polyzero.h

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

#ifndef POLYZERO_H
#define POLYZERO_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 polyzero.h
*
*   polynomial zero-finding support functions.
*
*   This file contains helper routines to polynomial roots-solving.
*
*   @ingroup Poly
*
*   @author Michael Roy <a href="mailto:[e-mail]">[e-mail]</a>
*
*   @version 1.0
*
*   @see    class Polynomial\n
*           laguerre.h\n
*           muller.h
*/

#include "polynomial.h"

/**
*   @defgroup PolyZero Polynomial roots-solving functions
*   @ingroup  Poly
*   Operations specifically related to the general polynomial roots-solving
*   problem.
*
*/
/** @addtogroup Poly
*   @{ */
/**
*   resolves closest zero by Newton-Raphson iterations.
*
*   This function requires an initial zero estimate that is fairly close.
*
*   @param[in]  p       polynomial being evaluated.  template parameter 
*                       \e T can be any IEEE floating-point type, real or 
*                       complex.
*   @param[in,out] z    on entry, the original estimate for zero. On exit, the 
*                       refined zero. Template parameter U can be IEEE 
*                       floating-point type, real or complex.
*   @param[out] pz      on exit, contains \f$p(z)\f$
*   @param[out] mpz     on exit, contains \f$|p(z)|\f$, the asolute value of 
*                       \f$p(z)\f$.
*   @param[in]  adaptive    indicates wether to use the classic Newton-Raphson 
*                       algorithm, or the adaptive version.  The adaptive 
*                       algorithm is more efficient for resolving multiple zeros, 
*                       but can be much is slower when the original estimate is
*                       very close to the actual zero.
*
*   @ingroup PolyZero
*   @see Polynomial
*/

template <class T, class U, class V>
bool newtonZero(const Polynomial<T>& p, U& z, U& pz, V& mpz, bool adaptive = false) {
    const int MAX_RATIO = 6;
    U ppz, z1, z2, dz, y, y1, y2;
    V my, my1, my2;
    int i, j = 0, m = p.degree();

    if (m < 1)
        return false;

    if (adaptive) {
        z1 = z;
        do {
            ++j;
            y = evalAndDerive(p, z1, ppz);
            my = abs(pz);
            if (ppz == U(0))
                break;
            z = z1;
            pz = y;
            mpz = my;
            dz = pz / ppz;

            z1 = z - dz;
            ++j;
            y1 = eval(p, z1);
            my1 = abs(y1);

            z2 = z1 - dz;
            ++j;
            y2 = eval(p, z2);
            my2 = abs(y2);

            i = 3;
            while(my2 < my1 && i++ < MAX_RATIO) {
                z1 = z2;
                y1 = y2;
                my1 = my2;
                z2 = z2 - dz;
                ++j;
                y2 = eval(p, z2);
                my2 = abs(y2);
            }
            if (_isnan(my1))
                break;
        } while (my1 < my);
    }
    else {
        ++j;
        y = evalAndDerive(p, z, ppz);
        my = abs(y);
        z1 = z;
        do {
            z = z1;
            pz = y;
            mpz = my;
            if (ppz == U(0))
                break;
            z1 -= pz / ppz;
            if (z1 == z)
                break;
            ++j;
            y = evalAndDerive(p, z1, ppz);
            my = abs(y);
            if (_isnan(my))
                break;
        } while(my < mpz);
    }
    //printf("j:%d, mpz:%le\n", j, mpz);
    return true;
}

/**
*   divides a polynomial \f$p(x)\f$ by \f$x\f$, until is has no null roots.
*
*   Deflates the polynomial <i>p</i> until it has no null roots, effectively 
*   dividing \f$p(x)\f$ by \f$x\f$ until the coeficient of degree zero
*   is not null.
*
*   @par Example:
*   @include nullzeros.cpp
*
*   @par Program output:
*   @include nullzeros.txt
*
*   @param[in]  p   polynomial being evaluated.  template parameter \e T can 
*                   be any IEEE floating-point type, real or complex.
*   @returns The number of null roots found and removed from the original 
*               polynomial.
*
*   @ingroup PolyZero
*   @see    Polynomial \n
*           Polynomial::shift()
*/

template<class T>
int removeNullZeros(Polynomial<T>& p) {
    int i, m = p.degree();
    for (i = 0; i < m && p[i] == T(0); ++i)
        ;
    p.shift(-i);
    return i;
}

/**
*   returns the Cauchy lower bound for a polynomial.
*
*   approximates the Cauchy lower bound for polynomial \a p.
*
*   This function expects the input polynomial to be of the form:
*
*   \f[p(x)=-|a_0|+|a_1|x+|a_2|x^2+...+|a_n|x^n\f]
*
*   where \f$(a_0,a_1,a_2,...,a_n)\f$ are real or complex coeficients of the
*   polynomial being examined.  The Cauchy lower bound is the value of the only
*   positive zero of \f$p(x)\f$.  This zero is found using Newton-Raphson 
*   method.
*
*   The optional parameter \a upperBound indicates a starting point from
*   where to start looking for this zero. If \a upperBound is omited, or 
*   null, the geometric mean of the zeros will be used as a starting point for
*   the search.
*
*   @par Example:
*   @code
#include <complex>
#include "polynomial.h"
#include "polyzero.h"

void main() {
    //  this example shows how to compute the Cauchy lower bound
    //  of a polynomial
    Polynomial< std::complex<double> > p;
    Polynomial<double> mp;

    //  initialize polynomial p
    ...

    //  compute Cauchy polynomial
    modulus(p, mp)
    mp[0] = -mp[0];
    
    //  compute Cauchy lower bound of polynomial p
    double c = cauchyLowerBound(mp, 0);

    // do something
    ...
}
*   @endcode
*
*   @param[in]  p           polynomial being evaluated.  Template parameter 
*                           \e T can be any non-complex IEEE floating-point 
*                           type.
*   @param[in]  upperBound  starting point to use for the search.
*   @returns The Cauchy lower bound.
*
*   @see    Polynomial\n
*           newtonZero()
*/

template<class T>
T cauchyLowerBound(const Polynomial<T>& p, const T& upperBound = T(0)) {
    static FloatSpecs<T> fpSpecs;
    T x, y, my;
    if (upperBound <= 0) {
        x = zerosGeometricMean(p);
        if (p[1] != 0)
            x = std::min(x, -p[0] / p[1]);
    }
    else
        x = upperBound;

    newtonZero(p, x, y, my);
    return x;
}

/**
*   returns the Cauchy upper bound for a polynomial
*
*   Returns the Cauchy upper bound for polynomial \a p.
*
*   The Cauchy upper bound \f$C_u\f$ is the result of the equation
*
*   \f[C_u=1+max_n(|a_0|,|a_1|,...,|a_n|)\f]
*
*   @param[in]  p   polynomial being evaluated.  template parameter 
*                   \e T can be any real or complex IEEE floating-point
*                   type.
*
*   @see    Polynomial\n
*           cauchyLowerBound()
*/

template<class T>
T cauchyUpperBound(const Polynomial<T>& p) {
    int m = p.degree();
    T x = abs(p[0]);
    for (int i = 1; i <= m; ++i)
        x = std::max(x, abs(p[i]));
    return x + T(1);
}

/**
*   returns the geometric mean for zeros of a polynomial.
*
*   Returns the geometric mean for zeros of a polynomial.
*
*   According to Vieta's theorem, the geometric mean \f$G\f$ for the zeros 
*   of a polynomial can be simplified to the result of the equation:
*
*   \f[G_p=\left[\frac{|a_0|}{|a_n|}\right]^{\frac{1}{n}}, p(x)=a_0+a_1{x}+a_2{x^2}+...+a_n{x^n}, a_n\neq 0\f]
*
*   @param[in]  p   polynomial being evaluated.  template parameter 
*                   \e T can be any real or complex IEEE floating-point
*                   type.
*
*   @see    Polynomial\n
*           cauchyLowerBound()
*/

template<class T>
T zerosGeometricMean(const Polynomial<T>& p) {
    const T SMALL = T(1.e-3);
    T m = Float(p.degree());
    return exp(log(abs(p[0]) - log(abs(p[p.degree()])))/m + SMALL);
}

/**
*   returns the geometric mean for zeros of a polynomial.
*
*   See zerosGeometricMean(const Polynomial<T>& p) for more details.
*/

template<class T>
T zerosGeometricMean(const Polynomial< std::complex <T> >& p) {
    const T SMALL = T(1.e-3);
    T m = Float(p.degree());
    return exp(log(abs(p[0]) - log(abs(p[p.degree()])))/m + SMALL);
}

/**
*   sorts the zeros of a polynomial.
*
*   Sorts the zeros of a polynomial by ascending order of their real 
*   part.
*
*   @param[in,out]  zeros   zeros of a polynomial.
*
*   @ingroup PolyZero
*   @see    Polynomial\n
*           laguerreZeros()\n
*           mullerZeros()\n
*           jenkinsTraubZeros()
*/

template<class T>
inline void sortZeros(std::vector<std::complex<T> >& zeros) {
    class local {
    public:
        static inline bool zeroSortPred(const std::complex<T>& z1, const std::complex<T>& z2) {
            return z1.real() < z2.real();
        }
    };
    std::sort(zeros.begin(), zeros.end(), local::zeroSortPred);
}

/**
*   returns the root of a degree 1 polynomial.
*
*   Returns the unique root of the polynomial \f$p(x)=ax+b\f$
*
*   @param[in]  a   the degree 1 coeficient.
*   @param[in]  b   the degree 0 coeficient.
*   @returns        The fraction \f$-\frac{b}{a}\f$
*
*   @ingroup PolyZero
*/
template <class T>
inline T solveDegree1(const T& a, const T& b) {
    return -b / a;
}

/**
*   returns the root of a degree 2 polynomial.
*
*   Returns the two roots of the polynomial \f$p(x)=ax^2+bx+c\f$.
*
*   The roots \a z1 and \a z2 are calculated as follows.
*
*   \f[S=b+\sqrt{b^2-4ac}\f]
*   \f[D=b-\sqrt{b^2-4ac}\f]
*   \f[q=\left\{ \begin{array}{ll} -S & \mbox{if} \left| S\right| >\left| D\right| \\ -D & \mbox{otherwise}\end{array}\right. \f]
*   \f[z_1=\frac{2c}{q}\f]
*   \f[z_2=\frac{q}{2a}\f]
*
*   This method avoids dividing by a very small number and thus reduces error.
*
*   @param[in]  a   the degree 2 coeficient.  Template parameter \a T can be 
*                   any floating-point or cvomplex data type.
*   @param[in]  b   the degree 1 coeficient.
*   @param[in]  c   the degree 0 coeficient.
*   @param[out] z1  the root with the smallest absolute value.  Template 
*                   parameter \a U should be a complex type for which
*                   operators *, -, / and functions abs() and sqrt() are 
*                   defined.
*   @param[out] z2  the root with the largest absolute value.

*   @ingroup PolyZero
*/

template <class T, class U>
void solveDegree2(const T& a, const T& b, const T& c, U& z1, U& z2) {
    U sq, sum, dif, q;
    sq = sqrt(U((b * b) - (Float(4) * a * c)));
    sum = b + sq;
    dif = sq - b;
    if (abs(dif) > abs(sum))
        q = dif;
    else 
        q = -sum;

    z1 = (c + c) / q;
    z2 = q / (a + a);

    if (z1 > z2)
        std::swap(z1, z2);
}
/**@} Poly */
#endif // POLYZERO_H