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

Go to the documentation of this file.

/*
 * lib/numerics/rootsolve/direct.cc 
 * 
 * Numerics library root solving: Directly solvable cases (i.e. where an 
 * explicit formula exitst). 
 * 
 * Copyright (c) 2004 by Wolfgang Wieser ] wwieser (a) gmx <*> de [
 * 
 * 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"


// THREAD SAFETY CHECK: PASSED. 

namespace NUM
{

// Initialize static const data: 
// Values smaller than this are treated as zero. 
const PolyRootSolver_Direct::Params PolyRootSolver_Direct::defaults = 
{
    INIT_FIELD(nearly_zero) 1e-10  // Rayshade supposes 1e-9
};


// Solve roots of polynomial of degree 1, i.e. 
//  0 = c[0] + c[1]*x
// Root returned in r[0]; number of roots returned as return value. 
int PolyRootSolver_Direct::Solve_1(const dbl *c,dbl *r)
{
    // Well, THIS is easy. 
    
    if(c[1]==0.0)  return(0);
    
    r[0]=-c[0]/c[1];
    
    return(1);
}


// Solve roots of polynomial of degree 2, i.e. 
//  0 = c[0] + c[1]*x + c[2]*x^2
// Roots returned in r[0..1]; number of roots returned as return value. 
int PolyRootSolver_Direct::Solve_2(const dbl *c,dbl *r)
{
    if(c[2]==0.0)
    {
        // In this case the polynomial is actually of degree 1. 
        return(Solve_1(c,r));
    }
    
    #if 0
        //---- First method: use normal form x^2 + p x + q = 0 ----
        // This can be found in any better maths formula book 
        // (e.g. Bronstein et al). 
        dbl p = -0.5*c[1] / c[2];
        dbl q =      c[0] / c[2];
        dbl D = sqr(p) - q;
        // If the discriminant is nearly zero, we have one (dbl) root: 
        if(fabs(D)<par.nearly_zero)
        {  r[0]=p;  return(1);  }
        // Otherwise we have 2 unless the discriminant is smaller than zero. 
        if(D>0.0)
        {
            D = sqrt(D);
            r[0]=p+D;
            r[1]=p-D;
            return(2);
        }
    #else
        //---- Second method: Do as I have learned it in school :) ----
        // This method runs somewhat faster on my machine even if I 
        // replace the 2 divisions st the beginning of the upper 
        // variant by one division and 2 multiplications. 
        // Furthermore, I (Wolfgang) implemented some additional 
        // nonstandard code to reduce leading digit cancellation. 
        dbl D = sqr(c[1]) - 4.0*c[0]*c[2];
        // If the discriminant is nearly zero, we have one (dbl) root: 
        if(fabs(D)<par.nearly_zero)
        {  r[0] = -0.5*c[1]/c[2];  return(1);  }
        // Otherwise we have 2 unless the discriminant is smaller than zero. 
        if(D>0.0)
        {
            D = sqrt(D);
            #if 0
                // This is the standard school-book way of calculting it: 
                dbl f = -0.5/c[2];
                r[0] = (c[1]+D)*f;
                r[1] = (c[1]-D)*f;
            #else
                // This variant will reduce the digit cancellation effect. 
                // Disavantage: one additional fabs() and 
                //              one additional division
                // Advantage: one fewer subtraction, one fewer multiplication
                D = -0.5*( c[1]<0 ? (c[1]-D) : (c[1]+D) );
                // Note: The case D==0 cannot happen here since D>0 
                //       and gets added to c[1] in a way which increases 
                //       the absolute value. 
                r[0] = c[0]/D;
                r[1] = D/c[2];
            #endif
            return(2);
        }
    #endif
    
    return(0);
}


// Solve roots of polynomial of degree 3, i.e. 
//  0 = c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3
// Roots returned in r[0..2]; number of roots returned as return value. 
// This uses the algorithm proposed by Cardano. 
int PolyRootSolver_Direct::Solve_3(const dbl *c,dbl *r)
{
    if(c[3]==0.0)
    {
        // In this case the polynomial is actually of degree 2. 
        return(Solve_2(c,r));
    }
    
    // This application of Cardano's formula is based on 
    // - code taken from Raycer 1.0 BETA, Copyright (C) 2002 Tim E. Madsen, 
    //   Reuben Lee, Justin Pahl, distributed under the GPL. 
    // - "Exact solutions of cubic polynomial equations" by 
    //   Stephen R. Schmitt (ported from Zeno-1.2 to C(++))
    // - Bronstein, et al. 
    
    // Normal form: x^3 + Ax^2 + Bx + C = 0. 
    dbl cc = 1/c[3];
    dbl A = c[2]*cc;
    dbl B = c[1]*cc;
    dbl C = c[0]*cc;
    // Substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0
    dbl sq_A = sqr(A);
    // (These was modified by me (Wolfgang) by multiplying/factoring 
    //  out FP constants and variables; p was changed into -par.)
    dbl p = (sq_A-3.0*B) / 9.0;
    dbl q = A * (sq_A-4.5*B) / 27 + 0.5*C;
    dbl cb_p = p*p*p;
    dbl D = sqr(q)-cb_p;
    // Substitution value (see above; moved here by Wolfgang): 
    dbl sub = A/-3.0;
    
    // This case is extremely unlikely (my test did not trigger it and 
    // I am hence not even sure if it works correctly), so I will leave 
    // it away via #if 0...#endif
    #if 1
    if(fabs(D)<par.nearly_zero)
    {
CritAssert(!"CHECK ME!\n");
        
        if(fabs(q)<par.nearly_zero)
        {
            // One triple solution
            r[0]=sub;
            return(1);
        }
        // One single and one dbl solution
        dbl u = cbrt(-q);
        r[0]=2.0*u+sub;
        r[1]=sub-u;
        return(2);
    }
    #endif
    
    if(D>0.0)
    {
        // One real solution. 
        // This case is probably the most likely one. 
        // There are 3 methods to compute this. 
        
        #if 0
            // Fastest but worst method. This is the original Rayshade algo. 
            D = sqrt(D);
            r[0] = cbrt(D-q) - cbrt(D+q) + sub;
        #elif 0
            // This relies on cbrt() to be well implemented but suffers 
            // somewhat from leading digit cancellation in the sqrt()+q part. 
            D = -cbrt( sqrt(q*q-cb_p) + q );
            r[0] = (D + p/D) + sub;
        #else
            // This is essentially the same as above but working around the 
            // leading digit cancellation. See Solve_2() for a 
            // similar trick. This showed to be the best method in my tests. 
            // The improvement here is much more than in the 2nd order case. 
            D = sqrt(q*q-cb_p);
            D = q<0 ? cbrt(D-q) : -cbrt(D+q) ;
            r[0] = (D + p/D) + sub;
        #endif
        return(1);
    }
    
    // Casus irreducibilis: three real solutions. 
    dbl phi = (1.0/3) * acos(q / sqrt(cb_p));
    dbl t = -2*sqrt(p);
    r[0] = t * cos(phi        ) + sub;
    r[1] = t * cos(phi+FPConst<dbl>::mul_pi_2_3) + sub;
    r[2] = t * cos(phi-FPConst<dbl>::mul_pi_2_3) + sub;
    return(3);
}


// Solve roots of polynomial of degree 4, i.e. 
//  0 = c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4
// Roots returned in r[0..3]; number of roots returned as return value. 
// This uses the algorithm proposed by Francois Vieta. 
int PolyRootSolver_Direct::Solve_4(const dbl *c,dbl *_r)
{
    if(c[4]==0.0)
    {
        // In this case the polynomial is actually of degree 2. 
        return(Solve_3(c,_r));
    }
    
    // This is based on the the implementation as found in 
    // Raycer 1.0 BETA, Copyright (C) 2002 Tim E. Madsen, 
    // Reuben Lee, Justin Pahl, distributed under the GPL. 
    // I (Wolfgang) spent several hours in optimizing this algorithm by 
    // pulling factors out and calculating better temporaries when 
    // hand-inlining the needed solvers for 2nd and 3rd order polynomials. 
    
    // Use normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0. 
    dbl f=1.0/c[4];
    dbl A = c[3]*f;
    dbl B = c[2]*f;
    dbl C = c[1]*f;
    dbl D = c[0]*f;

    // Substitute x = y - A/4 to eliminate cubic term: 
    //   x^4 + px^2 + qx + r = 0
    // Note: These constants are all correct (e.g. 3/256 = 0.011718750000..); 
    //       I (Wolfgang) computed them from Raycer-1.0 code. 
    dbl sq_A = A * A;
    dbl p = -0.375 * sq_A + B;
    dbl q = (0.125 * sq_A - 0.5*B)*A + C;
    //dbl r = (-0.01171875*sq_A + 0.0625*B)*sq_A - 0.25*A*C + D;
    dbl r = -0.01171875*sq_A*sq_A + 0.0625*sq_A*B - 0.25*A*C + D;
    
    dbl coeffs[4];
    int num;
    if(fabs(r)<par.nearly_zero)
    {
        // No absolute term: y(y^3 + py + q) = 0 
        coeffs[0] = q;
        coeffs[1] = p;
        coeffs[2] = 0;
        coeffs[3] = 1;
        num = Solve_3(coeffs,_r);
        
        // Resubstitute: 
        dbl sub = -0.25*A;
        for(int i=0; i<num; i++)  _r[i]+=sub;
        
        // Add the solution y=0: 
        _r[num++] = sub;
    }
    else
    {
        // Solve the resolvent cubic...
        #if 0
            coeffs[0] =  0.5*r*p - 0.125*q*q;
            coeffs[1] = -r;
            coeffs[2] = -0.5*p;
            coeffs[3] =  1;
            num = Solve_3(coeffs,_r);
            if(!num)  return(0);  // This should not happen... (3rd order!!)

            // ...and take the one real solution...
            dbl z = _r[0];
        #else
            // This is the same as above but "hand-inlined": 
            // This works much faster than the above version. 
            dbl z;  // "any real solution"
            do {
                dbl sq_p = sqr(p);
                dbl pp = (sq_p+12*r)/36.0;
                dbl qq = -0.0625*q*q - p*(sq_p-36*r)/216.0; // <-- not "pp"
                dbl cb_pp = pp*pp*pp;
                dbl D = sqr(qq)-cb_pp;
                dbl sub = p/6.0;  // <-- not "pp"
                
                if(fabs(D)<par.nearly_zero)
                {
                    if(fabs(qq)<par.nearly_zero)
                    {  z=sub;  break;  }
                    // Take one of these: 
                    //z=2.0*cbrt(-q)+sub;   // single solution
                    z=sub-cbrt(-qq);       // dbl solution
                    break;
                }
                if(D>0.0)
                {
                    D = sqrt(qq*qq-cb_pp);
                    D = qq<0 ? cbrt(D-qq) : -cbrt(D+qq) ;
                    z = (D + pp/D) + sub;
                    break;
                }
                // Take any one of these: 
                z = -2*sqrt(pp) * cos((1.0/3) * acos(qq/sqrt(cb_pp))        ) + sub;
                //z = -2*sqrt(pp) * cos((1.0/3) * acos(qq/sqrt(cb_pp))+M_2PI_3) + sub;
                //z = -2*sqrt(pp) * cos((1.0/3) * acos(qq/sqrt(cb_pp))-M_2PI_3) + sub;
                break;
            } while(0);
        #endif
        
        // ...to build two quadric equations
        dbl u = z*z - r;
        if(fabs(u)<par.nearly_zero)  u=0;
        else if(u>0)  u=sqrt(u);
        else return 0;
        
        dbl v = 2*z - p;
        if(fabs(v)<par.nearly_zero)  v=0;
        else if(v>0)  v=sqrt(v);
        else return 0;
        
        #if 0
            coeffs[0] = z-u;
            coeffs[1] = q<0 ? -v : v;
            coeffs[2] = 1;
            num=Solve_2(coeffs,_r);
            
            coeffs[0] = z+u;
            coeffs[1] = -coeffs[1];  // q<0 ? v : -v;
            //coeffs[2] = 1;   (already)
            num += Solve_2(coeffs,_r+num);
            
            // Resubstitute: 
            dbl sub = -0.25*A;
            for(int i=0; i<num; i++)  _r[i]+=sub;
        #else
            // This is the same as above but "hand-inlined": 
            // This works faster than the above version. 
            num=0;
            dbl c_1 = 0.5 * (q<0 ? -v : v);
            bool c_1_gn = c_1>0;
            dbl DD = sqr(c_1) - z;
            dbl sub = -0.25*A;  // (undo substitution)
            dbl D = DD + u;
            if(fabs(D)<par.nearly_zero)
            {  _r[num++] = sub - c_1;  }
            else if(D>0)
            {
                D = sqrt(D);
                #if 0  /* faster */
                    dbl sub2 = sub - c_1;
                    _r[num++] = sub2 - D;
                    _r[num++] = sub2 + D;
                #else  /* slightly more accurate */
                    D = c_1_gn ? (c_1+D) : (c_1-D);
                    _r[num++] = sub - D;
                    _r[num++] = (u-z)/D + sub;
                #endif
            }
            
            D = DD - u;
            if(fabs(D)<par.nearly_zero)
            {  _r[num++] = sub + c_1;  }
            else if(D>0) 
            {
                D = sqrt(D);
                #if 0  /* faster */
                    dbl sub2 = sub + c_1;
                    _r[num++] = sub2 - D;
                    _r[num++] = sub2 + D;
                #else  /* slightly more accurate */
                    D = c_1_gn ? (c_1+D) : (c_1-D);
                    _r[num++] = D + sub;
                    _r[num++] = (z+u)/D + sub;
                #endif
            }
        #endif
    }
    
    return(num);
}


int PolyRootSolver_Direct::solve(int order,const dbl *c,dbl *r)
{
    switch(order)
    {
        case 1:  return(Solve_1(c,r));
        case 2:  return(Solve_2(c,r));
        case 3:  return(Solve_3(c,r));
        case 4:  return(Solve_4(c,r));
        //default: run off...
    }
    
    return(-2);   // invalid order parameter
}


PolyRootSolver_Direct::PolyRootSolver_Direct() : 
    PolyRootSolver(),
    par(defaults)
{
    // empty
}

PolyRootSolver_Direct::~PolyRootSolver_Direct()
{
    // Not inline because virtual. 
}

}  // end of namespace NUM


#if 0
// Test routine: NOTE: Test routine is for old non-object oriented code; 
//                     needs to be changed for current version. 
int main()
{
    // This is nothing intelligent; just throw random havoc onto the 
    // routines. 
    dbl c[5],r[4],rB[4];
    int deg=2;
    #if 1
    fprintf(stderr,"DEG=%d\n",deg);
    for(int i=0; i<1000000; i++)
    {
        for(int j=0; j<=deg; j++)
            c[j]=(rand()%2000-1000)/100.0;
        int n=NUM::SolvePolyRoots_2(c,r);
        for(int j=0; j<n; j++)
        {
            dbl R=NUM::PolvEval(r[j],c,deg);
            if(fabs(R)>5e-13)  // 13
                fprintf(stderr,"[%3d] %g %g %g -> %g for %g (n=%d, j=%d)\n",i,
                    c[0],c[1],c[2],R,r[j],n,j);
        }
    }
    
    deg=3;
    fprintf(stderr,"DEG=%d\n",deg);
    for(int i=0; i<1000000; i++)
    {
        for(int j=0; j<=deg; j++)
            c[j]=(rand()%2000-1000)/100.0;
        int n=NUM::SolvePolyRoots_3(c,r);
        //int nB=NUM::SolvePolyRoots_3_B(c,rB);
        for(int j=0; j<n; j++)
        {
            dbl R=NUM::PolvEval(r[j],c,deg);
            if(fabs(R)>1e-9)
                fprintf(stderr,"[%6d ] %g %g %g %g -> %g for %g (n=%d,j=%d)\n",
                    i,c[0],c[1],c[2],c[3],R,r[j],n,j);
        }
        /*for(int j=0; j<nB; j++)
        {
            dbl R=NUM::PolvEval(rB[j],c,deg);
            if(fabs(R)>1e-9)
                fprintf(stderr,"[%6dB] %g %g %g %g -> %g for %g (n=%d,j=%d)\n",
                    i,c[0],c[1],c[2],c[3],R,rB[j],nB,j);
        }
        if(n!=nB)
        {
            fprintf(stderr,"[%6d*]: n=%d, nB=%d:",i,n,nB);
            for(int j=0; j<n; j++)  fprintf(stderr," %g",r[j]);
            fprintf(stderr," <->");
            for(int j=0; j<nB; j++)  fprintf(stderr," %g",rB[j]);
            fprintf(stderr,"\n");
        }*/
    }
    #endif
    
    deg=4;
    fprintf(stderr,"DEG=%d\n",deg);
    for(int i=0; i<1000000; i++)
    {
        for(int j=0; j<=deg; j++)
            c[j]=(rand()%2000-1000)/100.0;
        int n=NUM::SolvePolyRoots_4(c,r);
        //int nB=NUM::SolvePolyRoots_4_B(c,rB);
        for(int j=0; j<n; j++)
        {
            dbl R=NUM::PolvEval(r[j],c,deg);
            if(fabs(R)>1e-5)  // 8
                fprintf(stderr,"[%6d ] %g %g %g %g %g -> %g for %g (n=%d,j=%d)\n",
                    i,c[0],c[1],c[2],c[3],c[4],R,r[j],n,j);
        }
        /*
        for(int j=0; j<nB; j++)
        {
            dbl R=NUM::PolvEval(rB[j],c,deg);
            if(fabs(R)>1e-5)  // 8
                fprintf(stderr,"[%6dB] %g %g %g %g %g -> %g for %g (n=%d,j=%d)\n",
                    i,c[0],c[1],c[2],c[3],c[4],R,rB[j],nB,j);
        }
        if(n!=nB)
        {
            fprintf(stderr,"[%6d*]: n=%d, nB=%d:",i,n,nB);
            for(int j=0; j<n; j++)  fprintf(stderr," %g",r[j]);
            fprintf(stderr," <->");
            for(int j=0; j<nB; j++)  fprintf(stderr," %g",rB[j]);
            fprintf(stderr,"\n");
        }*/
    }
    
    return(0);
}
#endif

---