123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270 |
- // find the roots of a polynomial numerically
- #include "stdafx.h"
- #include "defs.h"
- #define YMAX 101
- #define DELTA 1.0e-6
- #define EPSILON 1.0e-9
- #define ABS(z) sqrt((z).r * (z).r + (z).i * (z).i)
- #define RANDOM (4.0 * (double) rand() / (double) RAND_MAX - 2.0)
- static struct {
- double r, i;
- } a, b, x, y, fa, fb, dx, df, c[YMAX];
- void
- eval_nroots(void)
- {
- int h, i, k, n;
- push(cadr(p1));
- eval();
- push(caddr(p1));
- eval();
- p2 = pop();
- if (p2 == symbol(NIL))
- guess();
- else
- push(p2);
- p2 = pop();
- p1 = pop();
- if (!ispoly(p1, p2))
- stop("nroots: polynomial?");
- // mark the stack
- h = tos;
- // get the coefficients
- push(p1);
- push(p2);
- n = coeff();
- if (n > YMAX)
- stop("nroots: degree?");
- // convert the coefficients to real and imaginary doubles
- for (i = 0; i < n; i++) {
- push(stack[h + i]);
- real();
- yyfloat();
- eval();
- p1 = pop();
- push(stack[h + i]);
- imag();
- yyfloat();
- eval();
- p2 = pop();
- if (!isdouble(p1) || !isdouble(p2))
- stop("nroots: coefficients?");
- c[i].r = p1->u.d;
- c[i].i = p2->u.d;
- }
- // pop the coefficients
- tos = h;
- // n is the number of coefficients, n = deg(p) + 1
- monic(n);
- for (k = n; k > 1; k--) {
- findroot(k);
- if (fabs(a.r) < DELTA)
- a.r = 0.0;
- if (fabs(a.i) < DELTA)
- a.i = 0.0;
- push_double(a.r);
- push_double(a.i);
- push(imaginaryunit);
- multiply();
- add();
- divpoly(k);
- }
- // now make n equal to the number of roots
- n = tos - h;
- if (n > 1) {
- sort_stack(n);
- p1 = alloc_tensor(n);
- p1->u.tensor->ndim = 1;
- p1->u.tensor->dim[0] = n;
- for (i = 0; i < n; i++)
- p1->u.tensor->elem[i] = stack[h + i];
- tos = h;
- push(p1);
- }
- }
- // divide the polynomial by its leading coefficient
- void
- monic(int n)
- {
- int k;
- double t;
- y = c[n - 1];
- t = y.r * y.r + y.i * y.i;
- for (k = 0; k < n - 1; k++) {
- c[k].r = (c[k].r * y.r + c[k].i * y.i) / t;
- c[k].i = (c[k].i * y.r - c[k].r * y.i) / t;
- }
- c[n - 1].r = 1.0;
- c[n - 1].i = 0.0;
- }
- // uses the secant method
- void
- findroot(int n)
- {
- int j, k;
- double t;
- if (ABS(c[0]) < DELTA) {
- a.r = 0.0;
- a.i = 0.0;
- return;
- }
- for (j = 0; j < 100; j++) {
- a.r = RANDOM;
- a.i = RANDOM;
- compute_fa(n);
- b = a;
- fb = fa;
- a.r = RANDOM;
- a.i = RANDOM;
- for (k = 0; k < 1000; k++) {
- compute_fa(n);
- if (ABS(fa) < EPSILON)
- return;
- if (ABS(fa) < ABS(fb)) {
- x = a;
- a = b;
- b = x;
- x = fa;
- fa = fb;
- fb = x;
- }
- // dx = b - a
- dx.r = b.r - a.r;
- dx.i = b.i - a.i;
- // df = fb - fa
- df.r = fb.r - fa.r;
- df.i = fb.i - fa.i;
- // y = dx / df
- t = df.r * df.r + df.i * df.i;
- if (t == 0.0)
- break;
- y.r = (dx.r * df.r + dx.i * df.i) / t;
- y.i = (dx.i * df.r - dx.r * df.i) / t;
- // a = b - y * fb
- a.r = b.r - (y.r * fb.r - y.i * fb.i);
- a.i = b.i - (y.r * fb.i + y.i * fb.r);
- }
- }
- stop("nroots: convergence error");
- }
- void
- compute_fa(int n)
- {
- int k;
- double t;
- // x = a
- x.r = a.r;
- x.i = a.i;
- // fa = c0 + c1 * x
- fa.r = c[0].r + c[1].r * x.r - c[1].i * x.i;
- fa.i = c[0].i + c[1].r * x.i + c[1].i * x.r;
- for (k = 2; k < n; k++) {
- // x = a * x
- t = a.r * x.r - a.i * x.i;
- x.i = a.r * x.i + a.i * x.r;
- x.r = t;
- // fa += c[k] * x
- fa.r += c[k].r * x.r - c[k].i * x.i;
- fa.i += c[k].r * x.i + c[k].i * x.r;
- }
- }
- // divide the polynomial by x - a
- void
- divpoly(int n)
- {
- int k;
- for (k = n - 1; k > 0; k--) {
- c[k - 1].r += c[k].r * a.r - c[k].i * a.i;
- c[k - 1].i += c[k].i * a.r + c[k].r * a.i;
- }
- if (ABS(c[0]) > DELTA)
- stop("nroots: residual error");
- for (k = 0; k < n - 1; k++) {
- c[k].r = c[k + 1].r;
- c[k].i = c[k + 1].i;
- }
- }
- #if SELFTEST
- static char *s[] = {
- "nroots(x)",
- "0",
- "nroots((1+i)*x^2+1)",
- "(-0.17178-0.727673*i,0.17178+0.727673*i)",
- "nroots(sqrt(2)*exp(i*pi/4)*x^2+1)",
- "(-0.17178-0.727673*i,0.17178+0.727673*i)",
- // "nroots(x^4+1)",
- // "(-0.707107+0.707107*i,-0.707107-0.707107*i,0.707107+0.707107*i,0.707107-0.707107*i)",
- };
- void
- test_nroots(void)
- {
- test(__FILE__, s, sizeof s / sizeof (char *));
- }
- #endif
|