/* gnubc program: lra */

"





			LEAST REMAINDER  ALGORITHM







If lra(m,n) is typed, the quotients and remainders for the
least remainder algorithm for m/n are printed, starting with 
m=n*q+r. 
The length of the algorithm is printed. (Type quit to quit)
If nicf(m,n) is typed, the nearest integer continued fraction of m/n
m/n=a_0-1/a_1-...-1/a_n is printed as (a_0,a_1,...,a_n)
"
/* lmodd(m,n) returns r, m=q*n+r, -n/2 < r <= n/2 */
define lmodd(m,n){
	return(m-n*lnearint(m,n))
}

/* rmodd(m,n) returns r, m=q*n+r, -n/2 <= r < n/2 */
define rmodd(m,n){
	return(m-n*rnearint(m,n))
}

/* lnearint(m,n) returns z, the nearest integer to m/n,  
 * z=t if m/n=t+1/2
 */
define lnearint(m,n){
	auto x,y,z
        if(n<0){
          m= -m
          n= -n
        }
	y=int(m,n)
	x=n*y
	z=2*(m-x)
	if (z>n){y=y+1}
	return(y)
}

/* rnearint(m,n) returns z, the nearest integer to m/n,  
 * z=t+1 if m/n=t+1/2
 */
define rnearint(m,n){
	auto x,y,z
	y=int(m,n)
	x=n*y
	z=2*(m-x)
	if (abs(z)>= abs(n)){y=y+1}
	return(y)
}

define lra(m,n){
	auto a,b,r,q,i
	a=m
	b=n
	q=lnearint(a,b)
	r=a-b*q
	print "\n"
	print a," = ",b,"*",q,"+",r,"\n"
	i=1
        while(r!=0){
		i=i+1
                a=b
                b=r
		q=lnearint(a,b)
		r=a-b*q
		print a," = ",b,"*",q,"+",r,"\n"
        }
	print "the least remainder algorithm for m/n with m = ",m, " n = ",n
        print " has length "; return(i)
}   


/* Calculates a[0]=a,a[1]=b,...,a[n], where 
 * a[n+1]=rnearint(a[]*a[n]/a[n-1]), a,b positive integers.
 * Mentioned by Edward Brisse, 19/12/97, in an email.
 */
define rbrisse(a,b,n){
auto i,t
print "a[0]=",a,"\n"
if(n==0)return
print "a[1]=",b,"\n"
scale=5
print "a[1]/a[0]=",b/a,"\n"
scale=0
if(n==1)return
for(i=2;i<=n;i++){
	t=rnearint(b^2,a)
	print "a[",i,"]=",t,"\n"
	scale=5
	print "a[",i,"]/a[",i-1,"]=",t/b,"\n"
	scale=0
	a=b
	b=t
}
return
}
/* Calculates a[0]=a,a[1]=b,...,a[n], where 
 * a[n+1]=lnearint(a[]*a[n]/a[n-1]), a,b positive integers.
 * Mentioned by Edward Brisse, 19/12/97, in an email.
 */
define lbrisse(a,b,n){
auto i,t
print "a[0]=",a,"\n"
if(n==0)return
print "a[1]=",b,"\n"
scale=5
print "a[1]/a[0]=",b/a,"\n"
scale=0
if(n==1)return
for(i=2;i<=n;i++){
	t=lnearint(b^2,a)
	print "a[",i,"]=",t,"\n"
	scale=5
	print "a[",i,"]/a[",i-1,"]=",t/b,"\n"
	scale=0
	a=b
	b=t
}
return
}

/* gnubc program: Near Euclid */

"





			EUCLID'S ALGORITHM:







If neuclid(m,n) is typed in (m>0,n>0), 

the quotients q[k], remainders r[k] for the

nearest integer Euclidean algorithm for m/n are printed, 

starting with r[0]=r[1]*q[1]+r[2]. 

Also the s[k] and t[k] are printed where 

 s[0]=1,s[1]=0, s[j]=s[j-2]-q[j-1]*s[j-1],
 t[0]=0,t[1]=1, t[j]=t[j-2]-q[j-1]*t[j-1], j=2,...,n+1

The length of the algorithm is printed. (Type quit to quit)

"

define nicf(m,n){
auto a,b,r,q,i
        if(n==0){
           print "n=0\n"
           return
        }
        a=m
        b=n
	q=lnearint(a,b)
	r=-a+b*q
	print "The nearest integer continued fraction for m/n is (",q,","
	i=1
        while(r){
           a=b
           b=r
           q=lnearint(a,b)
           print q
	   r=lmodd(a,b)
           if(r){
             print ","
           }else{
             print ")\n"
           }
           r=-r
	   i=i+1
        }
        print " and  has ", i," partial quotients\n"
        return
}