jenkinstraub.h

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

#ifndef JENKINSTRAUB_H
#define JENKINSTRAUB_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 jenkinstraub.h
*
*   Polynomial roots by Jenkins-Traub method.
*
*   @ingroup Poly
*   @author Michael Roy <a href="mailto:[e-mail]">[e-mail]</a>
*   @version 1.0
*/

#include "polyzero.h"

/**
*   Jenkins-Traub polynomial zeros evaluator.
*
*   returns all zeros of polynomial \f$P(x)\f$
*
*   @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 jenkinsTraubZeros(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;

    class local {
    public:
        static void noShift(const ComplexPolynomial& p, ComplexPolynomial& h) {
            const unsigned M = 5;
            int m = p.degree();
            Float ii = Float(1);
            Float dd = Float(m);
            h.resize(m);
            //  first evalutate H0 = p'
            for (int i = 0; i < m; ++i, ii += Float(1))
                h[i] = p[i + 1] * (ii / dd);

            // Hl+1(z) = 1/z * (Hl(z) - (Hl(0)/p(0)) * p(z));
            for (unsigned i = 0; i < M; ++i) {
                Complex k = -h[0] / p[0];
                for (int j = 0; j < m - 1; ++j)
                    h[j] = h[j + 1] + (k * p[j + 1]);
                h[m - 1] = k * p[m];
            }
        }

        static bool fixedShift(unsigned pass, const ComplexPolynomial& p, ComplexPolynomial& h, ComplexPolynomial& qp, Complex& s) {

            bool failed = false;
            int conv = 0;
            Complex hs, ps, t, tOld, z;
            ComplexPolynomial qh;
            ps = p.evalAndDeflate(s, qp);
            t = nextT(s, ps, h, hs, qh);

            for (unsigned i = 0; i < (10 * pass); ++i) {
                tOld = t;
                //  compute next h and new t
                nextH(p, h, qp, qh, hs, ps, t);
                t = nextT(s, ps, h, hs, qh);
                z = s + t;

                if (hs == Complex(0))
                    continue;

                // weak convergence test
                Float k = abs(tOld - t);
                Float l = Float(.5) * abs(z);
                if (!failed && abs(t - tOld) < Float(.5) * abs(z)) {
                    if (++conv >= 2) {
                        if (variableShift(p, h, qp, qh, z, t, hs, ps)) {
                            s = z;
                            return true;
                        }
                        conv = 0;
                        failed = true;
                    }
                }
                else
                    conv = 0;
            }
            // do a last attempt with final h polynomial
            if (variableShift(p, h, qp, qh, z, t, ps, hs)) {
                s = z;
                return true;
            }
            return false;
        }

        static bool variableShift(const ComplexPolynomial& p, ComplexPolynomial& h, ComplexPolynomial& qp, ComplexPolynomial& qh, Complex& s, Complex& t, Complex& ps, Complex& hs) {
            static FloatSpecs<Float> fpSpecs;
            bool b = false;
            int i;
            Float step = 0, mpOld = fpSpecs.INFINY, mp = 0, mt = fpSpecs.INFINY, mtOld = fpSpecs.INFINY;
            ComplexPolynomial _H = h, _QH = qh;
            Complex _s = s, _Hs = hs, _t = t;
            Float e;

            for (i = 1; i <= 12; ++i) {
                ps = evalAndDeflate(p, s, qp, e);
                mpOld = mp;
                mp = abs(ps);
                mtOld = mt;
                mt = abs(t);

                if (_isnan(mp)) {
                    //printf("a");
                    break;
                }

                if (mp < e) {
                    //printf("i:%d ", i);
                    return true;
                }

                if (i != 0) {
                    if(!b && mp > mpOld && step < Float(.05)) {
                        // Iteration has stalled. Probably a cluster of zeros. Do 5 fixed 
                        // shift steps into the cluster to force one zero to dominate
                        Float tp = std::max(step, fpSpecs.ETA);
                        b = true;
                        Complex r(1 + sqrt(tp), sqrt(tp));
                        s *= r;
                        ps = p.evalAndDeflate(s, qp);
                        for(unsigned j = 0; j < 5; ++j) {
                            t = nextT(s, ps, h, hs, qh);
                            nextH(p, h, qp, qh, hs, ps, t);
                        }
                        mp = fpSpecs.INFINY;
                    }
                    // Exit if polynomial value increase significantly
                    else if((mp * Float(0.1)) > mpOld) 
                        break;
                }

                t = nextT(s, ps, h, hs, qh);
                nextH(p, h, qp, qh, hs, ps, t);
                t = nextT(s, ps, h, hs, qh);
                if (hs != Complex(0)) {
                    step = abs(t) / abs(s);
                    s = s + t;
                }
            }
            h = _H;
            qh = _QH;
            s = _s;
            t = _t;
            hs = _Hs;
            ps = evalAndDeflate(p, s, qp);
            //printf("i:%d ", i);
            return false;
        }

        static Complex nextT(const Complex& s, const Complex& ps, const ComplexPolynomial& h, Complex& hs, ComplexPolynomial& qh) {
            static FloatSpecs<Float> fpSpecs;
            Complex r;
            Float e;
            hs = evalAndDeflate(h, s, qh, e);
            if (abs(hs) < fpSpecs.ARE * (10 * abs(h[0]))) {
            //if (abs(hs) < e) {
                hs = Complex(0);
                return Complex(0);
            }
            return -ps / hs;
        }

        static void nextH(const ComplexPolynomial& p, ComplexPolynomial& h, const ComplexPolynomial& qp, const ComplexPolynomial& qh, const Complex& hs, const Complex& ps, const Complex& t) {
            int m = qp.degree();
            int n = qh.degree();
            int i;
            if (hs != Complex(0)) {
                // Hl+1(z) = 1/(z-s) * Hl(z) - (Hl(s)/p(s)) * p(z))
                h.resize(m + 1);
                for(i = 0; i <= n; ++i)
                    h[i] = (t * qh[i]) + qp[i];
                for( ; i <= m; ++i)
                    h[i] = qp[i];
            }
            else
                h = qh;
        }
    };

    const bool ADAPTIVE = false;
    const int MAXIT1 = 2;
    const int MAXIT2 = 9;
    static const Complex SHIFT = exp(Complex(0, Float(-94) * acos(Float(-1)) / Float(180)));

    ComplexPolynomial p, qp, qp2, h;
    FloatPolynomial mP;
    Complex z, pz, z1, pz1;
    Float e, mpz;
    int i, n;

    p.resize(P.degree() + 1);
    for (i = 0; i < (int)p.size(); ++i)
        p[i] = Complex(P[i]);

    zeros.clear();

    modulus(p, mP);
    scalePoly(p, mP);

    i = removeNullZeros(p);
    zeros.assign(i, Complex(0));

    Complex shift = SHIFT;
    while (p.degree() > 0) {

        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], z, z1);
            zeros.push_back(z);
            zeros.push_back(z1);
            break;
        }

        modulus(p, mP);
        mP[0] = -mP[0];
        Float beta = cauchyLowerBound(mP) * Float(1.05);
        bool zeroFound = false;

        for (i = 1; i < MAXIT1 && !zeroFound; ++i) {
            // stage1, no shift
            local::noShift(p, h);
            for (unsigned j = 1; j < MAXIT2 && !zeroFound; ++j) {
                // stage2, fixed shift
                z = beta * shift;
                shift *= SHIFT;
                if (local::fixedShift(j, p, h, qp, z)) {
                    zeroFound = true;
                    break;
                }
            }
        }
        if (zeroFound) {
            //printf("\n");
            if (sizeof(T) == sizeof(U)) {
                z1 = conj(z);
                pz1 = eval(qp, conj(z), e);
                if (abs(pz1) < e) {
                    if (polish) {
                        newtonZero(P, z, pz, mpz, ADAPTIVE);
                        evalAndDeflate(p, z, qp);
                    }
                    zeros.push_back(z);
                    z = conj(z);
                    if (polish) {
                        newtonZero(P, z, pz, mpz, ADAPTIVE);
                    }
                    evalAndDeflate(qp, z, qp2);
                    zeros.push_back(z);
                    n = qp2.degree();
                    for (int i = 0; i <= n; ++i)
                        *(Float*)&(p[i]) = qp2[i].real();
                    p.resize(n + 1);
                }
                else {
                    if (polish) {
                        newtonZero(P, z, pz, mpz, ADAPTIVE);
                        evalAndDeflate(p, z, qp);
                    }
                    zeros.push_back(z);
                    n = qp.degree();
                    for (int i = 0; i <= n; ++i)
                        *(Float*)&(p[i]) = qp[i].real();
                    p.resize(n + 1);
                }
            }
            else {
                if (polish) {
                    newtonZero(P, z, pz, mpz, ADAPTIVE);
                    evalAndDeflate(p, z, qp);
                }
                zeros.push_back(z);
                *(ComplexPolynomial*)&p = qp;
            }
        }
        else 
            break;
    }
    sortZeros(zeros);
    return (int)zeros.size();
}

#endif //JENKINSTRAUB_H