\\ \\ Gentleman's Diary Problem \\ \\ x,y: positive integers \\ such that x+y = u^2, x^2+y^2 = v^2, x^3+y^3 = w^2 \\ \\ t=w/u ==> x^2-xy+y^2 = t^2 \\ \\ X=(-3)v^2/t^2 \\ Y=3(x^2-y^2)v/t^3 \\ \\ E: Y^2=x^3+8X^2+12X \\ \\ E_tors(Q)=Z/2ZxZ/2Z generator [0,0],[-2,0] \\ ={[0,0],[-2,0],[-6,0],O} \\ rank(E(Q))=1 generator [-3,-3] \\ E(Q)=Z+Z/2Z+Z/2Z \\ default(realprecision,1000); \\ \\ [x,y]-->[u,v,w] \\ gg(p)= { local(x,y,z,d); x=p[1]; y=p[2]; z=lcm(denominator(x),denominator(y)); x=x*z; y=y*z; d=gcd(gcd(x,y),z); \\ print("# d=",d,";"); if(z<0, d=-d); [x/d,y/d,z/d] } issquareint(x)= { (x >= 0 && denominator(x) == 1 && sqrtint(x)^2 == x); } gentleman(p,r)= { local(q,v,v2,w,tt,xx,yy,zz,x,y,z); \\ print("# p=",p,";"); q=gg(p); \\ print("# q=",q,";"); xx=q[1]; yy=q[2]; zz=q[3]; \\ print("# gcd(X,Y)=",gcd(xx,yy),";"); tt=floor(zz^(1/3)+0.01); if(tt^3 == zz, v2=xx/((-3)*tt); if(v2 > 0, \\ print("# tt=",tt,";"); if(issquareint(v2), v = sqrtint(v2); \\ print("# v=",v,";"); w = yy/((-3)*v); if(denominator(w) == 1, \\ print("# w=",w,";"); if(w > 0, x=(v2+w)/2; y=(v2-w)/2, x=(v2-w)/2; y=(v2+w)/2 ); if(issquareint(x) && issquareint(y), x=sqrtint(x); y=sqrtint(y); print("# q=",q,"=",r[1],r[2],"+",r[3],";"); \\ print("# [x,y]=",[x,y],";"); if(v2>=tt^2, print("[x,y]=",[x*(x+y),y*(x+y)]), if(y>0, print("[x,y]=",[x*(x-y),-y*(x-y)]) ) ) ) ) ), tt = -tt; v2 = -v2; \\ print("# tt=",tt,";"); if(issquareint(v2), v = sqrtint(v2); \\ print("# v=",v,";"); w = -yy/((-3)*v); if(denominator(w) == 1, \\ print("# w=",w,";"); if(w > 0, x=(v2+w)/2; y=(v2-w)/2, x=(v2-w)/2; y=(v2+w)/2 ); if(issquareint(x) && issquareint(y), x=sqrtint(x); y=sqrtint(y); print("# q=",q,"=",r[1],r[2],"+",r[3],";"); \\ print("# [x,y]=",[x,y],";"); if(v2>=tt^2, print("[x,y]=",[x*(x+y),y*(x+y)]), if(y>0, print("[x,y]=",[x*(x-y),-y*(x-y)]) ) ) ) ) ) ) ) } gdiary(n)= { local(e,p,q,g); e=ellinit([0,8,0,12,0]); \\ p1=[0,0]; \\ p2=[-2,0]; \\ p3=[-6,0]; g=[-4,-4]; for(i=1,n, p=ellpow(e,g,i); \\ print("i=",i,":"); \\ print(p); gentleman(p,[i,g,[0]]); q=elladd(e,p,[0,0]); \\ print(q); gentleman(q,[i,g,[0,0]]); q=elladd(e,p,[-2,0]); \\ print(q); gentleman(q,[i,g,[-2,0]]); q=elladd(e,p,[-6,0]); \\ print(q); gentleman(q,[i,g,[-6,0]]); ) } check(p)= { local(x,y,u2,v2,w2,u,v,w); x=p[1]; y=p[2]; u2=x+y; v2=x^2+y^2; w2=x^3+y^3; \\ print([u2,v2,w2]); u=sqrtint(u2); v=sqrtint(v2); w=sqrtint(w2); \\ print([u,v,w]); if(u^2==u2 && v^2==v2 && w^2==w2, [u,v,w], [] ) }