// Example code from paper 5:
// _Computing with the analytic Jacobian of a genus 2 curve_
// by Paul B. van Wamelen
//
print "Examples from: Computing with the analytic Jacobian of a genus 2 curve";
load "Paper05_funcs";
ei := GetEchoInput();
SetEchoInput(true);
// Section 2: Finding genus 2 CM curves defined over the rationals
// Subsection: CM by the maximal order of Q(Sqrt(-2 + Sqrt(2)))
SetDefaultRealFieldPrecision(100);
Q := RationalField();
P<x> := PolynomialRing(Q);
R<s> := NumberField(x^2 - 2);
PP<x> := PolynomialRing(R);
RF := NumberField(x^2 - (-2 + s));
F<t> := AbsoluteField(RF);
O := MaximalOrder(F);
D := Different(O);
bool, gen := IsPrincipal(D);
xi := 1/gen;
xi in R;
xi^2 in R;
Conjugates(xi);
G, Aut, tau := AutomorphismGroup(F);
G;
exists(cc){ tau(g) : g in G | Order(g) eq 2 };
Z := IntegerRing();
E := Matrix(Z, 4, 4, [ Trace(xi*cc(a)*b) : b in Basis(O), a in Basis(O) ]);
D, C := FrobeniusFormAlternating(E); D;
newb := ElementToSequence(Matrix(O, C) * Matrix(O, 4, 1, Basis(O)));
C := ComplexField(100);
SetKantPrecision(O, 100);
PMat := Matrix(C, 2, 4,
  [ Conjugate(b, 2) : b in newb ] cat [ Conjugate(b, 4) : b in newb ]);
tau := Submatrix(PMat, 1, 3, 2, 2)^-1 * Submatrix(PMat, 1, 1, 2, 2);
S := RosenhainInvariants(tau);
KC<xc> := PolynomialRing(C);
f := xc*(xc - 1) * &*{ xc - a : a in S };
IC := IgusaClebschInvariants(f);
IC := [ 1, IC[2]/IC[1]^2, IC[3]/IC[1]^3, IC[4]/IC[1]^5 ];
ICp := [ BestApproximation(Re(r), 10^50) : r in IC ];
Maximum([ Abs(IC[i] - ICp[i]) : i in [1..#IC] ]);
ICp;
C1 := HyperellipticCurveFromIgusaClebsch(ICp);
C2 := ReducedModel(C1 : Al := "Wamelen");
C2;
f := Evaluate(HyperellipticPolynomials(C2), xc);
A := AnalyticJacobian(f);
IsIsomorphicSmallPeriodMatrices(tau, SmallPeriodMatrix(A));
MA := EndomorphismRing(SmallPeriodMatrix(A));
Dimension(MA);
MAGens := SetToSequence(Generators(MA)); MAGens;
NF_MAGens := NumberField(MinimalPolynomial(MAGens[1]));
if NF_MAGens eq Rationals() then
  NF_MAGens := NumberField(MinimalPolynomial(MAGens[2]));
end if;
IsIsomorphic(F, NF_MAGens);


// Subsection: CM by the maximal order of Q(Sqrt(-5 + Sqrt(5)))
Q := RationalField();
P<x> := PolynomialRing(Q);
R<s> := NumberField(x^2 - 5);
PP<x> := PolynomialRing(R);
RF := NumberField(x^2 - (-5 + s));
F<t> := AbsoluteField(RF);
O := MaximalOrder(F);
ClassNumber(F);
G, Aut, tau := AutomorphismGroup(F);
G;
exists(cc){ tau(g) : g in G | Order(g) eq 2 };
G, classmap := ClassGroup(F);
D := Different(O);
D @@ classmap;
I1 := classmap(G.1);
IsPrincipal(I1*cc(I1));
U, Umap := UnitGroup(O); U;
_, iso := IsIsomorphic(F,RF);
Uplus, Upmap := UnitGroup(R);
normUgens := { Norm(iso(F!Umap(g))) @@ Upmap : g in Generators(U) };
subUplus := sub< Uplus | normUgens >;
cosetreps := Transversal(Uplus, subUplus);
ups := { (RF!Upmap(cr)) @@ iso : cr in cosetreps };
ups;
uplus := 1/2*(t^2 + 6);
_, b1 := IsPrincipal(D); F!b1;
_, b2 := IsPrincipal(D*I1*cc(I1)); F!b2;
b1 eq -cc(b1);
b2 eq -cc(b2);
tau1 := xi2tau(1/b1, O, cc, 100);
tau2 := xi2tau(uplus/b1, O, cc, 100);
tau3 := xi2tau(1/b2, I1, cc, 100);
tau4 := xi2tau(uplus/b2, I1, cc, 100);
IsIsomorphicSmallPeriodMatrices(tau1, tau2);
bool := IsIsomorphicSmallPeriodMatrices(tau3, tau4); bool;
S := RosenhainInvariants(tau1);
f := xc*(xc-1) * &*{ xc - a : a in S };
IC := IgusaClebschInvariants(f);
IC := [ 1, IC[2]/IC[1]^2, IC[3]/IC[1]^3, IC[4]/IC[1]^5 ];
ICp := [ BestApproximation(Re(r), 10^50) : r in IC ];
Maximum([ Abs(IC[i] - ICp[i]) : i in [1..#IC] ]);
ICp;
C1 := HyperellipticCurveFromIgusaClebsch(ICp);
C2 := ReducedModel(C1 : Al := "Wamelen"); C2;
S := RosenhainInvariants(tau3);
f := xc*(xc-1) * &*{ xc - a : a in S };
IC := IgusaClebschInvariants(f);
IC := [ 1, IC[2]/IC[1]^2, IC[3]/IC[1]^3, IC[4]/IC[1]^5 ];
ICp := [ BestApproximation(Re(r), 10^50) : r in IC ];
ICp;
Maximum([ Abs(IC[i] - ICp[i]) : i in [1..#IC] ]);
C1 := HyperellipticCurveFromIgusaClebsch(ICp);
C2 := ReducedModel(C1 : Al := "Wamelen"); C2;


// Section 3: Isogenies
K<x> := PolynomialRing(RationalField());
C<i> := ComplexField(100);
KC<xc> := PolynomialRing(C);
f1 := x^5 + 40*x^4 + 136*x^3 + 96*x^2 + 16*x;
f1C := Evaluate(f1, xc);
A1 := AnalyticJacobian(f1C);
f2 := x^6 - 8*x^5 + 22*x^4 - 16*x^3 - 36*x^2 + 64*x - 24;
f2C := Evaluate(f2, xc);
A2 := AnalyticJacobian(f2C);
bool, M, alpha := IsIsogenous(A1, A2); bool;
alpha[1,1];
Mlst := AnalyticHomomorphisms(SmallPeriodMatrix(A1), SmallPeriodMatrix(A2));
Mlst;
P1 := BigPeriodMatrix(A1);
P2 := BigPeriodMatrix(A2);
alst := [ Submatrix(P2*Matrix(C, M), 1, 1, 2, 2)
          * Submatrix(P1, 1, 1, 2, 2)^-1 : M in Mlst ];
SetDefaultRealFieldPrecision(100);
Cp := ComplexField();
_<x> := PolynomialRing(IntegerRing());
plst := [];
for alpha in alst do
  for i in [1..2] do
    for j in [1..2] do
      pol := MyPowerRelation(Cp!alpha[i,j], 8);
      if Degree(pol) gt 1 then
        Append(~plst, pol);
      end if;
    end for;
  end for;
end for;
plst;
F := NumberField(plst[1]);
for i in [2..#plst] do
  F2 := NumberField(plst[i]);
  Flst := CompositeFields(F, F2);
  _,ind := Min([ Degree(G) : G in Flst ]);
  F := Flst[ind];
end for;
DefiningPolynomial(F);
O := MaximalOrder(F);
_, O2 := OptimizedRepresentation(O);
pol := DefiningPolynomial(O2); pol;
Kp<xp> := PolynomialRing(Cp);
aroot := Roots(Evaluate(pol, xp))[1][1];
basis := [ aroot^i : i in [0..7] ];
elst := [ [ inbase(basis, l) : l in ElementToSequence(alpha) ] : alpha in alst ];
elst;
mat := Matrix(4, 28, [ &cat[ Remove(e[i], 1) : i in [1..4] ] : e in elst ]);
Nullspace(mat);
alpha := &+[ alst[i] : i in Support(Basis(Nullspace(mat))[1]) ]; alpha;

C1 := HyperellipticCurve(f1);
J1 := Jacobian(C1);
pts1 := RationalPoints(J1 : Bound := 500);
pts1;
P1 := pts1[5];
apol := ElementToSequence(P1)[1];
bpol := ElementToSequence(P1)[2];
divs1 := [ r[1] : r in Roots(Evaluate(apol, xc)) ];
divs1 := [ <d1, Evaluate(bpol, d1)> : d1 in divs1 ];
pt1 := &+[ ToAnalyticJacobian(d[1], d[2], A1) : d in divs1 ];
P2 := FromAnalyticJacobian(alpha*pt1, A2);
clst := Coefficients((xc - P2[1][1])*(xc - P2[2][1]));
coefflst := [ BestApproximation(Re(c), 10^40) : c in clst ];
Maximum([ Abs(coefflst[i] - clst[i]) : i in [1..#clst] ]);
xpol := K!coefflst; xpol;
pntx := Pi(C) + i*67/109;
pnty := Sqrt(Evaluate(f1C,pntx));
pt := ToAnalyticJacobian(pntx, pnty, A1);
P2 := FromAnalyticJacobian(alpha*pt, A2);
(8*(2 + 3*pntx))/(4 + 8*pntx + pntx^2) - (P2[1][1] + P2[2][1]);
(8*(2 + 3*pntx))/(4 + 8*pntx + pntx^2) - (P2[1][1] * P2[2][1]);
SetEchoInput(ei);