/* 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)
}
}