\\ \\ 8-descent.gp \\ \\ Find a rational point (x,y) of Elliptic curve E: y^2=x^3+ax^2+bx, \\ but not with determining rank(E). \\ \\ \\ find non-trivial point of E_n: y^2=x^3-n^2x congr(n,s1b,s2b,s3b,s4b)=ed(0,-n^2,2,s1b,s2b,s3b,s4b); \\ find non-trivial point of E2_n: y^2=x^3-3nx^2+2n^2x congr2(n,s1b,s2b,s3b,s4b)=ed(-3*n,-2*n^2,2,s1b,s2b,s3b,s4b); \\ find non-trivial point of E3_n: y^2=x^3+3nx^2+2n^2x congr3(n,s1b,s2b,s3b,s4b)=ed(3*n,-2*n^2,2,s1b,s2b,s3b,s4b); \\ find non-trivial point of E4_n: y^2=x^3+2nx^2-3n^2x congrpi3(n,s1b,s2b,s3b,s4b)=ed(2*n,-3*n^2,2,s1b,s2b,s3b,s4b); \\ find non-trivial point of E5_n: y^2=x^3-2nx^2-3n^2x congr2pi3(n,s1b,s2b,s3b,s4b)=ed(-2*n,-3*n^2,2,s1b,s2b,s3b,s4b); \\ ed(a,b,s1a,s1b,s2b,s3b,s4b)=eightdescent(a,b,s1a,s1b,s2b,s3b,s4b); \\ eightdescent(a,b,s1a,s1b,s2b,s3b,s4b) = { local(e); fordiv(b,d, if(issquarefree(d), e = b/d; \\ print([d,e]); for(s1=s1a,s1b, if(eightdescent2(d,a,b,e,s1,s2b,s3b,s4b)>0, return(1) ); if(eightdescent2(-d,a,b,-e,s1,s2b,s3b,s4b)>0, return(1) ) ) ) ); return(0); } eightdescent2(d,a,b,e,s,s2b,s3b,s4b) = { local(h0,f0,g0); local(kk); local(so,so2); local(k0,u0,p0,q0); so=solve3(d,a,e,s); \\ print(so) if(so!=[], \\ print("[d,h0,f0,g0]=",so); h0=so[2]; f0=so[3]; g0=so[4]; kk=g0*abs(a^2-4*b); \\ print("kk=",kk); fordiv(kk,k, if(issquarefree(k), for(s2=2,s2b, \\ print("k=",k); so2=solve6(k,d,g0,s2); if(so2!=[], \\ print("[k0,u0,p0,q0]=",so2); k0=so2[1]; u0=so2[2]; p0=so2[3]; q0=so2[4]; if(eightdescent3(d,a,b,e,s3b,s4b,h0,f0,g0,k0,u0,p0,q0)>0, return(1) ); if(eightdescent3(d,a,b,e,s3b,s4b,h0,f0,g0,k0,u0,p0,-q0)>0, return(1) ) ) ) ) ) ); return(0) } \\ \\ return 1 or 0 \\ eightdescent3(d,a,b,e,s3b,s4b,h0,f0,g0,k0,u0,p0,q0) = { local(so3); local(z1,z2,z3,z4,z5); local(u1,u2,u3,v1,v2,v3); local(zzz); local(k1,t1,r1,s1); local(c1,c2,c3,c4,c5); local(i); local(kkk); z1=k0^3*(a*g0*q0^2+d*f0*q0^2+f0*p0^2-2*h0*p0*q0); z2=4*u0*k0^3*(h0*q0-f0*p0); z3=2*k0^2*(f0*(2*u0^2*k0-d*g0*q0^2+g0*p0^2)-a*g0^2*q0^2); z4=-4*u0*g0*k0^2*(f0*p0+h0*q0); z5=g0^2*k0*(a*g0*q0^2+d*f0*q0^2+f0*p0^2+2*h0*p0*q0); \\ print("[z1,z2,z3,z4,z5]=",[z1,z2,z3,z4,z5]); if((z1 > 0 || z3 > 0 || z5 > 0), \\ zzz=factor2pol2(z1,z2,z3,z4,z5); zzz=factor2pol(z1,z2,z3,z4,z5); if(zzz!=[], u1=zzz[1]; u2=zzz[2]; u3=zzz[3]; v1=zzz[4]; v2=zzz[5]; v3=zzz[6]; kkk=u1^2*v3^2-u1*u2*v2*v3-2*u1*u3*v1*v3+u1*u3*v2^2+u2^2*v1*v3-u2*u3*v1*v2+u3^2*v1^2; \\ print("kkk=",kkk); fordiv(kkk,k1p, if(issquarefree(k1p), \\ print("[kkk,k1p]=",[kkk,k1p]); for(s3=2,s3b, so3=solve19(k1p,v1,v2,v3,s3); if(so3!=[], \\ print("[k1,t1,r1,s1]=",so3); k1=so3[1]; t1=so3[2]; r1=so3[3]; s1=so3[4]; c1=k1^3*(r1^2*u1+r1*s1*u2+s1^2*u3); c2=-2*k1^3*t1*(2*r1*u1+s1*u2); c3=k1^2*(4*k1*t1^2*u1+2*r1^2*u1*v1+2*r1*s1*u1*v2+s1^2*(u2*v2-2*u3*v1)); c4=-2*k1^2*t1*(2*r1*u1*v1+s1*(2*u1*v2-u2*v1)); c5=k1*(r1^2*u1*v1^2+r1*s1*v1*(2*u1*v2-u2*v1)+s1^2*(u1*v2^2-v1*(u2*v2-u3*v1))); \\ print("[c1,c2,c3,c4,c5]=",[c1,c2,c3,c4,c5]); for(s4=2,s4b, for(j=1,(s4-1), i=s4-j; if(eightdescent4(d,a,b,e,h0,f0,g0,k0,u0,p0,q0,k1,t1,r1,s1,v1,v2,c1,c2,c3,c4,c5,i,j)>0, return(1) ); if(eightdescent4(d,a,b,e,h0,f0,g0,k0,u0,p0,q0,k1,t1,r1,s1,v1,v2,c1,c2,c3,c4,c5,i,-j)>0, return(1) ) ) ) ) ) ) ) ) ); return(0) } \\ \\ return 1 or 0 \\ eightdescent4(d,a,b,e,h0,f0,g0,k0,u0,p0,q0,k1,t1,r1,s1,v1,v2,c1,c2,c3,c4,c5,i,j) = { local(p,q,r,s); local(x,y); local(zz); local(uu,vv,ww); zz=c1*i^4+c2*i^3*j+c3*i^2*j^2+c4*i*j^3+c5*j^4; if(zz > 0 && issquare(zz), \\ print("[z,i,j]=",[sqrtint(zz),i,j]); r=i^2*k1*r1-2*i*j*k1*t1+j^2*(r1*v1+s1*v2); s=s1*(i^2*k1-j^2*v1); \\ print("[r,s]=",[r,s]); p=2*u0*k0*r*s-p0*(g0*s^2+k0*r^2); q=q0*(g0*s^2-k0*r^2); \\ print("[p,q]=",[p,q]); uu=a*g0*q^2+d*f0*q^2+p*(f0*p-2*h0*q); vv=g0*(p^2-d*q^2); if((uu%k1)==0 && (vv%k1)==0, uu=uu/k1; vv=vv/k1; if(uu > 0 && vv > 0 && issquare(uu) && issquare(vv), x=d*uu/vv; ww=d*uu^2+a*uu*vv+e*vv^2; if(ww > 0 && issquare(ww), \\ print("x=",x); y=d*sqrtint(uu)*sqrtint(ww)/(sqrtint(vv)*vv); print("[x,y]=",[x,y]); return(1) ) ) ) ); return(0); } \\ sqrt of rational number \\ return sqrt(x) if x >= 0, -1 if x < 0 sqrtrational(x) = { local(n,d); n=numerator(x); d=denominator(x); if(issquare(n) && issquare(d), return(sqrtint(n)/sqrtint(d)) ); -1 } \\ factor z_1r^4+z_2r^3s+z_3r^2s^2+z_4rs^3+z_5s^4 \\ to (u_1r^2+u_2rs+u_3s^2)(v_1r^2+v_2rs+v_3s^2) \\ return [u1,u2,u3,v1,v2,v3] or [] \\ z_1=u_1v_1 \\ z_5=u_3v_3 \\ z_2=u_1v_2+v_1u_2 \\ z_4=u_3v_2+v_3u_2 \\ z_3=u_1v_3+u_2v_2+u_3v_1 \\ \\ return [u1,u2,u3,v1,v2,v3] or [] \\ factor2pol(z1,z2,z3,z4,z5) = { local(u2,v1,v2,v3); local(d,w,ww); if(z1!=0 && z5!=0, fordiv(z1,u1, v1=z1/u1; if(u1<=abs(v1), \\ print("[u1,v1]=",[u1,v1]); fordiv(z5,u3, v3=z5/u3; \\ print("[u3,v3]=",[u3,v3]); d=u1*v3-u3*v1; w=z3-u1*u3-u3*v1; ww=z3+u1*u3+u3*v1; if(d!=0, \\ case d!=0 u2=u3*z2-u1*z4; v2=v1*z4-v3*z2; if(u2%d==0 && v2%d==0, \\ print("[u2*d,v2*d]=",[u2,v2]); u2=u2/d; v2=v2/d; \\ print("[u2,v2]=",[u2,v2]); if(u2*v2==w, \\ print("[u1,u2,u3,v1,v2,v3]=",[u1,u2,u3,v1,v2,v3]); return([u1,u2,u3,v1,v2,v3]) ) ); u2=u3*z2+u1*z4; v2=-(v1*z4+v3*z2); if(u2%d==0 && v2%d==0, \\ print("[u2*d,v2*d]=",[u2,v2]); u2=u2/d; v2=v2/d; \\ print("[u2,v2]=",[u2,v2]); if(u2*v2==ww, \\ print("[u1,u2,u3,v1,v2,v3]=",[u1,u2,-u3,v1,v2,-v3]); return([u1,u2,-u3,v1,v2,-v3]) ) ), \\ case d==0 \\ print("d=0"); if(w!=0, fordiv(w,u2, v2=w/u2; \\ if(u2<=abs(v2) && u1*v2+u2*v1==z2, if(u1*v2+u2*v1==z2, return([u1,u2,u3,v1,v2,v3]) ) ) ); if(ww!=0, fordiv(ww,u2, v2=ww/u2; \\ if(u2<=abs(v2) && u1*v2+u2*v1==z2, if(u1*v2+u2*v1==z2, return([u1,u2,-u3,v1,v2,-v3]) ) ) ) ) ) ) ) ); [] } \\ factor z_1r^4+z_2r^3s+z_3r^2s^2+z_4rs^3+z_5s^4 \\ to (u_1r^2+u_2rs+u_3s^2)(v_1r^2+v_2rs+v_3s^2) \\ return [u1,u2,u3,v1,v2,v3] or [] factor2pol2(z1,z2,z3,z4,z5) = { local(u,v,ff,msize); ff=factor(z1*'r^4+z2*'r^3+z3*'r^2+z4*'r+z5); msize=matsize(ff); \\ print("poly=",ff); if(msize[1]==1 && ff[1,2]==4, u=ff[1,1]; u=u*u; v=u; \\ print("poly[u,v]=",[u,v]); return([polcoeff(u,2,'r),polcoeff(u,1,'r),polcoeff(u,0,'r), polcoeff(v,2,'r),polcoeff(v,1,'r),polcoeff(v,0,'r)]) ); if(msize[1]==2, u=ff[1,1]; if(ff[1,2]==2, u=u*u); v=ff[2,1]; if(ff[2,2]==2, v=v*v); if(ff[1,2]==3, v=u*v;u=u*u); \\ print("poly[u,v]=",[u,v]); return([polcoeff(u,2,'r),polcoeff(u,1,'r),polcoeff(u,0,'r), polcoeff(v,2,'r),polcoeff(v,1,'r),polcoeff(v,0,'r)]) ); if(msize[1]==3 && ff[1,2]==1 && ff[2,2]==1, u=ff[1,1]*ff[2,1]; v=ff[3,1]; if(ff[3,2]==2, v=v*v); \\ print("poly[u,v]=",[u,v]); return([polcoeff(u,2,'r),polcoeff(u,1,'r),polcoeff(u,0,'r), polcoeff(v,2,'r),polcoeff(v,1,'r),polcoeff(v,0,'r)]) ); if(msize[1]==4, u=ff[1,1]*ff[2,1]; v=ff[3,1]*ff[4,1]; \\ print("poly[u,v]=",[u,v]); return([polcoeff(u,2,'r),polcoeff(u,1,'r),polcoeff(u,0,'r), polcoeff(v,2,'r),polcoeff(v,1,'r),polcoeff(v,0,'r)]) ); \\ debug if(msize[1]>1, print("msize=",msize); print("ff=",ff) ); [] } \\ solve k_1t^2=v_1r^2+v_2rs+v_3s^2 with |r|+|s|=s3 \\ return [k_1,t_1,r_1,s_1] or [] solve19(k1,v1,v2,v3,s3) = { local(r); local(kt2,tt); \\ print("[k1,v1,v2,v3,s3]=",[k1,v1,v2,v3,s3]); for(s=1,(s3-1), r=s3-s; kt2=v1*r^2+v2*r*s+v3*s^2; if((kt2%k1)==0, if(kt2 < 0, k1=-k1); tt=kt2/k1; if(issquare(tt), return([k1,sqrtint(tt),r,s]) ) ); kt2 = v1*r^2-v2*r*s+v3*s^2; if((kt2%k1)==0, if(kt2 < 0, k1=-k1); tt=kt2/k1; if(issquare(tt), return([k1,sqrtint(tt),-r,s]) ) ) ); [] } \\ solve ku^2=g_0p^2-g_0dq^2 with |p|+|q|=s2 \\ return [k_0,u_0,p_0,q_0] or [] solve6(k,d,g0,s2) = { local(p); local(ku2,uu); \\ print("[k,d,g0,s2]=",[k,d,g0,s2]); for(q=1,(s2-1), p=s2-q; ku2=g0*(p^2-d*q^2); if((ku2%k)==0, if(ku2 < 0, k=-k); uu=ku2/k; if(issquare(uu),return([k,sqrtint(uu),p,q])) ) ); [] } \\ solve h^2=df^2+afg+eg^2 with |f|+|g|=s1 , g != 0 \\ return [d,h_0,f_0,g_0] or [] solve3(d,a,e,s1) = { local(f,h2); \\ print([d,a,e,s1]); for(g=1,(s1-1), f = s1-g; h2 = d*f^2+a*f*g+e*g^2; if(issquare(h2),return([d,sqrtint(h2),f,g])); h2=d*f^2-a*f*g+e*g^2; if(issquare(h2),return([d,sqrtint(h2),-f,g])) ); [] }