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

/*  gcd(m,n) for positive integers m and n */
/* Euclid's division algorithm is used. */

define gcd(m,n){
	auto a,b,c
	
	if(m==0){
		return(abs(n))
	}
	if(n==0){
		return(abs(m))
	}
	a=abs(m)         /* a=r[0] */ 
        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)
}   

define degreepi(a[],m){
	auto i
	for(i=m;i>=0;i--){
		if(a[i])
			break;
	}
	return(i)
}

/* prints polynomial in a[] to screen */
define printp(a[],n)
{
  auto i
  n=degreepi(a[],n)
  if (n>1){ 
	    for (i=n;i>1;i--) {
		if(a[i]^2!=1){
			print a[i]
		}else{
			if(a[i]==-1){
				print "-"
			}
		}
		print "x^",i
		if(a[i-1]>=0){
		"+"
		}
            }
  }
  if (n>0){
	i=1
	if(a[i]^2!=1){
	    print a[1]
	}else{
		if(a[i]==-1){
			print "-"
		}
	}
	print "x"
	if(a[0]>=0){
	 "+";
	}
  }
  print a[0], "\n";
	return(0)
}


define leadpi(a[],m){
	auto i
	for(i=m;i>=0;i--){
		if(a[i])
			break;
	}
	if(i>-1){
		return(a[i])
	}else{
		return(0)
	}
}

define addpi(a[],b[],m,n){
	auto i,j,k
	if(m<=n){
		for(i=0;i<=m;i++){
			api[i]=a[i]+b[i]
		}
		for(i=m+1;i<=n;i++){
			api[i]=b[i]
		}
		dpi=n
	}else{
		for(i=0;i<=n;i++){
			api[i]=a[i]+b[i]
		}
		for(i=n+1;i<=m;i++){
			api[i]=a[i]
		}
		dpi=m
	}
	return(0)
}

define divmodpi(a[],b[],m,n){
	auto i,j,g,gg,aa,d,f[],t
	
	if(m<n){
		for(i=0;i<=m;i++){
			modpi[i]=a[i]
		}
		degmodpi=m
		degdivpi=-1
		divpi[0]=0
		return(0)
	}else{
		g=b[n]
		d=m-n
		gg=g^(d+1)
		for(i=0;i<=m;i++){
			f[i]=gg*a[i]
		}
		for(i=0;i<d;i++){
			divpi[i]=0
		}
		while(1){
			j=degreepi(f[],m)
			if(j<n)
				break
			aa=leadpi(f[],m)
			if(g!=1){
				aa=aa/g
			}
			t=j-n
			divpi[j]=aa
			for(i=n;i>=0;i--){
				f[i+t]=f[i+t]-aa*b[i]
			}
		}
		for(i=0;i<n;i++){
			modpi[i]=f[i]
		}
	}
	t=degreepi(modpi[],n-1)
	degmodpi=t
	degdivpi=d
	return(0)
}

/* Evaluates a polynomial a[n]*x^n + ... + a[0] at x=b */
define evalpi(a[],n,b) 
{
	auto i,val
	val=a[n]
	for(i=n;i>0;i--){
		val=val*b+a[i-1]
	}
	return (val)
}

define derivpi(a[],n){
	auto i,t
	if(n==0){
		derpi[0]=0
		degderpi=-1
		return(0)
	}
	for(i=n-1;i>=0;i--){
		t=i+1
		derpi[i]=t*a[t]
	}
	degderpi=degreepi(derpi[],n-1)
	return(0)
}

define contentpi(a[],n){
	auto i
	g=a[0]
	for(i=1;i<=n;i++){
		g=gcd(g,a[i])
	}
	return(g)
}

define primitivepi(a[],n){
	auto i,g
	if(n){
		g=contentpi(a[],n)
		for(i=0;i<=n;i++){
			primpi[i]=a[i]/g
		}
	}else{
		primpi[0]=a[0]
	}
	return(0)
}

/* here it's assumed m>0 and n>0 */
/* returns deg(gcd(a,b)) */
/* prints the remainders if e=1 */
define gcdpi(a[],b[],m,n,e){
	auto ca,cb,d,i,j,g,t,prima[],primb[],count

	ca=contentpi(a[],m)
	cb=contentpi(b[],n)
	d=gcd(ca,cb)
	g=primitivepi(a[],m)
	for(i=0;i<=m;i++){
		prima[i]=primpi[i]
	}
	if(e){
		print "r[0]:"
		g=printp(prima[],m)
	}
	g=primitivepi(b[],n)
	for(j=0;j<=n;j++){
		primb[j]=primpi[j]
	}
	if(e){
		print "r[1]:"
		g=printp(primb[],n)
	}
	count=2
	while(1){
		g=divmodpi(prima[],primb[],m,n)
		for(i=0;i<=n;i++){
			prima[i]=primb[i]
		}
		m=n
		n=degmodpi
		if(n== -1){
			break
		}
		for(i=0;i<=n;i++){
			primb[i]=modpi[i]
		}
		g=primitivepi(primb[],n)
		for(i=0;i<=n;i++){
			primb[i]=primpi[i]
		}
		if(e){
			print "r[",count,"]:"
			g=printp(primb[],n)
		}
		count=count+1
	}
	return(degreepi(prima[],m))
}

/* n=number of terms in signs[], n>1 */
define sign_count(signs[],n){
	auto i,changes
	j=0
	for(i=0;i<n;i++){
		if(signs[i]==0){
			continue
		}else{
			s[j]=signs[i]
			j=j+1
		}
	}
	changes=0
	for(i=0;i<j;i++){
		if(s[i]*s[i+1]<0){
			changes=changes+1
		}
	}
	return(changes)
}

/* A.G. Akritas, Elements of computer algebra with applications */
/* p. 341 */
/* prints the sturm polynomials if e=1 */

define sturm(a[],m,b,e){
	auto bb[],d,i,j,g,t,prima[],primb[],count,n,changes
	g=derivpi(a[],m)
	n=degderpi
	deggcd=gcdpi(a[],derpi[],m,n,e)
	if(deggcd>0){
		print "a has a repeated root.\n"
		return
	}
	for(i=0;i<=n;i++){
		bb[i]=derpi[i]
	}

	g=primitivepi(a[],m)
	for(i=0;i<=m;i++){
		prima[i]=primpi[i]
	}
	if(e){
		print "sturm[0]:"
		g=printp(prima[],m)
	}
	signs[0]=evalpi(prima[],m,b)

	g=primitivepi(bb[],n)
	for(j=0;j<=n;j++){
		primb[j]=primpi[j]
	}
	if(e){
		print "sturm[1]:"
		g=printp(primb[],n)
	}
	signs[1]=evalpi(primb[],n,b)
	count=2
	while(n>0){
		g=divmodpi(prima[],primb[],m,n)
		for(i=0;i<=n;i++){
			prima[i]=primb[i]
		}
		m=n
		n=degmodpi
		for(i=0;i<=n;i++){
			primb[i]=modpi[i]
		}
		g=primitivepi(primb[],n)
		for(i=0;i<=n;i++){
			primb[i]=-primpi[i]
		}
		if(e){
			print "sturm[",count,"]:"
			g=printp(primb[],n)
		}
		signs[count]=evalpi(primb[],n,b)
		count=count+1
	}
	if(e){
		for(i=0;i<count;i++){
			print signs[i],","
		}
		print "\n"
	}
	changes=sign_count(signs[],count)
	return(changes)
}