libm/double/polrt.c

   1 /*                                                      polrt.c
   2  *
   3  *      Find roots of a polynomial
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * typedef struct
  10  *      {
  11  *      double r;
  12  *      double i;
  13  *      }cmplx;
  14  *
  15  * double xcof[], cof[];
  16  * int m;
  17  * cmplx root[];
  18  *
  19  * polrt( xcof, cof, m, root )
  20  *
  21  *
  22  *
  23  * DESCRIPTION:
  24  *
  25  * Iterative determination of the roots of a polynomial of
  26  * degree m whose coefficient vector is xcof[].  The
  27  * coefficients are arranged in ascending order; i.e., the
  28  * coefficient of x**m is xcof[m].
  29  *
  30  * The array cof[] is working storage the same size as xcof[].
  31  * root[] is the output array containing the complex roots.
  32  *
  33  *
  34  * ACCURACY:
  35  *
  36  * Termination depends on evaluation of the polynomial at
  37  * the trial values of the roots.  The values of multiple roots
  38  * or of roots that are nearly equal may have poor relative
  39  * accuracy after the first root in the neighborhood has been
  40  * found.
  41  *
  42  */
  43 \f
  44 /*                                                      polrt   */
  45 /* Complex roots of real polynomial */
  46 /* number of coefficients is m + 1 ( i.e., m is degree of polynomial) */
  47
  48 #include <math.h>
  49 /*
  50 typedef struct
  51         {
  52         double r;
  53         double i;
  54         }cmplx;
  55 */
  56 #ifdef ANSIPROT
  57 extern double fabs ( double );
  58 #else
  59 double fabs();
  60 #endif
  61
  62 int polrt( xcof, cof, m, root )
  63 double xcof[], cof[];
  64 int m;
  65 cmplx root[];
  66 {
  67 register double *p, *q;
  68 int i, j, nsav, n, n1, n2, nroot, iter, retry;
  69 int final;
  70 double mag, cofj;
  71 cmplx x0, x, xsav, dx, t, t1, u, ud;
  72
  73 final = 0;
  74 n = m;
  75 if( n <= 0 )
  76         return(1);
  77 if( n > 36 )
  78         return(2);
  79 if( xcof[m] == 0.0 )
  80         return(4);
  81
  82 n1 = n;
  83 n2 = n;
  84 nroot = 0;
  85 nsav = n;
  86 q = &xcof[0];
  87 p = &cof[n];
  88 for( j=0; j<=nsav; j++ )
  89         *p-- = *q++;    /*      cof[ n-j ] = xcof[j];*/
  90 xsav.r = 0.0;
  91 xsav.i = 0.0;
  92
  93 nxtrut:
  94 x0.r = 0.00500101;
  95 x0.i = 0.01000101;
  96 retry = 0;
  97
  98 tryagn:
  99 retry += 1;
 100 x.r = x0.r;
 101
 102 x0.r = -10.0 * x0.i;
 103 x0.i = -10.0 * x.r;
 104
 105 x.r = x0.r;
 106 x.i = x0.i;
 107
 108 finitr:
 109 iter = 0;
 110
 111 while( iter < 500 )
 112 {
 113 u.r = cof[n];
 114 if( u.r == 0.0 )
 115         {               /* this root is zero */
 116         x.r = 0;
 117         n1 -= 1;
 118         n2 -= 1;
 119         goto zerrut;
 120         }
 121 u.i = 0;
 122 ud.r = 0;
 123 ud.i = 0;
 124 t.r = 1.0;
 125 t.i = 0;
 126 p = &cof[n-1];
 127 for( i=0; i<n; i++ )
 128         {
 129         t1.r = x.r * t.r  -  x.i * t.i;
 130         t1.i = x.r * t.i  +  x.i * t.r;
 131         cofj = *p--;            /* evaluate polynomial */
 132         u.r += cofj * t1.r;
 133         u.i += cofj * t1.i;
 134         cofj = cofj * (i+1);    /* derivative */
 135         ud.r += cofj * t.r;
 136         ud.i -= cofj * t.i;
 137         t.r = t1.r;
 138         t.i = t1.i;
 139         }
 140
 141 mag = ud.r * ud.r  +  ud.i * ud.i;
 142 if( mag == 0.0 )
 143         {
 144         if( !final )
 145                 goto tryagn;
 146         x.r = xsav.r;
 147         x.i = xsav.i;
 148         goto findon;
 149         }
 150 dx.r = (u.i * ud.i  -  u.r * ud.r)/mag;
 151 x.r += dx.r;
 152 dx.i = -(u.r * ud.i  +  u.i * ud.r)/mag;
 153 x.i += dx.i;
 154 if( (fabs(dx.i) + fabs(dx.r)) < 1.0e-6 )
 155         goto lupdon;
 156 iter += 1;
 157 }       /* while iter < 500 */
 158
 159 if( final )
 160         goto lupdon;
 161 if( retry < 5 )
 162         goto tryagn;
 163 return(3);
 164
 165 lupdon:
 166 /* Swap original and reduced polynomials */
 167 q = &xcof[nsav];
 168 p = &cof[0];
 169 for( j=0; j<=n2; j++ )
 170         {
 171         cofj = *q;
 172         *q-- = *p;
 173         *p++ = cofj;
 174         }
 175 i = n;
 176 n = n1;
 177 n1 = i;
 178
 179 if( !final )
 180         {
 181         final = 1;
 182         if( fabs(x.i/x.r) < 1.0e-4 )
 183                 x.i = 0.0;
 184         xsav.r = x.r;
 185         xsav.i = x.i;
 186         goto finitr;    /* do final iteration on original polynomial */
 187         }
 188
 189 findon:
 190 final = 0;
 191 if( fabs(x.i/x.r) >= 1.0e-5 )
 192         {
 193         cofj = x.r + x.r;
 194         mag = x.r * x.r  +  x.i * x.i;
 195         n -= 2;
 196         }
 197 else
 198         {               /* root is real */
 199 zerrut:
 200         x.i = 0;
 201         cofj = x.r;
 202         mag = 0;
 203         n -= 1;
 204         }
 205 /* divide working polynomial cof(z) by z - x */
 206 p = &cof[1];
 207 *p += cofj * *(p-1);
 208 for( j=1; j<n; j++ )
 209         {
 210         *(p+1) += cofj * *p  -  mag * *(p-1);
 211         p++;
 212         }
 213
 214 setrut:
 215 root[nroot].r = x.r;
 216 root[nroot].i = x.i;
 217 nroot += 1;
 218 if( mag != 0.0 )
 219         {
 220         x.i = -x.i;
 221         mag = 0;
 222         goto setrut;    /* fill in the complex conjugate root */
 223         }
 224 if( n > 0 )
 225         goto nxtrut;
 226 return(0);
 227 }