/////////////////////////////////////////////////////////////////////////
//example4a.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 4
//(first part)
//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-162*x^2+x;

delta:=-1;

//The elliptic curves ...
E1:=PolToEllipticCurve(Kx!f);
E2:=PolToEllipticCurve(Kx!f*delta);

//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|[theta*delta]>;
S2small,mp2:=SelmerGroup(two2);

//And compute a structure in which both Selmer groups fit.
S:={LookupPrime(p[1]): p in Factorisation(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)>;

//Some stuff necessary to be able to compute the Norm map later on
Lrel:=NumberField(xK^2-delta);
L:=NumberField(x^2-delta);
LisLrel,LtoLrel:=IsIsomorphic(L,Lrel);
assert LisLrel;

EL:=BaseChange(E1,L);
twoL:=MultiplicationByMMap(EL,2);
muL:=IsogenyMu(twoL);
AL<thetaL>:=Codomain(muL);
AaL:=AbsoluteAlgebra(AL);
//We're lazy with the classgroup computation here. If you want, you can remove
//the "Bound:=300" and still arrive at the same (but now proven-to-be-correct)
//answer.
[ClassGroup(AaL[i]:Bound:=300):
    i in [1..NumberOfComponents(AaL)]];
//We help the package a bit by providing some auxiliary points.
//This saves SelmerGroup a lot of searching
p:=LookupPrime((1+L.1)*IntegerRing(L));
import "Algaem/descent.m":TryLocalAddPoint;
_:=LocalSelmerImageOfTwoTorsion(twoL,p);
[TryLocalAddPoint(twoL,p,x): x in
  [-1382*L.1 + 648,1/4*(616*L.1 + 687),-732*L.1 - 639]];

SL,mpL:=SelmerGroup(twoL);

//Unfortunately, @LtoLrel@@LtoLrel isn't always the identity.
//Conjugation of L/K is easy enough for Lrel.
Lrelcon:=hom<Lrel->Lrel|[-Lrel.1]>;
//To compute what this does to L.1, we consider two options and pick the
//nontrivial one
L1con:=[a:a in [LtoLrel(L.1)@@LtoLrel,Lrelcon(LtoLrel(L.1))@@LtoLrel]| a ne L.1];
//We should really end up with only one option
assert #L1con eq 1;
//That supposedly defines conjugation
Lcon:=hom<L->L|[L1con[1]]>;
//But better check that it really does, i.e., leaves K invariant.
assert k eq Lcon(k) where k:=(Lrel!K.1)@@LtoLrel;
//And, just to be sure, that conjugation is of order 2 (it is not the identity
//by picking L1con different from L.1)
assert Lcon(Lcon(L.1)) eq L.1;

ALcon:=hom<AL->AL|a:->AL![Lcon(b):b in Eltseq(a)]>;
ALnorm:=hom<AL->A|a:->A![Eltseq(LtoLrel(b))[1]:b in Eltseq(a*ALcon(a))]>;
NS:=sub<H1|[ALnorm(g@@mpL)@toH1:g in OrderedGenerators(SL)]>;

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

//Check that the Selmer group of E1 over Q is (Z/2Z)^3:
IsIsomorphic(S1,AbelianGroup([2,2,2]));

//Check that the Selmer group of E2 over Q is (Z/2Z)^2 and is generated by
//E2(Q):
IsIsomorphic(S2,AbelianGroup([2,2]));
S2 eq sub<H1|[toH1(A2toA(mu2(P))):P in [E2![0,0],E2![9/4,231/8]]]>;

//Check that E1(Q) does NOT land in S2, i.e.
//that E2 DOES visualise Sha[2](E1)
toH1(mu1(E1![0,0])) notin S2;

//However, the norm of S(E/Q(i)) equals S1 intersect S2, so gives no new info
NS eq S1 meet S2;
