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

Go to the documentation of this file.

/*
 * lib/numerics/rootsolve/sturm.cc 
 * 
 * Numerics library root solving: Polynomial root solver based on the 
 * Sturm sequence. 
 * 
 * Copyright (c) 2004 by Wolfgang Wieser ] wwieser (a) gmx <*> de [
 * 
 * This file is heavily based upon SturmSequences.c,
 * Using Sturm Sequences to Bracket Real Roots of Polynomial Equations
 * by D.G. Hook and P.R. McAree from "Graphics Gems", Academic Press, 1990
 * 
 * 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_Sturm::Params PolyRootSolver_Sturm::defaults = 
{
    // I am using a default of 1e-10 in simple.cc solvers so don't be more 
    // exact here. 
    INIT_FIELD(nearly_zero) 1.0e-10,  /*1e-12 */
    INIT_FIELD(relerror) 1.0e-14,
    INIT_FIELD(sbisect_maxit) 64,
    INIT_FIELD(modrf_maxit) 64,
};


/*
 * modrf
 *
 * Uses the modified regula-falsi method to evaluate the root
 * in interval [a,b] of the polynomial described in coef. The
 * root is returned is returned in *val. The routine returns zero
 * if it can't converge.
 */
int PolyRootSolver_Sturm::modrf(int ord,dbl *coef,dbl a,dbl b,
    dbl *val)
{
    // Evaluate polynomial at a and b: 
    dbl *scoef = coef;
    dbl *ecoef = &coef[ord];
    dbl fa,fb = fa = *ecoef;
    for(dbl *fp = ecoef - 1; fp >= scoef; fp--)
    {
        fa = a * fa + *fp;
        fb = b * fb + *fp;
    }
    
    // if there is no sign difference the method won't work
    if(fa * fb > 0.0)  return(0);
    
    if(fabs(fa) < par.relerror)  {  *val = a;  return(1);  }
    if(fabs(fb) < par.relerror)  {  *val = b;  return(1);  }
    
    dbl lfx = fa;
    dbl fx;
    for(int its = 0; its<par.modrf_maxit; its++)
    {
        dbl x = (fb*a - fa*b) / (fb-fa);
        // Evaluate at x: (x,order,coef)
        fx = *ecoef;
        for(dbl *fp = ecoef-1; fp >= scoef; fp--)
            fx = x*fx + *fp;
        
        if( (fabs(x) > par.relerror && fabs(fx/x) < par.relerror) || 
            fabs(fx) < par.relerror )
        {
            if(fabs(fx/x) < par.relerror)
            {  *val=x;  return(1);  }
        }
        
        if(fa*fx < 0)
        {
            b = x;
            fb = fx;
            if((lfx * fx) > 0)  fa *= 0.5;
        }
        else
        {
            a = x;
            fa = fx;
            if((lfx * fx) > 0)  fb *= 0.5;
        }
        
        if(fabs(b-a) < par.relerror)
        {  *val=x;  return(1);  }
        
        lfx = fx;
    }
    
    // This is actually no problem yet since we can still use the 
    // sbisect core code. 
    //fprintf(stderr, "modrf overflow %f %f %f\n", a, b, fx);
    return(0);
}
    

/*
 * modp
 *
 * Calculates the modulus of u(x) / v(x) leaving it in r, it
 * returns 0 if r(x) is a constant.
 * Note: this function assumes the leading coefficient of v 
 *       is 1 or -1
 */
int PolyRootSolver_Sturm::modp(poly *u,poly *v,poly *r)
{
#if !PolyRootSolver_Sturm_BufOnStack
    // NOTE: The number of coefficiants in r->coef[] used during 
    //       computation is u->ord+1 (i.e. order u->ord)!! 
    // (and NOT v->ord or dito -1, even if it seems so from below!! (Wolfgang)
    r->AllocCoeffs(u->ord);
#endif
        
    dbl *nr = r->coef;
    dbl *end = &u->coef[u->ord];
    dbl *uc = u->coef;
    while (uc <= end)
        *nr++ = *uc++;
    
    if(v->coef[v->ord] < 0.0)
    {
        for(int k = u->ord - v->ord - 1; k >= 0; k -= 2)
            r->coef[k] = -r->coef[k];
        
        for(int k = u->ord - v->ord; k >= 0; k--)
            for(int j = v->ord + k - 1; j >= k; j--)
                r->coef[j] = -r->coef[j] - r->coef[v->ord+k] * v->coef[j-k];
    }
    else
    {
        for(int k = u->ord - v->ord; k >= 0; k--)
            for(int j = v->ord + k - 1; j >= k; j--)
                r->coef[j] -= r->coef[v->ord+k] * v->coef[j-k];
    }
    
    int k = v->ord - 1;
    while(k >= 0 && fabs(r->coef[k]) < par.nearly_zero)  r->coef[k--]=0.0;
    
    r->ord = (k<0) ? 0 : k;
    
    return(r->ord);
}


/*
 * buildsturm
 *
 * Build up a sturm sequence for a polynomial in smat, returning
 * the number of polynomials in the sequence
 */
int PolyRootSolver_Sturm::buildsturm(int ord,poly *sseq)
{
    sseq[1].ord = ord-1;
#if !PolyRootSolver_Sturm_BufOnStack
    sseq[1].AllocCoeffs(ord-1);
#endif

    // calculate the derivative and normalise the leading coefficient.
    dbl f = fabs(sseq[0].coef[ord] * ord);
    dbl *fp = sseq[1].coef;
    dbl *fc = sseq[0].coef+1;
    for(int i=1; i<=ord; i++)
        *fp++ = *fc++ * i / f;

    // construct the rest of the Sturm sequence
    poly *sp;
    for(sp = sseq+2; modp(sp-2,sp-1,sp); sp++)
    {
        // reverse the sign and normalise
        f = -fabs(sp->coef[sp->ord]);
        for(fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)  *fp /= f;
    }
    
    sp->coef[0] = -sp->coef[0];  // reverse the sign
    
    return(sp-sseq);
}

/*
 * numroots
 *
 * Return the number of distinct real roots of the polynomial
 * described in sseq.
 */
int PolyRootSolver_Sturm::numroots(int np,poly *sseq,int *atneg,int *atpos)
{
    poly *s;
    
    // changes at positve infinity
    int atposinf=0;
    dbl lf = sseq[0].coef[sseq[0].ord];
    #if 0
        // Original code; see numchanges() for why this was replaced. 
        for(s = sseq+1; s <= sseq+np; s++)
        {
            dbl f = s->coef[s->ord];
            if(lf==0.0 || lf * f < 0)  atposinf++;
            lf=f;
        }
    #else
        for(s=sseq+1; lf==0.0 && s <= sseq+np; s++)  lf = s->coef[s->ord];
        for(; s<=sseq+np; s++)
        {
            dbl f = s->coef[s->ord];
            if(f==0)  continue;  // skip zeros
            if(f*lf<0)  ++atposinf;
            lf=f;
        }
    #endif
    *atpos = atposinf;

    // changes at negative infinity
    int atneginf=0;
    if(sseq[0].ord & 1)  lf = -sseq[0].coef[sseq[0].ord];
    else                 lf =  sseq[0].coef[sseq[0].ord];
    #if 0
        // Original code; see numchanges() for why this was replaced. 
        for(s = sseq+1; s <= sseq+np; s++)
        {
            dbl f = (s->ord & 1) ? -s->coef[s->ord] : s->coef[s->ord];
            if(lf==0.0 || lf * f < 0)  atneginf++;
            lf=f;
        }
    #else
        for(s=sseq+1; lf==0.0 && s <= sseq+np; s++)
            lf = (s->ord & 1) ? -s->coef[s->ord] : s->coef[s->ord];
        for(; s<=sseq+np; s++)
        {
            dbl f = (s->ord & 1) ? -s->coef[s->ord] : s->coef[s->ord];
            if(f==0)  continue;  // skip zeros
            if(f*lf<0)  ++atneginf;
            lf=f;
        }
    #endif
    *atneg = atneginf;
    
    #if 0
    // Sometimes, one may know the "left" border of interest is some 
    // fixed value. By substituting x one can then make the left border 
    // be 0. (In case you know the "right" and not the "left" border, 
    // simply mirror the polynomial or choose between posinf and neginf 
    // for the other border) 
    // Here is the code in case it is needed. (Just do not forget to 
    // remove the neginf calculation (or posinf for "right" border 
    // and update the return value).)
    
    // changes at 0
    int at0=0;
    lf = sseq[0].coef[0];
    #if 0
        for(poly *s = sseq+1; s <= sseq+np; s++)
        {
            dbl f = s->coef[0];
            if(lf==0.0 || lf * f < 0)  ++at0;
            lf=f;
        }
    #else
        for(s=sseq+1; lf==0.0 && s <= sseq+np; s++)  lf = s->coef[0];
        for(; s <= sseq+np; s++)
        {
            dbl f = s->coef[0];
            if(f==0)  continue;  // skip zeros
            if(f*lf<0)  ++at0;
            lf=f;
        }
    #endif
    *at00=at0;
    #endif
    
    return(atneginf-atposinf);
}


/*
 * numchanges
 *
 * Return the number of sign changes in the Sturm sequence in
 * sseq at the value a.
 */
int PolyRootSolver_Sturm::numchanges(int np,poly *sseq,dbl a)
{
    #if 0
        #warning "Buggy version enabled."
        
        // This is the original implementation. 
        int changes=0;
        dbl lf = PolvEval(a,sseq[0].coef,sseq[0].ord);
        for(poly *s = sseq+1; s<=sseq+np; s++)
        {
            dbl f = PolvEval(a,s->coef,s->ord);
            if(lf == 0.0 || lf*f < 0)  changes++;
            lf = f;
        }
    #else
        // It seems this has a bug since Bronstein, et al state 
        // that zero values are NOT to be counted as changes. 
        // Hence, I suggest the following: (Wolfgang)
        int changes=0;
        poly *s = sseq;
        dbl lf;  // <-- Ignore warning about uninitialized var "lf". 
        // Find first value which is not equal to 0: 
        for(; s<=sseq+np; s++)
        {
            lf = PolvEval(a,s->coef,s->ord);
            if(lf!=0.0)  break;
        }
        // Count number of sign changes: 
        for(; s<=sseq+np; s++)
        {
            dbl f = PolvEval(a,s->coef,s->ord);
            if(f==0.0)  continue;  // skip zeros
            if(f*lf<0)  ++changes;
            lf=f;
        }
    #endif
    
    return(changes);
}


/*
 * sbisect
 *
 * Uses a bisection based on the sturm sequence for the polynomial
 * described in sseq to isolate intervals in which roots occur,
 * the roots are returned in the roots array in order of magnitude.
 */
int PolyRootSolver_Sturm::sbisect(int np,poly *sseq,dbl min,dbl max,
    int atmin,int atmax,dbl *roots)
{
    if(atmin-atmax == 1)
    {
        // First try a less expensive technique (regula flasa): 
        if(modrf(sseq->ord,sseq->coef,min,max,&roots[0]))  return(1);
        
        // If we get here we have to evaluate the root the hard
        // way by using the Sturm sequence. 
        dbl mid;
        for(int its=0; its<par.sbisect_maxit; its++)
        {
            // Test in the middle; this will cut the interval in half for 
            // every iteration; number of iterations are guaranteed to 
            // be O(log(interval_size/accuracy)). 
            // NOTE: Using linear interpolation here is no good idea 
            //       (verified!). 
            mid = 0.5*(min+max);
            
            if( (fabs(mid) > par.relerror && fabs((max-min)/mid) < par.relerror) || 
                fabs(max-min) < par.relerror )
            {
                // Additional linear interpolation step at the end. (Wolfgang)
                // (Instead of: roots[0] = mid)
                dbl lv = PolvEval(min,sseq[0].coef,sseq[0].ord);
                dbl rv = PolvEval(max,sseq[0].coef,sseq[0].ord);
                roots[0] = min + (max-min)*lv/(lv-rv);
                return(1);
            }
            
            int atmid = numchanges(np,sseq,mid);
            
            // Note that we're only here for atmin-atmax=1. 
            Assert(atmid>=atmax && atmid<=atmin);
            // In case this assertion fails [and it is not due to a bug 
            // in numchanges()] we should probably return(0). 
            
            if(atmin==atmid)  min=mid;
            else              max=mid;
        }
        
        fprintf(stderr,"sbisect: overflow min %g max %g diff %g nroot %d\n",
            min,max,max-min,atmin-atmax);
        // This is the best guess we have: 
        // Additional linear interpolation step at the end. (Wolfgang)
        // (Instead of: roots[0] = mid)
        dbl lv = PolvEval(min,sseq[0].coef,sseq[0].ord);
        dbl rv = PolvEval(max,sseq[0].coef,sseq[0].ord);
        roots[0] = (lv*rv<0) ? ( min + (max-min)*lv/(lv-rv) ) : mid;
        
        return(1);
    }

    // more than one root in the interval, we have to bisect...
    int n1,n2;
    dbl mid;
    for(int its=0; its<par.sbisect_maxit; its++)
    {
        mid = 0.5*(min+max);
        if( (fabs(mid) > par.relerror && fabs((max-min)/mid) < par.relerror) || 
            fabs(max-min) < par.relerror )
        {
            // Additional linear interpolation step at the end. (Wolfgang)
            // (Instead of: roots[0] = mid)
            dbl lv = PolvEval(min,sseq[0].coef,sseq[0].ord);
            dbl rv = PolvEval(max,sseq[0].coef,sseq[0].ord);
            roots[0] = min + (max-min)*lv/(lv-rv);
            return(1);
        }
        
        int atmid = numchanges(np,sseq,mid);
        Assert(atmid>=atmax && atmid<=atmin);  // See same CritAssert() above. 
        
        // Number or roots in both intervals which should be there: 
        n1 = atmin-atmid;
        n2 = atmid-atmax;
        
        if(n1 && n2)
        {
            // n1, n2: Number of roots which were actually found: 
            n1=sbisect(np,sseq,min,mid,atmin,atmid,roots);
            n2=sbisect(np,sseq,mid,max,atmid,atmax,&roots[n1]);
            return(n1+n2);
        }
        
        if(!n1)  min=mid;
        else     max=mid;
    }
    
    fprintf(stderr, "sbisect: roots too close together\n");
    fprintf(stderr, 
        "sbisect: overflow min %g max %g diff %g nroot %d n1 %d n2 %d\n",
        min, max, max - min, atmin-atmax, n1, n2);
    
    // Additional linear interpolation step at the end. (Wolfgang)
    // (Instead of: roots[0] = mid)
    dbl lv = PolvEval(min,sseq[0].coef,sseq[0].ord);
    dbl rv = PolvEval(max,sseq[0].coef,sseq[0].ord);
    if(lv*rv<0)
    {  mid = min + (max-min)*lv/(lv-rv);  }
    
    #if 0
        // Fill up the roots with the best guess we have; this may be 
        // pretty stupid but at least we get the number of roots right. 
        int n2=atmin-atmax;
        for(n1=0; n1<n2; n1++)
            roots[n1] = mid;
        return(n2);
    #else
        // We could also just use one root which works since sbisect() 
        // now returns (properly) the number of roots found (Wolfgang). 
        roots[0]=mid;
        return(1);
    #endif
}


int PolyRootSolver_Sturm::solve(int order,const dbl *c,dbl *r)
{
    
    poly sseq[order+1];
#if PolyRootSolver_Sturm_BufOnStack
    dbl buf_on_stack[(order+1)*(order+1)];
    for(int i=0; i<=order; i++)
        sseq[i].coef = &buf_on_stack[i*(order+1)];
#else
    sseq[0].AllocCoeffs(order);
#endif
    
    // Set up the coefficients: 
    sseq[0].ord=order;
    for(int i=0; i<=order; i++)
        sseq[0].coef[i]=c[i];
    
    // build the Sturm sequence
    int np = buildsturm(order, sseq);
    Assert(np<=order);
    
    // get the number of real roots
    int atmin,atmax;
    int nroots = numroots(np,sseq,&atmin,&atmax);
    if(!nroots)  return(0);
    
    // calculate the bracket that the roots live in
    dbl min,max;
//  -> specify left border as 0 or -\infty (?), ...
dbl bracket_min=10,bracket_max=-10;
#define     MAXPOW      32  /* max power of 10 we wish to search to */

    
    if(bracket_min>bracket_max)
    {
        min = -1.0;
        int nchanges = numchanges(np,sseq,min);
        for(int i=0; nchanges!=atmin && i!=MAXPOW; i++)
        { 
            min *= 10.0;
            nchanges = numchanges(np,sseq,min);
        }
        if(nchanges != atmin)
        {
            fprintf(stderr,"sturm solve: unable to bracket all negetive roots\n");
            atmin = nchanges;
        }
        
        max = 1.0;
        nchanges = numchanges(np, sseq, max);
        for(int i=0; nchanges!=atmax && i!=MAXPOW; i++)
        {
            max *= 10.0;
            nchanges = numchanges(np,sseq,max);
        }
        if(nchanges != atmax)
        {
            fprintf(stderr,"sturm solve: unable to bracket all positive roots\n");
            atmax = nchanges;
        }
    }
    else
    {
        min=bracket_min;  // e.g. 0 if we need only positive roots
        max=bracket_max;
        
        atmin = numchanges(np,sseq,min);
        atmax = numchanges(np,sseq,max);
    }
    
    // Expected number of roots: 
    nroots = atmin-atmax;
    if(!nroots)  return(0);
    
    // perform the bisection.
    nroots=sbisect(np,sseq,min,max,atmin,atmax,r);
    return(nroots);
}


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

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

}  // end of namespace NUM


#if 0
// Test routine: 
/*
 make ADDFLAGS="-O0 -g -fno-omit-frame-pointer"
 gcc -o sturm -lm `find ../../.. -name "*.a"` `find ../../.. -name "*.a"` -pthread
*/

int main()
{
    // This is nothing intelligent; just throw random havoc onto the 
    // routines. 
    const int maxorder=11;
    dbl c[maxorder+1],r[maxorder];
    int cnt=0;
    
    NUM::PolyRootSolver_Sturm solver;
    
    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-3)  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);
    return(0);
}
#endif

---