legendre.h

#pragma once
#include <array>

namespace Faunus {
template <std::floating_point T, std::size_t max_order, bool use_table = false> class Legendre
{
  private:
    std::array<T, max_order + 1> y; 
    std::array<T, max_order + 1> P; 
  public:
    Legendre()
    {
        P[0] = 1.0;
        if constexpr (use_table) {
            for (std::size_t i = 1; i < max_order; ++i) {
                y[i] = 1.0 + 1.0 / T(i);
            }
        }
    }

    const auto& eval(T x)
    {
        if constexpr (max_order > 0) {
            P[1] = x;
            for (std::size_t i = 1; i < max_order; ++i) {
                if constexpr (use_table) {
                    P[i + 1] = ((y[i] + 1.0) * x * P[i] - P[i - 1]) / y[i];
                }
                else {
                    P[i + 1] = ((2.0 + 1.0 / T(i)) * x * P[i] - P[i - 1]) / (1.0 + 1.0 / T(i));
                }
            }
        }
        return P;
    }
};
#ifdef DOCTEST_LIBRARY_INCLUDED__
TEST_CASE_TEMPLATE("[Faunus] Legendre", LegendreType, Legendre<double, 3, false>,
                   Legendre<double, 3, true>)
{
    using doctest::Approx;
    LegendreType l;
    double x = 2.2;
    auto P = l.eval(x);
    CHECK_EQ(P[0], Approx(1));
    CHECK_EQ(P[1], Approx(x));
    CHECK_EQ(P[2], Approx(0.5 * (3 * x * x - 1)));
    CHECK_EQ(P[3], Approx(0.5 * (5 * x * x * x - 3 * x)));
}
#endif

} // namespace Faunus