/* bc file squareroot2 */
/* min(x,y) */
define min(x,y){
if(y<x)return(y)
return(x)
}
define rsv(x,y){
if(x>y)
return(1)
if(x==y)
return(0)
if(x<y)
return(-1)
}
define sign(a){
if(a>0) return(1)
if(a<0) return(-1)
return(0)
}
define abs(n){
if(n>=0) return(n)
return(-n)
}
define mod(a,b){
auto c
c=a%b
if(a>=0) return(c)
if(c==0) return(0)
return(c+b)
}
define mpower(a,b,c){
auto x,y,z
x=a
y=b
z=1
while(y>0){
while(y%2==0){
y=y/2
x=(x*x)%c
}
y=y-1
z=mod(z*x,c)
}
return(z)
}
/* Shanks-Tonelli algorithm.
* See "Square roots from 1; 24, 51, 10 to Dan Shanks",
* Ezra Brown, The College Mathematics Journal, 30, No. 2, 82-95, 1999
*/
define tonelli(a,p){
auto b,e,g,h,i,m,n,q,r,s,t,x,y,z
q=p-1
t=mpower(a,q/2,p)
if(t==q)return(0)/* a is a quadratic non-residue mod p */
s=q
e=0
while(s%2==0){s=s/2;e=e+1}
/* p-1=s*2^e */
x=mpower(a,(s+1)/2,p)
b=mpower(a,s,p)
if(b==1)return(x)
for(n=2;n>=2;n++){
t=mpower(n,q/2,p)
if(t==q) break
}
/* n is the least quadratic non-residue mod p */
g=mpower(n,s,p)
r=e
while(1){
y=b
for(m=0;m<=r-1;m++){
if(y==1)break
y=(y*y)%p
}
/*print "(m,x,b,g,r)=(",m,",",x,",",b,",",g,",",r,")\n"*/
/* ord_p(b)=2^m, m <= r-1, so m has decreased */
if(m==0)return(x)
z=r-m-1
h=g
for(i=0;i<z;i++) h=(h*h)%p
x=(x*h)%p
b=(b*h*h)%p
g=(h*h)%p
r=m
}
}
/* the congruence mx=p(mod n) */
define cong1(m,p,n){
auto a,b,x,y
a=gcd(m,n)
if(mod(p,a)>0){ "not soluble "; return(-1) }
b=gcd1(m,n)
y=n/a
p=p/a
x=mod(b*p,y)
return(x)
}
define gcd1(m,n){
auto a,b,c,h,k,l,q,t
if(n==0) return(sign(m))
a=m
b=abs(n)
c=mod(a,b)
if(c==0) return(0)
l=1
k=0
while(c>0){
q=(a-c)/b
a=b
b=c
c=mod(a,b)
h=l-q*k
l=k
k=h
}
return(k)
}
define gcd2(m,n){
auto a,b,c,h,k,l,q,t
if(n==0) return(0)
a=m
b=abs(n)
c=mod(a,b)
if(c==0) return(sign(n))
l=0
k=1
while(c>0){
q=(a-c)/b
a=b
b=c
c=mod(a,b)
h=l-q*k
l=k
k=h
}
q=(a-c)/b
if(n<0) return(-k)
return(k)
}
define gcd(m,n){
auto a,b,c
a=abs(m) /* a=r[0] */
if(n==0) return(a)
b=abs(n) /* b=r[1] */
c=a%b /* c=r[2]=r[0] mod(r[1]) */
while(c>0){
a=b
b=c
c=a%b /* c=r[j]=r[j-2] mod(r[j-1]) */
}
return(b)
}
/* inv(a,n) returns the inverse of a (mod n), n>0 if gcd(a,n)=1.
* uses the fact that |s[j]|<n in general.
*/
define inv(a,n){
auto s
if(gcd(a,n)>1){print "gcd(a,n)>1\n";return(0)}
s=gcd1(a,n)
if(s<0)return(s+n)
return(s)
}
/* the bth power of a, where a is an integer, b a positive integer.*/
define exp(a,b){
auto x,y,z
x=a
y=b
z=1
while(y>0){
while(y%2==0){
y=y/2
x=x*x
}
y=y-1
z=z*x
}
return(z)
}
pglobal[0]=2
pglobal[1]=3
pglobal[2]=5
pglobal[3]=7
pglobal[4]=11
pglobal[5]=13
pglobal[6]=17
pglobal[7]=19
pglobal[8]=23
pglobal[9]=29
pglobal[10]=31
pglobal[11]=37
pglobal[12]=41
pglobal[13]=43
pglobal[14]=47
pglobal[15]=53
pglobal[16]=59
pglobal[17]=61
pglobal[18]=67
pglobal[19]=71
pglobal[20]=73
pglobal[21]=79
pglobal[22]=83
pglobal[23]=89
pglobal[24]=97
pglobal[25]=101
pglobal[26]=103
pglobal[27]=107
pglobal[28]=109
pglobal[29]=113
pglobal[30]=127
pglobal[31]=131
pglobal[32]=137
pglobal[33]=139
pglobal[34]=149
pglobal[35]=151
pglobal[36]=157
pglobal[37]=163
pglobal[38]=167
pglobal[39]=173
pglobal[40]=179
pglobal[41]=181
pglobal[42]=191
pglobal[43]=193
pglobal[44]=197
pglobal[45]=199
pglobal[46]=211
pglobal[47]=223
pglobal[48]=227
pglobal[49]=229
pglobal[50]=233
pglobal[51]=239
pglobal[52]=241
pglobal[53]=251
pglobal[54]=257
pglobal[55]=263
pglobal[56]=269
pglobal[57]=271
pglobal[58]=277
pglobal[59]=281
pglobal[60]=283
pglobal[61]=293
pglobal[62]=307
pglobal[63]=311
pglobal[64]=313
pglobal[65]=317
pglobal[66]=331
pglobal[67]=337
pglobal[68]=347
pglobal[69]=349
pglobal[70]=353
pglobal[71]=359
pglobal[72]=367
pglobal[73]=373
pglobal[74]=379
pglobal[75]=383
pglobal[76]=389
pglobal[77]=397
pglobal[78]=401
pglobal[79]=409
pglobal[80]=419
pglobal[81]=421
pglobal[82]=431
pglobal[83]=433
pglobal[84]=439
pglobal[85]=443
pglobal[86]=449
pglobal[87]=457
pglobal[88]=461
pglobal[89]=463
pglobal[90]=467
pglobal[91]=479
pglobal[92]=487
pglobal[93]=491
pglobal[94]=499
pglobal[95]=503
pglobal[96]=509
pglobal[97]=521
pglobal[98]=523
pglobal[99]=541
pglobal[100]=547
pglobal[101]=557
pglobal[102]=563
pglobal[103]=569
pglobal[104]=571
pglobal[105]=577
pglobal[106]=587
pglobal[107]=593
pglobal[108]=599
pglobal[109]=601
pglobal[110]=607
pglobal[111]=613
pglobal[112]=617
pglobal[113]=619
pglobal[114]=631
pglobal[115]=641
pglobal[116]=643
pglobal[117]=647
pglobal[118]=653
pglobal[119]=659
pglobal[120]=661
pglobal[121]=673
pglobal[122]=677
pglobal[123]=683
pglobal[124]=691
pglobal[125]=701
pglobal[126]=709
pglobal[127]=719
pglobal[128]=727
pglobal[129]=733
pglobal[130]=739
pglobal[131]=743
pglobal[132]=751
pglobal[133]=757
pglobal[134]=761
pglobal[135]=769
pglobal[136]=773
pglobal[137]=787
pglobal[138]=797
pglobal[139]=809
pglobal[140]=811
pglobal[141]=821
pglobal[142]=823
pglobal[143]=827
pglobal[144]=829
pglobal[145]=839
pglobal[146]=853
pglobal[147]=857
pglobal[148]=859
pglobal[149]=863
pglobal[150]=877
pglobal[151]=881
pglobal[152]=883
pglobal[153]=887
pglobal[154]=907
pglobal[155]=911
pglobal[156]=919
pglobal[157]=929
pglobal[158]=937
pglobal[159]=941
pglobal[160]=947
pglobal[161]=953
pglobal[162]=967
pglobal[163]=971
pglobal[164]=977
pglobal[165]=983
pglobal[166]=991
pglobal[167]=997
/*
* the Brent-Pollard method for returning a proper factor of a composite n.
* See R. Brent, Bit 20, 176-184
*/
define brent(n){
auto a,y,r,g,x,i,k,f
a=1;y=0;r=1;g=1;prod=1
while(g==1){
x=y
for(i=1;i<=r;i++)y=(y*y+a)%n
k=0
f=0
while(f==0){
y=(y*y+a)%n
k=k+1
z=abs(x-y)
prod=prod*z
if(k%10==0){
g=gcd(prod,n)
prod=1
if (g>1)f=1
}
if(k>=r)f=1
}
r=2*r
if(g==n || r == 16384){
g=1
r=1
y=0
a=a+1
if(a==5){
return(0)
}
}
}
return(g)
}
define pollard(n){
auto i,p,t,b
b=t=2
p=1
for(i=2;i<=10^4;i++){
/*if(i%10==0){"i=";i}*/
t=mpower(t,i,n)/* now t=b ^(i!) */
p=gcd(n,t-1)
if(p>1){
if(p<n){
return(p)
}
}
if(p==n){
i=2;t=b=b+1
"b= ";b
if(b==5) {
return(0)}
}
}
return(0)
}
/*
* babydivide(n) returns n/m, where m is the part of n composed of primes
* qglobal[i],...,qglobal[jglobal-1] < 1000, where i is the value of global
* variable jglobal before babydivide(n) is called.
* qglobal[] and exponent array kglobal[] are global variables.
*/
define babydivide(n){
auto a,u,x,k,i
a=167
x=n
for(i=0;i<=a;i++){
k=0
if(x<pglobal[i])break
if(x%pglobal[i]==0){
while(x%pglobal[i]==0){
k=k+1
x=x/pglobal[i]
}
qglobal[jglobal]=pglobal[i]
kglobal[jglobal]=k
jglobal=jglobal+1
}
}
return(x)
}
/*
* n > 1 and odd, gcd(n,b)=1, or more generally, n does not divide b.
* If miller(n,b)=1, then n passes Miller's test for base b and n is
* either a prime or a strong pseudoprime to base b.
* if miller(n,b)=0, then n is composite.
*/
define miller(n,b){
auto a,s
s=(n-1)/2
a=s
while(a%2==0)a=a/2
b=mpower(b,a,n)
if(b==1)return(1)
while(a<=s){
if(b==n-1)return(1)
b=mod(b*b,n)
a=2*a
}
return(0)
}
/*
* n > 1 is distinct from pglobal[0],...,pglobal[4].
* if q(n)=1, then n passes Miller's test for bases pglobal[0],...,pglobal[4]
* and is likely to be prime.
* if q(n)=0, then n is composite.
*/
define q(n){
auto i
for(i=0;i<=4;i++){
if(miller(n,pglobal[i])==0){
return(0)
}
}
return(1)
}
/*
* n>1 is not divisible by pglobal[0],...,pglobal[167].
* v(n) returns a factor of n which is < 1,000,000 (and hence prime)
* or which passes Miller's test for bases pglobal[0],...,pglobal[4] and is
* therefore likely to be prime.
*/
define v(n){
auto f,x,y,b
b=1000
if(n<b*b){
return(n)
}
if(q(n)==1)return(n)
f=1
x=n
while(f!=0){
y=brent(x)
if(y==0){
y=pollard(x)
if(y==0){
print "Exiting factoring program with x=",x,"\n"
halt
}
}
x=min(x/y,y)
if(x<b*b)return(x)
if(q(x)==1)f=0
}
return(x)
}
/* A quasi-prime (q-prime) factor of n is > 1,000,000, passes Miller's test
* and is not divisible by pglobal[0],...,pglobal[167]. It is likely to be
* prime.
* primefactors(n) returns the number lglobal-t of q-prime factors of n.
* rglobal[] and lglobal are global variables.
* The prime and q-prime factors of n are qglobal[i],...,qglobal[jglobal-1]
* with exponents
* kglobal[i],...,kglobal[jglobal-1] where i is the value of the global
* variable jglobal before primefactors(n) is called.
*/
define primefactors(n){
auto b,p,x,k,t
b=1000
x=babydivide(n)
t=lglobal
while(x!=1){
k=0
p=v(x)
while(x%p==0){
k=k+1
x=x/p
}
if(p>b*b){
rglobal[lglobal]=p
lglobal=lglobal+1
}
qglobal[jglobal]=p
kglobal[jglobal]=k
jglobal=jglobal+1
}
return(lglobal-t)
}
/*
* Selfridge's test for primality - see "Prime Numbers and Computer
* Methods for Factorization" by H. Riesel, Theorem 4.4, p.106.
* input: n (q-prime)
* first finds the prime and q-prime factors of n-1. If no q-prime factor
* is present and 1 is returned, then n is prime. However if at least one
* q-prime is present and 1 is returned, then n retains "q-prime" status.
* If 0 is returned, then either n or one of the q-prime factors of n-1 is
* composite.
*/
define selfridge(n){
auto i,x,s,t,u
i=jglobal
u=primefactors(n-1)
cglobal=u+cglobal
/* cglobal,jglobal,lglobal are global variables. */
/* primefactors(n-1) returns jglobal-i primes and q-primes
qglobal[i],...,qglobal[jglobal-1] */
/* and q-primes rglobal[l-u],...,rglobal[lglobal-1], where u>=0. */
while(i<=jglobal-1){
for(x=2;x<n;x++){
s=mpower(x,n-1,n)
if(s!=1){
return(0)
}
t=mpower(x,(n-1)/qglobal[i],n)
if(t!=1){
i=i+1
break
}
}
if(x==n){ /* unlikely to be tested */
return(0)
}
}
return(1)
}
/*
* a factorization of n into prime and q-prime factors is first obtained.
* Selfridge's primality test is then applied to any q-prime factors;
* the test is applied repeatedly until either a q-prime is found to be
* composite (in which case the initial factorization is doubtful and
* an extra base should be used in Miller's test) or we run out of q-primes,
* in which case every q-prime factor of n is certified as prime.
* omega(n) returns the number of distinct prime factors of n.
*/
define omega(n){
auto i,s,t
jglobal=lglobal=0
cglobal=primefactors(n)
/* cglobal gives a progress count of the number of q-prime factors */
/* yet to be tested by Selfridges' method */
s=jglobal
/* s is the number of prime and q-prime factors of n */
if(cglobal==0){
return(s)
}
while(cglobal>0){
t=selfridge(rglobal[i])
if(t==0){
return(0)
}
i=i+1
cglobal=cglobal-1
}
return(s)
}
/* the Chinese remainder theorem for the congruences x=a(mod m)
* and x=b(mod n), m>0, n>0, a and b arbitrary integers.
* The construction of O. Ore, American Mathematical Monthly,
* vol.59,pp.365-370,1952, is implemented.
* Same as chinese() in gcd, but with printing suppressed, except that -1
* is returned if there is no solution.
*/
define chinesee(a,b,m,n) {
auto c,d,x,y,z
d = gcd(m,n)
if(mod(a-b,d)!=0){
return(-1)
}
x= m/d;y=n/d
z=(m*n)/d
c = mod(b*gcd1(x,y)*x + a*gcd2(x,y)*y,z)
return(c)
}
define chineseae(a[],m[],n){
auto i,m,x
m=m[0]
x=a[0]
for(i=1;i<n;i++){
x=chinesee(a[i],x,m[i],m)
if(x==-1)return(-1)
m=lcm(m[i],m)
}
return(x)
}
/* lcm(a,b) for any integers a and b */
define lcm(a,b){
auto g
g=gcd(a,b)
if(g==0) return(0)
return(abs(a*b)/g)
}
/* Solving x^2=a (mod p^n), p an odd prime not dividing a. */
define sqroot1(a,p,n){
auto i,k,t,u,v
t=mpower(a,(p-1)/2,p)
if (t!=1){
return (-1)
}
if (n==1){
return(tonelli(a,p))
}else{
u=sqroot1(a,p,n-1)
v=u^2-a
for(i=0;i< n-1;i++)v=v/p
k=cong1(2*u,-v,p)
return(u+k*exp(p,n-1))
}
}
/* Solving x^2=a (mod 2^n), a odd. */
define sqroot2(a,n){
auto i, u, v
if (n==1){/* x=1(mod 2) */
return(1)
}
if (n==2){
if(mod(a,4)!=1){
return(-1)
}else{
/* x=1(mod 2) */
return (1)
}
}
if(mod(a,8)!=1){
return (-1)
}
/* two solutions mod 2^{n-1} */
if (n==3){
return(1)
}else{
u=sqroot2(a,n-1)
v=u^2-a
for(i=0;i< n-1;i++)v=v/2
if ((v+1)%2 == 0){
u=u+2^(n-2)
}
return(u)
}
}
/* Solving x^2=a (mod p^n), p a prime possibly dividing a.
* If p does not divide a, the story is well-known:
* If p=2 and n=1, solution is x=1(mod 2) :number=1,modulus=2
* If p=2 and n=2, solution is x=1(mod 2) :number=1,modulus=2
* If p=2 and n>2, two solutions mod 2^(n-1):number=2,modulus=2^(n-1)
* If p>2, two solutions mod p^n: number=2,modulus=p^n
* If p divides a, there are two cases:
* (i) p^n divides a,
* (ii) p^n does not divide a.
* In case (i) there are two cases:
* (a) n=2m+1, (b) n=2m.
* In case (a), the solution is x=0 (mod p^(m+1)): number=1, modulus=p^(m+1)
* In case (b), the solution is x=0 (mod p^m): number=1, modulus=p^m
* In case (ii) gcd(a,p^n)=p^r, r<n. If r is odd, no solution.
* If r=2m, write x=(p^m)*X, d=(p^2m)*A, p not dividing A.
* Then x^2=a (mod p^n) becomes X^2=A (mod p^(n-2m)).
* If p is odd, this has 2 solutions X mod p^(n-2m)
* and we get two solutions x (mod p^(n-m)): number=2, modulus=p^(n-m)
* If p=2, we get
* 1 solution mod (2^(n-m)) if n-2m=1, viz. x=2^m (mod 2^(m+1);
* number=1, modulus=2^(m+1)
* 1 solution mod (2^(n-m-1) if n-2m=2 and A=1 (mod 4), viz. x=2^m (mod 2^(m+1);
* number=1, modulus=2^(m+1)
* 2 solutions mod (2^(n-m-1)) if n-2m>2 and A=1 (mod 8).
* number=2, modulus=2^(n-m-1)
*/
define sqroot3(a,p,n){
auto aa,d,m,r,s,t,u,v
d=a%p
if(d){
if(p==2){
t=sqroot2(a,n)
if(n==1 || n==2){
number=1
exponent=2
}else{
number=2
exponent=2^(n-1)
}
return(t)
}else
t=sqroot1(a,p,n)
exponent=p^n
number=2
return(t)
}else{
u=p^n
v=a%u
if(v==0){
number=1
if(n%2==0){
exponent=p^(n/2)
}else{
exponent=p^(1+n/2)
}
return(0)
}
r=0
aa=a
v=aa%p
while(v==0 && r<=n){
aa=aa/p
v=aa%p
r=r+1
}
if(r%2){
return(-1)
}else{
m=r/2
s=n-r
if(p==2){
t=sqroot2(aa,s)
if(s==1 || s==2){
exponent=2^(m+1)
number=1
}else{
number=2
exponent=2^(n-m-1)
}
if(t==-1){
return(-1)
}else{
return((2^m)*t)
}
}else{
t=sqroot1(aa,p,s)
exponent=p^(n-m)
number=2
if(t==-1){
return(-1)
}else
return((p^m)*t)
}
}
}
}
define solution(){
auto i,j,x
for(i=0; i<numprimefactors; i++){
a[i]=selection[i]*soln[i]
}
x=chineseae(a[],moduli[],numprimefactors)
if(2*x<=reduced_modulus){
reduced_solution[globalsqcount]=x
}else{
reduced_solution[globalsqcount]=x-reduced_modulus
}
globalsqcount=globalsqcount+1
}
/* the recursive trick in recur was supplied by Peter Adams, 13th April 2004. */
define recur(depth){
auto j,ns,s,t
if(depth==numprimefactors){
t=solution();
return;
}
ns=nsolns[depth];
for(j=0; j<ns; j++)
{
if(j%2){
selection[depth]=1
}else{
selection[depth]=-1
}
s=recur(depth+1)
}
}
/*
* r=sqroot(d,n,e) returns the solutions of x^2=d (mod n).
* r is the number of solutions (mod n).
* If e=1, we print the solutions (mod reduced_modulus) as
* reduced_solution[0],...,reduced_solution[globalsqcount-1].
* if e=0, solutions are not printed. Used eg. in cornacchia().
* If omega(n)>1, we use the Chinese remainder theorem after solving mod
* qglobal[i]^kglobal[i], i=0,...,e-1=omega(n)-1.
* The array solution[0],...,solution[numbr-1] consists of the solutions
* (mod n) in the range 0<=x<=n/2 and is used in cornacchia().
* numbr is a global variable.
*/
define sqroot(d,n,flag){
auto e,i,j,k,t,z
globalsqcount=0
reduced_modulus=1
if(n==1){
reduced_solution[0]=0
numbr=1
solution[0]=0
if(flag){
print "Solution of x^2=",d,"(mod ",n,") is x=0 (mod 1)\n"
print "The number of solutions (mod ",n,") is "
}
return(1)
}
numprimefactors=omega(n)
for (i=0;i<numprimefactors;i++){
globalsqc=sqroot3(d,qglobal[i],kglobal[i])
if(globalsqc==-1){
if(flag){
print "x^2=",d," mod(",n,") has no solution\n"
}
return(0)
}
if(number==1){
nsolns[i]=1
}else{
nsolns[i]=2
}
soln[i]=globalsqc
moduli[i]=exponent
reduced_modulus=reduced_modulus*moduli[i]
}
if(numprimefactors==1){
if(2*soln[0]>reduced_modulus){
reduced_solution[0]=reduced_modulus-soln[0]
reduced_solution[1]=-reduced_solution[0]
globalsqcount=2
}
if(2*soln[0]<reduced_modulus && soln[0]){
reduced_solution[0]=soln[0]
reduced_solution[1]=-soln[0]
globalsqcount=2
}
if(2*soln[0]==reduced_modulus){
reduced_solution[0]=soln[0]
globalsqcount=1
}
if(soln[0]==0){
reduced_solution[0]=soln[0]
globalsqcount=1
}
}else{
t=recur(0)
}
if(flag){
print "Solution of x^2=",d,"(mod ",n,") is \n"
for(j=0;j<globalsqcount;j++){
print reduced_solution[j]," (mod ",reduced_modulus,")\n"
}
}
/* now we omit the negative solutions and extend to solutions mod
n */
z=n/reduced_modulus
numbr=0
for(j=0;j<globalsqcount;j++){
for(k=0;k<z;k++){
temp=mod(reduced_solution[j]+k*reduced_modulus,n)
if(2*temp<=n){
solution[numbr]=reduced_solution[j]+k*reduced_modulus
numbr=numbr+1
}
}
}
if(flag){
print "The number of solutions (mod ",n,") is "
}
return(globalsqcount*z)
}
/* Returns the positive primitive solutions and their number (x,y) of
* a*x^2+b*y^2=m, x>=y.
* Here a>0,b>0,m>a+b, gcd(a,b)=1=gcd(a,m).
* If a=b=1, we get only solutions with 0<y<x.
* See A. Nitaj, 'L'algorithme de Cornacchia', Expositiones Mathematicae 13
* (1995), 358-365.
* Finished coding on 15th April 2004 after getting a better version of sqroot() in place.
*/
define cornacchia(a,b,m){
auto i,a1,a2,ma,mb,tt,t1,aa,bb,tee2,jj,r,t2,tee1,q1,qq,soln_number,temp,count
if(a<=0){
print "a<=0\n"
halt
}
if(b<=0){
print "b<=0\n"
halt
}
if(a+b>m){
print "a+b>m\n"
halt
}
temp=gcd(a,m)
if(temp>1){
print "gcd(a,m)>1\n"
halt
}
temp=gcd(a,b)
if(temp>1){
print "gcd(a,b)>1\n"
halt
}
a1=inv(a,m)
a2=a1*(-b)
ma=m/a
mb=m/b
n=sqroot(a2,m,0)
print "Solving "
if(a>1){
print a
}
print"x^2+"
if(b>1){
print b
}
print "y^2=",m,":\n"
if(n==0){
print "Number of primitive solutions="
return(0)
}
soln_number=0
count=0
for(i=0;i<numbr;i++){
tt=0
t1=1
aa=m
bb=solution[i]
tee2=0
jj=0
r=aa
t2=0
while(tee2<=0){
tee1=rsv(r*r,ma)
if(tee1<=0){
print "(x,y)=(",r,","
x_cornacchia[count]=r
if(t2<0){
t2= -t2
}
y_cornacchia[count]=t2
count=count+1
print t2,")\n"
soln_number=soln_number+1
break
}
jj=jj+1
if(jj==1){
r=bb
t2=1
tee2=rsv(t2,mb)
}
if(jj>1){
r=mod(aa,bb)
q1=aa/bb
aa=bb
bb=r
qq= -q1
t2=tt+(qq*t1)
tee2=rsv(t2*t2,mb)
tt=t1
t1=t2
}
}
}
if(a==1 && b==1){
print "Number of positive primitive solutions with y<=x is "
}else{
print "Number of positive primitive solutions is "
}
print soln_number,"\n"
return(soln_number)
}
/* Solving ax^2+bx+c=0 (mod n), where gcd(a,n)=1.
* k=quadratic(a,b,c,n,flag) returns number k of solutions.
* If k>0, we get a global array:
* quadratic_solution[0],...,quadratic_solution[k-1].
* First solve t^2=b^2-4ac (mod 4n), -n<t<=n.
* Note that t=b (mod 2) here.
* Then solve ax=(t-b)/2 (mod n), 0<=x<n.
* flag=1 prints the solution
*/
define quadratic(a,b,c,n,flag){
auto d,i,j,k,s,temp[]
d=b^2-4*a*c
s=sqroot(d,4*n,0)
if(s==0){
if(flag){
print "no solution\n"
}
return(0)
}
/* now we have x=solution[0],...,solution[numbr-1], 0<=x<=2n */
k=0
for(i=0;i<numbr;i++){
if(solution[i]<=n){
temp[k]=solution[i]
k=k+1
if(solution[i]>0 && solution[i]<n){
temp[k]=-solution[i]
k=k+1
}
}
}
for(j=0;j<k;j++){
quadratic_solution[j]=cong1(a,(temp[j]-b)/2,n)
if(flag){
print "quadratic_solution[",j,"]=",quadratic_solution[j],"\n"
}
}
return(k)
}
/* This checks John's conjecture, given a solution (x,y), x>=1,k-1>y>=1, of x^2+1=(k^2+1)(y^2+1)
he says there exist r,s,a,b, 1<=r<s,1<=a<b,gcd(r,s)=gcd(a,b)=1,r^2+s^2=k^2+1,a^2+b^2=y^2+1,
x=ar+bs, rb-sa=1 or -1. If it returns count=1, the conjecture holds.
If count=0, this means the construction only produces (X,Y) with X^2+Y^2=(k^2+1)(Y^2+1) but (X,Y) != (x,1)
or (x,-1).
*/
define testjpr(x,y,k,flag){
auto count1, count2,i,j,r[],s[],a[],b[],xtemp,pmtemp
count=0
count1=cornacchia0(1,1,k^2+1)
for(i=0;i<count1;i++){
s[i]=x_cornacchia[i]
r[i]=y_cornacchia[i]
}
count2=cornacchia0(1,1,y^2+1)
for(j=0;j<count2;j++){
b[j]=x_cornacchia[j]
a[j]=y_cornacchia[j]
}
for(i=0;i<count1;i++){
for(j=0;j<count2;j++){
xtemp=a[j]*r[i]+b[j]*s[i]
pmtemp=b[j]*r[i]-a[j]*s[i]
if(flag){
print "(X,Y,gcd(X,Y))=(",xtemp,",",pmtemp,",",gcd(xtemp,pmtemp),"):"
print "(a,b,r,s)=(",a[j],",",b[j],",",r[i],",",s[i],")"
}
if(x==xtemp && pmtemp*pmtemp==1){
print "(a,b,r,s)=(",a[j],",",b[j],",",r[i],",",s[i],")"
globala=a[j]
/*if(globala>1){
temp=(k^2+1)/2
print "\n\t2a(x-yk)-r=",2*a[j]*(x-y*k)-r[i],",a/2=",a[j]/2
if(temp>x){
print ",x<(k^2+1)/2\n"
}else{
print ",x>(k^2+1)/2!!!!\n"
}
}*/
globalb=b[j]
return
}
}
}
}
define testjprprint(x,y,k,flag){
auto count1, count2,i,j,r[],s[],a[],b[],xtemp,pmtemp
count=0
count1=cornacchia0(1,1,k^2+1)
for(i=0;i<count1;i++){
s[i]=x_cornacchia[i]
r[i]=y_cornacchia[i]
}
count2=cornacchia0(1,1,y^2+1)
for(j=0;j<count2;j++){
b[j]=x_cornacchia[j]
a[j]=y_cornacchia[j]
}
for(i=0;i<count1;i++){
for(j=0;j<count2;j++){
xtemp=a[j]*r[i]+b[j]*s[i]
pmtemp=b[j]*r[i]-a[j]*s[i]
if(flag==1){
print "(X,Y,gcd(X,Y))=(",xtemp,",",pmtemp,",",gcd(xtemp,pmtemp),"):"
print "(a,b,r,s)=(",a[j],",",b[j],",",r[i],",",s[i],")"
}
if(flag==0){
if(x==xtemp && pmtemp*pmtemp==1){
print "(a,b,r,s)=(",a[j],",",b[j],",",r[i],",",s[i],")"
globala=a[j]
globalb=b[j]
return
}
}
if(flag==2){
if(x==xtemp && pmtemp*pmtemp==1){
print r[i]," & ",s[i]
globala=a[j]
globalb=b[j]
return
}
}
}
}
}
define cornacchia0(a,b,m){
auto i,a1,a2,ma,mb,tt,t1,aa,bb,tee2,jj,r,t2,tee1,q1,qq,soln_number,temp,count,n
if(a<=0){
print "a<=0\n"
halt
}
if(b<=0){
print "b<=0\n"
halt
}
if(a+b>m){
print "a+b>m\n"
halt
}
temp=gcd(a,m)
if(temp>1){
print "gcd(a,m)>1\n"
halt
}
temp=gcd(a,b)
if(temp>1){
print "gcd(a,b)>1\n"
halt
}
a1=inv(a,m)
a2=a1*(-b)
ma=m/a
mb=m/b
n=sqroot(a2,m,0)
/*print "Solving "
if(a>1){
print a
}
print"x^2+"
if(b>1){
print b
}
print "y^2=",m,":\n"
if(n==0){
print "Number of primitive solutions="
return(0)
}*/
soln_number=0
count=0
for(i=0;i<numbr;i++){
tt=0
t1=1
aa=m
bb=solution[i]
tee2=0
jj=0
r=aa
t2=0
while(tee2<=0){
tee1=rsv(r*r,ma)
if(tee1<=0){
/*print "(x,y)=(",r,","*/
x_cornacchia[count]=r
if(t2<0){
t2= -t2
}
y_cornacchia[count]=t2
count=count+1
/*print t2,")\n"*/
soln_number=soln_number+1
break
}
jj=jj+1
if(jj==1){
r=bb
t2=1
tee2=rsv(t2,mb)
}
if(jj>1){
r=mod(aa,bb)
q1=aa/bb
aa=bb
bb=r
qq= -q1
t2=tt+(qq*t1)
tee2=rsv(t2*t2,mb)
tt=t1
t1=t2
}
}
}
/*if(a==1 && b==1){
print "Number of positive primitive solutions with y<=x is "
}else{
print "Number of positive primitive solutions is "
}*/
/*print soln_number,"\n"*/
return(soln_number)
}
/* This returns the number of Type 1 exceptional solutions for the equation x^2+1=(k^2+1)(y^2+1). */
/* Dujella's conjecture implies that at most one can be met. */
define rs(k){
auto count,count1,s[],r[],i,x,y
count1=0
count=cornacchia0(1,1,k^2+1)
for(i=0;i<count;i++){
s[i]=x_cornacchia[i]
r[i]=y_cornacchia[i]
if(r[i]>1 && s[i]>1){
if((s[i]-1)%r[i]==0){
y=(s[i]-1)/r[i]
x=r[i]+y*s[i]
print "k=",k,",r[",i,"]=",r[i],",s[",i,"]=",s[i],",x=",x,",y=",y,",gcd(x,y)=",gcd(x,y),",e=1\n"
count1=count1+1
if(y==2){
continue
}
}
if((s[i]+1)%r[i]==0){
y=(s[i]+1)/r[i]
x=r[i]+y*s[i]
print "k=",k,",r[",i,"]=",r[i],",s[",i,"]=",s[i],",x=",x,",y=",y,",gcd(x,y)=",gcd(x,y),",e=-1\n"
count1=count1+1
}
}
}
return(count1)
}
/* this returns the number of Type 1 exceptional solutions in the range [m,n]. */
define rstest(m,n){
auto k,count,count1
count=0
for(k=m;k<=n;k++){
count1=rs(k)
if(count1>1){
print "k=",k,",count1=",count1,">1\n"
return
}else{
if(count1==1){
count=count+1
}else{
/* print "k=",k,",count1=0\n"*/
}
}
}
return(count)
}
/* This returns the number of Type 1 exceptional solutions for the equation x^2+1=(k^2+1)(y^2+1). */
/* Dujella's conjecture implies that at most one can be met. */
define rs0(k){
auto count,count1,s[],r[],i,x,y,u
count1=0
count=cornacchia0(1,1,k^2+1)
for(i=0;i<count;i++){
s[i]=x_cornacchia[i]
r[i]=y_cornacchia[i]
if(r[i]>1 && s[i]>1){
if((s[i]-1)%r[i]==0){
y=(s[i]-1)/r[i]
x=r[i]+y*s[i]
u=sqrt(r[i])
if(x%y){
print "(",u,",",k/u,",",s[i],",e=1)\n"
}
count1=count1+1
if(y==2){
continue
}
}
if((s[i]+1)%r[i]==0){
y=(s[i]+1)/r[i]
x=r[i]+y*s[i]
u=sqrt(r[i])
if(x%y){
print "(",u,",",k/u,",",s[i],",e=-1)\n"
}
count1=count1+1
}
}
}
return(count1)
}
/* this returns the number of Type 1 exceptional solutions in the range [m,n]. */
define rstest0(m,n){
auto k,count,count1
count=0
for(k=m;k<=n;k++){
count1=rs0(k)
if(count1>1){
print "k=",k,",count1=",count1,">1\n"
return
}else{
if(count1==1){
count=count+1
}
}
}
return(count)
}
/*
This takes an integer which is expressible as the sum of two
relatively prime squares n=r^2+s^2, r<s and for each representation it calculates x < n/2 such that
x^2=-1 (mod n) and expresses x=ar+bs, 0<=a<=b,1<r<s, br-as=\eps= 1 or -1 and a<=r/2, b<=s/2.
Also xr=s\eps(mod n). We get a=0 only for the representation n=1^2+s^2.
See http://www.numbertheory.org/pdfs/hermite_serret.pdf */
define squares(n,flag){
auto i,l
r=sqroot(-1,n,1)
print r,"\n"
if(r==0){
print "no solution of x^2=-1 (mod n)\n"
return
}
for(i=0;i<numbr;i++){
l=euclid0(n,solution[i],flag)
print "length=",l,"\n"
}
}
define euclid0(m,n,flag){
auto a,b,c,r,q,i,j,l,s0,s1,s2,t0,t1,t2,s[],t[],rr[],qq[]
s0=1;s1=0
s[0]=1;s[1]=0
t0=0;t1=1
t[0]=0;t[1]=1
a=m
b=n
rr[0]=m
rr[1]=n
r=m%n
q=(a-r)/b
qq[1]=q
i=1
while(r>0){
i=i+1
s2=(-q)*s1+s0
s[i]=s2
t2=(-q)*t1+t0
t[i]=t2
a=b
b=r
r=a%b
q=(a-r)/b
qq[i]=q
rr[i]=b
s0=s1;s1=s2
t0=t1;t1=t2
}
s2=(-q)*s1+s0
t2=(-q)*t1+t0
i=i+1
s[i]=s2
t[i]=t2
rr[i]=0
print "the Euclidean algorithm for m/n with m = ",m, " and n = ",n,"\n"
l=i-1
if(flag){
for(j=0;j<=i;j++){
if(j>0 && j<i){
print "q[",j,"]=",qq[j],","
}
print "r[",j,"]=",rr[j],",s[",j,"]=",s[j],",t[",j,"]=",t[j],"\n"
}
print " has length ",l,"\n"
}
if(l%2==0){
c=l/2
a=abs(s[c]);b=abs(s[c+1]);r=abs(t[c]);s=abs(t[c+1]);
print "(a,b,r,s)=(",a,",",b,",",r,",",s,")\n"
print "x=ar+bs=",a*r+b*s,",br-as=",b*r-a*s,",",r,"^2+",s,"^2=",m,"\n"
}
return(l)
}
/* n divides x^2+1, 1<x<n/2. Then n is the sum of relatively prime squares n=r^2+s^2, r<s.
x=ar+bs, 0<=a<=b,1<r<s, br-as=\eps= 1 or -1 and a<=r/2, b<=s/2.
Also xr=s\eps(mod n). We get a=0 only for the representation n=1^2+s^2.
Returns the length 2i of the euclidean algorithm performed on (r0,r1)=(n,x). We exit when the first
remainder r_c <=sqrt(n), c>=1.
See http://www.numbertheory.org/pdfs/hermite_serret.pdf
*/
define euclid1(n,x){
auto a,b,c,i,r,s,z,e,r0,r1,r2
z=sqrt(n)
r0=n
r1=x
i=1
if(x<=z){
s=x
r=r0%r1
a=(x*r-s)/n
b=(r+s*x)/n
print "(a,b,r,s)=(",a,",",b,",",r,",",s,")\n"
print "x=ar+bs=",a*r+b*s,",br-as=",b*r-a*s,",",r,"^2+",s,"^2=",n,"\n"
return(2)
}
r2=r0%r1
while(r2>0){
i=i+1
if(r2<=z){
s=r2
r=r1%r2
e=parity(i+1)
a=(x*r-e*s)/n
b=(r*e+s*x)/n
print "(a,b,r,s)=(",a,",",b,",",r,",",s,")\n"
print "x=ar+bs=",a*r+b*s,",br-as=",b*r-a*s,",",r,"^2+",s,"^2=",n,"\n"
return(2*i)
}
r0=r1
r1=r2
r2=r0%r1
}
}
define squares1(n){
auto i,l
r=sqroot(-1,n,1)
print r,"\n"
if(r==0){
print "no solution of x^2=-1 (mod n)\n"
return
}
for(i=0;i<numbr;i++){
l=euclid1(n,solution[i])
print "length=",l,"\n"
}
}
define parity(e){
if(e%2){
return(-1)
}else{
return(1)
}
}