If the coefficients are real a little plot is shown.
Uses Newtons method via iterations, so it might come up short now and then.
If the degree of the polynomial is d, it should have d roots.
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
#macro cleanup(a,b)
redim as typeof(a) b(0)
scope
dim as long flag
for n1 as long=lbound(a) to ubound(a)
flag=0
for n2 as long=n1+1 to ubound(a)
if a(n1)=a(n2) then flag=1
next n2
if flag=0 then
redim preserve b(1 to ubound(b)+1)
b(ubound(b))=a(n1)
end if
next n1
end scope
#endmacro
Function decround ( a As Double, b As long) As Double
Dim y As Double
Dim i As Double
Dim r As Double
y = (Abs(a) - Int(Abs(a))) * (10 ^ b)
i = Int(y)
y = y - i
If y >= .5 Then i = i + 1
i = i / (10 ^ b)
r = Int(Abs(a)) + i
If a < 0 Then r = -r
decround = r
End Function
#define Cabs(z) sqr(z.re*z.re+z.im*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 lbound(coff) 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 Newton(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) ''x1=x0-(f(x0)/f1(x0))
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,byval range As Single,roots() As complex)
print "Message from iterator"
#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 complex start,root,er
Dim As complex grad(Ubound(pol))
Redim As complex lastroot()
differentiate(pol(),grad())
Dim As Integer count,flag,count2
Dim As Single stepsize=.5
start:
For x As Single=-range To range Step stepsize
For y As Single=-range To range Step stepsize
start=Type<complex>(x,y)
root=newton(pol(),grad(),start,iterate)
er=poly(pol(),root)
If cabs(er)> 1e-6 Then Goto jump
flag=1
count=count+1
For z2 As Integer=1 To count-1
If root=lastroot(z2) Then flag=0:Goto nxt
If cabs((root-lastroot(z2)))<1e-10 Then flag=0:Goto nxt
Next z2
If flag Then
count2=count2+1
Redim Preserve roots(1 to count2)
roots(count2)=root
End If
nxt:
Redim Preserve lastroot(1 to count)
lastroot(count)=root
jump:
Next y
Next x
If range>10 Then 'arbitary get out
Beep
Print "ERROR -- Range "; range;" Iteratations ";iterate
sort(roots)
Exit Sub
End If
If range>5 Then
stepsize=1
iterate=iterate+1000
End If
If Ubound(roots)< Ubound(pol) Then
range=range+1:iterate=iterate+10:Goto start
end if
sort(roots)
Print "Range ";range;" Iteratations ";iterate
print
for n as long=lbound(roots) to ubound(roots)
roots(n).re=decround(roots(n).re,5)
roots(n).im=decround(roots(n).im,5)
next
redim as complex outarray()
cleanup(roots,outarray)
redim roots(lbound(outarray) to ubound(outarray))
for n as long=lbound(outarray) to ubound(outarray)
roots(n)=outarray(n)
next
End Sub
Sub printroots(pol() As complex,roots() As complex)
Dim As complex er
Dim As Double er2
Print "ROOTS";
Print Tab(5);"Real part";Tab(30);"Imag part";Tab(57);"Error"
For z As Double=lbound(roots) 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
End Sub
redim shared as complex pL(),roots()
dim shared as string fn
do
print
dim as string s
dim as long p,ctr
dim as double x,y
redim roots(0)
input "Enter highest power of polynomial ",p
redim pL(p)
print "Enter the coefficients: real part first then comma then imaginary part"
print "If the imaginary part = 0 then you can press return after entering the real part"
print
dim as string msg1="Enter the complex coefficients of x^"
dim as string msg2="Enter the complex constant term "
ctr=p
do
if p>0 then
if p<>1 then print msg1;str(p); else print rtrim(msg1,"^") + " ";
else
print msg2;
end if
input" ", x,y
pL(ctr)=type(x,y)
s+="("+str(x)+","+str(y)+")*x^"+str(p)+" + "
p-=1
ctr-=1
loop until ctr=-1
print "The polynomium in ordered pairs (real,imaginary)"
s= rtrim(s," + ")
s=rtrim(s,"*x^0")
print s
fn=s
print
get_roots(pL(),20,2,roots())
printroots(pL(),roots())
print
print "Another polynomium? y/n "
s=input(1)
print
if lcase(s)="n" then exit do
loop
print "done"
'set up to sketch
dim as double sx=1e10,bx=-1e10,sy=1e10,by=-1e10
for n as long=lbound(pL) to ubound(pL)
if pL(n).im<>0 then end 'only plot with real coefficients
next
'to sketch, get an x domain.
for n as long=lbound(roots) to ubound(roots)
if bx<roots(n).re then bx=roots(n).re
if sx>roots(n).re then sx=roots(n).re
next
if sx=bx then
sx-=1
bx+=1
end if
'get the function y range rounded to six dec places
for n as double=sx to bx step(bx-sx)/100
var value=poly(pL(),type(n,0))
if by<value.re then by=value.re
if sy>value.re then sy=value.re
next n
by=decround(by,5)
sy=decround(sy,5)
'' plot part ===
dim as integer xres,yres
dim as double minx,maxx,miny,maxy,PLOT_GRADE=5000
screen 19,32
screeninfo xres,yres
#macro _window(topleftX,topleftY,bottomrightX,bottomrightY)
minx=(topleftX)
maxx=(bottomrightX)
miny=(bottomrightY)
maxy=(topleftY)
#endmacro
#macro drawaxis(colour)
line(0,(yres-(miny/(miny-maxy))*yres))-(xres,(yres-(miny/(miny-maxy))*yres)),(colour) 'x axis
line(((minx/(minx-maxx))*xres),0)-(((minx/(minx-maxx))*xres),yres),(colour) 'y axis
draw string (0,yres/2),str(minx),(colour)
draw string (xres-8-8*(len(str(maxx))),yres/2),str(maxx),(colour)
draw string(xres/2,0),str(maxy),(colour)
draw string(xres/2,yres-16),str(miny),(colour)
#endmacro
#macro sketch(_function,colour)
locate 4
print fn
For x As Double=minx To maxx Step (maxx-minx)/PLOT_GRADE
dim as double x1=(xres)*(x-minx)/(maxx-minx)
dim as double y1=(yres)*((_function)-maxy)/(miny-maxy)
if x=minx then Pset(x1,y1),colour else line -(x1,y1),colour
Next x
#endmacro
function f(x as double) as double'the polynomial wrapped in a function
return poly(pl(),type(x,0)).re
end function
_window(sx,by,bx,sy) 'set cartesian window area
sketch(f(x),rgb(150,150,150))
drawaxis(rgb(100,255,100))
sleep