1 /**********
2 Copyright 1990 Regents of the University of California. All rights reserved.
3 **********/
4
5 #include "spice.h"
6 #include <stdio.h>
7 #include <math.h>
8 #include "cktdefs.h"
9 #include "pzdefs.h"
10 #include "util.h"
11 #include "sperror.h"
12 #include "niext.h"
13
14 #define DEBUG(N) if (0)
15 /*if (Debug >= (unsigned) (N))
16 static unsigned int Debug = 0;
17 */
18
19 extern CKTpzTrapped;
20 double NIpzK;
21 int NIpzK_mag;
22
NIpzSym(set,new)23 NIpzSym(set, new)
24 PZtrial *set[3], *new;
25 {
26 #ifndef notdef
27 return NIpzSym2(set, new);
28 #else
29 double a, b, c, x0, x1;
30 double dx0, dx1;
31 int a_mag, b_mag, c_mag;
32
33 dx0 = set[1]->s.real - set[0]->s.real;
34 dx1 = set[2]->s.real - set[1]->s.real;
35
36 zaddeq(&a, &a_mag, set[1]->f_def.real, set[1]->mag_def,
37 -set[0]->f_def.real, set[0]->mag_def);
38 a /= dx0;
39 zaddeq(&b, &b_mag, set[2]->f_def.real, set[2]->mag_def,
40 -set[1]->f_def.real, set[1]->mag_def);
41 b /= dx1;
42 zaddeq(&c, &c_mag, b, b_mag, -a, a_mag);
43
44 /* XXX What if c == 0.0 ? */
45
46 x0 = (set[0]->s.real + set[1]->s.real) / 2.0;
47 x1 = (set[1]->s.real + set[2]->s.real) / 2.0;
48
49 c /= (x1 - x0);
50
51 new->s.real = - a / c;
52 c_mag -= a_mag;
53
54 new->s.imag = 0.0;
55
56 while (c_mag > 0) {
57 new->s.real /= 2.0;
58 c_mag -= 1;
59 }
60 while (c_mag < 0) {
61 new->s.real *= 2.0;
62 c_mag += 1;
63 }
64 new->s.real += set[0]->s.real;
65 #endif
66 }
67
NIpzComplex(set,new)68 NIpzComplex(set, new)
69 PZtrial *set[3], *new;
70 {
71 return NIpzSym2(set, new);
72 #ifdef notdef
73 NIpzMuller(set, new);
74 #endif
75 }
76
77
NIpzMuller(set,newtry)78 NIpzMuller(set, newtry)
79 PZtrial *set[3], *newtry;
80 {
81 SPcomplex A, B, C, D, E;
82 SPcomplex h0, h1;
83 SPcomplex j0, j1;
84 SPcomplex lambda_i, delta_i, ds;
85 double scale[3];
86 double q;
87 int mag[3], magx, min;
88 int i, j, error, total;
89
90 min = -999999;
91 j = 0;
92 total = 0;
93 for (i = 0; i < 3; i++) {
94 if (set[i]->f_def.real != 0.0 || set[i]->f_def.imag != 0.0) {
95 if (min < set[i]->mag_def - 50)
96 min = set[i]->mag_def - 50;
97 total += set[i]->mag_def;
98 j += 1;
99 }
100 }
101
102 magx = total / j;
103 if (magx < min)
104 magx = min;
105
106 DEBUG(2) fprintf(stderr, "Base scale: %d / %d (0: %d, 1: %d, 2: %d)\n",
107 magx, min, set[0]->mag_def, set[1]->mag_def, set[2]->mag_def);
108
109 for (i = 0; i < 3; i++) {
110 mag[i] = set[i]->mag_def - magx;
111 scale[i] = 1.0;
112 while (mag[i] > 0) {
113 scale[i] *= 2.0;
114 mag[i] -= 1;
115 }
116 if (mag[i] < -90)
117 scale[i] = 0.0;
118 else {
119 while (mag[i] < 0) {
120 scale[i] /= 2.0;
121 mag[i] += 1;
122 }
123 }
124 }
125
126 C_SUBEQ(h0,set[0]->s,set[1]->s);
127 C_SUBEQ(h1,set[1]->s,set[2]->s);
128 C_DIVEQ(lambda_i,h0,h1);
129
130 /* Quadratic interpolation (Muller's method) */
131
132 C_EQ(delta_i,lambda_i);
133 delta_i.real += 1.0;
134
135 /* Quadratic coefficients A, B, C (Note: reciprocal form of eqn) */
136
137 /* A = lambda_i * (f[i-2] * lambda_i - f[i-1] * delta_i + f[i]) */
138 C_MULEQ(A,scale[2] * set[2]->f_def,lambda_i);
139 C_MULEQ(C,scale[1] * set[1]->f_def,delta_i);
140 C_SUB(A,C);
141 C_ADD(A,scale[0] * set[0]->f_def);
142 C_MUL(A,lambda_i);
143
144 /* B = f[i-2] * lambda_i * lambda_1 - f[i-1] * delta_i * delta_i
145 + f[i] * (lambda_i + delta_i) */
146 C_MULEQ(B,lambda_i,lambda_i);
147 C_MUL(B,scale[2] * set[2]->f_def);
148 C_MULEQ(C,delta_i,delta_i);
149 C_MUL(C,scale[1] * set[1]->f_def);
150 C_SUB(B,C);
151 C_ADDEQ(C,lambda_i,delta_i);
152 C_MUL(C,scale[0] * set[0]->f_def);
153 C_ADD(B,C);
154
155 /* C = delta_i * f[i] */
156 C_MULEQ(C,delta_i,scale[0] * set[0]->f_def);
157
158 while (FABS(A.real) > 1.0 || FABS(A.imag) > 1.0
159 || FABS(B.real) > 1.0 || FABS(B.imag) > 1.0
160 || FABS(C.real) > 1.0 || FABS(C.imag) > 1.0) {
161 A.real /= 2.0;
162 B.real /= 2.0;
163 C.real /= 2.0;
164 A.imag /= 2.0;
165 B.imag /= 2.0;
166 C.imag /= 2.0;
167 }
168
169 /* discriminate = B * B - 4 * A * C */
170 C_MULEQ(D,B,B);
171 C_MULEQ(E,4.0 * A,C);
172 C_SUB(D,E);
173
174 DEBUG(2) fprintf(stderr, " Discr: (%g,%g)\n",D.real, D.imag);
175 C_SQRT(D);
176 DEBUG(1) fprintf(stderr, " sqrtDiscr: (%g,%g)\n",D.real, D.imag);
177
178 #ifndef notdef
179 /* Maximize denominator */
180
181 DEBUG(1) fprintf(stderr, " B: (%g,%g)\n",B.real, B.imag);
182 /* Dot product */
183 q = B.real * D.real + B.imag * D.imag;
184 if (q > 0.0) {
185 DEBUG(1) fprintf(stderr, " add\n");
186 C_ADD(B,D);
187 } else {
188 DEBUG(1) fprintf(stderr, " sub\n");
189 C_SUB(B,D);
190 }
191
192 #else
193 /* For trapped zeros, the step should always be positive */
194 if (C.real >= 0.0) {
195 if (B.real < D.real) {
196 C_SUB(B,D);
197 } else {
198 C_ADD(B,D);
199 }
200 } else {
201 if (B.real > D.real) {
202 C_SUB(B,D);
203 } else {
204 C_ADD(B,D);
205 }
206 }
207 #endif
208
209 DEBUG(1) fprintf(stderr, " C: (%g,%g)\n", C.real, C.imag);
210 C_DIV(C,-0.5 * B);
211
212 newtry->next = NULL;
213
214 DEBUG(1) fprintf(stderr, " lambda: (%g,%g)\n",C.real, C.imag);
215 C_MULEQ(newtry->s,h0,C);
216
217 DEBUG(1) fprintf(stderr, " h: (%g,%g)\n", newtry->s.real, newtry->s.imag);
218
219 C_ADD(newtry->s,set[0]->s);
220
221 DEBUG(1) fprintf(stderr, "New try: (%g,%g)\n",
222 newtry->s.real, newtry->s.imag);
223
224 return OK;
225
226 }
227
228
NIpzSym2(set,new)229 NIpzSym2(set, new)
230 PZtrial *set[3], *new;
231 {
232 double a, b, c, x0, x1;
233 double dx0, dx1, d2x, diff;
234 double disc;
235 int a_mag, b_mag, c_mag;
236 int tmag;
237 int error;
238 int disc_mag;
239 int new_mag;
240
241 error = OK;
242
243 /*
244 NIpzK = 0.0;
245 NIpzK_mag = 0;
246 */
247
248 /* Solve for X = the distance from set[1], where X > 0 => root < set[1] */
249 dx0 = set[1]->s.real - set[0]->s.real;
250 dx1 = set[2]->s.real - set[1]->s.real;
251
252 x0 = (set[0]->s.real + set[1]->s.real) / 2.0;
253 x1 = (set[1]->s.real + set[2]->s.real) / 2.0;
254 /* d2x = x1 - x0; */
255 d2x = (set[2]->s.real - set[0]->s.real) / 2.0;
256
257 zaddeq(&a, &a_mag, set[1]->f_def.real, set[1]->mag_def,
258 -set[0]->f_def.real, set[0]->mag_def);
259
260 tmag = 0;
261 R_NORM(dx0,tmag);
262 a /= dx0;
263 a_mag -= tmag;
264 R_NORM(a,a_mag);
265
266 zaddeq(&b, &b_mag, set[2]->f_def.real, set[2]->mag_def,
267 -set[1]->f_def.real, set[1]->mag_def);
268
269 tmag = 0;
270 R_NORM(dx1,tmag);
271 b /= dx1;
272 b_mag -= tmag;
273 R_NORM(b,b_mag);
274
275 zaddeq(&c, &c_mag, b, b_mag, -a, a_mag);
276
277 tmag = 0;
278 R_NORM(d2x,tmag);
279 c /= d2x; /* = f'' */
280 c_mag -= tmag;
281 R_NORM(c,c_mag);
282
283 if (c == 0.0 || (a == 0.0 || c_mag < a_mag - 40)
284 && (b = 0.0 ||c_mag < b_mag - 40)) {
285 /*fprintf(stderr, "\t- linear (%g, %d)\n", c, c_mag);*/
286 if (a == 0.0) {
287 a = b;
288 a_mag = b_mag;
289 }
290 if (a != 0.0) {
291 new->s.real = - set[1]->f_def.real / a;
292 a_mag -= set[1]->mag_def;
293 while (a_mag > 0) {
294 new->s.real /= 2.0;
295 a_mag -= 1;
296 }
297 while (a_mag < 0) {
298 new->s.real *= 2.0;
299 a_mag += 1;
300 }
301 new->s.real += set[1]->s.real;
302 } else
303 new->s.real = set[1]->s.real;
304 } else {
305
306 /* Quadratic power series about set[1]->s.real */
307 /* c : d2f/dx2 @ s1 (assumed constant for all s), or "2A" */
308
309 /* a : (df/dx) / (d2f/dx2) @ s1, or "B/2A" */
310 a /= c;
311 R_NORM(a,a_mag);
312 a_mag -= c_mag;
313
314 diff = set[1]->s.real - x0;
315 tmag = 0;
316 R_NORM(diff,tmag);
317
318 zaddeq(&a, &a_mag, a, a_mag, diff, tmag);
319
320 /* b : f(s1) / (1/2 d2f/ds2), or "C / A" */
321 b = 2.0 * set[1]->f_def.real / c;
322 b_mag = set[1]->mag_def - c_mag;
323 R_NORM(b,b_mag);
324
325 disc = a * a;
326 disc_mag = 2 * a_mag;
327
328 /* disc = a^2 - b :: (B/2A)^2 - C/A */
329 zaddeq(&disc, &disc_mag, disc, disc_mag, - b, b_mag);
330
331 if (disc < 0.0) {
332 /* Look for minima instead, but save radical for later work */
333 disc *= -1;
334 new_mag = 1;
335 }
336
337 if (disc_mag % 2 == 0)
338 disc = sqrt(disc);
339 else {
340 disc = sqrt(2.0 * disc);
341 disc_mag -= 1;
342 }
343 disc_mag /= 2;
344
345 if (new_mag != 0) {
346 if (NIpzK == 0.0) {
347 NIpzK = disc;
348 NIpzK_mag = disc_mag;
349 DEBUG(1) fprintf(stderr, "New NIpzK: %g*2^%d\n",
350 NIpzK, NIpzK_mag);
351 } else {
352 DEBUG(1) fprintf(stderr,
353 "Ignore NIpzK: %g*2^%d for previous value of %g*2^%d\n",
354 disc, disc_mag,
355 NIpzK, NIpzK_mag);
356 }
357 disc = 0.0;
358 disc_mag = 0;
359 }
360
361 /* NOTE: c & b get reused here-after */
362
363 if (a * disc >= 0.0) {
364 zaddeq(&c, &c_mag, a, a_mag, disc, disc_mag);
365 } else {
366 zaddeq(&c, &c_mag, a, a_mag, -disc, disc_mag);
367 }
368
369 /* second root = C / (first root) */
370 if (c != 0.0) {
371 b /= c;
372 b_mag -= c_mag;
373 } else {
374 /* special case */
375 b = 0.0;
376 b = 0;
377 }
378
379 zaddeq(&b, &b_mag, set[1]->s.real, 0, -b, b_mag);
380 zaddeq(&c, &c_mag, set[1]->s.real, 0, -c, c_mag);
381
382 while (b_mag > 0) {
383 b *= 2.0;
384 b_mag -= 1;
385 }
386 while (b_mag < 0) {
387 b /= 2.0;
388 b_mag += 1;
389 }
390
391 while (c_mag > 0) {
392 c *= 2.0;
393 c_mag -= 1;
394 }
395 while (c_mag < 0) {
396 c /= 2.0;
397 c_mag += 1;
398 }
399
400 DEBUG(1) fprintf(stderr, "@@@ (%.15g) -vs- (%.15g)\n", b, c);
401 /* XXXX */
402 if (b < set[0]->s.real || b > set[2]->s.real) {
403 /* b not in range */
404 if (c < set[0]->s.real || c > set[2]->s.real) {
405 DEBUG(1) fprintf(stderr, "@@@ both are junk\n");
406 if (CKTpzTrapped == 1)
407 new->s.real = (set[0]->s.real + set[1]->s.real) / 2.0;
408 else if (CKTpzTrapped == 2)
409 new->s.real = (set[1]->s.real + set[2]->s.real) / 2.0;
410 else if (CKTpzTrapped == 3) {
411 if (FABS(set[1]->s.real - c) < FABS(set[1]->s.real - b)) {
412 DEBUG(1) fprintf(stderr, "@@@ mix w/second (c)\n");
413 new->s.real = (set[1]->s.real + c) / 2.0;
414 } else {
415 DEBUG(1) fprintf(stderr, "@@@ mix w/first (b)\n");
416 new->s.real = (set[1]->s.real + b) / 2.0;
417 }
418 } else
419 ERROR(E_PANIC,"Lost numerical stability");
420 } else {
421 DEBUG(1) fprintf(stderr, "@@@ take second (c)\n");
422 new->s.real = c;
423 }
424 } else {
425 /* b in range */
426 if (c < set[0]->s.real || c > set[2]->s.real) {
427 DEBUG(1) fprintf(stderr, "@@@ take first (b)\n");
428 new->s.real = b;
429 } else {
430 /* Both in range -- take the smallest mag */
431 if (a > 0.0) {
432 DEBUG(1) fprintf(stderr, "@@@ push -- first (b)\n");
433 new->s.real = b;
434 } else {
435 DEBUG(1) fprintf(stderr, "@@@ push -- first (b)\n");
436 new->s.real = c;
437 }
438 }
439 }
440
441 }
442
443 new->s.imag = 0.0;
444
445 return error;
446 }
447