UNPKG

5.32 kBtext/coffeescriptView Raw
1# find the roots of a polynomial numerically
2
3
4
5NROOTS_YMAX = 101
6NROOTS_DELTA = 1.0e-6
7NROOTS_EPSILON = 1.0e-9
8NROOTS_ABS = (z) ->
9 return Math.sqrt((z).r * (z).r + (z).i * (z).i)
10
11# random between -2 and 2
12theRandom = 0.0
13
14NROOTS_RANDOM = ->
15 #theRandom += 0.2
16 #return theRandom
17 return 4.0 * Math.random() - 2.0
18
19class numericRootOfPolynomial
20 r: 0.0
21 i: 0.0
22
23nroots_a = new numericRootOfPolynomial()
24nroots_b = new numericRootOfPolynomial()
25nroots_x = new numericRootOfPolynomial()
26nroots_y = new numericRootOfPolynomial()
27nroots_fa = new numericRootOfPolynomial()
28nroots_fb = new numericRootOfPolynomial()
29nroots_dx = new numericRootOfPolynomial()
30nroots_df = new numericRootOfPolynomial()
31nroots_c = []
32for initNRoots in [0...NROOTS_YMAX]
33 nroots_c[initNRoots] = new numericRootOfPolynomial()
34
35Eval_nroots = ->
36 h = 0
37 i = 0
38 k = 0
39 n = 0
40
41 push(cadr(p1))
42 Eval()
43
44 push(caddr(p1))
45 Eval()
46 p2 = pop()
47 if (p2 == symbol(NIL))
48 guess()
49 else
50 push(p2)
51
52 p2 = pop()
53 p1 = pop()
54
55 if (!ispoly(p1, p2))
56 stop("nroots: polynomial?")
57
58 # mark the stack
59
60 h = tos
61
62 # get the coefficients
63
64 push(p1)
65 push(p2)
66 n = coeff()
67 if (n > NROOTS_YMAX)
68 stop("nroots: degree?")
69
70 # convert the coefficients to real and imaginary doubles
71
72 for i in [0...n]
73 push(stack[h + i])
74 real()
75 yyfloat()
76 Eval()
77 p1 = pop()
78 push(stack[h + i])
79 imag()
80 yyfloat()
81 Eval()
82 p2 = pop()
83 if (!isdouble(p1) || !isdouble(p2))
84 stop("nroots: coefficients?")
85 nroots_c[i].r = p1.d
86 nroots_c[i].i = p2.d
87
88 # pop the coefficients
89
90 moveTos h
91
92 # n is the number of coefficients, n = deg(p) + 1
93
94 monic(n)
95
96 for k in [n...1] by -1
97 findroot(k)
98 if (Math.abs(nroots_a.r) < NROOTS_DELTA)
99 nroots_a.r = 0.0
100 if (Math.abs(nroots_a.i) < NROOTS_DELTA)
101 nroots_a.i = 0.0
102 push_double(nroots_a.r)
103 push_double(nroots_a.i)
104 push(imaginaryunit)
105 multiply()
106 add()
107 NROOTS_divpoly(k)
108
109 # now make n equal to the number of roots
110
111 n = tos - h
112
113 if (n > 1)
114 sort_stack(n)
115 p1 = alloc_tensor(n)
116 p1.tensor.ndim = 1
117 p1.tensor.dim[0] = n
118 for i in [0...n]
119 p1.tensor.elem[i] = stack[h + i]
120 moveTos h
121 push(p1)
122
123# divide the polynomial by its leading coefficient
124
125monic = (n) ->
126 k = 0
127 t = 0.0
128 nroots_y.r = nroots_c[n - 1].r
129 nroots_y.i = nroots_c[n - 1].i
130 t = nroots_y.r * nroots_y.r + nroots_y.i * nroots_y.i
131 for k in [0...(n - 1)]
132 nroots_c[k].r = (nroots_c[k].r * nroots_y.r + nroots_c[k].i * nroots_y.i) / t
133 nroots_c[k].i = (nroots_c[k].i * nroots_y.r - nroots_c[k].r * nroots_y.i) / t
134 nroots_c[n - 1].r = 1.0
135 nroots_c[n - 1].i = 0.0
136
137# uses the secant method
138
139findroot = (n) ->
140 j = 0
141 k = 0
142 t = 0.0
143
144 if (NROOTS_ABS(nroots_c[0]) < NROOTS_DELTA)
145 nroots_a.r = 0.0
146 nroots_a.i = 0.0
147 return
148
149 for j in [0...100]
150
151 nroots_a.r = NROOTS_RANDOM()
152 nroots_a.i = NROOTS_RANDOM()
153
154 compute_fa(n)
155
156 nroots_b.r = nroots_a.r
157 nroots_b.i = nroots_a.i
158
159 nroots_fb.r = nroots_fa.r
160 nroots_fb.i = nroots_fa.i
161
162 nroots_a.r = NROOTS_RANDOM()
163 nroots_a.i = NROOTS_RANDOM()
164
165 for k in [0...1000]
166
167 compute_fa(n)
168
169 nrabs = NROOTS_ABS(nroots_fa)
170 if DEBUG then console.log "nrabs: " + nrabs
171 if ( nrabs < NROOTS_EPSILON)
172 return
173
174 if (NROOTS_ABS(nroots_fa) < NROOTS_ABS(nroots_fb))
175 nroots_x.r = nroots_a.r
176 nroots_x.i = nroots_a.i
177
178 nroots_a.r = nroots_b.r
179 nroots_a.i = nroots_b.i
180
181 nroots_b.r = nroots_x.r
182 nroots_b.i = nroots_x.i
183
184 nroots_x.r = nroots_fa.r
185 nroots_x.i = nroots_fa.i
186
187 nroots_fa.r = nroots_fb.r
188 nroots_fa.i = nroots_fb.i
189
190 nroots_fb.r = nroots_x.r
191 nroots_fb.i = nroots_x.i
192
193 # dx = nroots_b - nroots_a
194
195 nroots_dx.r = nroots_b.r - nroots_a.r
196 nroots_dx.i = nroots_b.i - nroots_a.i
197
198 # df = fb - fa
199
200 nroots_df.r = nroots_fb.r - nroots_fa.r
201 nroots_df.i = nroots_fb.i - nroots_fa.i
202
203 # y = dx / df
204
205 t = nroots_df.r * nroots_df.r + nroots_df.i * nroots_df.i
206
207 if (t == 0.0)
208 break
209
210 nroots_y.r = (nroots_dx.r * nroots_df.r + nroots_dx.i * nroots_df.i) / t
211 nroots_y.i = (nroots_dx.i * nroots_df.r - nroots_dx.r * nroots_df.i) / t
212
213 # a = b - y * fb
214
215 nroots_a.r = nroots_b.r - (nroots_y.r * nroots_fb.r - nroots_y.i * nroots_fb.i)
216 nroots_a.i = nroots_b.i - (nroots_y.r * nroots_fb.i + nroots_y.i * nroots_fb.r)
217
218 stop("nroots: convergence error")
219
220compute_fa = (n) ->
221 k = 0
222 t = 0.0
223
224 # x = a
225
226 nroots_x.r = nroots_a.r
227 nroots_x.i = nroots_a.i
228
229 # fa = c0 + c1 * x
230
231 nroots_fa.r = nroots_c[0].r + nroots_c[1].r * nroots_x.r - nroots_c[1].i * nroots_x.i
232 nroots_fa.i = nroots_c[0].i + nroots_c[1].r * nroots_x.i + nroots_c[1].i * nroots_x.r
233
234 for k in [2...n]
235
236 # x = a * x
237
238 t = nroots_a.r * nroots_x.r - nroots_a.i * nroots_x.i
239 nroots_x.i = nroots_a.r * nroots_x.i + nroots_a.i * nroots_x.r
240 nroots_x.r = t
241
242 # fa += c[k] * x
243
244 nroots_fa.r += nroots_c[k].r * nroots_x.r - nroots_c[k].i * nroots_x.i
245 nroots_fa.i += nroots_c[k].r * nroots_x.i + nroots_c[k].i * nroots_x.r
246
247# divide the polynomial by x - a
248
249NROOTS_divpoly = (n) ->
250 k = 0
251 for k in [(n - 1)...0]
252 nroots_c[k - 1].r += nroots_c[k].r * nroots_a.r - nroots_c[k].i * nroots_a.i
253 nroots_c[k - 1].i += nroots_c[k].i * nroots_a.r + nroots_c[k].r * nroots_a.i
254
255 if (NROOTS_ABS(nroots_c[0]) > NROOTS_DELTA)
256 stop("nroots: residual error")
257
258 for k in [0...(n - 1)]
259 nroots_c[k].r = nroots_c[k + 1].r
260 nroots_c[k].i = nroots_c[k + 1].i
261
262