/ray/src/lib/numerics/rootsolve/bairstow.cc

Go to the documentation of this file.

/*
 * lib/numerics/rootsolve/bairstow.cc 
 * 
 * Numerics library root solving: Polynomial root solver based on the 
 * Bairstow method and some other tricks. 
 * 
 * Copyright (c) 2004 by Wolfgang Wieser ] wwieser (a) gmx <*> de [
 * This file is heavily based upon code published by C. Bond, (c) 2002,2003. 
 * Thanks to him for making it available! 
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 2 as published by the Free Software
 * Foundation. (See COPYING.GPL for details.)
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 * 
 */ 

#include "rootsolve.h"

#include <stdio.h>  /* FIXME!! */


namespace NUM
{

// Initialize static const data: 
// Values smaller than this are treated as zero. 
const PolyRootSolver_Bairstow::Params PolyRootSolver_Bairstow::defaults = 
{
    INIT_FIELD(nearly_zero) 1e-10,
    INIT_FIELD(maxiter) 500,
    INIT_FIELD(iter_chg) 200,
    INIT_FIELD(initial_epsilon) 1e-10,  // originally 1e-15
    INIT_FIELD(disc_epsilon) 1e-10
};


// Extract individual real roots from list of quadratic factors. 
int PolyRootSolver_Bairstow::_QuadricRealRoots(dbl *a,int n,dbl *wr)
{
    int m=n;
    int numroots=0;
    while(m>1)
    {
        dbl b2 = -0.5*a[m-2];
        dbl c = a[m-1];
        
        // Comute the discriminant: 
        dbl disc = b2*b2-c;
        if(fabs(disc) <= par.disc_epsilon * (b2*b2+fabs(c)))
        {
            // One dbl real root. 
            wr[numroots++] = b2;
        }
        else if(disc>0.0)
        {
            // Two real roots. 
            disc = sqrt(disc);
            disc = b2<0 ? (b2-disc) : (b2+disc);
            wr[numroots++] = disc;
            wr[numroots++] = c/disc;
        }
        
        m -= 2;
    }
    
    if(m==1)
    {  wr[numroots++] = -a[0];  }
    
    return(numroots);
}


// Extract individual real or complex roots from list of quadratic factors. 
int PolyRootSolver_Bairstow::_QuadricRoots(dbl *a,int n,dbl *wr,
    dbl *wi)
{
    dbl sq,b2,c,disc;
    int m,numroots;
    
    // This routine should not be used as such. Should use numroots 
    // as counter for wr[] and wi[] indices and not store zero value 
    // pairs for "no root" indication. See _QuadricRealRoots()...
    CritAssert(0);  // not used / unimplemented
    
    m = n;
    numroots = 0;
    while (m > 1) {
        b2 = -0.5*a[m-2];
        c = a[m-1];
        // Comute the discriminant: 
        disc = b2*b2-c;
        if (fabs(disc) <= par.disc_epsilon*(b2*b2+fabs(c))) disc = 0.0;
        if (disc < 0.0)
        {
            // Two conjugated complex roots. 
            sq = sqrt(-disc);
            wr[m-2] = b2;
            wi[m-2] = sq;
            wr[m-1] = b2;
            wi[m-1] = -sq;
            numroots+=2;
        }
        else
        {
            sq = sqrt(disc);
            
            // Original code: 
            //  wr[m-2] = fabs(b2)+sq;
            //  if (b2 < 0.0) wr[m-2] = -wr[m-2];
            // This should be slightly more efficient (Wolfgang): 
            wr[m-2] = b2<0 ? (b2-sq) : (b2+sq);
            
            if(wr[m-2] == 0)
            {
                // Little bugfix: original code did "wr[m-1] == 0;"
                // (which is a no-op as statement). (Wolfgang)
                // This case should not happen since the polynomials 
                // may not have roots at x=0 for this routine anyways. 
                
                // So, use wr[m-2]=0 (already) and wr[m-1] = 0 and 
                // have them skipped by the output routine. 
                wr[m-1] = 0;
            }
            else
            {
                wr[m-1] = c/wr[m-2];
                numroots+=2;
            }
            
            // Real roots have zero imaginary part: 
            wi[m-2] = 0.0;
            wi[m-1] = 0.0;
        }
        m -= 2;
    }
    if (m == 1) {
       wr[0] = -a[0];
       wi[0] = 0.0;
       numroots++;
    }
    return numroots;
}


// Find a quadratic factor applying Bairstow's (= quadratic Newton) method. 
// Returns number of iterations needed and error estimate in *err_est. 
int PolyRootSolver_Bairstow::_FindQuadFact(dbl *a,int n,dbl *b,
    dbl *quad,int red_fact,dbl *err_est)
{
    // This routine had a severe problem when optimization is switched on 
    // using the GCC compiler on the i386 architecture and there are 
    // multiple roots such as in x^4-9x^3+29x^2-39x+18. This may be due 
    // to registers being more precise (80 bit) than normal dbl (64 bit). 
    // Well... hmmm... I played around a bit, converted the while() loop 
    // into a do..while() loop, re-arranged some statements and the 
    // problem magically disappeared. 
    // Oh-oh... there is a sleeping dragon inside this code! (Wolfgang)
    
    //fprintf(stderr,"FQ:");
    //for(int j=0; j<=n; j++)
    //  fprintf(stderr," %g,%g",a[j],b[j]);

    dbl c[n+1];
    c[0] = 1.0;
    dbl r = quad[1];
    dbl s = quad[0];
    dbl dr = 1.0;
    dbl ds = 0;
    dbl eps = par.initial_epsilon;
    
    int iter=1;
    
    do
    {
        if(iter>par.maxiter)  break;
        
        // Increase epsilon in case of convergence problems...
        if(!(iter % par.iter_chg))
        {  eps*=10.0;  }
        
        b[1] = a[1] - r;
        c[1] = b[1] - r;
        
        for(int i=2; i<=n; i++)
        {
            b[i] = a[i] - r * b[i-1] - s * b[i-2];
            c[i] = b[i] - r * c[i-1] - s * c[i-2];
        }
        
        dbl drn=b[n] * c[n-3] - b[n-1] * c[n-2];
        dbl dsn=b[n-1] * c[n-1] - b[n] * c[n-2];
        
        dbl dn=c[n-1] * c[n-3] - c[n-2] * c[n-2];
        // Work around zero slope at iteration point: 
        if(fabs(dn) < 1e-10)
        {  dn = dn<0.0 ? -1e-8 : 1e-8;  }
        
        // Some sort of underrelaxation in case it will not converge 
        // otherwise. (Wolfgang)
        if(red_fact)
        {
            switch(red_fact)
            {
                case 1:  dn*=0.999;  break;
                case 2:  dn*=1.1;  break;
                default: dn*=0.8;  break;
            }
        }
        
        dn = 1.0/dn;
        dr = drn * dn;
        ds = dsn * dn;
        r += dr;
        s += ds;
        
        ++iter;
        
        //fprintf(stderr,"d=%g  (%g,%g)  eps=%g\n",1.0/dn,dr,ds,eps);
    }
    while((fabs(dr)+fabs(ds)) > eps);
    
    quad[0] = s;
    quad[1] = r;
    *err_est = fabs(ds)+fabs(dr);
    
    //fprintf(stderr," -> %g,%g %g %d\n",quad[0],quad[1],*err_est,iter);
    
    return(iter);
}


// Top level routine to manage the determination of all roots of the given
// polynomial 'a', returning the quadratic factors (and possibly one linear
// factor) in 'x'.
void PolyRootSolver_Bairstow::_GetQuads(dbl *a,int n,dbl *quad,dbl *x)
{
    dbl tmp;
    if((tmp = a[0]) != 1.0)
    {
        a[0] = 1.0;
        for(int i=1; i<=n; i++)
            a[i] /= tmp;
    }
    if(n == 2)
    {
        x[0] = a[1];
        x[1] = a[2];
        return;
    }
    if(n == 1)
    {
        x[0] = a[1];
        return;
    }
    
    dbl b[n+1];
    dbl z[n+1];
    b[0] = 1.0;
    for(int i=0; i<=n; i++)
    {
        z[i] = a[i];
        x[i] = 0.0;
    }
    
    dbl err;
    int m = n;
    do 
    {
        int red_fact=0;
        
        if(n>m)
        {
            quad[0] = 3.14159e-1;
            quad[1] = 2.78127e-1;
        }
loop:
        int iter=_FindQuadFact(z,m,b,quad,red_fact,&err);
        if( err>1e-7 || iter>par.maxiter )
        {
            _DiffPoly(z,m,b);
            iter = _Recurse(z,m,b,m-1,quad,&err);
        }
        err=_Deflate(z,m,b,quad);
        if(err > 0.01)
        {
            fprintf(stderr,"Bairstow: Excessive error %g (iter=%d, n=%d) %s\n",
                err,iter,n,red_fact<3 && iter ? "[re-try]" : "[giving up]");
            if(red_fact<3 && iter)
            {
                // Try again with "under-relaxation" (Wolfgang). 
                ++red_fact;
                goto loop;
            }
            
            // What should we do...? Simply go on or store quad[] ?
            // No because this would end up in a potential endless loop. 
            // Hack: Set x values for imaginary root (skipped by 
            //       real quad solver). 
            // (Note: The ONE illegal solutioo which still remains is 
            //        probably the one at the end of this function: 
            //        "else x[0] = b[1];". 
            x[m-2] = 0;
            x[m-1] = 2;
            goto go_on;
        }
        if(err>1)  goto loop;
        x[m-2] = quad[1];
        x[m-1] = quad[0];
go_on:
        m -= 2;
        for(int i=0; i<=m; i++)
            z[i] = b[i];
    }
    while(m>2);
    
    if(m==2)
    {
        x[0] = b[1];
        x[1] = b[2];
    }
    else x[0] = b[1];
}


// Attempt to find a reliable estimate of a quadratic factor using modified
// Bairstow's method with provisions for 'digging out' factors associated
// with multiple roots.
//
// This resursive routine operates on the principal that differentiation of
// a polynomial reduces the order of all multiple roots by one, and has no
// other roots in common with it. If a root of the differentiated polynomial
// is a root of the original polynomial, there must be multiple roots at
// that location. The differentiated polynomial, however, has lower order
// and is easier to solve.
//
// When the original polynomial exhibits convergence problems in the
// neighborhood of some potential root, a best guess is obtained and tried
// on the differentiated polynomial. The new best guess is applied
// recursively on continually differentiated polynomials until failure
// occurs. At this point, the previous polynomial is accepted as that with
// the least number of roots at this location, and its estimate is
// accepted as the root.
// 
// Returns iteration count. 
int PolyRootSolver_Bairstow::_Recurse(dbl *a,int n,dbl *b,int m,
    dbl *quad,dbl *err)
{
    //              vvv--- originally 1e-16
    if(fabs(b[m]) < par.nearly_zero) --m;  // this bypasses roots at zero
    
    if(m == 2)
    {
        quad[0] = b[2];
        quad[1] = b[1];
        *err = 0;
        return(0);
    }
    
    dbl c[m+1];
    dbl x[n+1];
    c[0] = x[0] = 1.0;
    dbl rs[2];
    rs[0] = quad[0];
    rs[1] = quad[1];
    int iter=_FindQuadFact(b,m,c,rs,0,err);
    dbl tst = fabs(rs[0]-quad[0]) + fabs(rs[1]-quad[1]);
    if(*err < 1e-12)
    {
        quad[0] = rs[0];
        quad[1] = rs[1];
    }
// tst will be 'large' if we converge to wrong root
    if( (iter>5 && tst<1e-4) || (iter>20 && tst<1e-1) )
    {
        _DiffPoly(b,m,c);
        iter=_Recurse(a,n,c,m-1,rs,err);
        quad[0] = rs[0];
        quad[1] = rs[1];
    }
    
    return(iter);
}


// Differentiate polynomial 'a' returning result in 'b'. 
void PolyRootSolver_Bairstow::_DiffPoly(dbl *a,int n,dbl *b)
{
    dbl coef=n;  // divide by n so that leading coefficient stays 1. 
    b[0] = 1.0;
    for(int i=1; i<n; i++)
        b[i] = a[i]*(n-i)/coef;
}


// Deflate polynomial 'a' by division of 'quad'. Return quotient
// polynomial in 'b' and error (metric based on remainder). 
dbl PolyRootSolver_Bairstow::_Deflate(dbl *a,int n,dbl *b,
    dbl *quad)
{
    dbl c[n+1];
    
    dbl r = quad[1];
    dbl s = quad[0];
    
    b[1] = a[1] - r;
    c[1] = b[1] - r;
    
    for(int i=2; i<=n; i++)
    {
        b[i] = a[i] - r * b[i-1] - s * b[i-2];
        c[i] = b[i] - r * c[i-1] - s * c[i-2];
    }
    
    return( fabs(b[n])+fabs(b[n-1]) );
}


// Arguments as usual (e.g. see sturm.cc). 
// Returns number of roots found. 
int PolyRootSolver_Bairstow::solve(int n,const dbl *_c,dbl *r)
{
    int nroots=0;
    
    // This routine cannot deal with highest and lowest coefficiants 
    // being 0, so we need to deal with that here. 
    // First see the high ones: 
    while(n>=0 && fabs(_c[n])<par.nearly_zero)  --n;
    if(!n)  return(nroots);
    // Then the low ones: 
    if(fabs(_c[0])<par.nearly_zero)
    {
        // We have the trivial solution x=0. 
        r[nroots++]=0;
        
        // Divide away all zero solutions: 
        int i=1;
        while(i<=n && fabs(_c[i])<par.nearly_zero)  ++i;
        if(i>=n)  return(nroots);
        _c+=i;
        n-=i;
    }
    
    // These routines store the coefficients the other way round, 
    // so I will copy them: 
    dbl c[n+1];
    for(int i=0; i<=n; i++)  c[i]=_c[n-i];
    
    // initialize estimate for 1st root pair 
    dbl quad[2];
    quad[0] = 2.71828e-1;  // e/10
    quad[1] = 3.14159e-1;  // pi/10
    
    // Get the quadratic factors: 
    dbl x[n+1];
    _GetQuads(c,n,quad,x);
    //for(int i=0; i<n; i++)
    //{  fprintf(stderr,"%g ",x[i]);  }
    //fprintf(stderr,"\n");
    
    // Compute the (real) roots from the quadratic factors. 
    nroots+=_QuadricRealRoots(x,n,&r[nroots]);
    
    return(nroots);
}


PolyRootSolver_Bairstow::PolyRootSolver_Bairstow() : 
    PolyRootSolver(),
    par(defaults)
{
    // nothing to do
}

PolyRootSolver_Bairstow::~PolyRootSolver_Bairstow()
{
    // nothing to do
}

}  // end of namespace NUM


#if 0
// Test routine: 
int main()
{
    NUM::PolyRootSolver_Bairstow solver;
    
    {
        dbl c[10],r[9];
        int deg;
        deg=4; c[4]=1; c[3]=-9; c[2]=29; c[1]=-39; c[0]=18;
        
        int n=solver.solve(deg,c,r);
        fprintf(stderr,"%d roots:",n);
        for(int i=0; i<n; i++)
        {  fprintf(stderr," %.15g",r[i]);  }
        fprintf(stderr,"\n");
    }
    
    {
        dbl c[10],r[9];
        int deg;
        deg=3; c[3]=3.96; c[2]=2.69; c[1]=-5.78; c[0]=2.24;
        
        int n=solver.solve(deg,c,r);
        fprintf(stderr,"%d roots:",n);
        for(int i=0; i<n; i++)
        {  fprintf(stderr," %.15g",r[i]);  }
        fprintf(stderr,"\n");
    }
    
    // This is nothing intelligent; just throw random havoc onto the 
    // routines. 
    #if 1
    const int maxorder=11;
    dbl c[maxorder+1],r[maxorder];
    int cnt=0;
    for(int i=0; i<30000; i++)
    {
        int deg=rand()%(maxorder-3)+3;
        for(int j=0; j<=deg; j++)
            c[j]=(rand()%2000-1000)/100.0;
        int n=solver.solve(deg,c,r);
        for(int j=0; j<n; j++)
        {
            dbl R=NUM::PolvEval(r[j],c,deg);
            if(fabs(R)<1e-6)  continue;
            
            fprintf(stderr,"[%3d]",i);
            for(int k=0; k<=deg; k++)  fprintf(stderr," %g",c[k]);
            fprintf(stderr,"-> %g for %g (n=%d, j=%d)\n",R,r[j],n,j);
            ++cnt;
        }
    }
    fprintf(stderr,"cnt=%d\n",cnt);
    #endif
    
    return(0);
}
#endif

---