rem roots of polynomials - bairstow's method
rem rewrite of apple basic program
rem others at basicprograms.dns2go.com

cls
print
print "real and imaginary roots of polynomials"
print "using bairstow's method"
print "form is y=A0+A1*x+A2*x^2+A3*x^3+...."
print

dim a(22)
dim b(22)
dim e(22)
e1=.0001
e4=1e-10
k1=100

input "degree of polynomial: ";n$
if n$="q" then end
n=val(n$)
print
for i=0 to n
 print "coefficent of A(";i;")= ";
 input " ";a$
 if a$="q" then end
 atmp=val(a$)
 rem stored in array a in reverse order
 a(n-i+1)=atmp
next i

print
print "Roots:"

rem branch for special treatment of 1st and 2nd degree equations
if n<=2
	goto onetwo
end if

a(n+2)=0
n1=2*int((n+1)/2)
for m1=1 to n1/2
	p=1
 	q=1
	for k=1 to k1
  		for l=1 to k1
   			rem store coefficients in array b
    		for i=1 to n1+1
     			b(i)=a(i)
    		next i
			for j=n1-2 to n1-4 step -2
     			for i=1 to j+1
      				b(i+1)=b(i+1)-1*p*b(i)
      				b(i+2)=b(i+2)-1*q*b(i)
     			next i
    		next j
    		r0=b(n1+1)
    		r1=b(n1)
    		s0=b(n1-1)
    		s1=b(n1-2)
    		v0=-1*q*s1
    		v1=s0-1*s1*p
    		d0=v1*s0-1*v0*s1
    		if abs(d0)>=e4
				goto  jumpout
			end if
    		p=p+5
    		q=q+5
   		next l
		label jumpout
  		d1=s0*r1-1*s1*r0
  		d2=r0*v1-1*v0*r1
  		p1=d1/d0
  		q1=d2/d0
  		p=p+p1
  		q=q+q1
  		if abs(r0)>=e1
  			goto  aleap
  		end if	
  		if abs(r1)>=e1
    		goto  aleap
  		end if
  		e(m1)=1
  		goto  bleap
   		label aleap
  		if abs(p1)>=e1
  			goto  cleap
  		end if
  		if abs(q1)>=e1
  			goto  cleap
  		end if
  		e(m1)=2
  		goto  bleap
   		label cleap
  		if p=0
  			goto  dleap
  		end if
  		if abs(p1/p)>=e1
  			goto  eleap
  		end if
   		label dleap
  		if q=0 then
  			goto eleap
  		end if
  		if abs(q1/q)>=e1 then
    		goto eleap
  		end if
  		e(m1)=3
  		goto  bleap
   		label eleap
 	next k
 	e(m1)=4
  	label bleap
 	s=(-1*p)/2
 	t=s*s-1*q
 	if t<0 then
   		goto fleap
 	end if
 	t=t^.5
  	label kleap
 	print
 	print s+t
 	print s-1*t
 	goto  gleap
  	label fleap
 	t=(-1*t)^.5
  	label jleap
 	print
 	print s;" + I * ";t
 	print s;" - I * ";-1*t
  	label gleap
 	if e(m1)=4
  		goto endit
 	end if
 	for j=1 to n1-1
  		a(j+1)=a(j+1)-1*p*a(j)
  		a(j+2)=a(j+2)-1*q*a(j)
 	next j
 	n1=n1-2
 	if n1>1 then
  		goto hleap
 	end if
 	goto  endit
  	label hleap
 	if n1>3 then
  	goto ileap
 	end if
 	m1=m1+1
 	e(m1)=1
 	p=a(2)/a(1)
 	q=a(3)/a(1)
 	goto  bleap
  	label ileap
next m1

label onetwo
if n=2 then
	goto  lleap
end if
print -1*a(2)/a(1)
goto  endit

label lleap
a(3)=a(2)*a(2)-4*a(1)*a(3)
s=(-1*a(2))/2/a(1)
t=(abs(a(3))^.5)/2/a(1)
m1=4
e(4)=4
if sgn(a(3))<0 then
	goto  jleap
end if
goto  kleap

label endit
end