elltateform(E) =
{
	local(u,r,s,t);
	u=1/6; r=-1/12*E.b2;
	s=-1/2*E.a1; t=1/2*(1/12*E.a1*E.b2-E.a3); 
	return([u,r,s,t]);
}

ellcomposetrans(T12,T01) =
{
	local(u1,u2,r1,r2,s1,s2,t1,t2,u,r,s,t);
	u1 = T01[1]; r1 = T01[2]; s1 = T01[3]; t1 = T01[4];
	u2 = T12[1]; r2 = T12[2]; s2 = T12[3]; t2 = T12[4];
	u = u1*u2;
	r = u1^2*r2 + r1;
	s = u1*s2 + s1;
	t = u1^3*t2 + u1^2*s1*r2 + t1;
	return([u,r,s,t]);
}

ellinverttrans(T) = 
{
	local(u1,u2,r1,r2,s1,s2,t1,t2);
	u1 = T[1]; r1 = T[2]; s1 = T[3]; t1 = T[4];
	u2 = 1/u1;
	r2 = -r1/u1^2;
	s2 = -s1/u1;
	t2 = -(s1*r2*u1^2 + t1)/u1^3;
	return([u2,r2,s2,t2]);
}

vectoquad(V) =
{
	local(D,w,t,n,r,r1,r2,rootD0);
	if(V[3] % 4 == 2 || V[3] % 4 == 3, D = V[3]*4, D = V[3]);
	w = quadgen(D);
	t = trace(w); n = norm(w);
	r = V[3]/(t^2 - 4*n);
	r1 = numerator(r); r2 = denominator(r);
	if(!issquare(r1) || !issquare(r2), error("error in vectoquad: r nonsquare"));
	rootD0 = sqrtint(r1)/sqrtint(r2)*(2*w - t);
	return(V[1] + V[2]*rootD0);
}

quadtovec(quad,D) = 
{
	local(quadpol,n,t,w,r,V);
	if(type(quad) != "t_QUAD", 
		V = quad;
	,
		quadpol = component(quad,1);
		n = component(quadpol,1);
		t = -component(quadpol,2);
		if(D == 0, D = t^2 - 4*n);
		r = (t^2 - 4*n)/D;
		r1 = numerator(r); r2 = denominator(r);
		if(!issquare(r1) || !issquare(r2), error("error in quadtovec: r nonsquare"));
		w = [t/2,1/2*sqrtint(r1)/sqrtint(r2),D];
		V = [real(quad) + imag(quad)*w[1], imag(quad)*w[2],D];
	);
	return(V);
}

quadtovecpt(P,D) = 
{
	return([quadtovec(P[1],D),quadtovec(P[2],D)]);
}

elltwist(E0,D) =
{
	local(E1,E2,E3,T01,T12,T23,T03,L);
	T01 = elltateform(E0);
	E1 = ellchangecurve(E0,T01);
	if(issquare(D),
		T12 = [1,0,0,0];
	,
		T12 = [vectoquad([0,1,D]),0,0,0];
		E2 = ellinit([0,0,0,E1.a4/D^2,E1.a6/D^3]);
	);
	if(
		imag(E2.a1) != 0 || imag(E2.a2) != 0 || imag(E2.a3) != 0 ||
		imag(E2.a4) != 0 || imag(E2.a6) != 0
	,
		error("error in elltwist: imaginary parts not all zero")
	);
	L = [real(E2.a1),real(E2.a2),real(E2.a3),real(E2.a4),real(E2.a6)];
	E2 = ellinit(L);
	T23 = ellglobalred(E2)[2];
	E3 = ellchangecurve(E2,T23);
	L = [E3.a1,E3.a2,E3.a3,E3.a4,E3.a6];
	T03 = ellcomposetrans(T23,ellcomposetrans(T12,T01));
	return([T03,L]);
}

mwrankline = concat(componentfolder,"cremonaware/mwrank_ntl");

ellrank(E) =
{
	local(L);
	write("mwrank.in",[E.a1,E.a2,E.a3,E.a4,E.a6]);
	write("mwrank.in",[0,0,0,0,0]);
	system(mwrankline " -q -v 0 -o < mwrank.in > mwrank.out1");
	system("tail -2 mwrank.out1 > mwrank.out2");
	L = read("mwrank.out2");
	system("rm mwrank.in mwrank.out1 mwrank.out2");
	return(L);
}

applyellchangepoint(V,v,D) =
{
	local(Vout,j);
	Vout = vector(length(V),j,quadtovecpt(ellchangepoint(V[j],v),D));
	return(Vout);
}

twistgens(E0,D) =
{
	local(L,Etw,twtrans,torsdata,freedata,G);
	L = elltwist(E0,D);
	Etw = ellini(L[2]);
	twtrans = ellinverttrans(L[1]);
	torsdata = elltors(Etw);
	freedata = ellrank(Etw);
	G = [[torsdata[1],torsdata[2],applyellchangepoint(torsdata[3],twtrans,D)],[freedata[1],applyellchangepoint(freedata[2],twtrans,D)]];
	return(G);
}