polyzero.h

Go to the documentation of this file.

#ifndef POLYZERO_H
#define POLYZERO_H

#include "polynomial.h"

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;
}

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;
}

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;
}

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);
}

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);
}

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);
}

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);
}

template <class T>
inline T solveDegree1(const T& a, const T& b) {
    return -b / a;
}

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);
}
#endif // POLYZERO_H 

---