////////////////////////////////////////////////////////////////////////
//          FIRST CREATE SOME CURVES AND TEST ARITHMETIC              //
////////////////////////////////////////////////////////////////////////
p := 23;
FF := FiniteField(p);
PF<x> := PolynomialRing(FF);
////////////////////////////////////////////////////////////////////////
// g = 2
////////////////////////////////////////////////////////////////////////
// #0 (absolutely reducible)
C0 := SingularHyperellipticCurve(-(x^3+x+3)^2);
J0 := SingularPicardGroup(C0);
p0 := J0![x^2 + 3*x + 1,x];
(p^2-p+1)*p0 eq J0!0;
////////////////////////////////////////////////////////////////////////
// g = 3
////////////////////////////////////////////////////////////////////////
// #0 (absolutely reducible)
/*
// This absolutely reducible construction doesn't work for odd genus.
C0 := SingularHyperellipticCurve(-(x^4 + x + 2)^2);
J0 := SingularPicardGroup(C0);
p0 := J0![x^2 + 21*x + 2,x];
(p^2+1)*p0 eq J0!0;
*/
// #1
C1 := SingularHyperellipticCurve(x*(x^3+x+1)^2);
J1 := SingularPicardGroup(C1);
p1 := J1!C1![13,5,1];
q1 := (p-1)*p1;
(p^2+1)*(p-1)*p1 eq J1!0;
(p^2+1)*q1 eq J1!0;
// #2:
C2 := SingularHyperellipticCurve(-x*(x^3+x+1)^2);
J2 := SingularPicardGroup(C2); 
p2 := J2!C2![15,8,1];
q2 := (p+1)*p2; 
(p^2+1)*q2 eq J2!0;
// #3:
C3 := SingularHyperellipticCurve(x*(x^3+x+3)^2);
J3 := SingularPicardGroup(C3); 
p3 := Random(J3);
q3 := (p+1)*p3;
(p^2-p+1)*q3 eq J3!0;
////////////////////////////////////////////////////////////////////////
// g = 4 
////////////////////////////////////////////////////////////////////////
// #0
C0 := SingularHyperellipticCurve(-(x^5 + 16*x^2 + 14*x + 6)^2);
J0 := SingularPicardGroup(C0);
p0 := J0![x^2 + 15*x + 22,x];
(p^5+1)*p0 eq J0!0;
// #1 (1,1,2)-split
C1 := SingularHyperellipticCurve(-x*(x^4 + x + 2)^2);
J1 := SingularPicardGroup(C1);
(p^4-1)*Random(J1) eq J1!0;
// #2 (1,1,2)-split
C2 := SingularHyperellipticCurve(x*(x^4 + x + 2)^2);
J2 := SingularPicardGroup(C2);
(p^4-1)*Random(J2) eq J2!0;
////////////////////////////////////////////////////////////////////////
//       USE INDEX CALCULUS TO COMPUTE RELATIONS IN PICARD GROUP      //
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
// g = 6 or split of genus 7
////////////////////////////////////////////////////////////////////////
p := 17;
FF := FiniteField(p);
N1 := p^6-p^5+p^4-p^3+p^2-p+1;
N2 := p^6+p^5+p^4+p^3+p^2+p+1;
assert IsPrime(N1);
PF<x> := PolynomialRing(FF);
g := x^7 + 11*x^4 + 14*x^3 + 15*x^2 + 16*x + 10;
// #0
C0 := SingularHyperellipticCurve(3*g^2);
J0 := SingularPicardGroup(C0);
p0 := J0![x^2 + 14,x];
N1*p0 eq J0!0;
// Discrete log computation:
RelsC0 := [ ];
DivsC0 := [ Eltseq(p0) ];
while #RelsC0 lt #DivsC0 do
    SmoothRelations(J0,~DivsC0,~RelsC0,2,500);
    printf " #Rels/#Divs = %o/%o\n", #RelsC0, #DivsC0;
end while;
V := RSpace(Integers(),#DivsC0);
R := sub< V | [ &+[ t[2]*V.t[1] : t in rel ] : rel in RelsC0 ] >;
BasisMatrix(R);
// #1
C1 := SingularHyperellipticCurve(x*g^2);
J1 := SingularPicardGroup(C1);
p1 := Random(J1);
N1*(p+1)*p1 eq J1!0;
// Discrete log computation:
PntsC1 := [ ];
m := map< C1->ProjectiveSpace(FF,1) | [A.1,A.3] > where A := Ambient(C1);
for a in FF do
    pnts := RationalPoints((Codomain(m)![a,1])@@m);
    if #pnts ne 0 then Append(~PntsC1,pnts[1]); end if;
end for;
RelsC1 := [ ];
DivsC1 := [ Eltseq(J1!p) : p in PntsC1 ];
SmoothRelations(J1,~DivsC1,~RelsC1,2,1000);
degs := [ Degree(D[1]) : D in DivsC1 ]; 
scal := [ d eq 1 select 1 else 2^64 : d in degs ]; 
V := RSpace(Integers(),#DivsC1);
R := sub< V | [ &+[ t[2]*V.t[1] : t in rel ] : rel in RelsC1 ] >;
M := LLL(BasisMatrix(sub< V |
    [ &+[ scal[t[1]]*t[2]*V.t[1] : t in rel ] : rel in RelsC1 ] >));
Submatrix(M,1,1,#PntsC1,#PntsC1);
////////////////////////////////////////////////////////////////////////
// Larger prime: (p = 1 mod 3)
////////////////////////////////////////////////////////////////////////
p := 1123;
FF := FiniteField(p);
N0 := p^6+p^5+p^4+p^3+p^2+p+1;
N1 := p^6-p^5+p^4-p^3+p^2-p+1;
assert IsPrime(N1);
PF<x> := PolynomialRing(FF);
g := x^7 - x + 2;
// #0
C0 := SingularHyperellipticCurve(-g^2);
J0 := SingularPicardGroup(C0);
p0 := J0![x^4 + 962*x^3 + 560*x^2 + 629*x + 988,x];
N1*p0 eq J0!0;
// Discrete log computation:
RelsC0 := [ ];
DivsC0 := [ Eltseq(p0) ];
// while #RelsC0 lt #DivsC0 do
   SmoothRelations(J0,~DivsC0,~RelsC0,2,1000);
// printf " #Rels/#Divs = %o/%o\n", #RelsC0, #DivsC0;
// end while;
// V := RSpace(Integers(),#DivsC0);
// R := sub< V | [ &+[ t[2]*V.t[1] : t in rel ] : rel in RelsC0 ] >;
// BasisMatrix(R);
// #1
C1 := SingularHyperellipticCurve(-x*g^2);
J1 := SingularPicardGroup(C1);
p1 := Random(J1);
N1*(p+1)*p1 eq J1!0;
// Discrete log computation:
PntsC1 := [ ];
m := map< C1->ProjectiveSpace(FF,1) | [A.1,A.3] > where A := Ambient(C1);
for a in FF do
    pnts := RationalPoints((Codomain(m)![a,1])@@m);
    if #pnts ne 0 then Append(~PntsC1,pnts[1]); end if;
end for;
RelsC1 := [];
DivsC1 := [ Eltseq(J1!p) : p in PntsC1 ];
while #RelsC1 lt #DivsC1 do
    SmoothRelations(J1,~DivsC1,~RelsC1,2,1000);
    printf " #Rels/#Divs = %o/%o\n", #RelsC1, #DivsC1;
end while;
V := RSpace(Integers(),#DivsC1);
R := sub< V | [ &+[ t[2]*V.t[1] : t in rel ] : rel in RelsC1 ] >;
BasisMatrix(R);
// #2 (wrong twist)
C2 := SingularHyperellipticCurve(x*g^2);
J2 := SingularPicardGroup(C2);
p2 := Random(J2);
N0*(p-1)*p2 eq J2!0;
////////////////////////////////////////////////////////////////////////
// Even larger prime: (p = 2 mod 3)
////////////////////////////////////////////////////////////////////////
p := 2657;
FF := FiniteField(p);
N0 := p^6+p^5+p^4+p^3+p^2+p+1;
N1 := p^6-p^5+p^4-p^3+p^2-p+1;
assert IsPrime(N1);
PF<x> := PolynomialRing(FF);
g := x^7 - x + 6;
// #0
C0 := SingularHyperellipticCurve(3*g^2);
J0 := SingularPicardGroup(C0);
p0 := J0![x^2 + 2606*x + 1060,x];
N1*p0 eq J0!0;
// Discrete log computation:
RelsC0 := [ ];
DivsC0 := [ Eltseq(p0) ];
while #RelsC0 lt #DivsC0 do
    SmoothRelations(J0,~DivsC0,~RelsC0,2,500);
    printf " #Rels/#Divs = %o/%o\n", #RelsC0, #DivsC0;
end while;
V := RSpace(Integers(),#DivsC0);
R := sub< V | [ &+[ t[2]*V.t[1] : t in rel ] : rel in RelsC0 ] >;
BasisMatrix(R);
// #1
C1 := SingularHyperellipticCurve(x*g^2);
J1 := SingularPicardGroup(C1);
p1 := Random(J1);
N1*(p^2-1)*p1 eq J1!0;
// Discrete log computation:
PntsC1 := [ ];
m := map< C1->ProjectiveSpace(FF,1) | [A.1,A.3] > where A := Ambient(C1);
for a in FF do
    pnts := RationalPoints((Codomain(m)![a,1])@@m);
    if #pnts ne 0 then Append(~PntsC1,pnts[1]); end if;
end for;
RelsC1 := [ ];
DivsC1 := [ Eltseq(J1!p) : p in PntsC1 ];
while #RelsC1 lt #DivsC1 do
    SmoothRelations(J1,~DivsC1,~RelsC1,2,200);
    printf " #Rels/#Divs = %o/%o\n", #RelsC1, #DivsC1;
end while;
V := RSpace(Integers(),#DivsC1);
R := sub< V | [ &+[ t[2]*V.t[1] : t in rel ] : rel in RelsC1 ] >;
BasisMatrix(R);
////////////////////////////////////////////////////////////////////////
//    MULTIPLICATIVE GROUP AND GENERALIZED JACOBIANS ISOMORPHISMS     //
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
q := 1000099;
////////////////////////////////////////////////////////////////////////
ExtendToFq6 := true;
////////////////////////////////////////////////////////////////////////
Fq1 := FiniteField(q);
PFq1<x1> := PolynomialRing(Fq1);
ff := PrimitivePolynomial(Fq1,3);
Fq3<delta> := ext< Fq1 | ff >;
PFq3<x3> := PolynomialRing(Fq3);
Fq6<gamma> := ext< Fq3 | x3^2 - delta>;
PFq6<x6> := PolynomialRing(Fq6);
////////////////////////////////////////////////////////////////////////
if ExtendToFq6 then
    PFq1<x1> := PFq6;
    PFq3<x3> := PFq6;
end if;
////////////////////////////////////////////////////////////////////////
gg := (PFq3!Eltseq(ff)) div (x3 - delta);
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
Gm<T> := AffineSpace(Fq6,1);
EE := SingularCubicCurve([Fq3|0,delta,0,0,0]);
FF := SingularHyperellipticCurve(x3^2*(x3 + delta));
JF := SingularPicardGroup(FF);
////////////////////////////////////////////////////////////////////////
HH := SingularHyperellipticCurve(x1*Evaluate(ff,x1)^2);
HH := SingularHyperellipticCurve(x6*Evaluate(ff,x6)^2);
JH := SingularPicardGroup(HH);
////////////////////////////////////////////////////////////////////////
// ISOMORPHISMS
////////////////////////////////////////////////////////////////////////
iEEtoFF := map< Parent(EE!0)->FF |
    P :-> FF!Eltseq(P), Q :-> EE!Eltseq(Q) >;
iEEtoJF := map< Parent(EE!0)->JF |
    P :-> JF![x-SP[1],SP[2]] where SP := Eltseq(P), 
    Q :-> EE![-SQ[1](0),SQ[2](0)] where SQ := Eltseq(Q) >;
iFFtoJF := map< FF->JF | P :-> JF![x3-P[1],P[2]],
    Q :-> FF![-SQ[1](0),SQ[2](0)] where SQ := Eltseq(Q) >;
////////////////////////////////////////////////////////////////////////
// HOMOMORPHISMS
////////////////////////////////////////////////////////////////////////
mGmtoEE := map< Gm->EE | P :-> [ delta*(Z^2-1), gamma^3*Z*(Z^2-1) ]
    where Z := (1+T)/(1-T) where T := P[1] >;
mEEtoGm := map< EE->Gm | P :-> [ (y-gamma*x)/(y+gamma*x) ]
    where y := P[2] where x := P[1] >;
mGmtoFF := map< Gm->FF | P :-> FF![ delta*(Z^2-1), gamma^3*Z*(Z^2-1), 1 ]
    where Z := (1+T)/(1-T) where T := P[1] >;
////////////////////////////////////////////////////////////////////////
// The maps from HH -> EE and HH -> EE -> Gm induce group homomorphisms
// on the divisor class groups, and JH = Pic^0(HH) -> EE but the pullback
// from EE to the group Pic^0(HH) does not.
////////////////////////////////////////////////////////////////////////
mEEtoHH := map< EE->HH |
    P :-> HH(KK)![x+delta,y*Evaluate(gg,x+delta)]
    where KK := Parent(x) where x, y := Explode(Eltseq(P)) >;
mGmtoHH := map< Gm->HH |
    [ X, gamma^3*Z*(Z^2-1)*Evaluate(gg,X), 1 ]
    where X := delta*Z^2 where Z := (1+T)/(1-T) >;
mGmtoJH := map< Gm->JH | P :-> JH!mGmtoHH(P) >;
////////////////////////////////////////////////////////////////////////
// Here we construct the isomorphisms of divisor class groups.
////////////////////////////////////////////////////////////////////////
mJHtoJF := map< JH->JF |
    P :-> [ Evaluate(a,x3+delta) : a in [D[1],D[2]*InverseMod(gg,D[1])] ]
    where D := Eltseq(P) >;
////////////////////////////////////////////////////////////////////////
function mJFtoJH(P)
    u, v, w := Explode(Eltseq(P@@iFFtoJF));
    assert w eq 1;
    if u eq 0 then return JF![x,0]; end if;
    t := Fq3!(v/u);
    cx := PFq1!Eltseq(t);
    ax := cx^2-x1;
    ax /:= LeadingCoefficient(ax);
    bx := (cx*ff) mod ax;
    return JH![ax,bx];
end function;
////////////////////////////////////////////////////////////////////////
dlog := Random([1..q^2-q+1]);
g := PrimitiveElement(Fq6)^((q^3 - 1)*(q + 1));
h := g^dlog;
////////////////////////////////////////////////////////////////////////
Gg := Gm![g];
Pg := mGmtoEE(Gg);
Qg := iFFtoJF(mGmtoFF(Gg));
Rg := mJFtoJH(Qg);
////////////////////////////////////////////////////////////////////////
Gh := Gm![h];
Ph := mGmtoEE(Gh);
Qh := iFFtoJF(mGmtoFF(Gh));
Rh := mJFtoJH(Qh);
////////////////////////////////////////////////////////////////////////
assert g^dlog eq h;
assert dlog*Qg-Qh eq JF!0;
assert dlog*Rg-Rh eq JH!0;
////////////////////////////////////////////////////////////////////////
// These maps require the base extension of HH to Fq6.
////////////////////////////////////////////////////////////////////////
if ExtendToFq6 then
    Sg := JH!mGmtoHH(Gg);
    Sh := JH!mGmtoHH(Gh);
    dlog*Sg-Sh ne JH!0;
    mJHtoJF(dlog*Sg-Sh) eq JF!0;
end if;
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
//                         FROBENIUS MORPHISM                         //
////////////////////////////////////////////////////////////////////////
function Frobenius(P,q,i)
    J := Parent(P);
    C := Curve(J);
    f, g, r := HyperellipticPolynomials(C); 
    a1, b1 := Explode(Eltseq(P));
    acoeffs := Eltseq(a1); 
    bcoeffs := Eltseq(b1);
    P := Parent(a1);
    aq := P![ c^q^i : c in acoeffs ];
    bq := P![ c^q^i : c in bcoeffs ];
    return J![aq,bq];
end function;
////////////////////////////////////////////////////////////////////////
// Note that while Tr_Sg and Tr_Sh are defined over the prime field F_q,
// there images in JF are not.  Moreover, since the map is not defined
// over F_q its kernel is not stablized by Frobenius, so the image of
// Galois conjugates of dlog*Sg-Sh are no longer zero.
////////////////////////////////////////////////////////////////////////
if ExtendToFq6 then
    Tr_Sg := &+[ Frobenius(Sg,q,i) : i in [0..2] ];
    Tr_Sh := &+[ Frobenius(Sh,q,i) : i in [0..2] ];
    dlog*Tr_Sg-Tr_Sh ne JH!0;
    mJHtoJF(Frobenius(dlog*Sg-Sh,q,1)) ne JF!0;
    mJHtoJF(dlog*Tr_Sg-Tr_Sh) ne JF!0;
end if;
////////////////////////////////////////////////////////////////////////