/* gnubc program: gcd */
/*
		NUMBER THEORY ALGORITHMS
1.  sign(n);
2.  abs(n);
3.  mod(a,b), (b > 0): returns a(mod b);
4.  int(a,b):  returns the integer part of a/b, b nonzero;
5.  gcd(m,n);
6.  gcd1(m,n): gcd(m,n)=gcd1(m,n)m+gcd2(m,n)n;
7.  gcd2(m,n);
8.  lcm(m,n);
9.  inv(a,n): returns the inverse of a mod m;
10. cong(a,b,m): solves ax=b (mod m);
11. chinese(a,b,m,n): solves x=a(mod m),x=b(mod n);
12. gcda(m[],n): finds gcd(m[0],...,m[n-1]) and expresses it
    as a linear combination of m[0],...,m[n-1];
13. lcma(m[],n): finds lcm(m[0],...,m[n-1]);
14. chinesea(a[],m[],n): solves x=a[i](mod m[i]), i=1,...,n-1;
15. chineseab(a[],b[],m[],n): solves b[i]*x=a[i](mod m[i]),
    i=1,...,n-1;
16. mpower(a,b,c): returns a^b (mod c);
17. exp(a,b): returns a^b;
18. mthroot(a,b,m): returns the integer part of the mth-root of a/b. 
    Here a,b and m are positive integers, m > 1.
19. mthrootr(a,b,m,r): this gives the mth-root of a/b truncated
    to r decimal places, r>=0.
20. min(x,y): this returns x if x<=y, y if x>y.
21. max(x,y): this returns y if x<=y, x if x>y.
22. gcd3(a,b,c): returns gcd(a,b,c)
23. ceiling(a,b): returns least integer x>= a/b.
24. bezout(a,b): here a>=0 and b>=0.  This returns d=gcd(a,b) and global variables
    globalu and globalv, where d = globalu.a + globalv.b.
25. bezout1(a,b): here a and b are of arbitrary sign.  This returns d=gcd(a,b) and global variables
    globalu and globalv, where d = globalu*a + globalv*b.
26. testarrayvalue(arrayvalue[],count,value).
27. nearestsquare(d).
*/
/* sign of an integer a */
/* s(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)
}

/*NOTE: in bc we have */
/* a%b=m(a,b) if a>=0 or a<0 and b divides a */
/* a%b=m(a,b)-b if a<0, b>0, a not divisible by b */
/* a/b=[a/b] if a>=0 or a<0 and b divides a */
/* a/b=[a/b]+1 if a<0, b>0, a not divisible by b */
/* a=b(a/b)+a%b */

/* 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){
	if(b<0){
	a= -a
	b= -b
	}
	return((a-mod(a,b))/b)
}

/*  gcd(m,n) for any integers m and n */
/* Euclid's division algorithm is used. */
/* We use gcd(m,n)=gcd(m,|n|) */

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

/* gcd(a,b)=a*h(a,b)+b*k(a,b) for any integers a and b,
 * where h(a,b) and k(a,b) are calculated as follows:
 * Let q[1],...,q[n] be the quotients in Euclid's algorithm.
 * Then the recurrence relation 
 * s[0]=1, s[1]=0, s[j]=s[j-2]-q[j-1]*s[j-1] for j=2,...,n
 * is implemented. Then h(a,b)=s[n].
 * The recurrence relation 
 * t[0]=0, t[1]=1, t[j]=t[j-2]-q[j-1]*t[j-1] for j=2,...,n
 * is implemented. Then k(a,b)=t[n]. 
 */

/* gcd1(m,n)=h(m,n).
 * gcd2(m,n)=k(m,n).
 */

define gcd1(m,n){
	auto a,b,c,h,k,l,q,t
	if(n==0) return(sign(m))
        a=m                /* a=r[0] */ 
        b=abs(n)           /* b=r[1] */ 
        c=mod(a,b)         /* c=r[2] */ 
	if(c==0) return(0)
        l=1                /* l=s[0] */ 
        k=0                /* k=s[1] */ 
        while(c>0){
                q=(a-c)/b  /* q=q[j-1]=[r[j-2]/r[j-1]] */
                a=b
                b=c
                c=mod(a,b) /* c=r[j],r[j-2]=q[j-1]*r[j-1]+r[j]*/
                h=l-q*k    /* h=s[j],k=s[j-1], l=s[j-2] */ 
              	l=k 
              	k=h
        }         
	return(k)          /* k=s[j] */
}   


define gcd2(m,n){
        auto a,b,c,h,k,l,q,t,j
	if(n==0) return(0)
        a=m                /* a=r[0] */ 
        b=abs(n)           /* b=r[1] */ 
        c=mod(a,b)         /* c=r[2] */ 
	if(c==0) return(sign(n))
	l=0                /* l=t[0] */ 
	k=1                /* k=t[1] */ 
/*	print "     ","   r[0]=",a," t[0]=",0,"\n"
	print "q[1]=",(a-c)/b,"  r[1]=",b," t[1]=",1,"\n" */
/*j=2*/
        while(c>0){
                q=(a-c)/b  /* q=q[j-1]=[r[j-2]/r[j-1]] */
                a=b
                b=c
                c=mod(a,b) /* c=r[j],r[j-2]=q[j-1]*r[j-1]+r[j]*/
                h=l-q*k    /* h=t[j],k=t[j-1], l=t[j-2] */ 
/*              j=j+1
		print "q[",j-1,"]=",(a-c)/b,"  r[",j-1,"]=",b," t[",j-1,"]=",h,"\n" */
               	l=k
               	k=h
        }
		q=(a-c)/b
/*		print "        r[",j,"]=",c," t[",j,"]=",l-q*k,"\n" */
	if(n<0)	return(-k)
	return(k)          /* k=t[j] */
}  

/* 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 congruence mx=p(mod n) */

define cong(m,p,n){
        auto a,b,t,x,y,z
	a=gcd(m,n)
	if(mod(p,a)>0){
 "the number of solutions is "; return(0)
}
	b=gcd1(m,n)
	y=n/a
	p=p/a
	z=mod(b*p,y)
	print "The solution is x ="
	for(t=0;t<a;t++){print " ",z+t*y,","}
	print " mod ",n,"\n"
	print "or equivalently, x = ",z," mod "; return(y)
}

/* the congruence mx=p(mod n)
 * with printout suppressed-used in general form of the Chinese
 * remainder  theorem .
 */
define cong1(m,p,n){
        auto a,b,x,y
	a=gcd(m,n)
	if(mod(p,a)>0){return(-1)}
	b=gcd1(m,n)
	y=n/a
	p=p/a
	x=mod(b*p,y)
	return(x)
}

/* 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.
 */

define chinese(a,b,m,n) {
	auto c,d,x,y,z
	d = gcd(m,n)
	if(mod(a-b,d)!=0){
		"the number of solutions is "
		return(0)
	}
	x= m/d;y=n/d
	z=(m*n)/d
	c = mod(b*gcd1(x,y)*x + a*gcd2(x,y)*y,z)
	print "the solution mod(",z,") is "; return(c)
}

/* gcd(m[0],m[1],...,m[n-1])=b[0]*m[0]+...+b[n-1]*m[n-1] */

define gcda(m[],n){
	auto a,i,c,d,k,b[],t
	t=n
	n=n-1
	for(i=0;i<=n;i++) b[i]=m[i]
	for(i=1;i<=n;i++) b[i]=gcd(b[i],b[i-1])
	a=b[n]
	b[n]=m[n]
	k=1
	for(i=n;i>=1;i--) {
		c=b[i];d=b[i-1]
		b[i]=gcd1(c,d)*k
		if(i==1) break
		k=k*gcd2(c,d)
		b[i-1]=m[i-1]
	}
	b[0]=k*gcd2(c,d)
	for(i=0;i<n;i++){print "(",b[i],")*",m[i],"+"}
	print "(",b[n],")*",m[n]
	print " = "; return(a)
}

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

/* lcm(m[0],m[1],...,m[n-1]) */

define lcma(m[],n){
	auto i,b[]
	for(i=0;i<n;i++) b[i]=m[i]
	for(i=1;i<n;i++) b[i]=lcm(b[i],b[i-1])
	return(b[n-1])
}

/*
 * a^b (mod c), a,b,c integers, a,b>=0,c>=1
 */

define mpower(a,b,c){
	auto x,y,z
	x=a
	y=b
	z=1
/*	print z," ",x," ",y,"\n"*/
	while(y>0){
		while(y%2==0){
			y=y/2
			x=(x*x)%c
		}
		y=y-1
		z=mod(z*x,c)
	}
	return(z)
}

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

/* We use a discrete form of Newton's method to find [x<sup>1/m</sup>], 
 * the integer part of the mth root of a positive integer x. 
 * (See K.R. Matthews, Computing m-th roots, The College Mathematics Journal 19
 * (1988) 174-176 and H. Lüneburg's book On the Rational Normal Form of 
 * Endomorphisms 1987, B.I. Wissenschaftsverlag, Mannheim/Wien/Zürich.)
 * Note that if x is a positive rational, then [x^(1/m)]=[[x]^(1/m)].
 * The answer is given truncated to r decimal places.
 */
define mthrootr(a,b,m,r){
	auto s,x,y,i,z,d,f,t
	f=exp(10,r)
	a=a*exp(f,m)
	a=a/b
	if(a==0){
		if(r==0){
			"answer = ";return(0)
		}
		"answer = 0."
		for(i=1;i<r;i++)"0"
		return(0)
	}
	s=a
	while(s>1){
		s=s/2
		i=i+1
	}
	x=exp(2,1+(i/m))
	while(1==1){
		y=x
		z=exp(x,m-1)
		x=x+int(a-x*z,m*z)
		if(y<=x)break
	}
	if(r==0){
		 "answer = ";return(y)
	}
	d=y%f
	print "answer = ",y/f, "."
	t=r-length(d)
	for(i=1;i<=t;i++){
          print "0"
	}
         return(d)
}

/* We use a discrete form of Newton's method to find [x<sup>1/m</sup>], 
 * the integer part of the mth root of a positive integer x. 
 * (See K.R. Matthews, Computing m-th roots, The College Mathematics Journal 19
 * (1988) 174-176 and H. Lüneburg's book On the Rational Normal Form of 
 * Endomorphisms 1987, B.I. Wissenschaftsverlag, Mannheim/Wien/Zürich.)
 * Note that if x is a positive rational, then [x^(1/m)]=[[x]^(1/m)].
 */
define mthroot(a,b,m){
	auto s,x,y,i,z,d,j
	a=a/b
	if(a==0)return(0)
	s=a
	while(s>1){
		s=s/2
		i=i+1
	}
	x=exp(2,1+(i/m))
	/*print "x[0]=",x,"\n"*/
	j=1
	while(1==1){
		y=x
		z=exp(x,m-1)
		x=x+int(a-x*z,m*z)
		/*print "x[",j,"]=",x,"\n"*/
		j=j+1
		if(y<=x) break
	}
	return(y)
}

/* From the bc manual. */
/*
 * b(n,m) returns the binomial coefficient, the number of subsets each 
 * containing m elements from a set containing n elements.
 */

define binomial(n,m){
auto x,j
if(n<m)return(0)
if(m==0)return(1)
x=1
for(j=1;j<=m;j=j+1) x=x*(n-j+1)/j
return(x)
}

/* min(x,y) */

define min(x,y){
        if(y<x)return(y)
        return(x)
}

/* max(x,y) */

define max(x,y){
        if(y<x)return(x)
        return(y)
}

/* the Chinese remainder theorem for n simultaneous congruences
 * x=a[i] (mod m[i]), i=0,1,...,n-1, where m[i]>0 for all i and
 * a[0],...,a[n-1] are arbitrary integers.
 * the construction of O. Ore, American Mathematical Monthly,
 * vol.59,pp.365-370,1952, is implemented.
 */

/* 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(), but with printing suppressed, except that -1
 * is returned if here is no solution.
 */

define chinese1(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 chinesea(a[],m[],n){
	auto i,m,x
	m=m[0]
	x=a[0]
	for(i=1;i<n;i++){
		x=chinese1(a[i],x,m[i],m)
		if(x==-1)break
		m=lcm(m[i],m)
	}
	if(x==-1){
		print "there is no solution\n"
	}\
	else{
                print "HERE1\n"
		print "the solution mod(",m,") is "
                print "HERE2\n"
                return(x)
	}
}

/* the Chinese remainder theorem for n simultaneous congruences
 * a[i]x=b[i] (mod m[i]), i=0,1,...,n-1, where m[i]>0 for all i.
 * a[0],...,a[n-1] and b[0],...,b[n-1] are arbitrary integers.
 */

define chineseb(a[],b[],m[],n){
        auto c[],i,m,s,x
	m=m[0]
	x=cong1(a[0],b[0],m[0])
	for(i=1;i<n;i++){
		s=cong1(a[i],b[i],m[i])
		if(s==-1){x=-1;break}
		c[i]=s
                x=chinese1(c[i],x,m[i],m)
                if(x==-1)break
                m=lcm(m[i],m)
        }
        return(x)

}

define gcd3(a,b,c){
auto t
t=gcd(a,b)
t=gcd(t,c)
return(t)
}

/* ceiling function */
define ceiling(a,b){
auto x
x=int(a,b)
if(b*x==a){return(x)}\
else{return(x+1)}
}

/* From Henri Cohen' book, Alg. 1.3.6  13/07/2011 
 * This assumes a >=0 and b >= 0. 
 * returns d= gcd(a,b) and global variables globalu and globalv,
 * where d = globalu.a + globalv.b.
 */
define bezout(a,b){
auto d,v1,v3,t1,t3,q
   globalu=1
   d=a
   if(b==0){
      globalv=0
      return(a)
   }else{
      v1=0
      v3=b
   }
   while(v3>0){
      q=d/v3
      t3=d%v3
      t1=globalu-q*v1
      globalu=v1
      d=v3
      v1=t1
      v3=t3
   }
   globalv=(d-a*globalu)/b
   return(d)
}

/* Here a and b are any integers.
 * returns d= gcd(a,b) and global variables globalu and globalv,
 * where d = globalu.a + globalv.b.
 */
define bezout1(a,b){
auto u,v,absa,absb,d

   if(a<0){
     absa=-a
   }else{
     absa=a
   }
   if(b<0){
     absb=-b
   }else{
     absb=b
   }
   d=bezout(absa,absb)
   globalu=globalu*sign(a)
   globalv=globalv*sign(b)
   return(d)
}

define squaretest(n){
auto g
g=sqrt(n)
if(g*g==n){
   return(g)
}else{
   return(-1)
}
}

define minarray(a[],n){
auto l,k
l=0
for(k=1;k<n;k++){
   if(a[k]<a[l]){
     l=k
   }
}
return(a[l])
}

# this prints the AP an+b when a > 0.
define printexponent(a,b){
auto abs
  absb=abs(b)
  if(a==1){
    if(b==0){
      print "n"
    }
    if(b>0){
      print "n + ",b
    }
    if(b<0){
      print "n - ",absb
    }
  }
  if(a>1){
    if(b==0){
      print a,"n"
    }
    if(b>0){
      print a,"n + ",b
    }
    if(b<0){
      print a,"n -",absb
    }
  }
}

#define testarrayvalue(arrayvalue[],count,value){
#  flag=1
#  for(h=0;h<count;h++){
#     if(i$arrayvalue[$h]!=$value){
#       flag=0
#       break
#     }
#  }
#  return(flag)
#}

define rationalization(a,b,c,d,dd){
  auto aa,bb,cc,g
 aa=a*c-b*d*dd
 bb=b*c-a*d
cc=c^2-d^2*dd
g=gcd3(aa,bb,cc)
aa=aa/g
bb=bb/g
cc=cc/g
if(cc<0){
aa=-aa
bb=-bb
cc=-cc
}
print "(",aa,"+",bb,"sqrt(",dd,"))/",cc,"\n"
}

define nearestsquare(d){
auto x,y,z,t

   x=sqrt(d)
   y=d-x^2
   t=x+1
   z=t^2-d
   if(y<z){
     global_nearest_square_remainder=y
     return(x^2)
   }else{
     global_nearest_square_remainder=-z
     return(t^2)
   }
}