/////////////////////////////////////////////////////////////////////////
//example4b.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
//(second 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());
K:=NumberField(x^2+1);
Kx<xK>:=PolynomialRing(K);

f:=x^3-162*x^2+x;

delta:=K.1-2;

//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);
Aa,toAa1:=AbsoluteAlgebra(Codomain(mu1));
[ClassGroup(Aa[i]): i in [1..NumberOfComponents(Aa)]];

//We give some hint for local points to help the 2-Selmer computation
p:=LookupPrime((1+K.1)*IntegerRing(K));
import "Algaem/descent.m":TryLocalAddPoint;
_:=LocalSelmerImageOfTwoTorsion(two1,p);
[TryLocalAddPoint(two1,p,x): x in
  [-1382*K.1 + 648,1/4*(616*K.1 + 687),-732*K.1 - 639]];
S1small,mp1:=SelmerGroup(two1);
mu2:=IsogenyMu(two2:Fields:={Component(Aa,i):i in [1..NumberOfComponents(Aa)]});
A2:=Codomain(mu2);
A2toA:=hom<A2->A|[theta*delta]>;
//Note that we do not need the 2-Selmer group of E2. We just need the maps.

//We can just take H1 to be the ambient group taken by SelmerGroup
H1:=Codomain(mp1);
toH1:=mp1;

S1:=sub<H1|[(g@@mp1)@toH1:g in OrderedGenerators(S1small)]>;

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

EL:=BaseChange(E1,L);
twoL:=MultiplicationByMMap(EL,2);

//We feed in a representation of a field that is otherwise hard to find
L2:=NumberField( x^8 + 4*x^7 - 148*x^6 - 544*x^5 +
    10040*x^4 + 27728*x^3 - 347784*x^2 - 500400*x + 4903144);
muL:=IsogenyMu(twoL:Fields:={L,L2});
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
OL:=IntegerRing(L);
p:=LookupPrime(2*OL+(1+L.1)*OL);
_:=LocalSelmerImageOfTwoTorsion(twoL,p);
[TryLocalAddPoint(twoL,p,x): x in [
    1/250*(-2772402*L.1^3 + 810317*L.1^2 - 9387140*L.1 - 3840185),
    1/1250*(26389762*L.1^3 - 38936877*L.1^2 + 91237840*L.1 - 24174015),
    1/125*(-227934*L.1^3 - 2623011*L.1^2 + 537370*L.1 - 7067020),
    -60166*L.1^3 - 144620*L.1^2 - 101746*L.1 + 40893,
    80*L.1^3 + 60*L.1^2 + 120*L.1 - 55 
]];

SL,mpL:=SelmerGroup(twoL);

//Conjugation of L/K is easy enough for Lrel.
Lrelcon:=hom<Lrel->Lrel|[-Lrel.1]>;
//We pick a conjugate of L.1 over L that maps to
//Lrelcon(LtoLrel(L.1)).
L1con:=[a[1]:a in Roots(ext<L|>!MinimalPolynomial(L.1))|
   LtoLrel(a[1])eq Lrelcon(LtoLrel(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]]>;
//Just to be sure, check 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.
//Note that in this section
//E1 :      y_1^2=x^3-162x^2+x over K:=Q(i)
//E2 : (i-2)y_2^2=x^3-162x^2+x over K:=Q(i)
//(i.e., E3 from the article)
//L is the field over which E1 and E2 become isomorphic, i.e. K(sqrt(i-2))

//First, that E2(Q) has an image inside the 2-Selmer group of E1 over K
toH1(A2toA(mu2(E2![0,0]))) in S1;

//Which is not in the Norm of the Selmer group over L
//so this element does not come from E1(Q), i.e. represents non-trivial Sha.
toH1(A2toA(mu2(E2![0,0]))) notin NS;
