// This Magma script will be used to try and construct
// automorphisms of a canonically embedded
// non-hyperelliptic curve defined over a number field.

// One tool will be the fact that for a canonically embedded curve,
// of genus g every automorphism comes from a projective transformation 
// of P^(g-1). 

// Overall strategy:
// 1. Write down the "automorphism scheme" X of a curve C defined over
// a number field K.
// 2. Use Magma's automorphism group code over a finite field to find mod p points on X.
// 3. Use Hensel's lemma to lift those points to p-adic points of X.
// 4. Find points in K that reduce to those p-adic points of X.

bigprec := 200;
prec := 30;

// Helper function - This routine takes a number field K and a p-adic
// number a. It tries to find a simple element of K that has the same
// p-adic expansion as a.

function guess(K,a)
  ChangePrecision(~a,prec);
  OK := LLL(MaximalOrder(K));
  p := Prime(Parent(a));
  I := ideal<OK | p>;
  P := Factorization(I)[1][1];
  if Norm(P) ne p then
    printf "Inertial degree is > 1. This is a problem.\n";
    bad := 0;
    bad2 := 1/bad;
  end if;
  Kp, phi := Completion(K,P : Precision := prec);
  cflist := [ Integers()!phi(OK.i) : i in [1..Degree(K)]];
  Append(~cflist,Integers()!(a));
  Append(~cflist,p^prec);
  // Scaling factor.
  // This example is very similar to that given in Lat_LLLXGCD
  S := 10^15;
  X := RMatrixSpace(Integers(),Degree(K)+2,Degree(K)+3)!0;
  for i in [1..Degree(K)+2] do
    X[i][i+1] := -1;
  end for;
  for i in [1..Degree(K)+2] do
    X[i][1] := S*cflist[i];
  end for;
  L := LLL(X : Proof := false);
  v := L[1];
  denom := -v[Degree(K)+2];
  ret := (1/denom)*(&+[ v[i+1]*OK.i : i in [1..Degree(K)]]);
  return ret;
end function;

F := Rationals();
dval := 2;
R<a,b,c,d,e> := PolynomialRing(F,5);
f1 := -2705*a^2 + 1681*b^2 - 1967*b*c + 2048*c^2 - dval*d^2;
f2 := 73*a^2 - 41*b^2 + 64*b*c - 64*c^2 - dval*d*e;
f3 := -2*a^2 + b^2 - 2*b*c + 2*c^2 - dval*e^2;
D := Curve(ProjectiveSpace(F,4),[f1,f2,f3]);

p := 31;
RR<x> := PolynomialRing(Integers());
f := x^4 - 4*x^2 + 2;
K<z> := NumberField(f);

D2 := D;
gen := Genus(D2);
S := PolynomialRing(F,gen);
de := DefiningEquations(D2);
de := [S!de[i] : i in [1..#de]];

// Make the equations integral if necessary.

mindeg := Min([Degree(de[i]) : i in [1..#de]]);
maxdeg := Max([Degree(de[i]) : i in [1..#de]]);
// For k between mindeg and maxdeg, make a basis for the space
// of degree k polynomials vanishing on D.

polybasis := [];
ind := [mindeg..maxdeg];
for k in [mindeg..maxdeg] do
  kbasis := [];
  if (k eq mindeg) then
    kbasis := [ de[i] : i in [1..#de] | Degree(de[i]) eq mindeg ];
  end if;
  if (k gt mindeg) then
    oldbasis := polybasis[#polybasis];
    // Multiply the old basis vectors by the variables
    // in the coordinate ring. Then add any new basis.
    kbasis := [ oldbasis[i]*S.j : i in [1..#oldbasis], j in [1..gen]];
    kbasis := kbasis cat [ de[i] : i in [1..#de] | Degree(de[i]) eq k ];
  end if;
  Append(~polybasis,kbasis);
end for;

numvars := gen^2 + &+[ #polybasis[Degree(de[i])-mindeg+1] : i in [1..#de]] + 1;
printf "The automorphism scheme has %o variables.\n",numvars;
// Polynomial ring used to define automorphism scheme
T := PolynomialRing(F,numvars);
// Polynomial ring over that
T2 := PolynomialRing(T,gen);
imlist := [];
for i in [1..gen] do
  Append(~imlist, &+[ T2!(T.((i-1)*gen+j))*T2.j : j in [1..gen]]);
end for;
eqlist := [];
varcount := gen^2;
for i in [1..#de] do
  poly := Evaluate(de[i],imlist);
  for j in [1..#polybasis[Degree(de[i])-mindeg+1]] do
    varcount := varcount + 1;
    poly := poly - T2!(T.varcount)*(T2!polybasis[Degree(de[i])-mindeg+1][j]);
  end for;
  cofs, monlist := CoefficientsAndMonomials(poly);
  for k in [1..#monlist] do
    if cofs[k] ne 0 then
      Append(~eqlist,cofs[k]);
    end if;
  end for;
end for;
mat := Matrix(T,gen,gen,[T.i : i in [1..gen^2]]);
Append(~eqlist,Determinant(mat)*T.numvars - 1);
printf "The automorphism scheme has %o equations.\n",#eqlist;
I := ideal<T | eqlist>;
X := Scheme(AffineSpace(F,numvars),Basis(I));

// Pick a prime ideal of the ring of integers of F.
// Compute the automorphism group modulo that prime.

X2 := ChangeRing(X,GF(p));

S := PolynomialRing(GF(p),gen);
varlist := [ S.i : i in [1..gen]];
newdes := [ Evaluate(de[i],varlist) : i in [1..#de]];
Pr<a,b,c,d,e> := ProjectiveSpace(GF(p),gen-1);
D3 := Curve(Pr,newdes);
ptlist := Points(D3);
// We need enough points on D3 to be able to determine
// the linear representation of the automorphisms from them.

A := AutomorphismGroup(D3);
G, phi0 := PermutationRepresentation(A);
psi := Inverse(phi0);
// Find the matrices
matlist := [];
count := 0;
ordlist := [];
for g in Generators(G) do
  count := count + 1;
  printf "Finding matrix for automorphism %o.\n",count;
  mp := psi(g);
  Append(~ordlist,Order(mp));
  M := ZeroMatrix(GF(p),0,gen^2);
  done := false;
  it := 0;
  N := 0;
  while (done eq false) do
    it := it + 1;
    pt1 := ptlist[it];
    for j in [2..gen] do
      vec := ZeroMatrix(GF(p),1,gen^2);
      pt2 := mp(pt1);
      for k in [1..gen] do
        vec[1][k] := pt1[k]*pt2[j];
        vec[1][(j-1)*gen+k] := -pt1[k]*pt2[1];
      end for;
      M := VerticalJoin(M,vec);
    end for;
    N := NullSpace(Transpose(M));
    if (Dimension(N) eq 1) then
      done := true;
    end if;
  end while;
  mat := Matrix(GF(p),gen,gen,Eltseq(N.1));
  Append(~matlist,mat);
end for;
// Use the matrices to create p-adic points on the scheme X.
ptlist := [];
Qp := pAdicField(p,bigprec);
ptset := PointSet(X,Qp);
elt := O(Qp!p);
for i in [1..#matlist] do
  M := matlist[i];
  imlist := [ &+[ M[i][j]*S.j : j in [1..gen]] : i in [1..gen]];
  ptonX := [];
  for j in [1..gen^2] do
    Append(~ptonX,Qp!Eltseq(M)[j] + elt);
  end for;
  for j in [1..#de] do
    polyind := Degree(de[j]) - mindeg + 1;
    // if k = #polybasis[polyind] then
    // we have a kxk linear system to solve. We want to plug
    // the equations for the automorphism into de[j] and write
    // the result as a linear combination of the polynomials in kbasis.
    polymatch := Evaluate(de[j],imlist);
    monoset := &join [ Set(Monomials(S!polybasis[polyind][k])) : k in [1..#polybasis[polyind]]];
    monoseq := SetToSequence(monoset);
    k := #polybasis[polyind];
    Amat := ZeroMatrix(GF(p),k,#monoseq);
    for l in [1..k] do
      for i in [1..#monoseq] do
        Amat[l][i] := MonomialCoefficient(S!polybasis[polyind][l],monoseq[i]);
      end for;
    end for;
    Bmat := ZeroMatrix(GF(p),1,#monoseq); 
    for i in [1..#monoseq] do
      Bmat[1][i] := MonomialCoefficient(polymatch,monoseq[i]);
    end for;
    tst, sol := IsConsistent(Amat,Bmat);
    if (tst eq false) then
      printf "Ack! The polynomial isn't in the set of polynomials it needs to be in.\n";
      bad := 0;
      bad2 := 1/bad;
    end if;
    for i in [1..k] do
      Append(~ptonX,Qp!sol[1][i] + elt);
    end for; 
  end for;
  Append(~ptonX,Qp!(Determinant(M)^(-1)) + elt);
  Append(~ptlist,ptset!ptonX);
end for;
printf "Found %o points on X(F_p).\n",#ptlist;
des := DefiningEquations(X);
Qp := pAdicField(p,bigprec);
T3 := PolynomialRing(Qp,numvars);
des := DefiningEquations(X);
polylist := [ T3!des[i] : i in [1..#des]];
liftedpoints := [];
for i in [1..#ptlist] do
  printf "Calling LiftPoint for i = %o.\n",i;
  ptlift := LiftPoint(polylist,1,Coordinates(ptlist[i]),prec);
  Append(~liftedpoints,ptlift);
end for;
// Now, reconstruct algebraic numbers using the entries from these lists.
matlist := [ Matrix(Qp,gen,gen,[ liftedpoints[i][j] : j in [1..gen^2]]) : i in [1..#liftedpoints]];
matlist2 := [];
for i in [1..#matlist] do
  test := &and [ Valuation((matlist[i]^ordlist[i])[j][k]) ge prec : j in [1..gen], k in [1..gen] | j ne k];
  if (test eq false) then
    printf "Error. It looks like matrix %o doesn't have order %o.\n",i,ordlist[i];
  end if;
  ent := (matlist[i]^ordlist[i])[1][1];
  bol, elt := IsPower(ent,ordlist[i]);
  if (bol eq true) then
    mat := matlist[i]/elt;
  else
    printf "Error. We can't lift automorphism %o into GL.\n",i;
    mat := matlist[i];
  end if;
  Append(~matlist2,mat);
end for;

matlist3 := [];
SS<yy> := PolynomialRing(Qp);
for i in [1..#matlist] do
  printf "Lifting matrix %o to characteristic zero.\n",i;
  curmat := ZeroMatrix(K,gen,gen);
  ord := ordlist[i];
  g := yy^ord - 1;
  L := Roots(g);
  done := false;
  it := 1;
  while done eq false do
    for j in [1..gen] do
      for k in [1..gen] do
        curmat[j][k] := guess(K,matlist2[i][j][k]/L[it][1]);
      end for;
    end for;
    if curmat^ord eq ScalarMatrix(K,gen,1) then
      done := true;  
      printf "Done!\n";
      Append(~matlist3,curmat);
    else
      it := it + 1;
      if (it gt #L) then
        printf "Error! Failed to find any characteristic zero lifts.\n";
      end if;
    end if;
  end while;
end for;
autlist := [];
matset := { ScalarMatrix(K,gen,1) };
done := false;
while done eq false do
  newmatset := { matlist3[i]*m : m in matset, i in [1..#matlist3] };
  newmatset := newmatset join matset;
  if #newmatset eq #matset then
    done := true;
  else
    matset := newmatset;
  end if;
end while;
matlist4 := [];
for m in matset do
  found := false;
  for j in [1..#matlist4] do
    chk := m*matlist4[j]^(-1);
    if IsDiagonal(chk) and &and[ chk[i][i] eq chk[1][1] : i in [2..gen]] then
      found := true;
    end if; 
  end for;
  if (found eq false) then
    Append(~matlist4,m);
  end if;
end for;
autlist := [];
DK<a,b,c,d,e> := BaseChange(D2,K);
for i in [1..#matlist4] do
  mat := matlist4[i];
  matinv := mat^(-1);
  maplist1 := [ &+[ mat[j][k]*DK.k : k in [1..gen]] : j in [1..gen]];
  maplist2 := [ &+[ matinv[j][k]*DK.k : k in [1..gen]] : j in [1..gen]];
  aut := iso<DK->DK| maplist1, maplist2>;
  Append(~autlist,aut);
end for; 
matset := { ScalarMatrix(K,gen,1) };
done := false;
while done eq false do
  newmatset := { matlist3[i]*m : m in matset, i in [1..#matlist3] };
  newmatset := newmatset join matset;
  if #newmatset eq #matset then
    done := true;
  else
    matset := newmatset;
  end if;
end while;
matlist4 := [];
for m in matset do
  found := false;
  for j in [1..#matlist4] do
    chk := m*matlist4[j]^(-1);
    if IsDiagonal(chk) and &and[ chk[i][i] eq chk[1][1] : i in [2..gen]] then
      found := true;
    end if; 
  end for;
  if (found eq false) then
    Append(~matlist4,m);
  end if;
end for;
// matlist 4 is a list of genxgen matrices that correspond to the linear
// transformations found in the automorphism group of the curve.