null.c

Polynomial Root Finder is a reliable and fast C program (+ Matlab gateway) for finding all roots of

/******************************************************************/
/*                                                                */
/* file:          null.c                                          */
/*                                                                */
/* main function: null()                                          */
/*                                                                */
/* version:       1.3                                             */
/*                                                                */
/* author:        B. Frenzel                                      */
/*                                                                */
/* date:          Jan 19 1993                                     */ 
/*                                                                */
/* input:         p[]     coefficient vector of the original      */
/*                        polynomial                              */
/*                pred[]  coefficient vector of the deflated      */
/*                        polynomial                              */
/*                *n      the highest exponent of the original    */
/*                        polynomial                              */  
/*                flag    flag = TRUE  => complex coefficients    */
/*                        flag = FALSE => real    coefficients    */
/*                                                                */
/* output:        root[]  vector of determined roots              */
/*                *maxerr estimation of max. error of all         */
/*                        determined roots                        */
/*                                                                */
/* subroutines:   poly_check(),quadratic(),lin_or_quad(),monic(), */
/*                muller(),newton()                               */ 
/*                                                                */
/* description:                                                   */
/* main rootfinding routine for polynomials with complex or real  */
/* coefficients using Muller's method combined with Newton's m.   */
/*                                                                */
/* Copyright:                                                     */
/* Lehrstuhl fuer Nachrichtentechnik Erlangen                     */
/* Cauerstr. 7, 8520 Erlangen, FRG, 1993                          */
/* e-mail: [e-mail]                       */
/*                                                                */
/******************************************************************/
#define  NUL
#include "header.h"










/***** main function null() *****/
unsigned char null(dcomplex *p,dcomplex *pred,int *n,dcomplex *root,
                   double *maxerr,unsigned char flag)
/*dcomplex *p,     coefficient vector of the original polynomial   */
/*         *pred,  coefficient vector of the deflated polynomial   */
/*         *root;  determined roots                                */
/*int      *n;     the highest exponent of the original polynomial */
/*double   *maxerr; max. error of all determined roots             */ 
/*unsigned char flag;  flag = TRUE  => complex coefficients        */
/*                     flag = FALSE => real    coefficients        */
{
     double   newerr; /* error of actual root                      */
     dcomplex ns;     /* root determined by Muller's method        */
     int      nred,   /* highest exponent of deflated polynomial   */
              i;      /* counter                                   */
     unsigned char error; /* indicates an error in poly_check      */
     int      red,
              diff;   /* number of roots at 0                      */

     *maxerr = 0.;    /* initialize max. error of determined roots */
     nred = *n;       /* At the beginning: degree defl. polyn. =   */
                      /* degree of original polyn.                 */

                      /* check input of the polynomial             */
     error = poly_check(p,&nred,n,root);
     diff  = (*n-nred); /* reduce polynomial, if roots at 0        */
     p    += diff;      
     *n   =  nred;

     if (error)
          return error; /* error in poly_check(); return error     */ 

                        /* polynomial is linear or quadratic       */
     if (lin_or_quad(p,nred,root)==0) { 
          *n += diff;     /* remember roots at 0                   */
          *maxerr = DBL_EPSILON; 
          return 0;       /* return no error                       */
     }

     monic(p,n);          /* get monic polynom                     */

     for (i=0;i<=*n;i++)  pred[i]=p[i];  /* original polynomial    */
                           /* = deflated polynomial at beginning   */
                           /* of Muller                            */

     do {                  /* main loop of null()                  */ 
                           /* Muller method                        */
          ns = muller(pred,nred); 



                           /* Newton method                        */
          root[nred-1] = newton(p,*n,ns,&newerr,flag);

                           /* stores max. error of all roots       */
          if (newerr>*maxerr)  
               *maxerr=newerr;
                           /* deflate polynomial                   */
          red = poldef(pred,nred,root,flag);
          pred += red;        /* forget lowest coefficients        */
          nred -= red;        /* reduce degree of polynomial       */
     } while (nred>2); 
                              /* last one or two roots             */
     (void) lin_or_quad(pred,nred,root);
     if (nred==2) {
          root[1] = newton(p,*n,root[1],&newerr,flag);
          if (newerr>*maxerr)  
               *maxerr=newerr;
     }    
     root[0] = newton(p,*n,root[0],&newerr,flag);
     if (newerr>*maxerr)  
          *maxerr=newerr;

    *n += diff;              /* remember roots at 0               */
    return 0;                /* return no error                   */
}










/***** poly_check() check the formal correctness of input *****/
unsigned char poly_check(dcomplex *pred,int *nred,int *n,dcomplex *root)
/*dcomplex *pred,  coefficient vector of the original polynomial   */
/*         *root;  determined roots                                */
/*int      *nred,  highest exponent of the deflated polynomial     */
/*         *n;     highest exponent of the original polynomial     */
{
     int  i = -1, /* i stores the (reduced) real degree            */
          j;      /* counter variable                              */
     unsigned char
           notfound=TRUE; /* indicates, whether a coefficient      */
                          /* unequal zero was found                */

     if (*n<0) return 1;  /* degree of polynomial less than zero   */
                          /* return error                          */

     for (j=0;j<=*n;j++) {      /* determines the "real" degree of       */
          if(Cabs(pred[j])!=0.) /* polynomial, cancel leading roots      */
               i=j; 
     }
     if (i==-1) return 2;   /* polynomial is a null vector; return error */
     if (i==0) return 3;    /* polynomial is constant unequal null;      */
                            /* return error                              */

     *n=i;                  /* set new exponent of polynomial            */
     i=0;                   /* reset variable for exponent               */
     do {                   /* find roots at 0                           */
          if (Cabs(pred[i])==0.) 
               i++; 
          else 
               notfound=FALSE;
     } while (i<=*n && notfound);

          if (i==0) {            /* no root determined at 0              */
               *nred = *n;       /* original degree=deflated degree and  */
               return 0;         /* return no error                      */
          } else {               /* roots determined at 0:               */
               for (j=0;j<=i-1;j++) /* store roots at 0                  */ 
                    root[*n-j-1] = Complex(0.,0.); 
               *nred = *n-i;  /* reduce degree of deflated polynomial    */
               return 0;      /* and return no error                     */
          }
}










/***** quadratic() calculates the roots of a quadratic polynomial *****/
void quadratic(dcomplex *pred,dcomplex *root)
/*dcomplex *pred,  coefficient vector of the deflated polynomial   */
/*         *root;  determined roots                                */

{
     dcomplex discr,       /* discriminate                         */ 
              Z1,Z2,       /* numerators of the quadratic formula  */
              N;           /* denominator of the quadratic formula */
 
                           /* discr = p1^2-4*p2*p0                 */   
     discr   = Csub(Cmul(pred[1],pred[1]),
               RCmul(4.,Cmul(pred[2],pred[0])));
                           /* Z1 = -p1+sqrt(discr)                 */
     Z1      = Cadd(RCmul(-1.,pred[1]),Csqrt(discr));
                           /* Z2 = -p1-sqrt(discr)                 */
     Z2      = Csub(RCmul(-1.,pred[1]),Csqrt(discr));
                           /* N  = 2*p2                            */
     N       = RCmul(2.,pred[2]);
     root[0] = Cdiv(Z1,N); /* first root  = Z1/N                   */
     root[1] = Cdiv(Z2,N); /* second root = Z2/N                   */
}










/***** lin_or_quad() calculates roots of lin. or quadratic equation *****/
unsigned char lin_or_quad(dcomplex *pred,int nred,dcomplex *root)
/*dcomplex *pred,  coefficient vector of the deflated polynomial   */
/*         *root;  determined roots                                */
/*int      nred;   highest exponent of the deflated polynomial     */
{
     if (nred==1) {     /* root = -p0/p1                           */
          root[0] = Cdiv(RCmul(-1.,pred[0]),pred[1]);
          return 0;     /* and return no error                     */
     } else if (nred==2) { /* quadratic polynomial                 */
          quadratic(pred,root);
          return 0;        /* return no error                      */
     }

     return 1; /* nred>2 => no roots were calculated               */
}










/***** monic() computes monic polynomial for original polynomial *****/
void monic(dcomplex *p,int *n)
/*dcomplex *p;     coefficient vector of the original polynomial   */
/*int      *n;     the highest exponent of the original polynomial */
{
     double factor;   /* stores absolute value of the coefficient  */
                      /* with highest exponent                     */
     int    i;        /* counter variable                          */

     factor=1./Cabs(p[*n]);     /* factor = |1/pn|                 */
     if ( factor!=1.)           /* get monic pol., when |pn| != 1  */
         for (i=0;i<=*n;i++) 
               p[i]=RCmul(factor,p[i]);
}