/* BC program powerdd */ /* outputs global variables z1 and z2 */ define powerdd(a,b,d,n){ /* (a+b*sqrt(d))^n=z1+z2*sqrt(d) */ auto x1,x2,temp,temp1,temp2,temp3,temp4,temp5,y x1=a x2=b y=n z1=1 z2=0 while(y){ while(y%2==0){ y=y/2 temp=x1 temp1=x2*x2 temp2=x1*x1 temp3=d*temp1 x1=temp2+temp3 temp4=temp*x2 x2=2*temp4 } y=y-1 temp=z1 temp1=z2*x2 temp2=z1*x1 temp3=d*temp1 z1=temp2+temp3 temp4=temp*x2 temp5=z2*x1 z2=temp4+temp5 } return } /* returns (x+ysqrt(d))(u+vsqrt(d))^n, n >= 0. * for use in depicting the equivalence class containing (x,y) when we take (u,v) to be * (x1,y1) or (x1,-y1), where (x1,y1) is the fundamental solution of Pell's equaiton. */ define hyperbolan(x,y,u,v,d,n){ auto i,j for(i=0;i<=n;i++){ j=multrootd(x,y,u,v,d) print "(",globalxrootd,",",globalyrootd,")\n" x=globalxrootd y=globalyrootd } return } /* outputs global variables z1 and z2 */ define modpowerdd(a,b,d,n,m){ /* (a+b*sqrt(d))^n=z1+z2*sqrt(d) (mod m) */ auto x1,x2,temp,temp1,temp2,temp3,temp4,temp5,y x1=a x2=b y=n z1=1 z2=0 while(y){ while(y%2==0){ y=y/2 temp=x1 temp1=mod(x2*x2,m) temp2=mod(x1*x1,m) temp3=mod(d*temp1,m) x1=mod(temp2+temp3,m) temp4=mod(temp*x2,m) x2=mod(2*temp4,m) } y=mod(y-1,m) temp=z1 temp1=mod(z2*x2,m) temp2=mod(z1*x1,m) temp3=mod(d*temp1,m) z1=mod(temp2+temp3,m) temp4=mod(temp*x2,m) temp5=mod(z2*x1,m) z2=mod(temp4+temp5,m) } return } /* returns (x+ysqrt(d))(u+vsqrt(d)/2)^n, n >= 1. * for use in depicting the equivalence class containing (x,y) when we take (u,v) to be * (x1,y1) or (x1,-y1), where (x1,y1) is the fundamental solution of Pell's equaiton. */ define hyperbolas(x,y,u,v,d,n){ auto i,j,t print "(",x,",",y,")\n" for(i=1;i<=n;i++){ j=multrootd(x,y,u,v,d) t=2^i print "(",globalxrootd/t,",",globalyrootd/t,")\n" x=globalxrootd y=globalyrootd } return } /* outputs global variables z1 and z2 */ define modpowerdd(a,b,d,n,m){ /* (a+b*sqrt(d))^n=z1+z2*sqrt(d) (mod m) */ auto x1,x2,temp,temp1,temp2,temp3,temp4,temp5,y x1=a x2=b y=n z1=1 z2=0 while(y){ while(y%2==0){ y=y/2 temp=x1 temp1=mod(x2*x2,m) temp2=mod(x1*x1,m) temp3=mod(d*temp1,m) x1=mod(temp2+temp3,m) temp4=mod(temp*x2,m) x2=mod(2*temp4,m) } y=mod(y-1,m) temp=z1 temp1=mod(z2*x2,m) temp2=mod(z1*x1,m) temp3=mod(d*temp1,m) z1=mod(temp2+temp3,m) temp4=mod(temp*x2,m) temp5=mod(z2*x1,m) z2=mod(temp4+temp5,m) } return } /* returns an SFS in the Stolt class defined by (2x,2y) in relation to solving x^2-dy^2=4m, odd m. * Here x>0 and y>0 are even and (x,y) is an NSF for x^2-dy^2=4m. * Either (x,y) is an SFS or if x1+y1\sqrt(d)=(x+ysqrt(d))(x4-y4\sqrt(d)/2), then * (x1+y1\sqrt(d))/2 is an SFS. Here (x4+y4\sqrt(d))/2 is the least solution of * x4^2-dy4^2=4, with x4 and y4 odd. Here d=8k+5. We apply this to d=5 and 13. */ define hyperbolas1(x,y,x4,y4,d,m){ auto i,j,x1,y1,s for(i=0;i<=1;i++){ s=testsfs(x,y,x4,y4,d,m) if(i==1 && s==0){ print "something is wrong here!\n" return } if(s){ print "(",x,",",y,") is an SFS\n" return } j=multrootd(x,y,x4,-y4,d) x=globalxrootd/2 y=globalyrootd/2 print "(",x,",",y,")\n" } return } define testsfs(x,y,x4,y4,d,m){ auto s,g if(m>0){ if(y^2>0 && y^2*(x4+2)=4*m && y^2*(x4-2) 0 (sign=1), but (r,s/2^k) if N < 0, (sign=-1) if 2^k divides s. */ define hyperbolas1(x,y,u,v,d,n,k,sign){ auto i,j,t,r,s,z for(i=1;i<=n;i++){ j=multrootd(x,y,u,-v,d) r=globalxrootd/2 s=globalyrootd/2 z=2^k if(s%z==0){ print "i=",i,":" if(sign==1){ print "(",-r,",",-s/z,")\n" }else{ print "(",r,",",s/z,")\n" } } x=r y=s } return } /* This takes an ambiguous SFS xi=(x,y) of x^2-dy^2=4n and computes * xi.e^i for i=0,...,l. */ define case4kambig(x,y,u,v,d,k){ auto l,l1,l2,z,j l=lpower(d,k) z=2^k if(y%z==0){ print "i=0:(",x,",",y/z,")\n" } j=hyperbolas0(x,y,u,v,d,l-1,k) } /* This first computes rn+sn.sqrt(d)=(r0+s0.sqrt(dd))(x4+-y4.sqrt(dd)/2)^i, for 1<= i <= l, * where (x4,y4) is the fundamental solution of Pell's equation u^2-ddv^2=4 * and sign=1 or -1 acoording as n>0 or n<0, where (r0,s0) is a solution of x^2-dd.y^2=4n. * Here d=4^kdd, where 4 does not divide D. Here l is the rank of apparition. 19th January 2018. * This assumes (r0,s0) is not an ambiguous solution of x^2-dd*y^2=4n. */ define hyperbolas2(r0,s0,d,sign){ auto i,j,x,y,z,x4,y4,r,s,xtemp,ytemp,l,k,rem,two2k,rr,count rr=d%4 k=0 while(rr==0){ d=d/4 rr=d%4 rr=d%4 k=k+1 } dd=d print "dd=",dd,"\n" l=lpower(dd,k) print "l=",l,"\n" print "k=",k,"\n" z=fund4(dd,0) two2k=2^k print "two2k=",two2k,"\n" x4=globalx;y4=globaly print "(x4,y4)=(",x4,",",y4,")\n" rem=s0%two2k if(rem==0){ print "(r[0],s[0])=(",r0,",",s0/two2k,")\n" } r=r0;s=s0 count=1 while(1){ j=multrootd(r,s,x4,-y4,dd) x=globalxrootd/2; y=globalyrootd/2 r=x;s=y rem=y%two2k if(rem==0 && sign==1){ print "(-r[",(-1)^count,"],-s[",(-1)^count,"])=(",-x,",",-y/two2k,")\n" } if(rem==0 && sign==-1){ print "(r[",count*(-1)^count,"],s[",count*(-1)^count,"])=(",x,",",y/two2k,")\n" } if(count==1){ xtemp=r0;ytemp=s0 } count=count+1 if(count==l){ return } j=multrootd(xtemp,ytemp,x4,y4,dd) xtemp=globalxrootd/2 ytemp=globalyrootd/2 rem=ytemp%two2k if(rem==0){ print "(r[",count*(-1)^count,"],s[",count*(-1)^count,"])=(",xtemp,",",ytemp/two2k,")\n" } count=count+1 if(count==l){ return } } } /* This first computes rn+sn.sqrt(d)=(r0+s0.sqrt(dd))(x4+-y4.sqrt(dd)/2)^i, for 1<= i <= l, * where (x4,y4) is the fundamental solution of Pell's equation u^2-ddv^2=4 * and sign=1 or -1 acoording as n>0 or n<0, where (r0,s0) is a solution of x^2-dd.y^2=4n. * Here d=4^kdd, where 4 does not divide D. Here l is the rank of apparition. 19th January 2018. * This assumes (r0,s0) is an ambiguous solution of x^2-dd*y^2=4n. */ define hyperbolas22(r0,s0,d,sign){ auto i,j,x,y,z,x4,y4,r,s,xtemp,ytemp,l,k,rem,two2k,rr,count rr=d%4 k=0 while(rr==0){ d=d/4 rr=d%4 rr=d%4 k=k+1 } dd=d print "dd=",dd,"\n" l=lpower(dd,k) print "l=",l,"\n" print "k=",k,"\n" z=fund4(dd,0) two2k=2^k print "two2k=",two2k,"\n" x4=globalx;y4=globaly print "(x4,y4)=(",x4,",",y4,")\n" rem=s0%two2k if(rem==0){ print "(r[0],s[0])=(",r0,",",s0/two2k,")\n" } r=r0;s=s0 count=1 while(1){ if(count==1){ xtemp=r0;ytemp=s0 } if(count==l){ return } j=multrootd(xtemp,ytemp,x4,y4,dd) xtemp=globalxrootd/2 ytemp=globalyrootd/2 rem=ytemp%two2k if(rem==0){ print "(r[",count*(-1)^count,"],s[",count*(-1)^count,"])=(",xtemp,",",ytemp/two2k,")\n" } count=count+1 if(count==l){ return } } } /* This finds the fundamental solution (globalx,globaly) of Pell's equation x^2-dy^2=4, * where d>0 and is nonsquare. * We use the fact that a positive solution with gcd(x,y)=1 must be a convergent of * the continued fraction expansion of sqrt(d). If there is no such convergent, we use th * familiar midpoint method of finding the least solution of Pell's equations X^2-dY^2=1 * and then (globalx,globaly)=(2X,2Y). */ define fund4(d,printflag){ auto h,p,q,e,f,g,x,y,k,l,m,n,u,v,r,s,oldl,oldm,oldu,oldv,t,flag flag=0 if(d==2){ globalx=6 globaly=4 flag=1 } if(d==3){ globalx=4 globaly=2 flag=1 } if(d==5){ globalx=3 globaly=1 flag=1 } if(d==6){ globalx=10 globaly=4 flag=1 } if(d==7){ globalx=16 globaly=6 flag=1 } if(d==8){ globalx=6 globaly=2 flag=1 } if(d==10){ globalx=38 globaly=12 flag=1 } if(d==11){ globalx=20 globaly=6 flag=1 } if(d==12){ globalx=4 globaly=1 flag=1 } if(d==13){ globalx=11 globaly=3 flag=1 } if(d==14){ globalx=30 globaly=8 flag=1 } if(d==15){ globalx=8 globaly=2 flag=1 } x=sqrt(d) if(d==x^2+1 && d >17){ globalx=4*x^2+2 globaly=4*x flag=1 } if(printflag && flag){ print "globalx=",globalx,",globaly=",globaly,"\n" } if(flag){ return(d) } x=sqrt(d) p=0;q=1 g=x*x l=0;k=1;m=1;n=0 for(h=0;1;h++){ y=(x+p)/q u=k*y+l;v=n*y+m /*print "n=",n,",y=",y,",m=",m,",v=",v,"\n"*/ if(printflag){ print "A[",h,"]/B[",h,"]=",u,"/",v,"\n" } oldl=l;oldm=m oldu=u;oldv=v; l=k;m=n k=u;n=v f=p p=y*q-p e=q q=(d-(p*p))/q t=(-1)^(h+1) if(q==4 && t==1){ globalx=k globaly=n if(printflag){ print "globalx=",globalx,",globaly=",globaly,"\n" } return(4) } if(q==4 && t!=1){ globalx=(k^2+d*n^2)/2 globaly=k*n if(printflag){ print "globalx=",globalx,",globaly=",globaly,"\n" } return(-4) } if(p==f){/* P_h=P_{h+1}, even period 2h */ if(printflag){ print "P_",h,"=P_",h+1,"\n" } globalx=2*(m*(oldu+oldl)+(-1)^h) globaly=2*m*(oldv+oldm) if(printflag){ print "globalx=",globalx,",globaly=",globaly,"\n" } return(1) } if(q==e){/* Q_h=Q_{h+1}, odd period 2h+1 */ if(printflag){ print "Q_",h,"=Q_",h+1,"\n" } r=k*n+l*m s=n^2+m^2 globalx=2*(r*r+d*s*s) globaly=4*r*s if(printflag){ print "globalx=",globalx,",globaly=",globaly,"\n" } return(-1) } } } define lpower(d,k){ auto l,x1,x2,y1,y2,x,z x=fund4(d,0) x1=globalx y1=globaly x=fund4((4^k)*d,0) x2=globalx y2=globaly l=1 while(1){ z=powerdd(x1,y1,d,l) if(z1==2^(l-1)*x2){ return(l) } l=l+1 } }