laguerre.h

This library defines basic operation on polynomials, and contains also 3 different roots (zeroes)-fi

#ifndef LAGUERRE_H
#define LAGUERRE_H
/** 
*       This program has been developed by Micha雔 Roy.  
*	You have the author's permission to incorporate this code 
*       into your product, royalty free.  

*       The author specifically disclaims all warranties, express or
*       implied, and all liability, including consequential and
*       other indirect damages, for the use of this code, 
*       including liability for infringement of any proprietary
*       rights, and including the warranties of merchantability
*       and fitness for a particular purpose.  And he does not 
*       assume any responsibility for any errors which may 
*       appear in this code nor any responsibility to update it.
*
*   (c) 2005, Michael Roy.  [e-mail]
*/
/** 
*   @file laguerre.h
*
*   Polynomial roots by a modified Laguerre's method.
*
*   @ingroup Poly
*   @author Michael Roy <a href="mailto:[e-mail]">[e-mail]</a>
*   @version 1.0
*/

#include "polyzero.h"

/**
*   modified Laguerre polynomial zeros evaluator.
*
*   Returns all zeros of polynomial \f$P(x)\f$.
*
*   This is a very good roots evaluator, that can be pushed to the limits.  
*   It will quickly find the roots of even quite large polynomials 
*   (>1000 coeficients) in a reasonable amount of time.
*
*   It is also very precise and the roots usually will not need any further 
*   polishing.
*
*   @param[in]  P       polynomial to evaluate. \e T coeficients type can be 
*                       either real or complex.
*
*   @param[out] zeros   vector to hold zeros.
*
*   @param[in]  polish  if set to \e true, the zeros will be polished by 
*                       Newton-Raphson algorithm on the original polynomial 
*                       before deflation.  This 'quick and dirty' polishing 
*                       can lead to problems on large polynomials.
*
*   @returns    The number of zeros found.  This number can be checked against 
*               the degree of \f$P(x)\f$ to check for errors.
*
*   @ingroup PolyZero
*   @see class Polynomial \n
*              newtonZero()
*/

template<class T, class U>
int laguerreZeros(const Polynomial< T >& P, std::vector< std::complex < U > >& zeros, bool polish) {

    typedef U Float;
    typedef std::complex<Float> Complex;
    typedef Polynomial<Complex> ComplexPolynomial;
    typedef Polynomial<Float> FloatPolynomial;

    // local scope functions 
    class local {
    public:
        static Complex fejer(Float m, const Complex& pz, const Complex& ppz, const Complex& pppz) {
            Complex t1, t2, t3, z1, z2;
            t1 = pppz / (m * (m - 1));
            t2 = Float(2) * ppz / m;
            t3 = pz;
            solveDegree2(t1, t2, t3, z1, z2);
            return z1;
        }

        static Complex laguerreStep(Float m, const Complex& zs, const Complex& pz, const Complex& ppz) {
            Complex t;
            t = ((m - 2) * ppz) / (m * pz);
            t = (zs * t) + (m - 1);
            return zs / t;
        }
    };

    static FloatSpecs<Float> fpSpecs;
    const Float SMALL = Float(1.e-3);
    const Float BIG_ONE = Float(1.0001);
    const Float SML_ONE = Float(0.9999);
    const Float RCONST = Float(1.445);
    const Float ONEPQT = Float(1.25);
    const Float GAMA = Float(1);
    const Float THETA = Float(2);

    Polynomial<T> p = P;
    ComplexPolynomial q, q2;
    FloatPolynomial mP, d, p2, rem;
    int startd, iter, spiral, n;
    Complex spir, temp, r;
    Float sqm, m, ratio;

    Complex zs, l0, step, l, z0, z, pz, ppz, pppz, pz0;
    Float c, f, g, w, ml0, ub, lb, mstep, ml, mpz, mpz0;

    d.resize(3);

    modulus(p, mP);
    scalePoly(p, mP);
    mP[0] = -mP[0];
    
    n = removeNullZeros(p);
    zeros.assign(n, Complex(0));

    while (p.degree() > 0) {

        m = Float(p.degree());

        if (p.degree() == 1) {
            z = solveDegree1(p[1], p[0]);
            zeros.push_back(z);
            break;
        }
        if (p.degree() == 2) {
            solveDegree2(p[2], p[1], p[0], z0, z);
            if (polish) {
                newtonZero(P, z, pz, mpz, false);
                newtonZero(P, z0, pz, mpz, false);
            }
            zeros.push_back(z0);
            zeros.push_back(z);
            break;
        }

        //  get an initial estimate z, p(z)
        z = Complex(0);
        spir = Complex(-ONEPQT / 10, 1);
        do {
            pz = evalAndDerive(p, z, ppz, pppz);
            //  compute zs, fejer bound
            zs = local::fejer(m, pz, ppz, pppz);
            f = abs(zs);

            //  diro = l0
            l0 = local::laguerreStep(m, zs, p[0], p[1]);
            ml0 = abs(l0);

            if (_isnan(ml0)) {
                if (z == Complex(0))
                    z = -p[0];
                else
                    z = spir * z;
            }
            else
                break;
        } while (1);

        sqm = sqrt(m);
        g = zerosGeometricMean(p);

        //  Laguerre bound w = sqrt(m).|step|;
        w = sqm * ml0;

        //  get cauchy bound
        ub = std::min(g, BIG_ONE * std::min(f, w));
        c = cauchyLowerBound(mP, ub);

        //  compute upper and lower bound of smallest zero
        ub = std::min(RCONST * m * c, ub);
        lb = SML_ONE * c;

        //  init first pass
        f = ub;
        g = ub;
        step = l0;
        mstep = ml0;
        ratio = mstep / g;
        mpz = abs(pz);
        mpz0 = mpz;
        spiral = 0;
        startd = 0;
        iter = 0;
        while (1) {
            ++iter;

            if (ratio > THETA) {
                if (startd) {
                    // we've lost convergence, reduce step by half
                    ml *= Float(.5);
                    l *= Float(.5);
                    if (ml > fpSpecs.ARE * abs(z))
                        z = z0 + l;
                    else
                        break;
                }
                else {
                    // no convergence yet, try to catch a dip with an outward 
                    // spiral
                    if (spiral) {
                        z = spir * z;
                    }
                    else {
                        spiral = 1;
                        spir = Complex(-ONEPQT / m, 1);
                        ml = lb / (m * m);
                        z =  (step / mstep) * lb;
                    }

                    if (iter >= 25) {
                        //  we're spiraling too much, might as well let laguerre run its 
                        //  course (<40 iters), or we could be going for 100ks 
                        //  iterations !
                        //  CHECK: spiraling because the lower bound is too small ?

                        mpz0 = mpz;
                        pz0 = pz;
                        z0 = z;
                        startd = 1;
                        //printf("X");
                    }
                }
            }
            else {
                // converging, limit the step to keep converging
                startd = 1;
                if (ratio > GAMA && (spiral || lb <= g)) {
                    ratio = GAMA / ratio;
                    step /= ratio;
                    mstep /= ratio;
                }
                g = f;
                l = step;
                ml = mstep;
                z0 = z;
                mpz0 = mpz;
                z += l;
            }
            _clearfp();
            
            pz = evalDeriveAndDeflate(p, z, ppz, pppz, q);
            
            mpz = abs(pz);
            if (_isnan(mpz))
                mpz = fpSpecs.INFINY;

            if (mpz == 0)
                break;
            
            if (mpz >= mpz0 && startd) {
                ratio = THETA * BIG_ONE;
            }
            else {
                // get Fejer bound
                zs = local::fejer(m, pz, ppz, pppz);
                f = abs(zs);
                
                //  compute next laguerre step
                step = local::laguerreStep(m, zs, pz, ppz);
                mstep = abs(step);

                // Laguerre bound
                w = sqm * mstep;
                f = std::min(w, f);

                //  g is Fejer bound from preceding point (z0)
                ratio = mstep / g;

                if (mstep < fpSpecs.ARE * abs(z))
                    break;
            }
        }
        // store the zero, different strategies for real and complex-valued
        // coeficients
        if (sizeof(T) == sizeof(U)) {
            if (abs(z.imag()) > abs(z.real()) * 20 * fpSpecs.MRE) {
                if (polish) {
                    newtonZero(P, z, pz, mpz);
                    evalAndDeflate(p, z, q);
                }
                zeros.push_back(z);
                z = conj(z);
                if (polish)
                    newtonZero(P, z, pz, mpz);
                evalAndDeflate(q, z, q2);
                zeros.push_back(z);
                n = q2.degree();
                for (int i = 0; i <= n; ++i)
                    *(Float*)&(p[i]) = q2[i].real();
                p.resize(n + 1);
            }
            else {
                if (polish) {
                    newtonZero(P, z, pz, mpz);
                    evalAndDeflate(p, z, q);
                }
                zeros.push_back(z);
                n = q.degree();
                for (int i = 0; i <= n; ++i)
                    *(Float*)&(p[i]) = q[i].real();
                p.resize(n + 1);
            }
        }
        else {
            if (polish) {
                newtonZero(P, z, pz, mpz);
                evalAndDeflate(p, z, q);
            }
            zeros.push_back(z);
            *(ComplexPolynomial*)&p = q;
        }
        //printf("%d\n", iter);
        modulus(p, mP);
        mP[0] = -mP[0];
    }
    sortZeros(zeros);
    return (int)zeros.size();
}

#endif //LAGUERRE_H