©Set Complex format to Rectangular
©12may02/dburkett@infinet.com
local m,roots,eps,maxm,i,j,jj,ad,x,b,c,root
©------------------------------------------------------------
© Define local function to find 1 root with Laguerre'e method
Define root(a,x)=Func
©({coeffs},guess) polynomial root
©NOTE: {a} is reversed!
local m,maxit,mr,mt,epss,iter,j,abx,abp,abm,err,frac,dx,x1,b,d,f,g,h,sq,gp,gm,g2
dim(a)-1→m © Find order of polynomial
2⁻12→epss © Set estimated fractional error
8→mr © Set number of fractional values to break limit cycles
10→mt © Set number of steps to break limit cycle
mt*mr→maxit © Maximum number of iterations
{.5,.25,.75,.13,.38,.62,.88,1.}→frac © Fractions used to break limit cycle
for iter,1,maxit © Loop to find a root
a[m+1]→b
abs(b)→err
0→d
0→f
abs(x)→abx
for j,m,1,⁻1 © Evaluate the polynomial and first two derivatives
x*f+d→f
x*d+b→d
x*b+a[j]→b
abs(b)+abx*err→err
endfor
epss*err→err © Estimate roundoff error in polynomial evaluation
if abs(b)≤err then © Return root if error condition met
return x
else © Execute Laguerre's method
d/b→g
g^2→g2
g2-2.*f/b→h
√((m-1)*(m*h-g2))→sq
g+sq→gp
g-sq→gm
abs(gp)→abp
abs(gm)→abm
if abp < abm: gm→gp
if max(abp,abm)>0 then
m/gp→dx
else
ℯ^(ln(1+abx)+iter*)→dx
endif
endif
x-dx→x1
if x=x1:return x © Converged to root
if mod(iter,mt)≠0 then © Every mt steps take a fractional step to break (rare) limit cycle
x1→x
else
x-dx*frac[iter/mt]→x
endif
endfor
return "polyroot too many its" © May have this error for complex roots
EndFunc
©-----------------------------------------------------------
© Begin polyroot()
dim(a)-1→m © Find order of polynomial
6 - 120
To Next Page
Loading...
Need help?
General
