-
Notifications
You must be signed in to change notification settings - Fork 4
Expand file tree
/
Copy path120.bas
More file actions
220 lines (186 loc) · 3.46 KB
/
Copy path120.bas
File metadata and controls
220 lines (186 loc) · 3.46 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
'''This program will display the real and imaginary roots for a polynomial of the form f(x)=A0 + A1*x + A2*x^2 +...
''spagetti code in http://www3.sympatico.ca/ltoms/html/bairstow_s_method.html modified by:
''Adolfo Leon Sepulveda
''Aug-01-2004
''Version 1.0
''Test case:
''real and imaginary roots of polynomials
''using bairstow''s method
''form is y = A0 + A1*x + A2*x^2 + A3*x^3 +....
''degree of polynomial: 4
''coefficent of A(0)= 12
''coefficent of A(1)= -19
''coefficent of A(2)= 12
''coefficent of A(3)= -6
''coefficent of A(4)= 1
''Roots:
''0.5 + 1.6583124 i
''0.5 - 1.6583124 i
''4.0
''1.0
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)
Const err=.0001
Const cero=1e-19
Const iter=100
InputData a,n
print
print "Roots:"
Bairstow a,n,err,iter
End
sub InputData( ByRef a(), ByRef n)
local i,xa
input "degree of polynomial: ";n
print
for i=0 to n
print "coefficent of A(";i;")= ";
input " ";xa
rem stored in array a in reverse order
a(n-i+1)=xa
next i
end
Sub Bairstow(a,n,err,iter)
dim b(22)
dim e(22)
local p, q, p1, q1, r0, r1, v0, v1, s0, s1, d0, d1, d2, s, t
local m1, n1, k, i, h, j
local const cero = 1e-19
rem branch for special treatment of 1st and 2nd degree equations
if n<=2 then
OrderLessThan2 a, n
Exit Sub
endif
a(n+2)=0
n1=2*int((n+1)/2)
m1 = 1
While true
If m1 >= n1 / 2 Then Exit loop
p=1
q=1
For k=1 to iter
while true
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) - p*b(i)
b(i+2) = b(i+2) - q*b(i)
next i
next j
r0 = b(n1+1)
r1 = b(n1)
s0 = b(n1-1)
s1 = b(n1-2)
v0 = -q*s1
v1 = s0 - s1*p
d0 = v1*s0 - v0*s1
if abs(d0) >= cero then Exit loop
p=p+5
q=q+5
wend
d1 = s0*r1 - s1*r0
d2 = r0*v1 - v0*r1
p1=d1/d0
q1=d2/d0
p=p+p1
q=q+q1
if !(Abs(r0)>=err or abs(r1)>=err) then
e(m1)=1
Exit For
endif
if !(abs(p1)>=err or abs(q1)>=err) then
e(m1)=2
Exit For
endif
if p=0 then
if !(q=0) Then
if !(abs(q1/q)>=err) Then
e(m1)=3
Exit For
endif
endif
Else
if !(abs(p1/p)>=err) Then
if !(q=0) Then
if !(abs(q1/q)>=err) Then
e(m1)=3
Exit For
endif
endif
endif
endif
Next k
If k > iter Then
e(m1) = 4 ''Is Order 2
endif
While true
s = -p/2
t = s^2 - q
if !(t<0) then
t=Sqr(t)
PrintReal s,t
Else
t=Sqr(-1*t)
PrintImag s,t
endif
If e(m1) = 4 then Exit Sub
for j=1 to n1-1
a(j+1) = a(j+1) - p*a(j)
a(j+2) = a(j+2) - q*a(j)
next j
n1 = n1 - 2
if !(n1>1) Then Exit Sub
if n1>3 then Exit loop
m1 = m1 + 1
e(m1)=1
p=a(2)/a(1)
q=a(3)/a(1)
Wend
Wend ''m1
end
Sub OrderLessThan2(a(),n)
If n = 2 Then
Order2 a
Else
print -a(2)/a(1)
Endif
End
sub order2(a())
local s,t
a(3) = a(2)*a(2) - 4*a(1)*a(3)
s = -a(2) / (2*a(1))
t=Sqr(abs(a(3))) / (2*a(1))
if sgn(a(3))<0 then
PrintImag s,t
Else
PrintReal s,t
endif
End
Sub PrintImag(s,t)
Print
If Abs(s) > cero
print s;" + ";t;" i"
print s;" - ";t;" i"
Else
Print 0;" + ";t;" i"
print 0;" - ";t;" i"
endif
End
Sub PrintReal(s,t)
print
If Abs(s+t) > cero
print s+t
endif
If Abs(s-t) > cero
print s-t
endif
End
'