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"
#include "../core/iter_result.h"
 
 
namespace theoretica {
 
 
    template<typename RealFunction, typename Bracket = vec2>
    inline std::vector<Bracket> find_root_brackets(
        RealFunction f, real a, real b, unsigned int steps = 10) {
 
        std::vector<Bracket> res;
        const real dx = (b - a) / 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 iter_result<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, MathError::InvalidArgument);
            return iter_result<real>(ConvergenceStatus::InvalidInput);
        }
 
        if(f(a) * f(b) >= 0) {
            TH_MATH_ERROR("root_bisect", f(a) * f(b), MathError::InvalidArgument);
            return iter_result<real>(ConvergenceStatus::InvalidInput);
        }
 
        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, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(x_max - x_min) / 2.0);
        }
 
        return iter_result<real>(x_avg, iter, abs(x_max - x_min) / 2.0);
    }
 
 
    template<typename RealFunction>
    inline iter_result<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, MathError::InvalidArgument);
            return iter_result<real>(ConvergenceStatus::InvalidInput);
        }
 
        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, MathError::InvalidArgument);
            return iter_result<real>(ConvergenceStatus::InvalidInput);
        }
 
        // 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 {
                // Falls exactly on the zero
                return iter_result<real>(x_new, iter, 0.0);
            }
 
            eps *= 0.5;
            iter++;
        }
 
        if(abs(b - a) > 2 * tol) {
            TH_MATH_ERROR("root_itp", (a + b) / 2.0, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(b - a) / 2.0);
        }
 
        return iter_result<real>((a + b) / 2.0, iter, abs(b - a) / 2.0);
    }
 
 
    template<typename RealFunction>
    inline iter_result<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();
        real Df_x = 0;
        unsigned int iter = 0;
 
        while(abs(f_x) > tol && iter <= max_iter) {
 
            Df_x = Df(x);
 
            // Check for division by zero
            if (abs(Df_x) < MACH_EPSILON) {
                TH_MATH_ERROR("root_newton", Df_x, MathError::DivByZero);
                return iter_result<real>(ConvergenceStatus::Diverged, iter, f_x);
            }
 
            f_x = f(x);
            x = x - (f_x / Df_x);
            iter++;
        }
 
        if(iter > max_iter) {
            TH_MATH_ERROR("root_newton", x, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(f_x));
        }
 
        return iter_result<real>(x, iter, abs(f_x));
    }
 
 
    inline iter_result<real> root_newton(
        dual(*f)(dual), real guess = 0,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_NEWTON_ITER) {
 
        real x = guess;
        dual s = dual(inf(), 0.0);
        unsigned int iter = 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, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(s.Re()));
        }
 
        return iter_result<real>(x, iter, abs(s.Re()));
    }
 
 
    template <
        typename Type = real,
        typename ComplexFunction = std::function<complex<Type>(complex<Type>)>
    >
    inline iter_result<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(), MathError::NoConvergence);
            return iter_result<complex<Type>>(ConvergenceStatus::MaxIterations, iter, f_z.norm());
        }
 
        return iter_result<complex<Type>>(z, iter, f_z.norm());
    }
 
 
    template<typename RealFunction>
    inline iter_result<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, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(f_x));
        }
 
        return iter_result<real>(x, iter, abs(f_x));
    }
 
 
    inline iter_result<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, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(s.Re()));
        }
 
        return iter_result<real>(x, iter, abs(s.Re()));
    }
 
 
    template<typename RealFunction>
    inline iter_result<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, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(f_x));
        }
 
        return iter_result<real>(x, iter, abs(f_x));
    }
 
 
    inline iter_result<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, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(s.Re()));
        }
 
        return iter_result<real>(x, iter, abs(s.Re()));
    }
 
 
    template<typename RealFunction>
    inline iter_result<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, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(f_x));
        }
 
        return iter_result<real>(x, iter, abs(f_x));
    }
 
 
    template<typename RealFunction>
    inline iter_result<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, MathError::NoConvergence);
            return iter_result<real>(ConvergenceStatus::MaxIterations, iter, abs(f_x));
        }
 
        return iter_result<real>(x, iter, abs(f_x));
    }
 
 
    template<typename RealFunction>
    inline iter_result<real> root(
        RealFunction f, real a, real b,
        real tol = OPTIMIZATION_TOL,
        unsigned int max_iter = OPTIMIZATION_BISECTION_ITER) {
        
        return root_itp(f, a, b, tol, max_iter);
    }
 
 
    template<typename RealFunction, typename Vector = vec<real>>
    inline Vector roots(
        RealFunction f, real a, real b,
        real tol = OPTIMIZATION_TOL, real steps = 10) {
 
        if(steps == 0) {
            TH_MATH_ERROR("roots", steps, MathError::DivByZero);
            return {nan()};
        }
 
        // Find candidate intervals
        auto intervals = find_root_brackets(f, a, b, steps);
 
        Vector res;
        res.resize(intervals.size(), nan());
        size_t idx = 0;
 
        // 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[idx] = intervals[i][0];
                idx++;
                continue;
            }
 
            if(abs(f(intervals[i][1])) <= MACH_EPSILON) {
                res[idx] = intervals[i][1];
                idx++;
                continue;
            }
 
            // Approximate the roots using bisection inside each interval
            real r = root_bisect(f, intervals[i][0], intervals[i][1], tol);
 
            if (!is_nan(r)) {
                res[idx] = r;
                idx++;
            }
        }
 
        return res;
    }
 
 
    template<typename Field, typename Vector = vec<Field>>
    inline Vector roots_polynomial(
        const polynomial<Field>& p, real tolerance = OPTIMIZATION_TOL,
        unsigned int steps = 0) {
        
        // Real polynomial degree
        const unsigned int p_order = p.find_order();
        if (p_order == 0) {
            TH_MATH_ERROR("roots_polynomial", p.coeff[0], MathError::InvalidArgument);
            return {nan()};
        }
 
        // Normalize the polynomial by the leading coefficient to apply Cauchy's bound
        polynomial<Field> p_norm = p / p.coeff[p_order];
        Field a_max = 0;
 
        // Find the maximum absolute coefficient
        for (unsigned int i = 0; i < p_order - 1; ++i)
            a_max = max(a_max, abs(p_norm.coeff[i]));
 
        // The roots are bounded in absolute value by the maximum
        const Field M = 1 + a_max;
 
        return roots(p_norm, -M, M, tolerance, steps > 0 ? steps : 10 * p_order);
    }
 
}
 
#endif