/////////////////////////////////////////////////////////////////////////
//example3.m
//
//Nils Bruin, 11/09/2002
//
//Accompanying MAGMA script for article
//
//Nils Bruin, Visualising Sha[2] in Abelian surfaces
//
//This script verifies the computations referred to in Section 7, Example 3
//The script has been tested on Magma V2.9-12

//First load the package for computing Selmer groups
AttachSpec("Algaem/Algaem.spec");
_<x>:=PolynomialRing(Rationals());
//Oddly enough, we need the *number field* Q.
Q:=NumberField(x-1:DoLinearExtension);
K:=Q;
Kx<xK>:=PolynomialRing(K);

f:=x^3-22*x^2+21*x+1;

d:=-1;
a:=1;

//The elliptic curves ...
E1:=PolToEllipticCurve(Kx!f);
E2:=PolToEllipticCurve(Kx!d*Parent(x)!Reverse(Eltseq(Evaluate(f,x+a))));

//Getting the 2-Selmer groups ...
two1:=MultiplicationByMMap(E1,2);
two2:=MultiplicationByMMap(E2,2);
mu1:=IsogenyMu(two1);
A<theta>:=Codomain(mu1);

S1small,mp1:=SelmerGroup(two1);
Aa,toAa1:=AbsoluteAlgebra(Codomain(mu1));
mu2:=IsogenyMu(two2:Fields:={Component(Aa,i):i in [1..NumberOfComponents(Aa)]});
A2:=Codomain(mu2);
A2toA:=hom<A2->A|[d*Evaluate(f,a)/(theta-a)]>;
S2small,mp2:=SelmerGroup(two2);

//And compute a structure in which both Selmer groups fit.
S:={LookupPrime(p[1]): p in Factorisation(5*2*Disc*IntegerRing(BaseRing(E2)))}
  where Disc:=Discriminant(E2);
H1,toH1:=pSelmerGroup(Codomain(mu1),2,S);
S1:=sub<H1|[(g@@mp1)@toH1:g in OrderedGenerators(S1small)]>;
S2:=sub<H1|[(g@@mp2)@A2toA@toH1:g in OrderedGenerators(S2small)]>;

Slower:=sub<H1|OrderedGenerators(S1) cat OrderedGenerators(S2)>;
Supper:=sub<H1|OrderedGenerators(S1 meet S2)>;

//The hyperelliptic curve
C:=HyperellipticCurve(Evaluate(Parent(x)!f,d*x^2+a));
//And the algebra used to represent the fake-selmer group
AC<Theta>:=quo<Parent(x)|Evaluate(Parent(x)!f,d*x^2+a)>;

//We need a little function that is used to transform the curve internally
//in TwoSelmerGroupData. Unfortunately, the information returned about the fake
//selmer group is not transformed back, so we do this ourselves.
function ModelForTSG(C)
    g := Genus(C);
    // First simplify
    C1, m := SimplifiedModel(C);
    // Then make integral
    C1, m1 := IntegralModel(C1 : Reduce := true);
    m := m*m1;
    // Now, we have to do different things according to the parity of deg(f)
    f := HyperellipticPolynomials(C1);
    if IsOdd(Degree(f)) then
	a := LeadingCoefficient(f);
	C1, m1 := Transformation(C1, [a,0,0,1], a^g);
	m := m*m1;
    else // IsEven(Degree(f))
	a := Abs(LeadingCoefficient(f))/
             GCD(ChangeUniverse(Coefficients(f),Integers()));
	C1, m1 := Transformation(C1, [a,0,0,1], a^(g+1));
	m := m*m1;
    end if;
    return C1, m;
end function;

//so TwoSelmerGroupData will actually work with the model below.
Ctsg,toCtsg:=ModelForTSG(C);
Jtsg:=Jacobian(Ctsg);

//Compute the fake Selmer group (remove Bound:=300 to get the same, but
//guaranteed to be correct answer. Takes quite a while longer.)
_,SCrank,l:=TwoSelmerGroupData(Jtsg:Bound:=300);
//Now we have to figure out how to transform the fake Selmer group back
//to the model C.
br:=BaseRing(C);
dp1,dp2:=HyperellipticPolynomials(C);
brx:=Parent(dp1);
brxy:=quo<brxy|brxy.1^2+dp2*brxy.1-dp1> where brxy:=PolynomialRing(brx);
//Below is a generic point on Ctsg, expressed in a generic point 
//on C. [brx.1,brxy.1]
genpt:=toCtsg(C(brxy)![brx.1,brxy.1]);
//this had better still be an affine point
assert genpt[3] eq 1;
//If so, then we can easily compute the relation between Theta and the "Theta"
//for Ctsg.
Thetatsg:=Evaluate(Evaluate(genpt[1],1),Theta);
//and we can obtain generators of the fake selmer group of C
CS:=[Evaluate(p,Thetatsg):p in l[1]];
CSseq:=<CS,l[2],l[3]>;

//And compute the norm of the fake selmer group
Ax<xA>:=PolynomialRing(A);
_<ThetaR>:=quo<Ax|xA^2-(theta-a)/d>;
NS:=sub<H1|[Norm(Evaluate(p,ThetaR))@toH1: p in CS]>;

///////////////////////////////////////////////////////////////////////////
//Given the computations above, we now have computed the objects with which we
//can check the claims in Example 3 in Section 7 of the article.

//That the 2-Selmer group of E1 is (Z/2Z)^4:
IsIsomorphic(S1,AbelianGroup([2,2,2,2]));

//That the 2-Selmer group of E2 is (Z/2Z)^3:
IsIsomorphic(S2,AbelianGroup([2,2,2]));

//That the 2-Selmer group of E2 is generated by rational points on E2
sub<H1|[toH1(A2toA(mu2(P))):P in [E2![1,1],E2![2,7],E2![5,23]]]> eq S2;

//That the 2-Selmer rank of Jac_C is 5:
SCrank eq 5;

//In fact, the 2-rank of the fake-selmer group is 4 (see doc. why this is
//the correct way of computing the size of the fake selmer group)
#CSseq[1]-#CSseq[2]+#CSseq[3] eq 4;

//Two independent elements that are in S2 (and therefore come from E2(Q))
toH1(-theta*(1-2*theta)) in S2;
toH1(-theta*(1-8*theta)) in S2;

//but are not in the Norm of the Selmer group of Jac(C)
NS eq sub<H1|[toH1(21-theta)]>;
