/* gnubc program : patz */ /* Needs squareroot */ /* to run sswgeneral, one needs * squareroot2 powertest3 genfacs gcd posformrep reduceneg congruence powerdd */ /* Typing pn(a[],n) returns p[n], where * p[0]=a[0],p[1]=a[0]*a[1]+1,p[i+1]=a[i+1]*p[i]+p[i-1] if i >= 1. * Typing qn(a[],n) returns q[n], where * q[0]=1,q[1]=a[1],q[i+1]=a[i+1]*q[i]+q[i-1] if i >= 1. * Hence pn/qn is the value of the simple continued fraction * [a[0];a[1],...,a[n]]. */ define pn(a[],n){ auto x,y,z,i if (n == 0)return(a[0]) if(n==1)return(a[0]*a[1]+1) x=a[0];y=a[0]*a[1]+1 for(i=2;i<=n;i++){ z=a[i]*y+x x=y y=z } return(z) } define qn(a[],n){ auto x,y,z,i if(n==0)return(1) if(n==1)return(a[1]) x=1;y=a[1] for(i=2;i<=n;i++){ z=a[i]*y+x x=y y=z } return(z) } /* * This program finds the continued fraction expansion of a * quadratic irrational (u+tsqrt(d))/v, where d is not a square, * t,u,v integers, v nonzero. * After converting to (U+sqrt(D))/V, where V divides D-U^2, */ /* sign of an integer a */ /* sign(a)=1,-1,0, according as a>0,a<0,a=0 */ define sign(a){ if(a>0) return(1) if(a<0) return(-1) return(0) } /* absolute value of an integer n */ define abs(n){ if(n>=0) return(n) return(-n) } /* mod(a,b)=the least non-negative remainder when an integer a is divided by a positive integer b */ define mod(a,b){ auto c c=a%b if(a>=0) return(c) if(c==0) return(0) return(c+b) } /* int(a,b)=integer part of a/b, a, b integers, b != 0 */ define int(a,b){ auto c c=sign(b) a=a*c b=b*c return((a-mod(a,b))/b) } /* * This is a function for finding the period of the continued fraction * expansion of reduced quadratic irrational a=(u+sqrt(d))/v. * Here d is non-square, 1<(u+sqrt(d))/v, -1<(u-sqrt(d))/v<0. * The algorithm also assumes that v divides d-u*u and is based on K. Rosen, * Elementary Number theory and its applications, p.379-381 and Knuth's The art * of computer programming, Vol. 2, p. 359. 0 is returned if a is not reduced. */ define period(d,u,v,i){ auto r,s,t,f,j,k,temp f=sqrt(d) if(u<=0){"a is not reduced: u<=0:";return(0)} if(v<=0){"a is not reduced: v<=0:";return(0)} if(u+f0:";return(0)} k=d-u*u if(k%v!=0){"v does not divide d - u*u:";return(0)} s=v r=u a_surd_length=i /* i is created by surd(d,t,u,v) below and indexes the ith convergent of (u+t*sqrt(d))/v */ for(j=i;1;j++){ a=(f+u)/v /* print "a[",j,"]=",a,"\n" */ a_surd[j]=a u=a*v-u temp=v v=(d-u*u)/v u_surd[j+1]=u v_surd[j+1]=v if(u==r && v==s){ t=j+1-i a_surd_length=a_surd_length+t return(t) /* the continued fraction has period length j+1-i */ } } } /* * This function uses the continued fraction algorithm expansion in K. Rosen, * Elementary Number theory and its applications,p.379-381 and Knuth's * The art of computer programming, Vol. 2, p. 359. It locates the first reduced * complete quotient and then uses the function period(d,u,v,i,flag) * to locate the period of the continued fraction expansion. */ define surd(d,t,u,v,print_flag){ auto a,c,e,g,f,i,j,k,w,z,con_length c=d;g=u;e=v if(v==0){ print "v=" return(0) } f=sqrt(d) if(f*f==d){ print "d is the square of " return(f) } if(t==0){ print "t = " return(0) } z=sign(t) u=u*z v=v*z d=d*t*t w=abs(v) if((d-u*u)%w!=0){ d=d*v*v u=u*w v=v*w } f=sqrt(d) for(j=0;1;j++){ u_surd[j]=u v_surd[j]=v g=reducetest(u,v,d) if(g){/* start if(g) */ /* (u+sqrt(d))/v is reduced */ i=j pre_period_length=i l=period(d,u,v,i) if(l%2){ l=period(d,u,v,j+l) if(print_flag){ print "printing first and second period:\n" } } con_length=a_surd_length-1 for(k=0;k<=con_length;k++){ p_surd[k]=pn(a_surd[],k) q_surd[k]=qn(a_surd[],k) } if(print_flag){ for(k=0;k<=con_length;k++){ } print "\n" for(k=0;k<=con_length;k++){ if(k==pre_period_length || l%2 && k==pre_period_length+l){ print "\n period:\n" } print "(P[",k,"]+sqrt(",d,"))/Q[",k,"]=(",u_surd[k],"+sqrt(",d,"))/",v_surd[k],", " print "A[",k,"]/B[",k,"]=",p_surd[k],"/",q_surd[k],", " print "a[",k,"]=",a_surd[k],"\n" } print "cfrac has period length ",l,"\n" } return(l) }/* if(g) */ if(v>0){ a=int(f+u,v) } if(v<0){ a=int(f+u+1,v) } /*print "a[",j,"]=",a,"\n" */ a_surd[j]=a u=a*v-u v=int(d-u*u,v) } } /* This finds fundamental solutions (globalfundx[h], globalfundy[h]) for the diophantine equation * x^2-dy^2=n, d>0, d not a perfect square n nonzero and not equal to 1. * Returns the number of solution classes. */ define patz(d,n,print_flag){ auto a,b,k,q0,s,t,period_length,temp,g1_surd,g2_surd,f1,f2,solution_number,pos,pos1,pos2,flagpos,g solution_number=0 if(d<=1){ print "d<=1\n" return(-1) } g=sqrt(d) if(g^2==d){ print "d is a perfect square\n" return(-1) } if(n==1){ globalfundx[0]=1 globalfundy[0]=0 return(1) } q0=abs(n) s=sqroot(d,q0,0) /* s=0 means no solutions of x^2=d (mod |n|) */ /* If s>0, we get all solutions x, 0<=x<=|n|/2 in the form * solution[0],...,solution[numbr-1] */ if(s==0){ if(print_flag){ print "x^2=",d," mod(",n,") has no solution\n" } return(solution_number) } /* Now to test the cfrac of each of w[j]=(-solution[j]+\sqrt{D})/q0 to see if v_surd[j]=(-1)^j*sign(n) holds for some j in 1<=j<=t+l, where l is the period of w[j]=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]] and t=pre_period_length. */ for(k=0;kq0){ m=m-q00 } return(m) } /* binary0(a,b,c,n,print_flag) tests the diophantine equation * ax^2+bxy+cy^2=n for solubility. Here d=b^2-4ac>0 and is not a perfect * square. Also n is non-zero, gcd(a,b,c)=1 and gcd(a,n)=1. * In the case of solubility, the solution (x,y) in each class, * with least y is obtained. * With print_flag =1, we get verbose output. * See K.R. Matthews, 'The Diophantine equation ax^2+bxy+cy^2=N, * D=b^2-4ac>0', J. de The'orie des Nombres de Bordeaux 14 (2002) 257-270. * The test for necessity is simplified by the observation that * if Q_n=2(-1)^nN/|N| occurs in the cfrac of omega=(-n+sqrt(d))/q, * then Q_m=2(-1)^{m+1}N/|N| occurs in the cfrac of omega*=(-n-sqrt(d))/q. * creates global fund_x[] and fund_y[]. */ define binary0(a,b,c,n,print_flag){ auto no_of_solutions, d,k,m,mm,q,q0,q00,s,t,temp,x,y,ee,qdash,ndash,rdash,theta,result if(print_flag){ print "Solving " result=printform(a,b,c,n) } if(a==n){ fund_x[0]=1 fund_y[0]=0 mm=calculate_n(a,b,c,1,0,n) if(print_flag){ print "class n=",mm," (mod 2)\n" print "solution (1,0)\n" } return(1) } d=b^2-4*a*c q0=abs(n) q00=2*q0 q=2*a*q0 s=quadratic(a,b,c,q0,0) /* s=0 means no solutions of ax^2+bx+c=0 (mod |n|) */ /* If s>0, we get all solutions x, 0<=x<|n| in the form * quadratic_solution[0],...,quadratic_solution[s-1] */ if(s==0){ if(print_flag){ print "ax^2+bx+c=0 mod(n) has no solution\n" } return(0) } /* With m=2*a*quadratic_solution[k]+b, to test the cfrac of each of w[j]=(-m+\sqrt{d})/q to see if v_surd[j]=(-1)^j*2**sign(n) holds for some j in 1<=j<=t+l, where l is the period of w[j]=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]] and t=pre_period_length. */ no_of_solutions=0 for(k=0;k1. Again we assume d=b^2-4ac>0 * and that d is not a square. Also we assume gcd(a,b,c)=1. * The number r of solution classes is returned. * verbose output if print_flag=1. * Printing of solution classes is suppressed if gcd(a,n)>1. */ define binary1(a,b,c,n,print_flag){ auto aa,bb,cc,i,q00,mm,r,s,tmp1,tmp2,tmp3,tmp4,result if(print_flag){ print "Solving " result=printform(a,b,c,n) } if(gcd(a,n)==1){ r=binary0(a,b,c,n,print_flag) return(r) } /* Now we calculate aa,bb,cc, where ax^2+bxy+cy^2 is * transformed into aaX^2+bbXY+ccY^2 under the transformation * x=alpha*X+beta*Y, y=gamma*X+delta*Y. */ q00=2*abs(n) s=gauss(a,b,c,n) aa=gauss_m gauss_delta=gcd1(gauss_alpha,gauss_gamma) gauss_beta= -gcd2(gauss_alpha,gauss_gamma) tmp1=a*gauss_beta*gauss_beta tmp2=b*gauss_beta*gauss_delta tmp3=c*gauss_delta*gauss_delta tmp4=tmp1+tmp2 cc=tmp4+tmp3 tmp1=2*a*gauss_alpha*gauss_beta tmp2=2*c*gauss_gamma*gauss_delta tmp3=b*gauss_alpha*gauss_delta tmp4=b*gauss_beta*gauss_gamma bb=tmp1+tmp2+tmp3+tmp4 r=binary0(aa,bb,cc,n,print_flag) /* replace 0 by 1 to print fundamental solutions for reduced * case gcd(aa,n)=1. */ if(r==0){ return(0) } /* this creates arrays of fundamental solutions: * fund_x[0],...,fund_x[r-1] and * fund_y[0],...,fund_y[r-1]. */ for(i=0;i1 and produces * (x,y)=(gauss_alpha,gauss_gamma) and an m=gauss_m such that ax^2+bxy+cy^2=m, * gcd(m,n)=1. See L.-K. Hua 'Introduction to number theory', 311-312. */ define gauss(a,b,c,n){ auto absn,e,g,i,s,t,tmp1,tmp2,xx[],yy[] absn=abs(n) e=omega(absn) for(i=0;i0 and nonsquare, n nonzero, gcd(a,b,c)=1. * print_flag=1 gives verbose output. */ define binary(a,b,c,n,print_flag){ auto s,i,g,absn,p,u,v,delta,d,result g=gcd3(a,b,c) absn=abs(n) if(n%absn){ print "gcd(a,b,c) does not divide n\n" return(0) }else{ a=a/g b=b/g c=c/g n=n/g } d=b^2-4*a*c if(print_flag){ print "d = ",d,"\n" } if(b%2==0){ p=b/2 delta=p^2-a*c if(print_flag){ print "delta = ",delta,"\n" } } s=binary1(a,b,c,n,print_flag) if(print_flag){ print "The equation " result=printform(a,b,c,n) if(s==0){ print "has no primitive solution\n" }else{ if(s==1){ print "has one primitive solution class: " }else{ print "has ", s," primitive solution classes:\n"; } for(i=0;i0, we get all solutions x, 0<=x<=|n|/2 in the form * solution[0],...,solution[numbr-1] */ if(s==0){ return(solution_number) } /* Now to test the cfrac of each of w[j]=(-solution[j]+\sqrt{D})/q0 to see if v_surd[j]=(-1)^j*sign(n) holds for some j in 1<=j<=t+l, where l is the period of w[j]=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]] and t=pre_period_length. */ for(k=0;k4 && flag==0){ print "t=",i,",d=",d, ",N=",s,",sn=",sn,":\n" } } } define selin(t){ return(2*t^3-3*t^2-6*t+9) } define dujella_test(m,n,flag){ auto k,d,sn,temp for(k=m;k<=n;k++){ temp=k^2 d=temp+1 if(flag){ print "k=",k,":" } sn=patz0(d,temp,flag) if(sn>2 && flag==0){ print "k=",k,",d=",d, ",N=",temp,",sn=",sn,":\n" } } } /* define dujella_test1(m,n){ auto t,d,sn,temp for(t=m;t<=n;t++){ temp=f(2*t,2*t) d=temp^2+1 print "t=",t,":" sn=patz0(d,temp^2) print "sn=",sn,":\n" } }*/ define conjecture(d,n,s){ auto p,i,q[] i=0 p=s q[0]=p n[0]=n print "q[0]= ",s,"," print "n[0]=",n,"," while(n*n>d){ i=i+1 n=(p^2-d)/n n[i]=n print "n[",i,"]=",n,"," p=p%n q[i]=p print "q[",i,"]=",p,"," } if(n==1){ print "reached n[",i,"]=1\n" }else{ print "reached n[",i,"]=",n,"\n" } } define dujella_test2(k){ auto t,d,n,i,r d=k^2+1 n=d-1 print "k=",k,",d=k^2=",d,":" print "Solving congruence P^2=",1," mod(k^2)\n" print "1<=P<=n/2=",n/2,":\n" t=sqroot(1,n,0) t=t/2 print "We get ",t," solutions: " for(i=0;i=2 && f>1)||(sn>=4 && f==1)){ sn=patz0(kk,(k/f)^2,1) print "sn:",sn,"," print "k:",k,"," print "f:",f,"\n" } } } } return } define dujella_test4(k){ auto d,div2k,t,s,kk,divisor2k[],divkk,ss,twok,n,a,b,g,flag,p,q kk=k^2+1 twok=2*k t=omega(twok) div2k=divisors(qglobal[],kglobal[],t) for(s=0;s=b){ continue } if(gcd(a,b)==1){ /* g=binary1(a,0,-b,n,0)*/ /*print "processing a=",a,",d=",d,"\n"*/ g=rootdovera(a,k,d) if(g){ print "k=",k,"," print "d=",d,"," print "a=",a,"," print "b=",b,"," print "\n" print "=================\n" flag=1 } } } } return(flag) } /* for use when d=k^2+1 > 5 and a divides d. */ define rootda(d,a){ auto h,p,q,e,f,g,x,y,k,l,m,n,u,v,r,s x=sqrt(d) p=0;q=a g=x*x l=0;k=1;m=1;n=0 for(h=0;1;h++){ y=(x+p)/q print y,"," u=k*y+l;v=n*y+m /* u/v is the i-th convergent to sqrt(d) */ if(h>=0){ print "A[",h,"]/B[",h,"]=",u,"/",v,"\n" } l=k;m=n k=u;n=v f=p p=y*q-p e=q q=(d-(p*p))/q if(q==e){/* Q_h=Q_{h+1}, odd period 2h+1 */ print "Q_",h,"=Q_",h+1,"\n" return(2*h+1) } } } # a divides dd=kk^2+1, ab=kk^2+1,20){ if(flag){ for(i=0;i<=h;i++){ print "a[",i,"]=",a[i],",P[",i,"]=",pcap[i],",Q[",i,"]=",qcap[i],",A[",i,"]/B[",i,"]=",pu[i],"/",qv[i],"\n" } } return(flag) } l=k;m=n k=u;n=v p=y*q-p e=q q=(d-(p*p))/q if(dd==1){ t1=abs(e-q+2*p) if(t1==2*kk){ temp=u+l temq=v+m print "|Q[",h,"]-Q[",h+1,"]+2P[",h+1,"]|=",2*k,"=2k, A[",h,"]+A[",h-1,"]/B[",h,"]+B[",h-1,"]=",temp,"/",temq,"=p/q\n" print "solves ap^2-bq^2=pm2k/d=",a*(temp)^2-b*(temq)^2,"\n" print "x=",dd*(a*(temp)^2+b*(temq)^2)/2,", y=",dd*(temp)*(temq),"\n" flag=1 } t2=abs(e-q-2*p) if(t2==2*kk){ temp=u-l temq=v-m print "|Q[",h,"]-Q[",h+1,"]-2P[",h+1,"]|=",2*kk,"=2k, A[",h,"]-A[",h-1,"]/B[",h,"]-B[",h-1,"]=",temp,"/",temq,"=p/q\n" print "solves ap^2-bq^2=pm2k/d=",a*(u-l)^2-b*(v-m)^2,"\n" print "x=",dd*(a*(u-l)^2+b*(v-m)^2)/2,", y=",dd*(u-l)*(v-m),"\n" flag=1 } } if(q==twokoverd){ print "Q[",h+1,"]=",q,"=2k/d and A[",h,"]/B[",h,"]=",u,"/",v,"=p/q\n" print "solves ap^2-bq^2=pm2k/d=",a*u^2-b*v^2,"\n" print "x=",dd*(a*u^2+b*v^2)/2,", y=",dd*u*v,"\n" flag=1 } } } define dujella_cfrac_test(m,n){ auto k,t for(k=m;k<=n;k++){ t=dujella_test4(k) if(t==0){ continue } } } /* This lists the primitive fundamental solutions of x^2 - dy^2=n. */ define patztest(d,n){ auto h,g g=patz(d,n,0) print "Fundamental solutions:\n" for(h=0;h0, a!=0,n!=0, gcd(a,b,c)=1. */ define binarygen(a,b,c,n,print_flag){ auto absn,t,g,countbinarygen,temp,h,f,result g=gcd3(a,b,c) if(g>1){ print "gcd(a,b,c)>1\n" return(-1) } if(b^2-4*a*c<=0){ print "b^2-4ac<=0\n" return(-1) } absn=abs(n) if(absn!=1){ t=omega(absn) g=divisors(qglobal[],kglobal[],t) }else{ g=1 divisor[0]=1 } countbinarygen=0 for(f=0;f0){ for(h=0;h 0, a!=0, n!=0. */ define binarygenlist(a,b,c,n,print_flag){ auto h,g,d,bmod2,bigdelta,phi,psi,null # global globalx,globaly,globalr,globals g=binarygen(a,b,c,n,print_flag) if(g==-1){ return(-1) } if(g==0){ if(print_flag){ print "there are no solutions for ",a,"x^2+(",b,")xy+(",c,")y^2=",n,":\n" } return(0) } d=b^2-4*a*c bmod2=b%2 if(bmod2==0){ bigdelta=d/4 print "(b^2-4ac)/4=",bigdelta,"\n" }else{ print "b^2-4ac=",d,"\n" } if(bmod2){ null=fund4(d,0) phi=globalx;psi=globaly print "least positive solution (",phi,",",psi,") of Pell's equation\n" print "phi^2 - ",d,"psi^2=4\n"; }else{ null=fund1(bigdelta) phi=globalr;psi=globals print "least positive solution (",phi,",",psi,") of Pell's equation\n" print "phi^2 - ",bigdelta,"psi^2=1\n"; } if(print_flag){ if(g==1){ print "one solution class for ",a,"x^2+(",b,")xy+(",c,")y^2=",n,":\n" }else{ print g," solution classes for ",a,"x^2+(",b,")xy+(",c,")y^2=",n,"\n" } } for(h=0;h0 and is nonsquare. * We solve ax^2+bxy+cy^2=1. */ define generalizedpell1(a,b,c,print_flag){ auto d,i,k,q,period_length,temp,g1_surd,g2_surd,f1,f2,u,v,x,y,ss,delta,p,returnflag,t,result if(print_flag){ print "Solving " t=printform(a,b,c,1) } returnflag=0 if(a==1){ generalizedpellfundx=1 generalizedpellfundy=0 print "solution generalized pell equation " result=printform(a,b,c,1) if(b%2==0){ p=b/2 delta=p^2-a*c print "delta = ",delta,"\n" print "x = Phi - (",p,")Psi\n" print "y = Psi,\n" print "where Phi^2 - ",delta,"Psi^2 = 1\n" }else{ print "x = (Phi - (",b,")Psi)/2\n" print "y = Psi,\n" print "where Phi^2 - ",d,"Psi^2 = 4\n" } return(1) } d=b^2-4*a*c q=2*a /* to test the cfrac of w=(-b+\sqrt(d))/q to see if v_surd[j]=2(-1)^j holds for some j in 1<=j<=t+l, where l is the period of w=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]] and t=pre_period_length. also to test the cfrac of w*=(-b-\sqrt(d))/q to see if v_surd[j]=2(-1)^(j+1) holds for some j in 1<=j<=t+l, where l is the period of w*=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]] and t=pre_period_length. */ if(print_flag){ print "processing omega=(",-b,"+sqrt(",d,")/",q,"\n" } period_length=surd(d,1,-b,q,print_flag) temp=pre_period_length+period_length if(period_length%2){ temp=temp+period_length } for(i=1;i<=temp;i++){ temp1=(-1)^i*2 if(v_surd[i]==temp1){ returnflag=1 solution_flag=1 f1=q_surd[i-1] g1_surd=p_surd[i-1] if(print_flag){ print "solution (x,y)=(",g1_surd,",",f1,")\n" print "of generalized Pell equation " result=printform(a,b,c,1) } break } } if(print_flag){ print "processing omega*=(",-b,"-sqrt(",d,")/",q,"\n" } period_length=surd(d,-1,-b,q,print_flag) temp=pre_period_length+period_length if(period_length%2){ temp=temp+period_length } for(i=1;i<=temp;i++){ temp1=(-1)^(i+1)*2 if(v_surd[i]==temp1){ f2=q_surd[i-1] g2_surd=p_surd[i-1] if(print_flag){ print "solution (x,y)=(",g2_surd,",",f2,")\n" print "of generalized Pell equation " result=printform(a,b,c,1) } break } } if(d!=5 || (d==5 && a>0)){ if(returnflag==0){ return(0) } } if(d!=5 || (d==5 && a>0)){ if(f1<=f2){ u=g1_surd v=f1 }else{ u=g2_surd v=f2 } generalizedpellfundx=u generalizedpellfundy=v if(print_flag){ print "solution with least positive y (", u,",",v,") " print "of generalized Pell equation " result=printform(a,b,c,1) } } if(d==5 && a<0){ ss=pre_period_length if(ss>1){ tmp1=p_surd[ss-1]-p_surd[ss-2] tmp2=q_surd[ss-1]-q_surd[ss-2] }else{ tmp1=p_surd[0]-1 tmp2=q_surd[0] } returnflag=1 generalizedpellfundx=tmp1 generalizedpellfundy=tmp2 u=tmp1 v=tmp2 if(print_flag){ print "exceptional solution (", tmp1,",",tmp2,") with least positive y\n" print "of generalized Pell equation " result=printform(a,b,c,1) } } if(print_flag){ print "general solution of generalized pell equation:\n" if(b%2==0){ p=b/2 delta=p^2-a*c print "delta = ",delta,"\n" print "x=",u,"Phi - (",c*v+p*u,")Psi\n" print "y=",v,"Phi + (",p*v+a*u,")Psi\n" print "Phi^2 - ",delta,"Psi^2 = 1\n" }else{ print "d = ",d,"\n" print "x=",u,"Phi/2 - (",2*c*v+b*u,")Psi/2\n" print "y=",v,"Phi/2 + (",b*v+2*a*u,")Psi/2\n" print "Phi^2 - ",d,"Psi^2 = 4\n" } } return(returnflag) } /* Here a is nonzero, d=b^2-4ac>0 and is nonsquare. * We solve ax^2+bxy+cy^2=1. */ define generalizedpell(a,b,c,print_flag){ auto d,i,k,q,period_length,temp,g1_surd,g2_surd,f1,f2,u,v,x,y,ss,delta,p,returnflag,t,result if(print_flag){ print "Solving " t=printform(a,b,c,1) } returnflag=0 if(a==1){ generalizedpellfundx=1 generalizedpellfundy=0 print "solution generalized pell equation " result=printform(a,b,c,1) print "with smallest nonnegative y: (1,0)\n" return(1) } d=b^2-4*a*c q=2*a /* to test the cfrac of w=(-b+\sqrt(d))/q to see if v_surd[j]=2(-1)^j holds for some j in 1<=j<=t+l, where l is the period of w=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]] and t=pre_period_length. also to test the cfrac of w*=(-b-\sqrt(d))/q to see if v_surd[j]=2(-1)^(j+1) holds for some j in 1<=j<=t+l, where l is the period of w*=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]] and t=pre_period_length. */ if(print_flag){ print "processing omega=(",-b,"+sqrt(",d,")/",q,"\n" } period_length=surd(d,1,-b,q,print_flag) temp=pre_period_length+period_length if(period_length%2){ temp=temp+period_length } for(i=1;i<=temp;i++){ temp1=(-1)^i*2 if(v_surd[i]==temp1){ returnflag=1 solution_flag=1 f1=q_surd[i-1] g1_surd=p_surd[i-1] if(print_flag){ print "solution (x,y)=(",g1_surd,",",f1,")\n" print "of generalized Pell equation " result=printform(a,b,c,1) } break } } if(print_flag){ print "processing omega*=(",-b,"-sqrt(",d,")/",q,"\n" } period_length=surd(d,-1,-b,q,print_flag) temp=pre_period_length+period_length if(period_length%2){ temp=temp+period_length } for(i=1;i<=temp;i++){ temp1=(-1)^(i+1)*2 if(v_surd[i]==temp1){ f2=q_surd[i-1] g2_surd=p_surd[i-1] if(print_flag){ print "solution (x,y)=(",g2_surd,",",f2,")\n" print "of generalized Pell equation " result=printform(a,b,c,1) } break } } if(d!=5 || (d==5 && a>0)){ if(returnflag==0){ return(0) } } if(d!=5 || (d==5 && a>0)){ if(f1<=f2){ u=g1_surd v=f1 }else{ u=g2_surd v=f2 } generalizedpellfundx=u generalizedpellfundy=v if(print_flag){ print "solution with least positive y (", u,",",v,") " print "of generalized Pell equation " result=printform(a,b,c,1) } } if(d==5 && a<0){ ss=pre_period_length if(ss>1){ tmp1=p_surd[ss-1]-p_surd[ss-2] tmp2=q_surd[ss-1]-q_surd[ss-2] }else{ tmp1=p_surd[0]-1 tmp2=q_surd[0] } returnflag=1 generalizedpellfundx=tmp1 generalizedpellfundy=tmp2 u=tmp1 v=tmp2 if(print_flag){ print "exceptional solution (", tmp1,",",tmp2,") with least positive y\n" print "of generalized Pell equation " result=printform(a,b,c,1) } } if(print_flag){ print "general solution of generalized pell equation:\n" if(b%2==0){ p=b/2 delta=p^2-a*c print "delta = ",delta,"\n" print "x=",u,"Phi - (",c*v+p*u,")Psi\n" print "y=",v,"Phi + (",p*v+a*u,")Psi\n" print "Phi^2 - ",delta,"Psi^2 = 1\n" }else{ print "d = ",d,"\n" print "x=",u,"Phi/2 - (",2*c*v+b*u,")Psi/2\n" print "y=",v,"Phi/2 + (",b*v+2*a*u,")Psi/2\n" print "Phi^2 - ",d,"Psi^2 = 4\n" } } return(returnflag) } define printform(a,b,c,n){ print a,"x^2+(",b,")xy+(",c,")y^2=",n,"\n" } #define ssw2(a,b,c,d,e,f){ #auto gcount,h,dd,u,v,k,delta0,delta,gamma,zeta,epsilon,alpha,beta,psi,phi,p,aa,m,temp,temp1,x1,y1,null,twodelta0 # if(a==0){ # print "a=0\n" # return(-1) # } # gcount=0 # dd=b^2-4*a*c # if(dd<0){ # print "D<0\n" # return(-1) # } # if(dd==0){ # print "D=0\n" # return(-1) # } # print "D=",dd,"\n" # if(b%2){ # u=b*e-2*c*d # v=b*d-2*a*e # }else{ # u=(b/2)*e-c*d # v=(b/2)*d-a*e # } # alpha=-u # beta=-v # if(b%2==0){ # p=b/2 # delta0=p^2-a*c # twodelta0=2*delta0 # print "delta0 = ",delta0,"\n" # null=fund1(delta0) # phi=globalr # psi=globals # print "phi^2 - ",delta0,"psi^2 = 1\n" # } # if(b%2){ # null=fund4(dd,0) # phi=globalx # psi=globaly # print "phi^2 - ",dd,"psi^2 = 4\n" # } # print "phi=",phi,", psi=",psi,"\n" # print "alpha=",alpha,", beta=",beta,"\n" # if(b%2==0){ # k=-(a*u^2+b*u*v+c*v^2-d*twodelta0*u-e*twodelta0*v+f*twodelta0^2) # kk=-delta0*(a*e^2-b*e*d+c*d^2+4*f*delta0) # }else{ # k=-dd*(a*e^2-b*e*d+c*d^2+f*dd) # } # print "k=",k,"\n" # print "kk=",kk,"\n" # if(k==0){ # print "k=0\n" # r=-u # s=-v # if(r%dd==0 && s%dd==0){ # print "(x,y)=(",r/dd,",",s/dd,")\n" # gcount=1 # return(1) # }else{ # print "no integer solution\n" # return(0) # } # } # if(b%2){ # print dd,"x=X+(",alpha,")," # print dd,"y=Y+(",beta,")\n" # print "\n" # }else{ # print twodelta0,"x=X+(",alpha,")," # print twodelta0,"y=Y+(",beta,")\n" # print "\n" # } # # # print "First find the solution classes of ",a,"X^2+",b,"XY+",c,"Y^2=",k,"\n" # g=binarygenlist(a,b,c,k,0) # for(h=0;h=0){ print "A[",h,"]/B[",h,"]=",u,"/",v,"\n" }*/ oldl=l;oldm=m;oldn=n 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 if(p==f){/* P_h=P_{h+1}, even period 2h */ /* print "P_",h,"=P_",h+1,"\n"*/ globalr=oldn*(oldu+oldl)+(-1)^h globals=oldn*(oldv+oldm) /*print "globalr=",globalr,",globals=",globals,"\n"*/ return(2*h) } if(q==e){/* Q_h=Q_{h+1}, odd period 2h+1 */ /* print "Q_",h,"=Q_",h+1,"\n"*/ r=k*n+l*m s=n^2+m^2 globalr=r*r+d*s*s globals=2*r*s /*print "globalr=",globalr,",globals=",globals,"\n"*/ return(2*h+1) } } } /* 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 printlegendre(aa,m,parity){ if(parity==0){ if(aa==0){ print 2*m,"k\n" }else{ print "(",2*m,"k+",2*aa,")\n" } }else{ if(aa==0){ print "(",2*m,"k+1)\n" }else{ print "(",2*m,"k+",2*aa+1,")\n" } } } # This solves the diophantine equation ax^2+bxy+cy^2+dx+ey+f=0, where (a,b,c)!=(0,0,0). # It uses a method of John Robertson, communicated on March 25, 2015 to deal with the # general hyperbolic case. # needs squareroot2 powertest3 genfacs gcd posformrep reduceneg congruence powerdd define sswgeneral(a,b,c,d,e,f){ auto gg1,gg2,gg,gcount,h,dd,k,delta,gamma,zeta,epsilon,alpha,beta,psi,phi,x1,y1,null,twoa,absdd,twoadd,bigdelta,g dd=b^2-4*a*c print "D=",dd,"\n" twoa=2*a if(dd>0){ g=squaretest(dd) if(g<0 && a^2==1 && b==0 && d==0 && e==0 && f!=0){ if(a==-1){ a=1;c=-c;f=-f } g=patzgentest(-c,-f) t=fund1(-c) print "general solutions +-(x+ysqrt(",-c,"))(",globalr,"+",globals,"sqrt(",-c,")^n\n" return(g) } if(g<0 && d==0 && e==0 && f!=0){ g=binarygen(a,b,c,-f,0) for(h=0;h0){ g=squaretest(dd) } gg=gcd3(a,b,c) gcount=0 alpha=-b*e+2*c*d beta=-b*d+2*a*e print "alpha=",alpha,", beta=",beta,"\n" print dd,"x=X+(",alpha,")," print dd,"y=Y+(",beta,")\n" print "\n" k=-dd*(a*e^2-b*e*d+c*d^2+f*dd) print "solving (",a,")X^2+(",b,")XY+(",c,")Y^2=",k,"\n" if(k%gg){ print "no solution\n" return(0) } if(dd>0){ g=squaretest(dd) } absdd=abs(dd) if(dd<0 || (dd>0 && g<0)){ if(k==0){ print "k=0\n" if(alpha%absdd==0 && beta%absdd==0){ print "(x,y)=(",alpha/dd,",",beta/dd,")\n" return(1) }else{ print "no integer solution\n" return(0) } } } if(dd<0 && k){# ellipse a=a/gg;b=b/gg;c=c/gg;k=k/gg gcount=ellipse(a,b,c,k,dd,alpha,beta) return(gcount) } if(dd>0 && g<0 && gg>1 &&k!=0){# hyperbola gcount=hyperbolagg(a,b,c,d,e,f,dd,gg) return(gcount) } if(dd>0 && g<0 && gg==1 && k!=0){# hyperbola print "a=",a,",b=",b,",c=",c,",dd=",dd,",g=",g,",k=",k,",alpha=",alpha,",beta=",beta,"\n" gcount=hyperboladd(a,b,c,k,alpha,beta,dd) return(gcount) } gcount=ddsquare(a,b,c,k,g,twoa,dd,alpha,beta) return(gcount) } define hyperboladd(a,b,c,k,alpha,beta,dd){ auto gcount,null,phi,psi,g,bmod2,h,gamma,delta,epsilon,zeta,temp1,temp2,deltaover2,zetaover2,\ z1,z2,temp,x1,y1,tempx,tempy,bigdelta,phiover2 # globalx, globaly,globalbinarygenx,globalbinarygeny null=fund4(dd,0) phi=globalx psi=globaly # print "phi^2 - ",dd,"psi^2 = 4\n" # print "phi=",phi,", psi=",psi,"\n" print "First find the solution classes of ",a,"X^2+",b,"XY+",c,"Y^2=",k,"\n" bmod2=b%2 if(bmod2==0){ bigdelta=dd/4 phiover2=phi/2 } g=binarygenlist(a,b,c,k,0) for(h=0;h1){ print "There are ",gcount," solution classes\n" } if(gcount==1){ print "There is one solution class\n" } if(gcount==0){ print "There are no solution classes\n" } return(gcount) } #dd>0 && g>0 define ddsquare(a,b,c,k,g,twoa,dd,alpha,beta){ auto g1,g2,fourak,h,tcount,gg,t,g1timesg2,flag,xhplusalpha,yhplusalpha #global globallinearnx,globallinearny fourak=4*a*k print "k=",k,"\n" if(a!=0){ g1=gcd(twoa,b+g) g2=gcd(twoa,b-g) g1timesg2=g1*g2 if(fourak%g1timesg2){# g1g2 does not divide 4ak print "no solution\n" return(0) }else{ print "i.e., solving (",twoa/g1,"X+(",(b+g)/g1,")Y)(",twoa/g2,"X+(",(b-g)/g2,")Y)=",fourak/g1timesg2,"\n" } if(k){# g1g2 divides 4ak and k is nonzero print "This gives finitely many solutions on factorising the rhs\n" t=lineareqnn(twoa/g1,(b+g)/g1,twoa/g2,(b-g)/g2,0,0,fourak/g1timesg2) tcount=0 for(h=0;h1){ b=b/gg;c=c/gg;k=k/gg print "i.e., Y(",b,"X+(",c,")Y=",k,"\n" } } if(k){# k is nonzero and now gcd(b,c)=1 t=lineareqnn(0,1,b,c,0,0,k) tcount=0 for(h=0;h1){ print a,"u" } if(a==1){ print "u" } if(a==-1){ print "-u" } if(a<-1){ print a,"u" } if(b>1){ print "+",b,"v" } if(b==1){ print "+v" } if(b==-1){ print "-v" } if(b<-1){ print b,"v" } } define printauuplusbv(a,b){ if(a>1){ print a,"U" } if(a==1){ print "U" } if(a==-1){ print "-U" } if(a<-1){ print a,"U" } if(b>1){ print "+",b,"v" } if(b==1){ print "+v" } if(b==-1){ print "-v" } if(b<-1){ print b,"v" } } define lineareqn(a,b,c,d,e,f){ auto delta1,delta2,absdelta,delta delta1=e*d-b*f delta2=a*f-e*c delta=a*d-b*c absdelta=abs(delta) if(delta1%absdelta==0 && delta2%absdelta==0){ globallinearx=delta1/delta globallineary=delta2/delta return(1) }else{ return(0) } } define lineareqnn(a,b,c,d,r,s,n){ auto g,t,tcount,f,twog tcount=0 absn=abs(n) if(absn!=1){ t=omega(absn) g=divisors(qglobal[],kglobal[],t) }else{ g=1 divisor[0]=1 } for(f=0;f0 g=sqrt(v) if(g^2==v){ if(b==0){ if((g-d)%twoa==0){ print "solution x=",(g-d)/twoa," y, arbitrary\n" } if((-g-d)%twoa==0){ print "solution x=",-(g+d)/twoa,", y arbitrary\n" } return(-2) }else{#b!=0 t=gcd(twoa,b) tcount=0 if((-d+g)%t==0){ print "we get a line of solutions ",twoa/t,"x+(",b/t,")y=",(-d+g)/t,"\n" tcount=1 } if((-d-g)%t==0){ print "and a line of solutions ",twoa/t,"x+(",b/t,")y=",(-d-g)/t,"\n" tcount=tcount+1 } if(tcount>0){ return(-2) }else{ print "no solutions\n" return(0) } } }else{ print "no solutions\n" return(0) } }#end of v>0 }#end of v!=0 return }#end of if u==0 # we can now proceed with the knowledge that u is nonzero below. # solving (t+d)^2=uy+v, u nonzero absu=abs(u) if(flag){ print "absu=",absu,"\n" print "u=",u,"\n" print "v=",v,"\n" print "solve T^2=",v," (mod ",absu,")\n" } g=squareroot0(v,absu) if(g==0){ print "no solutions\n" return(0) } tcount=0 for(h=0;h1){ print tcount," families of solutions\n" }else{ print "one family of solutions\n" } }else{ print "no solutions\n" } }#end of define ddzero define reducetest(u,v,d){ auto f f=sqrt(d) if(v>0 && u>0 && u<=f && f