roots.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_ROOTS_H
#define THEORETICA_ROOTS_H
 
#include "../core/function.h"
#include "../calculus/deriv.h"
#include "../autodiff/dual.h"
#include "../autodiff/dual2.h"
#include "../algebra/vec.h"
#include "../complex/complex.h"
 
 
namespace theoretica {
 
 
    template<typename RealFunction, typename Vector = vec2>
    inline std::vector<Vector> find_root_intervals(
        RealFunction f, real a, real b, unsigned int steps = 10) {
 
        std::vector<Vector> res;
        const real dx = (b - a) / (real) steps;
 
        for (unsigned int i = 0; i < steps; ++i) {
            
            const real x1 = a + i * dx;
            const real x2 = a + (i + 1) * dx;
 
            if(f(x1) * f(x2) <= 0)
                res.push_back({x1, x2});
        }
 
        return res;
    }
 
 
    template<typename RealFunction>
    inline real root_bisect(
        RealFunction f, real a, real b,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_BISECTION_ITER) {
 
        if(a > b) {
            TH_MATH_ERROR("root_bisect", a, INVALID_ARGUMENT);
            return nan();
        }
 
        if(f(a) * f(b) >= 0) {
            TH_MATH_ERROR("root_bisect", f(a) * f(b), INVALID_ARGUMENT);
            return nan();
        }
 
        real x_avg = 0.0;
        real x_min = a;
        real x_max = b;
 
        unsigned int iter = 0;
 
        while((x_max - x_min) > (2 * tol) && iter <= max_iter) {
 
            x_avg = (x_max + x_min) / 2.0;
 
            if(f(x_avg) * f(x_min) > 0)
                x_min = x_avg;
            else
                x_max = x_avg;
 
            ++iter;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_bisect", x_avg, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x_avg;
    }
 
 
    template<typename RealFunction>
    inline real root_itp(
        RealFunction f, real a, real b, real tol = OPTIMIZATION_TOL,
        unsigned int n0 = 1, real k1 = 0.0) {
 
        // Default value for k1
        if (k1 == 0.0)
            k1 = 0.2 / (b - a);
 
        if(a > b) {
            TH_MATH_ERROR("root_itp", a, INVALID_ARGUMENT);
            return nan();
        }
 
        real y_a = f(a);
        real y_b = f(b);
 
        if(y_a * y_b >= 0) {
            TH_MATH_ERROR("root_itp", y_a * y_b, INVALID_ARGUMENT);
            return nan();
        }
 
        // Monotonicity of the function
        const int monotone = (y_a < y_b) ? 1 : -1;
 
        real x_t, x_new;
 
        const long int n_half = th::floor(th::log2((b - a) / tol));
        const long int n_max = n_half + n0;
 
        real eps = tol * pow(2.0, n_max);
        long int iter = 0;
 
        while((b - a) > (2 * tol) && iter <= n_max) {
 
 
            // Interpolation
            const real x_f = (a * y_b - b * y_a) / (y_b - y_a);
            const real x_half = 0.5 * (a + b);
 
 
            // Truncation
            const int sigma = sgn(x_half - x_f);
            const real delta = k1 * square(b - a);
 
            if (delta <= abs(x_half - x_f))
                x_t = x_f + sigma * delta;
            else
                x_t = x_half;
 
 
            // Projection
            const real r = eps - (b - a) / 2.0;
 
            if (abs(x_t - x_half) <= r)
                x_new = x_t;
            else
                x_new = x_half - sigma * r;
 
 
            // Update
            const real y_new = f(x_new);
 
            if (monotone * y_new > 0) {
                b = x_new;
                y_b = y_new;
            } else if (monotone * y_new < 0) {
                a = x_new;
                y_a = y_new;
            } else {
                return x_new;
            }
 
            eps *= 0.5;
            iter++;
        }
 
        if(abs(b - a) > 2 * tol) {
            TH_MATH_ERROR("root_itp", (a + b) / 2.0, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return (a + b) / 2.0;
    }
 
 
    template<typename RealFunction>
    inline real root_newton(
        RealFunction f, RealFunction Df, real guess = 0.0,
        real tol = OPTIMIZATION_TOL, unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
 
        real x = guess;
        real f_x = inf();
        unsigned int iter = 0;
 
        while(abs(f_x) > tol && iter <= max_iter) {
 
            f_x = f(x);
            x = x - (f_x / Df(x));
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_newton", x, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    inline real root_newton(
        dual(*f)(dual), real guess = 0,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
 
        real x = guess;
        unsigned int iter = 0;
 
        dual s = dual(inf(), 0.0);
 
        while(abs(s.Re()) > tol && iter <= max_iter) {
 
            // Compute the function and its derivative at the same time
            s = f(dual(x, 1.0));
 
            x = x - (s.Re() / s.Dual());
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_newton", x, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    template <
        typename Type = real,
        typename ComplexFunction = std::function<complex<Type>(complex<Type>)>
    >
    inline complex<Type> root_newton(
        ComplexFunction f, ComplexFunction Df, complex<Type> guess,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
 
        complex<Type> z = guess;
        complex<Type> f_z = Type(inf());
        unsigned int iter = 0;
 
        while(f_z * f_z.conjugate() > tol * tol && iter <= max_iter) {
 
            f_z = f(z);
            z = z - (f_z / df(z));
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_newton", z.Re(), NO_ALGO_CONVERGENCE);
            return complex<Type>(nan(), nan());
        }
 
        return z;
    }
 
 
    template<typename RealFunction>
    inline real root_halley(
        RealFunction f, RealFunction Df,
        RealFunction D2f, real guess = 0,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_HALLEY_ITER) {
 
        real x = guess;
        real f_x = inf();
        unsigned int iter = 0;
 
        while(abs(f_x) > tol && iter <= max_iter) {
 
            f_x = f(x);
            const real df_x = Df(x);
 
            x = x - (2 * f_x * df_x) / (2 * square(df_x) - f_x * D2f(x));
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_halley", x, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    inline real root_halley(
        dual2(*f)(dual2), real guess = 0,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_HALLEY_ITER) {
 
        real x = guess;
        dual2 s = dual2(inf(), 0.0, 0.0);
        unsigned int iter = 0;
 
        while(abs(s.Re()) > tol && iter <= max_iter) {
 
            // Compute the function value and the first and
            // second derivatives at the same time.
            s = f(dual2(x, 1, 0));
 
            const real f_x = s.Re();
            const real df_x = s.Dual1();
            const real d2f_x = s.Dual2();
 
            x = x - (2 * f_x * df_x) / (2 * square(df_x) - f_x * d2f_x);
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_halley", x, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    template<typename RealFunction>
    inline real root_steffensen(
        RealFunction f, real guess = 0,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_STEFFENSEN_ITER) {
 
        real x = guess;
        real f_x = inf();
        unsigned int iter = 0;
 
        while(abs(f_x) > tol && iter <= max_iter) {
 
            const real f_x = f(x);
            const real g_x = (f(x + f_x) / f_x) - 1;
 
            x = x - (f_x / g_x);
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_steffensen", x, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    template<typename RealFunction>
    inline real root_chebyshev(
        RealFunction f, RealFunction Df,
        RealFunction D2f, real guess = 0,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_CHEBYSHEV_ITER) {
 
        real x = guess;
        real f_x = inf();
        unsigned int iter = 0;
 
        while(abs(f_x) > tol && iter <= max_iter) {
 
            f_x = f(x);
            const real df_x = Df(x);
            const real u = f_x / df_x;
 
            x = x - u - square(u) * D2f(x) / (2.0 * df_x);
 
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_chebyshev", x, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    inline real root_chebyshev(
        dual2(*f)(dual2), real guess = 0,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_CHEBYSHEV_ITER) {
 
        real x = guess;
        dual2 s = dual2(inf(), 0.0, 0.0);
        unsigned int iter = 0;
 
        while(abs(s.Re()) > tol && iter <= max_iter) {
 
            s = f(dual2(x, 1.0, 0.0));
 
            const real f_x = s.Re();
            const real df_x = s.Dual1();
            const real u = f_x / df_x;
 
            x = x - u - square(u) * s.Dual2() / (2.0 * df_x);
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_chebyshev", x, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    template<typename RealFunction>
    inline real root_ostrowski(
        RealFunction f, RealFunction Df, real guess = 0.0,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_OSTROWSKI_ITER) {
 
        real x = guess;
        real f_x = inf();
        unsigned int iter = 0;
 
        while(abs(f_x) > tol && iter <= max_iter) {
 
            f_x = f(x);
            const real df_x = Df(x);
            const real u = f_x / df_x;
            const real f_xu = f(x - u);
 
            x = x - u - (f_xu / df_x) * (f_x / (f_x - 2 * f_xu));
 
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_ostrowski", x, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    template<typename RealFunction>
    inline real root_jarrat(
        RealFunction f, RealFunction Df, real guess = 0.0,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_JARRAT_ITER) {
 
        real x = guess;
        real f_x = inf();
        unsigned int iter = 0;
 
        while(abs(f_x) > tol && iter <= max_iter) {
 
            f_x = f(x);
            const real df_x = Df(x);
            const real u = f_x / df_x;
            const real f_xu = Df(x - 2.0 * u / 3.0);
 
            x = x - 0.5 * u + f_x / (df_x - 3 * f_xu);
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_jarrat", x, NO_ALGO_CONVERGENCE);
            return nan();
        }
 
        return x;
    }
 
 
    template<typename RealFunction>
    inline std::vector<real> roots(
        RealFunction f, real a, real b,
        real tol = OPTIMIZATION_TOL, real steps = 10) {
 
        if(steps == 0) {
            TH_MATH_ERROR("roots", steps, DIV_BY_ZERO);
            return {nan()};
        }
 
        // Find candidate intervals
        std::vector<vec2> intervals = find_root_intervals(f, a, b, steps);
 
        std::vector<real> res;
        res.reserve(intervals.size());
 
        // Iterate through all candidate intervals and refine the results
        for (unsigned int i = 0; i < intervals.size(); ++i) {
 
            // Check whether the extremes of the candidate intervals
            // happen to coincide with the roots
            if(abs(f(intervals[i][0])) <= MACH_EPSILON) {
                res.push_back(intervals[i][0]);
                continue;
            }
 
            if(abs(f(intervals[i][1])) <= MACH_EPSILON) {
                res.push_back(intervals[i][1]);
                continue;
            }
 
            // Approximate the roots using bisection inside each interval
            res.push_back(
                root_bisect(f, intervals[i][0], intervals[i][1], tol)
            );
        }
 
        return res;
    }
 
 
    template<typename Field>
    inline std::vector<Field> roots(
        const polynomial<Field>& p,
        real tolerance = OPTIMIZATION_TOL,
        unsigned int steps = 0) {
 
        // Effective order of the polynomial
        const unsigned int n = p.find_order();
        p /= p.coeff[n];
 
        // Absolute value of the highest coefficient
        Field a_hi = abs(p.coeff[n]);
        Field a_sum = 0;
 
        // Sum the absolute values of the lesser coefficients
        for (unsigned int i = 0; i < n; ++i)
            a_sum += abs(p.coeff[i]);
 
        // The roots are bounded in absolute value by the maximum
        const Field M = max(a_hi, a_sum);
 
        // Back track from the bounds to the first sign inversion ?
        
        return roots(
            [p](real x) { return p(x); },
            -M, M, tolerance, steps ? steps : (2 * n));
    }
 
}
 
#endif