tools.c

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

/******************************************************************/
/*                                                                */
/* file:          tools.c                                         */
/*                                                                */
/* version:       1.2                                             */
/*                                                                */
/* author:        B. Frenzel                                      */
/*                                                                */
/* date:          Jan 7 1993                                      */
/*                                                                */
/* function:      description:                                    */ 
/*                                                                */
/* hornc()        hornc() deflates the polynomial with coeff.     */
/*                stored in pred[0],...,pred[n] by Horner's       */
/*                method (one root)                               */
/*                                                                */
/* horncd()       horncd() deflates the polynomial with coeff.    */
/*                stored in pred[0],...,pred[n] by Horner's       */
/*                method (two roots)                              */
/*                                                                */
/* poldef()       decides whether to call hornc() or horncd()     */
/*                                                                */
/* fdvalue()      fdvalue() computes f=P(x0) by Horner's method;  */
/*                if flag=TRUE, it additionally computes the      */
/*                derivative df=P'(x0)                            */
/*                                                                */
/* Copyright:                                                     */
/* Lehrstuhl fuer Nachrichtentechnik Erlangen                     */
/* Cauerstr. 7, 8520 Erlangen, FRG, 1993                          */
/* e-mail: [e-mail]                       */
/*                                                                */
/******************************************************************/
#define  TOOLS
#include "header.h"










/***** Horner method to deflate one root *****/ 
void hornc(dcomplex *pred,int nred,dcomplex x0,unsigned char flag) 
/*dcomplex      *pred,  coefficient vector of the polynomial           */
/*              x0;     root to be deflated                            */
/*int           nred;   the highest exponent of (deflated) polynomial  */
/*unsigned char flag;   indicates how to reduce polynomial             */

{
     int      i;          /* counter                        */
     dcomplex help1;      /* help variable                  */

     if ((flag&1)==0)     /* real coefficients              */
          for(i=nred-1; i>0; i--) 
               pred[i].r += (x0.r*pred[i+1].r);
     else                 /* complex coefficients           */
          for (i=nred-1; i>0; i--) {
               CMUL(help1,pred[i+1],x0);
               CADD(pred[i],help1,pred[i]);
          }
}










/***** Horner method to deflate two roots *****/
void horncd(dcomplex *pred,int nred,double a,double b)
/*dcomplex *pred;     coefficient vector of the polynomial           */
/*double   a,         coefficients of the quadratic polynomial       */
/*         b;         x^2+ax+b                                       */
/*int      nred;      the highest exponent of (deflated) polynomial  */
{
     int i;        /* counter */

     pred[nred-1].r += pred[nred].r*a; 
     for (i=nred-2; i>1; i--)
          pred[i].r += (a*pred[i+1].r+b*pred[i+2].r);
}










/***** main routine to deflate polynomial *****/
int poldef(dcomplex *pred,int nred,dcomplex *root,unsigned char flag)
/*dcomplex *pred,     coefficient vector of the polynomial           */
/*         *root;     vector of determined roots                     */
/*int       nred;     the highest exponent of (deflated) polynomial  */
/*unsigned char flag;  indicates how to reduce polynomial            */
{
     double   a,   /* coefficients of the quadratic polynomial       */
              b;   /* x^2+ax+b                                       */
     dcomplex x0;  /* root to be deflated                            */


     x0 = root[nred-1];
     if (x0.i!=0.)             /* x0 is complex                      */
          flag |=2; 

     if (flag==2) {            /* real coefficients and complex root */
          a = 2*x0.r;          /* => deflate x0 and Conjg(x0)        */
          b = -(x0.r*x0.r+x0.i*x0.i); 
          root[nred-2]=Conjg(x0); /* store second root = Conjg(x0)   */
          horncd(pred,nred,a,b);
          return 2;            /* two roots deflated                 */
     } else {
          hornc(pred,nred,x0,flag); /* deflate only one root         */
          return 1;            /* one root deflated                  */
     }
}










/***** fdvalue computes P(x0) and optional P'(x0) *****/
void fdvalue(dcomplex *p,int n,dcomplex *f,dcomplex *df,dcomplex x0,
             unsigned char flag) 
/*dcomplex      *p,     coefficient vector of the polynomial P(x)   */
/*              *f,     the result f=P(x0)                          */
/*              *df,    the result df=P'(x0), if flag=TRUE          */
/*              x0;     polynomial will be computed at x0           */
/*int           n;      the highest exponent of p                   */
/*unsigned char flag;   flag==TRUE => compute P'(x0)                */
{
     int      i;     /* counter                                     */
     dcomplex help1; /* help variable                               */ 

     *f  = p[n];
     if (flag==TRUE) {              /* if flag=TRUE, compute P(x0)  */
          COMPLEXM(*df,0.,0.);      /* and P'(x0)                   */
          for (i=n-1; i>=0; i--) {
               CMUL(help1,*df,x0);  /* *df = *f   + *df * x0        */
               CADD(*df,help1,*f);
               CMUL(help1,*f,x0);   /* *f  = p[i] + *f * x0         */
               CADD(*f,help1,p[i]);
          }
     } else                       /* otherwise: compute only P(x0)  */
          for (i=n-1; i>=0; i--) {
               CMUL(help1,*f,x0);   /* *f = p[i] + *f * x0          */
               CADD(*f,help1,p[i]);
          }
}