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

/* From Otto Forster, Algorithmsche Zahlentheorie, p.30. */
define gcd_coeff(x,y){
auto temp,q,q21,q22,t21,t22,signx,signy

  signx=1; signy=1
  if(x<0){
     signx=-1
     x=-x
  }
  if(y<0){
    signy=-1
    y=-y
  }
  q11=q22=1; q12=q21=0
  while(y!=0){
     temp=y
     q=int(x,y)
     y=mod(x,y)
     x=temp
     t21=q21; t22=q22
     q21=q11-q*q21
     q22=q12-q*q22
     q11=t21; q12=t22
  }
  q11global=signx*q11; q12global=signy*q12
  return(x)
}

/* 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,temp
	t=n
	n=n-1
	for(i=0;i<=n;i++) b[i]=m[i]
	for(i=1;i<=n;i++) b[i]=gcd_coeff(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]
                temp=gcd_coeff(c,d)
		b[i]=q11global*k
		if(i==1) break
                k=k*q12global
		b[i-1]=m[i-1]
	}
        temp=gcd_coeff(c,d)
	b[0]=k*q12global
	for(i=0;i<n;i++){print "(",b[i],")*",m[i],"+"}
	for(i=0;i<n;i++){bglobal[i]=b[i]}
	print "(",b[n],")*",m[n]
        bglobal[n]=b[n]
	print " = ",a,"\n"; return(a)
}

/* lcm(a,b) for any integers a and b */

define lcm(a,b){
	auto g
	g=gcd_coeff(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,l
        l=m[0]
	for(i=1;i<n;i++){
            l=lcm(l,m[i])
        }
	return(l)
}

/* 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.
 */
 
define chinese_remainder_ore(a[],m[],n){
auto l,i,j,mmi[],g,sum,temp,t,d,dij
     t=n-1  
     for(j=1;j<=t;j++){
        for(i=0;i<j;i++){
           dij=gcd_coeff(m[i],m[j])
           d=mod(a[i]-a[j],dij)
           if(d!=0){
             print "no solution\n"
             return(-1)
           }
        }
     }
     l=lcma(m[],n)
     for(i=0;i<n;i++){
        mmi[i]=l/m[i]
     }
     g=gcda(mmi[],n)
     sum=a[0]*bglobal[0]*mmi[0]
     for(i=1;i<n;i++){
        temp=a[i]*bglobal[i]*mmi[i]
        sum=sum+temp
     }
     sum=mod(sum,l)
     return(sum)
}