home *** CD-ROM | disk | FTP | other *** search
/ Hacker Chronicles 2 / HACKER2.BIN / 737.POLY.CPP < prev    next >
Text File  |  1993-04-17  |  6KB  |  209 lines
  1. //  poly.cpp - a public domain polynomial root finder
  2. //             by D.A. Steinman (steinman@me.utoronto.ca)
  3. //
  4. //  As part of a project I was working on recently, I required a
  5. //  robust 4th order polynomial root solver.  Simple iterative
  6. //  methods (e.g. Bairstow) don't always work, and most good black
  7. //  box routines are written in Fortran and cost $$$.  However, since
  8. //  I was concerned only with the 4th order case, I knew it could be
  9. //  solved exactly.  The code which follows represents the fruit of my
  10. //  labours...
  11. //
  12. //  POLY solves the general polynomial:
  13. //
  14. //      a[n] x^n + a[n-1] x^n-1 ... a[1] x + a[0] = 0  n<=4
  15. //
  16. //  where the coefficients a[i] are real and, for numerical
  17. //  convenience, the leading term is normalized to unity by dividing
  18. //  each coefficient by a[n].  Polynomials of order 4 or less (i.e.
  19. //  n<=4) can be solved exactly, and this is what POLY will do for
  20. //  you.  No messy iterations.  No worries about multiple roots.  The
  21. //  routines for the cubic and quartic cases are based upon methods
  22. //  outlined in:
  23. //
  24. //      Tignol J-P.  Galois' theory of algebraic equations.  Longman
  25. //          Scientific & Technical, Harlow, 1988  (ISBN 0-470-20919-4)
  26. //
  27. //  POLY was written in C++ (for the complex variables), and was
  28. //  compiled using Turbo C++ v1.0; I haven't tested it for any other
  29. //  flavour of C++.  POLY *should* solve any quartic or lower order
  30. //  polynomial.  If it doesn't, feel free to find the bug and post a
  31. //  fixed version.
  32.  
  33. #include <stdio.h>
  34. #include <complex.h>
  35. #include <math.h>
  36. #include <stdlib.h>
  37.  
  38. #define N  5
  39. #define SIGN(a)  (a>=0 ? "+" : "-")
  40.  
  41. void poly1(double *, complex *);
  42. void poly2(double *, complex *);
  43. void poly3(double *, complex *);
  44. void poly4(double *, complex *);
  45.  
  46. void main(int argc, char *argv[])
  47. {
  48.     complex x[N],resid,tx;
  49.     double a[N];
  50.     int i,j,k,n,tr;
  51.  
  52.     /* print syntax message */
  53.     n = (--argc)-1;
  54.     if (n<1 || n>4) {
  55.         puts("");
  56.         puts("poly - zeros of 4th order (and lower) polynomials");
  57.         puts("       v1.00 public domain by D.A. Steinman");
  58.         puts("");
  59.         puts("syntax:  poly  a[n]..a[0] (n=[1..4])");
  60.         exit(1);
  61.     }
  62.  
  63.     /* initialize coefficients and roots */
  64.     for (i=0;i<=n;i++) a[i] = 0.0;
  65.     for (i=1;i<=n;i++) x[i] = complex(0.0);
  66.  
  67.     /* assign coefficients */
  68.     for (i=n;i>=0;i--) a[i]=atof(argv[n+1-i]);
  69.     if (fabs(a[n])==0.0) {
  70.         puts("\n*error* - a[n] must be non-zero");
  71.         exit(2);
  72.     }
  73.     for (i=0;i<=n;i++) a[i] /= a[n];
  74.  
  75.     /* write polynomial */
  76.     printf("\np(x):  ");
  77.     if (n>1) printf("x^%1d",n);
  78.     else printf("x");
  79.     for (i=n-1;i>1;i--) if (fabs(a[i])>0.0) printf(" %s %lg x^%1d",SIGN(a[i]),fabs(a[i]),i);
  80.     if (fabs(a[1])>0.0 && n>1) printf(" %s %lg x",SIGN(a[1]),fabs(a[1]));
  81.     if (fabs(a[0])>0.0) printf(" %s %lg",SIGN(a[0]),fabs(a[0]));
  82.     puts(" = 0\n");
  83.  
  84.     /* route solver based upon order */
  85.     switch(n) {
  86.       case 1: 
  87.         poly1(&a[0],&x[0]);              // linear
  88.         break;
  89.       case 2:
  90.         if (fabs(a[0])==0.0)
  91.             poly1(&a[1],&x[1]);          // degenerate linear
  92.         else
  93.             poly2(&a[0],&x[0]);          // quadratic
  94.         break;
  95.       case 3:
  96.         if (fabs(a[0])==0.0)
  97.             if (fabs(a[1])==0.0)
  98.                 poly1(&a[2],&x[2]);      // degenerate linear
  99.             else
  100.                 poly2(&a[1],&x[1]);      // degenerate quadratic
  101.         else
  102.             poly3(&a[0],&x[0]);          // cubic
  103.         break;
  104.       case 4:
  105.         if (fabs(a[0])==0.0)
  106.             if (fabs(a[1])==0.0)
  107.                 if (fabs(a[2])==0.0)
  108.                     poly1(&a[3],&x[3]);  // degenerate linear
  109.                 else
  110.                     poly2(&a[2],&x[2]);  // degenerate quadratic
  111.             else
  112.                 poly3(&a[1],&x[1]);      // degenerate cubic
  113.         else
  114.             poly4(&a[0],&x[0]);          // quartic
  115.         break;
  116.     }
  117.  
  118.     /* sort by descending real part */
  119.     for (i=1;i<=n;i++) for (j=i+1;j<=n;j++)
  120.         if (real(x[j]) > real(x[i])) {
  121.             tx = x[i];
  122.             x[i] = x[j];
  123.             x[j] = tx;
  124.         }
  125.  
  126.     /* display roots */
  127.     for (i=1;i<=n;i++) {
  128.         resid = a[0];
  129.         for (j=1;j<=n;j++) resid += a[j]*pow(x[i],double(j));
  130.         printf("x[%1d]: % lf%+lfi   p(x) = %lg\n",i,real(x[i]),imag(x[i]),abs(resid));
  131.     }
  132. }
  133.  
  134. void poly1(double *a, complex *x)
  135. {
  136.     x[1] = -a[0];
  137. }
  138.  
  139. void poly2(double *a, complex *x)
  140. {
  141.     complex t0;
  142.  
  143.     t0 = sqrt(complex(a[1]*a[1]-4*a[0]));
  144.     x[1] = -0.5*(a[1]-t0);
  145.     x[2] = -0.5*(a[1]+t0);
  146. }
  147.  
  148. void poly3(double *a, complex *x)
  149. {
  150.     complex y[N],z[N];
  151.     double p,q;
  152.     int i;
  153.  
  154.     /* eliminate x^2 term via y=x+a[2]/3 */
  155.     p = a[1] - a[2]*a[2]/3;
  156.     q = a[0] - a[2]*a[1]/3 + 2*a[2]*a[2]*a[2]/27;
  157.  
  158.     if (fabs(p)==0.0) {  // special case of p=0 ==> x^3 = -q
  159.        y[1] = pow(complex(-q),1.0/3.0);
  160.        for (i=2;i<=3;i++) y[i] = y[i-1]*0.5*complex(-1,sqrt(3.0));
  161.     } else if (fabs(q)==0.0) {  // special case of q=0 ==> x^2 = -p
  162.        y[1] = 0.0;
  163.        y[2] = sqrt(complex(-p));
  164.        y[3] = -sqrt(complex(-p));
  165.     } else {
  166.        z[1] = pow(q/2+sqrt(complex(p*p*p/27+q*q/4)),1.0/3.0);
  167.        for (i=2;i<=3;i++) z[i] = z[i-1]*0.5*complex(-1,sqrt(3.0));
  168.        for (i=1;i<=3;i++) y[i] = p/3/z[i]-z[i];
  169.     }
  170.     for (i=1;i<=3;i++) x[i] = y[i]-a[2]/3;
  171. }
  172.  
  173. void poly4(double *a, complex *x)
  174. {
  175.     complex ux[N],y[N],t0,t1;
  176.     double ua[N],p,q,r;
  177.     int i;
  178.  
  179.     /* eliminate x^3 term via y=x+a[3]/4 */
  180.     p = a[2] - 3*a[3]*a[3]/8;
  181.     q = a[1] - a[3]*a[2]/2 + a[3]*a[3]*a[3]/8;
  182.     r = a[0] - a[3]*a[1]/4 + a[3]*a[3]*a[2]/16 - 3*a[3]*a[3]*a[3]*a[3]/256;
  183.  
  184.     /* solve resolvent cubic */
  185.     ua[2] = p;
  186.     ua[1] = p*p/4-r;
  187.     ua[0] = -q*q/8;
  188.     poly3(ua,ux);
  189.  
  190.     if (fabs(q)==0.0) {  // special case of p=0 ==> x^3 = -q
  191.         t0 = sqrt(complex(p*p/4-r));
  192.         t1 = sqrt(-p/2+t0);
  193.         y[1] =  t1;
  194.         y[2] = -t1;
  195.         t1 = sqrt(-p/2-t0);
  196.         y[3] =  t1;
  197.         y[4] = -t1;
  198.     } else {
  199.         t0 = sqrt(ux[1]/2);
  200.         t1 = sqrt(-ux[1]/2-p/2-q/2/sqrt(2*ux[1]));
  201.         y[1] =  t0+t1;
  202.         y[2] =  t0-t1;
  203.         t1 = sqrt(-ux[1]/2-p/2+q/2/sqrt(2*ux[1]));
  204.         y[3] = -t0+t1;
  205.         y[4] = -t0-t1;
  206.     }
  207.     for (i=1;i<=4;i++) x[i] = y[i]-a[3]/4;
  208. }
  209.