\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ \\ Author: Sungkon Chang \\ The following is the implementation of the 3-descent on an elliptic curve y^2=x^3+b \\ Update: Nov 2003. \\ For more information, read the technical report \\ http://www.math.uga.edu/~schang/math/selmer3.pdf \\ \\ It is written for the Gp-Pari version 2.1.5; it may not be going to work on a higher version \\ as the newer version changed some output formats. \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ print("It is written for the Gp-Pari version 2.1.5; it may not work on a higher version"); print("Please, report bugs to schang@math.uga.edu"); { selmer(b,pN)= local(); b=makesixpowerfree(b); if(b==-2^2*3^3 || b==2^4, print("The following cases have not been implemented yet."); print("Case: b=-2^2*3^3, or 2^4 up to sixth power"), if(b==0,print("b shouldn't be 0"), Selmer3(b,pN)) ) } { makesixpowerfree(b)= local(k,t,T,n,sGn); t=abs(b); if(b>0,sGn=1,sGn=-1); T=factor(t); n=matsize(T)[1]; for(k=1,n, T[k,2]=T[k,2]%6); factorback(T)*sGn } { \\ b is sixth power free. Selmer3(b,pN)= local(n,m,S,k,j,dimm); dimm=0; if(biscube(4*b)>0,MGen=H1generator(b),MGen=H1generator2(b)); \\ It also generates VV,WW for A1 and A2. n=matsize(MGen)[1]; m=matsize(MGen)[2]; S=theS(b); setupfortestat3(b,pN); setupfortestatothers(b,S,pN); for(k=1,n, for(j=1,m, if( MGen[k,j]==1 && localcond(k,j,S)==1 && localat3(k,j,pN)==1, dimm=dimm+1, MGen[k,j]=0 ) ) ); k=0; while(dimm >= 3, dimm=dimm/3; k=k+1); k } { setupfortestatothers(b,S,pN)= local(q,k,kk); littlecubes=vector(length(S)-1); bigcubes=vector(length(S)-1); for(k=1,length(S)-1, q=S[k+1]; littlecubes[k]=cubesatother(A1,q); bigcubes[k]=cubesatother(A2,q) ) } { cubesatother(L,q)= local(k,T,P,N,geN,NN,limiT,thecoeff,thelist); if(L==A1,T=SA1,T=SA2); P=picktheprime(T,q); T=idealstar(L,P,0); N=length(T[2]); geN=vector(N); limiT=vector(N); for(k=1,N, if(T[2]%3==0, geN[k]=nfeltpow(L,T[3][k],3); limiT[k]=T[2][k]/3, geN[k]=T[3][k]; limiT[k]=T[2][k] ) ); NN=prod(k=1,N,limiT[k]); thelist=vector(NN); if(N==0,NN=0; thelist=[nfalgtobasis(L,1)],); thecoeff=vector(N); if(N > 0,thecoeff[1]=-1,); for(k=1,NN, thecoeff=carry3(thecoeff,limiT);\\ carry3 defined in localat3_2.gp thelist[k]=putiT(L,geN,thecoeff) ); thelist } { picktheprime(T,q)= local(k); k=1; while(k <= length(T) && T[k][1]!=q, k=k+1); T[k] } { localcond(kk,jj,S)= local(a1,a2,a22,P,n,q,k,thepower,sB); P=1; n=length(S); a1=totheA1(kk); a2=totheA2(jj); for(k=2,n, q=S[k]; sB=if(isbsq==0,thesqrt(b)+O(q^pN),polrootspadic(x^2-b,q,pN)[1]); thepower=testatA1_q(a1,q,sB); if(thepower==-1,P=0, P=testatA2_q(a2,q,thepower,sB) ) ); P } { totheA1(k)= local(n,a1,f,j); n=length(Pass1); f=A1.nf[1]; a1=prod(j=1,nn, Mod(VV[j]^(Pass1[k][j]),f) ); a1=component(a1,2) } { totheA2(j)= local(m,a2,f,k); m=length(Pass2); f=A2.nf[1]; a2=prod(k=1,mm, Mod(WW[k]^(Pass2[j][k]),f) ); a2=component(a2,2) } { totheA22(j)= local(m,a22,f,k); m=length(Pass2); f=A22.nf[1]; a22=prod(k=1,mmm, Mod(WWW[k]^(Pass2[j][k+mm]),f) ); a22=component(a22,2) } { testatA1_q(aa,q,sB)= local(k,w,ww,w0,pass,triv,f); triv=1+O(q^pN); f=A1.nf[1]; w0=Mod(4*b,f); aa=subst(aa,x,sB); w=Mod(aa,f); k=-1; pass=0; while(k < 2 && pass==0, k=k+1; ww=w/w0^k; ww=ww*triv; pass=wiscubeatq(A1,ww,q) ); if(pass==0,-1,k) } { testatA2_q(aa,q,thepower,sB)= local(k,w,P,pass,omega2_q,triv); triv=1+O(q^pN); omega2_q=find_omega2_q(q,sB); \\ omega_q in poly mod A2.nf[1] with coeff q-adic. w=aa/omega2_q^thepower; w=w*triv; wiscubeatq(A2,w,q) } { find_omega2_q(q,sB)= local(f,w); f=A2.nf[1]; w=Mod(x,f); w=(w^3-9*b*w-4*b)/(3*w^2-3*b);\\ w=sqrt(-3b) w=w/sB; w=(-1+w)/2; w^2 } { \\ w is a polmod mod Q-polynomial with coeff q-adic. wiscubeatq(L,w,q)= local(P,T,pi,n,ans,ww,k); if(L==A1,T=SA1,T=SA2); P=picktheprime(T,q); pi=P[2]; w=component(w,2); w=fitiT_q(w,q); w=nfalgtobasis(L,w); n=nfeltval(L,w,P); if(n%3!=0,ans=0, ww=nfeltpow(L,pi,n); w=nfeltdiv(L,w,ww); \\ w has valuation 0 at P. \\ w has valuation >=0 at other P' lying over q k=kthprime(S,q); if(L==A1,T=littlecubes[k],T=bigcubes[k]); ans=testit_q(L,w,P,T,q) ); ans } { kthprime(S,q)= local(k); k=1; while(S[k]!=q, k=k+1); k-1 } { testit_q(L,w,P,thelist,q)= \\ w has valuation 0 at P. local(N,counter,PP,ww); PP=0; N=length(thelist); w=makeintegeR_q(w,q); counter=0; while( counter < N && PP==0, counter=counter+1; ww=thelist[counter]-w; ww=nfeltval(L,ww,P); if(ww > 0, PP=1, ) ); PP } { \\ Reduce coeffs mod q, \\ and the result is not supposed to be the zero vector \\ because w has valuation 0 at P. makeintegeR_q(w,q)= local(k,ww,n); n=length(w); for(k=1,n, ww=w[k]; ww=Mod(w[k],q); ww=component(ww,2); w[k]=ww ); w; } { \\ ww is a polynomial in x with 3-adic coeff. \\ fitiT(ww) raitionalizes the coefficients of ww \\ by multiplying smallest 3^(3*t) for some t>=0. \\ So, its class up to cubes doesn't change. fitiT_q(ww,q)= local(n,m,mm,S,ind,cFF,preC,preCC,k); n=poldegree(ww); S=0; m=valuation(polcoeff(ww,0),q); preC=padicprec(polcoeff(ww,0),q); for(k=1,n, cFF=polcoeff(ww,k); preCC=padicprec(cFF,q); mm=valuation(cFF,q); if(mm=0,m=0, m=-m; if(m%3==0,,m=3*truncate(m/3)+3) ); ww=ww*q^m; for(k=0,n, cFF=polcoeff(ww,k); cFF=Mod(cFF,q^(preC+m));cFF=component(cFF,2); S=S+x^k*cFF ); S } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ H1_alt2.gp is copied here. \\ It requires the output of fullgenerator() for the fields A1 and A2; \\ Definitions: { thingsforH1(b)= local(cube4b,u,u1,u2,counter); checkgroundfields(b); if(isbsq==0,f1=x-thesqrt(b),f1=x^2-b); if(is4bcube==0 && isn3bsq==0, baseQ2=0; baseQ22=0; g2=t-thecubert(4*b); a2=Mod(t,g2); f2=(x-a2)-thesqrt(-3*b); G2=bnfinit(g2), ); if(is4bcube==0 && isn3bsq==1, baseQ2=0; baseQ22=1; g2=t-thecubert(4*b); a2=Mod(t,g2); trivialone=Mod(1,g2); g22=t^2+thecubert(4*b)*t+thecubert(4*b)^2; a22=Mod(t,g22); trivialone22=Mod(1,g22); f2=(x-thecubert(4*b))^2+3*b*trivialone; G2=bnfinit(g2); f22=(x-a22)^2+3*b; G22=bnfinit(g22), ); if(is4bcube==1 && isn3bsq==0, baseQ2=1; baseQ22=1; g2=t^3-4*b; a2=Mod(t,g2); f2=(x-a2)-thesqrt(-3*b); G2=bnfinit(g2), ); if(is4bcube*isn3bsq==1, baseQ2=1; baseQ22=1; g2=t^3-4*b; a2=Mod(t,g2); f2=(x-a2)^2+3*b; G2=bnfinit(g2), ); trivialone2=Mod(1,g2); if(isbsq==0 && is4bcube==0, baseQ3=0; baseQ33=1; g3=t-(thesqrt(b)+thecubert(4*b)); a3=Mod(t,g3); trivialone=Mod(1,g3); sb=thesqrt(b); f3=(x-thecubert(4*b))^2 +(3*b)*trivialone; G3=bnfinit(g3); omegA=(-1+(Mod(x,f3)-a3+sb)/sb)/2, ); if(isbsq==1 && is4bcube==0 && makesquarefree(b)==-3, baseQ3=-1; baseQ33=0; g3=(t-thecubert(4*b))^2-b; a3=Mod(t,g3); trivialone=Mod(1,g3); sb=a3-thecubert(4*b); f3=x-thecubert(4*b) -thesqrt(-3*b)*trivialone; G3=bnfinit(g3); omegA=(-1+(Mod(x,f3)-a3+sb)/sb)/2, ); if(isbsq==1 && is4bcube==0 && makesquarefree(b)!=-3, baseQ3=-1; baseQ33=1; u=thecubert(4*b); g3=(t-u)^2-b; a3=Mod(t,g3); trivialone=Mod(1,g3); \\sb=a3-u; sb=(a3^3+3*b*a3-u^3)/(3*a3^2+b); f3=((x-u)^2 +(3*b))*trivialone; G3=bnfinit(g3); omegA=(-1+(Mod(x,f3)-a3+sb)/sb)/2; g33=t^4+2*u*t^3+(3*u^2-2*b)*t^2+(2*u^3-2*b*u)*t+u^4+b*u^2+b^2; a33=Mod(t,g33); trivialonee=Mod(1,g33); sbb=(a33^3+3*b*a33-u^3)/(3*a33^2+b); omegAA=(a33-sbb)/u; nthreeb=sbb*(2*omegAA+1); G33=bnfinit(g33), ); \\ f33=(x-a33+sbb)^2+(3*b);f33=f33*trivialonee; G33=bnfinit(g33), ); if(isbsq==0 && is4bcube==1, baseQ3=1; baseQ33=1; g3=(t-thesqrt(b))^3-4*b; a3=Mod(t,g3); trivialone=Mod(1,g3); sb=thesqrt(b); f3=(x-a3+sb)^2+3*b; G3=bnfinit(g3); omegA=(-1+(Mod(x,f3)-a3+sb)/sb)/2, ); if(isbsq*is4bcube==1 && makesquarefree(b)==-3, baseQ3=2; baseQ33=0; g3=t^6-3*b*t^4-8*b*t^3+3*b^2*t^2-24*b^2*t+16*b^2-b^3; a3=Mod(t,g3); trivialone=Mod(1,g3); sb=(a3^3+3*b*a3-4*b)/(3*a3^2+b); f3=x-(thesqrt(-3*b) +a3-sb); G3=bnfinit(g3); omegA=(-1+(Mod(x,f3)-a3+sb)/sb)/2, ); if(isbsq*is4bcube==1 && makesquarefree(b)!=-3, baseQ3=2; baseQ33=1; g3=t^6-3*b*t^4-8*b*t^3+3*b^2*t^2-24*b^2*t+16*b^2-b^3; a3=Mod(t,g3); trivialone=Mod(1,g3); sb=(a3^3+3*b*a3-4*b)/(3*a3^2+b); f3=(x-a3+sb)^2+3*b; G3=bnfinit(g3); omegA=(-1+(Mod(x,f3)-a3+sb)/sb)/2, ); trivialone3=Mod(1,g3); if(is4bcube==0 && isn3bsq==0, baseQ4=0; baseQ44=0; g4=t-thesqrt(-3*b); a4=Mod(t,g4); f4=x-(a4+thecubert(4*b)), ); if(is4bcube==1 && isn3bsq==0, baseQ4=0; baseQ44=1; g4=t-thesqrt(-3*b); a4=Mod(t,g4); f4=(x-a4)^3-4*b, ); if(is4bcube==0 && isn3bsq==1, baseQ4=1; baseQ44=0; u=thecubert(4*b); g4=t^2+3*b; a4=Mod(t,g4); f4=(x-a4)-u; G4=bnfinit(g4), ); if(is4bcube==1 && isn3bsq==1, baseQ4=1; baseQ44=1; g4=t^2+3*b; a4=Mod(t,g4); f4=(x-a4)^3-4*b; G4=bnfinit(g4), ); trivialone4=Mod(1,g4); \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Global_A1(b); VV=fullgenerator(); nn=length(VV); Global_A2(b); WW=fullgenerator(); mm=length(WW); if(is4bcube==0, Global_A22(b); WWW=fullgenerator(); mmm=length(WWW), ); VV1=vector(nn);for(k=1,nn,VV1[k]=norm(Mod(VV[k],f1)) ); WW2=vector(mm);for(k=1,mm,WW2[k]=norm(Mod(WW[k]*trivialone2,f2)) ); if(is4bcube==0, WW22=vector(mmm);for(k=1,mmm,WW22[k]=norm(Mod(WWW[k]*trivialone22,f22)) ), ); VV3=vector(nn);for(k=1,nn,VV3[k]=subst(VV[k],x,sb) );\\ sb in Polmod if(is4bcube==0, VV33=vector(nn);for(k=1,nn,VV33[k]=subst(VV[k],x,sbb) ), ); WW3=vector(mm); triv=Mod(1,g3); for(k=1,mm, WW3[k]=norm(subst(WW[k],x,x+(a3-sb)*(omegA-1))*Mod(triv,f3) )); if(is4bcube==0, WW3=vector(mmm); triv=Mod(1,g3); for(k=1,mmm, WW3[k]=norm(subst(WWW[k],x,x+(a3-sb)*(omegA-1))*Mod(triv,f3) )), ); if(is4bcube==0, WW33=vector(mmm); triv=Mod(1,g33); for(k=1,mmm, WW33[k]=subst( WWW[k],x,u*omegAA-nthreeb ) ), ); WW4=vector(mm);for(k=1,mm,WW4[k]= norm( Mod(WW[k]*trivialone4,f4) ) ) } { H1generator(b)= local(vv,ww,pass1,counter1,pass2,counter2,dd,PassM,k,j); thingsforH1(b); vv=vector(nn);vv[1]=-1; ww=vector(mm);ww[1]=-1; counter1=0; counter2=0; pass1=vector(3^nn);pass2=vector(3^mm); for(k=1,3^nn, dd=coord_VV(); if( Test1(dd)==0,, counter1=counter1+1; pass1[counter1]=dd ) ); for(k=1,3^mm, dd=coord_WW(); if( Test2(dd)*Test4(dd)==0,, counter2=counter2+1; pass2[counter2]=dd ) ); Pass1=vector(counter1); Pass2=vector(counter2); for(k=1,counter1, Pass1[k]=pass1[k]); for(k=1,counter2, Pass2[k]=pass2[k]); PassM=matrix(counter1,counter2); for(k=1,counter1, for(j=1,counter2, if(Test3(Pass1[k],Pass2[j])==1,PassM[k,j]=1,) ) ); PassM } { H1generator2(b)= local(vv,ww,pass1,counter1,pass2,counter2,dd,PassM,k,j); thingsforH1(b); vv=vector(nn);vv[1]=-1; ww=vector(mm+mmm);ww[1]=-1; counter1=0; counter2=0; pass1=vector(3^nn);pass2=vector(3^(mm+mmm)); for(k=1,3^nn, dd=coord_VV(); if( Test1(dd)==0,, counter1=counter1+1; pass1[counter1]=dd ) ); for(k=1,3^(mm+mmm), dd=coord_WW(); if( Test22(dd)*Test42(dd)==0,,counter2=counter2+1; pass2[counter2]=dd ) ); Pass1=vector(counter1); Pass2=vector(counter2); for(k=1,counter1, Pass1[k]=pass1[k]); for(k=1,counter2, Pass2[k]=pass2[k]); PassM=matrix(counter1,counter2); for(k=1,counter1, for(j=1,counter2, if(Test32(Pass1[k],Pass2[j])==1,PassM[k,j]=1,) ) ); PassM } { \\ Requires ZZ in Polmod or to be rational numbers. tothebase(v)= local(P,k); P=1; for(k=1,length(ZZ), P=P*ZZ[k]^v[k]); P } { coord_VV()= local(); vv[1]=vv[1]+1; vv=carry(vv); vv } { coord_WW()= local(); ww[1]=ww[1]+1; ww=carry(ww); ww } { carry(u)= local(n,k); n=length(u); for(k=1,n, if(u[k]>2,u[k]=u[k]-3;u[k+1]=u[k+1]+1,) ); u } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { \\ Requires g2,g3,g4 the functions defining the ground field of \\ the norm condition. \\ And F2,F3,F4 the ground fields defined by the g's Test1(dd)= local(w,ZZ); ZZ=VV1; w=tothebase(dd); if(isbsq==0,1,wiscube1(w)) } { Test2(dd)= local(w,ZZ,L,g,k); if(baseQ2==0, ZZ=vector(mm); for(k=1,mm,ZZ[k]=component(WW2[k],2)); w=tothebase(dd); wiscube1(w), ZZ=WW2; L=G2; g=g2; w=tothebase(dd); if(isn3bsq==0,1,wiscube(w)) ) } { Test22(dd)= local(w1,w2,ZZ,L,g,k,dd1,dd2,P); dd1=vector(mm); for(k=1,mm,dd1[k]=dd[k]); dd2=vector(mmm); for(k=1,mmm,dd2[k]=dd[mm+k]); ZZ=WW2; L=G2; g=g2; w1=tothebase(dd1); P=wiscube(w1); ZZ=WW22; L=G22; g=g22; w2=tothebase(dd2); P=P*wiscube(w2); P } { Test4(dd)= local(w,ZZ,L,g,k); if(baseQ4==0, ZZ=vector(mm); for(k=1,mm, ZZ[k]=component(WW4[k],2)); w=tothebase(dd); wiscube1(w), ZZ=WW4; L=G4; g=g4; w=tothebase(dd); if(is4bcube==0,1,wiscube(w)) ) } { Test42(dd)= local(w,L,g,k); L=G4; g=g4; w=specialnorm(dd); wiscube(w) } { specialnorm(T)= local(k,u1,u2,u); u1=prod(k=1,mm, Mod(WW[k]^T[k],A2.nf[1]) ); u2=prod(k=1,mmm, Mod(WWW[k]^T[mm+k],A22.nf[1]) ); u=chinese(u1,u2); u=component(u,2); u=u*Mod(1,g4); u=Mod(u,(x-a4)^3-4*b); norm(u) } { Test3(ddd,eee)= local(w,ZZ,L,g,k); if(baseQ3==0, ZZ=vector(nn); for(k=1,nn, ZZ[k]=component(VV3[k],2)); w=tothebase(ddd); ZZ=vector(mm); for(k=1,mm,ZZ[k]=component(WW3[k],2)); w=w*tothebase(eee); wiscube1(w), L=G3; g=g3; w=1; ZZ=VV3; w=tothebase(ddd); ZZ=WW3; w=w*tothebase(eee); if(isn3bsq==0,1,wiscube(w)) ) } { Test32(ddd,eee)= local(w,ZZ,L,g,k,u1,u2,eee1,eee2,P,u,tt); u=thecubert(4*b); eee1=vector(mm);for(k=1,mm,eee1[k]=eee[k]); eee2=vector(mmm);for(k=1,mmm,eee2[k]=eee[mm+k]); ZZ=VV3; L=G3; g=g3; u1=tothebase(ddd); ZZ=WW3; u1=u1*tothebase(eee2); P=wiscube(u1); ZZ=WW; for(k=1,mm, ZZ[k]=Mod(WW[k],A2.nf[1])); tt=tothebase(eee1); tt=component(tt,2); tt=subst( tt,x,u+nthreeb ); ZZ=VV33; L=G33; g=g33; u2=tothebase(ddd); ZZ=WW33; u2=u2*tt*tothebase(eee2); P=P*wiscube(u2); P } { wiscube1(w)= local(ww,n,P); P=1; w=abs(w); ww=factor(w); n=matsize(ww)[1]; for(k=1,n, if(ww[k,2]%3==0,P=P*1,P=P*0) ); P } { \\ Requries g defining the ground field \\ and L the field. \\ w is polymod. wiscube(w)= local(P,J,n,V,m,k,gen); P=1; J=idealprincipal(L.nf,w); J_factor=idealfactor(L,J); n=matsize(J_factor)[1]; for(k=1,n, if(J_factor[k,2]%3==0, P=P*1; J_factor[k,2]=J_factor[k,2]/3,P=P*0) ); if(P==0,, J=putitback(J_factor); n=length(L.clgp[2]); V=bnfisprincipal(L,J); if(V[1]==vector(n)~, gen=L.zk*V[2]; gen=Mod(gen,L.nf[1]); gen=w/gen^3, P=0) ); n=length(L.fu); m=L.tu[1]; if(P==0,, V=bnfisunit(L,gen); for(k=1,n, if(V[k]%3==0,,P=0) ); if(P==1 && m==6, if(V[n+1]*2==Mod(0,6),,P=0), ) ); P } { \\ Requires the field L \\ and R an output by idealfactor, a matrix. putitback(R)= local(D,n,k); D=idealprincipal(L,1); n=matsize(R)[1]; for(k=1,n, D=idealmul(L,D,idealpow(L,R[k,1],R[k,2]))); D } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ theS_alt.gp is copied here. \\ The followings are miscellaneous (user-defined) functions: { \\ b is an integer. \\ 0 when it's a square; a positive number when it's not. bissquare(b)= local(sq,n,s,k); n=0; if(b<0,n=1, sq=factor(b); s=matsize(sq)[1]; for(k=1,s, n=n+sq[k,2]%2) ); n } { \\ b is an integer. \\ 0 when it's a cube; a positive number when it's not. biscube(b)= local(sq,n,s,bb,k); n=0; bb=abs(b); sq=factor(bb); s=matsize(sq)[1]; for(k=1,s, n=n+sq[k,2]%3); n } { makesquarefree(b)= local(a,n,sgnn,q,k); if(b>0,sgnn=1,sgnn=-1); a=factor(abs(b)); n=matsize(a)[1]; for(k=1,n, if(a[k,2]%2==0,, sgnn=sgnn*a[k,1]) ); sgnn } { \\ The following variables are certain indicators for fields involved. checkgroundfields(b)= local(); if(bissquare(-3*b)==0, isn3bsq=0, isn3bsq=1); if(biscube(4*b)==0, is4bcube=0, is4bcube=1); if(bissquare(b)==0, isbsq=0, isbsq=1) } { \\ a required to be a square (i.e., a>0) thesqrt(a)= local(F,n,P,k); P=1; F=factor(a); n=matsize(F)[1]; for(k=1,n,P=P*F[k,1]^(F[k,2]/2)); P } { \\a required to be a cube. thecubert(a)= local(F,n,P,sgnn,k); P=1; if(a>0,sgnn=1,sgnn=-1); a=abs(a); F=factor(a); n=matsize(F)[1]; for(k=1,n,P=P*F[k,1]^(F[k,2]/3)); P*sgnn } { \\ count the number of 1's counttheones(MM)= local(n,m,counter); counter=0; n=matsize(MM)[1]; m=matsize(MM)[2]; for(k=1,n,for(j=1,m, if(MM[i,j]==1,counter=counter+1,) ) ); if(counter==1,0,counter/3) } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ Assume: b is a sixth-power free integer. { \\ It computes the subset S of the rational primes. theS(b)= local(F,sizeofF,S,counter,n,q,qadic,f,SS,bb,k); bb=abs(b); F=factor(4*bb); sizeofF=matsize(F)[1]; S=vector(sizeofF+1); S[1]=3; counter=1; for(k=1,sizeofF, n=F[k,2];q=F[k,1]; if(n%6==0 || q==3,, f=x^2-4*b; qadic=matsize(factorpadic(f,q,10))[1]; if(qadic==1,, counter=counter+1; S[counter]=q) ) ); SS=vector(counter); for(k=1,counter, SS[k]=S[k]); SS; } { \\ K is a number field, and S the set of primes. \\ It finds the set of primes of K lying over S. primesoverS(K,S)= local(SS,SSS,numberofprimes,counter,n,k); numberofprimes=length(S); SS=vector(numberofprimes*3); counter=0; for(k=1,numberofprimes, A=idealprimedec(K,S[k]); n=length(A); for(j=1,n, counter=counter+1; SS[counter]=A[j]) ); SSS=vector(counter); for(k=1,counter,SSS[k]=SS[k]); SSS } { \\ A1=Q(sqrt(b)). \\ It computes the generators of the S-unit group: SunitA1[1], \\ and S-class group: SunitA1[5]. SgroupA1(b)= local(f); if(isbsq==0, f=x-thesqrt(b), f=x^2-b); A1=bnfinit(f); SA1=primesoverS(A1,S); SunitA1=bnfsunit(A1,SA1); SunitA1[5] } { SgroupA2(b)= local(f); if(isn3bsq==0 && is4bcube==0, f=x-(thesqrt(-3*b)+thecubert(4*b)),); if(isn3bsq==0 && is4bcube==1, f=(x-thesqrt(-3*b))^3-4*b,); if(isn3bsq==1 && is4bcube==0, f=(x-thecubert(4*b))^2+3*b,); if(isn3bsq*is4bcube==1, f= x^6+9*b*x^4-8*b*x^3+27*b^2*x^2+72*b^2*x +27*b^3+16*b^2,); A2=bnfinit(f); SA2=primesoverS(A2,S); SunitA2= bnfsunit(A2,SA2); SunitA2[5] } { \\ This auxilary field is considered when 4b is a cube. SgroupA22(b)= local(f,u); u=thecubert(4*b); f= x^4+2*u*x^3+(3*u^2+6*b)*x^2 +(6*b*u+2*u^3)*x+u^4+9*b^2-3*u^2*b; A22=bnfinit(f); SA22=primesoverS(A22,S); SunitA22= bnfsunit(A22,SA22); SunitA22[5] } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ Global_A1(b) and Global_A2(b) compute information \\ about the fields A1 and A2. \\ About A1: \\ S is the special set of rational primes. \\ A1 is renamed K for the purpose of running fullgenerator(). \\ T is a gp-description of the S-class group and generators. \\ Cl_S_gen is the set of generators of the class group of A1 \\ which is chosen for a good description of S-class group. \\ S_gen is the set of primes of A1 lying over S. \\ Sunits is the set of generators of the S-unit group \\ (modulo the fundamental unit group). \\ T[2] is a (diagonal) matrix that describes the relations of \\ the generators of the class group of A1 \\ which hold in the S-class group. \\ SClassno is the S-class number. { \\ M is the transition matrix \\ from the coordinate of to that of A1.clgp \\ e.g., SA1=[P1,P2,P3], and [C1,C2] are the generator of the Cl(K) \\ [1,0,0]~=P1, and M*[1,0,0]~=[a,b]~=C1^a C2^b. \\ bnfisprincipal solves: P1=C1^s C2^t. \\ Then the transpose [s,t]~ forms a column vector of M. Global_A1(b)= local(TT,n,a,k,j,m); S=theS(b); T=SgroupA1(b); K=A1; Cl_S_gen=T[3]; S_gen=SA1; Sunits=SunitA1[1]; TT=T[2]; if(TT==[],SClassno=1;S_N=0, S_N=0; SClassno=T[1]; N=length(K.clgp[2]); d=K.clgp[2]; a=1; while(a==1 && N>S_N, a=TT[N-S_N,N-S_N]; S_N=S_N+1); S_N=S_N-1; m=length(S_gen); M=matrix(N,m); for(k=1,m, TT=bnfisprincipal(K,S_gen[k])[1]; for(j=1,N,M[j,k]=TT[j]) ) ) } { Global_A2(b)= local(TT,n,a,k,j,m); S=theS(b); T=SgroupA2(b); K=A2; Cl_S_gen=T[3]; Sunits=SunitA2[1]; TT=T[2]; if(TT==[],SClassno=1;S_N=0, S_N=0; S_gen=SA2; K=A2; SClassno=T[1]; N=length(K.clgp[2]); d=K.clgp[2]; a=1; while(a==1 && N>S_N, a=TT[N-S_N,N-S_N]; S_N=S_N+1); S_N=S_N-1; m=length(S_gen); M=matrix(N,m); for(k=1,m, TT=bnfisprincipal(K,S_gen[k])[1]; for(j=1,N,M[j,k]=TT[j]) ) ) } { Global_A22(b)= local(TT,n,a,k,j,m); S=theS(b); T=SgroupA22(b); K=A22; Cl_S_gen=T[3]; Sunits=SunitA22[1]; TT=T[2]; if(TT==[],SClassno=1;S_N=0, S_N=0; S_gen=SA22; K=A22; SClassno=T[1]; N=length(K.clgp[2]); d=K.clgp[2]; a=1; while(a==1 && N>S_N, a=TT[N-S_N,N-S_N]; S_N=S_N+1); S_N=S_N-1; m=length(S_gen); M=matrix(N,m); for(k=1,m, TT=bnfisprincipal(K,S_gen[k])[1]; for(j=1,N,M[j,k]=TT[j]) ) ) } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ fullgenerator generates \\ a preimage of K/K^3 -> Cl_S(K)[3] \\ and the generators of the S-unit. \\ To run fullgenerator, one should run either Global_A1 or Global_A2 \\ so that it can use all the information. { \\ Here, S is a local variable. fullgenerator()= local(S,R,r,l,n,m,k); l=K.tu[1]; l=if(l%6==0,l=1,l=0); if(K.clgp[1]==1,allSunits(), S=preimage(); if(S==[1],allSunits();r=0, n=length(K.fu); m=length(Sunits); r=length(S); GNR=vector(l+n+m+r); if(l==1,GNR[1]=component(Mod(K.tu[2],K.nf[1])^2,2), ); for(k=1+l,n+l,GNR[k]=K.fu[k-l]); for(k=n+1+l,n+m+l,GNR[k]=Sunits[k-n-l]); for(k=n+m+1+l,n+m+r+l,GNR[k]=S[k-n-m-l]) ); GNR ) } { \\ S is a local variable. allSunits()= local(S,k,n,m); n=length(K.fu); m=length(Sunits); S=vector(l+n+m); for(k=1+l,n+l, S[k]=K.fu[k-l]); for(k=n+1+l,n+m+l, S[k]=Sunits[k-n-l]); if(l==1,S[1]=component(Mod(K.tu[2],K.nf[1])^2,2), ); S } { preimage()= if(Sclassno%3==0, generators(), [1]) } { generators()= local(Mx,R,k); Mx=vector(N-S_N); R=T[2]; counter=0; for(k=1,N-S_N,dd=R[k,k]; if(dd%3==0, counter=counter+1; B=idealpow(K,Cl_S_gen[k],dd); Mx[counter]=gen_element(B),) ); MM=vector(counter); for(k=1,counter,MM[k]=Mx[k]); MM } { \\ for given an ideal B, find a in K such that B = a A for some A in gen_element(B)= local(T,n,element,zerovector,v); zerovector=vector(N); v=coord_in_S(B); n=length(v); T=idealpow(K,S_gen[1],v[1]); for(k=2,n, T=idealmul(K,T,idealpow(K,S_gen[k],v[k])) ); T=idealmul(K,T,idealinv(K,B)); element=bnfisprincipal(K,T); if(element[1]==zerovector~, coord_x(element[2]), print("something is wrong")) } { \\ coord_x changes the coordinates in the integral basis to \\ an element of Q[x]/f. coord_x(u)=K.zk*u } { \\ Finding B in in terms of coordinates. \\ y is not a local variable; CL(K)-coord of B. coord_in_S(B)= local(); y=bnfisprincipal(K,B)[1]; matsolvemod(M,d,y) } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ localat3_2.gp copied here { \\ It computes: \\ (1) f1at3, f2at3: the irreducible poly over Q3. \\ (2) List1, List2: the cubes for the units. \\ (3) geNN: the generator of the image of the boundary map at q=3. setupfortestat3(b,pN)= findAat3(b,pN); geNN=findtheimageat3(b,pN); } { \\ localat3(...) =1 when it passed, and 0 when not. localat3(kk,jj,pN)= local(a1,a2,P,match,r1,r2); P=0; r1=0; if(length(geNN)==1, r2=2, r2=0); a1=totheA1(kk); a2=totheA2(jj); if(is4bcube==0, a22=totheA22(jj),); while(P==0 && r2<3, P=isinimage1(a1,r1,r2,pN)*isinimage2(a2,r1,r2,pN); if(is4bcube==0,P=P*isinimage22(a22,r1,r2,pN),); r1=r1+1; if(r1==3, r1=0; r2=r2+1,) ); P } \\\\\\\\\\\\\\\\\\\\\\\\\\\\ { \\ It finds a list of cubes of the units mod SA1[1]^((3e+1)/2), \\ and the same for A2. \\ List1 and List2 given in L.zk-basis. findAat3(b,pN)= local(TT,NN,k); if(is4bcube==0,LisT=vector(3),LisT=vector(2)); TT=SA1; NN=primesover3(TT); LisT[1]=vector(NN); for(k=1,NN, LisT[1][k]=cubesat3(A1,TT[k]); ); TT=SA2; NN=primesover3(TT); LisT[2]=vector(NN); for(k=1,NN, LisT[2][k]=cubesat3(A2,TT[k]); ); if(is4bcube==0, TT=SA22; NN=primesover3(TT); LisT[3]=vector(NN); for(k=1,NN, LisT[3][k]=cubesat3(A22,TT[k]); ), ) } { primesover3(T)= local(counter,n); counter=0; n=length(T); while(counter0, ind=1, if(rt2*rt3>0, ind=1, ind=0) ); ind+1 \\ dim E/3E = 1 + dim E[3](Q3). } { \\ The gp-command polrootspadic has deficiency: \\ If the polynomial has non-unit leading coefficient, \\ then often it cannot solve it. \\ So, one has to transform it to monic one. \\ T is a certain signature of this transformation. getsignature(b)= local(m,T); m=valuation(b,3); if(m==0, if(b%9==1 || b%9==8, T=[1,0,1], if(b%3==1,T=[0,1,0],T=[0,0,0]) ) ); if(m==1 || m==5, T=[0,0,0], ); if(m==2, if((b/3^2)%3==1, T=[0,1,1], T=[0,0,0]) ); if(m==3, if((b/3^3)%9==1 || (b/3^3)%9==8, T=[1,1,3], if((b/3^3)%9==2 || (b/3^3)%9==4, T=[1,1,2], T=[0,0,0] ) ) ); if(m==4, if((b/3^4)%3==1, T=[0,2,2], T=[0,0,0]) ); T } { randomsearch(n,T)= local(s,t,M,theEq,P); s=T[2]; t=T[3]; M=min(3*s,2*t); theEq=getequation(T,M); \\ equation to solve, written terms of x and k. \\ k is a unit parameter. P=findP(theEq,T,-1); if(n==2, P=findPP(P,theEq,T), P=[P[1]]); P } { \\ produces a transformed equation. getequation(T,M)= local(p,ss,tt,eqn); p=T[1]; ss=T[2]; tt=T[3]; if(p==0, eqn=3^(2*tt-M)*x^2 - 3^(3*ss-M)*r^3 - b/3^M, eqn=3^(2*tt-M)*r^2 - 3^(3*ss-M)*x^3 - b/3^M ); eqn } { \\ Finds a point not divisible by 3. findP(theEq,T,kk)= local(ind,eqn,solN,PP,inD); inD=0; while(inD==0, ind=0; while(ind==0, if(kk%3==2,kk=kk+2,kk=kk+1); eqn=subst(theEq,r,kk); solN=polrootspadic(eqn,3,pN); if(solN==[]~, , ind=1) ); solN=solN[1]; if(T[1]==0, PP=[kk*3^T[2],solN*3^T[3]], PP=[solN*3^T[2], kk*3^T[3]] ); if(torsioN(PP,b,pN)==0 && divisibility_by_3(PP)==0,inD=1,) ); [PP,kk] } { \\ 1 if divisible by 3, and 0 if not. divisibility_by_3(P)= local(w1,w2,indd); w1=Fmapat3(b,A1,P); w2=Fmapat3(b,A2,P); if(is4bcube==0,w22=Fmapat3(b,A22,P),); indd=cubetestat3(A1,w1)*cubetestat3(A2,w2); if(is4bcube==0,indd=indd*cubetestat3(A22,w22),); if(indd==1,1,0) } { \\ It spits 1 if P is 6-torsion, and 0 if not. torsioN(P,b,pN)= local(xx,yy,xxx,yyy,PP,ind); trivialone=1+O(3^pN); ind=0; if( P[1]^3-b==O(3^pN), ind=1, PP=multbytwo(b,P)); if(ind==0, xx=PP[1]*trivialone; yy=PP[2]*trivialone; xxx=xx*(xx^3+4*b); yyy=yy^2-xx^3-b; if(xxx==O(3^pN) && yyy==O(3^pN), ind=1,), ); ind } { \\ It finds a second generator. \\ P=[first point, kk] \\ where kk is the value of the parameter used \\ for the first point. findPP(P,theEq,T)= local(kk,Q,QQ,ind,thelimit,counter); thelimit=20; Q=P[1]; kk=P[2]; ind=0; counter=0; while(ind==0 && counter<=thelimit, counter=counter+1; QQ=findP(theEq,T,kk); kk=QQ[2]; QQ=QQ[1]; if(independency(Q,QQ)==1,ind=1,); ); if(ind==0, QQ=lastchance(Q,b,pN); if(independency(Q,QQ)==1,ind=1,) ); if(ind=1, [Q,QQ], print("It failed to find all generators");[Q] ) } { \\ Q and QQ are not 3-torsion. independency(Q,QQ)= local(v1,v2,v22,w1,w2,w22,ind,ss); v1=Fmapat3(b,A1,Q); v2=Fmapat3(b,A2,Q); if(is4bcube==0,v22=Fmapat3(b,A22,Q),); w1=Fmapat3(b,A1,QQ); w2=Fmapat3(b,A2,QQ); if(is4bcube==0,w22=Fmapat3(b,A22,QQ),); ind=cubetestat3(A1,v1*w1)*cubetestat3(A2,v2*w2); if(is4bcube==0, ind=ind*cubetestat3(A22,v22*w22), ); ss=cubetestat3(A1,v1^2*w1)*cubetestat3(A2,v2^2*w2); if(is4bcube==0,ss=ss*cubetestat3(A22,v22^2*w22),); ind=ind+ss; if(ind==0,1,0) } { \\ If nothing worked, use the torsion points. lastchance(Q,b,pN)= local(soln1,soln21,soln22,RR); soln1=polrootspadic(x^2-b,3,pN); soln21=polrootspadic(x^3+4*b,3,pN); soln22=polrootspadic(x^2+3*b,3,pN); if(soln1!=[]~,RR=[0,soln1[1]], RR=[soln21[1],soln22[1]] ); addpoints(b,Q,RR) } { theimage(MM,n)= local(k,thevector,PP,QQ,a1,a2,a22); thevector=vector(n); QQ=findP1(b,pN); for(k=1,n, PP=MM[k]; if(torsioN(PP,b,pN)==0, a1=Fmapat3(b,A1,PP); a2=Fmapat3(b,A2,PP); if(is4bcube==0, a22=Fmapat3(b,A22,PP); thevector[k]=[a1,a2,a22], thevector[k]=[a1,a2]), a1=Fmapat3fortor(b,A1,PP,QQ); a2=Fmapat3fortor(b,A2,PP,QQ); if(is4bcube==0, a22=Fmapat3fortor(b,A22,PP,QQ); thevector[k]=[a1,a2,a22], thevector[k]=[a1,a2]) ) ); thevector } \\\\\\\\\\\\\\\\\\\\\\\\\\ \\ The F-map at 3 { \\ When P is not a 2-torsion point. Fmapat3(b,L,P)= local(twoP); twoP=multbytwo(b,P); imG(b,L,twoP,P) } { multbytwo(b,P)= local(x0,y0,x3,y3,lam,mu); x0=P[1]; y0=P[2]; x3=(x0^4-8*b*x0)/(4*x0^3+4*b); lam=3*x0^2/(2*y0); mu=(-x0^3+2*b)/(2*y0); y3=-lam*x3-mu; [x3,y3] } { imG(b,L,Q,P)= local(ff,thex,they,ffat3); if(L==A1,ff=Y-x;ffat3=L.nf[1], they=(x^3-9*b*x-4*b)/(3*x^2-3*b);\\they=sqrt(-3*b) thex=they-x;\\ thex=-sqrt[3](4*b) lam=3*thex^2/(2*they); ff=Y-lam*(X-thex)-they; ffat3=L.nf[1] ); evaL(ffat3,ff,Q)/evaL(ffat3,ff,P) } { evaL(ffat3,ff,R)= local(xxx,yyy,zzz); xxx=R[1]; yyy=R[2]; zzz=subst(subst(ff,X,xxx),Y,yyy); Mod(zzz,ffat3) } { Fmapat3fortor(b,L,P,QQ)= local(PQQ); PQQ=addpoints(b,P,QQ); imG(b,L,PQQ,QQ) } { addpoints(b,P,QQ)= local(x1,y1,x2,y2,x3,y3,lam,mu); x1=P[1];y1=P[2]; x2=QQ[1];y2=QQ[2]; lam=(y2-y1)/(x2-x1); mu=(y1*x2-y2*x1)/(x2-x1); x3=lam^2-(x1+x2); y3=-lam*x3-mu; [x3,y3] } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { isinimage1(a1,r1,r2,pN)= local(P,w,aa,g1,g2,trivialone,ffat3); P=0; trivialone=1+O(3^pN); ffat3=A1.nf[1]; aa=Mod(a1*trivialone, ffat3*trivialone); if(length(geNN)==2, g1=geNN[1][1]*trivialone; g2=geNN[2][1]*trivialone; w=aa/(g1^r1*g2^r2), ); if(length(geNN)==1, g1=geNN[1][1]*trivialone; w=aa/g1^r1, ); cubetestat3(A1,w) } { isinimage2(a2,r1,r2,pN)= local(P,w,aa,g1,g2,trivialone,ffat3); P=0; trivialone=1+O(3^pN); ffat3=A2.nf[1]; aa=Mod(a2*trivialone, ffat3*trivialone); if(length(geNN)==2, g1=geNN[1][2]*trivialone; g2=geNN[2][2]*trivialone; w=aa/(g1^r1*g2^r2), ); if(length(geNN)==1, g1=geNN[1][2]*trivialone; w=aa/g1^r1, ); cubetestat3(A2,w) } { isinimage22(a22,r1,r2,pN)= local(P,w,aa,g1,g2,trivialone,ffat3); P=0; trivialone=1+O(3^pN); ffat3=A22.nf[1]; aa=Mod(a22*trivialone, ffat3*trivialone); if(length(geNN)==2, g1=geNN[1][3]*trivialone; g2=geNN[2][3]*trivialone; w=aa/(g1^r1*g2^r2), ); if(length(geNN)==1, g1=geNN[1][3]*trivialone; w=aa/g1^r1, ); cubetestat3(A22,w) } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { \\ w is in polmod over Q3 . \\ Some coeff of w should be outside O(3^3). \\ 1 if cube, and 0 if not. cubetestat3(L,w)= local(ans,T,n,k,ind); ind=1; if(L==A1,T=SA1,if(L==A2,T=SA2,T=SA22)); n=primesover3(T); for(k=1,n, ind=ind*cubetestatP(L,T,k,w)); ind } { cubetestatP(L,T,k,w)= local(e,P,j); P=T[k]; e=P[3]; if(L==A1,j=1,if(L==A2,j=2,j=3)); Pe=idealpow(L,P,ceil((3*e+1)/2)); Pe=idealhnf(L,Pe); pi=P[2]; w=component(w,2); w=fitiT(w,P); w=nfalgtobasis(L,w); n=nfeltval(L,w,P); if(n%3!=0,ans=0, ww=nfeltpow(L,pi,n); w=nfeltdiv(L,w,ww); \\ w has valuation 0 at P ans=testit(L,w,P,Pe,LisT[j][k]) ); ans } { testit(L,w,P,Pe,thelist)= \\ w has valuation 0 at P. local(k,N,counter,PP,ww,nn,e); ee=P[3]; PP=0; nn=ceil( (3*ee+1)/2 ); N=length(thelist); w=makeintegeR(w); counter=0; while( counter=nn, PP=1, ) ); PP } { \\ Reduce coeffs mod 9, \\ and the result is not supposed to be the zero vector \\ because w has valuation 0 at P. makeintegeR(w)= local(k,ww); n=length(w); for(k=1,n, ww=w[k]; ww=Mod(w[k],9); ww=component(ww,2); w[k]=ww ); w; } { \\ ww is a polynomial in x with 3-adic coeff. \\ fitiT(ww) raitionalizes the coefficients of ww \\ by multiplying smallest 3^(3*t) for some t>=0. \\ So, its class up to cubes doesn't change. fitiT(ww,P)= local(n,m,mm,S,ind,cFF,preC,preCC,k); n=poldegree(ww); S=0; m=valuation(polcoeff(ww,0),3); preC=padicprec(polcoeff(ww,0),3); for(k=1,n, cFF=polcoeff(ww,k); preCC=padicprec(cFF,3); mm=valuation(cFF,3); if(mm=0,m=0, m=-m; if(m%3==0,,m=3*truncate(m/3)+3) ); ww=ww*3^m; for(k=0,n, cFF=polcoeff(ww,k); cFF=Mod(cFF,3^(preC+m));cFF=component(cFF,2); S=S+x^k*cFF ); S }