rem roots of polynomials - bairstow's method 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=.00000000000000000001 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 a=val(a$) rem stored in array a in reverse order a(n-i+1)=a next i print print "Roots:" rem branch for special treatment of 1st and 2nd degree equations if n<=2 goto [onetwo] 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] p=p+5 q=q+5 next l [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] if abs(r1)>=e1 goto [aleap] e(m1)=1 goto [bleap] [aleap] if abs(p1)>=e1 goto [cleap] if abs(q1)>=e1 goto [cleap] e(m1)=2 goto [bleap] [cleap] if p=0 goto [dleap] if abs(p1/p)>=e1 goto [eleap] [dleap] if q=0 goto [eleap] if abs(q1/q)>=e1 goto [eleap] e(m1)=3 goto [bleap] [eleap] next k e(m1)=4 [bleap] s=(-1*p)/2 t=s*s-1*q if t<0 goto [fleap] t=t^.5 [kleap] print print s+t print s-1*t goto [gleap] [fleap] t=(-1*t)^.5 [jleap] print print s;" + I * ";t print s;" - I * ";-1*t [gleap] if e(m1)=4 goto [endit] 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 goto [hleap] goto [endit] [hleap] if n1>3 goto [ileap] m1=m1+1 e(m1)=1 p=a(2)/a(1) q=a(3)/a(1) goto [bleap] [ileap] next m1 [onetwo] if n=2 goto [lleap] print -1*a(2)/a(1) goto [endit] [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 goto [jleap] goto [kleap] function sgn(num) if num>0 then sgn=1 if num<0 then sgn=-1 if num=0 then sgn=0 end function [endit] end