##################################### # f for MapleV, releases 4, 5 and 5.1 ##################################### daat:=` Nov. 21, 2000`:vern:=6: Abel:=proc() global iee; #Does E stabilize to everywhere good or multiplicative reduction over #some abelian extension of Q? Answer is yes (true) when the list #nabel:=[..[p,n]..] is empty; [p,n] denotes a prime p where E is "non-abel." #and n is a code -- see Menu(abel) for details. iee:=1:abel(args);iee:=0: if nabel=NULL then RETURN():fi; if nabel=[] then lprint('j'=j,`is abelian`); else lprint('j'=j,`is non-abelian at p for reason n using the notation [p,n]:`); print(op(nabel)); fi; lprint(`For an explanation of the terminology use the command Menu(abel);`); NULL end: abel:=proc() global nabel, j, a5, a6, p, vj, x; local e,e1,e2,e3,e4,e5,i,nj,pr; nabel:=NULL:pr:=[]: if nargs=0 then if bouv then aM();RETURN();fi; j:=jay:pr:=badp: elif nargs=1 and type(args,rational) then j:=args: else if nargs=5 then j:=[args]: elif nargs=1 and type(args,list) then j:=args: else print(qts):RETURN(): fi: a5:=a6: for i to 5 do e.i:=a.i:a.i:=j[i]:od: a6:=a5:j:=jay: for i to 5 do a.i:=e.i:od: a6:=a5: fi; if j=0 then nabel:=[[3,-1]]: elif j=1728 then nabel:=[[2,-2]]: else nj:=numer(j):nabel:=[]: if pr=[] then pr:=factorset(nj*(nj-1728*denom(j))):fi; for p in pr do if p=1 then next:fi; vj:=val(j): if vj<0 then next:fi: if vj=0 then if modp(p,4)=3 and modp(val(j-1728),2)=1 then nabel:=[op(nabel),[p,1]]: fi; next; fi; if modp(p,6)=5 then if modp(vj,3)>0 then nabel:=[op(nabel),[p,2]]:fi; elif p<5 then x:=(abel.p)(j): if x<>[] then nabel:=[op(nabel),x]:fi; fi; od; fi; nabel end: abel2:=proc(j); if modp(vj,3)>0 then RETURN([2,6]):fi; if vj>0 and vj<12 and (vj<>6 or modp(j/64,8)<>1) then RETURN([2,7]):fi; [] end: abel3:=proc(j) global y; local u; if modp(val(j-1728),2)=1 then RETURN([3,3]):fi; u:=[]: if vj=2 then if modp(j/9,3)=2 then u:=[3,4]:fi; elif vj=3 then y:=modp(j/27,27): if y<>4 and y<>10 and y<>13 and y<>22 then u:=[3,5]:fi: fi; u end: Absc:=proc(x); print(absc(x)) end: absc:=proc() local x; global abb, abc, esqrt; #find point [x,y] on E with given abscissa x:=sex_(args); abb:=a1*x+a3:abc:=x^3+a2*x^2+a4*x+a6: if not essqr(abb^2+4*abc) then esqrt:=simplify(sqrt(abb^2+4*abc),sqrt,symbolic); fi; [x,normal(sex_((-abb+esqrt)/2))] end: Afac:=proc() local q; if nargs=1 then print(args=afac(args));q:=`Now the`; else q:=`The`; fi; lprint(cat(q,` remember table of afac is`)); print(op(4,eval(afac)));NULL end: afac:=proc(n) local f, u; global afac_; options `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`, remember; #this is modelled on Maple's ifactors with the extra option remember #and calls afactor which is an "enhanced" version of ifactor if not (type(n,rational) or type(n,facint)) then ERROR(`afac takes a rational or a factored integer`) fi; if op(4,eval(afac))<>NULL and type(op(4,eval(afac))[-n],list) then f:=afac(-n);RETURN([-f[1],f[2]]); fi; afac_ := afactor(n); if member(afac_,[0, 1, -1]) then RETURN([afac_, []]) fi; if type(afac_,`*`) and op(1,afac_) = -1 then u := -1; afac_ := -afac_ else u := 1 fi; f:=[u, map(proc (x) [op(x[1]), x[2]] end,convert(afac_,multiset))]; afac(-n):=[-u,f[2]]; f end: afactor:=proc(n) local a, b, i, sol, r, t1; global afacnops, afacprimes, `ifactor/bottom`, p, t; option remember; #Adapted from Maple's ifactor if type(n,'integer') then if 0 < n then sol := 1; r := n elif n < 0 then sol := -1; r := -n else RETURN(0); fi; elif type(n,'fraction') then RETURN(procname(op(1,n),args[2 ..-1])/procname(op(2,n),args[2 .. -1])) else ERROR(`invalid arguments`) fi; unprotect(`ifactor`); t1 := igcd(r,720720); while t1 <> 1 do sol := sol*`ifactor/ifact235`(t1); r := iquo(r,t1); t1 := igcd(t1,r) od; t1 := igcd(r,1165429511437014421182564255394118062787055181073690352\ 484362724429288655197089578084537426127020962982242734844013923456967013254066\ 860138774356183909763696443067069059862065192038988478418219082943879223019472\ 668430243781122990303492985397077416788992156275405498099765357945192479722213\ 857227212000608128560171197659349424196117167227902201722912212770434834579779\ 297150544465360850511025966392537774945829901979588813143194547543583825739867\ 652770670299983480797335874643414610343435150187989728154996577619208807489382\ 147598111758441747971820027670613338896899878156111058787897307947215801259320\ 492843992867678444319238233981852191131219418884757870660310007540404320811366\ 733424902110215685045241323250289709233); while t1 <> 1 do sol := sol*`ifactor/ifact1st`(t1); r := iquo(r,t1); t1 := igcd(t1,r) od; if r<>1 and afacflag then t:=op(op(4,eval(afac))); if nops(t)>afacnops then afacnops:=nops(t); for i to afacnops do a:=op(2,t[i]); if a[1]=-1 then next fi; for b in a[2] do if b[1]<1700 then next else afacprimes:=afacprimes union {b[1]} fi od; od; fi; for p in afacprimes while r<>1 do if vam(r)>0 then sol:=sol*``(p)^v;r:=zz fi; od; fi; if r <> 1 then if nargs = 1 then `ifactor/bottom` := readlib('`ifactor/morrbril`'); t1 := `ifactor/ifact0th`(r) else `ifactor/bottom` := readlib(`ifactor/`.(args[2])); t1 := `ifactor/ifact0th`(r,args[3 .. nargs]) fi; if t1 <> FAIL then sol:=sol*t1 else sol:=FAIL fi fi; protect(`ifactor`); sol end: AgM:=proc(x,y) local a; a:=agM(args); if a<>NULL then lprint(`The arithmetic-geometric mean of `,x,` and `,y); lprint(` to `,Digits,` digits is`); print(a); fi; NULL end: agM:=proc(x1,y1) local x, y, z; #Calculate Gauss's arithmetic-geometric mean of x1 and x2. Gauss actually #denoted this mean by agM. x:=evalf(x1);y:=evalf(y1); if not type(x,numeric) or not type(y,numeric) or x<=0 or y<=0 then lprint(`args of agM must be positive reals`,qtr);RETURN(); fi; x:=sort([x,y]):y:=x[1]:x:=x[2]: while (x-y)>10.^(-Digits+2)*x do z:=(x+y)/2:y:=sqrt(x*y):x:=z: od; (x+y)/2 end: aM:=proc(); print(`A non-sing. E/Q must be named using the apecs command ein`) end: aN:=proc(); print(`The Weierstrass coeff.'s a1,..,a6 must be rational`,qtr) end: ap0:=proc () local bqe, q; global iee, qqR, sB, tB; q:={op(QQ)} union ({op(RR)} minus uR); if nops(q)>0 then if nargs<2 then do_tor(args);fi; bqe:=iee;iee:=0;map(ap1_,q);iee:=bqe; if bouv then RETURN();fi; fi; if qqR = 1 then if 1 < sB and 0 < MM then lprint(` We must redo seek(sB) to replace members of RR by points \ with`); lprint(`smaller canon ht when possible, where sB`=sB); tB := sB; sB := 1;seek(tB); fi; qqR := -1; fi; end: ap1_:=proc() local api,h,ha,hflag,j,k,l,m,n,rR,Rz,Sz,z; global apa, i, iee, MM, my, ouP, PP, QQ, r4, rc, reg, RR, uR, x3, y3, yR; #(internal) Adjoin new point Rz and its siblings +-Rz+T, T torsion, to the #set uR; put it in RR if appropriate. When K_=1 #loop through yR to see if RR can be improved; carry any changes over to #isogenous curves if feasible z:=args; if z=[] or z=[`?`] or z=O then RETURN();fi; do_tor(1);api:=iee; if member(z,ouP) then fprint(z,` has finite order.`);RETURN();fi; if NN=0 then n:=max(10,nops(ouP)+1);ford(z,n);uR:=uR union {Qq}; if ord<=n then eprint(z,` is a torsion point of order`,ord); iee:=0;subg(op({PP} union {Qq}));iee:=api; PP:=op(GG);print_PP(); else if member(z,QQ) then apa:=5; elif member(neg(z),QQ) then apa:=6; else QQ:=[op(QQ),z];apa:=7; fi; fi; RETURN(); fi; if check_pt(z) then RETURN();fi; rR:=RR:ha:=0;Sz:=z; while nops(half(z))>0 do z:=q[1];ha:=ha+1;od; Rz:=z; if member(Rz,uR) then if member(Rz,RR) then apa:=1:RETURN():fi: if nops(yR)=0 or member(Rz,yR) or member(neg(Rz),yR) then fprint(Rz,` is redundant for RR of `,cur);RETURN(); fi; fi; avR(Rz); if bouu then if member(Rz,uR) then if vR intersect {op(yR)}={} then apa:=2:iee:=api:RETURN():fi: apa:=3: else apa:=4: fi: fi; uR:=uR union vR:yR:=[op(yR),Rz]:hflag:=1: while hflag=1 do hflag:=0:my:=nops(yR): for h from my by -1 to 1 while hflag=0 do Rz:=yR[h]:yR:=eau(yR,h):iee:=0:rR:=RR;j:=[op(RR),Rz];i:=MM+1; if bouu then j:=sort(j,proc(a,b):if ht_(a)0 then if ha>0 then fprint(Sz,` replaced by new point`,Rz); fprint(2^ha,Rz=Sz); fi; RR:=lsp;MM:=nops(RR):reg:=adet:hflag:=1: fprint(`New point `,Rz); else if lsp[i]<>Rz then RR:=eau(lsp,i); elif iee=1 and not member(Rz,RR) then lprint(Rz,` is redundant for RR of `,cur); fi; if bouu then if den>1 and not member(lsp[i],yR) then yR:=[op(yR),lsp[i]]: if lsp[i]<>Rz then hflag:=1:fi; fi; reg:=lin_(RR); fi; fi; nowRR(rR,cur); od; od; MM:=nops(RR):iee:=api: if bouu then if r1>=0 then if MM=r1 then rc:=7:r4:=r1:fi; if r4=-1 then r4:=r1:fi; fi; if qR>0 and RR<>rR and iflag=0 then apq();fi; if RR<>rR then upis();fi; fi; MM end: App:=proc() global iee; iee:=1:app(args);iee:=0: if apa=1 then print(apz,`is already in RR.`): elif apa=4 and member(apz,yR) then print(apz,`has been added to the list yR.`): lprint(`Now yR = `,yR): elif apa=5 then print(apz,`was already in the list QQ.`); elif apa=6 then print(`Th negative of`,apz,`was already in the list QQ.`); elif apa=7 then print(apz,`has been added to the list QQ.`); elif apa>1 and not member(apz,RR) then print(apz,`is redundant.`): fi: NULL end: app:=proc() global apa, apz; #Adjoin new point to uR as in ap1_. (This is the user-friendly version of the #internal procedure ap1_.) apa:=0: if sT=4 then print(qtt);RETURN();fi; apz:=sex_(args); if nargs=2 then apz:=[apz];fi; if not type(apz,list) then print(`wrong args in app`);RETURN();fi; if nops(apz)=1 then apz:=absc(apz[1]);fi; if not on(apz) then print(apz,noc);RETURN(); elif check_pt(apz) then RETURN(); fi; ap1_(apz);NULL end: apq:=proc() global iflag, StorE, storE, X, y, Y, z, zI; local a,b,cpq,d7,dpq,Iq,Jq,l,m,Mpq,npq,rr,vpq; #(internal) Transfer points of RR (see # for ap1_) to isogenous curves (when #not too much work). if nops(nisog)=0 or MM=0 then RETURN();fi; npq:=[]:dpq:=[]:vpq:=[]:iflag:=1: for Iq to nops(nisog) do if not has(visog[Iq],X) and abs(disog[Iq])<8 then tri(nisog[Iq]):fi; if has(visog[Iq],X) then y:=nisog[Iq]: npq:=[op(npq),y]:dpq:=[op(dpq),disog[Iq]]:vpq:=[op(vpq),visog[Iq]]: fi; od; if nops(npq)=0 then iflag:=0:RETURN();fi; rr:=RR:Mpq:=MM:cpq:=ncur: for Iq to nops(npq) do StorE:=storE;storE:=true:ein(npq[Iq]);storE(StorE):zI:=vpq[Iq]:d7:=dpq[Iq]: for Jq to Mpq do X:=rr[Jq][1]:Y:=rr[Jq][2]:z:=[zI[2],zI[3]]:X:='X':Y:='Y':ap1_(z): od; od; StorE:=storE;storE:=true:ein(cpq);iflag:=0:storE(StorE):NULL end: avR:=proc() local b, l, m, s, z, z1; global Rz, vR; Rz:=args;s:=Rz,seq(eadd0(Rz,z1),z1=ouP); b:=type([op(yo),op(Rz)],list(rational)); if b then vR:={s,seq(neg(z),z=[s])}; elif K_=m_ then vR:=map(qf_conv,{s,seq(neg(z),z=[s])}); else vR:=sex_({s,seq(neg(z),z=[s])}); fi; if b then m:=denom(Rz[1]); for z in vR do l:=denom(z[1]); if labs(z[1]) or (Rz[1]=z[1] and (abs(z[2])\ 0 then badmp:=[op(badmp),p]:badme:=[op(badme),ep]: if p<5 then aa:=emoe(p): else aa:=mods(power(-c6,1/2*p-1/2),p): fi; SFEm:=-SFEm*aa: else aa:=0: fi; badt:=[op(badt),aa]: od end: Bas:=proc() global iee; iee:=1:bas(args);iee:=0:NULL end: bas:=proc() local a,ape,b,B,i,MM1,nkur,w,x,y,z; global BASL, b2, baSM, basL, feJ, htM, iee, j, k, nSmL, niM, niv, pcn , pp, q, q1, qR, qqR, RRsm, shaf, tim, tom, top, ub_index, uS, v, X, Y, z1; #see Menu(bas) for description of this proc. ape:=iee: if nargs=0 then BASL:=2000:else BASL:=args:fi; if bouv then aM();RETURN();else do_tor();fi; if tim>0 then eprint(`clock started`);tom:=time():fi; isog(); for a in nisog do b:=eT[a]; if nops(b[6])=0 then next;fi; member(ncur,b[11],'i');B:=b[13][i]; if not has(B,X) then next;fi; for z in b[6] do X:=z[1];Y:=z[2];j:=eval(B);X:='X';Y:='Y'; if j[1]<>0 then ap1_([j[2],j[3]]);fi; od; od; if r4<0 or rc<0 or (rc<4 and qR<2) then if NNeven and dill>3 then a:=trunc(dill*1.5);else a:=NULL;fi; fprint(`We try to get rank without conj.:`);niv:=1:rkNC(a);iee:=ape; if rc>=0 then if rc<4 then fprint(`We'll make do with the conjectural rank r4 = `,r4); else fprint(`So the rank and "quality" of the rank are now r4,rc = `, r4,rc); fi; else fprint(`We must assume standard conj.'s to determine the rank.`); if NNeven then a:=dill+30;else a:=NULL;fi; if Nc<10^6 then rk(a); elif Nc<10^20 then rk1(a); else rk0(a); fi; fi; if rc<0 then fprint(`Rank not yet determined; apecs suggests discretionary use of Mest`); fprint(`& Seek. (Note that Seek1 can be much quicker when there are`); fprint(`sufficiently many integral points.)`); RETURN(qR); fi; fi; niv:=0: if tim<0 then tim:=0:RETURN(qR);fi; if r4=0 then pr4();qR:=5:qqR:=-1:shaf_tr();RETURN(5): elif qR>4 then eprint(`a basis has already been found:-`);bas1();eprint(op(RR)); fprint(`The regulator = `);eprint(reg);shaf_tr(); RETURN(qR); fi; ouu(200);iee:=ape;isog();MM1:=MM:iee:=ape:hdif_Sik(); if NNeven then two();else q1:=0:a:=50;fi; for feJ to 10 while (MM0) do if NNeven then if nops(dol[1])+nops(dol[2])>0 then rkNC(dill+30);fi; a:=trunc(evalf(20*sqrt(sB+2))); fi; ouu(sB+a);iee:=ape: if tim<0 then tim:=0:RETURN(qR);fi; od; if MM=r4 then pcn:=1:a:=ht_(RR[1]):pcn:=0: if MM=1 then a:=a/4:fi; baSM:=max(1,trunc(evalf(exp((a-lowB)*2)))); else if MM>MM1 then ouu(max(trunc(sB*1.2),sB+100));iee:=ape:fi: if baSM='baSM' then baSM:=BASL:fi; fi; ouu(min(baSM,BASL));iee:=ape; MM1:=max((BASL-sB)/4+1,100); while MM0 then i:=max(20,trunc(BASL/50)); eprint(`Resorting to seekb(i) where i`=i);iee:=0;seekb(i);iee:=ape; fi; if MM`,sB); else lprint(`To get the first successive minimum for h^, which is a \ minimum requirement`); lprint(`for this procedure, call Bas again with args`=baSM); fi; fi; if MM0 then qR:=5/2:fprint(`We have found the right number of independent points`); fprint(`but because yR isn't empty RR is not a basis.`); else qR:=3:shaf_tr(); fi; RETURN(qR); fi; MM1:=0; while MM1=0 do MM1:=1:B:=array(1..MM,1..MM):iee:=0;pcn:=1; for i to MM do h.i:=ht_(RR[i]):B[i,i]:=2*h.i:od; for i from 2 to MM while MM1=1 do for j to i-1 do x:=eadd(RR[i],RR[j]):y:=ht_(x); if ysB then RRsm:=b-1;a:=a*h.RRsm^(r4-RRsm);nSmL:=trunc(i+1);break;fi; a:=a*h.RRsm; od; if r4<9 then i:=MINK[r4];else i:=(4/Pi)^r4*GAMMA(r4/2+1)^2;fi; ub_index:=trunc(evalf(sqrt(i*det(B)/a/2^r4)+.01)); if ub_index=1 then fprint(`By a result of Siksek, RR is certified as a Mordell-Weil basis:`); bas1();print(` RR`=RR);qR:=5;shaf_tr();RETURN(qR); fi; if RRsm=r4 then if nops(yR)>0 then print(`!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!`);print(); print(`This appears to be an example where the succesive minima`); print(`don't form a Mordell-Weil basis!`); print(\ `please send details (at least the 5 Weierstrass coef.'s a1,..,a 6)`); print(`to connell@math.mcgill.ca`);print(qtr);qR:=7:RETURN(7); fi; if r4<5 then if iee=1 then bas1(); if r4=1 then lprint(`The point`,RR[1],` with height`);print(h1); lprint(`is a Mordell-Weil basis.`); lprint(`The regulator is twice the height:`); print(` reg`=reg); else lprint(`since rank<5, Minkowski's theorem applies:`); lprint(`the succesive minima form a Mordell-Weil basis;`); lprint(`in order they are`); print(RR);lprint(`with canonical heights, respectively,`); print(h.(1..MM));lprint(`and with regulator`);print(reg); fi; fi; qR:=5;shaf_tr();RETURN(5); fi; eprint(`upper bound for index of +tor in E(Q) = ub_index`=ub_index); qR:=4; fprint(`RR consists of successive minima, so is probably a basis. But \ since`); fprint(`rank r4>4, Minkowski's result doesn't apply.`); else qR:=3: if RRsm=1 then fprint(`RR[1] is the first of the successive minima.`); else fprint(`The first RRsm`=RRsm,` of the`,r4,` members of RR are \ successive minima.`); fi; fprint(`The next succ. min. will be found by Seek(`,nSmL,`)`); fi; fprint(`We now apply Siksek's filter to the primes`<=ub_index); i:=0; for z in PP do if z[3]=2 then i:=i+1;fi;od; if i=3 or sik_siev(2) then if NNeven then if Rflag = 0 then two() fi; a:=NULL;b:=1; for i in [PP] do if 1 < igcd(i[3],15) then next fi; if i[3]=2 then a:=a,i fi; if b0 then fprint(`Halving linear comb.'s of known points uncovered the poin\ t `,q[1]); ap1_(q[1]):ape:=1:MM1:=0;break: fi; od; od; if MM1=0 then eprint(`So we restart Siksek's sieve.`);next;fi; fi; pp:=3:fprint(`index not divisible by 2`);nkur:=0; while pp<=ub_index do if sik_siev(pp) then fprint(`index may be divisible by`,pp);nkur:=1: else fprint(`index not divisible by`,pp); fi; pp:=nextprime(pp); od; if nkur=0 then qR:=5:fprint(`Siksek's sieve certifies that RR =`); eprint(RR);fprint(`is a Mordell-Weil basis.`);bas1(); elif sB 1, so we try harder.`); if RRsm1. You can try Bas(args); with args >`,MM1); fi; od; if qR>=3 then shaf_tr(); if qR=3 and iee=1 then lprint(`The following independent point(s) possibly form a basis:`); print(op(RR)); fi; else eprint(`Apply Bas(xx); or Seek(xx); for some`); eprint(`xx> the current search limit sB = `,sB); fi; tim:=0:qR end: bas1:=proc(); if rc<4 and iee=1 then lprint(`(assuming the rank`=r4,` --- we used standard conj.'s to get t\ his value)`); fi end: bc:=proc() global a6, NN, ouP, PP, sT, tor_able, uR, yo; local i; #(internal) Note that b2,..,b8,c4,c6,DD,jay and various other things such as #the Fricke polynomials fr5,.. have already been defined in terms of the #weierstrass coefficients ai in the file apecs, and so they acquire #specific values when the values of the ai are specified. #trick: a dummy 6-th arg forces bc to consider simplifying args[1..5] #and to set sT and NN if nargs=5 then for i to 5 do a.i:=args[i];od; a6:=a5:tor_able:=true;yo:=[a.(1..5)] else for i to 5 do a.i:=sex_(args[i]):od; a6:=a5:yo:=[a.(1..5)]; if sex_(DD)=0 then sT:=4;NN:=0;uR:={};ouP:=[];PP:=NULL;tor_able:=false; else sT:=-1;NN:=-1; tor_able:=evalb(type(yo,list(rational)) or K_=T or K_=m_); fi; fi; yo end: bsGs:=proc() local N,r,ii,m,g1,n; #Shank's baby step Giant step method to determine the order of a point #g on E mod p (or on its unique twist mod p --- j<>0 or 1728 mod p). #Called by newp which in turn is called by emoe, which is the general #(internal) proc. to calc. the trace of Frob. of E mod p. n:=nnn;N:=[n]; if n=[] then r:=p-e;else r:=0;fi; for ii to es-1 do n:=eedp(n,g,tb2,tb4); if n=[] then if r=0 then r:=p-e+ii; else RETURN(p-e+ii-r); fi; fi; N:=[op(N),n]; od; g1:=malp(es,[g[1],-g[2]],tb2,tb4); if g1=[] then RETURN(es);fi; n:=g1; for ii to 2*es do if member(n,N,'j') then if r=0 then r:=p-e+j-1+ii*es; else RETURN(p-e+j-1+ii*es-r); fi; fi; n:=eedp(n,g1,tb2,tb4); od; [r] end: buc:=proc() global yo; #(internal) used by con0, ein, gcub and qoo (the latter two in biraF). if iee=1 then lprint(`The coordinates U,V of the original equation are related to the`); lprint(`coordinates X,Y of the Weierstrass equation by`); print(U=args[1]/args[2],V=args[3]/args[4]); print(X=args[5]/args[6],Y=args[7]/args[8]); lprint(`To transfer points between the original curve and the Wei\ erstrass form`); if nops(bic)>1 then lprint(`use the commands Trcw(u,v,n) and Trwc(x,y,n).`); lprint(`where the code number n = `,nops(bic)); else lprint(`use the commands Trcw(u,v) and Trwc(x,y).`); fi; lprint(`These commands remain available during this apecs session.`); # elif iee=0 then print(cus,cur=yo); fi end: check_pt:=proc(z) if K_=1 and not type(z,list(rational)) then eprint(z,` not defined over Q --- point discarded.`); RETURN(true); elif (K_=T and not type(z,list(ratpoly(rational,T)))) or (K_=m_ and not type(z,list(ratpoly(rational,w_)))) then eprint(z,` not defined over field with label K_`=K_,` --- point discarded`); RETURN(true); fi; false end: chco:=proc(x,y,z) global a6p, a4p, a3p, a2p, a1p, b8p, b6p, b4p, b2p; ; #(internal) used by tate, tate2 and tate3 a6p:=a6p+x*(x*(x+a2p)+a4p)-z*(z+x*a1p+a3p): a4p:=a4p+x*(3*x+2*a2p)-(z+x*y)*a1p-y*(a3p+2*z): a3p:=a3p+x*a1p+2*z: a2p:=a2p-y*(a1p+y)+3*x: a1p:=a1p+2*y: b8p:=b8p+x*(x*(x*(3*x+b2p)+3*b4p)+3*b6p): b6p:=b6p+x*(x*(4*x+b2p)+2*b4p): b4p:=b4p+x*(6*x+b2p): b2p:=b2p+12*x end: Cld:=proc() global iee; iee:=1:cld();iee:=0;NULL end: cld:=proc(); if type(yo,list(rational)) then cld_(); if uV>1 then if iee>0 then lprint(`Denom.'s of the Weier. coeff.'s cleared by uV`=uV); Ell(a.(1..5),K_);print(cus,yo); else ell(a.(1..5),K_); fi; else fprint(`All Weier. coeff.'s have denom. = 1 -- curve stays the same.`); fi; else lprint(`E must be defined /Q for the procedure cld.`); print(qts); fi; uV end: cld_:=proc() global a5, uV, a6; local i,j,k,p; #(internal) for weierstrass coefficients in Q, Q(sqrt(m_)) or Q(T), #find lcm of denom.'s (denoted uV) and transform coord.'s to get #integral coeff.'s a5:=a6:uV:=1: for p in factorset(mul(denom(a.i),i=1..5)) do for j to 5 do while modp(denom(a.j),p)=0 do for k to 5 do a.k:=p^k*a.k:od; a5:=a5*p:uV:=uV*p: od; od; od; a6:=a5 end: clk:=proc() global tim, tom; ; #(internal) used by Clock. if tim+tom0 then print(`Clock is operating,`);x:=tim+tom-time(): if x>0 then print(`remaining time`=x); else print(`but time is up. Call Clock(x) again if wanted.`);tim:=0: fi; else print(`Clock not operating.`); fi; elif nargs=1 then x:=args: if type(x,numeric) and x>0 then tim:=x:lprint(`The following apecs procedures will be halted after run\ ning `,x,` seconds :`); print(`Gcub,gcub,Ind=Bas,ind=bas,Quar,quar,Rk,rk,Sfe,sfe,Sha,sha`); lprint(`More precisely, the time limit will be imposed if and when the\ se procedures`); lprint(`become involved in evaluating infinite series or searching f\ or rational points.`); lprint(`N.B. Clock is stopped after each use - it must be re-called if wa\ nted again.`); else lprint(qts); fi; else lprint(qts); fi; end: clock:=proc(x) global tim; ; #see first few lines of Clock. if type(x,numeric) and x>0 then tim:=x: else lprint(qts); fi end: Comb:=proc(); comb(args); if xc=NULL then print(`suM = O`); elif xc<>'xc' then print('suM'=suM); fi end: comb:=proc() global xc, yc, pcn, suM; local i,n,z; #calculate n1z1+...+nkzk for integers ni (pos. or neg.) and points zi=[xi,yi] #on the elliptic curve; result denoted suM=[xc,yc]. This proc. even accepts #non-singular points zi satisfying a singular Weierstrass equation. if nargs=0 or modp(nargs,2)=1 then print(qts);xc:='xc':RETURN(); elif sT=4 and has(args,zs) then print(`The singular point`=zs,` is not allowed as an arg.`,qts); xc:='xc':RETURN(); fi; xc:=NULL:yc:=NULL:pcn:=1: for i to nargs/2 do n:=args[2*i-1]:z:=args[2*i]: if n=0 or z=[] then next:fi: if not on(z) and not hastype(z,float) then print(z,noc);xc:='xc':pcn:=0:RETURN(); fi; mult(n,z); if xc=NULL then xc:=xm:yc:=ym: else if xm=NULL then next;fi; eadd(xc,yc,xm,ym);xc:=x3:yc:=y3: fi; od; pcn:=0: suM:=[xc,yc] end: con0:=proc() local Weq; global ape, bec, bic, boc, i, iee, p_at_inf_c, p_at_inf_w, sBb, stoRE, storE, U_, V_, X_, Y_; ; #(internal) calculate minimal weierstrass form and store birational #transformations ape:=iee:iee:=0:stoRE:=storE:storE:=true: ein(a1p,a2p,a3p,a4p,a6p);storE(stoRE):iee:=ape: if DD=0 then print(`curve has genus 0`,qtr); else eprint(cus,cur=yo); X_:=simplify(X_,{czb});Y_:=simplify(Y_,{czb}); U_:=normal(U_):V_:=normal(V_):X_:=normal(X_):Y_:=normal(Y_): Weq:={Y^2=-a1*X*Y-a3*Y+X^3+a2*X^2+a4*X+a6}; boc:=[simplify(numer(U_),Weq),denom(U_),simplify(numer(V_),Weq),\ denom(V_),numer(X_),denom(X_),numer(Y_),denom(Y_)]: bic:=[op(bic),boc]:buc(op(boc));bec:=[op(bec),czb]:sBb:=0; p_at_inf_c:=[op(p_at_inf_c),[]];p_at_inf_w:=[op(p_at_inf_w),[]]; fi; U_:=X:V_:=Y:X_:=U:Y_:=V: end: con9:=proc(); fprint(`We now invoke the apecs command ein(`,a1,a2,a3,a4,a6,`)`); fprint(`to determine the minimal Weierstrass form.`); cld_();weier() end: diop:=proc() local hy,Lv,n; global Ah,b,Bh,bh,bbh,co_list,hu,i,iop,Lc,Ld,Lt,Lu,r3,R3,r4,R4; #(internal) called twice by eiop Lt:=args[1]:Lu:=args[2]:Lc:=Lt union Lu:Lv:=Lu; Ah:=args[3];Bh:=args[4];hu:=args[5];co_list:=[]; if iee>0 then lprint(`Pondering the homogeneous spaces`); print(cat(`d1*u^4 + `,Ah,`*u^2*v^2 + d2*v^4 = w^2`)); lprint(`where d1*d2`=Bh,` and d1 is from`); print(Lu); fi; while nops(Lv)>0 do bh:=Lv[1];bbh:=Bh/bh;hy:=mss();diyp(); Ld:=map(proc(e) e*bh/igcd(e,bh)^2 end,Lt):Lv:=Lv minus Ld: if hy>0 then Lu:=Lu minus Ld: if hy<4 then R4:=R4/nops(Lc):Lc:=Lc minus Ld: R4:=R4*nops(Lc):r4:=sam(R4); else R3:=R3/nops(Lt):Lt:= Lt union Ld:R3:=R3*nops(Lt):r3:=tam(R3): fi; else co_list:=[op(co_list),co_lust];n:=nops([co_lust])-1;i:=nops(h1_); if i=1 then fprint(`There is one associated quartic,`); if n=1 then fprint(`and this quartic is loc. sol.`); else fprint(`but this quartic is not loc. sol.`); fi; else fprint(`There are`,i,` associated quartics.`); i:=` loc. sol.`; if n=0 then i:=cat(`none is`,i); elif n=1 then i:=cat(`one is`,i); else i:=cat(n,` are`,i); fi; fprint(`Of these`,i); fi; fi; od; Lc:=Lt union Lu end: diyp:=proc() local l,z;global rflag,uflag; if ncur<>cn then l:=[]; for z in RR do l:=[op(l),tri(z,cn)];od; ein(cn); for z in l do ap1_(z);od; fi; rflag:=3;uflag:=2 end: Div:=proc(n) global i; if not type(n,posint) then print(qts);RETURN();fi; div(n);i:=al.n:print(cat(`al`,n)=i); if modp(n,2)=1 then lprint(cat(`= the usual division polynomial psi_`,n)); else lprint(cat(`= the usual division polynomial psi_`,n)/[2*Y+a1*X+a3]); fi; NULL end: div:=proc(n) global i, aln;local r,m; #calc. n-th div. poly. by recurrence if not type(n,posint) then print(qts);RETURN();fi; if n<3 then RETURN(1): elif n<5 then RETURN(al.n): fi; for i from trunc(aln/2)+3 to trunc(n/2)+2 do if i<>aln then div0(i):fi; od; aln:=n:div0(n) end: div0:=proc(i) global r, m; : r:=modp(i,4);m:=trunc(i/4): if r=1 then al.(4*m+1):=collect(sex_(al0*(al.(2*m+2))*(al.(2*m))^3- (al.(2*m-1))*(al.(2*m+1))^3),X): elif r=2 then al.(4*m+2):=collect(sex_((al.(2*m+1))*((al.(2*m+3))*(al.(2*m))^2- (al.(2*m-1))*(al.(2*m+2))^2)),X): elif r=3 then al.(4*m+3):=collect(sex_((al.(2*m+3))*(al.(2*m+1))^3- al0*(al.(2*m))*(al.(2*m+2))^3),X): else m:=m-1:al.(4*m+4):=collect(sex_((al.(2*m+2))*((al.(2*m+4))*(al.(2*\ m+1))^2-(al.(2*m))*(al.(2*m+3))^2)),X): fi end: do_tor:=proc() local e; global iee; if NN=-1 then if tor_able then fprint(`Must do tor first.`);fi; e:=iee;iee:=0;tor();iee:=e; if NN>0 and nargs=1 then eprint(`Order of torsion subgroup`=NN);fi; fi end: Eadd:=proc() global Edd; Edd:=0:eadd(args); if Edd=0 then if x3<>NULL then lprint([x1,y1],`+`,[x2,y2],`=`,[x3,y3]); else lprint([x1,y1],`+`,[x2,y2],`= O`); fi; fi end: eadd:=proc() local m; global Edd, x3, y3; #add two points on elliptic curve orgs(args); if pcn=0 then if not on([x1,y1]) and not hastype([x1,y1],float) then Edd:=1:print([x1,y1],noc);RETURN(); elif sT=4 and ([x1,y1]=zs or [x2,y2]=zs) then print(`The singular point`=zs,` is not allowed as an arg.`,qtr); Edd:=1:RETURN(); fi; if not on([x2,y2]) and not hastype([x2,y2],float) then Edd:=1:print([x2,y2],noc);RETURN(); fi; fi; if x1=NULL then x3:=x2:y3:=y2:RETURN([x2,y2]): elif x2=NULL then x3:=x1:y3:=y1:RETURN([x1,y1]): fi; if x1=x2 then m:=sex_(y1+y2+a1*x1+a3): if m=0 then x3:=NULL:y3:=NULL:RETURN([]): else m:=rationalize(simplify((3*x1^2+2*a2*x1+a4-a1*y1)/m)): fi; else m:=rationalize(simplify((y2-y1)/(x2-x1))): fi; x3:=sex_(m^2+a1*m-a2-x1-x2): y3:=sex_((x1-x3)*m-y1-a1*x3-a3):[x3,y3] end: eadd0:=proc(z1,z2) local m;global x3,y3; # internal: z1+z2 \neq O m:=(z1[2]-z2[2])/(z1[1]-z2[1]):x3:=m^2+a1*m-a2-z1[1]-z2[1]: y3:=(z2[1]-x3)*m-z2[2]-a1*x3-a3: [x3,y3] end: eao:=proc() global eae, eai, eaen; local b,c,d,n,p; #(internal) prolong table of Fourier coefficients up to next prime for n to nops(nisog) do if member(et1[nisog[n]],ft1,'b') and ft[b][4]>max(eai,ebi) then eae:=copy(ft[b][5]):eai:=ft[b][4]:RETURN(); fi; od; p:=nxtprm(eai): if not type(eae[p],integer) then eae[p]:=emoe(p):fi: #mest can insert eae[p] for prime p without calculating intermediate eae[n]. for n from eai+1 to p-1 do d:=ifactors(n)[2];c:=d[1][2];b:=d[1][1]; if nops(d)=1 then eaen:=eae[b]*eae[b^(c-1)]: if modp(Nc,b)>0 then eaen:=eaen-b*eae[b^(c-2)]:fi; else b:=b^c:eaen:=eae[b]*eae[n/b]: fi; eae[n]:=eaen: od; eai:=p end: eau:=proc(x,i); #(internal) delete i-th coordinate of list [op(x[1..(i-1)]),op(x[(i+1)..nops(x)])] end: ebo:=proc(pa) global am, pm, eae, ebe, ebu, ebi; local a,b,c,p,m; #internal - used by mest to extend table ebe am:=a^m:pm:=p^m:p:=1: while p1 then if n=1 then ntab:=[op(ntab),op(afac(zz)[2])]: else for p in badp do if p=2 or p=3 then next:fi; vam(zz); if v>0 then ntab:=[op(ntab),[p,v]]:fi; od; fi; fi; nP:=nops(ntab) end: efad:=proc(x,n) global um, ue;local z; #(internal) find biggest n-th power divisor efac(x);um:=1: for z to nP do ue[z]:=trunc(ntab[z][2]/n);um:=um*(ntab[z][1]^ue[z]): od; um:=expand(um) end: Ein0:=proc() global iee; iee:=1:ein0(args);iee:=0:NULL end: ein0:=proc() global z, U_, V_, X_, Y_, i, a6, bac, a3, a4, czb,a1p,a2p,a3p,a4p,a6p; ; #like ein except allow coefficient a0 of X^3 if nargs=1 then z:=args: else z:=[args]: fi; if nops(z)<>6 then print(`ein0 expects 6 args`,qtr):RETURN();fi; if z[1]=0 then U_:=X:V_:=Y:X_:=U:Y_:=V:print(`curve has genus 0`,qtr);RETURN(); fi; for i to 6 do a.(i-1):=z[i]:od: a6:=a5: bac:=a0*U^3+a2*U^2+a4*U+a6-V^2-a1*U*V-a3*V:a3:=a3*a0:a4:=a4*a0:a6:=a6*a0^2: if DD<>0 then if U_=X then czb:=bac:eprint(cus,czb=0);fi; con9();unk((Wu^2*XR+Wr)/(a0*uV^2),(Wu^3*YR+Wu^2*Ws*XR+Wt)/(a0*uV^3)); xnk((a0*UR*uV^2-Wr)/Wu^2,(a0*VR*uV^3-a0*Ws*UR*uV^2+Ws*Wr-Wt)/Wu^3); else a1p:=a1:a2p:=a2:a3p:=a3:a4p:=a4:a6p:=a6: fi; con0();NULL end: Ein:=proc() global iee; iee:=1:ein(args):iee:=0:NULL end: ein:=proc() global a, a1, a2, a3, a4, a6, alf, allP, bac, bec, bic, boc, bouu, bouv, cP, cplist, crem_dat, cur, curi, dij, dill, disog, Dn, Do, dolt, Dt, e, e4, e6, eae, eai,ebe, ebi, eT, et1, et2, etj, ewq, eyr, ft1, ft2, i, j, jm, jn, K_, Kod, kr, lds, lowB, misc, MM, Mp, Nc, Ncl, ncur, Ner, nisog, NN, Ny, oK_, ouP, p, pcn, PP, prank, p_at_inf_c, p_at_inf_w, Py, qqR, qR, q_1, r1, r4, raC, raF, rc, reg, Rflag, RKI, RKJ, RR, QQ, sB, sBb, se4, se6, SFE, shaf, stac, trou, tuQ, tW, u, U, uppB, uR, uuV, V, visog, Wr, Ws, Wt, Wu, yo, yR, yW, z, zst; local d,jk,S,t,x,y,ya,yy; #input (rational) weierstrass coefficients a1,..,a6 of elliptic curve oK_:=K_;K_:=1; if oin(args)>0 then go(psn);RETURN();fi; if K_<>1 or DD=0 then ell1(yo,K_);RETURN();fi; if not type(yo,list(rational)) then eprint(`ein expects rational args --- field label changed to K_=0.`); K_:=0;ell1(yo,0);RETURN(); fi; ewq:=0:Mp:=-1:boc:=[]:Rflag:=0:raC:=[]:uppB:=0:Ner:=NULL;q_1:=0; stac:=0:RKI:=[];RKJ:={};dolt:=0;prank:='prank';bouu:=true:bouv:=false: if member(yo,ft2,'zst') then stec();RETURN(cus,cur); #When K_=1, the ft2 entry is [a1,..,a6] rather than [a1,..,a6,K_] elif member(yo,et2,'ncur') then eon(); else cur:=new: fi; x:=yo: if cur=new then bac:=U^3+a2*U^2+a4*U+a6-V^2-a1*U*V-a3*V:cld_(); Do:=DD:weier();uuV:=u*uV:bc(a1p,a2p,a3p,a4p,a6p): if uuV>1 then boc:=[(Wu^2*X+Wr),uV^2,(Wu^3*Y+Wu^2*Ws*X+Wt),uV^3,(U*uV^2-Wr),Wu^2,(\ V*uV^3-Ws*U*uV^2+Ws*Wr-Wt),Wu^3]:curi:={};uR:={}; PP:=op(map(ein_tr,[PP]));QQ:=map(ein_tr,QQ); RR:=map(ein_tr,RR); NN_pair();oin_();zst:='zst'; if member(yo,ft2,zst) then stec();RETURN(cus,cur);fi; if member(yo,et2,'ncur') then eon();fi; fi; efac(DD);Dt:=ntab:Dn:=nP: else Wu:=1:Ws:=0:Wr:=0:Wt:=0:u:=1: fi; ill();eai:=1:eae:='eae':eae[1]:=1:dij:=2: eyr:=0:ebi:=1:ebe:='ebe':crem_dat:=0; if nops(boc)>0 then bec:=[bac]:bic:=[boc]:p_at_inf_c:=[[]];p_at_inf_w:=[[]]; fi; if cur=new then Nc:=-1:Kod:=[]:qR:=-1:rc:=-1:r4:=-1:r1:=-1:dill:=0: reg:=1:shaf:=0:SFE:=0:tuQ:=[]:sB:=1:sBb:=0:cP:=1:misc:=table(): Ncl:='Ncl':qqR:=1:cplist:=0;lowB:='lowB';uppB:='uppB'; if iee=1 then if uuV>1 then lprint(`By Laska's algorithm, minimal Weierstrass coefficients are`); print(yo);buc(op(boc));#t_ein:=[Wr/uV^2,Ws/uV,Wt/uV^3,Wu/uV]; else lprint(`Curve is Weierstrass minimal.`); fi; Tate(); else tate(); fi; ncur:=nops(et1)+1:nisog:=[]:disog:=[]:visog:=[]:et2:=[op(et2),yo]: for y while member(cat(alf[y],Nc),et1) do od; if pcn>0 then ya:=eT[pcn]:rc:=max(rc,ya[7]):r4:=max(r4,ya[8]):SFE:=ya[17]: r1:=ya[9]:ya:=[pcn,op(ya[11])]: for jm to nops(ya) do jn:=ya[jm]:nisog:=[op(nisog),jn]:disog:=[op(disog),0]: visog:=[op(visog),[]]:e:=eT[jn]:p:=[op(e[11]),ncur]:z:=[op(e[12]),0]: d:=[op(e[13]),[]]:eT[jn]:=uae(p,z,d,e,11): od; pcn:=0: fi; if y=nops(alf) then for z to 26 do alf:=[op(alf),cat(alf[nops(alf)-25],`&`)]: od; fi; cur:=cat(alf[y],Nc):et1:=[op(et1),cur]: tW:=[[],[]]:e4:=c4:e6:=c6:lds:=[]:se4:=1/2:x:=yo; if storE then yW:=1; else yW:=ncur; fi; while yW0 then e:=eT[jk][20]: if nops(tW[1])=0 then tW:=[[a,op(map(proc(jm) qfr(a*jm,2) end,e[1]))],[jk,op(e[2])]]: fi; e:=[[op(e[1]),a],[op(e[2]),ncur]]:eT[jk]:=uae(e,eT[jk],20): fi; od; bc(op(x)); for jk in oRR do if not type(jk,list(rational)) then next;fi; if on(jk) then ap1_(jk);fi; U:=jk[1]:V:=jk[2]:t:=bec:U:='U':V:='V': if member(0,t,'d') then trcw(jk,d);fi; od; eT[ncur]:=[Ncl,[Kod],-1,[],-1,[],rc,r4,r1,-1,nisog,disog,visog,0,[\ ],0,SFE,[],[1,0],tW,[cP,cplist],copy(misc)]:etj:=[op(etj),jay]: fprint(`This new curve has now been entered into the apecs catalog.`); fi; stac:=nops(ft1)+1:ft1:=[op(ft1),cur]:ft2:=[op(ft2),yo]: eprint(cus,cur=yo); cus,cur end: ein_tr:=proc(z) local t; t:=[(z[1]*uV^2-Wr)/Wu^2,(z[2]*uV^3-Ws*z[1]*uV^2+Ws*Wr-Wt)/Wu^3]; if nops(z)=3 then t:=[op(t),z[3]];fi; t end: eiop:=proc() local i,k,T_,uh,xh1; global bh,bbh,c,c_l,co_list,dill,dol,hx,L1,L3,L2,L4,Ld,Lt,Lu,lquar_,\ r3,R3,T,T1,T2,tut,u,uflag,um,x,xh,yh; #internal - used by rka and rkb to make the two calls to diop and #to find rational points on homogeneous spaces. uflag:=2:lquar_:=lquar;hx:=0; if Rflag<5 then diop(L1,L3,AA1,BB1,1);L1:=Lt:L3:=Lu:dol:=[co_list]; if r1=-1 or nops(RR)0 and iee>0 then lprint();lprint(`We now search for points up to naive height`,dil); lprint(`on the`,tut,` uncertain quartics associated to the 2-descent`); print(cur,` ---> `,naic); if q1<>0 then x:='x';lprint(`with respect to the shifted coordinate `,x-q1); fi; fi; for yh to dil while nops(dol[1])+nops(dol[2])>0 and nops(RR)dill then xh1:=1; else xh1:=dill+1; fi; for xh from xh1 to dil do if igcd(yh,xh)>1 then next;fi; T1:=xh/yh; for T_ to 2 do T1:=-T1: for uh to 2 do T:='T';co_list:=dol[uh]; if nops(co_list)=0 then next;fi; Lt:=L.uh;Lu:=L.(2+uh);c_l:={};T:=T1; for i to nops(co_list) do if essqr(co_list[i]) then c_l:=c_l union {i};T:='T';x:=co_list[i]; k:='k';c:=seq(coeff(x,T,5-k),k=1..5); T2:=[T1,esqrt];fprint(`Found point`,T2,` on`,x=Y^2);hx:=1; x:=co_list[1];bh:=coeff(x,T,4);bbh:=qfr(coeff(x,T,0),2); if i=1 then if uh=1 then um:=um1/(2*yh):x:=q1+bh*((um*xh)^2): else um:=(um1*um2)/(4*yh):x:=q1+((um*esqrt*yh^2/xh)^2): fi; ap1_(absc(x)): else writeto(garb);quar(c,op(T2));writeto(terminal);diyp(); fi; Ld:=map(proc(e) e*bh/igcd(e,bh)^2 end,Lt): R3:=R3/nops(Lt):Lt:= Lt union Ld:R3:=R3*nops(Lt):r3:=tam(R3): T:=T1;L.uh:=Lt;break; fi; od; od; x:=[];T:='T'; for i to nops(co_list) do if not member(i,c_l) then x:=[op(x),co_list[i]];fi; od; if uh=1 then dol:=[x,dol[2]]; else dol:=[dol[1],x]; fi; T:=T1; od; od; T:='T'; od; if hx=0 and nops(dol[1])+nops(dol[2])>0 and iee>0 then lprint(`No points found.`);lprint(); fi; dill:=max(dill,dil); pra();uflag:=0: if nops(dol[1])+nops(dol[2])=0 or nops(RR)=r4 then RETURN(); elif iee>0 and q1=q_1[nops(q_1)] then lprint(`The quartics a*T^4+..+e=Y^2 which are loc. sol. but for which`); lprint(`no rational point was found are listed (with other info) in dolt.`); lprint(`To continue the search use RkNC(args) where args is a number >`); lprint(`the present value`,dill,` of dill`); fi end: Ell:=proc() global iee; iee:=1:ell(args);iee:=0:NULL end: ell:=proc() global K_, oK_; #input (general) weierstrass coefficients oK_:=K_;K_:=0; if oin(args)>0 then RETURN();fi; if K_=1 then K_:=0:eprint(` K_ set to 0 --- K_=1 reserved for ein (see Menu(ell))`); fi; ell1() end: ell1:=proc() local i, y; global bouu,bouv,Dn,ft1,ft2,sB,sBb,stac,zst; #called both by ell and ein, the latter when K_<>1, i.e. the ground field <>Q #or when E is singular, i.e. DD=0. bouu:=false;bouv:=true;Dn:=0; ill();y:=yo;zst:='zst'; if K_<>1 then y:=[op(y),K_];fi; if member(y,ft2,zst) then if zst<>psn then stec():fi; eprint(cus,yo,` field label K_`=K_):RETURN(yo): fi; sB:=0;sBb:=0; ft1:=[op(ft1),CUR]:stac:=nops(ft1):ft2:=[op(ft2),y]: if sT=4 then eprint(`Curve is singular.`); else RR_update(); fi; yo end: Emod:=proc(p) global iee; iee:=1:emod(p);iee:=0: if card>0 then print(`card(E mod `.p.`) = `.card);fi; NULL end: emod:=proc(p) global card; #find order of group of points of cur mod p if bouv then card:=0;aM();RETURN();fi; if not type(p,integer) or p<=0 or not isprim(p) then print(`args of emod must be a prime number`);card:=0:RETURN(); fi; emoe(p);card end: Emods:=proc() global iee; ; iee:=1:emods(args);iee:=0:NULL end: emods:=proc() global z, ap; local p,p1,p2,x,y,a,d,j; #produce table of number of points etc. of cur mod p where p1<=p<=p2 if bouv then aM();RETURN();fi; if nargs=2 then p1:=args[1]:p2:=args[2]: else p1:=op(1,args):p2:=op(2,args): fi; p:=p1-1:a:=evalf(180./Pi):print(); if p2<11 then lprint(`p N ap theta DD j`); elif p2<101 then lprint(`p N ap theta DD j`); else lprint(`p N ap theta DD j`); fi; for x do p:=nxtprm(p): if p<=p2 then d:=modp(DD,p):emoe(p); if d=0 then j:=`-`: if spur=0 then z:=`singular-additive`: elif card=p then z:=`singular-split multiplicative`: else z:=`singular-nonsplit multiplicative`: fi; else j:=modp(c4^3/d,p): if modp(card,p)=1 then z:=supersingular: elif modp(card,p)=0 then z:=anomalous: else z:=NULL: fi; fi; ap:=p+1-card:y:=evalf(a*arccos(-ap/(2*sqrt(p)))): lprint(p,card,ap,evalf(y,4),d,j,z); else break; fi; od end: emoe:=proc(p) global spur, c6s, card; local u,x,a,b,j0; #internal - used by allp,eao,ebo,emod,emods,tate and tor #calculate spur=trace of Frobenius of E mod p and card=p+1-spur if p>100 and tFr3[1]=101 then frobr():fi; if p=2 then spur:=modp(a6,2)+modp(a1+a2+a3+a4+a6,2)-modp(a2+a4+a6,2)-modp(a3+a6,2): elif p=3 then spur:=0: for x to 3 do spur:=spur-mods(b2*x^2+(1-b4)*x+b6,3):od: elif p>=tFr3[1] and modp(DD,p)>0 then spur:=newp(p): elif modp(DD,p)=0 then member(p,badp,j0):spur:=badt[j0]: else j0:=modp(jay,p):c6s:=mods(c6,p): if modp(j0-1728,p)=0 then if modp(p,4)=3 then spur:=0: elif p=5 then spur:=2*mods(c4/2,5): elif p<100 then spur:=mods(froB0[p][j0]*power(c4,1/4*p-1/4),p): elif p<855 then spur:=mods(tFr[p][j0]*power(c4,1/4*p-1/4),p): elif p<1350 then spur:=mods(tFr2[p][j0]*power(c4,1/4*p-1/4),p): else spur:=mods(tFr3[p][j0]*power(c4,1/4*p-1/4),p): fi; elif j0=0 then if modp(p,6)=5 then spur:=0: elif p<17 then spur:=froB0[p][p][modp(c6,p)]: elif p<100 then spur:=mods(froB0[p][p]*power(c6s,1/6*p-1/6),p): elif p<855 then spur:=mods(tFr[p][p]*power(c6s,1/6*p-1/6),p): elif p<1350 then spur:=mods(tFr2[p][p]*power(c6s,1/6*p-1/6),p): else spur:=mods(tFr3[p][p]*power(c6s,1/6*p-1/6),p): fi; elif p<100 then spur:=froB0[p][j0]*mods(power(c6s,1/2*p-1/2),p): elif p<855 then spur:=tFr[p][j0]*mods(power(c6s,1/2*p-1/2),p): elif p<1350 then spur:=tFr2[p][j0]*mods(power(c6s,1/2*p-1/2),p): else spur:=tFr3[p][j0]*mods(power(c6s,1/2*p-1/2),p): fi; fi; card:=p+1-spur:spur end: eon:=proc() local i, rr, ux; global badmp, badp, cP, cplist, cur, dill, disog, Dn, dolt, Dt, hN, Kod, misc, MM, Nc, Nce, Ncl, nisog, NN, omega, omegA, Omega, pcn, PP, qqR, qR, rc, reg, RKI, RKJ, RR, r4, r1, sB, sBb, SFE, shaf, tuQ, tW, visog; #(internal)used by ein and stec: fetch data from catalog. cur:=et1[ncur]:ux:=eT[ncur]:hN:=ux[1]: if type(hN,integer) then if zst='zst' then Nc:=hN:badmp:=0:efac(DD);Dt:=ntab:Dn:=nP:tate();fi; uae(uae(Ncl,ux,1),eT,ncur): elif zst='zst' then Nc:=hN[1]:badp:=hN[2]:Nce:=hN[3]:Ncl:=hN: efac(DD,1):Dt:=ntab:Dn:=nP:bad(); fi; Kod:=op(ux[2]):qR:=ux[5]:MM:=nops(RR); rc:=ux[7]:r4:=ux[8]:r1:=ux[9]:dill:=ux[10]: if type(dill,list) then RKI:=dill[1];RKJ:=dill[2];dolt:=dill[3]; if type(dolt[nops(dolt)],integer) then dill:=dolt[nops(dolt)]; else dill:=0; fi; fi; nisog:=ux[11]:disog:=ux[12]: visog:=ux[13]:reg:=ux[14]:omegA:=ux[15]: if type(omegA,list) and nops(omegA)>0 then i:=omegA[1]; if nops(i)=3 then Omega:=i[1];omega:=copy(i[2]); elif nops(i)=2 then omegA:=omegA[2..nops(omegA)];om(); fi; fi; shaf:=ux[16]:SFE:=ux[17]:tuQ:=ux[18]:sB:=ux[19]: if type(sB,list) then sBb:=sB[2];sB:=sB[1]; else sBb:=0; fi; qqR:=sign(sB):sB:=abs(sB):tW:=ux[20]: misc:=copy(ux[22]):pcn:=0: if type(ux[21],list) then cP:=ux[21][1]:cplist:=ux[21][2]: else cP:=ux[21]:cplist:=0: fi; RR_update(); if RR<>ux[6] and qR>-1 then upis();fi end: eprint:=proc(); #Capitalized apecs commands set iee:=1 or higher; the default value is 0; if iee>0 then print(args); fi end: esr:=proc(r) local SS, z; global x, y; #used by tor: a new point is added to PP if of finite order #or to RR if of infinite order. x:='x':SS:=[isolve(alz=(4*r)^2)]: for z to nops(SS) while NNt then z2:=[x,y]: if not member(z2,RR) and not member([x,-y-a1*x-a3],RR) then RR:=[op(RR),z2]:MM:=nops(RR):qqR:=1: fi; else if NN=1 then PP:=Qq:NN:=ord: else syb({PP},{Qq});PP:=op(sy):NN:=nops(sy)+1: fi; fi end: essqr:=proc() local n,d,p,s,sol,t; global bsq, esqrt, ssn, ssd, ssq; #if args is in Q (resp. Q(T), Q(sqrt(m_)) then if args is a square in that #field, set bsq:=true and esqrt:=its square root; otherwise set bsq:=false. #If args is not in any of those fields then set bsq:=true and #esqrt:=sqrt(args) (simplified) bsq:=false:ssq:=sex_(args[1]); if ssq=0 then bsq:=true:esqrt:=0:RETURN(bsq);fi; ssn:=numer(ssq):ssd:=denom(ssq): if K_=m_ then if type(ssq,polynom(rational,w_)) then p:=sex_((s+t*w_)^2);n:=coeff(p,w_,0);d:=coeff(p,w_,1); sol:=solve({n=coeff(ssq,w_,0),d=coeff(ssq,w_,1)}); for p in sol while not bsq do assign(p); if type([s,t],list(rational)) then bsq:=true;esqrt:=s+t*w_;break; fi; unassign('s','t'); od; fi; elif type(ssq,rational) then if ssq>0 then if issqr(ssn) and issqr(ssd) then bsq:=true:esqrt:=isqrt(ssn)/isqrt(ssd): fi; fi; elif type(ssq,ratpoly(rational,T)) and has(ssq,T) then readlib(psqrt);n:=psqrt(ssn); if type(n,polynom(rational,T)) then d:=psqrt(ssd); if type(d,polynom(rational,T)) then bsq:=true;esqrt:=n/d; fi; fi; else bsq:=true;esqrt:=simplify(sqrt(ssq),sqrt,symbolic); fi; bsq end: Expr:=proc() global iee; iee:=1:expr(args):iee:=0:NULL end: expr:=proc() global a, x, bqe, iee, t, l, j; local rr,z; if bouv then aM();RETURN();fi; z:=args:a:=nargs: if z=[] then eprint(`[] is a torsion point of order 1.`):RETURN([],[[]]): fi; if a>1 then z:=[z]:fi; if not on(z) then print(z,` not on`,cur);RETURN();fi; do_tor(); if member([z[1],z[2]],ouP,'i') then x:=[PP]:z:=x[i]:eprint(z,` is a torsion point of order`,z[3]): RETURN([],[z]): fi; bqe:=iee:iee:=0:lin_([op(RR),z]);iee:=bqe: if adet<>0 then fprint(z,` is independent of RR and will now be adjoined to RR:`); ap1_(z);RETURN(z); else if nops(TP)=0 then t:=NULL:l:=NULL: else t:=cat(`+`,TP):l:=cat(` (`,TP,` is a torsion point)`): fi; if den>1 then fprint(` `,z,` is only a rational, not an integral lin. comb. of the \ points in RR:`); eprint(cat(den,z)=cat(seq(cat(num.j,lsp[j],`+`),j=1..(MM-1)),num.MM,\ lsp[MM],t,l)): fprint(`Since this implies that RR is not a basis, we attempt to improv\ e RR.`); rr:=RR;ap1_(z); if rr=RR then eprint(`RR remains`=RR);fi; j:='j':RETURN(z); else eprint(z=cat(seq(cat(num.j,lsp[j],`+`),j=1..(MM-1)),num.MM,\ lsp[MM],t,l)): fi; fi; [num.(1..MM)],TP end: extend_L1:=proc(m) local d, h, hx, i, j, m0, p, s, x, y; global hens_d, L1; m0:=abs(L1[nops(L1)][1]);s:=[]; for i to L4 do p:=L3[i]; if p>m then break;fi; y:=NULL; for j to nops(L1) do d:=L1[j][1]; if modp(d,p)=0 then next;fi; x:=d*p; if x<=m0 or member(x,s) then next;fi; if x>m then break;fi; h:=L1[j][2];s:=[op(s),x]; if s_1=1 then if L2[i]=1 then if h=1 then y:=y,[x,1]; else y:=y,[x,-1]; fi; elif h=1 then y:=y,[x,-1]; else y:=y,[x,hens_(x)]; fi; else if L2[i]=1 then if h=1 then y:=y,[x,1],[-x,h_s]; elif h_s=1 then y:=y,[x,-1],[-x,-1]; else y:=y,[x,-1],[-x,hens_(-x)]; fi; else if h=1 then if h_s=1 then y:=y,[x,-1],[-x,-1]; else y:=y,[x,-1],[-x,hens_(-x)]; fi; else hx:=hens_(x); if h_s=1 then y:=y,[x,hx],[-x,hx]; elif hx=1 then y:=y,[x,1],[-x,-1]; else y:=y,[x,hx],[-x,hens_(-x)]; fi; fi; fi; fi; od; L1:=sort([op(L1),y],proc(a,b) evalb(abs(a[1])<=abs(b[1])) end); od; hens_d:=[]; for x in L1 do if x[2]=1 then hens_d:=[op(hens_d),x[1]];fi;od; hens_d end: ey:=proc(rr,dd) local fy, n; global de, Di, eay, eoo, euy, eyo, eyr, eyy, hex, hyy, j, n1, ne, nf, q, r; #calculate to "dd digits" the cusp form (various cases when r<0) or L^(r)(1) #(when r>=0 and the parity of r is even/odd when the sign of the functional #equation is +/- :this proc. will give faulty answers for r of wrong parity) if bouv then aM();RETURN();fi; r:=rr:de:=dd:Di:=Digits:Digits:=max(Di,dd+6): if r<0 then fy:=1; elif rc>=0 then if r0 then if hyy=-1 then j:=cat(`to within +-10^`,-dd); else j:=cat(`using at most `,hyy,` terms.`); fi; eprint(`Calculating the `.q.`L series at s=1 `,j); fi; hex:=0:eyt(eyr); for ne to nops(nisog) while hex=0 do if member(et1[nisog[ne]],ft1,'j') then eyt(ft[j][13]);fi; od; if hex=1 then eyp(0); if r>=0 then eyr:=[r,Digits,dd,eyy]:fi; Digits:=Di:RETURN(eyy); fi; #I don't remember why I had the following. Anyway it can pop up #unwanted in the sequence Rk();Sha(); so I've remarked it out. # if iee>0 then # hex:=readstat(`To proceed with the calculation anyway, enter yes; \ #else enter no;`); # fi; fi; eyo:=evalf(10^(-dd)/40):eay:=evalf(eoo*Pi/sqrt(Nc)):euy:=1:ne:=1: if r>-3 then eyy:=0:n1:=1:fi; for n from n1 while ne<50 and hyy<>n do while n>eai do eao();od; if eae[n]<>0 then if r<0 then euy:=evalf(exp(-eay*n))*eae[n]: else euy:=eae[n]*eya(eay*n)/n: fi; eyy:=eyy+euy: if abs(euy)0 then eyp(n);fi; if tim>0 then clk(); if tim<0 then hyy:=n;RETURN();fi; fi; od; eyp(n);eyy:=2*eyy: if eyy=0 then euy:=-de;else euy:=evalf(log10(abs(eyy)));fi; if euy<-de+2 then eyy:=0;eprint(`Sum deemed to be 0.`); elif abs(eyy)<10^5*eyo then eprint(`Must restart calc. of series with greater accuracy.`); ey(rr,dd+4);RETURN(); fi; if fy=0 and eyy<>0 then print(`Apecs error: sum should be 0.`);RETURN(); elif fy=1 and euy<-de+1 then print(`Apecs error: sum should be non-zero.`);RETURN(); fi; if r>=0 and n<>hyy then hyy:=-1;eyr:=[r,Digits,dd,eyy]:fi; if n=hyy and iee>0 then lprint(`Series truncated at `,hyy,` terms, so accuracy lower than in\ dicated.`); fi; Digits:=Di:eyy end: eya:=proc(x) local a,b,d,m,n,t; global eye, eyi, ts1, ts1d, yia, yie, yoa, yoe, yio, s; #(internal) used by ey - part of evaluation of L series etc. a:=.4343*x;b:=evalf(r*log10(x));d:=Digits; if a+b>de+.5 then eye:=0;RETURN(0);fi; if r=0 then eye:=evalf(exp(-x)): else t:=evalf(-ln(x)): if r=1 then if x<3 then eye:=evalf(t-gamma):eyi:=x: for n while abs(eyi)>eyo do eye:=eye+eyi:eyi:=-eyi*x*n/((n+1)^2): od; else yia:=0:yie:=evalf(exp(-x)):yoa:=1:yoe:=x:eye:=yie/yoe:eyi:=0: for n while abs(eye-eyi)>eyo do yio:=yie:yie:=yie+n*yia:yia:=yio:yio:=yie:yie:=x*yie+n*yia:yia:=yio: yio:=yoe:yoe:=yoe+n*yoa:yoa:=yio:yio:=yoe:yoe:=x*yoe+n*yoa:yoa:=yio: eyi:=eye:eye:=yie/yoe: od; fi; else s:='s': if d>ts1d or r>=degree(ts1) then fprint(`Apecs must take a long time-out to calculate a more precise`); fprint(`Taylor series for GAMMA(1+s).`); ts1:=evalf(taylor(GAMMA(1+s),s,r+1));ts1d:=d: fi; eye:=evalf(sum((t^m)*coeff(ts1,s,r-m)/(m!),m=0..r)): if max(0,a-b)+de>d-2 then Digits:=trunc(a-b+de+4);fi; eyi:=((-1)^(1+r))*x: for n while abs(eyi)>eyo do eye:=eye+eyi:eyi:=-eyi*x*((n/(n+1))^r)/(n+1): od; fi; fi; Digits:=d;eye end: eyp:=proc(m) local x,y,z;global eyy,hyy; x:=eyy; if m=0 then z:=`Series already calculated: `;hyy:=-1: elif nargs=2 then z:=m;print(); else z:=cat(m,` terms give`);x:=2*x: fi; if x=0 then fprint(z,` 0.00`);RETURN(); else y:=evalf(log10(abs(x))); fi; if y<-de then fprint(z,evalf(x,2),` which is deemed to be 0.`); eyy:=0; else fprint(z,evalf(x,de+trunc(y)+1)); fi end: eyt:=proc(x) global hex, eyr, eyy; #(internal) used by ey - has series been done (in part) before? if type(x,list) and r=x[1] and Digits<=x[2] and de<=x[3] then hex:=1:eyr:=x:eyy:=x[4]; fi end: FnL:=proc() local i; global iee; if bouv then aM();RETURN();fi; if nargs=1 then i:=args:else i:=args[1]:fi; if rc>=0 and i sum=0;`); i:=readstat(`to continue anyway, type y; else n;`); if i=`n` then RETURN();fi; fi; iee:=1:fnL(args):iee:=0:NULL end: fnL:=proc() local d, i; global q, NN, hyy, niv, tim; #Calculate the i-th derivative of L(s) at s=1; i must have have the same #parity as SFE (the sign of the functional equation) if bouv then aM();RETURN();fi; if SFE=0 then eprint(`SFE must first be determined -- we call roha():`):roha(); fi; if nargs=1 then i:=args:d:=4;hyy:=-1; elif nargs=2 then i:=args[1]:d:=args[2]:hyy:=-1; elif nargs=3 then i:=args[1]:d:=args[2]:hyy:=args[3]; if not type(hyy,posint) then print(`Third arg must be a posint.`,qtr);RETURN(); fi; else print(qts);RETURN(); fi; if not type(d,posint) then print(`Second arg must be a posint.`,qtr);RETURN(); fi; if modp(i-(1-SFE)/2,2)=1 then if SFE=1 then q:=`even`:else q:=`odd`:fi; lprint(`Only programmed to cal. the i-th deriv. of L(s) at s=1 for i `,q); lprint(`because the sign of the functional equation is `,SFE); lprint(qtr);RETURN(); fi; ey(i,d); if tim<0 then if niv=1 then niv:=0:tim:=0:fi; RETURN(); elif iee>0 then if eyy=0 then q:=` `:else q:=` approx. `:fi; eyp(`L^`.i.`(1)/`.i.`! is`.q,2); fi; eyy end: Ford:=proc() local i,Q; ford(args); if ord=0 then RETURN();fi; if ord>ort then Q:=Qq,etc:else Q:=Qq,O:fi; for i to ord do lprint(i,`*z = `,Q[i]); od end: ford:=proc() global x, y, ord, ort, Qq; local x2,y2,z,Q1; # list multiples z,2z,..up to tz of point z=[x,y] # where t=min{order(z),args[3]}. ort:=args[nargs]: if nargs=3 then x:=args[1]:y:=args[2]: elif nargs=2 and type(args[1],list) then if args[1]=[] then ord:=1;Qq:=O;RETURN(); else x:=args[1][1]:y:=args[1][2]: fi; else ord:=0:lprint(\ `ford takes args z,n where z is a point on E (z=x,y or [x,y] or [x,y,d])`); lprint(`and n is a pos. integer.`,qtr);RETURN(); fi; z:=[x,y]: if not on(z) and not hastype(z,float) then ord:=0:print(z,noc);RETURN(); elif sT=4 and z=zs then print(`The singular point`=zs,` is not allowed as an arg.`,qts); ord:=0;RETURN(); fi; x2:=x:y2:=y:z:=simplify(normal(-y-a1*x-a3)):Q1:=[x,y,0]: for ord from 2 to ort do if x<>x2 or x=x2 and y2<>z then eadd(x,y,x2,y2):x2:=x3:y2:=y3: Q1:=Q1,[x3,y3,0]: else Qq:=[x,y,ord]: for z from 2 to ord-1 do Qq:=Qq,[op(1..2,Q1[z]),ord/igcd(ord,z)]: od: break; fi; Qq:=Q1: od end: fprint:=proc(); #(internal) like eprint: lprint(args) when the flag iee has value >0 if iee>0 then lprint(args); fi end: frobr:=proc(); #(internal) read in trace of Frobenius tables fprint(`Momentary delay while reading in tables Frobe1, Frobe2 and Frobe3.. `); read Frobe1:read Frobe2:read Frobe3: fprint(`The tables have been read.`) end: Genj:=proc() global iee; iee:=1:genj(args):iee:=0:NULL end: genj:=proc() local j,k;global a1, a2, a3, a4, a6, a5; #Input "generic" elliptic curve with jay=args j:=args[1]; if nargs=2 then k:=args[2]; else k:=1; fi; if j=0 then ein(A27,k); elif j=1728 then ein(A32,k); else a1:=1:a2:=0:a3:=0:a4:=-36/(j-1728):a6:=-1/(j-1728):a5:=a6: if type(j,rational) then ein(a.(1..5),k); else ell(a.(1..5),k); fi; fi end: Go:=proc() global iee; iee:=1:go(args);iee:=0:NULL end: go:=proc() local a,i,n,q,y; #point to curve #args in stack, or list stack if no args n:=nops(ft1):a:=args: if n>0 then if nargs=0 then for i to n do y:=ft2[i]; if nops(y)=5 then y:=[op(y),1];fi; if i=stac then q:=`present curve`: else q:=``: fi; if ft1[i]=CUR then fprint(i,y[1..5],`field label K_`=y[6],q); else fprint(i,ft1[i],q); fi; od; elif nargs=1 and type(a,posint) and a<=n then if a=stac then eprint(`curve remains `.cur); elif nops(ft2[a])=6 then ell(ft2[a]): else ein(ft2[a]); fi; else print(`improper args in command go`); fi; else eprint(`stack is empty`); fi; NULL end: Half:=proc() global iee; iee:=1:print(half(args)):iee:=0:NULL end: half:=proc() local a, qu, t, w, z; global q, x, y; #Calc. list q of points P such that 2*P=args. z:=args:q:=[]; if nargs>1 then z:=[z];fi; if not on(z) then print(z,` not on`,cur);RETURN();fi; if sT=4 and z=zs then print(`The singular point`=zs,` is not allowed as an arg.`,qts); RETURN(); fi; if z=[] then x:='x'; if bouu or sT=0 or sT=1 then do_tor();q:=[]; if NN<1 then RETURN();fi; for z in [PP] do if z[3]=2 then q:=q,z:fi: od; elif type(convert(prac,RootOf),polynom(algnum,x)) then for z in roots(convert(prac,RootOf)) do q:=q,absc(z[1]); od; else for z in [solve(prac,x)] do q:=q,absc(z); od; fi; else x:=z[1];y:=z[2];q:=NULL;qu:=convert(QuarQ,RootOf); if type(qu,polynom(algnum,M)) then w:=[]; for z in roots(qu) do w:=[op(w),z[1]]; od; else w:=[solve(QuarQ,M)]; fi; for t in w do if (sT=0 or sT=1) and not type(t,ratpoly(rational,T)) then next;fi; a:=(t^2-b2/4-x)/2:q:=q,[a,t*(a-x)-y-a1*(a+x)/2-a3]; od; fi; q:=sex_([q]) end: Hdif_Sik:=proc(); lprint(`Siksek's lower bound for h^-h/2`=hdif_Sik()) end: hdif_Sik:=proc() global D_v,D1_v,e_v,f,f1,g,g1,i,j,k,mU,lowB,p,r,r_,r1_,s,t,u,x; if lowB<>'lowB' then RETURN(lowB);fi; if not type(yo,list(rational)) or sT=4 then RETURN(` undefined --- E must be /Q and nonsingular`); fi; lowB:=0;tate();raf();x:='x';f:=4*prac;g:=x^4-b4*x^2-2*b6*x-b8; u:='u';x:=1/u;f1:=expand(u^4*f);g1:=expand(u^4*g):x:='x';u:=x; r:=[op(raF),fsolve(g),fsolve(f+g),fsolve(f-g),fsolve(diff(f,x)),\ fsolve(diff(g,x))]; r_:=NULL; for j in r do if j>=-1 and j<=1 then r_:=r_,j;fi; od; r:=[fsolve(f1),fsolve(g1),fsolve(f1+g1),fsolve(f1-g1),\ fsolve(diff(f1,x)),fsolve(diff(g1,x))]; r1_:=NULL; for j in r do if j>=-1 and j<=1 then r1_:=r1_,j;fi; od; for i in cplist do p:=2;vam(i[2]); if zz>1 then mU:=1/3: elif i[2]<0 or i[2]=2 then mU:=1/4; else mU:=(1-1/4^v)/3: fi; p:=i[1];e_v:=p^(2*trunc(vam(DD)/2)); lowB:=lowB+evalf(mU*log(e_v)); od; r:=raF[1]: if DD<0 then if r<=-1 then D_v:=[[-1,1]]:D1_v:=[[-1,1/r],[0,1]]: elif r<=1 then D_v:=[[r,1]]:D1_v:=[[0,1]]: else D_v:=[]:D1_v:=[[0,1/r]]: fi; else s:=raF[2];t:=raF[3]; if r<=-1 then if t<=1 then if s>=-1 then D_v:=[[-1,s],[t,1]]:D1_v:=[[-1,1/r],[0,1]]: elif t>-1 then D_v:=[[t,1]]:D1_v:=[[1/s,1/r],[0,1]]: else D_v:=[[-1,1]]:D1_v:=[[-1,1/t],[1/s,1/r],[0,1]]: fi; else if s<-1 then D_v:=[]:D1_v:=[[1/s,1/r],[0,1/t]]: elif s<1 then D_v:=[[-1,s]]:D1_v:=[[-1,1/r],[0,1/t]]: else D_v:=[[-1,1]]:D1_v:=[[-1,1/r],[0,1/t],[1/s,1]]: fi; fi; else if s<1 then D_v:=[[r,s]]:D1_v:=[[0,1/t]]: elif r<=1 then D_v:=[[r,1]]:D1_v:=[[0,1/t],[1/s,1]]: else D_v:=[]:D1_v:=[[0,1/t],[1/s,1/r]]: fi; fi; fi; i:=NULL; for j in D1_v do i:=i,j[1],j[2]; for k in [r1_] do if k>=j[1] and k<=j[2] then i:=i,k;fi; od; od; k:=NULL; for j in {i} do x:=j;k:=k,max(f1,g1);od; t:=min(k);i:=NULL; for j in D_v do i:=i,j[1],j[2]; for k in [r_] do if k>=j[1] and k<=j[2] then i:=i,k;fi; od; od; if i<>NULL then k:=NULL; for j in {i} do x:=j;k:=k,max(f,g);od; t:=min(k,t); fi; lowB:=(evalf(log(t)/3)-lowB)/2 end: Hdif_Sil:=proc(); if bouv then aM();RETURN();fi; lprint(`Silverman's lower and upper bounds for h^ - h/2 are [lowB_Sil,uppB] =`); print(hdif_Sil());NULL end: hdif_Sil:=proc() global hb2, logaj, muE, uppB, lowB_Sil, i, p, j; ; if bouv then aM();RETURN();fi; if b2=0 then hb2:=0; else hb2:=(loga(b2/12)+log(2))/2: fi; logaj:=loga(jay):muE:=(loga(DD)+logaj)/12+hb2: uppB:=evalf(muE+1.07): lowB_Sil:=-logaj/24-muE-.973: if jay<>0 then for i from 2 to Dn do if Dt[i][2]=0 then next:fi; p:=Dt[i][1]:vam(c4):j:=3*v-Dt[i][2] : if j>=0 then next; elif j=-1 and v=0 then lowB_Sil:=lowB_Sil+log(p)/12: else lowB_Sil:=lowB_Sil+j*log(p)/24: fi; od; fi; lowB_Sil:=evalf(lowB_Sil):[lowB_Sil,uppB] end: Heeg:=proc() global iee, ubnh; iee:=1;print(heeg(args));iee:=0; if hht>0 then print(`is the approx. canonical height of a Heegner point.`); hdif_Sik();ubnh:=evalf(2.01*(hht-lowB),5): lprint(`Using hdif_Sik (see Menu), the naive height of this point`); lprint(`(log(max(|numer|,denom))) is at most`); print(` ubnh`=ubnh); lprint(`and therefore will be found by Seek(x,1) where`); print(` x`=trunc(exp(ubnh))); lprint(`The purpose of the 2nd arg. in seek is to terminate seek`); lprint(`as soon as a point is found -- see Menu(seek).`); fi; NULL end: heeg:=proc() local n, q; global hht, iee; hht:=0; if nargs>1 or (nargs=1 and not type(args,posint)) then RETURN(`Improper args -- see Menu(heeg).`); fi; q:=` -- must be 1 for the proc. heeg.`; if K_<>1 then RETURN(cat(`The field label K_`,q));fi; if SFE=0 then if iee=1 then Roha();iee:=1;else roha();fi; fi; if SFE=1 then RETURN(cat(`The rank of E(Q) is even`,q));fi; if rc>=0 and r4>1 then RETURN(cat(`The rank of E(Q)`=r4,q));fi; tor();om(); if MM>1 then RETURN(cat(`The rank of E(Q) >= |RR|`=MM,q));fi; if rc=-1 then eprint(`The rank of E(Q) has not been determined; at least we know`); eprint(`it is odd. We proceed on the assumption that it is 1.`) fi; if nargs=1 then n:=args else n:=50000;fi; fnL(1,4,n);fprint(); if hyy>0 then eprint(`Not enough terms in the series to estimate L'(1).`); RETURN(cat(`Try Heeg(args) with args >= ,say,`,n+50000)); fi; if eyy<.001 then if eyy<-.001 then q:=` < 0.`; elif rc=-1 then RETURN(`Computed value of L'(1) approx. 0 -- seems rank of E(Q)>=3.`); else q:=` approx. 0.`; fi; eprint(cat(`Something wrong -- computed value of L^{1}(1)`,q)); cat(`Try Heeg(args) with args >`,n); else hht:=evalf(eyy*NN^2/2/Omega/cP,4); fi end: hens_:=proc(d) local hex, vh; if nrs(d,bb2)=-1 then RETURN(-1);fi; hex:=1;set_(d,0,AA1,0,BB1/d); for vh in SS while hex=1 do if (vh>2 and modp(d,vh)>0 and L(d,vh)=1) or \ (vh=2 and (modp(d,8)=1 or modp(d+4*AA1,8)=1)) then next; fi; hex:=-1+2*Zp_sol(vh); od; hex end: Ht:=proc() global iee; iee:=1:print(ht(args));iee:=0;NULL #print is needed when ht(args) is remembered by ht__ end: ht:=proc() local ape, hz; global iee; #Calculate the canonical height if bouv then aM();RETURN();fi; if nargs>1 then hz:=[args[1],args[2]]: elif type(args,rational) then hz:=absc(args): else hz:=args: fi; if hz<>[] and not type(hz,list(rational)) then print(hz,` not rational`,qtr);RETURN(); fi; if not on(hz) then print(hz,noc);RETURN();fi; ape:=iee:iee:=iee*(1-pcn):do_tor();iee:=ape; ht__(stac,hz,Digits); end: ht_:=proc() ht__(stac,args,Digits); end: ht__:=proc() local a, b, BB, c, ght, hz, j, k, lL, m, n, T, w, x1, x2, z; global b2p, b4p, b6p, b8p, dec, Di, HH, iee, Kud, lht, p, sing, sw, u, x, y; option remember; hz:=args[2]; if member(hz,ouP) or hz=[] then ght:=0:fprint(hz,` has finite order; global height = 0`);RETURN(0); fi; Di:=Digits:x:=hz[1]:y:=hz[2]: HH:=max(4,abs(b2),2*abs(b4),2*abs(b6),abs(b8)): dec:=evalf(5*Di/3+1.5+.75*log(7+4*log(HH)/3)): Digits:=trunc(evalf(Di+2+(log(dec))/log(10))): b2p:=b2-12:b4p:=b4-b2+6:b6p:=b6-2*b4+b2-4:b8p:=b8-3*b6+3*b4-b2+3: if abs(x)>=1/2 then T:=1/x:b:=1: else T:=1/(1+x):b:=0: fi; lL:=evalf(-(log(abs(T)))/2):m:=0:sw:=0: for n from 0 to dec do if b=1 then w:=evalf(T*(T*(T*(T*b6+2*b4)+b2)+4)): z:=evalf(T*(T*(T*(-b8*T-2*b6)-b4))+1): if abs(w)<=2*abs(z) then m:=evalf(m+(.25^n)*log(abs(z))):T:=w/z: else m:=evalf(m+(.25^n)*log(abs(z+w))):T:=w/(z+w):b:=1-b:sw:=sw+1: fi; else w:=evalf(T*(T*(T*(T*b6p+2*b4p)+b2p)+4)): z:=evalf(T*(T*(T*(-b8p*T-2*b6p)-b4p))+1): if abs(w)<=2*abs(z) then m:=evalf(m+(.25^n)*log(abs(z))):T:=w/z: else m:=evalf(m+(.25^n)*log(abs(z-w))):T:=w/(z-w):b:=1-b:sw:=sw+1: fi; fi; od; lL:=lL+(m/8):eprint(`number of switches in Silverman's algorithm = `.sw); eprint(`Table of local heights ((log D)/12 term omitted)`); eprint(` and Kodaira type if P is singular mod p:`);eprint(`arch`,lL); lht:=lL:ght:=lL:a:=x*(3*x+2*a2)+a4-a1*y:b:=2*y+a1*x+a3: c:=x*(x*(x*(3*x+b2)+3*b4)+3*b6)+b8:x1:=denom(x):k:=0:j:='j':Kud:=[Kod]: for j to Dn do if Dt[j][2]=0 then next;fi; p:=Dt[j][1]:x2:=vam(x1):x1:=zz:BB:=val(b):val(a);k:=k+1:sing:=Kud[k]: if v>0 and BB>0 then val(c4); if v>0 then val(c); if v=1/2 or v>=3*BB then lL:=-BB/3: else lL:=-v/8: fi; else z:=Dt[j][2]: if BB=1/2 then lL:=-z/8: else u:=min(BB,z/2):lL:=u*(u-z)/(2*z): fi; fi; else lL:=x2/2: fi; lL:=evalf(lL*log(p)): eprint(p,lL,sing); lht:=lht,lL:ght:=ght+lL: od; lL:=evalf((log(x1))/2): eprint(`the sum of the remaining local heights= `,lL); ght:=ght+lL:Digits:=Di: for z in ouP do ht__(stac,eadd(hz,z),Digits):=ght;od; if pcn=0 then eprint();eprint(`global height of`,hz,` is`,ght);fi; ght end: ian:=proc() global a, ncir; #internal - used by Dat, dat,We and we to check ID of curve if (nargs<>1 and nargs<>5) or (nargs=1 and nops(args[1])<>1 and \ nops(args[1])<>5) then print(`args must be a catalog name (eg A11) or a catalog number (eg 1)`); print(`or a sequence or list of five Weierstrass coefficients`); lprint(`possibly augmented by the field code K_ -- see Menu(ell)`); RETURN(2); fi; a:=args:ncir:='ncir': if nargs=1 then if type(a,integer) then if 1<=a and a<=nops(et1) then ncir:=a:fi; elif nops(a)=5 then member(a,et2,ncir): else member(a,et1,ncir); fi; else member([a],et2,ncir); fi; if not type(ncir,integer) then print(`curve not in apecs catalog`);RETURN(2); elif ncir=ncur then oan(); fi; 0 end: ill:=proc() local dd,x; global a5, bouu, bouv, f4, i, jey, K_, natsing, nats, p, raF, rc, sC4, sC6, xs, y, yo, ys, zs; #internal - used by Dat,ein,and ell1 to print curve data a5:=a6;yo:=[a.(1..5)]; for i from 2 by 2 to 8 do sB.i:=sex_(b.i):od; sC4:=sex_(c4):sC6:=sex_(c6): dd:=sex_((sC4^3-sC6^2)/1728): if dd<>0 then jey:=normal(sex_(sC4^3/dd)):fi; if iee>0 then if K_<>1 then print(`Field label K_`=K_);fi; print(`b's = `,sB2,sB4,sB6,sB8);print(`c's = `,sC4,sC6); if type(dd,integer) and dd<>0 then print('DD'=dd,` =`,afactor(dd)); else print('DD'=dd); fi; if dd=0 then print(`jay undefined`); elif type(jey,rational) then if member(jey,CMj,'y') then p:=CMd[y]:f4:=` CM by `.p;else f4:=NULL;fi; print(`jay = `,jey,` denom(jay) = `,afactor(denom(jay)),f4); else lprint(`jay = `,jey); fi; print_PP(); fi; if dd=0 then rc:=-1;raF:=[]; if K_=1 then K_:=0;fi; if sC4<>0 then xs:=sex_((18*sB6-sB2*sB4)/sC4):natsing:=node: nats:=[solve(x^2+a1*x-3*xs-a2,x)]: else xs:=-sB2/12:natsing:=cusp:nats:=[-a1/2]: fi; ys:=sex_(-(a1*xs+a3)/2):zs:=[xs,ys]: if iee>0 then lprint(`The singularity is at `,zs,` and is a`,natsing); if nops(nats)>0 then lprint(`Tangent slope(s) there = `,op(nats)); fi; fi; fi; if K_=1 and dd<>0 then bouu:=true;bouv:=false; else bouv:=true;bouu:=false; fi end: Inch:=proc(x); inch(x); if cur=x then print(`The catalog names `.y,x.` have been interchanged.`); print(`The catalog numbers remain the same.`); fi end: inch:=proc(x) global y, cur, et1, ft1; : #Interchange name of present curve with that of args; catalog numbers of #the two curves remain the same (used to make apecs notation agree with #other, e.g. Antwerp IV). if bouv then aM();RETURN();fi; if member(x,et1,'i') then y:=eT[i][1]: if type(y,list) then y:=y[1]:fi; if y<>Nc then print(`The two curves have different conductors - name switch disallowed.`); else y:=cur:cur:=x:et1:=uae(y,et1,i):et1:=uae(x,et1,ncur): if member(x,ft1,j) then ft1:=uae(y,ft1,j):fi; ft1:=uae(x,ft1,stac): fi; else print(`args must be name in apecs catalog like D162`); fi end: Isog:=proc() global iee, ncir, e7; ; if bouv then aM();RETURN();fi; iee:=2:isog();iee:=0:ncir:=ncur:e7:=qR:isug(nisog,disog) end: isog:=proc() global isn, iwog, k, uflag, vflag, storE, nisw, cisw, yw, ny, i, e2, e4, e6, e1, e3, iew, iee, qR, yo, nc, n, d, v, nb, xa, na, da, va, zw, oc, ne, l, vaj, eT, nisog, disog, visog, ftw, lcyi, cyi, p, j, y, z, hex, isl, ism, w, u, t, s, e, f; local a, g, h, q, STorE, x; #Find all curves isogenous to the present curve. if bouv then aM();RETURN();fi; if qR>0 then eprint(`isog already done`);RETURN();fi; isn:=ncur:iwog:={}:ttw(1);uflag:=1:vflag:=1:STorE:=storE:storE:=true: if atw<>1 then nisw:=otc[11]:cisw:=[]: for yw to nops(nisw) do ny:=nisw[yw]:x:=et2[ny];for i to 5 do g.i:=x[i]:od;e2:=g1^2+4*g2; e4:=g1*g3+2*g4;e6:=g3^2+4*g5;e1:=e2^2-24*e4:e3:=-e2^3+36*e2*e4-216*e6: iew:=iee:iee:=0:qR:=0:ein(0,0,0,-27*e1*(atw^2),-54*e3*(atw^3));iee:=iew: cisw:=[op(cisw),ncur]:yo:=[a1,a2,a3,a4,a6]: fprint(cat(`twist `,et1[ny],`*`,atw,` = `,cur)=yo); od; nisw:=[ntc,op(nisw)]:cisw:=[isn,op(cisw)]: for yw to nops(cisw) do nc:=cisw[yw]:x:=eT[nc]:n:=x[11]:d:=x[12]:v:=x[13]:nb:=nisw[yw]: xa:=eT[nb]:na:=xa[11]:da:=xa[12]:va:=xa[13]: for zw to nops(cisw) do oc:=cisw[zw]: if oc=nc then next;fi; ne:=nisw[zw]:member(ne,na,'j');l:=-abs(da[j]):vaj:=va[j]: if member(oc,n,'j') then if d[j]=0 then d:=uae(l,d,j):fi; else n:=[op(n),oc]:d:=[op(d),l]:v:=[op(v),[]]: fi; od; eT[nc]:=uae(n,d,v,x,11): od; nisog:=n:disog:=d:visog:=v:qR:=0: fi; ftw:=0:lcyi:=[]: if qR=0 then cyi:=1: for i to nops(nisog) do d:=abs(disog[i]): if isprim(d) then cyi:=cyi*d:fi; od; if member(-3,disog) then ftw:=1:fi; eprint(cat(cus,cur));itog(); else lcyi:=lcyo: if qfr(Nc,2)=Nc then lcyi:=lcyi minus {13}:fi; itog();p:=1:x:=NULL: for i to 6 do p:=nxtprm(p): if modp(cyi,p)>0 then for j to nops(lcyi) do if modp(lcyi[j],p)=0 then x:=x,lcyi[j]:fi; od; fi; od; lcyi:=lcyi minus {x}: fi; while nops(iuog)>0 do ein(iuog[1]);itog(); od; ein(isn); if atw=1 then for i in iwog do if i=ncur then x:=nisog:y:=disog:z:=visog: else x:=eT[i][11];y:=eT[i][12];z:=eT[i][13]; fi; hex:=0; for j in iwog do if i<>j then if not member(j,x) then x:=[op(x),j]:y:=[op(y),0]:z:=[op(z),[]]:hex:=1; fi; fi; od; if hex=1 then if i=ncur then nisog:=x:disog:=y:visog:=z: else eT[i]:=uae(x,y,z,eT[i],11); fi; fi; od; isl:=1: while isl=1 do isl:=0: for n in iwog do if n=ncur then x:=nisog:y:=disog: else x:=eT[n][11]:y:=eT[n][12]: fi; ism:=0: for i to nops(x) do w:=x[i]: if w=n then next;fi; z:=y[i]: if w=ncur then u:=nisog:v:=disog: else u:=eT[w][11]:v:=eT[w][12]: fi; if z=0 then isl:=1:ism:=1:member(n,u,'k'); if v[k]<>0 then ism:=1:y:=uae(-abs(v[k]),y,k):fi; else for j to nops(u) do t:=u[j]:s:=v[j]: if s=0 or t=n or t=w then next;fi; k:='k':q:=-abs(z*s):member(t,x,k); if y[k]=0 or q>y[k] then y:=uae(q,y,k):isl:=1:ism:=1:g:='g': if t=ncur then e:=nisog:f:=disog: else e:=eT[t][11]:f:=eT[t][12]: fi; member(n,e,g);eT[t]:=iyog(t,uae(q,f,g)); fi; od; fi; od; if ism=1 then iyog(n,y);fi; od; od; fi; if nops(nisog)=0 and iee=1 then print(`There are no curves Q-isogenous to `,cur); fi; uflag:=0:vflag:=0:store(STorE): upis();apq();nisog,disog end: Isom:=proc() global iee; iee:=1;isom(args);iee:=0;NULL end: isom:=proc() local a,b,d,f,hex,i,ia1,ia2,ia3,ie,iK,iK1,iK2,iQ,ist,j,k,n1,n2,q,rr,sc4,sc6,sd; global iee,isom_list,r,s,t,taU,u; if nargs=0 or nargs>2 then print(qts);RETURN();fi; n1:=args[1]; if nargs=2 then n2:=args[2];else n2:=stac;fi; if not type([n1,n2],list(posint)) or n1>nops(ft1) or n2>nops(ft1) then print(qts);RETURN(); fi; ist:=stac:ie:=iee:iee:=0; go(n1);ia1:=a1;ia2:=a2;ia3:=a3;a:=sex_(c4); b:=sex_(c6);d:=sex_(DD); j:=sex_(a^3/d); if type(K_,integer) and issqrfree(K_) then iK1:=K_; #elif K_=0 and type([a1,a2,a3,a4,a6],list(rational)) then iK1:=1; elif K_=T then iK1:=T; else iK1:=0; fi; if n1<>n2 then go(n2);sc4:=sex_(c4);sc6:=sex_(c6); sd:=sex_(DD); if sex_(j-sc4^3/sd)<>0 then go(ist);iee:=ie: fprint(`The curves are not even twists of each other --- their j's are \ unequal.`); RETURN(false); fi; if type(K_,integer) and issqrfree(K_) then iK2:=K_; #elif K_=0 and type([a1,a2,a3,a4,a6],list(rational)) then iK2:=1; elif K_=T then iK2:=T; else iK1:=0; fi; iK:={iK1,iK2}; if iK={T} or iK={1,T} then iK:=T; elif member(0,iK) or member(T,iK) then iK:=0; elif iK1=1 or iK2=iK1 then iK:=iK2; elif iK2=1 then iK:=iK1;iK2:=iK;go(cur,iK);n2:=stac; else iK:=0; fi; else iK:=iK1:sc4:=a;sc6:=b;sd:=d; fi; i:=sex_(d/sd);f:=1;isom_list:=[]; q:=`The curves are twists, but are not isomorphic over `; if iK=T then iQ:=Q(T);readlib(iroot); s:=normal(i);t:=factors(denom(s));s:=factors(numer(s)); rr:=iroot(s[1],12);hex:=false; if rr^12=s[1] then j:=iroot(t[1],12); if j^12=t[1] then rr:=rr/j;hex:=true;s:=s[2];t:=t[2];fi; fi; for j to nops(s) while hex do if modp(s[j][2],12)=0 then rr:=rr*s[j][1]^(s[j][2]/12); else hex:=false; fi; od; for j to nops(t) while hex do if nops(t[j])=0 then break; elif modp(t[j][2],12)=0 then rr:=rr/t[j][1]^(t[j][2]/12); else hex:=false; fi; od; if not hex then go(ist);iee:=ie: lprint(cat(q,iQ));RETURN(false); else rr:=[rr,-rr]; fi; elif iK<>0 then if iK<>1 and type (i,rational) then f:=1+w_;fi; r:=roots(X^12-i*f^12); if iK=1 then iQ:=Q;else iQ:=Q(sqrt(iK));fi; if nops(r)=0 then go(ist);iee:=ie: fprint(cat(`Ratio of DD's not a 12-th power in `,iQ)); fprint(cat(q,iQ));RETURN([]); else rr:=[]; for i to nops(r) do rr:=[op(rr),r[i][1]]; od; fi; elif type(i,algnum) then rr:=[solve(X^12-i)]; iQ:=`the indicated field.`; else go(ist);iee:=ie: lprint(`Isom is only programmed for E defined over Q or Q(sqrt(m)`); lprint(`or Q(T), i.e., when the field label K_ is a nonzero`); lprint(`squarefree integer or T. For general fields,`); lprint(`see Proposition 4.4.2 in the Elliptic Curve Handbook.`); RETURN(); fi; d:=NULL;taU:=NULL; for i to nops(rr) do u:=rr[i]/f; if sex_(a-sc4*u^4)<>0 or sex_(b-sc6*u^6)<>0 then next; fi; s:=simplify((u*a1-ia1)/2);r:=simplify((u^2*a2-ia2+s*ia1+s^2)/3); t:=(u^3*a3-ia3-r*ia1)/2; j:=sex_([r,s,t,u]);isom_list:=[op(isom_list),j];k:=length(j); if d=NULL or k E2`;fi; lprint(cat(`The `,q,` defined over `,iQ)); if n1<>n2 then if n2=ist then lprint(`where E1 has stack number`=n1,` and E2 is the present curve.`); else lprint(`where E1, E2 have stack numbers`=n1,n2); fi; fi; lprint(`are listed in isom_list. The one of minimal Maple length is`); print(` taU`=taU); fi; fi; isom_list end: isoX:=proc() local c_r; global yoX, f, g0, h0, a1c, a3c, b2c, c, a1b, a3b, b2b, a1o, a3o, b2o, ci, cpow, XI, etA, X, Y, f1, g1, h1, v; : #(internal)called by tri (by isog no longer): when we know the equations of the #isogeny from curve #nb to #ne and #nc=the twist nb*a, calculate the isog. #from #nc to #oc=ne*a by "twisting" these equations c_r:=proc(x) local y;y:=simplify(x);sign(x)*abs(x)^(1/3);end: yoX:=a1,a2,a3,a4,a6;f:=vaj[1]:g0:=vaj[2]:h0:=vaj[3]: bc(op(et2[nc]));a1c:=a1:a3c:=a3:b2c:=b2: if jay=0 then c:=c6: # elif jay=1728 then c:=c4: else c:=c6/c4: fi; bc(op(et2[nb]));a1b:=a1:a3b:=a3:b2b:=b2: if jay=0 then c:=c_r(c/c6): elif jay=1728 then c:=simplify((c/c4)^(1/2)): else c:=c*c4/c6: fi; bc(op(et2[oc])):a1o:=a1:a3o:=a3:b2o:=b2: if jay=0 then ci:=c6: elif jay=1728 then ci:=c4: else ci:=c6/c4: fi; bc(op(et2[ne])):cpow:=c^(degree(f,X)):XI:='XI':etA:='etA': if jay=0 then ci:=c_r(ci/c6): elif jay=1728 then ci:=simplify((ci/c4)^(1/2)): else ci:=ci*c4/c6: fi; X:=XI/c-b2b/12:Y:=etA/c^(3/2)-(a1b*X+a3b)/2:f1:=cpow*f: g1:=(g0+b2/12)*ci-b2o/12: h1:=(h0+(a1*g0+a3)/2)*ci^(3/2)-(a1o*g1+a3o)/2: X:='X':Y:='Y':XI:=X+b2c/12:etA:=Y+(a1c*X+a3c)/2:f1:=sort(collect(f1,X)): bc(yoX);v:=normal(expand([f1,g1,h1])) end: isprim := proc (n) local btor, nr, p, r; #options system,`Copyright 1993 by Waterloo Maple Software`,remember; #The point is to remove "remember" from the options list. #None of the subprocedures called have the option remember. if n < 2 then false elif member(n,p_up_to_100) then true #isprime/w replaced by p_up_to_100; funny they use "has" --- #according to my tests "member(n,`isprime/w`)" is faster. elif igcd(2305567963945518424753102147331756070,n) <> 1 then false elif n < 10201 then true elif igcd(849696948923341811053233990918734996592606258664893273661154542\ 634220389327076939090906947730950913750978691711866802886149933382509768238672\ 298373796296306675767413112673657893644078815718696989373063311306647862044862\ 494925732402262739543736363903875260816675866125595683463069722044751229884822\ 222855006268378634251996022599630131594564447006472069662175047724452891592786\ 7113,n) <> 1 then false elif n < 1018081 then true else nr := igcd(408410100000,n-1); nr := igcd(nr^5,n-1); r := iquo(n-1,nr); btor := modp(power(2,r),n); if `isprime/cyclotest`(n,btor,2,r) = false or irem(nr,3) = 0 and `isprime/cyclotest`(n,btor,3,r) = false or irem(nr,5) = 0 and `isprime/cyclotest`(n,btor,5,r) = false or irem(nr,7) = 0 and `isprime/cyclotest`(n,btor,7,r) = false then RETURN(false) fi; for p from 3 while (numtheory[jacobi])(p^2-4,n) <> -1 do od; evalb(`isprime/TraceModQF`(p,n+1,n) = [2, p]) fi end: isq:=proc(f) global n; local d; #internal - used by hens to check if f is a p-adic square, also whether Hensel #works for one of w,x,y in w^2=bh*x^4+Ah*x^2*y^2+bbh*y^4 if f=0 then RETURN(1);fi; n:=vam(f): if modp(n,2)=0 and ((p=2 and modp(zz,8)=1) or (p>2 and L(zz,p)=1)) then RETURN(1); fi; if xh=0 then RETURN(0);fi; d:=4*bh*xh^3+2*Ah*xh: if d<>0 and n>2*vam(d) then RETURN(2):fi; d:=2*Ah*xh^2+4*bbh: if d<>0 and n>2*vam(d) then RETURN(3):fi; 0 end: isug:=proc(x,y) local i,n; #internal - printout for Isog; also used by Dat. n:=nops(x): if n=0 then if e7<0 then lprint(`isogenous curves`.qD,` Isog`); else print(et1[ncir],` possesses no rational isogenies`); fi; else print(` isogenous curve|degree of isogeny`); for i to n do print(cat(et1[x[i]],`|`,y[i])); od; fi end: itog:=proc() local fc_,p_; global iwog, RR1, iew, iee, auw, ftw, qR, cyc, i, gen, z, hen, q, di, dg, cyi, d, a, a1, a2, a3, a4, a6, ix, x, u, t, w, vis, x7, y7, p, sol, som, nom, trou, fc, lc, s, nni, iuog; ; #internal - called 3 times by isog iwog:=iwog union {ncur}: if qR<1 then if NN<1 then RR1:=RR:iew:=iee:iee:=0:auw:=atw:tor();iee:=iew: if auw=1 and atw<>1 and member(-3,disog) then ftw:=1:fi; fi; qR:=1: if NN>1 then cyc:=2: for i in [PP] do if i[3]=NN then cyc:=1:gen:=i:break; elif i[3]=NN/2 then gen:=i: fi; od; if cyc=2 then z:=mult(NN/4,gen): for i in [PP] do if i[3]=2 and i<>z then hen:=i:break; fi; od; fi; q:=`cyclic`;z:=gen: if cyc=2 then q:=cat(`bi-`,q):z:=z,` and `,hen:fi; fprint(`the torsion group is `,q,` generated by `,z); fprint(`We now calculate isogenies:`); for di from 2 to NN do if modp(NN,di)=0 then dg:=mult(NN/di,gen): if cyc=2 then if modp(di,2)=1 then next:fi; jac([op(eadd(dg,hen)),ilcm(2,NN/(2*igcd(NN/2,NN/di)))]); if di=2 then next:fi; fi; jac(dg); fi; od; else fprint(cur,` has no points of finite order.`); fi; cyi:=1: for i to nops(nisog) do d:=abs(disog[i]): if isprim(d) then cyi:=cyi*d:fi; od; if member(jay,ljay,'j') then a:=qfr(lC[j]*c4*c6,2):a1:=0:a2:=0:a3:=0:a4:=-27*nC4[j]*(a^2): a6:=-54*nC6[j]*(a^3):velu2(dI[j],0); else if member(3*cyi,lcyi) or ftw=1 then for ix in roots(al3) do x:=ix[1]:u:=x*(x*(4*x+b2)+2*b4)+b6: if type(u,square) then next;fi; t:=x*(6*x+b2)+b4:w:=u+x*t:vis:=X-x:x7:=X+t/vis+u/vis^2: y7:=Y-((u+t*vis/2)*(2*Y+a1*X+a3+a1*vis)-u*vis*a1/2)/vis^3: velu2(3,-1); od; fi; q:=` which we now calculate.`: for p in [5,7,13] do if not member(p*cyi,lcyi) then next;fi; x:='x':som:=map(z->z[1],roots(fr.p)): nom:=nops(som); for i to nops(disog) do if abs(disog[i])=p then nom:=nom-1:fi; od; if nom=0 then next;fi; if nom<0 then print(QT,` itog - nom<0`);trou:=[op(trou),ITOG_NOMn,cur,p]: RETURN(); fi; if modp(NN,p)=0 then if nom=1 then next; else nom:=nom-1: fi; fi; if nom>2 or (nom=2 and p>5) then print(QT,` itog - nom`);trou:=[op(trou),ITOG-NOM,cur,p,nom]: RETURN(); elif nom=1 then fprint(cur,` has one isogeny of degree `,p,q); else fprint(cur,` has two isogenies of degree 5,`,q); fi; lker(p);p_:=p; for fc_ in lz7 while nom>0 do fc:=collect(fc_,X):lc:=coeff(fc,X,lkd4): for i to 3 do sy.i:=coeff(fc,X,lkd4-i)/lc:od; s:=sy1^2-2*sy2:t:=6*s-b2*sy1+lkd4*b4: w:=10*(-sy1^3+3*sy1*sy2-3*sy3)+2*b2*s-3*b4*sy1+lkd4*b6: nni:=nops(nisog):velu2(p,2);fc:=1:p:=p_; if nops(nisog)>nni then nom:=nom-1:fi; od; od; fi; fi; iuog:={op(nisog)} minus iwog end: iyog:=proc(m,a) global disog, eT; ; #internal - called twice by isog if m=ncur then disog:=a; else eT[m]:=uae(a,eT[m],12): fi end: jac:=proc() global ape, iee, jargs, diso, eur; ; ape:=iee:iee:=-1:jargs:=args:diso:=jargs[3]:velu(args);eur:=cur:ein(pc); iee:=ape: if nops(jargs)=2 then jargs:=[op(jargs),diso]:fi; fprint(cur,`/<`,jargs,`> = `,eur);apq() end: lat:=proc() local d, i, j, Le, Ls, n, x; Le:=1,args[1]:d:=args[2]: if NN>2 then for j to NN-1 do if modp(PP[j][3],2)=0 then x:=4*(PP[j][1]-d): if x<>0 then Le:=Le,qfr(x,2): fi; fi; od; fi; for j to MM do x:=RR[j][1]-d: Le:=Le,qfr(x*denom(x)^2,2): od; Le:={Le}; do n:=nops(Le);Ls:=NULL; for i to n-1 do if Le[i]=1 then next;fi; for j from i+1 to n do Ls:=Ls,Le[i]*Le[j]/igcd(Le[i],Le[j])^2; od; od; Le:=Le union {Ls}; if nops(Le)=n then break;fi; od; Le end: Lin:=proc() local j, l; global iee, TP; iee:=1:lin(args);iee:=0: if adet>0 then if adet=infinity then print(`Result inconclusive: the points may or may not be indep. mod torsion.`); RETURN(); elif nops(lsp)>1 then print(`The points are independent --- no non-trivial comb. of them \ is a torsion point.`); if adet<>1 then lprint(`Grammian height-pairing det`=adet);fi; else print(lsp,` is a torsion-free point.`); if adet<>1 then lprint(`with canon. height`=adet);fi; fi; RETURN(): elif adet<0 or K_<>1 then RETURN(); fi; if nops(lsp)=1 then print(op(lsp),` is a torsion point.`);RETURN();fi; lprint(`The points satisfy the relation`);j:='j': if nops(TP)=0 then TP:=O:l:=NULL: else l:=` a torsion point.` fi; print(cat(seq(cat(num.j,lsp[j],`+`),j=1..(i-1)),cat(-den,lsp[i]))=TP,l): NULL end: lin:=proc() local a, aqe, b, B, c, d, h, hn, j, l, llc, lld, lle, l1, l2, m, n, Q, x; global adet, den, i, iee, k, lsp, pcn, TP, trou; aqe:=iee: if lin_flag='lin_flag' then adet:=-1;n:=nargs;lsp:=args: if n=0 or not type(args[1],list) then print(`args for lin must be a sequence or list of one or more points`); print(`in the format [x,y]`,qtr);RETURN(-1); fi; if not tor_able then print(`The procedure lin needs to know the torsion subgroup.`); tor_fail();RETURN(-1); fi; do_tor(); if n=1 then if lsp=[] then lsp:=[[]];fi; if lsp=[[]] then adet:=0;RETURN(0);fi; if type(lsp[1],list) then n:=nops(lsp);lsp:=op(lsp);fi; fi; a:=lsp;lsp:=[lsp]; for h in lsp do if check_pt(h) then print(`args not defined over present ground field`,qtr); RETURN(-1); fi; od; if n=1 then if on(a) then if member(a,ouP) then adet:=0;RETURN(0); elif bouu then iee:=0;adet:=2*ht_(a);iee:=aqe; else adet:=1; fi; RETURN(adet); else print(a,noc,qtr);adet:=-1;RETURN(-1); fi; fi; a:=NULL; for x in lsp do if x=[] or member(x,ouP) or member(x,[PP]) then lsp:=[x];adet:=0:RETURN(0); fi; if not type(x,list) or nops(x)<>2 then print(qts);adet:=-1;RETURN(-1); fi; x:=sex_(x); if not on(x) then print(x,noc,qtr);adet:=-1;RETURN(-1); fi; a:=a,x: od; lsp:=[a]; else lsp:=args;n:=nops(lsp);a:=op(lsp); if n=1 then if a=[] or member(a,ouP) then adet:=0; elif K_=1 then adet:=2*ht_(a); else adet:=1; fi; RETURN(adet); fi; fi; adet:=1: if bouv then with(combinat,subsets);B:=subsets({seq(i,i=1..n)});B[nextvalue]();l:={PP}; hn:=[];llc:={[]} union {op(ouP)};lld:={}; while not B[finished] do a:=[]:b:=B[nextvalue](); for j in b do a:=eadd(a,lsp[j]);od; if member(a,llc) then eprint(`The points are dependent (mod torsion): sum_{i in b} P_i = a, where`); if a=[] then j:=NULL; elif member(a,ouP) then j:=` is a torsion point`; else j:=` = a different lin. comb. of the P_i + a torsion point`; fi; eprint(` b`=b,`and a`=a,j);adet:=0;RETURN(0); else llc:=llc union {a};lld:=lld union {a}; for c in ouP do d:=eadd(a,c);llc:=llc union {d}; if half(d)<>[] then hn:=[op(hn),[a,b,c,d,q]]:fi; od; fi; od; if nops(hn)=0 then RETURN(1);fi; lle:={}; for l1 in llc do for l2 in lld do lle:=lle union {sub(l1,l2)}; od; od; for h in hn do for Q in h[5] do if member(Q,lle) then fprint(`The points are dependent (mod torsion):`); eprint(`sum_{i in b} P_i + c = 2*Q, where`); eprint(` b`=h[2],` c`=h[3],` (a torsion point) and`); eprint(` Q`=Q,` is a lin comb. of the points + torsion`); adet:=0;RETURN(0); fi; od; od; adet:=infinity;RETURN(adet); else B:=array(1..n,1..n):iee:=0:pcn:=1: for i to n do h.i:=ht_(lsp[i]):B[i,i]:=2*h.i:od; for i from 2 to n do for j to i-1 do B[i,j]:=ht_(eadd(lsp[i],lsp[j]))-h.i-h.j:B[j,i]:=B[i,j]: od; od; for i to n do x:=B[i,i]:adet:=adet*x: if abs(x)0 then iee:=aqe:pcn:=0:RETURN(adet): fi: for j from 2 to i-1 do for k to j-1 do x:=B[k,j]:B[k,j]:=0: for l from j+1 to i do B[k,l]:=B[k,l]-x*B[j,l]: od; od; od; for j to i-1 do num.j:=loo(B[j,i]);den.j:=k; od; den:=ilcm(den.(1..(i-1))): TP:=mult(-den,lsp[i]): for j to i-1 do num.j:=(den*num.j)/den.j:TP:=eadd(TP,mult(num.j,lsp[j])):od; iee:=aqe:pcn:=0: if nops(TP)>0 then if not member([TP[1],TP[2]],ouP) then print(QT,` lin`);trou:=[op(trou),LIN,cur,lsp]:RETURN(); else TP:=neg(TP); fi; fi; adet end: lin_:=proc() global adet, lin_flag; lin_flag:=0;lin(args):lin_flag:='lin_flag';adet end: lker:=proc(d8) global lz7, i, z7, j, lf, lkd4, lz6, lg, q, d, trou, k, hex; ; #(internal) make list lz7 of all factors of the poly. al.d8 such that the roots #of the factor are the x-coord.'s of the nonzero points of an isogeny of #degree d8 whose transf. equations haven't yet been stored in visog #--- called by tri and by itog lz7:=[]: if d8=2 then for i in [PP] do if i[3]=2 then z7:=X-i[1]: for j to nops(visog) do if nops(visog[j])>0 and visog[j][1]=z7 then next;fi; lz7:=[op(lz7),z7]: od; fi; od; else div(d8);lf:=factors(al.d8)[2]: z7:=1:lkd4:=iquo(d8,2):lz6:=[]:lg:=[]: for i in lf do q:=i[1]:d:=degree(q): if d>lkd4 then next; elif d=lkd4 then lz6:=[op(lz6),q]: else lg:=[op(lg),q]:z7:=z7*q: fi; od; z7:=collect(z7,X): if degree(z7,X)=lkd4 then lz6:=[op(lz6),z7]:z7:=1:fi; if nops(lz6)>2 or (nops(lz6)=2 and d8>5) then print(QT,` lker`);trou:=[op(trou),LKER,ncur,d8]:RETURN(); fi; if z7<>1 then eprint(`Note al5 has 4 linear factors!`); member(X-mult(2,absc(-lg[1][2]))[1],lg,'j'); if j=2 then k:=lg[3]*lg[4]: elif j=3 then k:=lg[2]*lg[4]: else k:=lg[2]*lg[3]: fi; lz6:=[op(lz6),lg[1]*lg[j],k]: fi; for z7 in lz6 do hex:=0: for i in visog while hex=0 do if has(i,X) and expand(z7)=expand(i[1]) then hex:=1:fi; od; if hex=0 then lz7:=[op(lz7),z7]:fi; od; fi end: loga:=proc(x); log(max(1,abs(x))) end: loo:=proc(x) global k,top; for k while abs(x-trunc(x*k+sign(x)*.5)/k)>.001/k do:od: top:=trunc(x*k+sign(x)*.5) end: malp:=proc(nn,zm,w2,w4) local n,xm,ym,u,v,z; xm:=NULL:ym:=NULL:n:=nn: if n=0 or zm=[] then RETURN([]):fi; u:=modp(zm[1],p):v:=modp(zm[2],p): while u<>NULL do if modp(n,2)=1 then if xm=NULL then xm:=u:ym:=v: else z:=eedp([xm,ym],[u,v],w2,w4); if z<>[] then xm:=z[1]:ym:=z[2]:else xm:=NULL:fi: fi; fi; if n=1 then break;fi; n:=trunc(n/2):z:=eedp([u,v],[u,v],w2,w4); if z<>[] then u:=z[1]:v:=z[2]:else u:=NULL:fi: od; if xm=NULL then ym:=NULL:fi: [xm,ym] end: maz:=proc() local i,n,ra,uy;global trou,ue; #internal - calculate simplified version of Mazur's upper bound for the #rank when E has a point of order p and reduction at p is good or mult. if NN=1 then RETURN(-1);fi; ra:=-1;ue:=2*nops([Kod])-nops(badmp)-1: if NNeven and (member(2,badmp) or modp(Nc,2)=1) then ra:=ue;fi; if NNeven and modp(Nc,4)>0 then n:=0: if NN=2 then if type(PP[1],integer) then n:=1:fi; else for i to nops([PP]) do if PP[i][3]=2 and type(PP[i][1],integer) then n:=n+1: fi; od; fi; if n=0 then eprint(`In Mazur's terminology, ker(2) is divisible by mu2 but not by Z/2Z`); else ra:=uaa(2)-u2: fi; if ra<0 then print(QT,` maz`);trou:=[op(trou),MAZ,cur,2]:fi; fi; for uy from 3 by 2 to 7 do if modp(NN,uy)=0 and modp(Nc,uy^2)>0 then if ra<0 then ra:=uaa(uy): else ra:=min(ra,uaa(uy)): fi; if ra<0 then print(QT,` maz-1`);trou:=[op(trou),MAZ1,cur,uy]:fi; fi; od; ra end: Mest:=proc() local nL,z:global iee; if bouv then aM();RETURN();fi; iee:=1:mest(args):iee:=0:NULL end: mest:=proc() global eai, eae, x, mstr, mstR, nL, mL, lamb, ebu, r4, rc, trou; local d,i: #upper bound for rank by Mestre's method #(analog of Weil's exact formula, assuming R.H. etc. for L) if bouv then aM();RETURN();fi; for i in nisog do if i=ncur then next:fi; if member(et1[i],ft1,'j') then d:=` from stack data of isogenous curve `; if ft[j][4]>eai then eai:=ft[j][4]:eae:=copy(ft[j][5]); fprint(`Trace of Frobenius table copied`); fprint(d,ft1[j]); fi; x:=ft[j][20]: if x[2]>mstr then mstr:=x[2]:mstR:=x[3]: fprint(`Borrowed result Mest(`.mstr.`)`=mstR); fprint(d,ft1[j]); fi; fi; od; if nargs=0 then if mstr=0 then fprint(`Mest was not yet called for `.cur):RETURN(): else nL:=50: fi; else nL:=args: fi; x:='x':d:=Digits:Digits:=10:mL:=trunc(evalf(nL+.5)): if mL<51 then i:=max(1,trunc((mL-16)/5)):mL:=5*(i+4): elif mL<101 then i:=1+trunc((mL+9)/10):mL:=10*(i-1): elif mL<1001 then i:=7+trunc((mL+24)/25):mL:=25*(i-7): elif mL<5001 then i:=37+trunc((mL+99)/100):mL:=100*(i-37): fi; if nL<>mL then fprint(`Parameter adjusted to `.mL);fi; if mL>mstr then mstr:=mL: else eprint(`Mest previously calculated with parameter `.mstr); eprint(`giving the upper bound `,mstR);RETURN(mstR): fi; lamb:=ln(mL): if mL<5001 then ebu:=evalf(ln(Nc)/2)-I_F[i]: else ebu:=evalf(ln(Nc/(2*Pi))/2)-mi(lamb): fi; ebo(mL):mstR:=evalf(ebu*Pi^2/(4*lamb)): eprint(`Mestre's u.b. for rank is `,mstR);x:=trunc(mstR); if r4<>-1 and r4<=x then if rc>-1 then fprint(`Rank was determined earlier to be `,r4); elif r4ln(1000); for smaller l #the value (with ln(2Pi) added on) is looked up in the list I_F z:=evalf(Int(odly(x/l)/(exp(x)-1)-exp(-x)/x,x=0..l))+Ei(-l) end: Minf:=proc(); `Apecs is not pointing to a quadratic field --- use Qfin(m) first.` end: minf:=proc(); `Call qfin(m) first.` end: mss:=proc() local deff,h1,ie,j,qu,rr,s,sfleg; global cd,ci,co,co_lust,ds1,ds2,f_,g,h,h1_,iee,lquar,sfds,\ t,T,t0,u,u0,u2,x0,y,y0; #2nd 2-descent on ds1*X^4+Ah*X^2+ds2=Y^2 ds1:=bh;ds2:=bbh;sfleg:=3;rr:=RR; co_lust:=ds1*T^4+Ah*T^2+ds2; fprint(`2nd 2-descent on `); eprint(co_lust=Y^2,` :---`); q_form(ds1,Ah,ds2);u0:=u1;y0:=y1; if essqr(u0) then mss_();RETURN(4);fi; s:=factorset((Ah^2-4*Bh)*denom(u0)); #Putting minus {1} in above line makes no diff to sfds even when args of #factorset = 1 or -1; however we need minus {-1} in Release 3. if nops(s)>15 then fprint(`Have to abandon 2-nd descent on this quartic because set of`); fprint(`divisors > 2^15. Sorry.`); RETURN(-1); fi; sfds:=map(proc(y) local i;product(y[i],i=1..nops(y)) end,powerset(s)); if ds2<0 then deff:=sign(u0);else deff:=0;fi; if ds1<0 and deff<0 then # fprint(`Homog. space with coeff.'s`,ds1,Ah,ds2,` not sol.`); # fprint(` by 2nd 2-descent at infinity.`); RETURN(1); elif ds1<0 or deff>0 then h:=sfds; else h1:=map(proc(j) -j end,sfds); if deff=-1 then h:=h1;else h:=h1 union sfds;fi; fi; h1_:={};cd:=ilcm(denom(u0),denom(y0))^2;u2:=u0*cd; for j in h do if nr_test(j,-j*ds1)=1 and nr_test(j*u2,-j*ds2*u2)=1 then h1_:=h1_ union {j}; fi; od; if nops(h1_)=0 then sfleg:=2; else t:='t';f_:=cd*(u0*t^2-2*y0*t+Ah+ds1*u0); fi; for j in h1_ do if j=1 then t:=(T^2+ds1)/2/T; g:=collect(normal(4*T^2*f_),T); else q_form(j,0,-j*ds1);t:=u1; if essqr(j*f_) then u:=f_/(cd*(t^2-ds1));essqr(ds1*u^2+Ah*u+ds2);y0:=esqrt; essqr(u);mss_();RETURN(4); fi; t0:=u1;x0:=y1: t:=(t0*T^2-2*x0*T+j*t0)/(T^2-j); g:=collect(normal(j*(T^2-j)^2*f_),T); fi; h:=ilcm(op(map(denom,[coeffs(g,T)])))^2;g:=g*h;ci:=NULL; for ie from 0 to 4 do ci:=ci,coeff(g,T,4-ie);od; h:=ofad(igcd(ci),2)^2;g:=g/h;co:=NULL; for ie from 1 to 5 do co:=co,ci[ie]/h;od; # fprint(`We seek point on quartic`); # eprint(g=Y^2); ie:=iee;iee:=0;lquar:=2;rr:=RR:qu:=quar(co);lquar:=lquar_;iee:=ie: if type(qu,list) then fprint(`Found point on a descendant quartic.`); # if qu=[] then # fprint(`Found point at infinity on the quartic;`); # elif qu=[0] then # fprint(`Leading coeff. S1 os quartic = 0 --- called ein.`); # else # fprint(`Found point`,qu,` on the quartic;`); # fi; nowRR(rr,cur);RETURN(4) elif qu=FAIL then # fprint(`The quartic is loc. sol. but no points found.`); sfleg:=0;co_lust:=co_lust,g; else # eprint(`Quartic not loc. sol. at p`=P__); fi; od; if sfleg>0 then fprint(`Conclusion: 2nd 2-descent proves that there is no rational point.`); fprint(); fi; sfleg end: mss_:=proc() local r,um,x,y; fprint(`Found the point`,[esqrt,y0]); if hu=1 then um:=um1/2:x:=q1+bh*((um*esqrt)^2):y:=bh*esqrt*y0*(um^3): else um:=(um1*um2)/4:x:=q1+((um*y0/esqrt)^2): y:=bh*esqrt*y0*(um^3)*(bbh/(esqrt^4*bh)-1): fi; r:=RR;ap1_([x,y-(a1*x+a3)/2]) end: Mult:=proc() local n,z; n:=args[1]: if nargs=2 then z:=args[2]:else z:=[args[2],args[3]]:fi; if not on(z) and not hastype(z,float) then print(z,noc);RETURN();fi; if sT=4 and z=zs then print(`The singular point`=zs,` is not allowed as an arg.`,qts); RETURN(); fi; print(cat([n],z)=factor(simplify(mult(args)))); NULL end: mult:=proc() global zm, xm, ym, nN, Nz; local n,u,v; if nargs=3 then zm:=[args[2],args[3]]: elif nargs=2 and type(args[2],list) then zm:=args[2]: else print(`Improper args in mult`,qtr);RETURN(); fi; n:=args[1]:xm:=NULL:ym:=NULL: if n=0 or zm=[] then RETURN([]):fi: u:=zm[1]:v:=zm[2]:nN:=n: if pcn=0 then if not on(zm) and not hastype(zm,float) then print(zm,noc);RETURN(); elif sT=4 and z=zs then print(`The singular point`=zs,` is not allowed as an arg.`,qts); RETURN(); fi; fi; if nops(zm)=3 then n:=mods(n,zm[3]):fi; if n<0 then n:=-n:v:=-v-a1*u-a3:fi; while n>=1 and u<>NULL do if modp(n,2)=1 then if xm=NULL then xm:=u:ym:=v: else eadd(xm,ym,u,v);xm:=x3:ym:=y3: fi; fi; n:=trunc(n/2):eadd(u,v,u,v);u:=x3:v:=y3: od; Nz:=[xm,ym]: if nops(zm)=3 and xm<>NULL then Nz:=[op(Nz),zm[3]/igcd(zm[3],nN)]: fi; Nz end: Neg:=proc() global Edd; ; Edd:=0:neg(args); if Edd=0 then lprint(`-`,pz,` = `,nz);fi; NULL end: neg:=proc() global pz, nz, Edd; ; if nargs>1 then pz:=[args[1..nargs]]: else pz:=args: fi; if pz=[] then nz:=[]:RETURN([]):fi; if not on(pz) and not hastype(pz,float) then Edd:=1:print(pz,noc);RETURN(); elif sT=4 and pz=zs then print(`The singular point`=zs,` is not allowed as an arg.`,qtr); Edd:=1:RETURN(); fi; nz:=[pz[1],sex_(-pz[2]-a1*pz[1]-a3),op(pz[3..nops(pz)])] end: newp:=proc() global p, j0, q, b, tb2, tb4,e,es,e2,g,d2,nnn,h; local a,f,i,j,o,s,N,r,t,n; p:=args:j0:=modp(jay,p): if modp(j0-1728,p)=0 then if modp(p,4)=3 then RETURN(0): else RETURN(mods(binomial((p-1)/2,(p-1)/4)*power(-c4/48,(p-1)/4),p)); fi; elif j0=0 then if modp(p,6)=5 then RETURN(0): else RETURN(mods(binomial((p-1)/2,(p-1)/3)*power(-c6/864,(p-1)/6),p)): fi; elif member(jay,CMj,j) and L(CMd[j],p)=-1 then RETURN(0): fi; a:=NN;e:=trunc(evalf(2*sqrt(p)));es:=ceil(sqrt(e+1));e2:=2*e; q:=iquo(p,2):d2:=1;f:=a; for i from -MM while f<=e2 do if i<0 then h:=RR[-i][1]:else h:=i:fi; o:=oux/4 mod p: if o=0 then a:=ilcm(2,a):f:=ilcm(2,f):next:fi: tb2:=b2*o mod p;tb4:=b4*o^2 mod p; g:=[h*o,o^2] mod p:nnn:=malp(p-e,g,tb2,tb4): r:=bsGs(); if mods(power(o,q),p)=1 then if type(r,list) then RETURN(p+1-r[1]):fi; a:=ilcm(r,a): else if type(r,list) then RETURN(r[1]-p-1):fi; fi; f:=ilcm(r,f); od; mods(p+1-a*modp((2*p+2)/a,f/a),f) end: NN_pair:=proc() local rr, z; global NNeven, ouP, QQ, RR, sT, uR; if sT=4 then RETURN();fi; if NN>=0 then ouP:=[seq([z[1],z[2]],z=[PP])]:sT:=modp(nops(ouP)+1,2); NNeven:=evalb(sT=0); uR:=uR union {op(ouP)};rr:=RR;RR:=[]; for z in rr do uR:=uR union avR(z);RR:=[op(RR),Rz];od; QQ:=[op({op(QQ)} minus uR)]; fi end: nowRR:=proc() local r;global MM; r:=args[1];MM:=nops(RR); if iee>0 and nops(RR)>0 and not r=RR then if nargs=1 then print(`Now RR`=RR); else print(`Now RR of`,args[2]=RR); fi; fi end: nr3:=proc(x) global y; ; if modp(x,p)=0 then y:=3-L(x/p,p): else y:=2-L(x,p): fi end: nre:=proc(x) global y; ; if modp(x,2)=0 then y:=1+modp(x,16)/2: else y:=modp(x,8): fi end: nrs:=proc(ay,bz) global nr, p; if ay>0 or bz>0 then if ay<>1 and bz<>1 then nr:=nr2[nre(ay),nre(bz)]: for p in SSo while nr=1 do if modp(p,4)=1 then nr:=nrp[nr3(ay),nr3(bz)]: else nr:=nrq[nr3(ay),nr3(bz)]: fi; od; else nr:=1: fi; else nr:=-1: fi; nr end: nr_test:=proc() local i;global nr,p,s1,s2; for i to 2 do s.i:=qfr(args[i],2);od; if s1<0 and s2<0 then nr:=-1; else nr:=nr2[nre(s1),nre(s2)]; for p in factorset(s1*s2) minus {1,2} while nr=1 do if modp(p,4)=1 then nr:=nrp[nr3(s1),nr3(s2)]; else nr:=nrq[nr3(s1),nr3(s2)]; fi; od; fi; nr end: nxtprm:= proc (n) local i,j; if n < 2 then 2 else if irem(n,2) = 0 then j:=n+1 else j:=n+2 fi; for i from j by 2 while not isprim(i) do od; i fi end: oan:=proc() local a; global eT; if nops(RKI)>0 then a:=[RKI,RKJ,dolt]; else a:=dill; fi; eT[ncur]:=[Ncl,[Kod],NN,[PP],qR,RR,rc,r4,r1,a,nisog,disog,visog,reg,\ omegA,shaf,SFE,tuQ,[sB*qqR,sBb],tW,[cP,cplist],copy(misc)] end: odly:=proc(x); #Odlyzko's function used in mest (1-x)*cos(Pi*x)+sin(Pi*x)/Pi end: ofac:=proc() global zz, nP, ntab, p;local n; n:=nargs:zz:=abs(args[1]): if zz<2 then nP:=0:ntab:=[];RETURN();fi; if n=1 then ntab:=afac(zz)[2]: else for p in badp do vam(zz); if v>0 then ntab:=[op(ntab),[p,v]]:fi; od; fi; nP:=nops(ntab) end: ofad:=proc(x,n) global um, ue;local z; ofac(x);um:=1: for z to nP do ue[z]:=trunc(ntab[z][2]/n);um:=um*(ntab[z][1]^ue[z]): od; um:=expand(um) end: oin:=proc() local a, b, d, d1, e, hex, i, j, k, nergs, not_old_cur, p, z; global allP, aln, a5, a6, badp, bec, bic, cur, curi, ebi, ergs, et1, et2, etj, ft, ft1, ft2, ft3, iee, K_, mstr, MM, ncur, new_bec, new_bic, NN, omegA, ouP, oRR,p_at_inf_c, p_at_inf_w, PP, psn, q1, QQ, raF, RR, seekn, sT, stac, s1R,uR, uV, Weq, x, yo, yR, zst; ergs:=args:nergs:=nargs;a5:=a6;curi:={};psn:=stac;not_old_cur:=true; if nergs=2 then K_:=args[2]:ergs:=args[1]: # if ergs=[a1,a2,a3,a4,a6] then ergs:=cur;fi; elif nergs=6 then K_:=args[6]:ergs:=ergs[1..5]: # if [ergs]=[a1,a2,a3,a4,a6] then ergs:=cur:nergs:=2;fi; #causes trouble when 2 gcubs in a row, first with RR.. elif nergs<>1 and nergs<>5 then print(qts):RETURN(2): fi; if nergs<3 then if type(ergs,list) then if nops(ergs)=5 then ergs:=op(ergs): elif nops(ergs)=6 then K_:=ergs[6]:ergs:=op(1..5,ergs): else print(qts):RETURN(2): fi; elif type(ergs,posint) and ergs<=nops(et1) then ergs:=op(et2[ergs]): elif ergs=cur then if nops(ft1)=0 then print(qts);RETURN(2); else ergs:=a.(1..5); fi; elif member(ergs,et1,'i') then ergs:=op(et2[i]): else print(qts):RETURN(2): fi; fi; if type(stac,posint) then if [ergs]=yo and K_=oK_ then not_old_cur:=false;fi; #RETURN(1);fi; can lead to ERROR in seq. rk<--quar<--con9<--con0... if sT=4 then i:=[zs,natsing,nats,raF]; else if bouv then i:=[[allP,omegA],bec,bic,boc,raF,aln,[al.(5..aln)],sT,[QQ,RR],p_at_inf_c,\ p_at_inf_w]; if sT=0 or sT=1 then i:=[op(i),NN,[PP],uR,[sB,sBb],uV,s1R]; else i:=[op(i),[sB,sBb]]; fi; else oan(); i:=[ncur,[aln,[al.(5..aln)]],Dn,eai,copy(eae),[Ner],badmp,[badme,\ badt,bade,Ncl,SFEm,crem_dat],tt,Dt,yR,uR,eyr,bec,bic,[QQ,q1,RKI,RKJ,q_1],\ Rflag,raC,raQ,[raF,mstr,mstR,lowB,uppB],uppB,p_at_inf_c,p_at_inf_w]: #i[21]=uppB doesn't need to be remembered --- available if Rflag>0 and NNeven then i:=[op(i),L1,L2,L3,L4,AA1,BB1,ic,AA2,BB2,dol,dolt,um1,R3,um2,r3,SSo,big_flag]: fi; fi; fi; ft:=uae(i,ft,stac); fi; if type(K_,integer) and K_<>1 and issqrfree(K_) then qfin(K_);bc(ergs,1); if not type(yo,list(polynom(rational,w_))) then hex:=0; for i to 5 while hex=0 do a:=a.i; if nops(a)>2 then hex:=1; elif nops(a)=2 then if type(a,`*`) then if has(a,RootOf) then a_.i:=op(1,a)*w_; else hex:=1; fi; elif type(a,`+`) then b:=op(2,a); if nops(b)>2 then hex:=1: elif nops(b)=2 then a_.i:=op(1,a)+op(1,b)*w_; elif has(b,RootOf) then a_.i:=op(1,a)+w_; else hex:=1; fi; else hex:=1; fi; else if has(a,RootOf) then a_.i:=w_; else a_.i:=a; fi; fi; od; if hex=0 then bc(a_.(1..5),1);fi; if not type(yo,list(polynom(rational,w_))) then K_:=0; fprint(`Not all a.i are in `.Q(sqrt(m_))); fprint(`so field label K_ reset to 0. See Menu(ell);.`); fi; fi; else bc(ergs,1); fi; ergs:=a.(1..5); #Yes, bc 3 lines up is necessary when K_=1: mss calls quar, quar #finds rat. point on quartic which is birat. equiv. with original E; but #quar calls con0 which calls weier which can leave a1,.. with temp. #values, not those of E. zst:='zst':cur:='cur': if K_=1 then if DD=0 then K_:=0;hex:=[ergs,0]; else hex:=yo; fi; else hex:=[ergs,K_]; fi; if not member(hex,ft2,'stac') then if not storE then if nops(ft1)>FT1 then ft1:=ft1[1..FT1];ft2:=ft2[1..FT1]; fi; if nops(et1)>Ncur1 then et1:=et1[1..Ncur1]:et2:=et2[1..Ncur1]:etj:=etj[1..Ncur1]: fi; fi; stac:=nops(ft2)+1; fi; if psn>0 then ft3:=[seq(z[1..5],z=ft2)]; else ft3:=[]; fi; if not_old_cur then allP:={}:bec:=new_bec;new_bec:=[];bic:=new_bic;new_bic:=[]; raF:=[];aln:=4;oRR:=RR;p_at_inf_c:=[];p_at_inf_w:=[]; uR:={};ebi:=1:q1:=[]:mstr:=0;yR:=[];QQ:=[];seekn:=0; fi; if K_=1 and member(yo,et2,'ncur') then i:=eT[ncur];NN:=i[3];PP:=op(i[4]);RR:=i[6]; NN_pair(); if not member(yo,ft3) then RETURN(0);fi; else ncur:='ncur':badp:=0:PP:=NULL;RR:=[];ouP:=[];omegA:=[]: fi; oin_(); if K_=m_ then if not type(yo,list(polynom(integer,w_))) then cld_();fprint(`Denominators of the Weier. coeff.'s have been clear\ ed by substituting`); fprint(`X|-->X/uV^2, Y|-->Y/uV^3 in the Weier. equation, where uV\ `=uV); fprint(`so that now the Weierstrass coefficients are`); bc(a.(1..5),1);eprint(a.(1..5)); QQ:=map(z->[z[1]*uV^2,z[2]*uV^3],QQ); PP:=op(map(z->[z[1]*uV^2,z[2]*uV^3,z[3]],[PP])); if nops(bic)>0 then bic:=subs(X=X/uV^2,Y=Y/uV^3,bic); Weq:=Y^2+a1*X*Y+a3*Y=X^3+a2*X^2+a4*X+a6; bic:=map(z->[op(simplify(z[1..4],{Weq},[Y,X])),uV^2*z[5],z[6],\ uV^3*z[7],op(z[8..nops(z)])],bic); fi; curi:={};oin_(); fi; else if K_=T then if not type(yo,list(ratpoly(rational,T))) or not has(yo,T) then fprint(`The Weier. coeff.'s must generate Q(T) to have K_ = T.`); fprint(`Reseting K_ to 0.`);K_:=0; fi; fi; if K_=T then d:=1;d1:=1; for i to 5 do if a.i<>0 then e:=primpart(denom(a.i),T); if e=1 then d1:=d1*denom(a.i); else d:=d*e; fi; fi; od; if d*d1<>1 then eprint(\ `In order to qualify for field label K_ = T, E must be defined over Z[T].`); fprint(`Must clear denom.'s of Weier. coeff.'s.`); uV:=1: if d<>1 then d:=factors(d)[2]: for i to nops(d) do p:=d[i][1]: for j to 5 do while rem(denom(a.j),p,T)=0 do for k to 5 do a.k:=normal(p^k*a.k):od; a5:=a5*p:uV:=uV*p: od; od; od; d1:=1: for i to 5 do a.i:=normal(a.i):d1:=d1*denom(a.i);od; fi; d:=afac(d1)[2]: for i to nops(d) do p:=d[i][1]: for j to 5 do while modp(denom(content(a.j,T)),p)=0 do for k to 5 do a.k:=p^k*a.k:od; a5:=a5*p:uV:=uV*p: od; od; od; a6:=a5:bc(a.(1..5),1); fprint(`Now the coeff.'s are`,yo); fprint(`and x,y coordinates refer to these new coeff.'s. The trans.`); fprint(`to original coord.'s is [x,y]|->[x/uV^2,y/uV^3] where`); eprint(` uV`=uV); QQ:=map(z->[z[1]*uV^2,z[2]*uV^3],QQ); PP:=op(map(z->[z[1]*uV^2,z[2]*uV^3,z[3]],[PP])); if nops(bic)>0 then bic:=subs(X=X/uV^2,Y=Y/uV^3,bic); Weq:=Y^2+a1*X*Y+a3*Y=X^3+a2*X^2+a4*X+a6; bic:=map(z->[op(simplify(z[1..4],{Weq},[Y,X])),uV^2*z[5],z[6],\ uV^3*z[7],op(z[8..nops(z)])],bic); fi; curi:={};oin_(); fi; if sT=-1 then #the following is done since seekT expects s1R and sT to be there, without #summoning tor or torT; sR=0 or 1 at this point means that we're #returning to a previous curve in the stack x:=X;PP:=NULL;s1R:=0; for z in {solve(prac,X)} do if type(z,polynom(rational,T)) then s1R:=normal(z,expanded); p:=[s1R,-normal((a1*s1R+a3)/2,expanded),2];PP:=PP,p; fi; od; NN:=0;NN_pair();NN:=-1; fi; fi; fi; 0; end: oin_:=proc() local curi3, curn, f, hex, i, j, k, n, p, q, u; global curi, MM, ncur, NN, PP, QQ, uR; if member(yo,ft3,curn) then curi:={2};fi; if member(yo,et2,curi3) then curi:=curi union {3}; fi; if sex_(DD)=0 then curi:={};fi; if member(2,curi) then hex:=nops(ft3); while member(yo,eval(ft3[curn..hex]),'i') do curn:=i+curn;i:=curn-1;f:=ft[i]; if nops(ft2[i])=6 then k:=ft2[i][6];q:=f[9][1];j:=f[8]; if j=0 or j=1 then p:=f[13];u:=f[14];else p:=[];u:={};fi; if nops(f)>12 then n:=f[12];else n:=-1;fi; else k:=1;u:=f[12];q:=f[16][1]; fi; uR:=uR union {oin2(k,u)};QQ:=[op({op(QQ)} union {oin2(k,q)})]; if k<>1 then oin1(k,n,p,q,f[9][2]);fi; od; fi; if member(3,curi) then j:=eT[curi3];oin1(1,j[3],j[4],[],j[6]); fi; MM:=nops(RR);NN_pair() end: oin1:=proc(k,n,p,q,r) local e, s, t, z; global iee, MM, NN, PP, QQ, RR; if n>0 and NN=-1 then if (K_=T and k=T) or (K_=m_ and k=m_) then PP:=op(p);NN:=n; elif tor_able and (K_=k or K_=1) then PP:=oin2(k,p);NN:=nops([PP])+1 fi; fi; if NN<1 and nops([PP])>0 then s:={oin2(k,p)} minus {PP}; if nops(s)>0 then e:=iee;iee:=0;subg(PP,op(s));iee:=e;PP:=op(GG); fi; #in this case NN remains -1 fi; s:={oin2(k,q)};t:=[oin2(k,r)]; if nops(t)>nops(RR) then QQ:=[op({op(RR)} union {op(QQ)} union s)];RR:=t;MM:=nops(t); else QQ:=[op({op(t)} union {op(QQ)} union s)]; fi; end: oin2:=proc(k,p) local q, z; if k=K_ then q:=op(p); else q:=NULL; for z in p do if type(z,list(rational)) then q:=q,z;fi; od; fi; q end: Om:=proc(); if om()<>NULL then lprint(`Lattice basis omega[1],omega[2]:`); print('omega[1]'=omega[1],'omega[2]'=omega[2]); lprint(`Factor Omega in Birch, Swinnerton-Dyer conjecture:`); print('Omega'=Omega); if Omega<>omega[1] then lprint(`i.e, 2*omega[1] since DD>0.`); else lprint(`i.e., omega[1] since DD<0.`); fi; fi; NULL end: om:=proc() local a, b, e1, pi; global omega, omegA, Omega; #For E/|R, calc. lattice basis omega[1], omega[2] (omega is a table) #where omega[1] is real and positive, and omega[2]/omega[1] has positive #imag. part and real part=0 (when DD>0) or -1/2 (when DD<0). In the former #case Omega (the factor in the Birch, Swinnerton-Dyer conj.) = 2*omega[1], #and in the latter case Omega=omega[1]. Set omegA:=[[Omega,omega,Digits],..] #where ,.. represents the sequence (usually empty) of p-adic periods #calc. by omp, each a list of length>3. Thus type(omegA,list), nops(omegA)>0 #and nops(omegA[1])=3 signal that om has been called previously. #Return value = Omega. if not type(yo,list(numeric)) or sT=4 then lprint(`For om, Weier. coeff.'s a1,..,a6 must be real and DD<>0`); lprint(qtr); RETURN(); fi; if type(omegA,list) and nops(omegA)>0 and nops(omegA[1])=3 then if omegA[1][3]>=Digits then RETURN(Omega); else omegA:=omegA[2..nops(omegA)]; fi; else omegA:=[]; fi; raf();pi:=evalf(Pi); if DD>0 then Omega:=2*pi/agM(sqrt(raF[3]-raF[1]),sqrt(raF[3]-raF[2])): omega[1]:=Omega/2; omega[2]:=I*pi/agM(sqrt(raF[3]-raF[1]),sqrt(raF[2]-raF[1])): else e1:=raF[1]:a:=3*e1+b2/4;b:=sqrt((e1*(6*e1+b2)+b4)/2): Omega:=2*pi/agM(2*sqrt(b),sqrt(2*b+a)): omega[1]:=Omega; omega[2]:=-Omega/2+I*pi/agM(2*sqrt(b),sqrt(2*b-a)); fi; omegA:=[[Omega,copy(omega),Digits],op(omegA)]:Omega end: Omp:=proc() global iee; ; iee:=1:omp(args);iee:=0:NULL end: omp:=proc() local d;global nmp, imp; if bouv then aM();RETURN();fi; if nargs=0 then print(`Omp expects args - see Menu(omp)`);RETURN();fi; if nargs=1 then if type(jay,integer) then fprint(`jay = `,jay,` is integral`); else nmp:=args[1]:d:=afac(denom(jay))[2]; for imp to nops(d) do ump(d[imp][1],nmp);od; fi; else ump(args); fi end: On:=proc(); print(on(args)) end: on:=proc() local x,y,z; z:=args: if z=[] then RETURN(true):fi; if nargs>1 then z:=[args]:fi; if hastype(z,float) then RETURN(true);fi; x:=z[1]:y:=z[2]: evalb(sex_(y^2+a1*x*y+a3*y-(x^3+a2*x^2+a4*x+a6))=0) end: oon:=proc() global etag; ; if (nargs<>1 and nargs<>5) or (nargs=1 and nops(args[1])<>1 and \ nops(args[1])<>5) then print(`args must be a catalog name (eg A11) or a catalog number (eg 1)`); print(`or a sequence or list of five Weierstrass coefficients`);etag:=2: fi end: orgs:=proc() local ergs;global x1, y1, x2, y2, p; ergs:=sex_([args]); if nargs<4 then if ergs[1]=[] then x1:=NULL:y1:=NULL: else x1:=ergs[1][1]:y1:=ergs[1][2]: fi; if ergs[2]=[] then x2:=NULL:y2:=NULL: else x2:=ergs[2][1]:y2:=ergs[2][2]: fi; else x1:=ergs[1]:y1:=ergs[2]:x2:=ergs[3]:y2:=ergs[4]: fi; if nargs=3 then p:=ergs[3]:fi; if nargs=5 then p:=ergs[5]:fi end: ouu:=proc() global asB, d, dmin, d2, g, h, iee, iuu, m1, M1, m2, M2, m3, M3, m4, M4, m5, M5, m6, M6, nou, nwec, RR, searchl, sB, seekn, sflag, U, x; local e, f, i, j, k, l, L, M, N, r, rr, t, u; #internal - used by bas, rk, seek, seek1 and search #search for rational points [x,y] with max{|numer(x)|,denom(x)}<=L #where L is args possibly adjusted because of the location of the real point(s) #[e_i,f_i] of order 2 and, when e:=max{e_i}>0, #also integral points with e1. L:=trunc(evalf(args)):dmin:=1:nou:=0:g:=trunc(max(op(raf()))):rr:=RR; if L<=sB and sB>1 then RETURN();fi; if bouu then do_tor(1); for d from 2 to 3 do if emod(d)=1 then eprint(cur,` has no integral points since it has no`,d,`-adic integral\ points.`); dmin:=dmin*d: fi; od; if sflag=1 and dmin>1 then RETURN():fi; if sflag<2 then if qR<5 then M:=` not in the subgroup generated by RR and ouP;`; else M:=` ;`; fi; eprint(`Searching for [x,y] on `.cur.M); fprint(`our starting list RR = `,RR); if NNeven then ouu2(L):RETURN():fi; fi; else nwec:=stac;e:=iee:iee:=0;cld();iee:=e; if uV>1 then g:=trunc(max(op(raf()))):fi; fi; M:=L:N:=L: if b6<0 then e:=min(op(map(abs,raF))):f:=trunc(e):L:=L+f:M:=min(L,1+L/e): asB:=sB+f: else asB:=sB: fi; if sflag=1 then eprint(`Integral points where |x| is between `.asB.` and `.L); M:=1: else eprint(`Rational points where x has naive exp. height between `.asB.` & `.L); if seekn>0 then eprint(`or until number of points in RR is`,seekn);fi; fi; if bouv and uV>1 and iee>0 then lprint(`Because we had to clear denominators (temporarily) by the substit\ ution`); lprint(`x=X/uV^2, y=Y/uV^3, where uV`=uV,`,`); lprint(`these limits actually refer to X, not x.`); lprint(`However any points found will be reported in the notation [x,\ y].`); if sflag=1 then lprint(`The points [x,y] for which [X,Y] is integral are stored in the list`); lprint(`searchl, while IP is the list of integral [x,y] that are known.`); lprint(`(searchl contains just one of a pair +-P, while IP contains both.)`); fi; fi; for d from dmin by dmin while d^2<=M do d2:=d^2: for iuu to 2+sign(DD) do h.iuu:=trunc(d2*raF[iuu]+.01*((-1)^iuu)):od; m1:=max(d2,asB,h1): if DD>=0 then M1:=min(L,h2):m2:=max(d2,asB,h3):M2:=L:else M1:=L:fi; m3:=max(-L,h1):M4:=min(-d2,-asB): if DD>=0 then M3:=min(M4,h2):m4:=max(-L,h3):else M3:=M4:fi; if d2>=asB then M6:=d2-1:m5:=max(-M6,h1): if DD>=0 then M5:=min(M6,h2):m6:=max(-M6,h3):else M5:=M6:fi; fi; l:=[]: for iuu by (3-sign(DD))/2 to 6 do if iuu<5 or d2>=asB then if m.iuu>M.iuu then next:fi: h:=min(abs(m.iuu),abs(M.iuu)):x:=[m.iuu,M.iuu,h]: if m.iuu<0 and M.iuu>0 then l:=[[0,M.iuu,0],[m.iuu,0,0],op(l)]: elif nops(l)=0 then l:=[x]: else if m.iuu=0 or M.iuu=0 or h<=l[1][3] then l:=[x,op(l)]: elif h>=l[nops(l)][3] then l:=[op(l),x]: else for k from 2 to nops(l) do if h=0 then t:=l[iuu][2]:u:=1: else t:=k:k:=l[iuu][2]:u:=-1: fi; for h from k by u to t do if igcd(h,d)=1 and issqr(oux) then ouw();fi; if seekn>0 and seekn<=MM then if sflag=0 then sB:=max(abs(h),d^2)-1:else sflag:=0;fi; seekn:=0;RETURN(); fi; od; od; od; if bouu and sflag<2 then if nou=0 then fprint(`No points added to RR in this range.`); elif rr=RR then fprint(`Only redundant points found in this range.`); fi; nou:=0: if g>L then d:=1:d2:=1: eprint(`We look for integral points with x between`,g+sB,` and`,g+N); fprint(); for h from g+sB to g+N do if issqr(oux) then ouw():fi; od: if nou=0 then fprint(`None found.`);fi: fi: fi: if tim>0 then clk(); fi; if sflag=0 then sB:=N:else sflag:=0:fi; if bouv then if uV>1 then e:=iee;iee:=0;go(nwec);iee:=e; if sflag=0 then sB:=N:else sflag:=0:fi; fi; if nops(searchl)=0 then fprint(`No points in this range.`); else e:=iee;map(ap1_,searchl);iee:=e; if rr=RR then fprint(`Only redundant points found in this range.`);fi; fi; fi end: ouu2:=proc() local l,K,L,M,ouu3,R,raG,S,s1,s2,s3,s4,ssB,T; global bb1, bb2, d, d2, fleg, h, hens_d, h_s, iee, L1, L2, L3, L4, Lt1, nou, nP, ntab, qb, qd, qs, Rflag, sB, seekn, sflag, s_1, z; K:=args:L:=trunc(evalf(sqrt(K))):ouu3:=iee:fleg:=1:raG:=map(x->x-q1,raF): if DD<0 then T:=0:else member(0,raG,'T'):fi; if modp(Rflag,2)=1 then l:=L1 union L3: elif Rflag<4 then q1AB(): if big_flag then Lt1:=prod(bb2,bb1,AA1,SS1,L1) minus L1; hens(L1,Lt1,AA1,BB1):L3:=Lu minus L1;l:=L1 union L3:Rflag:=1 else L2:=NULL; for l to L4 do L2:=L2,hens_(L3[l]);od; L2:=[L2]; if bb2<0 or (AA1<=0 and bb1>0) then s_1:=1;L1:=[[1,1]]; else s_1:=2;h_s:=hens_(-1);L1:=[[1,1],[-1,h_s]]; fi; l:=extend_L1(K);Rflag:=4; fi; else l:=extend_L1(K); fi; if sB=1 then s1:=min(30,abs(q1*denom(q1)));d2:=1;d:=1; for h from max(trunc(raF[1]),-s1) to s1 do if issqr(oux) then ouw();fi; od; nou:=0; fi; if ouu3=1 then ssB:=trunc(evalf(sqrt(sB)+.99)): if sflag=0 then lprint(`Search for rational points with x = q1 + (s/d)^2*b where \ q1 = `.q1.` ,`); lprint(`max(s*sqrt(|b|),d) lies between `.ssB.` and `.L); else lprint(`Integral points where x = q1 +s^2*b where q1 = `.q1.` ,`); lprint(`s*sqrt(|b|) lies between `.ssB.` and `.L); fi; if nops(l)<200 then lprint(`and b is a Hensel tested divisor of `.BB1,`: one of `,op(l)); else lprint(`and b is a Hensel tested divisor of `.BB1); fi; fi; if sflag=1 then M:=1:else M:=L:fi: for qb in l do R:=evalf(sqrt(K/abs(qb))):S:=max(1,trunc(evalf(sqrt(sB/abs(qb))+.99))): for qd to M do d2:=qd^2: if d2>sB then s1:=1:else s1:=S:fi; if T=0 then if qb<0 then s2:=-1:else s2:=R:fi; elif T=1 then if qb<0 then s2:=-1: else s2:=min(R,evalf(sqrt(raG[2]/qb))): if K>raG[3] then s3:=max(S,trunc(evalf(sqrt(raG[3]/qb)+.99))):fleg:=2: fi; fi; elif T=2 then if qb<0 then s2:=min(R,evalf(sqrt(raG[1]/qb))): else s1:=max(s1,trunc(evalf(sqrt(raG[3]/qb)+.99))):s2:=R; fi; else if qb<0 then s1:=max(s1,trunc(evalf(sqrt(raG[2]/qb))+.99)): s2:=min(R,evalf(sqrt(raG[1]/qb))): else s2:=R: fi; fi; for s4 to fleg do if s4=2 then s1:=s3:s2:=R:fi: if s2=-1 then next;fi; for qs from s1 to s2 do if igcd(qs,qd)=1 then h:=q1*d2+qb*qs^2: #The previous version had the following which replaced "to s2" in the #qs loop. I don't remember what class of examples this served, but #its no good when q1 is large and negative: # if qs>s2 and ((abs(h)>K and sign(h)=sign(qb)) or \ # (fleg=2 and s4=1)) then break;fi; if issqr(16*oux) then z:=absc(h/d2):d:=qd:ouw(): if seekn>0 and seekn<=MM then if sflag=0 then z:=h/d2-q1;sB:=max(numer(z),denom(z))-1; else sflag:=0; fi; seekn:=0;RETURN(); fi; fi; fi; od; od; fleg:=1: od; od; iee:=ouu3: if sflag=0 then sB:=K: if nou=0 then fprint(`No points added to RR in this range.`);fi; else sflag:=0: if nou=0 then fprint(`None found.`);fi; fi end: ouw:=proc() local nP,P,Q: global nou, searchl, x; x:=h/d2:P:=[x,(isqrt(16*oux)/(4*d*d2)-a1*x-a3)/2]:nP:=neg(P); if bouv and uV>1 then P:=[P[1]/uV^2,P[2]/uV^3];nP:=[nP[1]/uV^2,nP[2]/uV^3];fi; if bouu and sflag<2 then nou:=1:ap1_(P); else searchl:=[op(searchl),P]:fprint(` P`=P,` -P`=nP): fi end: Pad:=proc(a,pp,n) local s; pad(args); if p>0 then if ad=0 then print(0); else lprint(`The `.p.`-adic expansion of`,a,`to`,n,`digits`); s:=cat(p,`^`,ad[1],`*`):print(s,ad[2..nops(ad)]); fi; else print(`improper args in pad`,qtr); fi end: pad:=proc(aD,pD,nD) global p, ad, aE, aF, za, va, zz, i, j; p:=0: if type(aD,rational) and isprim(pD) and nD>0 then p:=pD: if aD=0 then ad:=0: else aE:=denom(aD):aF:=aE*aD:vam(aF);za:=zz:va:=v:ad:=[va-vam(aE)]: zz:=modp(za/zz,p^nD): for i to nD do j:=modp(zz,p);ad:=[op(ad),j]:zz:=(zz-j)/p:od; fi; ad; fi end: pra:=proc() local x; global iee, MN, r4, r3, rc, SFE, qR, uflag, niv, tim, trou; MN:=max(MM,MN): if r1>-1 and (r4=-1 or r4>r1) then r4:=r1;fi; if MN=r1 then r4:=r1:r3:=r1:rc:=7:SFE:=(-1)^modp(r4,2): if r1=0 then qR:=5:x:=iee:iee:=0;iee:=x;fi; pr4(); if uflag<>2 then uflag:=0:upis();fi; RETURN(): fi; x:=max(r3,MN): if x>r3 then r3:=x:fprint(`Lower bound for rank is `,x,` because we have that number`); fprint(`of independent point(s) of infinite order`); if MN>MM then fprint(`on an isogenous curve.`);fi; fi; if rflag<3 then if (SFE=1 and modp(r3,2)=1) or (SFE=-1 and modp(r3,2)=0) then r3:=r3+1: eprint(`Lower bound r3 for rank increased to `.r3.` by conj'd parity of rank.`) fi; fi; if r3MN or r42 then uflag:=0:upis();fi; elif r4>=0 then print(QT,` pra`);trou:=[op(trou),PRA,cur,yo]: fi end: print_PP:=proc() local q, q1, q2; NN_pair(); if iee=0 or DD=0 then RETURN();fi; q1:=`All we know at present is that the order of the torsion subgroup`; q2:=`is odd, possibly 1.`; if NN=-1 and (K_<>m_ or old_NN='old_NN') then q:=`Tor has not been called --- `; if K_=T then if sT=1 then print(q);print(q1);print(q2);RETURN(); else print(q);q:=`However, the pts. of order 2 are known, and`; fi; else print(cat(q,`torsion subgroup unknown.`));RETURN(); fi; fi; if K_=m_ and NN=-1 then RETURN();fi; if NN=0 then if nops(ouP)=0 then print(q1);print(q2); RETURN(); else q:=`Now the`; print(`The order of the torsion group is (possibly a multiple of)`,\ nops(ouP)+1); fi; elif NN=1 then print(`The torsion group is trivial.`);RETURN(); else q:=`The`; print(`The order of the torsion group is`,NN); fi; lprint(cat(q,` list of torsion points is`));print(O,PP); lprint(`The non-O point(s) are in the notation [x,y,d], d=order; they form`); lprint(`the sequence PP, while in the notat. [x,y] they form the list ouP.`); if NN=0 then lprint(`(possibly incomplete).`); lprint(`Note that NN is not the group order --- it remains 0.`); fi end: pr4:=proc() local nergs,q;global prank,trou; if iee=0 or bouv then RETURN();fi; nergs:=nargs;fprint(); if prank<>[rc,r1,r3,r4] then if rc>3 then lprint(`The rank r4 unconditionally = `,r4,` -- no conj.'s assumed,`); lprint(`so the quality index rc = 7.`); pr4qR(); elif rc=-1 then lprint(`The rank r4 has not been determined, so the quality index rc remains \ at -1.`); if r4>0 then lprint(`And the value r4`=r4,` is only an upper bound for the rank.`); if r4=r1 then lprint(`This upper bound does not depend on any conjectures (r4=r1).`); else lprint(`This upper bound depends on one or more of the standard conjectures.`); fi; elif r4=0 then print(QT.` pr4: r4=0 and rc=-1`);trou:=[op(trou),PR4]; fi; else print(`Rank r4`=r4,` with quality index rc`=rc); lprint(`Since rc<4 (see Menu(Rk)) the value of r4 depends on`); lprint(`one or more of the standard conjectures.`); pr4qR(1); fi; if rc<4 and r1>r4 then lprint(`An upper bound for the rank not depending on any conjectures i\ s r1`=r1); fi; prank:=[rc,r1,r3,r4]; fi; if nergs>0 and rc<4 then lprint(); lprint(`The options to possibly determine the rank using no conj.'s are:`); if NNeven then lprint(`--- call RkNC(args) with the integer args bigger than`,dil); fi; lprint(`--- (for E "not too big") call Crem`); if r1>0 then lprint(`--- call Seek(args) with args bigger than`,sB); fi; if rc=-1 then lprint(`Otherwise allow one or more conj.'s to be used by calling`); lprint(`--- Rk0, Rk1, or Rk, possibly with args >`,dil); lprint(`--- Mest or Bas , possibly with (appropriately large) args`); lprint(`See Menu(Rk0), etc.`); fi; fi end: pr4qR:=proc() local q,r; if r4>0 then if qR>2 then if nargs>0 and qR>5/2 then r:=`, subject to the conjectures assumed,`; else r:=NULL; fi; if qR<3 then q:=` not `; elif qR=3 then q:=possibly; elif qR=4 then q:=probably; else q:=definitely; fi; fprint(`Since qR`=qR,`Bas has determined that`,r); fprint(`RR is`,q,`a basis`); fprint(`(For qR see Menu(Bas)).`); elif MN=r4 then if MM=r4 then q:=NULL;else q:=` of an isogenous curve `;fi; fprint(`Moreover, RR`,q,`contains the right number of indep. points.`); fi; fprint(); fi end: prod:=proc(b_,b,a,s,L) local j,k;global d,Ld,Lt,m,x,y; m:='m';x:=map(proc(y) product(y[m],m=1..nops(y)) end,powerset(s)); if b_>=0 and (a>0 or b<=0) then x:=x union map(proc(j) -j end,x); fi; Lt:=L;x:=x minus Lt: while nops(x)>0 do d:=x[1];Ld:=map(proc(k) d*k/igcd(d,k)^2 end,Lt); x:=x minus Ld; if nrs(d,b_)=1 then Lt:=Lt union Ld;fi; od; Lt end: pythy := proc () local consol, i, eqn1, a, b, c , ab, bc, ac, sa, sb, sc, A, B, C, solx, soly, solz, S; option `Copyright 1993 by Waterloo Maple Software`; #adapted from isolve/pyth #the call to mysolve takes us to code of Dave Rusin in the file legen eqn1 := igcd(args[1..3]); a := args[1]/eqn1; b := args[2]/eqn1; c := args[3]/eqn1; sa := sqrprt_(a); sb := sqrprt_(b); sc := sqrprt_(c); A := a/sa; B := b/sb; C := c/sc; ab := igcd(A,B); ac := igcd(A,C); bc := igcd(B,C); A := A*bc/ab/ac; sa := sqrt(sa)*ab*ac; B := B*ac/bc/ab; sb := sqrt(sb)*bc*ab; C := C*ab/bc/ac; sc := sqrt(sc)*bc*ac; S := zip(proc (X, Y) X/Y end,mysolve(A,B,C),[sa, sb, sc]); S := map(proc (n, d) n*d end,S,ilcm(op(map(denom,S)))); consol := igcd(op(S)); S:=map(proc(X,Y) X/Y end,S,consol) end: q_form:=proc(d1,A,d2) local T,v1;global q,u,u1,v,y1; u:='u';v:='v';q:=d1*u^2+A*u*v+d2*v^2; T:=pythy(-1,d1,4*d1^2*d2-d1*A^2);v1:=2*d1*T[3]; u1:=(T[2]-A*T[3])/v1;y1:=T[1]/v1; [u1,y1] end: q1AB:=proc() local i; global AA1,AA2,bb2,BB2,BB1,bb1,big_flag,L1,L3,L4,SS,SSo,SS1,SS2,um1,um2; #called by rka and ouu2 if Rflag=0 then if nargs=0 then two():fi; BB1:=(8*q1)*(6*q1+b2)+8*b4:AA1:=12*q1+b2: ofad(igcd(BB1,AA1^2),4);BB1:=BB1/(um^4):AA1:=AA1/(um^2):um1:=um: BB2:=AA1^2-4*BB1:bb1:=qfr(BB1,2): L3:=map(i->i[1],ntab);L4:=nops(L3);SS1:={op(L3)}; L1:=lat(bb1,q1);ofad(igcd(BB2,4*(AA1^2)),4);BB2:=BB2/(um^4): AA2:=-2*AA1/(um^2):um2:=um:bb2:=qfr(BB2,2):SS2:={op(map(i->i[1],ntab))}; if L4<15 and nops(SS2)<15 then big_flag:=true;else big_flag:=false;fi; SS:=SS1 union SS2 union {2}:SSo:=SS minus {2}; fi end: qfr:=proc(x,e) global y; ofad(x,e):y:=x/(um^e) end: rac:=proc() global x, raC, raQ; ; if raC=[] then x:='x':raC:=[solve(prac)]: raQ:=map(proc(x) if type(x,rational) then x;else NULL;fi;end,raC): fi end: Raf:=proc(); raf(args); if sT=4 then if natsing=node then lprint(`E is singular -- has node with x-coord.=raF[2]=raF[3] where`); else lprint(`E is singular -- has cusp with x-coord.=the common entry in`); fi; elif type(yo,list(numeric)) then lprint(`x-coord.'s of real point(s) of order 2:`); else lprint(`x-coord.'s of points of order 2:`); fi; print(` raF`=raF) end: raf:=proc() global x, raF; x:='x': if type(yo,list(numeric)) then if raF=[] then raF:=[fsolve(prac)]:fi; if (DD>=0 and nops(raF)<>3) or (DD<0 and nops(raF)<>1) then raF:=[]:Digits:=Digits+10:raf(); fi; else raF:=simplify([solve(prac,x)]): fi; raF:=sort(raF) end: Rk:=proc() global iee,prank; iee:=1:prank:='prank';rk(args);iee:=0:NULL end: rk:=proc() local ape,drb,i,j,jrb,MM1,mrb,q,q2,q3,rb,rki,srb,stor,x; global arc,ar4,dil,epe,go2,iee,niv,prank,rc,rflag,rflG,Rflag,RKI,RKJ,r3,r4, storE,tim,tom,trou,uflag; if bouv then aM();RETURN();fi; if nargs=1 then dil:=args: if not type(dil,posint) then print(`Args in rk must be a positive integer`,qtr);rflag:=0;RETURN(); fi; dil:=max(dil,2); else dil:=30: fi; #Note dil is now used when NN is odd: seek(dil) is called in rka if qR<1 then eprint(`Before proceeding with rk we complete the isogeny network.`); isog(); fi; upis(); if rc>-1 then fprint(`Rank already determined to be `.r4); if rc=2 then fprint(`using (possibly all of) the T, B-SD, RH conjectures.`); elif rc=3 then fprint(`using the T-p conjecture, but no other conjectures.`); else fprint(`without using any of the standard conjectures.`); if r4>0 then if rc=7 then if MM=r4 then q:=NULL;else q:=` of an isogenous curve `;fi; fprint(`Moreover, RR`,q,`contains the right number of indep. points.`); fi; fprint(`To find a Mordell-Weil basis try the command Bas.`); fi; rflag:=0;RETURN(); fi; q:=not NNeven; if iee>0 then q3:=`Consider using Seek, Bas or Crem ---`; if rflag=3 then q2:=`conjecture-independent upper bound for the rank`; if r1=-1 then lprint(`r1 has its default value -1 which means that we have no`); lprint(cat(q2,`.`)); if q then lprint(`Only Crem can determine a value r1>=0 ---`); fi; else lprint(A,q2,` is r1 =`,r1); if q then lprint(q3);fi; fi; fi; if q then fprint(q3); fprint(`RkNC can do no more (2-descent involving the cubic field Q(x) where`); fprint(`[x,y] is a point of order 2 is not programmed in apecs).`); fi; fi; if q then rflag:=0;RETURN(); elif NNeven and dil<=dill then fprint(`Consider calling RkNC(args) with args >`,dill); rflag:=0;RETURN(); fi; fi; go2:=stac; if SFE=0 then ape:=iee:iee:=0:roha();iee:=ape: fi; if not member(ncur,RKJ) then if rflag<3 and iee=1 then print(` `);lprint(`The sign of the functional eqn. is SFE`=SFE); if SFE=1 then q:=even:else q:=odd:fi; lprint(`which conjecturally implies that the rank r4 is `.q): lprint(`This information will be used if necessary.`):print(` `); fi; arc:=rc:ar4:=r4:r4:=r1:r3:=MN;rc:=-1;epe:=iee; drb:={ncur} union {op(nisog)} minus RKJ;rki:=NULL;srb:=NULL; writeto(garb);iee:=0;stor:=storE:storE:=true; for jrb in drb do ein(jrb);mrb:=Rflag;Rflag:=0;Rflag:=mrb; if NNeven and nops(q_1)=3 then rki:=jrb,rki;else srb:=srb,jrb;fi; od; writeto(terminal);iee:=epe; rki:=[op(sort([rki])),op(sort([srb]))];rflG:=rflag; for i in rki do if member(i,RKJ) then next;fi; rflag:=3;rka(i);r3:=MN;pra(); if rc>-1 then break fi; od; if r4=-1 then r4:=ar4;fi; if ar4>=0 then if (r4=0) or (r4>ar4 and rc>=0) then print(QT,` rk, arc`);trou:=[op(trou),RK,cur]:rflag:=0;RETURN(); fi; r4:=min(r4,ar4);rc:=max(rc,arc); fi; x:=RKI;j:=RKJ;writeto(garb);iee:=0; for i in rki do ein(i);RKJ:=j;RKI:=x;od; writeto(terminal);iee:=epe;go(go2); #now all isog. E are in stack and with same RKI and RKJ if iee=0 then print(cus,cur=yo);fi; prank:='prank';upis();store(stor); if rc>-1 then rflag:=0;RETURN(); else rflag:=rflG; fi; fi; if rflag<3 then if r4>-1 and modp(r4+(1-SFE)/2,2)=1 then r4:=r4-1: fprint(`Parity of SFE reduces u.b. for rank by 1 to r4`=r4); if r4<0 then print(QT,` rk`);trou:=[op(trou),` r4`]:RETURN();fi; fi; pra(); if r3=r4 then rflag:=0;RETURN();fi; # if rflag=2 then pr4(1);rflag:=0;RETURN();fi; fi; if tim>0 and niv=1 then eprint(`clock started`);tom:=time():fi; if MM=0 and r4<>0 then ouu(50); if MM>0 then pra();fi; fi; # upis();pr4(1);rflag:=0;RETURN(); if rflag<2 then rb:=evalf(ln(Nc)^2):drb:=[1,3/2,2,3]:mrb:=[1,3/2,4/3,3/2,4/3]: srb:=[100,200,300,500,700]:arc:=rc:ar4:=r4; for jrb to 5 do rb:=evalf(rb*mrb[jrb]): if rb>mstr then if iee>0 then if jrb=1 then lprint(`We now assume the T and B-SD conjectures and the R.H. for L`); lprint(`and calculate Mestre's upper bound for the rank`); lprint(`using the parameter ln(Nc)^2 = `,rb); elif jrb<5 then lprint(`We redo Mestre's u.b. with parameter`); print(drb[jrb],`*ln(Nc)^2 = `,rb); else lprint(`Final try at Mestre's u.b. with parameter`); print(`4*ln(Nc)^2 = `,evalf(rb)); fi; fi; rc:=-1:uflag:=1:mest(rb);uflag:=0:rc:=2:pra(); if r3=r4 then rflag:=0;RETURN();fi; fi; MM1:=MM:ouu(srb[jrb]); if MM>MM1 then pra();fi; if r3=r4 then rflag:=0;RETURN();fi; if tim<0 then if niv=1 then niv:=0:tim:=0:fi; rc:=arc:r4:=ar4:rflag:=0; RETURN(); fi; od; if rflag=0 then if tim<0 then if niv=1 then niv:=0:tim:=0:fi; rc:=arc:r4:=ar4: RETURN(); fi; pra(); if r3=r4 then RETURN();fi; fprint(`We now determine the smallest i such that L^{i}(1)<>0:`);x:=0; for rb from r3 by 2 to r4 do ape:=iee:i:=rc:rc:=-1:fnL(rb):rc:=i:iee:=ape: if eyy<>0 then x:=1;break;fi; od; if x=0 then print(QT,` rk`);trou:=[op(trou),` r4/fnL`]:rc:=-1:RETURN();fi; r4:=rb:r3:=r4:pra(); else rc:=-1:pr4(1); fi; fi; if rflag<3 then i:=MN+modp(MN+(1-SFE)/2,2); else i:=MN; fi; if not type(r3,integer) or r33 or (rc>-1 and (rflag<2 or (rflag=2 and rc=3))) then RETURN();fi; fprint(`The search for points on homogeneous spaces is resumed.`); for i in RKI do ape:=iee;iee:=0;ein(i);iee:=ape;Rflag:=5;rkb(); if rc>3 then if stac<>go2 then go(go2);fi; break; fi; od; if stac<>go2 then go(go2);fi; upis();pr4(1);rflag:=0;NULL end: rka:=proc() local ape, aq1, brc, br1, count, c_ur, idolt, ir1, iRKI, iRKJ, jrb, point, q, q1_, Q_1, rkj, sta_c: global cn, dolt, fleg, iee, Lt, Lt1, L2, L3, L4, L22, MN, naic, niv, nuic, pcn, prank, q1, q_1, rc, rflag, Rflag, RKI, RKJ, r1, r3, R3, r4, R4, tom: #for E with ncur=i, get best estimates for r1 (locally denoted ir1) #without using any conjectures: first apply maz, which no longer contains #SFE adjustments, and when NN is even apply 2-descent ape:=iee:iee:=0:brc:=rc:br1:=r1:iRKI:=RKI:iRKJ:=RKJ:rc:=arc:r4:=ar4: ein(args):RKI:=iRKI:RKJ:=iRKJ:sta_c:=stac;c_ur:=cur: iee:=ape:cn:=ncur:rc:=brc:r1:=br1:r4:=r1:r3:=MN:dolt:=0:Rflag:=0: prank:=[rc,r1,r3,r4]: fprint(cus,E=cur,` The torsion subgroup has order `.NN): ir1:=maz():niv:=niv+1: if ir1>=0 then fprint(`Mazur's upper bound for the rank is `.ir1): if r1>-1 and r10 then eprint(`This gives the upper bound r1`=r1):fi: r4:=r1: elif ir10 then eprint(`This gives the improved upper bound r1`=r1):fi: r4:=r1: fi: if rkc() then RETURN():fi: fi: if not NNeven then RKJ:=RKJ union {ncur}: if r1=-1 or r1>0 then seek(dil): if MM>MN then MN:=MM:fprint(`Lower bound for rank improved to MN`=MN): fi: fi: RETURN(): fi: fleg:=0:q1AB():Rflag:=2:idolt:=[]:count:=1:rkj:=ncur:Q_1:=q_1: if big_flag then for q1_ to nops(Q_1) while r4=-1 or r31 then Rflag:=0:q1AB(1):Rflag:=2: fi: aq1:=q1:brc:=rc:br1:=r1: pcn:=ncur:ape:=iee:iee:=0:point:=absc(q1): ell(0,AA2,0,BB2,0):tor(): ein(yo);naic:=cur;nuic:=ncur;ell(0,AA2,0,BB2,0); iee:=ape:RKI:=iRKI:RKJ:=iRKJ:rc:=brc:r1:=br1:r4:=r1: if count>1 then fprint(): fprint(`The previous attempt at descent didn't work, so we try another:`): fi: count:=count+1: eprint(`We use the 2-isogeny E-->E':=E/{O,P} where P`=point): eprint(`that is, `,c_ur,`-->`,naic): eprint(`Initial image in Q*/Q*^2 for E`=L1): L2:=lat(bb2,0):L22:=L2: eprint(`Initial image in Q*/Q*^2 for E' is `,L2): go(sta_c):RKI:=iRKI:RKJ:=iRKJ:q1:=aq1:L2:=L22: jrb:=r4:r1:=br1:r4:=r1:rc:=brc:prank:=[rc,r1,r3,r4]:q_1:=Q_1: Lt1:=prod(bb2,bb1,AA1,SS1,L1) minus L1: Lt:=prod(bb1,bb2,AA2,SS2,L2) minus L2: R3:=nops(L1)*nops(L2):r3:=max(r3,tam (R3)): R4:=nops(L1 union Lt1)*nops(L2 union Lt): if r4<0 then r4:=sam(R4): else r4:=min(r4,sam(R4)) fi: pra(): if rkc() then rflag:=rflG:q1:=q_1[1]: eprint(`We didn't even have to resort to Hensel's lemma.`): RETURN(): fi: fprint(`The quadratic norm residue test leaves the following candidates`): fprint(`to be tested for local solvability by Hensel's lemma:`): eprint(`for E: `,Lt1):eprint(`for E': `,Lt): hens(L1,Lt1,AA1,BB1):L3:=Lu minus L1:Rflag:=3: eprint(`For E hensel gives L3`=L3):pra(): if r3=r4 then q1:=q_1[1]:RETURN():fi: hens(L2,Lt,AA2,BB2):L4:=Lu minus L2: eprint(`For E' hensel gives L4`=L4):pra(): if rkc() then rflag:=rflG:q1:=q_1[1]:RETURN():fi: Rflag:=-1:eiop():Rflag:=5: if rkc() then rflag:=rflG:q1:=q_1[1]:RETURN():fi: idolt:=[op(idolt),[dol,q1,L1,L2,L3,L4,naic]]:rkj:=rkj,nuic: od: q_1:=Q_1:q1:=q_1[1]:dolt:=idolt:RKI:=[op(RKI),ncur]:RKJ:=RKJ union {rkj}: rkb():pr4(): else fprint(`Can't do 2-descent because at least one of initial lists of`): fprint(`candidates for the Selmer group has > 2^14 members, which is`): fprint(`beyond the present capacity of Maple. Sorry. But we carry on`): fprint(`as best we can.`): jrb:=L4+nops(SS2)-1: if r1=-1 or jrb`): eprint(`unconditional upper bound on rank r1`=r1): fi: if modp(r1,2)<>(1-SFE)/2 then jrb:=jrb-1:fi: if r4=-1 or r4>jrb then r4:=jrb: if iee>0 and r43 then RETURN();fi; fi; od end: rkc:=proc() local a; global epe; if r3=r4 then RETURN(true);fi; a:= modp(MN+(1-SFE)/2,2); if r3+a=r4 then if epe>0 then lprint();epe:=epe-1; lprint(`At this point we have the conjecture-free u.b. for the rank r1`=r4); lprint(`the conj.-free l.b r3`=r3,` and the sign of the conj. functional`); lprint(`equation SFE`=SFE,`.`); if rflG=3 then if NNeven then lprint(`Thus the rank probably`=r4,` but we continue the 2-descent process `); else lprint(`Thus the rank probably`=r4,` but we continue`); fi; lprint(`in the hope of confirming this value of r4 without any conj.'s.`); lprint();RETURN(false); else lprint(`Thus the T-p (parity of rank) conj. implies that the rank r4`=r4); lprint(`If desired, the command RkNC() or RkNC(args) with args >`,max(dill,30)); lprint(`may be able to confirm this value of r4 without using any conj.'s.`); lprint(); fi; fi; if rflG<3 then RETURN(true);fi; fi; false end: Rk0:=proc() global iee, prank, rflag; iee:=1;rflag:=2;prank:='prank';rk(args);iee:=0;rflag:=0;NULL end: rk0:=proc() global rflag; #invoke rk(args) with rflag=2 which means calculate the rank r4 allowing #use of T-p conj. (parity of r4 given by sign SFE of functional equation) #but do not allow use of Mestre's upper bound for r4, nor the conj. #r4 = the least r for which L^(r)(1)<>0. rflag:=2:rk(args);rflag:=0:print(`rc,r4,r1`=rc,r4,r1);NULL end: Rk1:=proc() global iee, prank, rflag; iee:=1:rflag:=1:prank:='prank';rk(args);iee:=0:rflag:=0:NULL end: rk1:=proc() global rflag; #invoke rk(args) with rflag=1 which means calculate the rank r4 allowing #use of T-p conj. (parity of r4 given by sign SFE of functional equation) #and Mestre's upper bound for r4, but do not allow use of the conj. #r4 = the least r for which L^(r)(1)<>0 if rk reaches that point rflag:=1:rk(args);rflag:=0:print(`rc,r4,r1`=rc,r4,r1);NULL end: RkNC:=proc() global iee, prank, rflag; iee:=1:rflag:=3:prank:='prank';rk(args);iee:=0:rflag:=0:NULL end: rkNC:=proc() global rflag; #invoke rk(args) with rflag=3 which means attempt to calculate the rank r4 #without using any conj.'s rflag:=3:rk(args);rflag:=0:print(`rc,r4,r1`=rc,r4,r1);NULL end: Roha:=proc() global iee; if bouv then aM();RETURN();fi; if jay=1728 or jay=0 then lprint(`Birch, Stephens, Liverance formula for j`=jay,` gives`); else lprint(`By results of Rohrlich and Halberstadt:`); lprint(`_______________________________________________________`); lprint(` bad prime p || Kodaira type || local root number at p`); lprint(`_____________||______________||______________________________`); fi; iee:=1:roha();iee:=0; NULL end: roha:=proc() global cas2,cas3,roh, p, iR, vj, ro, trou, SFE, eT; local K,Ko,rc4,rc6,rDD,x; if bouv then aM();RETURN();fi; if jay=1728 then x:=qfr(27*c4,4): if modp(x,4)=0 then x:=-x/4:fi; roh:=sign(-x): if member(modp(x,16),{1,3,11,13}) then roh:=-roh:fi; for p in factorset(abs(x)) do if modp(p,4)=3 and modp(x,p^2)=0 and modp(x,p^3)>0 then roh:=-roh:fi; od; elif jay=0 then roH0(qfr(-54*c6,6)): else tate():roh:=SFEm:Ko:=[Kod]:cas2:=-1:cas3:=-1; for iR to nops(badp) do p:=badp[iR]:K:=Ko[iR]; if member(p,badmp) then fprint(` `,p,` || `,K,` || `,-badt[iR]);next: fi: vj:=val(jay): if vj<0 then if p=2 then vam(-c6):ro:=mods(zz,4):cas2:=0; else ro:=mods(p,4): fi; if p=3 then cas3:=0:fi; else vj:=vam(DD); if p<5 then rDD:=zz:ro:=-1:vj:=[vam(c6),vj];rc6:=zz; vj:=[vam(c4),op(vj)];rc4:=zz; fi; if p=2 then if K=II then if vj[3]=4 then if vj[1]=4 then ro:=mods(rc4,4);cas2:=1; else cas2:=2; if vj[1]=5 then ro:=1;fi; fi; elif vj[3]=6 then if vj[1]=4 then cas2:=3; if vj[2]=7 then ro:=1:fi; elif vj[1]=5 then ro:=-mods(rc4,4);cas2:=5; else ro:=mods(rc6,4);cas2:=6; fi; else cas2:=4; if modp(rc6,8)=5 or modp(rc6-5*rc4,8)=0 then ro:=1;fi; fi; elif K=III then if vj[3]=4 then if vj[1]=4 then cas2:=7; if modp(rc4*rc6,8)=3 then ro:=1;fi; else cas2:=8; if modp(rc6,8)=5 then ro:=1;fi; fi; elif vj[3]=6 then cas2:=9;x:=modp(rc4-c6/32,16); if x=7 or x=11 then ro:=1;fi; elif vj[3]=8 then cas2:=10; if modp(rc6,8)=3 or modp(2*rc4+rc6,8)=7 then ro:=1;fi; elif vj[2]=8 then cas2:=11; if abs(mods(2*rc6+rc4,8))=1 then ro:=1;fi; else cas2:=12; if modp(rc4,8)<4 then ro:=1;fi; fi; elif K=IV then if vj[1]=4 then cas2:=13;else cas2:=14;fi; elif K=`IV*` then if vj[1]=4 then cas2:=15;else cas2:=16;fi; elif K=`III*` then if vj[3]=10 then cas2:=17;x:=modp(rc4-2*rc6,64); if x=3 or x=19 then ro:=1;fi; elif vj[3]=12 then cas2:=18; if (modp(rc4,4)=1 and abs(mods(rc6,8))=1)\ or (modp(rc4,4)=3 and modp(rc6,8)<4) then ro:=1;fi; elif vj[3]=14 then cas2:=19; if (modp(rc4,4)=1 and modp(rc4*rc6,8)>4)\ or (modp(rc4,4)=3 and abs(mods(rc4*rc6,8))=1) then ro:=1;fi; elif vj[2]=11 then cas2:=20; if modp(2*rc6+rc4,8)<4 then ro:=1;fi; else cas2:=21; if modp(rc4,8)>4 then ro:=1;fi; fi; elif K=`II*` then if vj[3]=11 then cas2:=22; if modp(rc6,8)<4 then ro:=1;fi; elif vj[3]=12 then cas2:=23; if modp(rc6,4)=1 then ro:=1;fi; else cas2:=24; if modp(rc6,4)=1 then ro:=1;fi; fi; elif K=`I*0` then if vj[3]=8 then if vj[2]=6 then cas2:=25; if modp(2*rc6+rc4,16)=3 then ro:=1;fi; else cas2:=27; if vj[1]=6 then ro:=1;fi; fi; elif vj[3]=9 then cas2:=26; if modp(rc6,8)=7 or modp(2*rc6+rc4,32)=11 then ro:=1;fi; elif vj[1]=6 then cas2:=28; if modp(rc4*rc6,4)=3 then ro:=1;fi; else cas2:=29; if modp(rc6,4)=1 then ro:=1;fi; fi; elif K=`I*1` then if vj[1]=4 then cas2:=30; if modp(2*rc6+rc4,32)=23 then ro:=1;fi; else cas2:=31; if modp(2*rc4+rc6,8)=3 then ro:=1;fi; fi; elif K=`I*2` then if vj[3]=10 then cas2:=32;ro:=1; elif vj[3]=12 then cas2:=33; if vj[2]=10 then ro:=1;fi; else cas2:=34; if modp(rc4,16)=11 or modp(rc4+4*rc6,16)=3 then ro:=1;fi; fi; elif K=`I*3` then if vj[2]=6 then cas2:=35;ro:=1; elif vj[2]=10 then cas2:=36; if member(modp(rc4+4*rc6,16),{9,13}) then ro:=1;fi; else cas2:=37; if member(modp(rc4+c6/256,16),{5,9}) then ro:=1;fi; fi; elif K=`I*4` then cas2:=38; if vj[3]=14 then cas2:=39; if modp(rDD-rc6,4)=0 then ro:=1;fi; fi; elif K=`I*5` then cas2:=40; if modp(rDD,4)=3 then ro:=1;fi; else cas2:=25+vj[3]; if modp(rc6,4)=3 then ro:=1;fi; fi; elif p=3 then if K=II then if vj[3]=3 then if vj[2]=3 then cas3:=1; if member(modp(rc6,9),{4,7,8}) then ro:=1;fi; else cas3:=2; if modp(rc4+rc6,3)=0 then ro:=1;fi; fi; elif vj[3]=4 then cas3:=3;ro:=1; else cas3:=4; if modp(rc6,3)=2 then ro:=1;fi; fi; elif K=III then ro:=1; if vj[2]=3 then cas3:=5;else cas3:=6;fi; elif K=IV then if vj[3]=5 then cas3:=7; if modp(rDD-rc6,3)=0 then ro:=1;fi; elif vj[3]=6 then cas3:=8; if modp(rc4,3)=2 then ro:=1;fi; else cas3:=9; if modp(rc6,3)=2 then ro:=1;fi; fi; elif K=`IV*` then if vj[3]=9 then if vj[2]=6 then if vj[1]=4 then cas3:=10; if member(modp(rc6,9),{4,8}) then ro:=1;fi; else cas3:=11; if modp(rc6,9)<3 then ro:=1;fi; fi; else cas3:=12; if modp(rc6,3)=2 then ro:=1;fi; fi; elif vj[3]=10 then cas3:=13; if abs(mods(rc6,9))=2 then ro:=1;fi; else cas3:=14; if modp(rc6,3)=1 then ro:=1;fi; fi; elif K=`III*` then ro:=1; if vj[2]=6 then cas3:=15;else cas3:=16;fi; elif K=`II*` then if vj[3]=11 then cas3:=17; if modp(rc6,3)=1 then ro:=1;fi; elif vj[3]=12 then cas3:=18;ro:=1; else cas3:=19; if modp(rc6,3)=1 then ro:=1;fi; fi; elif K=`I*0` then if vj[1]=2 then cas3:=20;else cas3:=21;fi; fi; else x:=igcd(vj,12):ro:=1: if ((x=2 or x=6) and modp(p,4)=3) or (x=3 and modp(p,8)>4) or (\ x=4 and modp(p,3)=2) then ro:=-1: fi; fi; fi; fprint(` `,p,` || `,K,` || `,ro);roh:=roh*ro: od; fi; if jay<>1728 and jay<>0 and iee=1 then print(`The product of the above times -1 (the root number at infinity\ ) ====>`); fi; eprint(`The root number SFE, conjecturally the sign of the functional eq\ uation, `=roh); if SFE<>0 then if SFE=roh then eprint(`which agrees with the previously determined value.`); else print(QT,` roha`);trou:=[op(trou),ROHA,cur,SFE,roh]; fi; else SFE:=roh: if nops(nisog)>0 and eT[nisog[1]][17]=0 then fprint(`Since SFE is an isogeny invariant, catalog data for all isoge\ nous curves`); fprint(`is being updated.`); for x in nisog do eT[x]:=[op(eT[x][1..16]),roh,op(eT[x][18..22])]:od: fi; fi; roh end: roH0:=proc(a) global p, roh; ; p:=2:vam(a): if (modp(v,2)=0 and modp(zz,4)=3) or (v=4 and modp(zz,4)=1) then roh:=-1: else roh:=1: fi; p:=3:vam(a): if modp(v,3)=2 or (modp(v,3)=0 and abs(mods(zz,9))=2) or (modp(v,3)=0 and\ abs(mods(zz,9))=1 and modp(zz*(-1)^v,3)=2) then roh:=-roh: fi; for p in badp do if p>3 and modp(p,3)=2 then roh:=-roh:fi; od end: RR_update:=proc() global MM, QQ; QQ:=[op({op(QQ)} minus uR)]; ap0(1);MM:=nops(RR); end: rsqr:=proc(x) global y; y:=x*denom(x)^2: if issqr(y) then isqrt(y)/denom(x):else -1:fi end: sam:=proc(x) global v, r1; local y; y:=x: for v from 0 while y>=8 do y:=y/2: od; if r1=-1 then r1:=v: else r1:=min(r1,v): fi; if (SFE=-1 and modp(v,2)=0) or (SFE=1 and modp(v,2)=1) then if rflag<3 then v:=v-1:fi; if v0 then lprint(`The variable searchl is a list of all points found.`): lprint(`Note that only P, and not -P, is listed.`); else lprint(`There are no points in this range.`); fi; NULL end: search:=proc() global sB, searchl, seekn, sflag; local a, tB, z; if bouv and not type(yo,list(rational)) then aN();RETURN();fi; tB:=sB:searchl:=[]:seekn:=0: if nargs=2 then sB:=args[1]:a:=args[2]: else sB:=1:a:=args: fi: if sB<=a then sflag:=2:ouu(a):sflag:=0:fi; sB:=tB:searchl end: Seek:=proc() global iee; iee:=1;seek(args);pr4();iee:=0;NULL end: seek:=proc() local l, x; global searchl, seekn; searchl:=[]; if bouv then if K_=T then if sB=1 then eprint(`Seek already done.`); else seekT(); fi; RETURN(); elif not type(yo,list(rational)) then aN();RETURN(); fi; fi; seekn:=0; if nargs=0 then l:=sB; elif nargs=1 and type(args,posint) then l:=args; elif nargs=2 and type([args],list(posint)) then l:=args[1];seekn:=args[2]; else print(qts);RETURN(); fi; if seekn>0 and MM>=seekn then eprint(`The 2nd arg requests that seek be terminated when RR contains`); eprint(seekn,`points; but RR already contains`,MM);seekn:=0;RETURN(); fi; if K_=1 and NN>0 and rc>=0 and r4<=seekn then if iee=0 then seekn:=0;RETURN();fi; print(`Rank is believed to be `,r4); x:=readstat(` --- to continue or not type y; or n;`); if x=`n` then seekn:=0;RETURN();fi; fi; if l<=sB then fprint(cur,` has already been searched up to sB`=sB); else ouu(l); fi; seekn:=0;NULL end: Seek1:=proc() global x, iee; if bouv and not type(yo,list(rational)) then aN();RETURN();fi; if bouu and rc>=0 and r4=0 then x:=readstat(`Rank is 0 --- to continue or not type y; or n;`); if x=`n` then RETURN();fi; fi; iee:=1:seek1(args):iee:=0: if dmin=1 then lprint(`List of integral points found so far: IP=`); fi; IP end: seek1:=proc() local sC, z; global IP, l, sB, searchl, sflag; if bouv and not type(yo,list(rational)) then aN();RETURN();fi; sC:=sB:IP:=[]:sflag:=1:searchl:=[]; if nargs=1 then sB:=1:l:=args: else sB:=args[1]:l:=args[2]: fi; seek(l):sflag:=0:sB:=sC:IP:=[]; if dmin>1 then RETURN();fi; sC:={op(searchl)} union {op(map(neg,searchl))}; if bouu then sC:=sC union uR;fi; for z in sC do if type(z,list(integer)) then IP:=[op(IP),z]:fi; od; IP:=sort(IP,proc(a,b);if a[1]<=b[1] then true:else false:fi:end): nops(IP) end: Seekb:=proc() local a, b, l, x; global iee; if bouv then lprint(`Field label must be K_=1 for this procedure (hence E is defined /Q)`); lprint(qtr);RETURN(); elif nops(bec)=0 then lprint(`bec is empty (no birat. equiv. curves listed)`,qtr);RETURN(); fi; if nargs=1 and type(args,posint) then l:=args; else print(qts);RETURN(); fi; if NN>0 and rc>=0 and r4=0 then x:=readstat(`Rank is 0 --- to continue or not type y; or n;`); if x=`n` then RETURN();fi; fi; if l<=sBb then lprint(`bec has already been searched up to sBb`=sBb); else iee:=1:a:=nops(bec); if a=1 then a:=NULL;b:=curve;else b:=curves;fi; eprint(\ `Searching for rational points with naive height between `,sBb,` and `,l); eprint(`on the `,a,` birat. equiv. `,b,` in bec.`); seekb(l);iee:=0; fi; NULL end: seekb:=proc() local d,h,l,m;global e,iee,sBb; if bouv then lprint(`Field label must be K_=1 for this procedure (hence E is defined /Q)`); lprint(qtr);RETURN(); elif nops(bec)=0 then fprint(`bec is empty (no birat. equiv. curves listed)`,qtr);RETURN(); fi; if nargs=1 and type(args,posint) then l:=args; else print(qts);RETURN(); fi; if l<=sBb then RETURN();fi; e:=iee;iee:=0; if sBb=0 then seekb_(0,1);fi; for d to l do if drr and e>0 then lprint(`The point`=[h/d,j[1]],` on birat. equiv. curve in \ bec[`.i.`] makes`); print(` RR`=RR); fi; od; od end: sex_:=proc(z) if K_=m_ then qf_conv(z); else simplify(expand(z));fi end: Sfe:=proc() global iee, tim, SFE, trou; local x,y; if bouv then aM();RETURN();fi; iee:=1:sfe(args);iee:=0: if tim<0 then if niv=0 then tim:=0:fi; SFE:=sfo:RETURN(); fi; if SFE=1 then x:=`even.`:y:=`+.`: elif SFE=-1 then x:=`odd.`:y:=`-.`: else RETURN(): fi; lprint(`Under the Taniyama conj., the sign of the functional equation is`,y); lprint(`Hence under the Birch, Swinnerton-Dyer conjecture, the parity`); lprint(`of the rank is`,x); if nops(badmp)=nops(badp) then lprint(`This is opposite to the parity of the number of split multiplica\ tive`);lprint(`primes, as it should be for a semi-stable curve.`); elif sfo<>0 then if SFE=sfo then lprint(`This agrees with the earlier determination of SFE.`); else print(QT,` sfe`);trou:=[op(trou),sFE,cur,SFE,sfo];RETURN(); fi; fi; if pcn=0 and rc>=0 then lprint(`This agrees with the value `.r4.` of the rank already found.`); fi; if sfo=0 and nops(nisog)>0 then lprint(`Since L(s) is an isogeny invariant, the catalog entries`); lprint(`of all isogenous curves are being updated with SFE.`); lprint(`Also straightforward cases of twists are updated.`); fi end: sfe:=proc() global hyy, sfo, SFE, sgi, tom, sfi, n1, eyy, ey1, tim, sn1, ey2, sn2, x, trou, eT, nc, e; local i,j; #determine the sign of functional equation analytically, compare with other #determination (if any), and transfer this information to isogenous curves if bouv then aM();RETURN();fi; if nargs=0 then hyy:=-1: else hyy:=args: fi; sfo:=SFE:SFE:=0:sgi:=-1: if tim>0 and niv=0 then eprint(`clock started`);tom:=time():fi; for sfi from 2 while SFE=0 do if sfi>2 then eprint(`It's necessary to calculate more terms in the series.`); if hyy>0 then SFE:=sfo;RETURN();fi; sgi:=-3:n1:=sn1:eyy:=ey1/2: fi; ey1:=ey(sgi,sfi): if tim<0 then if niv=0 and iee=0 then tim:=0:fi; RETURN(); fi; sn1:=n:eyy:=ey2/2:n1:=sn2:ey2:=ey(sgi-1,sfi): if tim<0 then if niv=0 then tim:=0:fi; RETURN(); fi; sn2:=n:x:=4.0*ey1/ey2: if abs(x^2-1)<.1 then SFE:=trunc(x*1.1):fi; od; if pcn=0 and rc>=0 and (1-SFE)/2<>modp(r4,2) then print(QT,` sfe`);trou:=[op(trou),sFE,cur,r4,SFE];RETURN(); fi; if sfo=0 then for i to nops(nisog) do eT[nisog[i]]:=uae(SFE,eT[nisog[i]],17): od; elif sfo<>SFE then print(QT,` sfe2`);trou:=[op(trou),sFE2,cur,sfo,SFE]:RETURN(); else fprint(`The newly calculated "analytic" value of SFE agrees with`); fprint(`the previous value`,SFE); fi; SFE end: Sha:=proc() global iee; if bouv then aM();RETURN();fi; if shaf>0 then print(`Sha already done with result shaf = `,shaf); lprint(`(To redo, set shaf:=0: then call Sha; again.`); else iee:=1:sha(args);iee:=0:NULL; fi end: sha:=proc() local dad,hux,qRo,r4o,rco,sr; global tom, fyy, hyy, tim, x, shaf, p, q,qR,r4,rc; #determine the order of the Shafarevich-Tate group analytically if bouv then aM();RETURN();fi; if shaf >0 then RETURN();fi; if nargs>0 then if SFE=0 then roha();fi; if modp(MM-(1-SFE)/2,2)=1 then fprint(`Proposed basis RR has wrong parity --- clashes with SFE.`); RETURN(); fi; do_tor();qRo:=qR;qR:=3;r4o:=r4;r4:=MM;rco:=rc:rc:=max(rc,0); sha();qR:=qRo:r4:=r4o:rc:=rco; if qR<3 and shaf>0 then lprint(` shaf`=shaf,` is reset to shaf = 0 because of the uncertainty`); lprint(`that RR is a basis --- call Bas`);shaf:=0; fi; RETURN(); fi; if qR<3 then eprint(`Must find a Mordell-Weil basis first --- use the apecs command Bas`); RETURN(); fi; shaf_tr(); if tim>0 then eprint(`clock started`);tom:=time():fi; hux:=1;dad:=4; while hux=1 do hux:=0;hyy:=-1:ey(r4,dad); if tim<0 then tim:=0:RETURN();fi; x:=evalf(eyy/fyy):shaf:=trunc(x+.1):p:=2: sr:=isqrt(shaf); if sr/50 or nops(nisog)=0 or qR<3 then RETURN();fi; for i in nisog do x:=eT[i]; if x[16]=0 or x[5]<3 or nops(x[15])=0 or nops(x[15][1])=0 then next;fi; w:=x[21]; if type(w,list) then w:=w[1];fi; z:=x[14]*x[15][1][1]*x[16]*w/(x[3]^2*fyy); z1:=trunc(z+0.5); fprint(`Borrowing the earlier result shaf[`,et1[i],`]=`,x[16]); fprint(`the conj. order of the Shaf.-Tate group of the present curve is`); eprint(z); if abs(z-z1)<.01 and issqr(z1) then shaf:=z1:fprint(`which rounds to `,z1=cat(isqrt(z1),`^2`)); else fprint(`Since this is not approx. an integer square,`); fprint(`we have not yet found a basis; we set qR=5/2.`); qR:=5/2: fi; break; od end: sik_siev:=proc() local i,n,z;global p,pp,s,RS,Ps,n_rel,p_vots,nRS,s_rel; pp:=args; if modp(NN,pp)>0 then s:=0;RS:=RR; else if NN=2 then RS:=[op(RR),PP]; else if modp(NN,pp^2)=0 then s:=2;else s:=1;fi; for Ps in PP do if Ps[3]=pp^s then RS:=[op(RR),Ps];break;fi; od; fi; fi; nRS:=nops(RS);p:=2; if modp(NN,4)=0 then for z in PP do if denom(z[1])>1 then p:=3;break;fi; od; fi; p_vots:=[];n_rel:=0;n:=0;s_rel:=vector(nRS); while nops(p_vots)0 or modp(card,pp^2)=0\ or sik_siew()=[] do p:=nextprime(p); od; p:=nextprime(p);z:=nops(p_vots); if z>n_rel then n_rel:=z;n:=0;else n:=n+1;fi; od; if n_rel[] then break;fi; od; if y<>[] then mi[i]:=1;mj:=y;z:=y; for j to pp-2 do mj:=mj,eadp(y,z);od; mj:=[mj]; for z from i+1 to nRS do w:=mulp(l,RS[z]); if w<>[] then member(w,mj,'j');mi[z]:=j;fi; od; if n_rel=0 then s_rel[i]:=copy(mi);p_vots:=[i]; else for j from i to nRS do if mi[j]>0 then if member(j,p_vots) then for k from j+1 to nRS do mi[k]:=mi[k]-mi[j]*s_rel[j][k] mod pp; od; mi[j]:=0; else for k from j to nRS do mi[k]:=mi[k]/mi[j] mod pp;od; s_rel[j]:=copy(mi);p_vots:=[op(p_vots),j];break; fi; fi; od; fi; fi; y end: sqr:=proc(a,s,pn) global S2, S1; ; #Newton's method for the square root of a mod p^n that is s mod p S2:=s:S1:=S2+pn: while S1<>S2 do S1:=S2:S2:=modp((S2+a/S2)/2,pn): od end: sqrprt_ := proc (n) local i, f, m; option `Copyright 1993 by Waterloo Maple Software`; #from isolve/pyth/sqrprt m := 1; f := afac(n); for i in op(2,f) do m := m*i[1]^(i[2]-`mod`(i[2],2)) od; m end: stec:=proc() local i,z; global AA1, AA2, allP, aln, bade, badme, badp, badmp, badt, BB2, BB1, bec, bic, big_flag, boc, crem_dat, Dn, dol, dolt, Dt, eae, eai, ebi, eyr, ic, K_, L1, L2, L3, L4, lowB, mstr, mstR, nats, natsing, Nc, Nce, Ncl, ncur, Ner, omega, omegA, Omega, p_at_inf_c, p_at_inf_w, q1, q_1, r3, R3, raC, raF, raQ, Rflag, RKI, RKJ, RR, sB, sBb, SFEm, s1R, SSo, sT, stac, tt, um1, um2, uppB, uV, yR, zs; #internal - retrieve curve data from the stack # if stac=zst then eprint(cus,cur=yo);RETURN();fi; #Looks reasonable, but there can be prob.'s with new oin, etc. stac:=zst:z:=ft[zst]: if nops(z)=4 then zs:=z[1];natsing:=z[2];nats:=z[3];raF:=z[4];sT:=4;RETURN();fi; if ft1[zst]=CUR then allP:=z[1][1]:omegA:=z[1][2]; if type(omegA,list) and nops(omegA)>0 then Omega:=omegA[1][1];omega:=copy(omegA[1][2]); fi; bec:=z[2]:bic:=z[3]:boc:=z[4]:raF:=z[5]:aln:=z[6]:sT:=z[8]: RR:=z[9][2];p_at_inf_c:=z[10];p_at_inf_w:=z[11];K_:=ft2[zst][6]; for i from 5 to aln do al.i:=z[7][i-4]:od: if sT=0 or sT=1 then #now done in oin: NN:=z[12]:PP:=op(z[13]);uR:=z[14]:QQ:=z[9][1]: # " " " oin/NN_pair: ouP:=[seq([i[1],i[2]],i=[PP])]: #But note we still want these in ft so oin can look them up sB:=z[15][1]:sBb:=z[15][2];uV:=z[16]:s1R:=z[17]:RR_update(); else sB:=z[12][1];sBb:=z[12][2]; fi: RETURN(): fi; ncur:=z[1]:aln:=z[2][1]:K_:=1: for i from 5 to aln do al.i:=z[2][2][i-4]:od: Dn:=z[3]:eai:=z[4]:eae:=copy(z[5]):Ner:=op(z[6]):badmp:=z[7]:badme:=z[8][1]: badt:=z[8][2]:bade:=z[8][3]:Ncl:=z[8][4]:Nc:=Ncl[1]:badp:=Ncl[2]: Nce:=Ncl[3]:SFEm:=z[8][5]:crem_dat:=z[8][6];tt:=z[9]:Dt:=z[10]:ebi:=1: yR:=z[11]:eyr:=z[13]:bec:=z[14]:bic:=z[15]: q1:=z[16][2]:RKI:=z[16][3]:RKJ:=z[16][4]; q_1:=z[16][5]:Rflag:=z[17]:raC:=z[18]:raQ:=z[19]:raF:=z[20][1]:mstr:=z[20][2]: mstR:=z[20][3]:lowB:=z[20][4]:uppB:=z[20][5]: p_at_inf_c:=z[22];p_at_inf_w:=z[23]; #z[21] not used -- available if Rflag>0 and NNeven then L1:=z[24]:L2:=z[25]:L3:=z[26]:L4:=z[27]:AA1:=z[28]:BB1:=z[29]:ic:=z[30]: AA2:=z[31]:BB2:=z[32]:dol:=z[33]:dolt:=z[34]:um1:=z[35]:R3:=z[36]: um2:=z[37]:r3:=z[38]:SSo:=z[39]:big_flag:=z[40]; fi; if nops(boc)>0 and not member(bac,bec) then bec:=[op(bec),bac]:bic:=[op(bic),boc]:sBb:=0; p_at_inf_c:=[op(p_at_inf_c),[]];p_at_inf_w:=[op(p_at_inf_w),[]]; fi; eon();eprint(cus,cur=ft2[stac]); end: Store:=proc(); store(); if storE then print(`Normal storage of all new curves in the catalog and the stack.`); else print(`Only the present curve (and any isogenous curves found) will be \ added`); lprint(`to the catalog and stack.`); fi; NULL end: store:=proc() global Ncur1, FT1, storE; if nargs=0 then Ncur1:=nops(et1):FT1:=nops(ft1);storE:=not storE; elif not type(storE,boolean) or storE=FAIL then store(); else storE:=args; fi end: Sub:=proc(); orgs(args); if neg(x2,y2)=NULL then RETURN();fi; Eadd([x1,y1],nz);NULL end: sub:=proc(); orgs(args); if neg(x2,y2)=NULL then RETURN();fi; eadd([x1,y1],nz) end: Subg:=proc() global x; ; subg(args); if diso=0 then RETURN();fi; x:=op(GG): if diso>Nsubg then if bouv then print(`The group has order >`,Nsubg,`possibly infinite.`); else print(`The group is infinite.`); fi; lprint(`The set GG contains the non-zero elements calc.'d so far,`); lprint(`in the notation [x,y,0].`); else print(`The group has order `.diso); if diso=2 then print(`The non-zero element is `.x); elif diso>2 then print(`The non-zero elements are `.x); fi; fi end: subg:=proc() global g, gg, ig, diso, GG; local i; #Calculate the subgroup GG generated by the point(s) contained in args. #The order of GG is denoted diso. g:=args:gg:=NULL:ig:=1: if nargs>1 then g:=[g]: elif type(g,set) then g:=[op(g)]: fi; while ig<=nops(g) do if not type(g[ig],list) then if igy[1] then eadd(x[1],x[2],y[1],y[2]);ford(x3,y3,z);sy:=sy union {Qq}: fi; od; od end: Szp:=proc(); if bouv then aM();RETURN();fi; szp();lprint(`Szpiro ratio = `,szr) end: szp:=proc() global szr; ; if bouv then aM();RETURN();fi; szr:=evalf(log(abs(DD))/log(Nc)) end: tam:=proc(x) global v; ; for v from 0 while 2^vr3 and iee>0 and rflag<3 then lprint(`Parity conjecture used to improve lower bound for the rank.`); fi; fi; v end: Tate:=proc() local x,y; if bouv then aM();RETURN();fi; tate();lprint(`Conductor Nc = `.Nc.` = `,Ncl[2],`^`,Ncl[3]); if nops([Kod])>1 then lprint(`The Kodaira types at the bad primes are`);print(Kod); if Ner<>NULL then lprint(`The corresponding Neron symbols are`);print(Ner); fi; else lprint(`The Kodaira type at the bad prime is`);print(Kod); if Ner<>NULL then lprint(`The corresponding Neron symbol is`);print(Ner); fi; fi; lprint(`Product of the local Tamagawa numbers cP = `.cP); end: tate:=proc() global p, Ncl, Nc, Kod, ep, Ner, cP, tt, v, c, x, q, q2, y, z, a, b, g, p3, m2, u,cplist; local aa,m,B,it,Nce,N; if bouv then aM();RETURN();fi; if Nc>0 and type(badmp,list) and cplist<>0 then RETURN();fi; bad();Nce:=[]: if Nc>0 then for m to nops(badp) do p:=badp[m]:Nce:=[op(Nce),vam(Nc)]:od: Ncl:=[Nc,badp,Nce]: if cplist<>0 then RETURN();fi; fi; Nc:=1:Kod:=NULL:cP:=1:cplist:=[];Nce:=[]; for m to nops(badp) do p:=badp[m]:aa:=badt[m]:ep:=bade[m]: if aa<>0 then Nc:=Nc*p:Nce:=[op(Nce),1]:Ner:=Ner,B.ep: if aa=-1 then Kod:=Kod,`I `.ep.`(non-split)`:ep:=2-modp(ep,2): else Kod:=Kod,`I `.ep.`(split)`: fi; if ep>1 then cP:=cP*ep:cplist:=[op(cplist),[p,ep]]:fi; if p=2 then ep:=2*ep:fi; tt:=igcd(tt,ep*(p-aa)): else if p=2 then tate1();tate2();Nc:=2^f:Nce:=[f]:Kod:=Kod,KK:Ner:=Ner,reN: elif p=3 then tate1();tate3();Nc:=Nc*(3^f):Nce:=[op(Nce),f]:Kod:=Kod,KK:Ner:=Ner,reN: else if c4=0 then v:=3:else vam(c4);fi; Nc:=Nc*(p^2):Nce:=[op(Nce),2]: if ep=2 then Kod:=Kod,II:Ner:=Ner,C1:c:=1: elif ep=3 then Kod:=Kod,III:Ner:=Ner,C2:c:=2: elif v>2 and ep=9 then Kod:=Kod,`III*`:Ner:=Ner,C7:c:=2: elif v>2 and ep=10 then Kod:=Kod,`II*`:Ner:=Ner,C8:c:=1: else tate1();chco(mods(-b2p/12,p),0,mods((a1p*b2p/24)-(a3p/2),p)); if ep=4 then Kod:=Kod,IV:Ner:=Ner,C3:c:=2+L((a3p^2+4*a6p)/p^2,p): else chco(0,mods(-a1p/2,p),mods(-a3p/2,p^2)); if ep=6 then Kod:=Kod,`I*0`:Ner:=Ner,C4:x:='x': c:=1+nops([msolve(x^3+(a2p/p)*x^2+(a4p/p^2)*x+(a6p/p^3),p)]): if c=4 then c:=-4:fi; else if v>2 then Kod:=Kod,`IV*`:Ner:=Ner,C6:chco(p*mods(-a2p/(3*p),p),0,0); c:=2+L((a3p^2+4*a6p)/p^4,p): else chco(p*mods((a2p*a4p-9*a6p)/(2*p*(3*a4p-(a2p^2))),p),0,0); Kod:=Kod,`I*`.(ep-6):Ner:=Ner,C5.`,`.(ep-6): q:=iquo(p,2):q2:=q+q^2:x:=a3p/p^2:y:=a4p/p^3:z:=a6p/p^4: a:=a1p/p:b:=a2p/p:g:=1:p3:=0:N:=p^(ep+5): for m2 from 2 while p3=0 do u:=modp(x^2+4*z,p): if u=0 then y:=mods(y-x*a*q,N):z:=mods((z-q2*x^2)/p,N): u:=mods(y^2-4*b*z,p): if u=0 then c:=mods(-1/(2*b),p): z:=mods(g*(c*y)^3+(z+c*(y^2)*(1+b*c))/p,N): u:=y*(3*g*(c^2)*y+(1+2*b*c)/p): x:=mods(x+a*c*y,N):g:=g*p: b:=mods(b+3*c*y*g,N):y:=mods(u,N): else p3:=1:c:=3+L(u,p): fi; else p3:=1:c:=3+L(u,p): fi; od; if c=4 and modp(ep,2)=0 then c:=-4:fi; fi; fi; fi; fi; fi; tt:=igcd(tt,c*p): if c<>1 then cP:=cP*abs(c):cplist:=[op(cplist),[p,c]]:fi; fi; od; Nc:=expand(Nc):Ncl:=[Nc,badp,Nce]:Nc end: tate1:=proc() global a1p, a2p, a3p, a4p, a6p, b2p, b4p, b6p, b8p; ; a1p:=a1:a2p:=a2:a3p:=a3:a4p:=a4:a6p:=a6: b2p:=b2:b4p:=b4:b6p:=b6:b8p:=b8 end: tate2:=proc() global KK, reN, f, c, m, v, p3, trou; local a,b,u,x,y,z,g; if modp(a1p,2)=0 then chco(modp(a4p,2),0,modp(a4p*a2p+a6p,2)); else chco(modp(a3p,2),0,modp(a3p+a4p,2)); fi; tateb(2); if p3=0 then chco(0,modp(a2p,2),2*modp(a6p/4,2)); b:=a2p/2:y:=a4p/4:z:=a6p/8: if modp(b*y+z,2)=1 then KK:=`I*0`:reN:=C4: f:=ep-4:c:=2-modp(b+y,2): if c=4 then c:=-4:fi; elif modp(b+y,2)+modp(b+z,2)>0 then if modp(y,2)=1 then chco(2,0,0); fi; a:=a1p/2:b:=a2p/2:g:=1:x:=a3p/4:y:=a4p/8:z:=a6p/16: for m from 2 while p3=0 do if modp(x,2)=0 then if modp(z,2)=0 then x:=x/2:z:=z/2: else x:=(x/2)+1:y:=y-a:z:=(z+1)/2-x: fi; if modp(y,2)=0 then if modp(z,2)=0 then y:=y/2:z:=z/2: else z:=(z+b)/2+(y/2)+g:y:=(y/2)+b+3*g:x:=x+a:b:=b+6*g: fi; g:=2*g: else v:=2*m-2:p3:=1: fi; else v:=2*m-3:p3:=1: fi; od; f:=ep-4-v:KK:=`I*`.v:reN:=C5.`,`.v: if modp(z,2)=0 then if modp(bade[1],2)=0 then c:=-4:else c:=4:fi; else c:=2: fi; else if modp(b,2)=1 then chco(2,0,0); fi; if modp(a3p/4,2)=1 then KK:=`IV*`:f:=ep-6:reN:=C6: if modp(a6p/16,2)=1 then c:=1: else c:=3: fi; else if modp(a6p/16,2)=1 then chco(0,0,4); fi; if modp(a4p,16)>0 then KK:=`III*`:f:=ep-7:reN:=C7:c:=2: elif modp(a6,64)>0 then KK:=`II*`:f:=ep-8:reN:=C8:c:=1: else print(QT,` tate2`);trou:=[op(trou),TATE2,`curve not 2-minimal`]: fi; fi; fi; fi end: tate3:=proc() global KK, f, reN, c, m, p3, v, trou; local a,b,g,N,u,x,y,z; x:=modp(b2p,3): if x>0 then u:=mods(-x*b4p,3): else u:=mods(-b6p*(b4p+1),3): fi; chco(u,0,mods(a3p+u*a1p,3)); tateb(3); if p3=0 then chco(0,mods(a1p,3),mods(a3p,9)); b:=a2p/3:y:=a4p/9:z:=a6p/27:x:=mods(-b*z+(b*y)^2-y,3): if x<>0 then KK:=`I*0`:f:=ep-4:reN:=C4: if mods(b,3)=0 and mods(y,3)=-1 and mods(z,3)=0 then c:=-4: else c:=(3-x)/2: fi; elif modp(b,3)+modp(y,3)>0 then if modp(y,3)>0 then chco(3*mods(b*y,3),0,0); fi; x:=a3p/9:y:=a4p/27:z:=a6p/81:a:=a1p/3:b:=a2p/3:g:=1:N:=3^(ep+5): for m from 2 while p3=0 do u:=mods(x^2+z,3): if u=0 then y:=mods(y-x*a,N):z:=mods((z-2*x^2)/3,N):u:=mods(y^2-b*z,3): if u=0 then u:=b^2:z:=mods(g*((b*y)^3)+(z+b*(y^2)*(1+u))/3,N): u:=mods(y*(3*g*u*y+(1+2*u)/3),N):x:=mods(x+a*b*y,N): b:=mods(b*(1+9*g*y),N):y:=u:g:=3*g: else p3:=1:v:=2*m-2: fi; else p3:=1:v:=2*m-3: fi; od; f:=ep-4-v:KK:=`I*`.v:reN:=C5.`,`.v: if u>0 then if modp(ep,2)=0 then c:=-4: else c:=4:fi; else c:=2: fi; else u:=mods(z,3): if u<>0 then chco(-3*u,0,0); fi; u:=mods(b6p/81,3): if u<>0 then KK:=`IV*`:f:=ep-6:c:=2+u:reN:=C6: else u:=mods(a3p/9,3): if u<>0 then chco(0,0,9*u); fi; if modp(a4p,81)>0 then KK:=`III*`:f:=ep-7:reN:=C7:c:=2: elif modp(a6p,729)>0 then KK:=`II*`:f:=ep-8:reN:=C8:c:=1: else print(QT,` tate3`);trou:=[op(trou),TATE3,`curve not 3-minimal`]: fi; fi; fi; fi end: tateb:=proc(z) global p2, p3, KK, c, f, reN; local y; p2:=z^2:p3:=z*p2:y:=z-1: if modp(a6p,p2)>0 then KK:=II:c:=1:f:=ep:reN:=C1: elif modp(b8p,p3)>0 then KK:=III:c:=2:f:=ep-1:reN:=C2: elif modp(b6p,p3)>0 then KK:=IV:f:=ep-2:reN:=C3: if z=2 then c:=3-2*modp(a6p/4,2): else c:=2+mods(b6p/9,3): fi; else p3:=0: fi end: Tor:=proc() global iee; if sT=4 then print(qtt);RETURN();fi; iee:=1; if NN>=0 then lprint(`Tor already done with the following results.`); print_PP(); else tor(args); fi; iee:=0;NULL end: tor:=proc() local a, b, c, d, d4, d6, hox, i, itN, j, k, l, n, nbec, qq, r, rr1, SS, tf, ubt, uE, UV, vbt, z; global alz, dol, disw, iee, L2, naic, NN, NNeven, nuic, ouP,p, PP, pto, p_at_inf_c, p_at_inf_w, QQ, RR, sy, t, tN, tt, x, y; if sT=4 then eprint(qtt);NN:=0;RETURN();fi; if NN>=0 then RETURN(NN);fi; qq:=QQ; if bouv then if K_=m_ then a:=torq();print_PP();RETURN(a); elif K_=T then torT();RETURN(); elif type([a.(1..5)],list(rational)) then if member(yo,et2,j) then k:=eT[j];rr1:=op(k[6]);NN:=k[3];PP:=op(k[4]); fi; if NN=-1 then fprint(`We find the min. Weier. equation/Z, calc. the torsion group, then use`); fprint(`trwc (see Menu(trwc)) to trans. the group back to the present curve.`); a:=stac;c:=iee;iee:=0;UV:=-V^2-a1*U*V-a3*V+U^3+a2*U^2+a4*U+a6; ein(cur);naic:=cur;nuic:=ncur;member(UV,bec,nbec);j:=nops(bec)-1; if nops(p_at_inf_c)=j then p_at_inf_c:=[op(p_at_inf_c[1..(nbec-1)]),[[0,1,0]],op(p_at_inf_c[nbec..j])]; p_at_inf_w:=[op(p_at_inf_w[1..(nbec-1)]),[[]],op(p_at_inf_w[nbec..j])]; fi; tor();k:=NULL;rr1:=NULL; for j in [PP] do k:=k,[op(trwc(j,nbec)),j[3]];od; for j in RR do rr1:=rr1,trwc(j,nbec);od; go(a);PP:=k;iee:=c; fi; NN:=nops([PP])+1;QQ:=qq;RR:=[rr1];print_PP();RETURN(NN); elif K_=0 and has([a.(1..5)],T) and type([a.(1..5)],list(ratpoly\ (rational,T))) then print(`We change the field label K_ to T.`); ell(cur,T);torT();RETURN(): else z:=NN;tor_fail(); if z=0 then print_PP();RETURN(0);fi; x:=X;r:=map(z->z[1],factors(prac)[2]);n:=0; for z in r do if degree(z,X)<>1 then next;fi; a:=[op(absc(-coeff(z,X,0)/coeff(z,X,1))),2]; if member(a,[PP]) then next; else PP:=PP,a;n:=n+1; fi; od; fprint(` `); if n>0 then fprint(`However we did notice point(s) of order 2.`); fi; if n1 then b:=[args[1..nargs]]: else b:=args: fi; if b<>NULL then if nops(b)<2 or nops(b)>3 or not on(b) or not type(b[1],rational)\ or not type(b[2],rational) then b:=NULL:fi; fi; NN:=1: if tt=1 then print_PP();RETURN(1); elif tt=0 then tt:=5040; fi; NN_pair(); if modp(DD,2)=0 then fprint(`From calculations done in Tate's algorithm we see that NN di\ vides `,tt); else tt:=gcd(tt,2*(3-modp(a6,2)-modp(a1+a2+a3+a4+a6,2)+modp(a2+a4+a6,2)\ +modp(a3+a6,2))): if tt=1 then print_PP();RETURN(1); else fprint(`From calculations done in Tate's algorithm and reduc. \ mod 2, we see that NN`); fprint(`divides `,tt); fi; fi; disw:=disog:ttw(); if qR=0 or atw<>1 then n:=nops(disw):l:=map(proc(p) abs(p) end,disw):ubt:=max (1,op(l)):vbt:=1: if n=0 then ubt:=1: elif n=1 and ubt>7 then ubt:=1: elif n=2 and ubt>9 then ubt:=5: elif n=3 then if ubt=2 then vbt:=2: elif ubt=27 then ubt:=9: fi; elif n=5 then if ubt=18 then ubt:=6: elif ubt=4 then vbt:=2: elif ubt=25 then ubt:=5: fi; elif n=7 then if ubt=16 then ubt:=8: elif ubt=8 then vbt:=0: for p to 7 do if l[p]=2 then vbt:=vbt+1: if vbt>1 then break;fi; fi; od; elif ubt<8 then vbt:=2: fi; fi; tf:=igcd(tt,ubt*vbt): if tf1 and a<5 do pto:=nxtprm(pto): if modp(DD,pto)>0 then a:=a+1:emoe(pto);tt:=igcd(tt,card): fi; od; if tt=1 then fprint(`In fact reduction mod the first few primes shows that NN = 1.`); else if tf>tt then fprint(`Reduction mod odd good primes up to`,pto,`improves this to`); fprint(` NN divides `,tt); fi; t:=tt: if b<>NULL then x:=b[1]:y:=b[2]:ess(); fi; d4:=8*b4:d6:=16*b6:x:='x':alz:=x^3+b2*x^2+d4*x+d6: if (modp(a1,2)=1 or modp(a3,2)=0) and modp(tt,2)=0 then esr(0);sy:=[PP]: fi; if NN[] then x:=q[1][1];y:=q[1][2];ess();hox:=1:tN:=t/NN:break; fi; od; od; while modp(tN,2)=0 do tN:=tN/2:t:=t/2;od; torN(3); if modp(tt,9)=0 and modp(NN,9)>0 and modp(NN,3)=0 then t:=t*9;tN:=tN*9; fi; for itN from 5 by 2 to 9 do torN(itN); od; if NN0 then RETURN();fi; div(j); for z in roots(al.j) do x:=z[1]; if issqr(4*prac) then y:=(isqrt(4*prac)-a1*x-a3)/2:ess(); fi; od; while modp(tN,j)=0 do tN:=tN/j;t:=t/j;od end: tor_fail:=proc() global NN; if iee>0 then lprint(`Tor is only programmed to calc. the torsion group for E with field`); lprint(`label K_ in the following cases (for an explan. of K_ see Menu(Ell)):`); lprint(` --- E defined over Q, K_ = 1 or 0`); lprint(` --- E " " Q(T), K_=T (the var. must be T)`); lprint(` --- E " " Q(sqrt(m_)), K_=m_ (see Menu(Qfin))`); lprint(`To indicate that Tor does not apply to this E,`); lprint(`NN is given the symbolic value 0.`); else print(`apecs not programmed to determine the torsion subgroup for this E.`); fi; NN:=0 end: Trans:=proc(r,s,t,u) global iee; iee:=1;trans(args);iee:=0;NULL end: trans:=proc(r,s,t,u) local a,b,e,nn,ou,pp,q,rr,w; global bec,bic,iee,K_,ouP,p_at_inf_c,p_at_inf_w,PP,RR,RR_supp; w:=sex_(u); if w=0 then print(` u=0`,qts);RETURN();fi; if nargs=5 then K_:=args[5]; elif K_=1 then K_:=0; fi; b:=[a1,a2,a3,a4,a6];a:=RR;e:=iee;iee:=0;nn:=NN;ou:=ouP;pp:=PP; q:=sex_([(a1+2*s)/w,(a2-s*a1+3*r-s^2)/w^2,(a3+r*a1+2*t)/w^3,(a4-s*a3+2*r*a2\ -(t+r*s)*a1+3*r^2-2*s*t)/w^4,(a6+r*a4+r^2*a2+r^3-t*a3-t^2-r*t*a1)/w^6]); bec:=[op(bec),sort(U^3+q[2]*U^2+q[4]*U+q[5]-V^2-q[1]*U*V-q[3]*V,[U,V])]; bic:=[op(bic),[X-r,u^2,Y-s*(X-r)-t,u^3,u^2*U+r,1,u^3*V+s*u^2*U+t,1,[]]]; p_at_inf_c:=[op(p_at_inf_c),[]];p_at_inf_w:=[op(p_at_inf_w),[]]; ell(q,K_); if e=1 then iee:=1;We();iee:=1;else we();fi; bec:=[op(bec),sort(sex_(U^3+b[2]*U^2+b[4]*U+b[5]-V^2-b[1]*U*V-b[3]*V),\ [U,V])]; bic:=[op(bic),[u^2*X+r,1,u^3*Y+s*u^2*X+t,1,U-r,u^2,V-s*(U-r)-t,u^3,[]]]; p_at_inf_c:=[op(p_at_inf_c),[]];p_at_inf_w:=[op(p_at_inf_w),[]]; #The reason for using trant is we don't want to invoke #tor which can in some cases clear denom.'s, hence muck up the transf. #For instance, we don't want to use trcw when K_=T. if nn>0 then PP:=op(map(trant,[pp]));ouP:=map(trant,ou); fi; if nops(a)>0 then a:=map(trant,a); if nops(RR)=0 then RR:=a;nowRR([]); else RR_supp:=[op({op(a)} minus {op(RR)})]; if nops(RR_supp)>0 then fprint(`The list RR_supp`=RR_supp,` contains points transferred from`); fprint(`RR of the previous E to the present which are not in present RR.`); fi; fi; fi; NULL end: trant:=proc(p) local h,t,z;global U,V; h:=bic[nops(bic)][5..8]; if nops(p)=3 then t:=p[3];else t:=NULL;fi; U:=p[1];V:=p[2];z:=eval([h[1]/h[2],h[3]/h[4],t]); U:='U';V:='V';z end: Tri:=proc() global iee; if bouv then aM();RETURN();fi; iee:=1:tri(args);iee:=0: if tflag=0 then if diso>2 then lprint(`The non-O points of the kernel of the isogeny are the points whose`); lprint(`X-coordinates are roots of `,vo[1],` ---`,diso-1,` points in all.`); else print(`The kernel is O`,absc(-coeff(vo[1],X,0))); fi; lprint(`The equations of the isogeny are:`); print(`X'` = vo[2]);print(`Y'` = vo[3]); elif trz='trz' then RETURN(); elif trz=[] then lprint(tz,`|--->`,O); else lprint(tz,`|--->`,trz); fi; NULL end: tri:=proc() local ape,nic,rr,stuRe; global trz, tflag, oc, tz, trI, i, iflag, x, q, trip, a5, oldn, vo, d7, diso, storE, tWc, nc, vri, iee, itr, nb, eTt, vaj, ne, x7, y7, z7, vis, pcn, l, mri, uV, Wu, vflag, nisog, disog, visog, hex, nl, j, e, n, d, k, puh, pth, pvh, trIo, oco, tzo, oldo, voo, d7o, d8o, ncurn, rvX, rvY, tfleg, X, Y, e11, e13, a, b, U, V; ; trz:='trz':tflag:=1:ape:=iee; if bouv then aM();RETURN(); elif nargs=0 then print(qts);RETURN(); elif nargs=1 then oc:=args:tflag:=0:tz:=NULL: elif nargs=2 then oc:=args[2]:tz:=args[1]: else oc:=args[nargs]:tz:=[args[1..(nargs-1)]]: fi; trI:='trI': if not type(oc,integer) then i:='i': if member(oc,et1,i) or member(oc,et2,i) then oc:=i: else print(qts);tflag:=1:RETURN(); fi; fi; if oc<1 or oc>nops(et1) or oc=ncur then print(qts);tflag:=1:RETURN(); fi; if not member(oc,nisog,trI) then if qR<0 then x:=eT[oc][1]: if type(x,list) then x:=x[1]:fi; if x=Nc then q:=` known to be `:else q:=` `;fi; else q:=` `: fi; print(`This curve is not`.q.`isogenous to`,cur,qtr);tflag:=1:RETURN(); fi; if tz=[] then trz:=[]:RETURN([]);fi; if tflag=1 and (not type(tz,list) or nops(tz)<2 or nops(tz)>3 or not \ on(tz)) then print(pno,qtr);RETURN(); fi; trip:=[PP]:a5:=a6:oldn:=ncur:vo:=visog[trI]:d7:=disog[trI]:diso:=abs(d7): if not has(vo,X) then stuRe:=storE:storE:=true: if nops(tW[2])>0 and not member(oc,tW[2]) and jay<>1728 and jay<>66^3 then tWc:=tW:nc:=ncur:vri:=[]:iee:=0: for itr to nops(tWc[1]) while vri=[] do nb:=tWc[2][itr]:eTt:=eT[nb]: if has(eTt[13],X) then ein(oc);tw(tWc[1][itr]); if member(ncur,eTt[11],'k') then vaj:=eTt[13][k]: if has(vaj,X) then ne:=ncur:vri:=isoX():fi; fi; fi; od; ein(nc);iee:=ape: if vri<>[] then fprint(`Equations for isogeny`,cur,`--->`,et1[nisog[trI]]); fprint(`have been obtained from a twist.`); fprint(isog2);fprint(isog3);print(``);vo:=vri: fi; fi; if not has(vo,X) then if iflag=1 then storE(stuRe):RETURN():fi; if diso=13 then q:=`..few minutes `: elif diso=7 then q:=`..few seconds `: else q:=`..momentary `: fi; fprint(cat(q,`delay --- equations for isogeny have to be calculated.`)); fprint(isog2);fprint(isog3);print(``); if nops(vo)>0 and d7<>-3 then x7:=X:y7:=Y:z7:=[X,Y]:vis:=1:pcn:=1: for l to nops(vo)-5 do mri:=trip[vo[l]]:eadd(z7,mri); x7:=x7+x3-mri[1]:y7:=y7+y3-mri[2]:vis:=vis*(X-mri[1]): if modp(mri[3],2)=1 then mri:=neg(mri):eadd(z7,mri); x7:=x7+x3-mri[1]:y7:=y7+y3-mri[2]: fi; od; pcn:=0:uV:=vo[l]:x7:=x7*(uV^2):Wu:=vo[l+4]:x7:=(x7-vo[l+1])/(Wu^2): z7:={Y^2+a1*X*Y+a3*Y=X^3+a2*X^2+a4*X+a6}: x7:=simplify(x7,z7,[Y,X]): y7:=(y7*(uV^3)-Wu^2*x7*vo[l+2]-vo[l+3])/Wu^3: y7:=simplify(y7,z7,[Y,X]): vo:=[vis,x7,y7]: elif isprim(diso) then lker(diso):vflag:=2:iee:=0: for vis in lz7 do velu(vis):nic:=ncur: if not member(oldn,nisog) then nisog:=[op(nisog),oldn]:disog:=[op(disog),-diso]: visog:=[op(visog),[]]: fi; ein(oldn): if nic=oc then vo:=[vis,x7,y7]:break;fi; if not member(nic,nisog) then nisog:=[op(nisog),nic]:disog:=[op(disog),d7]: visog:=[op(visog),[]]: fi; od; vflag:=0:iee:=ape: else l:=[]:hex:=0: for i to nops(nisog) do if isprim(abs(disog[i])) then l:=[op(l),[nisog[i]]]:fi; od; for i while hex=0 do nl:=[]: for j in l do e:=eT[j[i]]:n:=e[11]:d:=e[12]: for k to nops(n) do e:=n[k]: if member(e,j) or e=ncur or not isprim(abs(d[k])) then next:fi; nl:=[op(nl),[op(j),e]]: if e=oc then hex:=1:fi; od; od; l:=nl: od; nl:=[]: for j in l do if j[i]=oc then nl:=[op(nl),j]:fi; od; hex:=10^5: for j in nl do member(j[1],nisog,'s'): if has(visog[s],X) then d:=0:puh:=[1]:else d:=1:puh:=[-1]:fi; for k to i-1 do e:=eT[j[k]]:member(j[k+1],e[11],'s'): if has(e[13][s],X) then puh:=[op(puh),1]: else d:=d+abs(e[12][s]):puh:=[op(puh),-1]: fi; od; if d=0 then RETURN();fi; y:=tW[1]:z:=tW[2]:hex:=0: for j to nops(y) while hex=0 do ntc:=z[j]:otc:=eT[ntc]: if otc[5]>=0 and not member(ntc,nisog) then atw:=y[j]:disw:=otc[12]: if nargs>0 and not has(otc[13],X) then for k from j+1 to nops(y) do n:=z[k]:o:=eT[n]: if o[5]>=0 and not member(n,nisog) and has(o[13],X) then atw:=y[k]:disw:=o[12]:ntc:=n:otc:=o:hex:=1:break: fi; od; fi; fprint(cur,` is the twist by `,atw,` of `,et1[ntc],` whose isogeny network`); fprint(`is known to apecs. This may shorten the task at hand.`);break; fi; od end: Tw:=proc(a) global iee; iee:=1:tw(a);iee:=0:NULL end: tw:=proc(a) local j;global twb, l, x, y, r; if bouv or not type(a,rational) then #see Elliptic Curve Handbook, p.410 ell(a1,sex_(a*a2+(a-1)*a1^2/4),a3,sex_(a^2*a4+(a^2-1)*a1*a3/2),\ sex_(a^3*a6+(a^3-1)*a3^2/4),K_); else twb:=qfr(a*denom(a)^2,2):eprint(`twist of `.cur.` by `.twb); if twb=1 then RETURN(cus,cur); elif twb=-1 and sex_(DD)<>0 and jay=1728 then fprint(`Twisting by -1 when jay = 1728 gives the same curve.`); RETURN(cus,cur); fi; if bouu and member(twb,tW[1],j) then fprint(`This twist was already calculated:`);ein(tW[2][j]); RETURN(cus,cur); fi; l:=V^2:x:=-27*c4*(a^2):y:=-54*c6*(a^3):r:=U^3+x*U+y: fprint(`Initial Weierstrass form of `.a.`*`.cur.` is `);eprint(l=r); ein(0,0,0,x,y);cus,cur; fi end: two:=proc() local i;global q1,q_1; if NN=2 then q1:=PP[1]:q_1:=[q1]; else q_1:=[]: for i in [PP] do if i[3]=2 then q_1:=[op(q_1),i[1]];fi; od; q_1:=sort(q_1);q1:=q_1[1]; fi end: uaa:=proc(uo) local ui,uu;global u2; #internal - used by maz ui:=ue:u2:=0: for uu to nops(badmp) do if modp(badme[uu],uo)>0 then if modp(badmp[uu],uo)<>1 then ui:=ui-1:fi; if uo=2 and modp(badmp[uu],4)<>1 then u2:=1:fi; fi; od; ui end: uae:=proc() local n,p,t; n:=nargs-2:p:=args[nargs]: if p+n>nops(args[n+1]) then t:=NULL; else t:=op(args[n+1][(p+n)..nops(args[n+1])]); fi; [op(args[n+1][1..(p-1)]),args[1..n],t] end: ump:=proc() global p, nnp, hex, i, e, q, u2, pn, omegA, v1, v2, nnq, e1, hix, x, y, s, eb, w, A1, b1, r, z, wi; ; p:=args[1]:nnp:=args[2]:hex:=0: for i to nops(omegA) do e:=omegA[i]: if nops(e)<4 then next;fi; if e[1]=p then if e[2]>=nnp then fprint(p,`-adic period already calculated:`);hex:=1:nnp:=e[2]:q:=e[3]: u2:=e[4]:pn:=p^nnp: else omegA:=eau(omegA,i): fi; break; fi; od; if hex=0 then v1:=vam(c4):v2:=vam(DD); if v2<=3*v1 then print(`j integral at`,p,qtr);RETURN(); elif v1>0 then print(`Not yet programmed for additive Kodaira type I*.`);RETURN(); elif p=2 then print(`Not yet programmed for p = 2.`);RETURN(); fi; nnq:=v2-3*v1+nnp:pn:=p^nnq:e:=modp(-b2/4-2*(18*b6-b2*b4)/c4,p):e1:=e+1: hix:=0: if NNeven then x:=[PP]: for i to nops(x) do y:=x[i]: if y[3]>2 then next; elif modp(y[1]-e,p)=0 then hix:=2:e:=y[1]:break; fi; od; fi; while hix=0 do while e<>e1 do e1:=e:e:=modp((e*(e*(8*e+b2))-b6)/(2*(e*(6*e+b2)+b4)),pn): od; if e=0 then hix:=2: else e1:=vam(e)+nnp-nnq: if e1>0 then pn:=pn*p^e1:nnq:=nnq+e1: else hix:=1: fi; fi; od; if iee=1 then if hix=2 then print(`root e1`=e); else e1:=pad(e,p,nnp):s:=cat(`root e1 = `,p,`^`,e1[1],`*`): print(s,e1[2..nops(e1)]); fi; fi; eb:=12*e+b2:hix:=0:w:=-eb/4:nnq:=nnq+v2: while hix=0 do w:=sqr(e*(eb+b2)+2*b4,w,pn):A1:=modp((-eb/4+w)/4,pn): b1:=modp(w/2,pn):vam(A1-b1); if v>nnq-nnp then nnq:=nnp+v:pn:=pn*p^v: else hix:=1: fi; od; r:=modp(A1-b1,pn):x:=modp(r/4,pn): while A1<>b1 do z:=modp(b1*sqr(A1/b1,1,pn),pn):A1:=modp((A1+b1+2*z)/4,pn):b1:=z: r:=modp(A1-b1,pn):x:=x*((1+sqr(1+r/x,1,pn))/2)^2: od; u2:=modp(1/(4*A1),pn):r:=2*u2*x:wi:=r/(1+r):u2:=pad(u2,p,nnp): q:=pad((1-sqr(1-wi^2,1,pn^2))/wi,p,nnp): omegA:=[op(omegA),[p,nnp,q,u2]]: fi; fprint(`The Henniart-Mestre algorithm for`,p,`^`,nnp=p^nnp,`yields`); s:=cat(`Tate's q`=p,`^`,q[1],`*`):eprint(s,q[2..nops(q)]); s:=cat(` u^2`=p,`^`,u2[1],`*`):eprint(s,u2[2..nops(u2)]) end: unk:=proc() global X, Y, vnk, ynk, XR, YR, U_, V_; ; X:=args[1]:Y:=args[2]:vnk:=U_:ynk:=V_:X:='X':Y:='Y':XR:=X:YR:=Y: U_:=sex_(vnk):V_:=sex_(ynk):XR:='XR':YR:='YR' end: upis:=proc() local x, y, z;global eT, MN, n, qR, rc, r1, r4; #(internal) used by eon, isog, pra, rk, crem and mest: on the basis of new #value(s) for rc,r4,r1, update the corresponding values of isogenous curves. if uflag=1 then RETURN();fi; MN:=MM;x:=[rc,r4,r1]: for n in nisog do if n=ncur then next:fi; y:=eT[n][6..9]; MN:=max(MN,nops(y[1])); rc:=max(rc,y[2]); if y[3]>-1 then if r4=-1 then r4:=y[3];else r4:=min(r4,y[3]);fi; fi; if y[4]>-1 then if r1=-1 then r1:=y[4];else r1:=min(r1,y[4]);fi; fi; od; if rc>3 then r1:=r4:fi; if MN=r1 then rc:=7; if r1=0 then qR:=5; for n in nisog do #Aug.5/97 if n=ncur then next:fi; eT[n]:=uae(5,eT[n],5); od; fi; fi; if iee>0 and x<>[rc,r4,r1] then lprint(`Values of rc,r4,r1 were `,op(x)); lprint(`now updated to `,rc,r4,r1); fi; for n in nisog do if n=ncur then next:fi; y:=eT[n]:x:=y[7..9]: if x<>[rc,r4,r1] then eT[n]:=uae(rc,r4,r1,eT[n],7); fprint(`Values of rc,r4,r1 for curve `,et1[n],` were `,op(x)); fprint(`now updated to `,rc,r4,r1); #no need to update stack since rc, r4, r1 are stored only in cat fi; od; if MM0 do v:=v-1: od; od end: val:=proc(xx) global v, v1; ; if xx=0 then v:=1/2: else vam(denom(xx));v1:=v:vam(numer(xx));v:=v-v1: fi end: vam:=proc(yy) global zz, v; ; zz:=yy: for v from 0 while modp(zz,p)=0 do zz:=zz/p: od; v end: Velu:=proc() global iee; print(\ `We use formulas derived from Velus' to calculate the isogenous curve.`); iee:=1:velu(args);iee:=0:NULL end: velu:=proc() global ergs, vis, n, dp, iee, z7, diso, vflag, veg, vig, x, d, vog, z1, d2p, d3p, p2, aav, ts, t, w, x7, y7, X, x1, ts2; local i,j,y,z,S; if vflag<>2 then if vflag<>1 then ergs:=args: if nargs=1 and type(ergs,polynom) and type(ergs,polynom(anything,X)) then vis:=collect(ergs,X):n:=degree(vis):dp:=diff(vis,X): if resultant(vis,dp,X)=0 then if dp=0 then print(qts): else print(`args has multiple factor`,qtr); fi; iee:=0:RETURN(): fi; z7:=1: if bouv then x:='x'; for z in roots(prac) do if rem(vis,X-z[1],X)=0 then z7:=z7*(X-z[1]):fi; od: else do_tor(); for z in [PP] do if z[3]=2 and rem(vis,X-z[1],X)=0 then z7:=z7*(X-z[1]):fi; od: fi; z7:=collect(z7,X): if degree(z7,X)=1 then vis:=quo(vis,z7,X):fi; if degree(vis)=0 then diso:=2: else diso:=2*n+1:div(diso): if rem(al.diso,vis,X)<>0 then diso:=diso-1:div(diso): if rem(al.diso,vis,X)<>0 then lprint(\ `The roots of the poly. must be the X-coordinates, unduplicated, of all the`); lprint(\ `points in a cyclic subgroup`,qtr); iee:=0:RETURN(): fi; fi; fi; vflag:=3:vis:=collect(vis*z7,X): else if bouv then if nargs>3 or (nargs>1 and not type(args[1],algebraic)) then lprint(` For an irrational curve the kernel must be specified by a single`); lprint(`generator--a point of finite order in one of the forms [x,y] or`); lprint(`[x,y,n] (n=order) or x,y or x,y,n or simply x.`); print(qts):RETURN(); elif nargs=1 and type(ergs,algebraic) then ergs:=absc(ergs): fi; fi; subg(ergs); if diso=0 or GG=NULL then RETURN();fi; if diso>Nsubg then eprint(\ `Proposed kernel has order > the allowed limit Nsubg whose current value`); fprint(`is`,Nsubg,` This value can be reset at any time.`); RETURN(); fi; if vog=0 then eprint(`The degree of the isogeny is `.diso);fi; if diso=1 then RETURN();fi; if bouu then do_tor();S:=GG:veg:=0: for i to diso-1 do z:=GG[i]: if veg2 and member(z,S) then x:=z[1]:S:=S minus {[x,-z[2]-a1*x-a3,z[3]]}: fi; od; if veg[] and not member(z[1],S) do S:=S union {z[1]}:z:=eadd(z,z1): od; fi; if vflag=3 then vflag:=0 else vis:=1: for z in S do vis:=vis*(X-z):od; vis:=collect(vis,X): fi; fi; dp:=diff(vis,X);d2p:=diff(dp,X):d3p:=diff(d2p,X):n:=degree(vis,X): for i to 3 do s.i:=(-1)^i*coeff(vis,X,n-i):od; p2:=s1^2-2*s2:aav:=4*X^3+b2*X^2+2*b4*X+b6:ts:=diff(aav,X)/2: t:=6*p2+b2*s1+n*b4:w:=10*(s1^3-3*s1*s2+3*s3)+2*b2*p2+3*b4*s1+n*b6: x7:=(2*n+1)*X-2*s1-ts*dp/vis+aav*(dp^2-d2p*vis)/vis^2: y7:=(2*n+1)*Y+n*a3+a1*s1-(3*a1*X^2+6*(2*Y+a3)*X+b2*\ Y+2*a2*a3-a1*a4)*dp/vis+(7*a1*X^3+(18*Y+a1*b2+9*a3)*X^2+(3*b2*Y+2*a1^2\ *a3+a1*a4+6*a2*a3)*X+(3*b4*Y+a1*a3^2+3*a3*a4-2*a1*a6))*(dp^2-d2p*vis)/vis^2-\ aav*(Y+(a1*X+a3)/2)*(d3p*vis^2-3*d2p*dp*vis+2*dp^3)/vis^3: if type(diso,even) then if bouv then i:=GG:else i:=[PP]:fi; for z in i do if z[3]>2 then next:fi; X:=z[1]: if vis=0 then x1:=X:break:fi; X:='X': od; if X<>'X' then ts2:=-ts/2:t:=t+ts2:w:=w+x1*ts2:X:='X':x7:=x7+ts2/(X-x1): y7:=y7-ts2*(a1*(X-x1)+(a1*x1+a3)/2+Y)/(X-x1)^2: fi; fi; x7:=normal(x7):y7:=normal(y7): fprint(`The isogenous curve,which now becomes the current apecs curve,`); fprint(`is as follows:`);velu1(1) end: velu1:=proc(vs) local disu; global a4, a6, disog, diso, eT, j, jV, nisog, pc, pcn, trou, us, uV, visog, wel, wIS, Wr, Ws, Wt, Wu, x, x7, y, y7, z; #called by velu and by velu2; the latter is only called by itog (3 times) if bouu then pcn:=ncur:fi; pc:=ncur:disu:=diso; if vs<>0 then a4:=a4-5*t:a6:=a6-b2*t-7*w: if bouu then cld_():fi; fi; if bouu then weier():wIS:=[uV,Wr,Ws,Wt,Wu]:ein(a1,a2,a3,a4,a6); pcn:=0:uV:=wIS[1]:Wr:=wIS[2]:Ws:=wIS[3]:Wt:=wIS[4]:Wu:=wIS[5]: else ell(a1,a2,a3,a4,a6,K_);RETURN(); fi; diso:=disu; if vs=1 or vs=-1 then x7:=(x7*uV^2-Wr)/Wu^2:y7:=(y7*uV^3-Wu^2*x7*Ws-Wt)/Wu^3: fi; if ncur=pc then x:=isqrt(diso):eprint(`the isogeny is just multiplication by `.x); else x:=eT[pc]: if not member(ncur,x[11],'jV') then x:=uae([op(x[11]),ncur],x,11):jV:=nops(x[11]): fi; if vs=0 or vs=2 then us:=-1:z:=uae([],x[13],jV): else us:=vs:z:=uae([collect(vis,X),x7,y7],x[13],jV): fi; y:=uae(us*diso,x[12],jV):x:=uae(y,z,x,12):eT[pc]:=x:j:='j': if not member(pc,nisog,j) then nisog:=[op(nisog),pc]:j:=nops(nisog):disog:=[op(disog),0]: visog:=[op(visog),[]]: fi; y:=disog[j]: if y=0 then disog:=uae(-diso,disog,j):wel:=1:fi; if abs(disog[j])<>diso then print(QT,` velu1`);trou:=[op(trou),VELU1,pc,ergs,cur]: fi; upis(); eprint(`To return to the previous curve use the apecs command go(psn);.`); fi; ncur end: velu2:=proc(d,vs) global diso, ape, iee, wel, eur, cyi; #called only by itog (3 times) diso:=d:ape:=iee:iee:=-1:wel:=0:velu1(vs);eur:=cur:ein(pc);iee:=ape: if wel=1 or (vs<>0 and ftw=1) then fprint(cur,`/subgroup of order `,d,` = `,eur);apq(); if wel=1 then cyi:=cyi*d:fi; fi end: Ver:=proc(); print(`Version `.vern,daat); lprint(`To see the latest version number, and the changes made since this \ version`); lprint(`see the file "README"; the file "install.txt" reminds you how to`); lprint(`update your apecs program:`); lprint(` ftp.math.mcgill.ca`); lprint(` anonymous`); lprint(` your_e_address`); lprint(` cd pub/apecs`); lprint(` get README`); lprint(` get install.txt`); print(); lprint(`>>>>>HIT SEMICOLON ; FOR NEXT SCREEN`);readstat(); lprint(`Please let me know about any bugs; also suggestions for additions a\ nd`); lprint(`improvements are always welcome. My e-mail and ordinary mail addres\ ses are:`); lprint(` connell@math.mcgill.ca`); lprint(` Ian Connell`); lprint(` Mathematics Dept.`); lprint(` McGill University`); lprint(` 805 Sherbrooke W.`); lprint(` Montreal, Quebec`); lprint(` Canada H3A 2K6`) end: ver:=proc(); `Version `.vern,daat end: We:=proc() global wec; if nargs=0 then if bouv then lprint(`The present curve is defined over the field with label K_`=K_); lprint(`(not a catalog curve) by the Weierstrass equation`); fi; wec:=yo: else if ian(args)=2 then RETURN();fi; lprint(`Catalog name and number: `.(et1[ncir]),ncir):wec:=et2[ncir]: fi; print(sort(Y^2+wec[1]*X*Y+wec[3]*Y,[Y,X])=sort(X^3+wec[2]*X^2+wec[4]*X+wec[5],X)) end: we:=proc() global wec; if nargs=0 then wec:=yo else if ian(args)=2 then RETURN();fi; wec:=et2[ncir]: fi end: weier:=proc() global Wu, i, p, d5, d4, d6, taU, Ws, Wr, Wt, a4p, a6p, a3p, a2p, a1p, u, ewq; local c1,c2,c3,x; efad(igcd(c6^2,DD),12);Wu:=1: for i to nP do p:=ntab[i][1]:d5:=ue[i]:Wu:=Wu*p^d5: if p=2 then d4:=c4/2^(4*d5):d6:=c6/2^(6*d5): if modp(d4,2)=1 and modp(d6,4)=3 then taU:=-(d6+1)/4: elif modp(d4,16)=0 and (modp(d6,32)=0 or modp(d6,32)=8) then taU:=d6/8: else Wu:=Wu/2:taU:=0: fi; elif p=3 then if val(c6)=6*d5+2 then Wu:=Wu/3:fi; fi; od; d4:=c4/Wu^4:d6:=c6/Wu^6: c1:=modp(d4,2):Ws:=(Wu*c1-a1)/2: c2:=mods(-d6-c1,3):x:=c1+4*c2:Wr:=(Wu^2*c2-a2+Ws*a1+Ws^2)/3: c3:=modp(taU+c1*c2,2):Wt:=(Wu^3*c3-a3-Wr*a1)/2: a4p:=(x^2-d4-24*c1*c3)/48: a6p:=(36*x*(c1*c3+2*a4p)-x^3-d6-216*c3)/864: a3p:=c3:a2p:=c2:a1p:=c1: u:=Wu+abs(Ws)+abs(Wr)+abs(Wt):ewq:=-Wr/(Wu^2): end: xnk:=proc() global X, Y, UR, VR, X_, Y_; ; X:=args[1]:Y:=args[2]:UR:=X_:VR:=Y_: X_:=sex_(X):Y_:=sex_(Y): UR:='UR':VR:='VR':X:='X':Y:='Y' end: Ypecs:=proc() global pcn, iee; pcn:=1:Zpecs();iee:=0:NULL end: ypecs:=proc() global pcn; pcn:=1:zpecs() end: Zpecs:=proc() global iee; print(`Please wait while catalog is being updated.`); iee:=1:zpecs() end: zpecs:=proc() global ans, pcn; oan(); if nops(trou)>0 then print(`there was trouble: trou = `,trou);ans:=readstat(`Do you wish to s\ tore eT anyway? Type yes; or no; `): if ans=no then if pcn=0 then quit; else pcn:=0:NULL;RETURN(); fi; fi; fi; save et1,et2,etj,eT,cat;eprint(`Catalog is now updated on disk.`); if pcn=0 then eprint(`***************`); eprint(`| Au Revoir |`); eprint(`***************`);quit; fi; pcn:=0:NULL end: