// Fermat Twists

load "TwistsCommon.m";    // load code common to Fermat and Klein cases

Q:=Rationals();
RT<T>:=PolynomialRing(Q);
Rxx<x,y,z>:=PolynomialRing(Q,3);
FM<i>:=NumberField(T^2+1);

G<s,t,u,zz>:=FreeGroup(4);
FG,Fpi:=quo<G|s^3, t^4, u^4, (s^2*t^3)^2, (s^2*u^3)^3, u^3*t^3*s*u*s^2*t, t^3*s^2*u^2*s*t^3, zz^2, zz*s*zz*s^2, zz*t*zz*t^3, zz*u*zz*u>;
FG0:=sub<FG|s,t,u>;
FGP,Fmu:=PermutationGroup(FG);
Fphi:=map<G->FGP|x:->Fmu(Fpi(x))>;
FG0P := sub<FGP|Fphi(s),Fphi(t),Fphi(u)>;
assert IdentifyGroup(FGP) eq <192,956>;
assert IdentifyGroup(FG0P) eq <96,64>;

// Over K the Jacobian of every twist is isogenous to the cube of the elliptic curve below
Fec := EllipticCurve("64a4");

// Data format is <n, H0 generators, h, ID(H), ID(H0), f(x,y,z), abs(disc(f)), g(t), h(t)>
//   n = # twist, f(x,y,z)=0 is the curve equation, g(t) defines the base field k (either Q or a quadratic field), where g(t)=t for Q
//   the splitting field of h(t) is the field K=L over which all the endomorphisms/homomorphisms are defined (it necessarily contains the field k defined by g(t)

Ftwists:=<
<1, [Id(G)], Id(G), <2, 1>, <1, 1>, x^4+y^4+z^4, 2^40, T, T^2+1>,
<2, [Id(G)], u^2*t, <2, 1>, <1, 1>, x^4-6*x^2*y^2+y^4-2*z^4, 2^67, T, T^2+1>,
<3, [Id(G)], t^3*u*t*u, <2, 1>, <1, 1>, 4*x^4-y^4-z^4, 2^58, T, T^2+1>,
<4, [t^2], t, <4, 1>, <2, 1>, 12*x^4+40*x^3*y-100*x*y^3-75*y^4-2*z^4, 2^67*5^18, T^2+5, T^4-T^2-1>,
<5, [t^2], t^3*u*t*u, <4, 2>, <2, 1>, 9*x^4+9*y^4-4*z^4, 2^58*3^36, T, T^4-T^2+1>,
<6, [t^2], u^2*t, <4, 2>, <2, 1>, 9*x^4-54*x^2*y^2+9*y^4-2*z^4, 2^67*3^36, T, T^4-T^2+1>,
<7, [t^2], Id(G), <4, 2>, <2, 1>, 9*x^4+y^4+z^4, 2^40*3^18, T, T^4-T^2+1>,
<8, [t^2], u, <4, 2>, <2, 1>, 9*x^4-4*y^4+z^4, 2^58*3^18, T, T^4-T^2+1>,
<9, [u^2*t], Id(G), <4, 2>, <2, 1>, 2*x^4+36*x^2*y^2+18*y^4+z^4, 2^76*3^18, T, T^4-T^2+1>,
<10, [u^2*t], t^3*u*t*u, <4, 2>, <2, 1>, 9*x^4+18*x^2*y^2+y^4-2*z^4, 2^67*3^18, T, T^4-T^2+1>,
<12, [s], u^2*t, <6, 1>, <3, 1>, x^4+4*x^3*y+12*x^2*y^2-12*x^2*y*z-6*x^2*z^2+36*x*y*z^2+6*y^4-36*y^2*z^2-12*y*z^3+9*z^4, 2^67*3^36, T, T^3-3*T-4>,
<11, [s], Id(G), <6, 2>, <3, 1>, 5*x^4+8*x^3*y-4*x^3*z+6*x^2*y^2+12*x^2*z^2+12*x*y*z^2+4*x*z^3+2*y^4+4*y^3*z+6*y^2*z^2+4*y*z^3+2*z^4, 2^40*7^18, T, T^6+5*T^4+6*T^2+1>,
<30, [t^3*u*t*u^3], t*u, <8, 1>, <4, 1>, x^4-10*x^3*z+30*x^2*z^2-2*y^4-100*z^4, 2^49*5^18, T^2+5, T^8-2*T^4-4>,
<24, [t], t^3*u*t*u, <8, 2>, <4, 1>, 3*x^4-4*x^3*y+12*x^2*y^2+4*x*y^3+3*y^4-5*z^4, 2^40*5^27, T, T^8+15*T^4+25>,
<27, [t^3*u*t*u^3], u^2*t, <8, 2>, <4, 1>, 12*x^4+40*x^3*y-100*x*y^3-75*y^4+10*z^4, 2^67*5^27, T, T^8+15*T^4+25>,
<29, [t], Id(G), <8, 2>, <4, 1>, 3*x^4+4*x^3*y+12*x^2*y^2-4*x*y^3+3*y^4+20*z^4, 2^58*5^27, T, T^8+15*T^4+25>,
<21, [u^2*t,t^2], t^3*u*t*u^3, <8, 3>, <4, 2>, 11*x^4+12*x^3*y+54*x^2*y^2-12*x*y^3+11*y^4-2*z^4, 2^85*5^18, T, T^4-6*T^2+10>,
<22, [t^2,u^2], u^2*t, <8, 3>, <4, 2>, x^4+5*x^3*y-25*x*y^3-25*y^4+z^4, 2^22*5^18, T, T^4-6*T^2+10>,
<28, [u^2*t,t^2], u^2, <8, 3>, <4, 2>, 19*x^4-12*x^3*y+6*x^2*y^2+12*x*y^3+19*y^4+2*z^4, 2^85*5^18, T, T^4-6*T^2+10>,
<18, [t], u^2, <8, 3>, <4, 1>, 9*x^4-18*x^2*y^2-12*x*y^3-2*y^4+12*z^4, 2^58*3^27, T, T^4-6*T^2+12>,
<19, [t], t^3*u*t*u^3, <8, 3>, <4, 1>, 9*x^4-18*x^2*y^2-12*x*y^3-2*y^4-3*z^4, 2^40*3^27, T, T^4-6*T^2+12>,
<25, [t^3*u*t*u^3], u, <8, 3>, <4, 1>, 9*x^4+3*y^4-4*z^4, 2^58*3^27, T, T^4-6*T^2+12>,
<26, [t^3*u*t*u^3], Id(G), <8, 3>, <4, 1>, 9*x^4+3*y^4+z^4, 2^40*3^27, T, T^4-6*T^2+12>,
<16, [t^3*u*t*u], u^2*t, <8, 3>, <4, 1>, x^4-6*x^2*y^2+y^4-6*z^4, 2^67*3^9, T, T^4-6*T^2+12>,
<20, [t^3*u*t*u], u, <8, 3>, <4, 1>, 3*x^4-4*y^4+z^4, 2^58*3^9, T, T^4-6*T^2+12>,
<31, [t^3*u*t*u], Id(G), <8, 3>, <4, 1>, 3*x^4+y^4+z^4, 2^40*3^9, T, T^4-6*T^2+12>,
<14, [t^3*u*t*u], t, <8, 4>, <4, 1>, 3*x^3*y-3*x*y^3-2*z^4, 2^31*3^18, T^2+2, T^8+9>,
<13, [t^2,u^2], Id(G), <8, 5>, <4, 2>, 9*x^4+25*y^4+z^4, 2^40*3^18*5^18, T, T^8-3*T^6+8*T^4-3*T^2+1>,
<15, [t^2,u^2], u, <8, 5>, <4, 2>, 9*x^4+25*y^4-4*z^4, 2^58*3^18*5^18, T, T^8-3*T^6+8*T^4-3*T^2+1>,
<17, [u^2*t,t^2], t^3*u*t*u, <8, 5>, <4, 2>, x^4+30*x^2*y^2+25*y^4-18*z^4, 2^67*3^18*5^18, T, T^8-3*T^6+8*T^4-3*T^2+1>,
<23, [u^2*t,t^2], Id(G), <8, 5>, <4, 2>, 2*x^4+60*x^2*y^2+50*y^4+9*z^4, 2^67*3^18*5^18, T, T^8-3*T^6+8*T^4-3*T^2+1>,
<32, [s,u^2*t], Id(G), <12, 4>, <6, 1>, 4*x^3*y-3*x^2*z^2+12*x*y^2*z-2*y^4-2*y*z^3, 2^31*3^27, T, T^6+2*T^3+2>,
<33, [t^3*u*t*u,u^2], u*t, <16, 6>, <8, 2>, x^4-30*x^2*y^2-80*x*y^3-55*y^4-2*z^4, 2^67*5^9, T^2+5, T^8-2*T^4+5>,
<38, [t^3*u*t*u^3,u^2*t], u, <16, 7>, <8, 3>, x^4-12*x^2*y^2-32*x*y^3-28*y^4+z^4, 2^76, T, T^8-6*T^4-8*T^2-1>,
<41, [t*u^2*t*u*t], t*u*t*u^2, <16, 7>, <8, 1>, 2*x^3*y-x*y^3-z^4, 2^31, T, T^8-8*T^4-2>,
<37, [u^2*t*u], u, <16, 8>, <8, 1>, x^3*y+2*x*y^3+z^4, 2^31, T, T^8-2>,
<45, [t^3*u*t*u^3,t], u, <16, 8>, <8, 4>, x^4+10*x^3*y+30*x^2*y^2-100*y^4-10*z^4, 2^49*5^27, T, T^8-10*T^4-100>,
<40, [t,u^2], u*t*u, <16, 11>, <8, 3>, 4*x^4+4*x^3*y+6*x^2*y^2-2*x*y^3+y^4-2*z^4, 2^49*3^18, T,  T^8-2*T^4+9>,
<43, [t,u^2], Id(G), <16, 11>, <8, 3>, 4*x^4+4*x^3*y+6*x^2*y^2-2*x*y^3+y^4+2*z^4, 2^49*3^18, T, T^8-2*T^4+9>,
<44, [t^3*u*t*u^3,u^2*t], Id(G), <16, 11>, <8, 3>, 5*x^4-8*x^3*y+12*x^2*y^2+16*x*y^3+20*y^4+2*z^4, 2^85*3^18, T, T^8-2*T^4+9>,
<36, [t^3*u*t*u,u^2], Id(G), <16, 11>, <8, 2>, 9*x^4+5*y^4+z^4, 2^40*3^18*5^9, T, T^8+5*T^4+25>,
<34, [t^3*u*t*u,u^2], u, <16, 11>, <8, 2>, 9*x^4+5*y^4-4*z^4, 2^58*3^18*5^9, T, T^8+5*T^4+25>,
<42, [u*t*u,t^3], u^2, <16, 11>, <8, 2>, 7*x^4+8*x^3*y+6*x^2*y^2+8*x*y^3+7*y^4+10*z^4, 2^67*3^18*5^9, T, T^8+5*T^4+25>,
<35, [t^3*u*t*u,u^2], u^2*t, <16, 13>, <8, 2>, x^4+5*x^3*y-25*x*y^3-25*y^4+2*z^4, 2^31*5^18, T, T^8-8*T^4+25>,
<46, [t,u*t*u], Id(G), <16, 13>, <8, 2>, 19*x^4+32*x^3*y+21*x^2*y^2+8*x*y^3+3*y^4+2*y^3*z+6*y^2*z^2+8*y*z^3+4*z^4, 2^58*5^18, T, T^8-8*T^4+25>,
<39, [t^3*u*t*u^3,t], Id(G), <16, 13>, <8, 4>, x^4-2*x^3*y+6*x^2*y^2+4*x*y^3+4*y^4+3*z^4, 2^40*3^27, T,  T^8+12*T^4+9>,
<48, [s,u^2], u^2*t, <24, 12>, <12, 3>, x^4-3*x^3*z-12*x^2*y*z+16*x*y^3-x*z^3+9*y^4+12*y^3*z+6*y^2*z^2, 2^22*3^36, T, T^4-16*T-24>,
<47, [s,u^2], Id(G), <24, 13>, <12, 3>, 3*x^4+4*x^3*y+4*x^3*z+6*x^2*y^2+6*x^2*z^2+8*x*y^3+12*x*y*z^2+5*y^4+4*y^3*z+12*y^2*z^2+z^4, 2^40*7^18, T, T^6-T^4-2*T^2+1>,
<50, [t*u^2*t*u*t,u^2], Id(G), <32, 7>, <16, 6>, 4*x^4-8*x^3*y+12*x^2*y^2+2*y^4+5*z^4, 2^58*5^18, T, T^8-10*T^4+20>,
<52, [u,u^2*t*u^3*t], u^2*t, <32, 11>, <16, 2>, x^4-30*x^2*y^2-80*x*y^3-55*y^4-2*z^4, 2^67*5^9, T, T^8-2*T^4+5>,
<49, [u,t^3*u*t], Id(G), <32, 34>, <16, 2>, 2*x^4+3*y^4+z^4, 2^49*3^9, T, T^16-4*T^12+6*T^8+20*T^4+1>,
<51, [t^3*u*t*u,u^2,t], u, <32, 43>, <16, 13>, 3*x^4-36*x^2*y^2-96*x*y^3-84*y^4+2*z^4, 2^85*3^18, T, T^8+6*T^4-9>,
<53, [t*u^2*t*u*t,u^2], u, <32, 43>, <16, 6>, 9*x^4+36*x^3*y-24*x*y^3-4*y^4-10*z^4, 2^67*3^18*5^18, T, T^8-10*T^4+45>,
<54, [u*t*u,u^2,t], Id(G), <32, 49>, <16, 13>, x^4+x^3*y+24*x^2*y^2+67*x*y^3+79*y^4+2*z^4, 2^31*3^18*7^18, T, T^8+8*T^4+9>,
<55, [s,t], Id(G), <48, 48>, <24, 12>, 3*x^4+2*x^3*z+6*x^2*y*z+12*x*y^3-30*x*y^2*z+2*x*z^3-27*y^4+38*y^3*z+18*y^2*z^2-10*y*z^3, 2^40*3^18*7^18, T, T^6-T^4+5*T^2+1>,
<56, [t,u], Id(G), <64, 134>, <32, 11>, x^4-42*x^2*y^2-168*x*y^3-203*y^4+z^4, 2^67*7^9, T, T^8-4*T^4-14>,
<57, [s,u], u^2*t, <96, 64>, <48, 3>, 6*x^3*z+3*x^2*y^2+27*x^2*z^2+6*x*y^3+18*x*y*z^2+4*y^4+2*y^3*z+18*y*z^3-36*z^4, 2^49*3^36, T, T^12+6*T^4+4>,
<58, [s,u], Id(G), <96, 72>, <48, 3>, x^3*y-x^3*z+3*x^2*y^2+28*x*y^3-84*x*y^2*z+84*x*y*z^2+28*y^4-56*y^3*z+98*z^4, 2^31*7^18, T, T^12+T^8-2*T^4-1>,
<59, [s,t,u], Id(G), <192, 956>, <96, 64>,44*x^4+120*x^3*y+36*x^3*z+60*x^2*y*z+9*x^2*z^2-200*x*y^3+x*z^3-150*y^4-15*y^2*z^2, 2^58*3^18*7^18, T, T^12+48*T^4+64>>;
//4*x^4-8*x^3*y+4*x^3*z+36*x^2*y*z-3*x^2*z^2+24*x*y^3+24*x*y^2*z+3*x*z^3+10*y^4+8*y^3*z+9*y^2*z^2+2*y*z^3, 2^31*23^18, T, T^12+4*T^10+3*T^8-4*T^6-21*T^4+24*T^2+1>>;

FG0PC,FG0PCphi := PCGroup(FG0P);
Frep := AbsolutelyIrreducibleModulesSchur(FG0PC,Rationals());

// given a triple of nonnegative integers v computes the expectation of a1^v[1]*a2^v[2]*a3^v[3] over the Sato-Tate group corresponding to H,H0,j
// j=1 means over the base field k, j=2 means base change to the compositum of k with the CM field

function Fexp(G0,H,H0,j,v)
    assert j in [1,2];
    assert #v eq 3;

    R<u,ubar,a1,a2,a3,a1bar,a2bar,a3bar>:=PolynomialRing(CyclotomicField(4),8);
    RR:=quo<R|u*ubar-1>;    

    A1:=a1*u+a1bar*ubar;
    A2:=a2*u^2+a2bar*ubar^2+a1*a1bar;
    A3:=a3*u^3+a3bar*ubar^3+a1bar*a2*u+a1*a2bar*ubar;
    grancoef:=Coefficient(Coefficient(RR!(A1^(v[1])*A2^(v[2])*A3^(v[3])),u,0),ubar,0);

    e:=0;
    for h in H0 do
        A:=Coefficients(CharacteristicPolynomial(Representation(Frep[6])(FG0PCphi(h))));
        B:=Coefficients(CharacteristicPolynomial(Representation(Frep[9])(FG0PCphi(h))));
        gcoef:=Evaluate(grancoef,a1,A[3]);
        gcoef:=Evaluate(gcoef,a1bar,B[3]);
        gcoef:=Evaluate(gcoef,a2,A[2]);
        gcoef:=Evaluate(gcoef,a2bar,B[2]);
        gcoef:=Evaluate(gcoef,a3,A[1]);
        gcoef:=Evaluate(gcoef,a3bar,B[1]);
        e +:= gcoef;
    end for;

    if j eq 1 then
        if (v[1] eq 0) and (v[3] eq 0) then
            return  ExactQuotient(e+obar(H,H0,2)*3^(v[2])+obar(H,H0,4)*(-1)^(v[2])+obar(H,H0,8)+obar(H,H0,6)*0^(v[2]), #H);
        else
            return ExactQuotient(e,#H);
        end if;
    else
        return ExactQuotient(e,#H0);
    end if;
end function;

function Fverify()
    print "Verifying ST/ST^o ~ Gal(K/k)...";
    verify_groups(Ftwists,Fphi,FM);
    print "Verifying Jac(C)~E^3 at totally split primes of K...";
    verify_splits(Ftwists,Fec);
    print "Bulding dictionary of STgroups...";
    STdict,STgroups:=STdictionary([<Ftwists,Fphi,Fexp>]);
    print "Verifying isogenies between families with the same ST group...";
    verify_isogenies(Ftwists,Fphi,Fexp,STdict);
    verify_minimal_fields_of_definition(Ftwists);
    return true;
end function;

print "Loaded Ftwists,Fphi,Fexp, use Fverify() to run verification.";
