polynomial.hpp

// Copyright Take Vos 2022.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt)
 
#pragma once
 
#include "../utility/utility.hpp"
#include "../container/container.hpp"
#include "../macros.hpp"
#include <numbers>
#include <array>
 
hi_export_module(hikogui.numeric.polynomial);
 
hi_export namespace hi::inline v1 {
 
template<typename T>
hi_force_inline constexpr lean_vector<T> solvePolynomial(T const& a, T const& b) noexcept
{
    if (a != 0) {
        return {static_cast<T>(-(b / a))};
    } else if (b == 0) {
        // Any value of x is correct.
        return {T{0}};
    } else {
        // None of the values for x is correct.
        return {};
    }
}
 
template<typename T>
hi_force_inline constexpr lean_vector<T> solvePolynomial(T const& a, T const& b, T const& c) noexcept
{
    if (a == 0) {
        return solvePolynomial(b, c);
    } else {
        auto const D = b * b - T{4} * a * c;
        if (D < 0) {
            return {};
        } else if (D == 0) {
            return {static_cast<T>(-b / (T{2} * a))};
        } else {
            auto const Dsqrt = sqrt(D);
            return {static_cast<T>((-b - Dsqrt) / (T{2} * a)), static_cast<T>((-b + Dsqrt) / (T{2} * a))};
        }
    }
}
 
template<typename T>
hi_force_inline lean_vector<T> solveDepressedCubicTrig(T const& p, T const& q) noexcept
{
    constexpr T oneThird = T{1} / T{3};
    constexpr T pi2_3 = (T{2} / T{3}) * std::numbers::pi_v<T>;
    constexpr T pi4_3 = (T{4} / T{3}) * std::numbers::pi_v<T>;
 
    auto const U = oneThird * acos(((T{3} * q) / (T{2} * p)) * sqrt(T{-3} / p));
    auto const V = T{2} * sqrt(-oneThird * p);
 
    auto const t0 = V * cos(U);
    auto const t1 = V * cos(U - pi2_3);
    auto const t2 = V * cos(U - pi4_3);
    return {static_cast<T>(t0), static_cast<T>(t1), static_cast<T>(t2)};
}
 
template<typename T>
hi_force_inline lean_vector<T> solveDepressedCubicCardano(T const& p, T const& q, T const& D) noexcept
{
    auto const sqrtD = sqrt(D);
    auto const minusHalfQ = T{-0.5} * q;
    auto const v = cbrt(minusHalfQ + sqrtD);
    auto const w = cbrt(minusHalfQ - sqrtD);
    return {static_cast<T>(v + w)};
}
 
template<typename T>
hi_force_inline lean_vector<T> solveDepressedCubic(T const& p, T const& q) noexcept
{
    constexpr T oneForth = T{1} / T{4};
    constexpr T oneTwentySeventh = T{1} / T{27};
 
    if (p == 0.0 && q == 0.0) {
        return {T{0}};
 
    } else {
        auto const D = oneForth * q * q + oneTwentySeventh * p * p * p;
 
        if (D < 0 && p != 0.0) {
            // Has three real roots.
            return solveDepressedCubicTrig(p, q);
 
        } else if (D == 0 && p != 0.0) {
            // Has two real roots, or maybe one root
            auto const t0 = (T{3} * q) / p;
            auto const t1 = (T{-3} * q) / (T{2} * p);
            return {static_cast<T>(t0), static_cast<T>(t1), static_cast<T>(t1)};
 
        } else {
            // Has one real root.
            return solveDepressedCubicCardano(p, q, D);
        }
    }
}
 
template<typename T>
hi_force_inline constexpr lean_vector<T> solvePolynomial(T const& a, T const& b, T const& c, T const& d) noexcept
{
    if (a == 0) {
        return solvePolynomial(b, c, d);
 
    } else {
        auto const p = (T{3} * a * c - b * b) / (T{3} * a * a);
        auto const q = (T{2} * b * b * b - T{9} * a * b * c + T{27} * a * a * d) /
            (T{27} * a * a * a);
 
        auto const b_3a = b / (T{3} * a);
 
        auto r = solveDepressedCubic(p, q);
        for (auto &x: r) {
            x -= b_3a;
        }
        
        return r;
    }
}
 
} // namespace hi::inline v1