/* gnubc program : surd */
"
This program finds the continued fraction expansion of a
quadratic irrational (u+tsqrt(d))/v, where d is not a square,
t,u,v integers, v nonzero.
After converting to (U+sqrt(D))/V, where V divides D-U^2,
typing surd(d,t,u,v,1) prints a[i] and (P[i[+sqrt(D))/Q[i],
typing surd(d,t,u,v,0) prints only the a[i].
"
/* sign of an integer a */
/* sign(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)
}

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

/*
 * This is a function for finding the period of the continued fraction
 * expansion of reduced quadratic irrational a=(u+sqrt(d))/v.
 * Here d is non-square, 1<(u+sqrt(d))/v, -1<(u-sqrt(d))/v<0.
 * The algorithm also assumes that v divides d-u*u and is based on K. Rosen,
 * Elementary Number theory and its applications, p.379-381 and Knuth's The art
 * of computer programming, Vol. 2, p. 359.  0 is returned if a is not reduced.
 */
define period(d,u,v,i,flag){
	auto a,r,s,t,f,j,k,temp
	f=sqrt(d)
	if(u<=0){"a is not reduced: u<=0:";return(0)}
	if(v<=0){"a is not reduced: v<=0:";return(0)}
	if(u+f<v){"a is not reduced: a<1:";return(0)}
	if(u+v<= f){"a is not reduced: a'<=-1:";return(0)}
	if(f<u){"a is not reduced: a'>0:";return(0)}
	k=d-u*u
	if(k%v!=0){"v does not divide d - u*u:";return(0)}
	s=v
	r=u
/* i is created by surd(d,t,u,v,flag) below and indexes the ith 
   convergent of (u+t*sqrt(d))/v */
	print "period:\n"
	for(j=i;1;j++){
		a=(f+u)/v
                print "a[",j,"]=",a,"\n"
		u=a*v-u
		temp=v
		v=(d-u*u)/v
		if(flag){
	           print "        u[",j+1,"]=",u, ", v[",j+1,"]=", v,"\n"
		}
		if(u==r){if(v==s)return(j+1-i)}
/* the continued fraction has period length j+1-i */
	}
}

/*
 * This function uses the continued fraction algorithm expansion in K. Rosen,
 * Elementary Number theory and its applications,p.379-381 and Knuth's
 * The art of computer programming, Vol. 2, p. 359. It locates the first reduced
 * complete quotient and then uses the function period(d,u,v,i,flag)
 * to locate the period of the continued fraction expansion.
 */

define surd(d,t,u,v,flag){
	auto a,c,e,g,f,i,j,w,z
	c=d;g=u;e=v
	if(v==0){"v=";return(0)}
	f=sqrt(d)
	if(f*f==d){"d is the square of ";return(f)}
	if(t==0){"t = ";return(0)}
	z=sign(t)
	u=u*z
	v=v*z
	d=d*t*t
	w=abs(v)
	if((d-u*u)%w!=0){
		d=d*v*v
		u=u*w
		v=v*w
	}
	f=sqrt(d)
	if(flag){print "D=",d," "}
	for(j=0;1;j++){
		if(flag){
			print "      u[",j,"]=",u, ", v[",j,"]=", v,"\n" 
		}
		if(v>0){
			if(u>0){
				if(u<= f){
					if(f<u+v){
						if(v-u <=f){
				 /* (u+sqrt(d))/v is reduced */
							i=j
							/*if (i)
								print ";"*/
							l=period(d,u,v,i,flag)
print "\nThe continued fraction for "
if(t>1) {
	if(e^2!=1){
		print "(",g,"+",t,"*sqrt(",c,"))/",e
	}else{
		if(e==1){
			print g,"+",t,"*sqrt(",c,")"
		}else{
			print "-(",g,"+",t,"*sqrt(",c,"))"
		}
	}
}
if(t==1) {
	if(e^2!=1){
		print "(",g,"+sqrt(",c,"))/",e
	}else{
		if(e==1){
			print g,"+sqrt(",c,")"
		}else{
			print "-(",g,"+sqrt(",c,"))"
		}
	}
}
if(t== -1) {
	if(e^2!=1){
		print "(",g,"-sqrt(",c,"))/",e
	}else{
		if(e==1){
			print g,"-sqrt(",c,")"
		}else{
			print "-(",g,"-sqrt(",c,"))"
		}
	}
}
if(t< -1) {
	if(e^2!=1){
		print "(",g,t,"*sqrt(",c,"))/",e
	}else{
		if(e==1){
			print g,t,"*sqrt(",c,")"
		}else{
			print "-(",g,t,"*sqrt(",c,"))"
		}
	}
}
print " has period length ";return(l)
						}
					}
				}
			}
		}
		if(v>0)a=int(f+u,v)
		if(v<0)a=int(f+u+1,v)
                print "a[",j,"]=",a,"\n"
		u=a*v-u
		v=int(d-u*u,v)
	}
}