1 |
|
2 |
|
3 |
|
4 |
|
5 | NROOTS_YMAX = 101
|
6 | NROOTS_DELTA = 1.0e-6
|
7 | NROOTS_EPSILON = 1.0e-9
|
8 | NROOTS_ABS = (z) ->
|
9 | return Math.sqrt((z).r * (z).r + (z).i * (z).i)
|
10 |
|
11 |
|
12 | theRandom = 0.0
|
13 |
|
14 | NROOTS_RANDOM = ->
|
15 |
|
16 |
|
17 | return 4.0 * Math.random() - 2.0
|
18 |
|
19 | class numericRootOfPolynomial
|
20 | r: 0.0
|
21 | i: 0.0
|
22 |
|
23 | nroots_a = new numericRootOfPolynomial()
|
24 | nroots_b = new numericRootOfPolynomial()
|
25 | nroots_x = new numericRootOfPolynomial()
|
26 | nroots_y = new numericRootOfPolynomial()
|
27 | nroots_fa = new numericRootOfPolynomial()
|
28 | nroots_fb = new numericRootOfPolynomial()
|
29 | nroots_dx = new numericRootOfPolynomial()
|
30 | nroots_df = new numericRootOfPolynomial()
|
31 | nroots_c = []
|
32 | for initNRoots in [0...NROOTS_YMAX]
|
33 | nroots_c[initNRoots] = new numericRootOfPolynomial()
|
34 |
|
35 | Eval_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 |
|
59 |
|
60 | h = tos
|
61 |
|
62 |
|
63 |
|
64 | push(p1)
|
65 | push(p2)
|
66 | n = coeff()
|
67 | if (n > NROOTS_YMAX)
|
68 | stop("nroots: degree?")
|
69 |
|
70 |
|
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 |
|
89 |
|
90 | moveTos h
|
91 |
|
92 |
|
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 |
|
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 |
|
124 |
|
125 | monic = (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 |
|
138 |
|
139 | findroot = (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 |
|
194 |
|
195 | nroots_dx.r = nroots_b.r - nroots_a.r
|
196 | nroots_dx.i = nroots_b.i - nroots_a.i
|
197 |
|
198 |
|
199 |
|
200 | nroots_df.r = nroots_fb.r - nroots_fa.r
|
201 | nroots_df.i = nroots_fb.i - nroots_fa.i
|
202 |
|
203 |
|
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 |
|
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 |
|
220 | compute_fa = (n) ->
|
221 | k = 0
|
222 | t = 0.0
|
223 |
|
224 |
|
225 |
|
226 | nroots_x.r = nroots_a.r
|
227 | nroots_x.i = nroots_a.i
|
228 |
|
229 |
|
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 |
|
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 |
|
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 |
|
248 |
|
249 | NROOTS_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 |
|