muller.h

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

#ifndef MULLER_H
#define MULLER_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 muller.h
*
*   Polynomial roots by Muller's method.
*
*   @ingroup Poly
*   @author Michael Roy <a href="mailto:[e-mail]">[e-mail]</a>
*   @version 1.0
*/

#include "polyzero.h"

/**
*   Muller'method polynomial roots 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 mullerZeros(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;

    const bool ADAPTIVE = false;
    static FloatSpecs < Float > fpSpecs;
    const int MAXIT          = 128;
    static const Float SHIFT = acos(Float(-1)) * Float(7 / 180.);
    const Float GROW         = Float(1.05);

    Polynomial < T > p = P, f;
    ComplexPolynomial q, q2;
    FloatPolynomial mp;

    Complex z0, zx, z, z1, z2, pz, pz1, pz2, bestZ, zE, shift;
    Complex h, h1, h2, d, d1, d2, ac, sq, b, D, E, x, sm, df, px;
    Float e, a = 0, da = 1, r, mpz, mpx, g;
    int i, j = 0, conv = 0, n;
    bool started;

    zeros.clear();

    modulus(p, mp);
    scalePoly(p, mp);
    f = p;
    i = removeNullZeros(p);
    zeros.assign(i, Complex(0));

    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);
            if (polish) {
                newtonZero(P, z, pz, mpz, ADAPTIVE);
                newtonZero(P, z1, pz, mpz, ADAPTIVE);
            }
            zeros.push_back(z);
            zeros.push_back(z1);
            break;
        }

        modulus(p, mp);
        mp[0] = -mp[0];
        Float c = cauchyLowerBound(mp, Float(0));
        r = c;
        g = zerosGeometricMean(p);
        g = (g - c) / MAXIT;
        for (i = 1; i <= MAXIT; ++i) {

            shift = exp(Complex(0, a));
            z2 = Complex(0, -r) * shift;
            z1 = 0;
            z = (r) * shift;
            pz = eval(p, z);
            pz1 = eval(p, z1);
            pz2 = eval(p, z2);
            if (abs(pz2) < abs(pz)) {
                std::swap(pz, pz2);
                std::swap(z, z2);
            }

            a += SHIFT;
            r += g;

            e = abs(pz) * (1 + fpSpecs.MRE), mpx, mpz;
            started = false;

            h1 = z1 - z2;
            h2 = z - z1;

            if (h1 == Complex(0) && h2 == Complex(0) || (h1 + h2) == Complex(0))
                continue;
            d1 = (pz1 - pz2) / h1;
            d2 = (pz - pz1) / h2;
            d = (d2 - d1) / (h2 + h1);

            px = pz;
            mpx = mpz = abs(pz);

            j = 0;
            conv = 0;
            while (1) {
                if (++j > 3)
                    started = true;

                if (!started || conv == 0) {
                    b = d2 + (h2 * d);
                    ac = pz * d;
                    ac = ac + ac;
                    ac = ac + ac;
                    D = sqrt((b * b) - ac);

                    sm = b + D;
                    df = b - D;
                    E = (abs(df) > abs(sm)) ? df : sm;
                    if (E == Complex(0)) {
                        //printf("e");
                        break;
                    }
                    h = (pz + pz) / E;
                }
                else {
                    //  if we've lost convergence, reduce step by half until we converge again
                    //printf("h");
                    h = h * Float(.5);
                }

                x = z - h;

                if (x == z) {
                    if (j == 1)
                        evalAndDeflate(p, z, q);
                    //printf("a");
                    break;
                }

                px = evalAndDeflate(p, x, q, e);
                mpx = abs(px);

                if (_isnan(mpx)) {
                    //printf("b");
                    break;
                }
                else if (mpx <= (mpz + mpz))
                    conv = 0;
                else
                    ++conv;
                if (conv > 5) {
                    //printf("g");
                    conv = 0;
                    break;
                }

                if (!started || conv == 0) {
                    conv = 0;

                    mpz = mpx;

                    z2 = z1;
                    z1 = z;
                    z = x;

                    if (mpz < e)
                        break;

                    pz2 = pz1;
                    pz1 = pz;
                    pz = px;

                    h1 = h2;
                    d1 = d2;
                    h2 = z - z1;
                    if (h2 == Complex(0)) {
                        //printf("c");
                        break;
                    }
                    d2 = (pz - pz1) / h2;
                    if ((h2 + h1) == Complex(0)) {
                        //printf("d");
                        break;
                    }
                    d = (d2 - d1) / (h2 + h1);
                }
            }

            //printf("j:%d%s", j, (mpz < e) ? "\n" : "f ");
            if (mpz < e || mpz == 0) {
                if (sizeof(T) == sizeof(U)) {
                    z1 = conj(z);
                    pz1 = eval(q, conj(z), e);
                    if (abs(pz1) < e) {
                        //if (abs(z.imag()) > abs(z.real()) * 20 * fpSpecs.MRE) {
                        if (polish) {
                            newtonZero(P, z, pz, mpz, ADAPTIVE);
                            evalAndDeflate(p, z, q);
                        }
                        zeros.push_back(z);
                        z = conj(z);
                        if (polish) {
                            newtonZero(P, z, pz, mpz, ADAPTIVE);
                        }
                        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, ADAPTIVE);
                            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);
                    }
                    break;
                }
                else {
                    if (polish) {
                        newtonZero(P, z, pz, mpz, ADAPTIVE);
                        evalAndDeflate(p, z, q);
                    }
                    zeros.push_back(z);
                    * (ComplexPolynomial *) & p = q;
                    break;
                }
            }
        }
        if (i >= MAXIT)
            break;
    }
    sortZeros(zeros);
    return (int)zeros.size();
}

#endif //MUELLER_H