#!/bin/sh ############################################################################### ############################################################################### ## ## thesis.sh ## ## Nils Bruin ## ## This archive consists of a concatanation of ASCII files. On a UNIX system, ## put this file in the directory you want to unpack it in and type ## ## "sh thesis.sh" ## ## If you unpack it on a non-UNIX system, delete the lines of the form ## "cat ... # UNIX extraction command" ## and save the lines between the designators ## "#--- (start) ---#" and ## "#--- (end) ---#" ## in the file . ## ############################################################################### ############################################################################### #!/bin/sh # this archive contains: # -rw-r--r-- 1 nbruin staff 4400 Jun 21 16:42 array.g # -rw-r--r-- 1 nbruin staff 2336 Jun 21 16:42 domain.g # -rw-r--r-- 1 nbruin staff 27776 Jun 22 15:59 ellcalc.g # -rw-r--r-- 1 nbruin staff 854 Jun 21 16:42 eq238E5.g # -rw-r--r-- 1 nbruin staff 2824 Jun 21 16:42 eq238E7.g # -rw-r--r-- 1 nbruin staff 2959 Jun 21 16:42 eq238E9.g # -rw-r--r-- 1 nbruin staff 892 Jun 21 16:42 eq245E3.g # -rw-r--r-- 1 nbruin staff 1229 Jun 21 16:42 eq245E6.g # -rw-r--r-- 1 nbruin staff 1095 Jun 21 16:42 eq245E7.g # -rw-r--r-- 1 nbruin staff 945 Jun 21 16:42 eq245E8.g # -rw-r--r-- 1 nbruin staff 1337 Jun 21 16:42 eq245E9.g # -rw-r--r-- 1 nbruin staff 777 Jun 21 16:42 eq283E2.g # -rw-r--r-- 1 nbruin staff 777 Jun 21 16:42 eq283E4.g # -rw-r--r-- 1 nbruin staff 890 Jun 21 16:42 eq283E7.g # -rw-r--r-- 1 nbruin staff 4997 Jun 21 17:23 index.txt # -rw-r--r-- 1 nbruin staff 14095 Jun 22 15:33 isogdesc.g # -rw-r--r-- 1 nbruin staff 2864 Jun 21 16:49 latcalc.g # -rw-r--r-- 1 nbruin staff 16491 Jun 21 16:49 loccalc.g # -rw-r--r-- 1 nbruin staff 4332 Jun 22 13:53 matalg.g # -rw-r--r-- 1 nbruin staff 3341 Jun 21 16:47 powseries.mpl # It is produced by a very simple shar that does not protect special # symbols, but does minimise work for manual extraction. # If possible, treat this data as binary. # #--- index.txt (start) ---# cat << "#--- index.txt (end) ---#" > index.txt # UNIX extraction command ############################################################################## ############################################################################## ## ## index.txt ## ## Nils Bruin, 21 june, 1999 ## ## this file describes the contents of the files in this archive. ## ############################################################################## ############################################################################## The files in this archive perform calculations used in the following publications. [thesis] Nils Bruin, Chabauty Methods and Covering Techniques applied to Generalised Fermat equations. PhD-thesis, Universiteit Leiden, 1999 [MI 99-14] Nils Bruin, Chabauty methods using elliptic curves. Report MI 99-14, Universiteit Leiden, 1999 [MI 99-15] Nils Bruin, Chabauty methods using covers on curves of genus 2. Report MI 99-15, Universiteit Leiden, 1999 ############################################################################## In order to use these files: 1) Make sure that you have access to KASH (these scripts were tested on KASH 2.0, but probably work (with possibly small modifications) on a far wider range of KASH versions). KASH can be obtained in executable form for a variety of platforms free of charge from http://www.math.tu-berlin.de/~kant/kash.html ftp://ftp.math.tu-berlin.de/pub/algebra/Kant/Kash/ 2) Make sure that all files from this archive are in the current directory and start KASH. 3) Type (for example) ``Read("eq283E2.g");''. You might want to type ``LogTo("eq283E2.log");'' have a file with the session transcript to examine later. The files do not require further input. In order to understand the output, you might want to examine the scripts themselves. It has been a lot of work to write these programs. Feel free to use them for scientific purposes, but please contact me on how to include proper references if you obtain publishable results by them. ############################################################################## The following files refer to curves related to x^2+y^8=z^3 ([thesis, Section 4.6]). eq283E2.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E2 eq283E4.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E4 eq283E7.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E7 ############################################################################## The following files refer to curves related to x^8+y^3=z^2 ([thesis, Section 4.7]). eq283E5.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E5 eq283E7.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E7 eq283E9.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E9 ############################################################################## The following files refer to curves related to x^2+-y^4=z^5 ([thesis, Section 4.8] and [MI 99-14, Section 7]). eq283E3.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E3 eq283E6.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E6 eq283E7.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E7 eq283E8.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E8 eq283E9.g: Does 2-isogeny descent and Elliptic curve Chabauty on curve E9 ############################################################################## The following file is a Maple V-script and computes the power series mentioned in [thesis, section 4.4] and [MI 99-14, section 3]. To see the results, make sure you have access to Maple V (Waterloo Software) and type "read`powseries.mpl`;" in Maple. powseries.mpl: Computes the power series for Theta. ############################################################################## The following files are used by the eq*.g files described above. These are libraries of functions. You might want to experiment with them to test their correctness. At the moment, I do not intend to make a software package out of them, so a user manual is not provided. The files eq*.g should give enough examples to get started (and, of course, the code itself documents exactly what it does :-)). array.g: Implements arbitrary dimensional arrays in KASH domain.g: Some functions to facilitate writing generic code for Z, number fields and finite fields. matalg.g: Some functions for calculating with matrices (especially linear subspaces spanned by the row vectors of a matrix) latcalc.g: Some routines to do calculations with sublattices of ZZ^r (has nothing to do with standard lattice functions of KASH) loccalc.g: Routines to do computations in local fields (with finite precision) ellcalc.g: Routines to do computations on elliptic curves, including EC Chabauty for elliptic degree 2 covers. isogdesc.g: 2-isogeny descent for elliptic curves over number fields. ############################################################################## # end of index.txt #--- index.txt (end) ---# #--- array.g (start) ---# cat << "#--- array.g (end) ---#" > array.g # UNIX extraction command ############################################################################### ############################################################################### ## ## array.g ## ## Nils Bruin, 5 april 1999 ## ## Some routines to work with arrays of arbitrary dimensions. ## ############################################################################### ############################################################################### Tail:=L->L{[2..Length(L)]}; #(returns the "tail" of a list: all but the first element) ############################################################################### ArrInit:=function(dim,L) # initialises an "array" structure: a 0-base indexed list-of-lists with # fixed depth. local R,l; if dim<1 then Error("We don't do arrays of nonpositive dimension"); elif dim=1 then return rec(dim:=1, L:=L); else R:=rec(dim:=dim,L:=[]); for l in L do Add(R.L,ArrInit(dim-1,l).L); od; fi; return R; end; ############################################################################### ArrSetElt:=function(arr,v,val) # sets element indexed by v to val. if Length(v)<>arr.dim then Error("IndexVector dimension mismatch"); fi; if arr.dim=1 then arr.L[v[1]+1]:=val; else ArrSetElt(rec(dim:=arr.dim-1,L:=arr.L[v[1]+1]),Tail(v)); fi; end; ############################################################################### ArrElt:=function(arr,v) # returns entry indexed by v. local R,i; if Length(v)<>arr.dim then Error("Index vector dimension mismatch"); fi; if v[1]<0 or v[1]>=R.dims[1] then Error("Index out of range"); fi; R:=arr.L; for i in [1..arr.dim] do R:=R[v[i]+1]; od; return R; end; ############################################################################### ListAugFront:=function(v,L) # returns a list with v added in front of L. It is undefined if entries # from L are duplicated or identically appended. local R; R:=[v]; Append(R,L); return R; end; ############################################################################### ArrLookup:=function(arr,val) # returns index of val in arr or false, if it does not occur. If val occurs # more than once, the lexicographically smallest index is returned (i.e. the # index that occurs first if arr is written out as list of lists). local i, t, v; if arr.dim=1 then for i in [1..Length(arr.L)] do if arr.L[i]=val then return [i-1]; fi; od; else for i in [1..Length(arr.L)] do t:=ArrLookup(rec(dim:=arr.dim-1,L:=arr.L[i]),val); if t<>false then return ListAugFront(i-1,t); fi; od; fi; return false; end; ############################################################################### ArrLookupAll:=function(arr,val) # returns a list of all indices of val in arr. local L,i,t; L:=[]; if arr.dim=1 then for i in [1..Length(arr.L)] do if arr.L[i]=val then Add(L,[i-1]); fi; od; else for i in [1..Length(arr.L)] do t:=ArrLookupAll(rec(dim:=arr.dim-1,L:=arr.L[i]),val); Append(L,List(t,u->ListAugFront(i-1,u))); od; fi; return L; end; ############################################################################### ArrApply:=function(arr,fun) # returns a new array obtained by applying fun to the entries of arr. if arr.dim=1 then return rec(dim:=1,L:=List(arr.L,fun)); else return rec(dim:=arr.dim,L:=List(arr.L, l->ArrApply(rec(dim:=arr.dim-1,L:=l),fun).L)); fi; end; ############################################################################### ArrToList:=arr->arr.L; # convert array to a list structure (well, actually it already is. We only # have to strip off the record structure). We lose dimension information # this way. If entries of the array are lists themselves, it is not so easy to # recapture this info from the list structure. ############################################################################### ArrAddRow:=function(arr1,arr2) # if dimension of arr1=1 then add arr2 as last entry, otherwise arr2 should # be an array of one dimension lower and the entries of arr2 are added as an # extra row to arr1 (the are added identically; NOT copied, so the effect is # ciomparable of "Add"ing a list to another list) if arr1.dim=1 then Add(arr1.L,arr2); elif arr1.dim=arr2.dim+1 then Add(arr1.L,arr2.L); else Error("Dimension mismatch"); fi; end; #--- array.g (end) ---# #--- domain.g (start) ---# cat << "#--- domain.g (end) ---#" > domain.g # UNIX extraction command ############################################################################### ############################################################################### ## ## domain.g ## ## Nils Bruin, 4 april 1999 ## ## Some routines to facilitate writing routines for generic domains. ## Z, Q, Orders and finite fields are mainly supported now. ## Main reason is to get representatives for 0 and 1 right. ## ############################################################################### ############################################################################### # Some missing type tests. IsFFElt:=x->TYPE(x)="KANT finite field elt"; IsFF:=x->TYPE(x)="KANT finite field"; ############################################################################### EltRing:=function(x) # returns ring of definition. Note that Z is NOT an order according to # KASH. Thus, Z and Q are handled differently. Perhaps Z should be returned # for rational numbers as well. This is not so critical since KASH automatically # casts between integers and rationals. (Not so for orders !) if IsInt(x) then return Z; elif IsRat(x) then return Q; elif IsElt(x) then return EltOrder(x); elif IsFFElt(x) then return FFEltFF(x); else Error("Ring type not recognised"); fi; end; ############################################################################### RingZero:=function(R) # returns neutral additive element of R. if IsFF(R) then return 0*FFGenerator(R); elif IsOrder(R) then return Elt(R,0); else return 0; fi; end; ############################################################################### RingOne :=function(R) # returns neutral multiplicative element of R. if IsFF(R) then return FFGenerator(R)^0; elif IsOrder(R) then return Elt(R,1); else return 1; fi; end; ############################################################################### RingForce:=function(R,x) # forces x in R. Only works if KASH can figure out the cast itself. # Not usable for lifting or descending from an order to Z. In those cases, # no error is generated, but just the original element is returned. This routine # should therefore not be considered user interface. Either put up with it or # don't use it. Do not delete, as other routines obviously will use it. return x+RingZero(R); end; #--- domain.g (end) ---# #--- ellcalc.g (start) ---# cat << "#--- ellcalc.g (end) ---#" > ellcalc.g # UNIX extraction command ############################################################################### ############################################################################### ## ## ellcalc.g ## ## Nils Bruin, 4 april 1999 ## ## Various routines to do calculations on elliptic curves of the form ## gamma*Y^2=X^3+a2*X^2+a4*X+a6 ## over number fields and finite fields and for computations ## at localisations of primes of good reduction. ## ############################################################################### ############################################################################### Read("domain.g"); Read("array.g"); Read("loccalc.g"); Read("latcalc.g"); ############################################################################### EllDisc:=function(ec) # discriminant of curve (gamma should be included in some way, so it is # linearly) return ec.gamma*(16*ec.a2^2*ec.a4^2-64*ec.a4^3-64*ec.a2^3*ec.a6+ 288*ec.a2*ec.a4*ec.a6-432*ec.a6^2); end; ############################################################################### EllJay:=function(ec) # jay invariant of curve (gamma is unimportant for this, as it is a geometric # invariant. return 4096*ec.gamma*(ec.a2^2-3*ec.a4)^3/EllDisc(ec); end; ############################################################################### EllInit:=function(gamma,a1,a2,a3,a4,a6) # initialise elliptic curve data structure. (other fields may be added later) local ec,R; R:=EltRing(gamma+a1+a2+a3+a4+a6); if RingForce(R,a1)<>RingZero(R) or RingForce(R,a3)<>RingZero(R) then Error("Can only handle a1=a3=0"); fi; ec:=rec(a2:=RingForce(R,a2),a4:=RingForce(R,a4), a6:=RingForce(R,a6),gamma:=RingForce(R,gamma),R:=R); if EllDisc(ec)=RingZero(R) then Error("Singular curve"); fi; return ec; end; ############################################################################### EllPointOnCurve:=function(ec,G) # test whether a point lies on the curve. return ec.gamma*G[2]^2*G[3]-(G[1]^3+ec.a2*G[1]^2*G[3]+ ec.a4*G[1]*G[3]^2+ec.a6*G[3]^3)=RingZero(ec.R); end; ############################################################################### EllPointNormal:=function(ec,G) # puts a projective point in normal form. (ec is not really important. In fact, # only the "R" field is used to force the entries of G in the right ring. # Normal form does NOT mean that all entries are integral if ec is over an # order or Z. G:=List(G,u->RingForce(ec.R,u)); if G[3]<>RingZero(ec.R) then G:=[G[1]/G[3],G[2]/G[3],RingOne(ec.R)]; elif G[2]<>RingZero(ec.R) then G:=[G[1]/G[2],RingOne(ec.R),RingZero(ec.R)]; elif G[1]<>RingZero(ec.R) then G:=[RingOne(ec.R),RingZero(ec.R),RingZero(ec.R)]; else Error("Point is not projective"); fi; return G; end; ############################################################################### #neutral element of elliptic curve. EllZero:=ec->[RingZero(ec.R),RingOne(ec.R),RingZero(ec.R)]; ############################################################################### EllNeg:=function(ec,G) # returns inverse of a point. if not(EllPointOnCurve(ec,G)) then Error("Point not on curve"); fi; return EllPointNormal(ec,[G[1],-G[2],G[3]]); end; ############################################################################### EllAdd:=function(ec,G1,G2) # returns the sum of two points # The formulas are just a matter of calculating. if not(EllPointOnCurve(ec,G1)and EllPointOnCurve(ec,G2)) then Error("Points not on curve"); fi; G1:=EllPointNormal(ec,G1);G2:=EllPointNormal(ec,G2); if G1=EllZero(ec) then return G2; elif G2=EllZero(ec) then return G1; elif G1[1]=G2[1] then if G1[2]=-G2[2] then return EllZero(ec); else return EllPointNormal(ec, [2*(8*ec.a6*G1[1]+2*G1[1]^2*ec.a4-ec.a4^2-G1[1]^4+4*ec.a6*ec.a2)* (ec.a6+G1[1]^3+ec.a2*G1[1]^2+ec.a4*G1[1]), -G1[2]*(5*G1[1]^4*ec.a4+20*ec.a6*G1[1]^3-8*ec.a6^2-ec.a4^3+ 4*ec.a6*ec.a2*ec.a4-5*ec.a4^2*G1[1]^2+G1[1]^6-4*ec.a6*ec.a4*G1[1]+ 20*ec.a6*ec.a2*G1[1]^2-2*ec.a2*G1[1]*ec.a4^2+8*ec.a2^2*G1[1]*ec.a6+ 2*ec.a2*G1[1]^5), -8*(ec.a6+G1[1]^3+ec.a2*G1[1]^2+ec.a4*G1[1])^2]); fi; else return EllPointNormal(ec, [(G1[1]-G2[1])*(G1[1]^2*G2[1]+ec.a4*G1[1]+2*ec.a2*G1[1]*G2[1]+ G1[1]*G2[1]^2-2*ec.gamma*G1[2]*G2[2]+2*ec.a6+ec.a4*G2[1]), -2*G1[2]*ec.a2*G1[1]*G2[1]-G1[2]*G2[1]^3-4*G1[2]*ec.a6- 3*G1[2]*G1[1]*G2[1]^2-G1[2]*ec.a4*G1[1]-3*G1[2]*ec.a4*G2[1]- 2*G1[2]*ec.a2*G2[1]^2+G2[2]*G1[1]^3+3*G2[2]*G1[1]^2*G2[1]+ 2*G2[2]*ec.a2*G1[1]*G2[1]+4*G2[2]*ec.a6+3*G2[2]*ec.a4*G1[1]+ G2[2]*ec.a4*G2[1]+2*G2[2]*ec.a2*G1[1]^2, (G1[1]-G2[1])^3]); fi; end; ############################################################################### EllMul:=function(ec,n,G) # returns n*G. local Gt; if n=0 then return EllZero(ec); elif n<0 then return EllNeg(ec,EllMul(ec,-n,G)); elif n mod 2 = 0 then Gt:=EllMul(ec,n/2,G); return EllAdd(ec,Gt,Gt); else Gt:=EllMul(ec,(n-1)/2,G); return EllAdd(ec,EllAdd(ec,Gt,Gt),G); fi; end; ############################################################################### EllGenInit:=function(ec,Gs,Rk) # routine to register "generators" in ec. For most purposes, these should # be Mordell-Weill pseudo-generators. Rk should reflect that the # first Rk entries of Gs are "free" and that the others are torsion. # This is not checked. "Rk" does not make sense for curves over finite fields, # but EllRed does inherit it. if IsBound(ec.Gs) then Error("Gs already defined"); fi; if not(ForAll(Gs,u->EllPointOnCurve(ec,u))) then Error("One of points not on curve"); else ec.Gs:=List(Gs,u->EllPointNormal(ec,u)); ec.Rk:=Rk; fi; end; ############################################################################### EllVecToPnt:=function(ec,v) # calculates the specified linear combination of the preregistered Gs. local G,i; if Length(v)<>Length(ec.Gs) then Error("Vector of wrong dimension"); fi; G:=EllZero(ec); for i in [1..Length(v)] do G:=EllAdd(ec,G,EllMul(ec,v[i],ec.Gs[i])); od; return G; end; ############################################################################### EllYsqr:=function(ec,xc) # Given an x-coordinate, returns y^2 belonging to it. Useful for checking # if an x-coordinate belongs to a rational point. (over whatever ring) return (xc^3+ec.a2*xc^2+ec.a4*xc+ec.a6)/ec.gamma; end; ############################################################################### EllPhiInit:=function(ec,c1,c2) # initialises a degree 2 cover E -> P_1. If phi is even, then the ill-defined # point (G_phi) should be [0:1:0]. local Gphi; if IsBound(ec.c) then Error("Phi already initialised"); fi; Gphi:=EllPointNormal(ec,[c1[2]*c2[3]-c1[3]*c2[2], c1[3]*c2[1]-c1[1]*c2[3], c1[1]*c2[2]-c1[2]*c2[1]]); if (not(EllPointOnCurve(ec,Gphi))) or (Gphi[2]=RingZero(ec.R) and Gphi<>EllZero(ec)) then Error("Cover should be of degree 2 and, if even, have Gphi=[0,1,0]."); fi; ec.c:=[c1,c2]; ec.Gphi:=Gphi; end; ############################################################################### EllPointRed:=function(ec,Plc,G) # calculates the reduction of a point at a place which should be of good red. local e,p,Gr; e:=Minimum(PlaceVal(Plc,G[1]),PlaceVal(Plc,G[2]),PlaceVal(Plc,G[3])); p:=PlacePrime(Plc); # A little trickery to still use EllPointNormal, while we don't have # the reduced ec at our disposal. (calling EllRed would be messy since # EllPointRed is also used there. This would create an ugly cross # dependence. Gr:=List(PlacePointNormal(Plc,G),u->EltToFFE(u,p)); return EllPointNormal(rec(R:=IdealResidueField(p)),Gr); end; ############################################################################### EllRed:=function(ec,Plc) # looks up the reduction of ec or computes it. The reduced curve inherits # Gs, such that EllPointRed and EllVecToPnt commute. # If ec.Phi is present, it is also inherited and should be of good reduction local Prm, r, redec, G; if not(IsBound(ec.red)) then ec.red:=[]; fi; Prm:=PlacePrime(Plc); for r in ec.red do if Prm=PlacePrime(r[1]) then return r[2]; fi; od; redec:=EllInit(EltToFFE(ec.gamma,Prm),0,EltToFFE(ec.a2,Prm),0, EltToFFE(ec.a4,Prm),EltToFFE(ec.a6,Prm)); Add(ec.red,[Plc,redec]); if IsBound(ec.Gs) then EllGenInit(redec,List(ec.Gs,G->EllPointRed(ec,Plc,G)),ec.Rk); fi; if IsBound(ec.c) then EllPhiInit(redec,List(ec.c[1],u->EltToFFE(u,Prm)), List(ec.c[2],u->EltToFFE(u,Prm))); if redec.Gphi=EllZero(redec) and ec.Gphi<>EllZero(ec) then Error("Cover has bad reduction"); fi; fi; return redec; end; ############################################################################### EllSize:=function(ec) # Count points on a curve over a finite field. Very inefficient implementation, # but since we only need curves over very small fields, this is not a problem. local N,c,S,w,r,i; if not(IsFF(ec.R)) then Error("ec not over a finite field"); fi; if IsBound(ec.size) then return ec.size; fi; # Get size of field. N:=FFSize(ec.R);N:=N[1]^N[2]; # count points with X=infty, X=0 c:=EllYsqr(ec,0); if c=RingZero(ec.R) then S:=2; elif FFEltLog(c) mod 2=0 then S:=3; else S:=1; fi; # use multiplicative structure to generate rest of X candidates. w:=FFPrimitiveElt(ec.R); r:=RingOne(ec.R); for i in [1..N-1] do c:=EllYsqr(ec,r); if c=RingZero(ec.R) then S:=S+1; elif FFEltLog(c) mod 2 =0 then S:=S+2; fi; r:=r*w; od; ec.size:=S; return S; end; ############################################################################### EllTau:=function(ec,G) #involution associated to the installed cover. return EllNeg(ec,EllAdd(ec,ec.Gphi,G)); end; ############################################################################### EllPhi:=function(ec,G) #calculate value of installed cover at a point. (takes care of G=Gphi) local res; if G=ec.Gphi then if ec.Gphi=EllZero(ec) then if ec.c[2][1]=RingZero(ec.R) then return [RingOne(ec.R),RingZero(ec.R)]; else return [ec.c[1][1]/ec.c[2][1],RingOne(ec.R)]; fi; else G:=EllTau(ec,G); fi; fi; res:=[ec.c[1][1]*G[1]+ec.c[1][2]*G[2]+ec.c[1][3]*G[3], ec.c[2][1]*G[1]+ec.c[2][2]*G[2]+ec.c[2][3]*G[3]]; if res[2]=RingZero(ec.R) then return [RingOne(ec.R),RingZero(ec.R)]; else return [res[1]/res[2],RingOne(ec.R)]; fi; end; ############################################################################### EllSpan:=function(ec,Gs) # calculates the "span" of Gs (Gs should generate torsion only): # returns a list L=[M,arr], where M is a list of lists, the rows of which # generate the Z-relations of Gs and arr is an array with as entries the # linear combination of Gs corresponding to the index. # M is in row HNF. local arr,G,i,v,M,R; if Length(Gs)=1 then arr:=ArrInit(1,[EllZero(ec)]); G:=Gs[1]; i:=1; while G<>EllZero(ec) do ArrAddRow(arr,G); G:=EllAdd(ec,G,Gs[1]); i:=i+1; od; return([[[i]],arr]); else R:=EllSpan(ec,Tail(Gs)); arr:=ArrInit(Length(Gs),[ArrToList(R[2])]); G:=ArrApply(R[2],u->EllAdd(ec,u,Gs[1])); v:=ArrLookup(G,EllZero(ec)); i:=1; while v=false do ArrAddRow(arr,G); G:=ArrApply(G,u->EllAdd(ec,u,Gs[1])); v:=ArrLookup(G,EllZero(ec)); i:=i+1; od; M:=[[i]]; Append(M[1],v); Append(M,List(R[1],u->ListAugFront(0,u))); return([M,arr]); fi; end; ############################################################################### EllGroup:=function(ec) # Returns an array refrecting the group of points on ec # The data is stored for future reference. (EllKernel is precomputed as a # side-effect) local tmp; if not(IsFF(ec.R)) then Error("Only applicable to elliptic curves over finite fields"); fi; if IsBound(ec.group) then return ec.group; fi; tmp:=EllSpan(ec,ec.Gs); ec.kerlat:=tmp[1]; ec.group:=tmp[2]; return ec.group; end; ############################################################################### EllKernel:=function(ec) # Returns a matrix the rows of which generate the Z-linear relations of Gs. # This describes the kernel of reduction if Gs comes from reduction on a # characteristic 0 curve. The data is stored for future ref. # (EllGroup is precomputed as a side-effect) local tmp; if not(IsFF(ec.R)) then Error("Only applicable to elliptic curves over finite fields"); fi; if IsBound(ec.kerlat) then return Mat(Z,ec.kerlat); fi; tmp:=EllSpan(ec,ec.Gs); ec.kerlat:=tmp[1]; ec.group:=tmp[2]; return Mat(Z,ec.kerlat); end; ############################################################################### EllFibers:=function(ec) # Returns the fibers of the rational (wrt. the primitive field) points # under Phi. Returned as a list of vectors wrt. Gs. Positions 1 .. p-1 contain # the fibers of 1 .. p-1. 0 is encoded as p and infinity as p+1. # Note that these vectors should be interpreted "mod EllKernel". local tmp,p,i; if not(IsFF(ec.R)) then Error("Only applicable to elliptic curves over finite fields"); fi; if IsBound(ec.fibers) then return ec.fibers; fi; tmp:=ArrApply(EllGroup(ec),G->EllPhi(ec,G)); p:=FFSize(ec.R)[1]; ec.fibers:=[]; for i in [1..p-1] do ec.fibers[i]:=ArrLookupAll(tmp,[RingForce(ec.R,i),RingOne(ec.R)]); od; ec.fibers[p]:=ArrLookupAll(tmp,[RingZero(ec.R),RingOne(ec.R)]); ec.fibers[p+1]:=ArrLookupAll(tmp,[RingOne(ec.R),RingZero(ec.R)]); return ec.fibers; end; ############################################################################### EllRatFibers:=function(ec,PlcList) #intersects EllFiber for the places in PlcList (which should be over one # rational prime p) to determine a candidate for # (Phi(E(L)) cap P_1(Q)) mod p. # and returns this. local p,L,ec1,fibers1,M1,ec2,fibers2,M2,M3,R,i,Plc1,Plc2; p:=PlcList[1].primebelow; if ForAny(PlcList,u->p<>u.primebelow) then Error("Places should all be above same rational prime."); fi; L:=Set([1..p+1]); for Plc1 in PlcList do ec1:=EllRed(ec,Plc1); fibers1:=EllFibers(ec1); M1:=EllKernel(ec1); for Plc2 in PlcList do ec2:=EllRed(ec,Plc2); fibers2:=EllFibers(ec2); M2:=EllKernel(ec2); M3:=LatSum(M1,M2); #looping through L would suffice, but L is changed in the loop, which #makes it dangerous to use L. #We cast out i from L if none of the elements in the fibers agree mod #the sum of the kernels of reduction. for i in [1..p+1] do if not(ForAny(fibers1[i], u->ForAny(fibers2[i], v->LatContainsVec(M3,u-v)))) then SetRemove(L,i); fi; od; od; od; # decode results to return an understandable answer. R:=[]; if p+1 in L then Add(R,[1,0]); fi; if p in L then Add(R,[0,1]); fi; for i in [1..p-1] do if i in L then Add(R,[i,1]); fi; od; return R; end; ############################################################################### EllAddPadicSpec:=function(ec,Plc,prec,G1,G2) # returns the sum of two points with finite p-adic precision. # G1 and G2 should be reduced (integral and primitive at Plc) # G1 and G2 should be distant and not near infty. local G; G:=[(G1[1]*G2[3]-G2[1]*G1[3])*(G1[1]^2*G2[1]*G2[3]+ec.a4*G1[1]*G1[3]*G2[3]^2+ 2*ec.a2*G1[1]*G2[1]*G1[3]*G2[3]+G1[1]*G2[1]^2*G1[3]-2*ec.gamma*G1[2]*G2[2] *G1[3]*G2[3]+2*ec.a6*G1[3]^2*G2[3]^2+ec.a4*G2[1]*G1[3]^2*G2[3]), -2*G1[2]*ec.a2* G1[1]*G2[1]*G1[3]*G2[3]^2-G1[2]*G2[1]^3*G1[3]^2-4*G1[2]*ec.a6*G1[3]^2*G2[3]^3-3 *G1[2]*G1[1]*G2[1]^2*G1[3]*G2[3]-G1[2]*ec.a4*G1[1]*G1[3]*G2[3]^3-3*G1[2]*ec.a4* G2[1]*G1[3]^2*G2[3]^2-2*G1[2]*ec.a2*G2[1]^2*G1[3]^2*G2[3]+G2[2]*G1[1]^3*G2[3]^2 +3*G2[2]*G1[1]^2*G2[1]*G1[3]*G2[3]+2*G2[2]*ec.a2*G1[1]*G2[1]*G1[3]^2*G2[3]+ 4*G2[2]*ec.a6*G1[3]^3*G2[3]^2+3*G2[2]*ec.a4*G1[1]*G1[3]^2*G2[3]^2+ G2[2]*ec.a4*G2[1]*G1[3]^3*G2[3]+2*G2[2]*ec.a2*G1[1]^2*G1[3]*G2[3]^2, (G1[1]*G2[3]-G2[1]*G1[3])^3]; if Minimum(List(G,u->PlaceVal(Plc,u)))<>0 then Error("Cannot do this addition (G1, G2 and infty should be distant)"); fi; return List(G,u->u mod PlaceGetIdeal(Plc,prec)); end; ############################################################################### EllAddPadic:=function(ec,Plc,prec,G1,G2) # Do an addition with limited p-adic precision. If G1 and G2 reduce to the same # point, then a third point, Gt, is chosen and we use G1+G2=((G1+Gt)+G2)-Gt. # Gt is chosen from Gs. If a suitable point cannot be found, an error is # generated. local ecr,G1r,G2r,i,Gt,R; ecr:=EllRed(ec,Plc); G1r:=EllPointRed(ec,Plc,G1); G2r:=EllPointRed(ec,Plc,G2); G1:=PlacePointNormal(Plc,G1); G2:=PlacePointNormal(Plc,G2); if G1r=EllZero(ecr) or G2r=EllZero(ecr) then Error("Cannot add points in kernel of reduction"); elif G1r<>G2r then return EllAddPadicSpec(ec,Plc,prec,G1,G2); fi; i:=1; Gt:=EllMul(ecr,-2,G1r); while i<=Length(ecr.Gs) and ( ecr.Gs[i]=EllZero(ecr) or ecr.Gs[i]=G1r or ecr.Gs[i]=Gt or EllMul(ecr,2,ecr.Gs[i])=Gt) do i:=i+1; od; if i>Length(ecr.Gs) then Error("Could not find a suitable translation to do addition."); else Gt:=PlacePointNormal(Plc,ec.Gs[i]); R:=EllAddPadicSpec(ec,Plc,prec,G1,Gt); R:=EllAddPadicSpec(ec,Plc,prec,R,G2); R:=EllAddPadicSpec(ec,Plc,prec,R,[Gt[1],-Gt[2],Gt[3]]); return R; fi; end; ############################################################################### EllMulPadic:=function(ec,Plc,prec,n,G) # determines n*G with finite P-adic precision. using repeated doubling. # Same restrictions as for EllAddPadic. local Gt; if n=0 then return EllZero(ec); elif n=1 then return PlacePointNormal(Plc,G); elif n<0 then Gt:=EllMulPadic(ec,Plc,prec,-n,G); return([Gt[1],-Gt[2] mod PlaceGetIdeal(Plc,prec),Gt[3]]); elif n mod 2 = 0 then Gt:=EllMulPadic(ec,Plc,prec,n/2,G); return EllAddPadic(ec,Plc,prec,Gt,Gt); else Gt:=EllMulPadic(ec,Plc,prec,(n-1)/2,G); return EllAddPadic(ec,Plc,prec,EllAddPadic(ec,Plc,prec,Gt,Gt),G); fi; end; ############################################################################### EllVecToPadicPnt:=function(ec,Plc,prec,v) # calculates the specified linear combination of the preregistered Gs. local G,i,Gn,e; if Length(v)<>Length(ec.Gs) then Error("Vector of wrong dimension"); fi; G:=EllZero(ec); for i in [1..Length(v)] do Gn:=ec.Gs[i]; e:=Minimum(List(Gn,u->PlaceVal(Plc,u))); Gn:=List(Gn,u->u/Plc.unif^e mod PlaceGetIdeal(Plc,prec)); if v[i]<>0 then if G=EllZero(ec) then G:=EllMulPadic(ec,Plc,prec,v[i],Gn); else G:=EllAddPadic(ec,Plc,prec,G,EllMulPadic(ec,Plc,prec,v[i],Gn)); fi; fi; od; return G; end; ############################################################################### EllKerZ:=function(ec,Plc) # computes the Z-coordinates of the kernel of reduction of Gs at Plc. # This data is stored in the reduced curve, which is strictly incorrect, # since it has nothing to do with the finite field and cannot be computed # from EllRed alone. However, no-one should notice this data here and it # saves a whole new data structure for "localisations" of the curve. local ecr,Gk; ecr:=EllRed(ec,Plc); if IsBound(ecr.kerZ) then return ecr.kerZ; fi; Gk:=List(MatToRowList(EllKernel(ecr)), u->EllVecToPadicPnt(ec,Plc,2,u)); ecr.kerZ:=List(Gk,u->u[1]/u[2] mod PlaceGetIdeal(Plc,2)); return ecr.kerZ; end; ############################################################################### EllZB:=function(ec,PlcList) local p,L,Plc,Mr,row; if IsBound(ec.PlcList) and PlcList=ec.PlcList then return ec.ZB; fi; p:=PlcList[1].primebelow; if ForAny(PlcList,Plc->Plc.primebelow<>p) or ForAny(PlcList,Plc->Plc.ramdeg>1) then Error("Places should all be of same characteristic and unramified."); fi; ec.PlcList:=PlcList; L:=MatId(Z,Length(ec.Gs)); for Plc in PlcList do L:=LatSect(L,EllKernel(EllRed(ec,Plc))); od; ec.M:=L; L:=[]; for Plc in PlcList do Mr:=LatRelBasis(ec.M,EllKernel(EllRed(ec,Plc))); row:=MatToColList(MatMove(Mr,ec.R)* MatTrans(Mat(ec.R,[EllKerZ(ec,Plc)])))[1]; Add(L,List(row{[1..ec.Rk]},u->(u/p) mod PlacePrime(Plc))); od; ec.PlcList:=PlcList; ec.ZB:=Mat(ec.R,L); return ec.ZB; end; ############################################################################### MatDecompose:=function(M) #Given a matrix over an order O, returns Degree(O) matrices over Q that form #the decomposition of M over the Z-basis of O. local O,L; if not(IsElt(MatElt(M,1,1))) then Error("Matrix should be over an order"); fi; O:=EltOrder(MatElt(M,1,1)); L:=List(MatToRowList(M),row->List(row,EltToList)); return List([1..OrderDeg(O)],i->Mat(Z,List(L,row->List(row,a->a[i])))); end; ############################################################################### EllLocalIndexOne:=function(ec,PlcList) # Assuming that the Gs span a subgroup of finite index in the Mordell-Weil # group, test if the kernel is generated by the appropriate linear # combinations of the Gs. (EllZB should not have nontrivial F_p relations, # where F_p is the primitive field associated to the places in PlcList) local M,t; M:=[]; for t in MatDecompose(EllZB(ec,PlcList)) do Append(M,MatToRowList(t)); od; return MatRows(MatKernel(MatTrans(Mat(FF(PlcList[1].primebelow),M))))=0; end; ############################################################################### EllTheta:=function(ec,PlcList,G) # Computes the "magic" Theta power series that bounds the number of # points with rational image under phi near G. local phival,c,Gf1,Gf2,Gf3,x,y,fpx,ZB,quad,C,L,ref,rat,aux,l,k,i,j,c; G:=EllPointNormal(ec,G); phival:=EllPhi(ec,G); c:=ec.c; if G=ec.Gphi and G<>EllZero(ec) then Error("Cannot do theta around G=Gphi if G<>infinity"); fi; if phival[2]=0 then Print("Point maps to infinity. Taking 1/phi\n"); c:=[c[2],c[1]]; phival:=0; else L:=EltToList(phival[1]); if ForAny(L{[2..Length(L)]},u->u<>0) then Error("Phi(G) not rational"); fi; phival:=L[1]; fi; # First we compute the power series. This will either be a # square approximation or a linear one (denoted by the boolean "quad") Gf1:=c[1][2]*c[2][3]-c[1][3]*c[2][2]; Gf2:=c[1][3]*c[2][1]-c[1][1]*c[2][3]; Gf3:=c[1][1]*c[2][2]-c[1][2]*c[2][1]; ZB:=MatToRowList(EllZB(ec,PlcList)); if ec.Gphi=EllZero(ec) and G=EllZero(ec) then quad:=true; C:=Gf2/ec.gamma/c[2][1]^2; elif ec.Gphi=EllZero(ec) and G[2]=0 then quad:=true; x:=G[1]/G[3]; fpx:=3*x^2+2*ec.a2*x+ec.a4; C:=-fpx*Gf2/ec.gamma/(c[2][1]*x+c[2][3])^2; elif ec.Gphi<>EllZero(ec) and G=EllZero(ec) then quad:=false; C:=Gf3/c[2][2]^2; else quad:=false; x:=G[1]/G[3]; y:=G[2]/G[3]; fpx:=3*x^2+2*ec.a2*x+ec.a4; C:=(fpx*(x*Gf3-Gf1)-2*ec.gamma*y*(y*Gf3-Gf2))/ (ec.gamma*(c[2][1]*x+c[2][2]*y+c[2][3])^2); fi; # L is either going to be a list of matrices (describing quadr. forms) # or a list of row matrices (describing the linear approximation. The only # thing we have to do is multiply C with ZB in the proper way. if quad then L:=List([1..Length(PlcList)],k-> Mat(ec.R,List([1..ec.Rk],i->List([1..ec.Rk],j-> (C*ZB[k][i]*ZB[k][j]) mod PlacePrime(PlcList[k]))))); else L:=MatToRowList(C*ZB); L:=List([1..Length(PlcList)],k-> Mat(ec.R,[List(L[k],u->u mod PlacePrime(PlcList[k]))])); fi; # we decompose the matrices in L wrt. a Z-basis. C:=List(L,MatDecompose); #a zero row of right dimension for comparison. c:=Mat(Z,List(MatToRowList(L[1]),row->List(row,u->0))); ref:=[]; rat:=[];aux:=[]; for l in C do # if we do not have a row describing the rational value, take it if ref=[] then ref:=l[1]; else # otherwise we have an equation. Add(rat,l[1]-ref); fi; # the components on non-rational basis vectors (that are non-zero) give # us equations as well. Append(aux,Filtered(Tail(l),u->u<>c)); od; # We put the equations coming from inertness at the end. They depend on # choice of basis much more than other equations. Append(rat,aux); # we return either a list of matrices or one matrix over the primitive finite # field. if quad then return List(rat,mat->Mat(FF(PlcList[1].primebelow),MatToRowList(mat))); else return Mat(FF(PlcList[1].primebelow),List(rat,mat->MatToRowList(mat)[1])); fi; end; ############################################################################### QuadHasSols:=function(M) # given a list of matrices over a finite field, describing quadratic forms, # determines if these simultaneously represent 0 in a non-trivial way. local F,p,idx1,i,mc,v,n; F:=FFEltFF(MatElt(M[1],1,1)); p:=FFSize(F)[1]; if FFSize(F)[2]>1 then Error("Non-primitive fields not implemented"); fi; n:=MatRows(M[1]); if ForAll(M,m->MatElt(m,1,1)=RingZero(F)) then return true; fi; for idx1 in [2..n] do mc:=List([1..n],i->0); mc[idx1]:=1; for i in [1..p^(idx1-1)] do v:=Mat(F,[mc]); if ForAll(M,m->MatElt(v*m*MatTrans(m),1,1)=RingZero(F)) then return true; fi; MultiCounterInc(mc,p); od; od; return false; end; ############################################################################### EllThetaTest:=function(ec,PlcList,G) # Tests whether EllTheta around G proves that G is only point with rational # phi in its fiber of reduction. local L,tht,m; Print("Computing Theta^G for G=",G,"...\n"); L:=List(PlcList,Plc->List(EllPhi(EllRed(ec,Plc),EllPointRed(ec,Plc,G)), u->FFEToElt(u,PlacePrime(Plc)))); if ForAny(L,u->u<>L[1]) then Print("No points near G have rational image\n"); return true; fi; tht:=EllTheta(ec,PlcList,G); if TYPE(tht)= "KANT matrix" then Print("G is only point in fiber if the following matrix has maximal rank mod ", PlcList[1].primebelow,"\n",tht,"\n\n"); return MatRows(MatKernel(MatTrans(tht)))=0; else Print( "G is only point in fiber if the intersection of the following quadratic\n", "forms has no rat. points mod ",PlcList[1].primebelow,"\n"); for m in tht do Print(m,"\n\n"); od; return not(QuadHasSols(tht)); fi; end; ############################################################################### EllChab:=function(ec,PlcList,fiblist) # local imglist,imglistred,t; Print("Testing whether fiblist covers rational image mod ", PlcList[1].primebelow,"...\n"); imglist:=List(fiblist,G->EllPhi(ec,G)); imglistred:=[]; for t in imglist do if t[2]=RingZero(ec.R) or PlaceVal(PlcList[1],t[1])<0 then Add(imglistred,[1,0]); else Add(imglistred, [EltMove(t[1] mod PlacePrime(PlcList[1]),Z)mod PlcList[1].primebelow,1]); fi; od; Print("phi(fiblist) mod ",PlcList[1].primebelow,"=",Set(imglistred),"\n"); Print("Local considerations give\n", "phi(E(L)) cap PP_1(Q) mod ",PlcList[1].primebelow," subset of ", Set(EllRatFibers(ec,PlcList)),"\n"); if not(SetIsEqual(Set(imglistred),Set(EllRatFibers(ec,PlcList)))) then Error("fiblist does not cover potential non-empty fibers"); else Print("Tested OK.\n"); fi; Print("Testing if points are unique in reduction fiber ...\n"); if ForAll(fiblist,G->EllThetaTest(ec,PlcList,G)) then Print("All fibers tested OK. Phi(E(K)) cap PP_1(Q)=\n",imglist,"\n"); return true; else Print("Some fiber not decisive. Try another prime.\n"); return false; fi; end; #--- ellcalc.g (end) ---# #--- eq238E5.g (start) ---# cat << "#--- eq238E5.g (end) ---#" > eq238E5.g # UNIX extraction command Read("isogdesc.g"); O:=OrderMaximal(x^3-2); beta:=Elt(O,[0,1,0]); u3:=1+beta; u2:=beta; p2:=u2*O; Plc2:=PlaceInit(p2); p5:=List([beta+2,beta^2-2*beta-1],u->Ideal(5,u)); Plc5:=List(p5,u->PlaceInit(u)); ec:=EllInit(Elt(O,-6),0,0,0,-1,0); Gs:=[[-2,1,1],[1,0,1],[0,0,1]]; Gps:=[[-2-2*beta,-(4+2*beta+4*beta^2)/3,1],[0,0,1]]; EllGenInit(ec,[Gs[1],EllPsip(ec,Gps[1]),Gs[2],Gs[3]],2); EllPhiInit(ec,[beta,-6*beta,beta],[2,0,-1]); fiblist:=List([[0,0,0,0],[-1,0,1,0],[0,1,0,0]],v->EllVecToPnt(ec,v)); PlcList:=Plc5; TorPlc:=[Plc5[1],PlaceInit(Ideal(11,beta-7))]; GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #--- eq238E5.g (end) ---# #--- eq238E7.g (start) ---# cat << "#--- eq238E7.g (end) ---#" > eq238E7.g # UNIX extraction command Read("isogdesc.g"); Oeq:=Order(x^12 + 6*x^10 + 39*x^8 + 64*x^6 + 15*x^4 - 6*x^2 - 3); OrderAutomorphisms(Oeq,[Elt(Oeq,[0,-1,0,0,0,0,0,0,0,0,0,0])]); CFmat:=MatTrans(Mat(Z, [[1832,0,0,0,0,0,0,0,0,0,0,0], [3050,0,2296,0,-19120,0,-13784,0,-2110,0,-372,0], [-52,0,2750,0,17748,0,12504,0,1892,0,334,0], [0,1832,0,0,0,0,0,0,0,0,0,0], [-336,0,-6998,0,3632,0,4908,0,740,0,150,0], [-1636,0,-6032,0,8568,0,8816,0,1324,0,256,0], [-2175,1855,-665,-7003,12028,-18296,8846,-9406,1347,-1471,239,-227], [194,1316,2124,-1636,-10690,-14836,-8706,-9376,-1316,-1448,-242,-244], [-635,-4149,1143,-5539,18434,4970,13144,5048,2001,771,353,147], [-635,4149,1143,5539,18434,-4970,13144,-5048,2001,-771,353,-147], [-194,464,-2124,4386,10690,32584,8706,21880,1316,3340,242,578], [-2317,-2511,-5539,-7663,4970,15660,5048,13754,771,2087,147,389]])); O:=Order(Oeq,CFmat,1832); if OrderDisc(OrderMaximal(O))=OrderDisc(O) then OrderSetMaximal(O); fi; sigma:=u->EltAutomorphism(u,2); alpha:=EltMove(Elt(Oeq,[0,1,0,0,0,0,0,0,0,0,0,0]),O); up2:=Elt(O,[0,1,0,0,1,0,1,0,0,0,0,0]); p2:=up2*O; up3:=Elt(O,[1,0,1,0,0,0,0,0,0,0,-1,0]); p3:=up3*O; uq3:=Elt(O,[1,0,0,0,0,-1,0,0,1,1,0,0]); q3:=uq3*O; Plc2:=PlaceInit(p2); p31:=List([alpha-5,alpha+12,alpha+5,alpha-12,alpha^2-5*alpha-13],u->Ideal(31,u)); Plc31:=List(p31,u->PlaceInit(u)); gamma7:=Elt(O,[0,0,0,-1,-1,0,-1,0,0,0,1,0]); ec:=EllInit(gamma7,0,2,0,2,0); Gs:=List([[[0,-1,6,-4,0,-1,-2,0,-4,0,-1,4],[8,-2,23,-12,-2,0,-6,7,-18,0,-6,14]], [[-1,-1,0,0,0,1,-1,1,0,0,0,-1],[0,0,0,1,1,-1,1,0,0,0,-1,-1]], [[1,0,0,-1,0,0,1,-1,-1,1,0,-1],[0,0,0,2,0,2,0,0,0,-4,0,-2]], [0,0]],u->[Elt(O,u[1]),Elt(O,u[2]),1]); Gps:=List([[[1,-1,-1,0,1,-2,0,1,1,1,-1,0],[0,-2,0,0,0,0,-4,0,0,2,-2,2]], [0,0]],u->[Elt(O,u[1]),Elt(O,u[2]),1]); EllGenInit(ec,[Gs[1],Gs[2],Gs[3],EllPsip(ec,Gps[1]),Gs[4]],4); EllPhiInit(ec, [Elt(O,[4,-4,-4,-2,0,6,-2,8,4,-6,10,-8]),-12, Elt(O,[2,1,2,5,2,-2,1,0,1,-3,-1,-3])], [Elt(O,[0,4,0,4,0,0,0,-4,4,0,0,4]),0, Elt(O,[4,-3,0,5,2,-8,-1,-4,3,5,-5,5])]); fiblist:=List([[0,0,0,0,0],[0,0,0,1,0]],v->EllVecToPnt(ec,v)); PlcList:=Plc31; TorPlc:=[Plc31[1],PlaceInit(Ideal(43,alpha-6))]; GlobalCertAddHint(ec,Plc2, List([[0,0,-7,-1,8,8,-4,-2,5,1,4,2],[1,1,-5,0,3,6,2,-5,5,3,-3,2], [-7,2,-4,5,-2,2,13,-7,-5,-3,2,-5],[11,12,4,3,-6,-6,9,-5,1,-7,6,5], [3,-4,6,5,0,6,1,7,-7,-1,4,-5],[-15,8,-6,-5,-6,2,-15,7,5,1,-4,-3]], u->Elt(O,u)), List([[18,-58,-48,-14,20,-28,-86,-20,-56,56,-2,54], [-78,-74,-8,-54,-28,-8,90,56,-28,48,6,26]],u->Elt(O,u))); GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #--- eq238E7.g (end) ---# #--- eq238E9.g (start) ---# cat << "#--- eq238E9.g (end) ---#" > eq238E9.g # UNIX extraction command Read("isogdesc.g"); Oeq:=Order(x^12 + 6*x^10 + 39*x^8 + 64*x^6 + 15*x^4 - 6*x^2 - 3); OrderAutomorphisms(Oeq,[Elt(Oeq,[0,-1,0,0,0,0,0,0,0,0,0,0])]); CFmat:=MatTrans(Mat(Z, [[1832,0,0,0,0,0,0,0,0,0,0,0], [3050,0,2296,0,-19120,0,-13784,0,-2110,0,-372,0], [-52,0,2750,0,17748,0,12504,0,1892,0,334,0], [0,1832,0,0,0,0,0,0,0,0,0,0], [-336,0,-6998,0,3632,0,4908,0,740,0,150,0], [-1636,0,-6032,0,8568,0,8816,0,1324,0,256,0], [-2175,1855,-665,-7003,12028,-18296,8846,-9406,1347,-1471,239,-227], [194,1316,2124,-1636,-10690,-14836,-8706,-9376,-1316,-1448,-242,-244], [-635,-4149,1143,-5539,18434,4970,13144,5048,2001,771,353,147], [-635,4149,1143,5539,18434,-4970,13144,-5048,2001,-771,353,-147], [-194,464,-2124,4386,10690,32584,8706,21880,1316,3340,242,578], [-2317,-2511,-5539,-7663,4970,15660,5048,13754,771,2087,147,389]])); O:=Order(Oeq,CFmat,1832); if OrderDisc(OrderMaximal(O))=OrderDisc(O) then OrderSetMaximal(O); fi; sigma:=u->EltAutomorphism(u,2); alpha:=EltMove(Elt(Oeq,[0,1,0,0,0,0,0,0,0,0,0,0]),O); up2:=Elt(O,[0,1,0,0,1,0,1,0,0,0,0,0]); p2:=up2*O; up3:=Elt(O,[1,0,1,0,0,0,0,0,0,0,-1,0]); p3:=up3*O; uq3:=Elt(O,[1,0,0,0,0,-1,0,0,1,1,0,0]); q3:=uq3*O; Plc2:=PlaceInit(p2); p31:=List([alpha-5,alpha+12,alpha+5,alpha-12,alpha^2-5*alpha-13],u->Ideal(31,u)); Plc31:=List(p31,u->PlaceInit(u)); gamma9:=Elt(O,[0,0,0,1,1,-1,1,1,0,-1,1,-1]); ec:=EllInit(gamma9,0,2,0,2,0); Gs:=List([[[1/3,-2/3,-4/3,-1/3,0,2/3,-1/3,1/3,1/3,1/3,2/3,-1],[-4/3,-2/3,0,0,0,2/3,4/3,-2/3,-4/3,2/3,0,4/3]], [[-7/3,1/3,-4,2/3,-2/3,2/3,2/3,-5/3,2/3,2/3,5/3,-2],[28/3,-2,16/3,12,-2/3,-2,2/3,26/3,10,-4,-4,2/3]], [[-1,2,0,0,-1,-2,1,-1,-1,1,0,1],[0,-4,-2,-1,-5,1,-5,2,4,-2,5,-1]], [0,0]],u->[Elt(O,u[1]),Elt(O,u[2]),1]); Gps:=List([[[2,-2,0,0,0,2,-2,2,0,0,0,-2],[-2,0,-4,2,-4,0,-2,2,4,-4,4,-2]], [0,0]],u->[Elt(O,u[1]),Elt(O,u[2]),1]); EllGenInit(ec,[Gs[1],Gs[2],Gs[3], EllAdd(ec,EllMul(ec,2,EllAdd(ec,Gs[1],Gs[2])), EllNeg(ec,EllPsip(ec,Gps[1]))), Gs[4]],4); EllPhiInit(ec, [Elt(O,[0,2,6,0,-6,10,-10,6,-2,-10,2,2]),-12, Elt(O,[0,7,2,1,-2,2,-3,2,1,-7,3,-1])], [Elt(O,[-4,0,0,0,-4,0,0,0,0,-4,4,4]),0, Elt(O,[0,3,2,-1,0,2,-1,0,-1,-3,1,-1])]); fiblist:=List([[0,0,0,0,0],[0,0,0,-1,0]],v->EllVecToPnt(ec,v)); PlcList:=Plc31; TorPlc:=[Plc31[1],PlaceInit(Ideal(43,alpha-6))]; GlobalCertAddHint(ec,Plc2, List([[-6,6,6,-1,-3,-7,5,-6,4,4,5,-1],[3,-1,5,0,-5,-4,-6,7,7,1,-3,-2], [-1,-13,4,6,0,-4,-6,3,6,-4,5,8],[3,-7,6,-4,0,6,6,-7,-6,8,-1,4], [1,1,-2,6,-6,8,-6,-3,8,6,3,-4],[1,-15,2,0,-6,4,16,7,-4,0,7,0]],u->Elt(O,u)), List([[50,46,-32,22,36,20,110,-32,-60,-4,-38,-46], [-114,-30,40,10,48,-12,2,16,16,52,-6,-22]],u->Elt(O,u))); GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #--- eq238E9.g (end) ---# #--- eq245E3.g (start) ---# cat << "#--- eq245E3.g (end) ---#" > eq245E3.g # UNIX extraction command Read("isogdesc.g"); O:=OrderMaximal(x^4-5*x^2+5); beta:=Elt(O,[0,1,0,0]); u5:=beta; u2:=1+beta-beta^2; p2:=u2*O; Plc2:=PlaceInit(p2); p41:=List([beta-3,beta-18,beta-23,beta-38],u->Ideal(41,u)); Plc41:=List(p41,u->PlaceInit(u)); ec:=EllInit(beta^2+beta-1,0,-5,0,5,0); Gs:=[[10+5*beta-2*beta^2-beta^3,20+15*beta-5*beta^2-4*beta^3,1],[beta^2,0,1]]; Gps:=[[5-4*beta^2,0,1],[0,0,1]]; EllGenInit(ec,[Gs[1],Gs[2],EllPsip(ec,Gps[1])],1); EllPhiInit(ec,[6*beta-2*beta^3,0,3*beta^3-10*beta],[0,0,5]); fiblist:=List([[0,0,0]],v->EllVecToPnt(ec,v)); PlcList:=Plc41; TorPlc:=[Plc41[1],Plc41[3]]; GlobCert:=GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #--- eq245E3.g (end) ---# #--- eq245E6.g (start) ---# cat << "#--- eq245E6.g (end) ---#" > eq245E6.g # UNIX extraction command Read("isogdesc.g"); O:=OrderMaximal(x^4-x^3+x^2-x+1); zeta:=Elt(O,[0,1,0,0]); u5:=1+zeta+zeta^2+zeta^3; u2:=2; p2:=u2*O; Plc2:=PlaceInit(p2); p31:=List([zeta-15,zeta-23,zeta-27,zeta-29],u->Ideal(31,u)); Plc31:=List(p31,u->PlaceInit(u)); ec:=EllInit(2*(zeta^2-zeta),0,-5,0,5,0); Gs:=[[1-3*zeta^2+3*zeta^3,1+3*zeta-2*zeta^2+zeta^3,1], [1-zeta^2+zeta^3,zeta^3,1],[2-zeta^2+zeta^3,0,1]]; Gps:=[[-3+4*zeta^2-4*zeta^3,0,1]]; EllGenInit(ec,[Gs[1],Gs[2],Gs[3],EllPsip(ec,Gps[1])],2); EllPhiInit(ec,[Elt(O,1),0,-2*zeta^2+zeta-2],[Elt(O,1),0,zeta^3+zeta^2-zeta-1]); fiblist:=List([[0,0,0,0],[1,0,1,1]],v->EllVecToPnt(ec,v)); PlcList:=Plc31; TorPlc:=[Plc31[1],PlaceInit(Ideal(41,zeta-4))]; GlobalCertAddHint(ec,Plc2,[4-2*zeta,4-2*zeta^2],[]); GlobCert:=GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #brk> GlobCert.lcts[2].gens; #[ [1, 0, -3, 3], [1, 0, -1, 1], [2, 0, -1, 1], [44, -82, -24, -112], # [60, -8, 14, -48] ] #brk> GlobCert.lcts[2].gensp; #[ [-7431, 2472, -2916, 8508] ] #--- eq245E6.g (end) ---# #--- eq245E7.g (start) ---# cat << "#--- eq245E7.g (end) ---#" > eq245E7.g # UNIX extraction command Read("isogdesc.g"); O:=OrderMaximal(x^4-x^3+x^2-x+1); zeta:=Elt(O,[0,1,0,0]); u5:=1+zeta+zeta^2+zeta^3; u2:=2; p2:=u2*O; Plc2:=PlaceInit(p2); p31:=List([zeta-15,zeta-23,zeta-27,zeta-29],u->Ideal(31,u)); Plc31:=List(p31,u->PlaceInit(u)); ec:=EllInit(-2*(zeta^2-zeta),0,-5,0,5,0); Gs:=[[2+zeta^2-zeta^3,1-zeta-zeta^3,1], [-1+3*zeta^2-3*zeta^3,1-7*zeta+8*zeta^2-4*zeta^3,1],[2-zeta^2+zeta^3,0,1]]; Gps:=[[-3+4*zeta^2-4*zeta^3,0,1]]; EllGenInit(ec,[Gs[1],Gs[2],Gs[3],EllPsip(ec,Gps[1])],2); EllPhiInit(ec,[Elt(O,1),0,-2*zeta^2+zeta-2],[Elt(O,1),0,zeta^3+zeta^2-zeta-1]); fiblist:=List([[0,0,0,0],[1,1,0,1],[1,-1,0,1]],v->EllVecToPnt(ec,v)); PlcList:=Plc31; TorPlc:=[Plc31[1],PlaceInit(Ideal(41,zeta-4))]; GlobalCertAddHint(ec,Plc2,[4+10*zeta+8*zeta^2+8*zeta^3,-4+2*zeta^2+8*zeta^3],[]); GlobCert:=GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #--- eq245E7.g (end) ---# #--- eq245E8.g (start) ---# cat << "#--- eq245E8.g (end) ---#" > eq245E8.g # UNIX extraction command Read("isogdesc.g"); O:=OrderMaximal(x^5-2); alpha:=Elt(O,[0,1,0,0,0]); u5:=alpha+1; u2:=alpha; p2:=u2*O; Plc2:=PlaceInit(p2); p151:=List([alpha-22,alpha-25,alpha-49,alpha-90,alpha-116],u->Ideal(151,u)); Plc151:=List(p151,u->PlaceInit(u)); ec:=EllInit(Elt(O,1),0,-5,0,5,0); Gs:=[[1-2*alpha+alpha^3,3-2*alpha^2-2*alpha^3+alpha^4,1],[0,0,1]]; Gps:=[[5,-20,1]]; EllGenInit(ec,[Gs[1],EllPsip(ec,Gps[1]),Gs[2]],2); EllPhiInit(ec,[-alpha^3,2*alpha^3,0],[Elt(O,4),0,Elt(O,-5)]); fiblist:=List([[0,0,0],[0,1,0],[-2,0,0]],v->EllVecToPnt(ec,v)); PlcList:=Plc151; TorPlc:=[Plc151[3],PlaceInit(Ideal(3,alpha-2))]; GlobalCertAddHint(ec,Plc2,[],[]); GlobCert:=GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #--- eq245E8.g (end) ---# #--- eq245E9.g (start) ---# cat << "#--- eq245E9.g (end) ---#" > eq245E9.g # UNIX extraction command Read("isogdesc.g"); O:=OrderMaximal(x^5-2); alpha:=Elt(O,[0,1,0,0,0]); u5:=alpha+1; u2:=alpha; p2:=u2*O; Plc2:=PlaceInit(p2); p151:=List([alpha-22,alpha-25,alpha-49,alpha-90,alpha-116],u->Ideal(151,u)); Plc151:=List(p151,u->PlaceInit(u)); ec:=EllInit(alpha^3+alpha^2-1,0,-5,0,5,0); Gs:=[[1+alpha+alpha^2-alpha^3+alpha^4,1-5*alpha+6*alpha^2-4*alpha^3+2*alpha^4,1], [0,0,1]]; Gps:=[[1+8*alpha-4*alpha^3,-8*alpha^3+8*alpha^4-4*alpha+20-4*alpha^2,1]]; EllGenInit(ec,[Gs[1], EllAdd(ec,EllPsip(ec,Gps[1]),EllNeg(ec,Gs[1])),Gs[2]],2); EllPhiInit(ec,[2*(6*alpha^4-4*alpha^3+5*alpha^2+11*alpha-13),12, -10*(3*alpha^4-alpha^3+4*alpha-6)], [-4*alpha^3+2*alpha^2+2*alpha+2,-12,-5*(3*alpha^4-3*alpha^3+alpha^2-4)]); fiblist:=List([[0,0,0]],v->EllVecToPnt(ec,v)); PlcList:=Plc151; TorPlc:=[PlaceInit(Ideal(3,alpha-2))]; GlobalCertAddHint(ec,Plc2,[],[]); GlobCert:=GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; #phi has bad reduction at TorPlc so we remove phi temporarily to avoid an error #note that phi is not important for bounding the torsion anyway. c:=ec.c; Unbind(ec.c); Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); ec.c:=c; EllChab(ec,PlcList,fiblist); #--- eq245E9.g (end) ---# #--- eq283E2.g (start) ---# cat << "#--- eq283E2.g (end) ---#" > eq283E2.g # UNIX extraction command Read("isogdesc.g"); O:=OrderMaximal(x^2-3); alpha:=Elt(O,[0,1]); u3:=alpha; u2:=1+alpha; p2:=u2*O; Plc2:=PlaceInit(p2); p23:=[(-3*alpha-2)*O,(3*alpha-2)*O]; Plc23:=List(p23,u->PlaceInit(u)); ec:=EllInit(2,0,0,0,-2*alpha,0); Gs:=[[0,0,1]]; Gps:=[[2+2*alpha,4+4*alpha,1],[0,0,1]]; EllGenInit(ec,[EllPsip(ec,Gps[1]),Gs[1]],1); EllPhiInit(ec,[1,0,0],[0,0,1]); fiblist:=List([[0,0],[0,1],[1,0]],v->EllVecToPnt(ec,v)); PlcList:=Plc23; TorPlc:=[Plc23[1],PlaceInit(Ideal(13,alpha-9))]; GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #--- eq283E2.g (end) ---# #--- eq283E4.g (start) ---# cat << "#--- eq283E4.g (end) ---#" > eq283E4.g # UNIX extraction command Read("isogdesc.g"); O:=OrderMaximal(x^2-3); alpha:=Elt(O,[0,1]); u3:=alpha; u2:=1+alpha; p2:=u2*O; Plc2:=PlaceInit(p2); p23:=[(-3*alpha-2)*O,(3*alpha-2)*O]; Plc23:=List(p23,u->PlaceInit(u)); ec:=EllInit(alpha,0,0,0,-2*alpha,0); Gs:=[[0,0,1]]; Gps:=[[6+6*alpha,24+16*alpha,1],[0,0,1]]; EllGenInit(ec,[EllPsip(ec,Gps[1]),Gs[1]],1); EllPhiInit(ec,[1,0,0],[0,0,1]); fiblist:=List([[0,0],[0,1]],v->EllVecToPnt(ec,v)); PlcList:=Plc23; TorPlc:=[Plc23[1],PlaceInit(Ideal(13,alpha-4))]; GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #--- eq283E4.g (end) ---# #--- eq283E7.g (start) ---# cat << "#--- eq283E7.g (end) ---#" > eq283E7.g # UNIX extraction command Read("isogdesc.g"); O:=OrderMaximal(x^4-3); beta:=Elt(O,[0,1,0,0]); u3:=beta; u2:=1+beta; p2:=u2*O; Plc2:=PlaceInit(p2); p13:=List([beta-3,beta+2,beta+3,beta-2],u->Ideal(13,u)); Plc13:=List(p13,u->PlaceInit(u)); ec:=EllInit(-(beta^2+beta),0,2,0,2,0); Gs:=[[-3-2*beta-2*beta^2-beta^3/3,2+5*beta+5*beta^2/3+2*beta^3,1],[0,0,1]]; Gps:=[[-2*beta-2*beta^2,2+4*beta+2*beta^2,1],[0,0,1]]; EllGenInit(ec,[Gs[1],EllPsip(ec,Gps[1]),Gs[2]],2); EllPhiInit(ec,[beta,0,beta],[0,0,1]); fiblist:=List([[0,0,0],[1,-1,0]],v->EllVecToPnt(ec,v)); PlcList:=Plc13; TorPlc:=[Plc13[1],PlaceInit(Ideal(11,beta-4))]; GlobalCert(ec,Gs,Gps,SquareBasis(PlaceSupport(EllDisc(ec)))); for Plc in PlcList do ecp:=EllRed(ec,Plc); Print("Make sure that [E(K):Gs] is prime to ", EllSize(ecp)/MatDet(EllKernel(ecp)),"\n"); od; Print("Torsion bounded by ", EllBoundTor(ec,TorPlc),"\n"); EllChab(ec,PlcList,fiblist); #--- eq283E7.g (end) ---# #--- isogdesc.g (start) ---# cat << "#--- isogdesc.g (end) ---#" > isogdesc.g # UNIX extraction command ############################################################################### ############################################################################### ## ## isogdesc.g ## ## Nils Bruin, 4 april 1999 ## ## Routines to do a 2-isogeny descent on elliptic curves of the form ## gamma*Y^2=X^3+a2*X^2+a4*X. ## ############################################################################### ############################################################################### Read("ellcalc.g"); Read("loccalc.g"); ############################################################################### EllPsi:=function(ec,G) # 2-isogeny to 2-isogeneous curve. if G[3]=0 then return EllZero(ec); else return EllPointNormal(ec, [ec.gamma*G[2]^2*G[3],G[2]*(ec.a4*G[3]^2-G[1]^2),G[1]^2*G[3]]); fi; end; ############################################################################### EllPsip:=function(ec,G) # dual 2-isogeny (to curve itself) if G[3]=0 then return EllZero(ec); else return EllPointNormal(ec, [2*ec.gamma*G[2]^2*G[3], G[2]*((ec.a2^2-4*ec.a4)*G[3]^2-G[1]^2),8*G[1]^2*G[3]]); fi; end; ############################################################################### EllMu:=function(ec,G) # group homomorphism E(K)/psi'(E'(K)) -> K^*/(K^*)^2 if G[3]=0 then return RingOne(ec.R); elif G[1]=0 then return ec.a4; else return ec.gamma*G[1]/G[3]; fi; end; ############################################################################### EllHomSpcCfs:=function(ec,delta) # Returns[f0,f2,f4], where Y^2=f0*X^4+f2*X^2+f4 is the homogeneous space # belonging to given delta. if ec.a6<>RingZero(ec.R) then Error("(0,0) not on curve"); fi; return [ec.gamma^2*delta^3,ec.gamma*delta^2*ec.a2,delta*ec.a4]; end; ############################################################################### LocCertInit:=function(Plc,ec,Gs,addX,Gsp,addXp) # Initialises a LocCert data-structure. This is a record containing data # describing E/psi'E and E'/psi E locally at Plc. # Gs, Gsp - global points on ec and 2-isogeneous curve resp. # addX, addXp - X - coordinates of additional local points. # # returns a record with fields: # Plc - Exactly given Plc. # ec, ecp - the two isogeneous curves (ec is identical to given ec) # Gs, Gsp - given Gs and Gsp # gens, gensp - X-coordinates of (local) points independent mod image(psi'), # mod image(psi) resp. # muimg,muimgp - PlaceSqrRep of EllMu of gens, gensp. local ecp,Cert,Certp,muimg,muimgp,G,t,xc; if ec.a6<>RingZero(ec.R) then Error("2-Isog. certificates only valid for a6=0"); fi; ecp:=EllInit(ec.gamma,0,-2*ec.a2,0,ec.a2^2-4*ec.a4,0); Cert:=[]; #initialise muimg muimg:=MatBasis(Mat(MatRing,[PlaceSqrRep(Plc,EllMu(ec,[0,1,0]))])); #loop through known generators and add to certificate if relevant for G in Gs do if not(EllPointOnCurve(ec,G)) then Error("Point in Gs not on curve"); fi; t:=MatBasis(MatUnion(muimg,Mat(MatRing,[PlaceSqrRep(Plc,EllMu(ec,G))]))); if MatRows(t)>MatRows(muimg) then Add(Cert,G[1]/G[3]); muimg:=t; fi; od; #loop through additional points (check if they really are X of local points) #add if relevant. for xc in addX do if PlaceIsSqr(Plc,EllYsqr(ec,xc)) then t:=MatBasis(MatUnion(muimg,Mat(MatRing,[PlaceSqrRep(Plc,EllMu(ec,[xc,0,1]))]))); if MatRows(t)>MatRows(muimg) then Add(Cert,xc); muimg:=t; fi; else Print("WARNING: ",xc," not X for a local point; skipping...\n"); fi; od; #same for ecp. Certp:=[]; muimgp:=MatBasis(Mat(MatRing,[PlaceSqrRep(Plc,EllMu(ecp,[0,1,0]))])); for G in Gsp do if not(EllPointOnCurve(ecp,G)) then Error("Point in Gsp not con curve"); fi; t:=MatBasis(MatUnion(muimgp,Mat(MatRing,[PlaceSqrRep(Plc,EllMu(ecp,G))]))); if MatRows(t)>MatRows(muimgp) then Add(Certp,G[1]/G[3]); muimgp:=t; fi; od; for xc in addXp do if PlaceIsSqr(Plc,EllYsqr(ecp,xc)) then t:=MatBasis(MatUnion(muimgp,Mat(MatRing,[PlaceSqrRep(Plc,EllMu(ecp,[xc,0,1]))]))); if MatRows(t)>MatRows(muimgp) then Add(Certp,xc); muimgp:=t; fi; else Print("WARNING: ",xc," not X for a local point; skipping...\n"); fi; od; return rec( Plc:=Plc, ec:=ec, ecp:=ecp, Gs:=Gs, Gsp:=Gsp, gens:=Cert, gensp:=Certp, muimg:=muimg, muimgp:=muimgp ); end; ############################################################################### LocCertIsComplete:=function(lct) # Test whether a LocCert is complete, i.e. describes the complete local # images of mu. (i.e. dimensions should add to 2+valuation(2) ) local dimsum,expdim; dimsum:=Length(lct.gens)+Length(lct.gensp); expdim:=2+IdealDegree(PlacePrime(lct.Plc))*PlaceVal(lct.Plc,2); if dimsum < expdim then return false; elif dimsum = expdim then return true; else Error("Dimension mismatch"); fi; end; ############################################################################### HomSpcLocalPoint:=function(Plc,cfs) # Returns an X such that cfs[1]*X^4+cfs[2]*X^2+cfs[3] is a local square # or "false" if it does not exist. local A,B,C,skew,e,L,tmp; # first, shift such that leading coeff. has valuation 0 or 1 e:=2*Floor(PlaceVal(Plc,cfs[1])/2); C:=cfs[1]/Plc.unif^e; A:=cfs[2]/Plc.unif^e; B:=cfs[3]/Plc.unif^e; # flatten Newton polygon (using scalar mult. of X) skew:=Minimum(Floor(PlaceVal(Plc,B)/4),Floor(PlaceVal(Plc,A)/2)); A:=A/Plc.unif^(2*skew); B:=B/Plc.unif^(4*skew); # make sure coeffs. are integral. e:=2*Floor(Minimum(PlaceVal(Plc,C),PlaceVal(Plc,A),PlaceVal(Plc,B))/2); C:=C/Plc.unif^e;A:=A/Plc.unif^e;B:=B/Plc.unif^e; # fill data structure for recursive search HSL:=rec(O:=PlaceOrder(Plc), C:=C, A:=A, B:=B, val2:=PlaceVal(Plc,2),val3:=PlaceVal(Plc,3), valA:=PlaceVal(Plc,A), valB:=PlaceVal(Plc,B), valC:=PlaceVal(Plc,C), oneonly:=true, Plc:=Plc, bounds:=IdealRepBounds(PlacePrime(Plc))); # call filtering routine. (and return found point, if any). L:=HomSpcLocalPoints_rec(0,0); if L<>[] then return L[1][1]*Plc.unif^skew; fi; # reverse coefficients tmp:=HSL.C;HSL.C:=HSL.B;HSL.B:=tmp; tmp:=HSL.valC;HSL.valC:=HSL.valB;HSL.valB:=tmp; # search around X=0 (was: X=infinity before reversal) L:=HomSpcLocalPoints_rec(0,1); if L<>[] then return Plc.unif^skew/L[1][1]; else return false; fi; end; ############################################################################### LocCertComplete:=function(lct) # try to complete a LocCert, if necessary. local Plc,Rep,mc,t,L,delta; if LocCertIsComplete(lct) then Print("Local Certificate already is complete.\n"); return lct; fi; Print("Not enough points. Trying to complete certificate...\n"); #loop through all deltas (mod squares) Plc:=lct.Plc; Rep:=PlaceSqrRepGens(Plc); mc:=MultiCounterInit(Length(Rep)); while MultiCounterInc(mc,1) do #note that PlaceSqrRep(delta)=mc by construction. t:=MatBasis(MatUnion(lct.muimg,Mat(MatRing,[mc]))); #if value might increase dimension of muimg, then investigate. if MatRows(t)>MatRows(lct.muimg) then delta:=Product([1..Length(mc)],i->Rep[i]^mc[i]); L:=HomSpcLocalPoint(Plc,EllHomSpcCfs(lct.ec,delta)); # if point is found, then add to data structures. if L<>false then Add(lct.gens,lct.ec.gamma*delta*L^2); lct.muimg:=t; fi; fi; t:=MatBasis(MatUnion(lct.muimgp,Mat(MatRing,[mc]))); #if value might increase dimension of muimgp, then investigate. if MatRows(t)>MatRows(lct.muimgp) then delta:=Product([1..Length(mc)],i->Rep[i]^mc[i]); L:=HomSpcLocalPoint(Plc,EllHomSpcCfs(lct.ecp,delta)); if L<>false then Add(lct.gensp,lct.ecp.gamma*delta*L^2); lct.muimgp:=t; fi; fi; od; if LocCertIsComplete(lct) then Print("Local Certificate is now complete.\n"); return lct; else Error("Could not complete certificate."); fi; end; ############################################################################### DescInfMuImg:=function(ec) # returns a list [L,Lp], where L gives the real embeddings (as in EltCon) # where EllMu(ec,G) is positive (due to real considerations). Lp gives the # same for the isogeneous curve. Note that L and Lp form a partition of the # real embeddings. (some basic real algebra) local L1,L2,i; L1:=[];L2:=[]; for i in [1..OrderSig(ec.R)[1]] do if Re(EltCon(ec.a2^2-4*ec.a4,i))<0 or (Re(EltCon(ec.gamma*ec.a2,i))<0 and Re(EltCon(ec.a4,i))>0) then Add(L1,i); else Add(L2,i); fi; od; return([L1,L2]); end; ############################################################################### SqrSpan:=function(L) # returns the nontrivial multiplicative combinations of L. local mc,R; R:=[]; mc:=List([1..Length(L)],i->0); while MultiCounterInc(mc,1) do Add(R,Product([1..Length(L)],i->L[i]^mc[i])); od; return R; end; ############################################################################### GlobalCertHints:=[]; GlobalCertLookUpHint:=function(ec,plc) local L; for L in GlobalCertHints do if L[1]=ec and PlacePrime(plc)=PlacePrime(L[2]) then return L[3]; fi; od; return false; end; GlobalCertAddHint:=function(ec,plc,AddX,AddXp) if GlobalCertLookUpHint(ec,plc)<>false then Error("Hint already present for this curve and place"); fi; Add(GlobalCertHints,[ec,plc,[AddX,AddXp]]); end; GlobalCert:=function(ec,Gs,Gsp,K2S) # does a 2-isogeny descent on ec and returns a certificate of it. # K2S should span K(S,2). # Gs and Gsp might be modified: if [0,0,1] on E and E' help, then this is added. # furthermore, it is checked if Gs and Gsp span E(K)/psi'E'(K) and # E'(K)/psi E(K). A data structure describing all this is returned, containing: # ec, ecp, K2S, # lcts - A list of the LocCert's used in the descents. # gens, gensp - the independent parts from Gs and Gsp (if useful, points with # X=0 are added to these) # muimg, muimgp - image as F2-basis on K2S (note that real considerations are # also in here. These cannot be found in lcts.) local ecp,M1,M2,gens,gensp,G,L,i,lcts,lct,Loc,Plc,S; if ec.a6<>RingZero(ec.R) then Error("2-isogeny only valid for a6=0"); fi; Print("Doing a 2-isogeny descent to determine selmer-rank...\n"); gens:=[]; M1:=[]; #if [0,0,1] is not in Gs, then add it (if an equivalent generator is # present, then this will not matter in gens, since gens is built up # iteratively) if ForAll(Gs,u->u[1]<>0) and PolyIsIrreducible(x^2-ec.a4) then Print("a4 is not a square: adding (0,0) to generators of ec.\n"); Add(Gs,EllPointNormal(ec,[0,0,1])); fi; #take independent G's in Gs for G in Gs do if not EllPointOnCurve(ec,G) then Error("Point not on curve."); fi; Add(M1,EllMu(ec,G)); if ForAll(SqrSpan(M1),u->PolyIsIrreducible(x^2-u)) then Add(gens,G); else Unbind(M1[Length(M1)]); fi; od; # same for isog. curve. ecp:=EllInit(ec.gamma,0,-2*ec.a2,0,ec.a2^2-4*ec.a4,0); gensp:=[]; M2:=[]; if ForAll(Gsp,u->u[1]<>0) and PolyIsIrreducible(x^2-ecp.a4) then Print("a2^2-4*a4 is not a square: adding (0,0) to generators of ec.\n"); Add(Gsp,EllPointNormal(ecp,[0,0,1])); fi; for G in Gsp do if not EllPointOnCurve(ecp,G) then Error("Point not on curve."); fi; Add(M2,EllMu(ecp,G)); if ForAll(SqrSpan(M2),u->PolyIsIrreducible(x^2-u)) then Add(gensp,G); else Unbind(M2[Length(M2)]); fi; od; # This gives us a lower bound on the rank. This might be negative if # very few points are given and a4 or a2^2-4*a4 is a square. Print("Rank lower bound due to given generators:", Length(M1)+Length(M2)-2,"\n"); M1:=MatId(MatRing,Length(K2S)); M2:=MatId(MatRing,Length(K2S)); Print("A priori rank upper bound:",MatRows(M1)+MatRows(M2)-2,"\n"); #let's intersect at infinite places. L:=DescInfMuImg(ec); for i in L[1] do M1:=MatIntersect(M1,MatKernel(Mat(MatRing,List(K2S,u->InfSqrRep(i,u))))); od; for i in L[2] do M2:=MatIntersect(M2,MatKernel(Mat(MatRing,List(K2S,u->InfSqrRep(i,u))))); od; Print("Rank upper bound after real considerations:", MatRows(M1)+MatRows(M2)-2,"\n"); #It seems like a good idea to use all places of bad reduction (plus some, #if K2S needs a larger support (due to class group). S:=List(Factor(Product(K2S*EllDisc(ec))),u->PlaceInit(u[1])); lcts:=[]; for Plc in S do Print("pondering on ",PlacePrime(Plc)," (above ",Plc.primebelow,")...\n"); i:=GlobalCertLookUpHint(ec,Plc); if i=false then i:=[[],[]]; else Print("Using hints for additional local points.\n"); fi; lct:=LocCertComplete(LocCertInit(Plc,ec,Gs,i[1],Gsp,i[2])); Add(lcts,lct); Loc:=Mat(MatRing,List(K2S,u->PlaceSqrRep(Plc,u))); M1:=MatIntersect(M1,MatKernel(Loc*MatTrans(MatOrth(lct.muimg)))); M2:=MatIntersect(M2,MatKernel(Loc*MatTrans(MatOrth(lct.muimgp)))); Print("Rank upper bound:",MatRows(M1)+MatRows(M2)-2,"\n"); od; if Length(gens)=MatRows(M1) then Print("gens generate image of mu^(\psi').\n"); elif Length(gens)>MatRows(M1) then Error("too many gens !"); else Print("Still room left on ec\n"); fi; if Length(gensp)=MatRows(M2) then Print("gensp generate image of mu^(\psi).\n"); elif Length(gensp)>MatRows(M2) then Error("too many gens !"); else Print("Still room left on ecp\n"); fi; return rec( ec:=ec, ecp:=ecp, K2S:=K2S, muimg:=M1, muimgp:=M2, lcts:=lcts, gens:=gens, gensp:=gensp); end; ############################################################################### EllBoundTor:=function(ec,PlcList) local Bound,Plc; Bound:=0; for Plc in PlcList do if Plc.ramdeg>=Plc.primebelow-1 then Error("Hey, we're lazy here. We don't do highly ramified primes."); fi; Bound:=Gcd(Bound,EllSize(EllRed(ec,Plc))); od; return Bound; end; # end of isogdesc.g #--- isogdesc.g (end) ---# #--- latcalc.g (start) ---# cat << "#--- latcalc.g (end) ---#" > latcalc.g # UNIX extraction command ############################################################################### ############################################################################### ## ## latcalc.g ## ## Nils Bruin, 5 april 1999 ## ## Routines to compute intersections and sums of lattices (subspaces of Z^r) ## Lattices are represented as rowspans of matrices. ## ############################################################################### ############################################################################### LatHermite:=M->MatTrans(MatHermiteColLower(MatTrans(M))); # HermiteRowUpper. This is missing from KASH, which is a little strange since # MatKernel does work with respect to rows. ############################################################################### LatSum:=function(M1,M2) # Compute the joint span of M1 and M2 in the following way: # stack M1 and M2, compute Hermite normal form. The non-zero rows # generate the Z-span of M1 and M2. local L,R,ZeroRow,v; L:=MatToRowList(M1); ZeroRow:=List(L[1],u->0); Append(L,MatToRowList(M2)); R:=[]; for v in MatToRowList(LatHermite(Mat(Z,L))) do if v<>ZeroRow then Add(R,v); fi; od; return Mat(Z,R); end; ############################################################################### LatSectExt:=function(M1,M2) # Compute intersection of 2 lattices and return representations with # respect to the bases for M1 and M2. in the following way: # Stack M1 and -M2 and compute the kernel M3 of this matrix. The first # Rows(M1) columns span the linear combinations of M1 that are also # expressible with respect to M2 (thus, the intersection). The same for the # last Rows(M2) columns of M3. local L,M3,I1,I2; L:=MatToRowList(M1); Append(L,MatToRowList(M2)); M3:=MatKernel(Mat(Z,L)); L:=MatToColList(M3); I1:=MatTrans(Mat(Z,L{[1..MatRows(M1)]})); I2:=MatTrans(Mat(Z,L{[MatRows(M1)+1..Length(L)]})); return([I1,I2]); end; ############################################################################### LatContainsVec:=function(M,v) # Determines if v is in M by seeing if adding v to M changes the lattice #(ie. if the HNF changes). local Ma,Mp; Mp:=LatHermite(M); Ma:=LatSum(Mp,Mat(Z,[v])); return (Ma=Mp); end; ############################################################################### LatSect:=function(M1,M2) local L; L:=LatSectExt(M1,M2); return LatHermite(L[1]*M1); end; ############################################################################### LatRelBasis:=function(M1,M2) # Expresses the basis of M1 with respect to a sublattice M2. M2 should be of # maximal rank in ambient space. local R,L; R:=M1*M2^(-1); L:=MatToRowList(R); if ForAny(Flat(L),u->Denominator(u)<>1) then Error("Can only calculate relative bases for sublattices"); fi; return Mat(Z,L); end; # end of latcalc.g #--- latcalc.g (end) ---# #--- loccalc.g (start) ---# cat << "#--- loccalc.g (end) ---#" > loccalc.g # UNIX extraction command ############################################################################### ############################################################################### ## ## loccalc.g ## ## Nils Bruin, 5 april 1999 ## ## Routines to do calculations in "local fields" ## ## We introduce a new data structure, a "place", which is a place holder for ## information at a prime ideal and we give some routines to do calculations ## with them ## ## Beware that this file sets ILS_zero_is_precise=true, which changes the ## semantics of IsLocSqr to a less intuitive one, but increases the performance ## of FindLocRatX. ## ############################################################################### ############################################################################### Read("matalg.g"); MatRing:=FF(2); #some instructions for KASH 1.8 compatibility #FFEltFF:=FFEFF; #FFEltLog:=function(t) # local w,t,u,l; # w:=FFPrimitiveElt(FFEltFF(t)); # if t=w-w then # Error("FFEltLog(0) undefined"); # fi; # u:=w^0;l:=0; # while u<>t do # u:=w*u;l:=l+1; # od; # return(l); #end; ############################################################################### IdealRepBounds:=function(I) # returns a list to be used in combination with the MultiCounter routines to # enumerate a full set of representatives of O/I local B; B:=IdealBasisUpperHNF(I)[2]; return List([1..MatRows(B)],i->MatElt(B,i,i)-1); end; ############################################################################### PlaceInit:=function(prime) # returns a place-structure belonging to the given prime ideal. local u,primebelow,ramdeg; u:=IdealIsPrincipal(prime); if u=false then Error("Cannot work with non-principal ideals yet"); fi; primebelow:=IdealMin(prime); ramdeg:=Valuation(prime,primebelow); return rec(unif:=u, filtr:=[prime], val4p1:=Valuation(prime,4)+1, primebelow:=primebelow, ramdeg:=ramdeg); end; ############################################################################### PlaceGetIdeal:=function(place,i) # returns the i-th power of the prime ideal belonging to place. When called # twice, it will return the same ideal (ie. a pointer to the same structure) # This mains that a basis for (O/p^i)^* remains fixed. local j; for j in [Length(place.filtr)+1..i] do place.filtr[j]:=place.filtr[1]^j; od; return place.filtr[i]; end; ############################################################################### PlacePrime:=(place)->place.filtr[1]; # shorthand for PlaceGetIdeal(place,1) ############################################################################### PlaceValZeroVal:=32767; PlaceVal:=function(place,elt) # valuation of elt at prime belonging to place. # if elt=0 then a preset "infinity" value is returned. if 2*elt=elt then return PlaceValZeroVal; else return Valuation(PlacePrime(place),elt); fi; end; ############################################################################### PlaceOrder:=(place)->IdealOrder(place.filtr[1]); # underlying order of place ############################################################################### PlaceSqrRep:=function(place,u) # represents the class of u in O_p^*/(O_p^*)^2 as a vector of an F2-vector # space. After place=PlaceInit(prime), the choice of basis is fixed, but there # is no guarantee that two runs PlaceInit's result in the same choice of # basis. The result is returns as a list of integers from {0,1}. The first # component always is the parity of the valuation. local e,v,w,i,res; #we make sure that u is a local integer by multiplying with a square. #(in fact, we make sure that u is an integer at all places of same #characteristic. This is necessary for "mod" in kash) #WARNING! patch for older versions of KASH. Denominators are not handled #correctly by those versions, so we remove any. u:=u*EltDen(u)^2; # end of patch v:=Valuation(place.primebelow,EltDen(u)); if v>0 then u:=u*place.primebelow^(2*Ceil(v/2)); fi; e:=Valuation(place.filtr[1],u); u:=u/place.unif^e; v:=MatToRowList(EltRayResidueRingRep(u,PlaceGetIdeal(place,place.val4p1)))[1]; w:=RayResidueRing(PlaceGetIdeal(place,place.val4p1))[2]; res:=[e mod 2]; for i in [1..Length(w)] do if (w[i] mod 2)=0 then Add(res,v[i] mod 2); fi; od; return res; end; ############################################################################### PlaceSqrRepGens:=function(place) # returns a lift of a basis for K_p^*/(K_p^*)^2 (elements of O_K) local L,u; L:=[place.unif]; for u in RayResidueRingCyclicFactors(PlaceGetIdeal(place,place.val4p1)) do if u[2] mod 2 = 0 then Add(L,u[1]); fi; od; return L; end; ############################################################################### InfSqrRep:=function(i,u) # The same as PlaceSqrRep, but at the i-th real place of EltOrder(u) local i; if i>OrderSig(EltOrder(u))[1] then return [0]; elif Re(EltCon(u,i))<0 then return [1]; else return [0]; fi; end; ############################################################################### PlaceIsSqr:=function(place,u) # Determines whether u is a square in the completion at place. if u=0 then return true; else return ForAll(PlaceSqrRep(place,u),v->v=0); fi; end; ############################################################################### ILS_zero_is_precise:=false; IsLocSqr:=function(place,level,u) # Tries to answer the question whether u is a local square given only # limited precision. Thus, u is only determined up to prime^level. # given this, the routine returns # -1 if the class u mod p^level does not contain local squares # 0 if the class u mod p^level contains both local squares and non-squares # 1 if the class u mod p^level consists of local squares. # if ILS_zero_is_precise is false, then the behaviour is changed for u=0. # This is considered to be an exact value (level=infinity), so 1 is returned. # This behaviour is useful for FindLocalRatX. local e,v,w,i,lvl; # handle u=0 separately if u=0 then if ILS_zero_is_precise then return 1; else return 0; fi; fi; # too close to 0 means: not enough precision to make a decision. e:=Valuation(place.filtr[1],u); if e>=level then # too close to 0 means: not enough precision to make a decision. return 0; elif e mod 2=1 then # odd valuation means: certainly not a square. return -1; else #by now, we divide by a square to make u a local unit. u:=u/place.unif^e; #of other precision lvl:=Minimum(level-e,place.val4p1); #look at the exponents in the multiplicative group v:=MatToRowList(EltRayResidueRingRep(u,PlaceGetIdeal(place,lvl)))[1]; w:=RayResidueRing(PlaceGetIdeal(place,lvl))[2]; #if one of the exponents on a generator of even order is odd, then #we don't have a square. for i in [1..Length(w)] do if (w[i] mod 2)=0 and (v[i] mod 2)=1 then return -1; fi; od; #otherwise, we do have a square on this level. It depends on the level #if this can change for higher levels. if lvl>=place.val4p1 then return 1; fi; return 0; fi; end; ############################################################################### IntIsLocSqr:=function(p,level,u) local e,lvl; if u=0 then if ILS_zero_is_precise then return 1; else return 0; fi; fi; e:=Valuation(p,u); if e>=level then return 0; elif e mod 2 =1 then return 1; else u:=u/p^e; if p=2 then lvl:=Minimum(level-e,3); if (u mod 2^lvl)=1 then if lvl=3 then return 1; else return 0; fi; else return -1; fi; else if (FFEltLog(FFElt(FF(p),u)) mod 2)=0 then return 1; else return -1; fi; fi; fi; end; ############################################################################### FFEToInt:=function(elt) # This routine lifts an element from a primitive finite field # to a non negative integer. local F,sz,i; F:=FFEltFF(elt); sz:=FFSize(F); if sz[2]>1 then Error("Only implemented for primitive fields"); fi; for i in [0..sz[1]-1] do if FFElt(F,i)=elt then return i; fi; od; Error("Couldn't find integral representative"); end; #FFEToInt := function (fe) # # suggested routine by KASH-team. Should be more efficient, but causes SEGV # after a while # #local str, feZ; # str := SPrint(fe); # #Print (str,"\n"); # feZ := SScan(str, "%Z")[1]; # return feZ; #end; ############################################################################### FLRX:=rec(); FLRX_rec:=function(center,level) # recursive part of FindLocalRatX. Shouldn't be called directly. local L,i,cnd,value,sqrdtl; L:=[]; # loop through candidates around center. for i in [0..FLRX.p-1] do cnd:=center+Elt(FLRX.O,i*FLRX.p^(level/FLRX.step)); # determine the value value:=Eval(FLRX.poly,cnd); # and determine the nature of this value. sqrdtl:=IsLocSqr(FLRX.Plc,level+FLRX.step,value); if sqrdtl=1 then # whole environment of cnd suffices, apparently. Add(L,[cnd,level+FLRX.step]); elif sqrdtl=0 then # not enough precision yet. Look more closely. Append(L,FLRX_rec(cnd,level+FLRX.step)); fi; # sqrdtl=-1 means: no points to be found around this cnd. od; return L; end; ############################################################################### ILS_zero_is_precise:=true; FindLocalRatX:=function(plc,poly) # Tries to find rational values for X such that poly(X) is a local square at # plc. Returns a list in the form [...,[u,e],...], such that u mod prime^e is # a suitable environment for X. If e<0, then 1/X should be in u mod prime^(-e) # routine fails if poly has a zero in Q_p that is not in Q. # # If ILS_zero_is_precise=true then zeroes in Q are handled correctly, # because in FLRX_rec, only exact values are fed to IsLocSqr, which will # then consider 0 to be a square regardless of level. local L,L1,i; #initialise data structure FLRX:=rec( O:=PlaceOrder(plc), poly:=poly, Plc:=plc, step:=plc.ramdeg, p:=plc.primebelow); #and search away from infinity L:=FLRX_rec(0,0); #exchange 0 and infinity (i.e. reverse polynomial) FLRX.poly:=Poly(PolyAlg(poly),Reversed(PolyToList(poly))); #and search around 0 L1:=FLRX_rec(0,FLRX.step); for i in [1..Length(L1)] do L1[i][2]:=-L1[i][2]; od; Append(L,L1); return L; end; ############################################################################### HSL:=rec(); HomSpcLocalPoints_rec:=function(base,level) local L, mc, vallevel, cnd, value, sqrdtl; L:=[]; mc:=MultiCounterInit(Length(HSL.bounds)); vallevel:=level+Minimum(2*HSL.val2+HSL.valC, HSL.val2+HSL.val3+HSL.valC+level, HSL.valC+3*level, HSL.val2+HSL.valA, HSL.valA+level); # vallevel:=Minimum( # level+HSL.valC+Minimum(2*HSL.val2,HSL.val2+HSL.val3+level,3*level), # HSL.valA+Minimum(HSL.val2+level,2*level)); while MultiCounterInc(mc,HSL.bounds) do cnd:=(base+Elt(HSL.O,mc)*HSL.Plc.unif^level) mod PlaceGetIdeal(HSL.Plc,level+1); value:=HSL.C*cnd^4+HSL.A*cnd^2+HSL.B; sqrdtl:=IsLocSqr(HSL.Plc,vallevel,value); if sqrdtl=1 then Add(L,[cnd,level+1]); if HSL.oneonly then return L; fi; elif sqrdtl=0 then Append(L,HomSpcLocalPoints_rec(cnd,level+1)); if HSL.oneonly and L<>[] and (L[1][1]<>0 or Length(L)>1) then return L; fi; fi; od; return L; end; ############################################################################### HomSpcLocalPoints:=function(plc,cfs,oneflag) local skew,L,a,b,c,F,cnd,prelt,i,prm,t,F0; c:=cfs[1];a:=cfs[2];b:=cfs[3]; if 2*a=a then skew:=Floor(PlaceVal(plc,b)/4); elif 2*b=b then skew:=Floor(PlaceVal(plc,a)/2); else skew:=Floor(Minimum(PlaceVal(plc,a)/2,PlaceVal(plc,b)/4)); fi; HSL:=rec(O:=PlaceOrder(plc), C:=c,A:=a/plc.unif^(2*skew),B:=b/plc.unif^(4*skew), val2:=PlaceVal(plc,2),val3:=PlaceVal(plc,3), valC:=PlaceVal(plc,c), oneonly:=oneflag, Plc:=plc, bounds:=IdealRepBounds(PlacePrime(plc))); if 2*a=a then HSL.valA:=1000000; else HSL.valA:=PlaceVal(plc,a)-2*skew; fi; if 2*b=b then HSL.valB:=1000000; else HSL.valB:=PlaceVal(plc,b)-4*skew; fi; if HSL.val2=0 and PlaceVal(plc,PolyDisc(HSL.C*x^4+HSL.A*x^2+HSL.B))=0 then prm:=PlacePrime(plc); F:=IdealResidueField(prm); a:=EltToFFE(HSL.A,prm);b:=EltToFFE(HSL.B,prm);c:=EltToFFE(HSL.C,prm); L:=[]; if FFEltLog(b) mod 2=0 then Add(L,[Elt(HSL.O,0),1]); fi; cnd:=FFElt(F,1); prelt:=FFPrimitiveElt(F); F0:=prelt-prelt; for i in [1..FFSize(F)[1]^FFSize(F)[2]-1] do t:=c*cnd^4+a*cnd^2+b; if t=F0 or FFEltLog(t) mod 2 = 0 then Add(L,[FFEToElt(cnd,prm),1]); fi; cnd:=cnd*prelt; od; else L:=HomSpcLocalPoints_rec(0,0); fi; return List(L,u->[u[1]*plc.unif^skew mod PlaceGetIdeal(plc,u[2]+skew),u[2]+skew]); end; ############################################################################### SelectPrimPlaces:=function(O,nplaces,D) local p,L,fct,f; p:=2; L:=[]; while Length(L)EltToFFE(u,prm))); L:=[]; Print(plc.primebelow,"\n"); for v in V do Print(FFEToInt(EltToFFE(v,prm))," "); od; Print("\n",FFEToInt(EltToFFE(v0,prm)),"\n"); for xc in [0..plc.primebelow-1] do xm:=FFElt(F,xc); t:=Eval(polym,xm); Add(L,t=F0 or FFEltLog(t) mod 2 =0); if (t=F0 or FFEltLog(t) mod 2 =0) then Print("1"); else Print("0"); fi; od; Print("\n"); od; Print("0\n"); end; WriteGenInfo:=function(FILE,poly,Bound,PlcList) ECHOon(FILE); PrintGenInfo(poly,Bound,PlcList); ECHOoff(); end; CommLine:=[]; DoGen:=function(poly,Bound,Nfilters) local PlcList,O; O:=PolyAlgCoef(PolyAlg(poly)); Print("determining places to filter with ...\n"); PlcList:=SelectPrimPlaces(O,Nfilters,2*PolyDisc(poly)); Print("sieving ...\n"); WriteGenInfo("/tmp/dogen.in",poly,Bound,PlcList); Exec("/home/nbruin/cover/filter /tmp/dogen.out"); Read("/tmp/dogen.out"); CommLine:=List(CommLine{[1..Length(CommLine)-1]}, u->Elt(O,u-Bound)); return Filtered(CommLine, u->not(PolyIsIrreducible(x^2-Eval(poly,u)))); end; ############################################################################### PlaceSupport:=function(object) #Returns a list of the places at which "object" is not a unit. #object should be something for which Factor returns an ideal factorisation #(i.e. an Elt or an Ideal) or a list of such things. if IsList(object) then object:=Product(object); fi; return List(Factor(object),u->PlaceInit(u[1])); end; ############################################################################### PlacePointNormal:=function(Plc,G) # brings a point in padic normal form. local e; e:=Maximum(List(G,u->Valuation(Plc.primebelow,EltDen(u)))); G:=List(G,u->u*Plc.primebelow^e); e:=Minimum(List(G,u->PlaceVal(Plc,u))); G:=List(G,u->u/Plc.unif^e); return G; end; ############################################################################### SquareBasis:=function(Plcs) local O,L; O:=EltOrder(Plcs[1].unif); Print("Computing class group...\n"); if OrderClassGroup(O)[1]mod 2 = 0 then Error("Can't do SquareBasis for even class number yet."); fi; L:=[OrderTorsionUnit(O)]; Print("Computing fundamental units...\n"); Append(L,OrderUnitsFund(O)); Append(L,List(Plcs,u->u.unif)); return L; end; # end of loccalc.g #--- loccalc.g (end) ---# #--- matalg.g (start) ---# cat << "#--- matalg.g (end) ---#" > matalg.g # UNIX extraction command ############################################################################### ############################################################################### ## ## matalg.g ## ## for KASH version 1.8 ## ## Version 1.0, January 27, 1998 ## ## Nils Bruin, nbruin@wi.leidenuniv.nl ## ## Description: Extended support for linear algebra. ## ## Convention: Bases of subspaces are represented as row vectors of a matrix. ## ## (C) This code is regarded as a scientific publication. Work (partially) ## based on results obtained using this code should contain an appropriate ## reference. You are free to copy this file unaltered for non-commercial ## purposes. ## This code is provided as is, with no claim of fitness for any particular ## purpose. Although great care has been taken to eliminate all bugs, there is ## no guarantee that results produced are correct. ## ## Note that similar restrictions hold for KASH, the platform this file is ## designed for. ############################################################################### ############################################################################### # # Be sure to set MatRing to the ring over which you want to do linear algebra. # Some routines need this. # if not(IsBound(MatRing)) then MatRing:=Z; fi; ############################################################################### # # MatOrth(A) # # Returns a basis of the orthogonal space ############################################################################### MatOrth:=A->MatKernel(MatTrans(A)); ############################################################################### # # MatUnion(A,B) # # Returns a matrix with the rows of A and B. # This routine needs MatRing being set correctly. ############################################################################### MatUnion:=function(A,B) local L; if MatRows(A)=0 then return B; elif MatRows(B)=0 then return A; fi; L:=MatToRowList(A); Append(L,MatToRowList(B)); return(Mat(MatRing,L)); end; ############################################################################### # # MatBasis(A) # # Constructs a matrix with same rowspace of maximal rank ############################################################################### MatBasis:=A->MatOrth(MatOrth(A)); ############################################################################### # # MatIntersect # # Returns a basis of the intersection of A and B. # This routine needs MatRing being set correctly. ############################################################################### MatIntersect:=function(A,B) return(MatOrth(MatUnion(MatOrth(A),MatOrth(B)))); end; ############################################################################### # # MatCatDiag(A,B) # # Makes a matrix [[A,O],[O,B]], where O are zero-matrices of appropriate size # This routine needs MatRing being set correctly. ############################################################################### MatCatDiag:=function(A,B) local Ac,Bc,l,Al,t; Ac:=MatCols(A); Bc:=MatCols(B); Al:=[]; for l in MatToRowList(A) do Append(l,List([1..Bc],u->0)); Add(Al,l); od; for l in MatToRowList(B) do t:=List([1..Ac],u->0); Append(t,l); Add(Al,t); od; return(Mat(MatRing,Al)); end; ############################################################################### # # MatIsZero(A) # # Determines whether A is a zero matrix # This routine needs MatRing being set to the ring of definition of A. ############################################################################### MatIsZero:=function(A) local RingZero; RingZero:=MatElt(Mat(MatRing,[[0]]),1,1); if MatRows(A)=0 or MatCols(A)=0 then return true; fi; return(ForAll(Flat(MatToRowList(A)),u->u=RingZero)); end; ############################################################################### # # SpecSolution(A,b) # # Determines a special solution such that xA=b or returns false if it doesn't # exist. ############################################################################### SpecSolution:=function(A,b) local M,bmat,Mt,Msquare; M:=MatBasis(A); Msquare:=MatUnion(M,MatOrth(M)); bmat:=Mat(MatRing,[b]); if MatIsZero(bmat*MatTrans(MatOrth(M))) then Mt:=MatTrans(M); return(bmat*Mt*(M*Mt)^(-1)*A*Mt*(M*Mt)^(-1)); else return(false); fi; end; #--- matalg.g (end) ---# #--- powseries.mpl (start) ---# cat << "#--- powseries.mpl (end) ---#" > powseries.mpl # UNIX extraction command # Start of powseries.mpl # This file is a Maple V script that computes the power series mentioned in # [thesis, section 4.4] and [MI 99-14, section 3]. with(student): ecXYD:=gamma*Y^2*D-X^3-a2*X^2*D-a4*X*D^2-a6*D^3; XYDtoZW:={X=Z,Y=1,D=W}; ecZW:=simplify(subs(XYDtoZW,ecXYD)); Witer:=-(ecZW-gamma*W)/gamma; W0:=series(Z^3,Z,2); for it from 1 to 7 do W0:=series(subs(W=W0,Witer),Z,15); od; Gz:=simplify(subs(XYDtoZW,W=W0,[X,Y,D])); phihom:=[c[1,1]*X+c[1,2]*Y+c[1,3]*D,c[2,1]*X+c[2,2]*Y+c[2,3]*D]; Gxy:=[x,RootOf(subs(X=x,D=1,ecXYD),Y),1]; lineeq:=lambda*Gz+Gxy; lambda0:=-2*Gxy[2]; for e from 1 to 4 do Ceq:=collect(coeff(series(subs(equate([X,Y,D],lineeq),lambda=lambda0+C*Z^e,ecXYD), Z,e+1),Z,e),C); if degree(Ceq,C)<>1 then ERROR(`something wrong`); fi; Cval:=simplify(-coeff(Ceq,C,0)/coeff(Ceq,C,1)); lambda0:=lambda0+Cval*Z^e; od; lambda0:=series(lambda0+O(Z^e),Z,e); GzPlusGxy:=subs(Gxy[2]=y,simplify(map(T->series(T,Z,90), subs(lambda=lambda0,[lambda*Gz[1]+Gxy[1], -(lambda*Gz[2]+Gxy[2]),lambda*Gz[3]+Gxy[3]])))); Gx0:=[RootOf(subs(Y=0,D=1,ecXYD),X),0,1]; lineeq:=lambda*Gz+Gx0; lambda0:=(3*Gx0[1]^2+2*a2*Gx0[1]+a4)/gamma*Z; for e from 2 to 4 do Ceq:=collect(coeff(series(subs(equate([X,Y,D],lineeq),lambda=lambda0+C*Z^e,ecXYD), Z,e+2),Z,e+1),C); if degree(Ceq,C)<>1 then ERROR(`something wrong`); fi; Cval:=simplify(-coeff(Ceq,C,0)/coeff(Ceq,C,1)); lambda0:=lambda0+Cval*Z^e; od; lambda0:=series(lambda0+O(Z^e),Z,e); GzPlusGx0:=subs(Gx0[1]=x,simplify(map(T->series(T,Z,90), subs(lambda=lambda0,[lambda*Gz[1]+Gx0[1], -(lambda*Gz[2]+Gx0[2]),lambda*Gz[3]+Gx0[3]])))); Gphi:=[c[1,2]*c[2,3]-c[1,3]*c[2,2], c[1,3]*c[2,1]-c[1,1]*c[2,3], c[1,1]*c[2,2]-c[1,2]*c[2,1]]; phiGz:=simplify(series(series(subs(equate([X,Y,D],Gz),phihom[1]/phihom[2]),Z,30),Z,2)); phiGzPlusGxy:=simplify(series(series(subs(equate([X,Y,D],GzPlusGxy), phihom[1]/phihom[2]),Z,30),Z,2)); phiGzsym:=simplify(series(series(subs(equate([X,Y,D],Gz),c[1,2]=0,c[2,2]=0, phihom[1]/phihom[2]),Z,30),Z,3)); phiGzPlusGx0sym:=simplify(series(series(subs(equate([X,Y,D],GzPlusGx0),c[1,2]=0,c[2,2]=0, phihom[1]/phihom[2]),Z,30),Z,3)); tphiGz:=c[1,2]/c[2,2]+Gphi[3]/c[2,2]^2*Z; tphiGzPlusGxy:=(c[1,1]*x+c[1,2]*y+c[1,3])/(c[2,1]*x+c[2,2]*y+c[2,3])+ ( (3*x^2+2*a2*x+a4)*(x*Gphi[3]-Gphi[1])-2*gamma*y*(y*Gphi[3]-Gphi[2]) )/(c[2,1]*x+c[2,2]*y+c[2,3])^2/gamma*Z; tphiGzPlusGx0:=(c[1,1]*x+c[1,3])/(c[2,1]*x+c[2,3])+ (3*x^2+2*a2*x+a4)*(x*Gphi[3]-Gphi[1])/(gamma*(c[2,1]*x+c[2,3])^2)*Z; tphiGzsym:=c[1,1]/c[2,1]+(Gphi[2]/gamma/c[2,1]^2)*Z^2; tphiGzPlusGx0sym:=(c[1,1]*x+c[1,3])/(c[2,1]*x+c[2,3])- (3*x^2+2*a2*x+a4)*Gphi[2]/((c[2,1]*x+c[2,3])^2*gamma)*Z^2; sernot0:=(S)->evalb(convert(S,polynom)<>0); if sernot0(simplify(series(phiGz-tphiGz,Z))) or sernot0(simplify(series(phiGzPlusGxy-tphiGzPlusGxy,Z))) or sernot0(simplify(series(phiGzPlusGx0-subs(y=0,tphiGzPlusGxy),Z))) or sernot0(simplify(series(phiGzsym-tphiGzsym,Z))) or sernot0(simplify(series(phiGzPlusGx0sym-tphiGzPlusGx0sym,Z))) then ERROR(`forms are not correct`); fi; termZphiGz:=coeff(tphiGz,Z,1); termZphiGzGxy:=coeff(tphiGzPlusGxy,Z,1); termZ2phiGz:=coeff(tphiGzsym,Z,2); termZ2phiGzGx0:=coeff(tphiGzPlusGx0sym,Z,2); # save termZphiGz,termZphiGzGxy,termZ2phiGz,termZ2phiGzGx0,`termZphi.mpl`; # end of powseries.mpl #--- powseries.mpl (end) ---#