However, this time they seem to have survived the onslaught.
Here is an assemblage which returns the N roots of a polynomium degree N.
The Eigenvalue sub can be used for, well, -- er-- Eigenvalues, if required.
By extracting the roots of a polynomium via Eigenvalues, the quality of the Eigenvalue can then easily be checked by testing the Zeros of the Polynomium.
Have fun!
Code: Select all
type complex
as double re,im
end type
declare sub balance(A()as double,d() as double)
declare sub GET_HESSENBERG(A() as double)
declare sub ITERATE_BY_QR(A() as double,wr() as double,wi() as double)
declare sub GET_EIGENVALUES(A() as double,Eigen() as complex)
declare sub makecompanion(polynomium() as double,mat() as double)
declare Function Cpoly(coff() As complex,number As complex)As complex
declare sub equalize(m1() as double,m2() as double)
declare Function decround ( a As Double, b As Integer) As Double
declare operator +(n1 as complex,n2 as complex) as complex
declare operator -(n1 as complex,n2 as complex) as complex
declare operator *(n1 as complex,n2 as complex) as complex
declare operator /(n1 as complex,n2 as complex) as complex
sub printout(bu()as double) 'only used for de-bugging.
for x as integer=1 to ubound(bu)
for y as integer=1 to ubound(bu)
print bu(x,y);
next y
print
next x
end sub
' USAGE
dim n as integer
input "Enter degree of polynomium";n
print
dim as double poly(n)
for z as integer=n to 0 step -1
if z<>0 then
print "Enter coefficient of x^";str(z);" ";
input poly(z)
else
print "Enter constant term ";
input poly(z)
end if
next z
dim as double val1=poly(n)
for z as integer=n to 0 step -1
poly(z)=poly(z)/val1
next z
print
redim companion(0,0) as double
makecompanion(poly(),companion())
''printout(companion())
redim as complex roots(0)
GET_EIGENVALUES(companion(),roots())
print " ROOTS"
print
print "REAL PART IMAGINARY PART"
print
for x as integer=1 to n
print roots(x).re,roots(x).im
next x
Print "Press a key"
sleep
dim as complex polyC(n),result
for z as integer=0 to n
polyc(z).re=poly(z)
polyc(z).im=0
next z
print "Check the roots by evaluating polynomium at the above roots (6 decimal places)"
print
print " REAL RETURN IMAGINARY RETURN"
print
for z as integer=1 to n
result=cpoly(polyC(),roots(z))
print "Root";z,decround(result.re,6),decround(result.im,6)
next z
print
print "DONE"
sleep
'START PROCEDURES
#macro EXC(m)
d(m) = 1# * j
IF j <> m THEN
FOR ii as integer = 1 TO k
f = A(ii, j)
A(ii, j) = A(ii, m)
A(ii, m) = f
NEXT ii
FOR ii as integer = l TO n
f = A(j, ii)
A(j, ii) = A(m, ii)
A(m, ii) = f
NEXT ii
END IF
#endmacro
sub balance(A() as double,d() as double)
dim n as integer = ubound(A)
dim as double b = 2
dim as integer low,ihi,noconv
dim as integer i,j,k,l,m
dim as double b2,c,f,g,r,s
dim as double zero=0
b2 = b * b
l = 1
k = n
lab1: FOR j = k TO 1 STEP -1
r = ZERO
FOR i = 1 TO j - 1
r = r + ABS(A(j, i))
NEXT i
FOR i = j + 1 TO k
r = r + ABS(A(j, i))
NEXT i
IF (r = 0) THEN
m = k
Exc(k)
k = k - 1
GOTO lab1
END IF
NEXT j
lab2: FOR j = l TO k
c = ZERO
FOR i = l TO j - 1
c = c + ABS(A(i, j))
NEXT i
FOR i = j + 1 TO k
c = c + ABS(A(i, j))
NEXT i
IF (c = 0) THEN
m = l: Exc(l)
l = l + 1
GOTO lab2
END IF
NEXT j
low = l
ihi = k
FOR i = 1 TO k
d(i) = 1#
NEXT i
lab3: noconv = 0
FOR i = l TO k
c = ZERO
r = c
FOR j = l TO i - 1
c = c + ABS(A(j, i))
r = r + ABS(A(i, j))
NEXT j
FOR j = i + 1 TO k
c = c + ABS(A(j, i))
r = r + ABS(A(i, j))
NEXT j
g = r / b
f = 1#
s = c + r
lab4: IF c < g THEN
f = f * b
c = c * b2
GOTO lab4
END IF
g = r * b
lab5: IF c >= g THEN
f = f / b
c = c / b2
GOTO lab5
END IF
IF (c + r) / f < .95 * s THEN
g = 1# / f
d(i) = d (i) * f
noconv = 1
FOR j = l TO n
A(i, j) = A(i, j) * g
NEXT j
FOR j = 1 TO k
A(j, i) = A(j, i) * f
NEXT j
END IF
NEXT i
IF noconv = 1 THEN GOTO lab3
end sub
sub GET_HESSENBERG(A() as double)
dim n as integer=ubound(a)
dim as integer m,j,i
dim as double y,x,ZERO=0
IF n > 2 THEN
FOR m = 2 TO n - 1
x = ZERO
i = m
FOR j = m TO n
IF ABS(A(j, m - 1)) > ABS(x) THEN
x = A(j, m - 1)
i = j
END IF 'IF Abs
NEXT j 'FOR j= m TO n
IF i <> m THEN
FOR j = m - 1 TO n
y = A(i, j)
A(i, j) = A(m, j)
A(m, j) = y
NEXT j
FOR j = 1 TO n
y = A(j, i)
A(j, i) = A(j, m)
A(j, m) = y
NEXT j
END IF 'IF i <> m
IF x <> ZERO THEN
FOR i = m + 1 TO n
y = A(i, m - 1)
IF y <> ZERO THEN
y = y / x
A(i, m - 1) = y
FOR j = m TO n
A(i, j) = A(i, j) - y * A(m, j)
NEXT j
FOR j = 1 TO n
A(j, m) = A(j, m) + y * A(j, i)
NEXT j
END IF 'IF y
NEXT i 'FOR i
END IF 'IF x
NEXT m 'FOR m
END IF 'if n>2
end sub
sub ITERATE_BY_QR(A() as double,wr() as double,wi() as double)
dim as double y,r,u,v,mmin
dim as string res
dim sign as double
dim as double ZERO=0
dim n as integer=ubound(A)
dim as double anorm,itsmx,t,its,s,x,p,q,nn,z,w
dim as integer i,j,l,m,k
dim as double aa,bb
#macro _SIGN(aa,bb)
IF bb < 0 THEN
Sign = -ABS(aa)
ELSE
Sign = ABS(aa)
END IF
#endmacro
itsmx = 50
anorm = ABS(A(1, 1))
FOR i = 2 TO n
FOR j = i - 1 TO n
anorm = anorm + ABS(A(i, j))
NEXT j
NEXT i
nn = n
t = ZERO
label4: its = 0
label2: FOR l = nn TO 2 STEP -1
s = ABS(A(l - 1, l - 1)) + ABS(A(l, l))
IF (s = 0!) THEN s = anorm
IF ((ABS(A(l, l - 1)) + s) = s) THEN GOTO label3
NEXT l
l = 1
label3: x = A(nn, nn)
IF (l = nn) THEN
wr(nn) = x + t
wi(nn) = ZERO
nn = nn - 1
ELSE
y = A(nn - 1, nn - 1)
w = A(nn, nn - 1) * A(nn - 1, nn)
IF (l = nn - 1) THEN
p = .5# * (y - x)
q = p * p + w
z = SQR(ABS(q))
x = x + t
IF q >= ZERO THEN
aa = z: bb = p: _SIGN(aa,bb)
z = p + Sign
wr(nn) = x + z
wr(nn - 1) = wr(nn)
IF z <> ZERO THEN wr(nn) = x - w / z
wi(nn) = ZERO
wi(nn - 1) = ZERO
ELSE
wr(nn) = x + p
wr(nn - 1) = wr(nn)
wi(nn) = z
wi(nn - 1) = -z
END IF
nn = nn - 2
ELSE
IF its = itsmx THEN
PRINT " Pause"
PRINT " Too many iterations!"
INPUT res
END IF
IF (its = 10) OR (its = 20) THEN
t = t + x
FOR i = 1 TO nn
A(i, i) = A(i, i) - x
NEXT i
s = ABS(A(nn, nn - 1)) + ABS(A(nn - 1, nn - 2))
x = .75# * s
y = x
w = -.4375# * s * s
END IF
its = its + 1
FOR m = nn - 2 TO 1 STEP -1
z = A(m, m)
r = x - z
s = y - z
p = (r * s - w) / A(m + 1, m) + A(m, m + 1)
q = A(m + 1, m + 1) - z - r - s
r = A(m + 2, m + 1)
s = ABS(p) + ABS(q) + ABS(r)
p = p / s
q = q / s
r = r / s
IF m = 1 THEN GOTO label1
u = ABS(A(m, m - 1)) * (ABS(q) + ABS(r))
v = ABS(p) * (ABS(A(m - 1, m - 1)) + ABS(z) + ABS(A(m + 1, m + 1)))
IF u + v = v THEN GOTO label1
NEXT m
label1: FOR i = m + 2 TO nn
A(i, i - 2) = ZERO
IF i <> (m + 2) THEN A(i, i - 3) = ZERO
NEXT i
FOR k = m TO nn - 1
IF k <> m THEN
p = A(k, k - 1)
q = A(k + 1, k - 1)
r = ZERO
IF k <> (nn - 1) THEN r = A(k + 2, k - 1)
x = ABS(p) + ABS(q) + ABS(r)
IF x <> ZERO THEN
p = p / x
q = q / x
r = r / x
END IF
END IF
aa = SQR(p * p + q * q + r * r): bb = p: _SIGN(aa,bb)'GOSUB 2900
s = Sign
IF s <> ZERO THEN
IF k = m THEN
IF l <> m THEN A(k, k - 1) = -A(k, k - 1)
ELSE
A(k, k - 1) = -s * x
END IF
p = p + s
x = p / s
y = q / s
z = r / s
q = q / p
r = r / p
FOR j = k TO nn
p = A(k, j) + q * A(k + 1, j)
IF k <> (nn - 1) THEN
p = p + r * A(k + 2, j)
A(k + 2, j) = A(k + 2, j) - p * z
END IF
A(k + 1, j) = A(k + 1, j) - p * y
A(k, j) = A(k, j) - p * x
NEXT j
IF nn < k + 3 THEN
mmin = nn
ELSE
mmin = k + 3
END IF
FOR i = 1 TO mmin
p = x * A(i, k) + y * A(i, k + 1)
IF k <> (nn - 1) THEN
p = p + z * A(i, k + 2)
A(i, k + 2) = A(i, k + 2) - p * r
END IF
A(i, k + 1) = A(i, k + 1) - p * q
A(i, k) = A(i, k) - p
NEXT i
END IF
NEXT k
GOTO label2:
END IF
END IF
IF nn >= 1 THEN GOTO label4:
end sub
sub GET_EIGENVALUES(MAT() as double,Eig() as complex)
dim n as integer=ubound(MAT)
redim Eig(n)
dim as double A(n,n)
equalize(A(),MAT())
dim as double wr(n),wi(n),d(n)
balance(a(),d())
GET_HESSENBERG(A())
ITERATE_BY_QR(A(),wr(),wi())
FOR i as integer = 1 TO n
Eig(i).re= wr(i):Eig(i).im= wi(i)
NEXT i
end sub
sub makecompanion(polynomium() as double,mat() as double)
dim as integer di=ubound(polynomium)
dim as double one,null=0
one=1:
redim mat (1 to di,1 to di) as double
for a as integer=1 to di
mat(1,di)=null-polynomium(0)
for b as integer=1 to di
if a=b+1 then mat(a,b)=one
if b>1 then
mat(b,di)=null-polynomium(b-1)
end if
next b
next a
end sub
Function Cpoly(coff() As complex,number As complex)As complex
Dim count As Integer 'evaluates the polynomial
Dim pol As complex
Dim deg As Integer=Ubound(coff)
dim temp as complex
For count = 1 To DEG + 1
temp.re=1:temp.im=0
for k as integer=1 to count-1
temp=temp*number
next k
pol = pol + coff(count-1) * temp
Next count
Cpoly = pol
End Function
sub equalize(m1() as double,m2() as double)
for x as integer=1 to ubound(m1)
for y as integer=1 to ubound(m1)
m1(x,y)=m2(x,y)
next y
next x
end sub
Function decround ( a As Double, b As Integer) 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
operator +(n1 as complex,n2 as complex) as complex
dim as complex n
n.re=n1.re+n2.re
n.im=n1.im+n2.im
return n
end operator
operator -(n1 as complex,n2 as complex) as complex
dim as complex n
n.re=n1.re-n2.re
n.im=n1.im-n2.im
return n
end operator
operator *(n1 as complex,n2 as complex) as complex
dim as complex n
n.re=n1.re*n2.re - n1.im*n2.im
n.im=n1.im*n2.re + n1.re*n2.im
return n
end operator
operator /(n1 as complex,n2 as complex) as complex
dim as complex n
dim as double d
d=n2.re^2+n2.im^2
n.re=(n1.re*n2.re+n1.im*n2.im)/d
n.im=(n1.im*n2.re - n1.re*n2.im)/d
return n
end operator