/* bc program padic */
 * mult(a[],b[],m,n,p) forms the product of a=a[0]+a[1]p+...+a[m]p^m 
 * and b=b[0]+b[1]p+...+b[n]p^n
 * where a[m] and b[n] are nonzero.
 * The product is given by the global variable prod[]
 * prod[0]+prod[1]p+...+prod[l]p^l, where prod[l] is nonzero.
 * If a or b is zero, we return 0 at the start, otherwise return l..
 * The program is an adaption of one in i1.c
 * from http://www.numbertheory.org/calc/
 */
define mult(a[],b[],m,n,p){
	auto c,i,t,j,bk,k,l
	if((m==0 && a[0]==0) || (n==0 && b[0]==0)){
		prod[0]=0
		print prod[0]
		l=0
		print "\n"
		return(l)
	}
        c=0
	l=m+n+1
        for(k=0;k<=m;k++){
		prod[k]=0
	}
        for(k=0;k<=n;k++){
		bk=b[k]
                if(bk){
			c=0
        		for(j=0;j<=m;j++){
				t=a[j]*bk+prod[j+k]+c
				prod[j+k]=t%p
				c=t/p
			}
			prod[m+k+1]=c
                }else{
			prod[m+k+1]=0
		}
	}
	if(c==0){
		l=m+n
	}
	return(l)
}

/*
 * rsv(a[],b[],m,n) is an adaption of one in i1.c
 * from http://www.numbertheory.org/calc/
 * it returns 1,0,-1 according as a[] >,=,< b[].
 */
define rsv(a[],b[],m,n){
	auto j
	if(m>n){
		return(1)
	}
	if(m<n){
		return(-1)
	}
	j=m
	while((a[j]==b[j]) && j>0){
		j=j-1
	}
	if(a[j]>b[j]){
		return(1)
	}
	if(a[j]<b[j]){
		return(-1)
	}
	return(0)
}

/*
 * sub(a[],b[],m,n,p) is an adaption of one in i1.c
 * from http://www.numbertheory.org/calc/
 * Here a[]>= b[] and returns the l and the global array 
 * diff[]=a[]-b[]=diff[0]+diff[1]*p+...+diff[l]*p^l.
 */
define sub(a[],b[],m,n,p){
	auto c,d,f,k,j,l,t
	if(n==0 && b[0]==0){
	      for(j=0;j<=m;j++){
                 diff[j]=a[j]
	      }
              l=m
	      return(l)
	}
	c=1
	l=m
        for(j=0;j<=n;j++){
		t=a[j]-b[j]+p+c-1
                diff[j]=t%p
                c=t/p
        }
	if(m>n){
		for(j=n+1;j<=m;j++){
			t=a[j]+p+c-1
			diff[j]=t%p
			c=t/p
		}
	}
	d=0
	k=m
	while (k >= 0)
        {
                if (diff[k] != 0)
                {
                        d = k
                        k = -1
                }else{
                        k=k-1
		}
        }
        f = m - d;
	if(f){
		l=d
	}
	return(l)
}

/* base(b,n) gives the base b expansion of n:
 * base[]=baseb[0]+baseb[1]*b+ ...+baseb[i]*b^i.
 * The integer i is returned, along with the global
 * variable baseb[].
 */
define base(b,n){
	auto i,q,t,x
	i=0
	if(n<0){"Try again, n<";return(0)}
	if(b<=1){"Try again, b<=";return(1)}
	x=n
	while(n>=b){
		q=n/b
		t=n-q*b
		baseb[i]=t
		n=q
		i=i+1
	}
	baseb[i]=n
	return(i)
}
/*
 * twoadicsqrt(a,n) n>=3, returns the first n binary digits of a
 * 2adic sqroot x of a positive integer a=8k+1. Here x=1 or 5 (mod 8).
 * See http://www.numbertheory.org/courses/MP313/lectures/lecture23/page3.html
 * and http://www.numbertheory.org/courses/MP313/solns/soln3/page1.html.
 */
define twoadicsqrt(a,n){
	auto i,j,k,l,x,y,z,b[]
	if(n<1){
		print "n<1\n"
                return
	}
	if(a%8!=1){
		print "a is not of the form 8k+1\n"
                return
	}
	if(a<1){
		print "a<1\n"
                return
	}
	b[0]=1; 
	if(n==1){
          print "1 is the first digit of the 2adic square root of "
          return(a)
        }
	b[1]=0; 
	if(n==2){
          print "10 are the first digits of the 2adic square root of "
          return(a)
        }
	x=base(2,a)
        l=0
	for(k=3;k<=n;k++){
	   y=mult(b[],b[],l,l,2)	
	   z=rsv(prod[],baseb[],y,x)
	   if(z<0){
		j=sub(baseb[],prod[],x,y,2)
	   }else{
		j=sub(prod[],baseb[],y,x,2)
	   }
	   if(diff[k]){
              b[k-1]=1
              if(k==3){
                l=2
              }else{
                  l=k-1
              }
           }
	}
	for(i=0;i<n;i++){
           if(i && i%20==0){
              print "\n"
           }
           print b[i]
        }
        print "\n"
        print "are the first ",n," digits of a 2adic square root of "
        return(a)
}

/*
 * padicsqrt(a,p,n) p an odd prime, returns the first n p-adic digits of a
 * padic sqroot x of a positive integer a. Here x=a (mod p).
 * See http://www.numbertheory.org/courses/MP313/lectures/lecture23/page3.html
 * and http://www.numbertheory.org/courses/MP313/solns/soln3/page3.html.
 */
define padic(a,p,n){
auto j,l,k,u,x,y,z,b[],sign
	u=tonelli(a,p)
	print "u=",u,"\n"
	x=base(p,a)
	sign=1
        l=0
	b[0]=u
	for(k=1;k<=n;k++){
	   y=mult(b[],b[],l,l,p)	
	   z=rsv(prod[],baseb[],y,x)
	   if(z<0){
		sign= -1
		j=sub(baseb[],prod[],x,y,p)
	   }else{
		j=sub(prod[],baseb[],y,x,p)
	   }
	   if(diff[k]){
              b[k]=cong1(sign*2*u,-diff[k],p)
           }else{
		b[k]=0
	   }
	   l=k
	}
	for(i=0;i<n;i++){
           if(i && i%20==0){
              print "\n"
           }
           print b[i]
        }
        print "\n"
        print "are the first ",n," digits of a padic square root of "
        return(a)
}

/* a[]~a[0]+a[1]p+...+a[n-1]p^(n-1)=a
 * b[]~b[0]+b[1]p+...+b[n-1]p^(n-1)=b
 * prodp[]~prodp[0]+prodp[1]p+...+prodp[n-1]p^(n-1)=prodp
 * where prodp=a*b(mod p^n)
 */
define multp(a[],b[],n,p){
	auto c,i,t,j,bk,k,flag
	flag=0
	for(i=0;i<n;i++){
		if(a[i]!=0){
			flag=1
			break
		}
	}
	if(flag==0){
		for(i=0;i<n;i++){
			prodp[i]=0
		}
		return
	}
	flag=0
	for(i=0;i<n;i++){
		if(b[i]!=0){
			flag=1
			break
		}
	}
	if(flag==0){
		for(i=0;i<n;i++){
			prodp[i]=0
		}
		return
	}
        c=0
        for(k=0;k<n;k++){
		prodp[k]=0
	}
        for(k=0;k<n;k++){
		bk=b[k]
                if(bk){
			c=0
        		for(j=0;j<n;j++){
				t=a[j]*bk+prodp[j+k]+c
				prodp[j+k]=t%p
				c=t/p
			}
			prodp[n+k+1]=c
                }else{
			prodp[n+k+1]=0
		}
	}
	return
}

/* a[]~a[0]=a[1]p+...a[n]p^n, 1<b<p, 0<=a[i]<p.
 * ratio[]~ratio[0]+ratio[1]+...+ratio[l], where ratio=int(a/b).
 * l is returned.
 */
define divp(a[],n,b,p){
        auto i,k,l,t
	if(p==2){
		if(n==0){
			ratiop[0]=0
			return(0)
		}
		for(i=0;i<n;i++){
			ratiop[i]=a[i+1]
		}
		return(n-1)
	}
	for(i=0;i<=n;i++){
		ratiop[i]=a[i]
	}
     	t = ratiop[n]
	l=n
     if(n!=0){
          for(k=n; k >= 1; k--){
               ratiop[k]=t/b
               t  = ratiop[k-1] + ((t%b)*p)
          }
          if (ratiop[n]==0){
		l=n-1
          }
     }
     ratiop[0] = t / b
     return(l)
}

define parity(a[],n,p){
	auto i,s

	if(p==2){
		t=a[0]%2
		return(t)
	}
	s=0
	for(i=0;i<=n;i++){
		t=a[i]%2
		s=s+t
	}
	return(s%2)
}

/* a[]~a[0]+a[1]p+...+a[k-1]p^(k-1) (mod p^k), 0<=a[i]<p.
 * b[]~b[0]+b[1]p+...+b[m]p^m, 0<=b[i]<p.
 * exponentiationp[]~exponentiationp[0]+exponentiationp[1]p+...exponentiationp[k-1]p^(k-1) is
 * a[]^b[](mod p^k).
 */
define mpowerp(a[],b[],m,k,p){
	auto i,x[],y[],z[],one[],u
	one[0]=1
	for(i=0;i<k;i++){
		x[i]=a[i]	
	}
	for(i=0;i<=m;i++){
		y[i]=b[i]
	}
	for(i=0;i<k;i++){
		z[i]=0
	}
	z[0]=1
	while(m>0 || (m=0 && y[0]>0)){
		t=parity(y[],m,p)
		while(t==0){
			m=divp(y[],m,2,p)
			for(i=0;i<=m;i++){
				y[i]=ratiop[i]
			}
			t=parity(y[],m,p)
			u=multp(x[],x[],k,p)
			for(i=0;i<k;i++){
				x[i]=prodp[i]	
			}
		}
		m=sub(y[],one[],m,0,p)
		for(i=0;i<=m;i++){
			y[i]=diff[i]
		}
		u=multp(z[],x[],k,p)
		for(i=0;i<k;i++){
			z[i]=prodp[i]	
		}
	}
	for(i=0;i<k;i++){
		exponentiationp[i]=z[i]	
	}
	return
}