Here is a rough and ready method to get the roots of a polynomial (real or complex coefficients).
You can tweak two parameters if necessary.
It is not 100% guaranteed for every polynomial, but then again neither is the eigenvalue method, which I assume gsl uses.
Code: Select all
Type complex
As Double re,im
End Type
Operator +(Byval n1 As complex,Byval n2 As complex) As complex
Return Type<complex>(n1.re+n2.re,n1.im+n2.im)
End Operator
Operator -(Byval n1 As complex,Byval n2 As complex) As complex
Return Type<complex>(n1.re-n2.re,n1.im-n2.im)
End Operator
Operator *(Byval n1 As complex,Byval n2 As complex) As complex
Return Type<complex>(n1.re*n2.re - n1.im*n2.im,n1.im*n2.re + n1.re*n2.im)
End Operator
Operator /(Byval n1 As complex,Byval n2 As complex) As complex
Dim As Double d = n2.re*n2.re+n2.im*n2.im
Return Type<complex>((n1.re*n2.re+n1.im*n2.im)/d,(n1.im*n2.re - n1.re*n2.im)/d)
End Operator
Operator *(Byval n1 As Double,Byval n2 As complex) As complex
Return Type<complex>(n1*n2.re,n1*n2.im)
End Operator
Operator =(Byval n1 As complex,Byval n2 As complex) As Integer
Return n1.re=n2.re And n1.im=n2.im
End Operator
#define Cabs(z) sqr(z.re*z.re+z.im*z.im)
#define conjugate(z) type<complex>(z.re,-z.im)
Function poly(coff() As complex,Byval z As complex) As complex
Dim As Integer deg=Ubound(coff)
Dim As complex p=coff(deg)
For i As Integer=deg-1 To 0 Step -1
p=p*z+coff(i)
Next i
Return p
End Function
Sub differentiate(pol() As complex,Dpol() As complex)
Redim Dpol(Ubound(pol)-1)
For count As Integer=0 To Ubound(Dpol)
Dpol(count)=(count+1)*pol(count+1)
Next count
End Sub
Dim Shared As Integer E
Function Newton1(f0() As complex,f1() As complex,Byval start As complex,Byval iterations As Integer) As complex
Dim As complex x1=start,x2
E=iterations
For z As Integer=0 To iterations
x2=x1-poly(f0(),x1)/poly(f1(),x1)
If x1=x2 Then E=z: Exit For
x1=x2
Next z
Return x2
End Function
Function Newton(f0() As complex,f1() As complex,f2() As complex,Byval start As complex,Byval iterations As Integer) As complex
Dim As complex x1=start,x2
E=iterations
Dim As complex numerator,denominator,d0,d1,d2
For z As Integer=0 To iterations
d0=poly(f0(),x1)
d1=poly(f1(),x1)
d2=poly(f2(),x1)
numerator=d0*d1
denominator=d1*d1-.5*d0*d2
x2=x1-numerator/denominator
If x1=x2 Then E=z: Exit For
x1=x2
Next z
Return x2
End Function
Sub get_roots(pol() As complex,Byval iterate As Integer=1,Byval range As Single=1,roots() As complex)
#macro sort(a)
Scope
Dim As Integer n=Ubound((a))
For p1 As Integer = 1 To n - 1
For p2 As Integer = p1 + 1 To n
If a(p1).re>=a(p2).re Then Swap a(p1),a(p2)
Next p2
Next p1
End Scope
#endmacro
Dim As Integer real=1
Dim As complex cj
For z As Integer=0 To Ubound(pol)
If pol(z).im<>0 Then real=0:Exit For
Next z
Dim As Single mag
For z As Integer=0 To Ubound(pol)
mag=mag+Abs(pol(z).re)
Next z
mag=mag/Ubound(pol)
Dim As complex start,root,er,centre
Dim As complex grad(Ubound(pol))
Redim As complex lastroot()
differentiate(pol(),grad())
Dim As complex grad2(Ubound(grad))
differentiate(grad(),grad2())
Dim As Integer count,flag,count2,deg=Ubound(pol)/10+1,ecount,rootflag=1
Dim As Double stepsize=1,small=1/(10^deg)
small=.00000001
start:
Randomize
'do
For x As Double=centre.re-range To centre.re+range Step stepsize
For y As Double=centre.im-range To centre.im+range Step stepsize
start=Type<complex>(x,y)
If rootflag=1 Then root=newton(pol(),grad(),grad2(),start,iterate)
If rootflag=2 Then root=newton1(pol(),grad(),start,iterate)
er=poly(pol(),root)
If cabs(er)> small Then Goto jump '1e-6
If Abs(root.im)<1e-12 Then root.im=0
flag=1
If e=iterate Then
ecount=ecount+1
centre=Type<complex>(10*Rnd*1-10*Rnd*1,10*Rnd*1-10*Rnd*1)
End If
count=count+1
For z2 As Integer=1 To Ubound(lastroot)
If cabs((root-lastroot(z2)))<1e-10 Then flag=0:Exit For
Next z2
If flag =1 Then
count2=count2+1
Redim Preserve roots(count2)
roots(count2)=root
Redim Preserve lastroot(count2)
lastroot(count2)=root
If Ubound(roots)= Ubound(pol) Then Goto oot
centre=Type<complex>(root.re,root.im)
End If
nxt:
jump:
Next y
Next x
If range>29 Then 'arbitary get out
Beep
Print "ERROR -- Range "; range;" Iteratations ";iterate
sort(roots)
Exit Sub
End If
If range>23 Then rootflag=2
If range>20 Then
centre=Type<complex>(10*Rnd*1-10*Rnd*1,10*Rnd*1-10*Rnd*10)
End If
If range>10 Then
If Ubound(roots) Then centre=Type<complex>(roots(Ubound(roots)-1).re,roots(Ubound(roots)-1).im)
End If
If range>5 Then
stepsize=.9
End If
If Ubound(roots)< Ubound(pol) Then
range=range+1:iterate=iterate+10:Goto start
End If
oot:
sort(roots)
' Print "Range ";range;" Iteratations ";iterate;" r ";rootflag
End Sub
Sub printroots(pol() As complex,roots() As complex)
Dim As complex er
Dim As Double er2
Print "ROOTS";
Print Tab(5);"Real";Tab(30);"Imag";Tab(57);"Error"
For z As Double=1 To Ubound(roots)
er=poly(pol(),roots(z)):er2=cabs(er)
Print Tab(0);z & ") ";Tab(5);roots(z).re;Tab(30);roots(z).im & "i";Tab(55); Using "###.###";er2;
Next z
print
Print string(60,"_")
print
End Sub
Redim As complex answer()
Dim As complex pol1(5)={ (-1),(0),(0),(0),(0),(1)} '-1 + x^5
get_roots(pol1(),,,answer())
printroots(pol1(),answer())
'complex quadratic polynomial
dim as complex pol2(...)={ (-1,1),(8.8,-7),(3,-3)}
get_roots(pol2(),,,answer())
printroots(pol2(),answer())
dim as complex pol3(...)={(0),(5),(0),(-20),(0),(16)} ' 5*x-20*x^3 +16*x^5
get_roots(pol3(),,,answer())
printroots(pol3(),answer())
Sleep
