/ray/src/lib/numerics/rootsolve/rootsolve.h

Go to the documentation of this file.

/*
 * lib/numerics/rootsolve/rootsolve.h
 * 
 * Numerics library root solving functions. 
 * 
 * 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.
 * 
 */

#ifndef _LIB_NUMERICS_ROOTSOLVE_ROOTSOLVE_H_
#define _LIB_NUMERICS_ROOTSOLVE_ROOTSOLVE_H_ 1

#include <lib/numerics/num_math.h>
#include <lib/salloc.h>


namespace NUM  // numerics
{

class PolyRootSolver
{
    private:
    public:
        PolyRootSolver() {}
        virtual ~PolyRootSolver() {}
        
        virtual int solve(int order,const dbl *c,dbl *r);
};


class PolyRootSolver_Direct : public PolyRootSolver
{
    public:
        struct Params
        {
            dbl nearly_zero;
        };
        Params par;
        
        static const Params defaults;
        
    public:
        PolyRootSolver_Direct();
        ~PolyRootSolver_Direct();
        
        int solve(int order,const dbl *c,dbl *r);
        
        int Solve_1(const dbl *c,dbl *r);  // Simple
        int Solve_2(const dbl *c,dbl *r);  // Midnight formula
        int Solve_3(const dbl *c,dbl *r);  // Cardano's formula
        int Solve_4(const dbl *c,dbl *r);  // Vieta's "formula"
};


class PolyRootSolver_Sturm : public PolyRootSolver
{
    public:
        struct Params
        {
            dbl nearly_zero;
            dbl relerror;
            int sbisect_maxit;
            int modrf_maxit;
        };
        Params par;
        
        static const Params defaults;
        
    private:
#define PolyRootSolver_Sturm_BufOnStack 1
        
        struct poly
        {
            int ord;        
            dbl *coef;   
            
            inline poly() : ord(0),coef(NULL) {}
            
#if !PolyRootSolver_Sturm_BufOnStack
            inline void AllocCoeffs(int order)
                {  coef=REALLOC(coef,order+1);  }
            inline ~poly()  {  coef=FREE(coef);  }
#else
            inline ~poly() {}
#endif
        };
    
        int modrf(int ord,dbl *coef,dbl a,dbl b,dbl *val);
        int modp(poly *u,poly *v,poly *r);
        int buildsturm(int ord,poly *sseq);
        int numroots(int np,poly *sseq,int *atneg,int *atpos);
        int numchanges(int np,poly *sseq,dbl a);
        int sbisect(int np,poly *sseq,dbl min,dbl max,int atmin,
            int atmax,dbl *roots);
    public:
        PolyRootSolver_Sturm();
        ~PolyRootSolver_Sturm();
        
        int solve(int order,const dbl *c,dbl *r);
};


class PolyRootSolver_Bairstow : public PolyRootSolver
{
    public:
        struct Params
        {
            dbl nearly_zero;
            int maxiter;
            int iter_chg;
            dbl initial_epsilon;
            dbl disc_epsilon;
        };
        Params par;
        
        static const Params defaults;
        
    private:
        int _QuadricRealRoots(dbl *a,int n,dbl *wr);
        int _QuadricRoots(dbl *a,int n,dbl *wr,dbl *wi);
        int _FindQuadFact(dbl *a,int n,dbl *b,dbl *quad,
            int red_fact,dbl *err_est);
        void _GetQuads(dbl *a,int n,dbl *quad,dbl *x);
        int _Recurse(dbl *a,int n,dbl *b,int m,dbl *quad,dbl *err);
        void _DiffPoly(dbl *a,int n,dbl *b);
        dbl _Deflate(dbl *a,int n,dbl *b,dbl *quad);
        
    public:
        PolyRootSolver_Bairstow();
        ~PolyRootSolver_Bairstow();
        
        int solve(int order,const dbl *c,dbl *r);
};

}  // end of namespace NUM

#endif  /* _LIB_NUMERICS_ROOTSOLVE_ROOTSOLVE_H_ */

---