// basic definitions
C10<z>:=CyclotomicField(10);
R<X,Y,Z>:=PolynomialRing(C10,3);
GalAction:=hom<C10 -> C10 | z^3>; // this is the automorphism \sigma : z_{10} \mapsto z_{10}^3 of Q(\zeta_{10})

// Equations of the curve
f2:=X^5*Z + 11*Z^6 - Y^5*Z; // the polynomial defining the curve in weighted projective space (together with XY=Z^2)
h:=X*Y-Z^2; // the equation of Q_0
I:=ideal<R | h>; // the ideal generated by the conic Q_0

// useful functions
function GalActionMatrix(M) // given a 3x3 matrix M, compute \sigma(M), where \sigma is the generator of Gal(Q(z_10)/Q) that sends z_10 to z_10^3
	N:=Matrix(C10,3,3,[]);
	for i in [1..3] do
		for j in [1..3] do
			N[i][j]:=GalAction(M[i][j]);
		end for;
	end for;
	return N;
end function;

function GalActionPoly(f) // given a polynomial f with coefficients in C10, compute \sigma(f)
	coeff:=Coefficients(f);
	monom:=Monomials(f);
	for i in [1..#coeff] do
		coeff[i]:=GalAction(coeff[i]);
	end for;
	return Polynomial(coeff,monom);
end function;

function EvaluatePolynomialInMatrix(p,M) // given a polynomial p in 3 variables and a 3x3 matrix M, computes p \circ M (x,y,z)
	transform:=Matrix(R,M)*Matrix(R,3,1,[X,Y,Z]);
	return Evaluate(p,[transform[1][1],transform[2][1],transform[3][1]]);
end function;

function poly_trace(f) // computes the trace of the polynomial f, that is, the polynomial whose coefficients are the traces of the coefficients of f
	coeff:=Coefficients(f);
	monom:=Monomials(f);
	for i in [1..#coeff] do
		coeff[i]:=Trace(coeff[i]);
	end for;
	return Polynomial(coeff,monom);
end function;


// the core stuff starts here
function compute_twist()
alpha:=Matrix(C10,3,3,[1,(z-1)^2,2*(z-1),z^2*(1-z)^2,z^4,-2*z^3*(1-z),z*(1-z),(z-1)*(-z^2),-z^2+z*(1-z)*(z-1)]); // the automorphism
lambda:=(-z^2+z*(1-z)*(z-1))/(z^4-(z-1)^2*z^2*(1-z)^2); // normalization factor for the embedding
alphaNorm:=lambda*alpha; // normalized matrix in GL_3

// Galois conjugates of alphaNorm
alphaNorm3:=GalActionMatrix(alphaNorm);  
alphaNorm9:=GalActionMatrix(alphaNorm3);
alphaNorm7:=GalActionMatrix(alphaNorm9);

// these are all automorphisms of order 3:
// alphaNorm^3;
// alphaNorm3^3;
// alphaNorm7^3;
// alphaNorm9^3;

// define the cocycle
cocyc1:=DiagonalMatrix(C10,3,[1,1,1]); // the values of the cocycle. cocyc(n) is the value of \xi(\tau), where \tau(z_10)=z_10^n
cocyc3:=alphaNorm;
cocyc9:=alphaNorm*alphaNorm3;
cocyc7:=alphaNorm*alphaNorm3*alphaNorm9;

// We check that cocyc(n) are automorphisms. The reason to do this explicitly is that just evaluating does not work:
"What's f \circ \xi_{\sigma}?", EvaluatePolynomialInMatrix(f2,cocyc3);

// however, they really are automorphisms, because we need to remember that the right condition is that we find the same polynomial up to a multiple of the conic:
"NormalForm after quotienting out by the ideal generated by the conic:";
NormalForm(EvaluatePolynomialInMatrix(f2,cocyc1),I);
NormalForm(EvaluatePolynomialInMatrix(f2,cocyc3),I);
NormalForm(EvaluatePolynomialInMatrix(f2,cocyc9),I);
NormalForm(EvaluatePolynomialInMatrix(f2,cocyc7),I);

// Set up an auxiliary matrix c (your M_0) to be used for Hilbert 90
c1:=DiagonalMatrix(C10,3,[1+z,1+z^3,1+z^7]); // a "sufficiently generic" matrix
c3:=GalActionMatrix(c1); // and its Galois conjugates
c9:=GalActionMatrix(c3);
c7:=GalActionMatrix(c9);

// Compute M as in Hilbert 90
M:=cocyc1*c1+cocyc3*c3+cocyc7*c7+cocyc9*c9; // a matrix given by Hilbert 90

// Now let's twist the curve. Compute first the twist of the polynomial f
Twistf:=EvaluatePolynomialInMatrix(f2,M); // this is f \circ M
SigmaTwistf:=GalActionPoly(Twistf); // this is \sigma(f \circ M); if f \circ M had rational coefficients, we would have Twistf=SigmaTwistf

// However this is not the case:
"The difference between \sigma(f \circ M) and f \circ M is:", SigmaTwistf-Twistf;

// Let's compute the twist of the conic:
twistConic:=EvaluatePolynomialInMatrix(h,M);

"New conic: ", twistConic;

I2:=ideal<R|twistConic>;

// the difference between Twistf and SigmaTwistf belongs to the ideal generated by the twisted conic
"Is the difference between f \circ M and \sigma(f \circ M) in the ideal generated by the conic?", Twistf-SigmaTwistf in I2;

// this implies that also the difference between the twisted polynomial and its (normalized) trace is in the ideal generated by the twisted conic
"Is the difference between f \circ M and its normalized trace in the ideal generated by the conic?", Twistf-1/4*poly_trace(Twistf) in I2;

// let's get an affine model of the twist and test that everything goes as planned. First we work over Q(\zeta_10), where this should be the same curve as before.

S<x,y,t>:=PolynomialRing(C10,3);
affineConic:=S!Evaluate(twistConic,[x,y,1]);
affinePolynomial:=t^2-Evaluate(1/4*poly_trace(Twistf),[x,y,1]);

CTwist:=Curve(AffineSpace(C10,3),[affineConic, affinePolynomial]);



return CTwist, twistConic, Twistf;
end function;






procedure compute_twist_and_check()
CTwist, twistConic, Twistf := compute_twist();
// The original curve
S<x,y>:=PolynomialRing(C10,2);
p:=y^2-x*(x^10+11*x^5-1);
C:=Curve(AffineSpace(C10,2),[p]);

"Original curve:", C;
"Twisted curve:", CTwist;

// Let's check that they really are isomorphic
checkIsomorphism:=IsIsomorphic(ProjectiveClosure(C),ProjectiveClosure(CTwist));
"Is CTwist isomorphic to C?", checkIsomorphism;


// Now let's work over the rationals. An interesting property of the twisted curve is that it has 12 rational automorphisms (in particular, one of order 3), whereas the original curve only had 4.

S<x,y,t>:=PolynomialRing(Rationals(),3);
affineConic:=S!Evaluate(twistConic,[x,y,1]);
affinePolynomial:=t^2-Evaluate(1/4*poly_trace(Twistf),[x,y,1]);

CTwist:=Curve(AffineSpace(Rationals(),3),[affineConic, affinePolynomial]);
G:=AutomorphismGroup(CTwist);
"How many rational automorphisms does CTwist have?", Order(G);
end procedure;

"Compute twist";
time(compute_twist());

"Compute twist, check isomorphism, and compute rational automorphism group";
time(compute_twist_and_check());
