// Klein Twists

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

Q:=Rationals();
RT<T>:=PolynomialRing(Q);
Rx<x,y,z>:=PolynomialRing(Q,3);
KM<a>:=NumberField(T^2+T+2);
RX<X,Y,Z>:=PolynomialRing(KM,3);
Q2<d>:=NumberField(T^2-2);
RU<U,V,W>:=PolynomialRing(Q2,3);

G<s,t,u,zz>:=FreeGroup(4);
KG,Kpi:=quo<G|s^3, t^2, u^7, (s*t)^2, (s^3*t*u)^3, (s^2*t^2*u)^3, (s^2*t*u^2)^3, (s*t^2*u)^3, zz^2, zz*s*zz*s^2, zz*t*zz*t, zz*u*zz*u>;
KG0:=sub<KG|s,t,u>;
KGP,Kmu:=PermutationGroup(KG);
Kphi:=map<G->KGP|x:->Kmu(Kpi(x))>;
KG0P := sub<KGP|Kphi(s),Kphi(t),Kphi(u)>;
assert IdentifyGroup(KGP) eq <336,208>;
assert IdentifyGroup(KG0P) eq <168,42>;

// Over K the Jacobian of every twist is isogenous to the cube of the elliptic curve below
Kec := EllipticCurve("49a3");

// Data format is <n, H0 generators, h, ID(H), ID(H0), f(x,y,z), abs(disc(f)), g(t), h(t), ST, ST_k(a))>
//   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

Ktwists:=<
<1, [Id(G)], Id(G), <2,1>, <1,1>, x^4+6*x^3*z-3*x^2*y^2+3*x^2*y*z-3*x^2*z^2+6*x*y^3+3*x*y^2*z+3*x*y*z^2+y^4-3*y^2*z^2+6*y*z^3+z^4, 7^34, T, T^2+7, "">,
<2, [t], Id(G), <4,2>, <2,1>, 3*x^4+28*x^3*y+105*x^2*y^2-21*x^2*z^2+196*x*y^3+147*y^4+147*y^2*z^2-49*z^4, 2^36*7^43, T, T^4-3*T^2+4, "">,
<12, [u*s*t*u^6,s*u*t*u^6*s^2], zz, <4,2>, <4,2>, X^4+9*a*X^2*Y^2+6*a*X^2*Z^2+9*Y^4+18*a*Y^2*Z^2+4*Z^4, 2^144*3^36*7^14, T^2+7, T^8+12*T^6+65*T^4+48*T^2+16, "">,
<3, [s], t, <6,1>, <3,1>, x^4+3*x^3*y-9*x^3*z+9*x^2*y^2-6*x^2*z^2+18*x*y^3+3*x*y^2*z-3*x*y*z^2+y^4+4*y^3*z-3*y^2*z^2+7*y*z^3, 5^18*7^25, T, T^3+T^2+2*T+3, "2450ba1">,
//<3, [s], t, <6,1>, <3,1>, x^4-6*x^3*y+7*x^3*z+3*x^2*y^2+21*x^2*y*z+3*x^2*z^2+18*x*y^3-9*x*y*z^2+9*y^4+21*y^3*z-9*y^2*z^2-3*z^4, 3^36*7^25, T, T^3-3*T-5>,
<4, [s], Id(G), <6,2>, <3,1>, x^3*y+y^3*z+z^3*x, 7^7, T, T^6+T^5+T^4+T^3+T^2+T+1, "">,
<6, [u^2*t*u^3*t*u^2], u^5*t*u^2, <8,1>, <4,1>, x^4+3*x^2*y^2-3*x^2*z^2+2*y^3*z+3*y^2*z^2+2*y*z^3, 2^36*7^7, T^2+1, T^8+2*T^7-14*T^4+16*T+4, "64a4">,
//<6, [u^2*t*u^3*t*u^2], u^5*t*u^2, <8,1>, <4,1>, 2*x^4-3*x^3*y-12*x^2*y^2+3*x^2*y*z-3*x^2*z^2-11*x*y^3+15*x*y^2*z-15*x*y*z^2-7*y^4+y^3*z-2*y*z^3+z^4, 7^7*23^18, T^2+23, T^8-4*T^7+14*T^5-28*T^3-7*T^2+24*T-145>,
//<6, [u^2*t*u^3*t*u^2], u^5*t*u^2, <8,1>, <4,1>, U^4-6*U^2*V^2+(24*d+12)*U*V*W^2-7*V^4+(4*d+9)*W^4, 2^144*7^32, T^4+1, T^8+7>,
<7, [u^2*t*u^3*t*u^2], Id(G), <8,3>, <4,1>, 12*x^4-80*x^3*y+60*x^2*y^2-24*x^2*z^2-104*x*y^3+24*x*y*z^2+83*y^4+36*y^2*z^2-2*z^4, 2^99*7^25, T, T^4-4*T^2-14, "">,
<13, [s*u,t*u], zz, <12,3>, <12,3>, 3*X^4+(-18*a+12)*X^3*Y+(12*a+4)*X^3*Z+(-27*a+36)*X^2*Y^2+(9*a+6)*X^2*Z^2+36*X*Y^2*Z+27*Y^4+(54*a-36)*Y^3*Z+(-54*a+36)*Y^2*Z^2+(18*a-12)*Y*Z^3+(-3*a+2)*Z^4, 2^108*3^90*7^14, T^2+7, T^6-147*T^2+343, "">,
<8, [s,t], Id(G), <12,4>, <6,1>, 7*x^3*z+3*x^2*y^2-6*x*y*z^2+2*y^3*z-4*z^4, 2^18*3^27*7^7, T, T^6+3*T^5+2*T^4-T^3+4*T^2+5*T+1 /*T^3-T^2+5*T+1*/, "36a1">,
//<8, [s,t], Id(G), <12,4>, <6,1>, 1587*x^4-1449*x^2*y^2+4347*x^2*y*z-483*x^2*z^2+3381*x*y^3-3381*x*y*z^2+3381*x*z^3-441*y^4+2646*y^3*z-4263*y^2*z^2+882*y*z^3-49*z^4, 3^36*7^43*23^36, T, T^6+T^5+2*T^4+5*T^3+6*T^2+8>,
<5, [u], Id(G), <14,1>, <7,1>, x^3*y-21*x^2*z^2+x*y^3-42*x*y*z^2-147*x*z^3+2*y^4+21*y^3*z+63*y^2*z^2-196*z^4, 7^25, T, T^7+7*T^3-7*T^2+7*T+1, "784h1">,
<10, [u^2*t*u^3*t*u^2, u^5*t*u^2], Id(G), <16,7>, <8,3>, x^4+3*x^2*y^2-3*x^2*z^2+2*y^3*z+3*y^2*z^2+2*y*z^3, 2^36*7^7, T, T^8+2*T^7-14*T^4+16*T+4, "64a4">,
//<10, [u^2*t*u^3*t*u^2, u^5*t*u^2], Id(G), <16,7>, <8,3>, 2*x^4-3*x^3*y-12*x^2*y^2+3*x^2*y*z-3*x^2*z^2-11*x*y^3+15*x*y^2*z-15*x*y*z^2-7*y^4+y^3*z-2*y*z^3+z^4, 2^7*23^18, T, T^8-4*T^7+14*T^5-28*T^3-7*T^2+24*T-145>,
//<10, [u^2*t*u^3*t*u^2, u^5*t*u^2], Id(G), <16,7>, <8,3>,  U^4-6*U^2*V^2+(24*d+12)*U*V*W^2-7*V^4+(4*d+9)*W^4, 2^144*7^32, T^2-2, T^8+7>,
<14, [s*u*s*t,s*u^6*s^2*t*u^2], zz, <24,12>, <24,12>, (3*a-2)*X^4+(30*a-20)*X^3*Y+(90*a-60)*X^2*Y^2+(9*a+6)*X^2*Z^2+(150*a-100)*X*Y^3+60*X*Y^2*Z+(12*a+4)*X*Z^3+(150*a-25)*Y^4+(-45*a+60)*Y^2*Z^2+(30*a-20)*Y*Z^3+3*Z^4, 2^108*3^54*5^36*7^14, T^2+7, T^8+16*T^6+10*T^5+45*T^4+10*T^3+146*T^2-60*T+134 /* T^4 + 2*T^3 + 6*T^2 - 6 */, "">,
<9, [u,s], Id(G), <42,1>, <21,1>, 2*x^3*y+y^3*z+z^3*x, 2^9*7^7, T, T^7-2, "">,
//<11, [t,u,s], Id(G), <336,208>, <168,42>, 4*x^3*y-2*x^3*z+3*x^2*y^2+x*y^3-2*x*z^3-2*y^3*z+3*y^2*z^2-2*y*z^3, 2^18*5^9*7^7, T, T^8-4*T^7-35*T^4-50*T+25>>;
<11, [t,u,s], Id(G), <336,208>, <168,42>, 2*x^3*y-2*x^3*z-3*x^2*z^2-2*x*y^3-2*x*z^3-4*y^3*z+3*y^2*z^2-y*z^3, 2^18*7^25, T, T^8+4*T^7+21*T^4+18*T+9, "144b1">>;

// 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 Kexp(G0,H,H0,j,v)
    assert j in [1,2];
    assert #v eq 3;

    KK<s7>:=QuadraticField(-7);
    R<u,ubar,a1,a2,a3,a1bar,a2bar,a3bar>:=PolynomialRing(KK,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);

    if IsDivisibleBy(#H0,7) then g7 := [g:g in H0|Order(g) eq 7][1]; end if;

    e := 0;
    for h in H0 do
        if Order(h) eq 1 then A:=[-1,3,-3,1]; B:=A; end if;
        if Order(h) eq 2 then A:=[-1,-1,1,1]; B:=A; end if;
        if Order(h) eq 3 then A:=[-1,0,0,1]; B:=A; end if;
        if Order(h) eq 4 then A:=[-1,1,-1,1]; B:=A; end if;
        if Order(h) eq 7 and IsConjugate(G0,h,g7) then A:=[-1,-1/2+s7/2,1/2+s7/2,1]; B:=[-1,-1/2-s7/2,1/2-s7/2,1]; end if;
        if Order(h) eq 7 and IsConjugate(G0,h,g7^3) then A:=[-1,-1/2-s7/2,1/2-s7/2,1]; B:=[-1,-1/2+s7/2,1/2+s7/2,1]; end if;

        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:=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 Kverify()
    print "Verifying ST/ST^o ~ Gal(K/k)...";
    verify_groups(Ktwists,Kphi,KM);
    print "Verifying Jac(C)~E^3 at totally split primes of K...";
    verify_splits(Ktwists,Kec);
    verify_minimal_fields_of_definition(Ktwists);
    return true;
end function;

print "Loaded Ktwists,Kphi,Kexp, use Kverify() to run verification.";

RAB<A,B>:=PolynomialRing(Rx,2);
HKpoly:=A*x^4+7*B*x^3*z+3*x^2*y^2-3*A^2*x^2*z^2-6*B*x*y*z^2-5*A*B*x*z^3+2*y^3*z+3*A*y^2*z^2+2*A^2*y*z^3-4*B^2*z^4;

    