// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// Polynomial.hpp
// Polynomial template class
//
// Internal representation: ascending order (c[0] = constant term, c[n] = leading coefficient)
// Design:
//   - Usable with T = double, float, Int, Float, Rational, Complex<T>, etc.
//   - O(n) evaluation via Horner's method
//   - Polynomial division (quotient and remainder)
//   - Euclidean GCD
//   - Differentiation / integration
//   - Special polynomials: Chebyshev, Legendre, Hermite, Laguerre

#pragma once

#include <vector>
#include <numeric>
#include <iosfwd>
#include <sstream>
#include <string>
#include <cmath>
#include <concepts>
#include <algorithm>
#include <initializer_list>
#include <cassert>

namespace sangi {

// ================================================================
// Concepts
// ================================================================

template<typename T>
concept PolynomialCoeff = requires(T a, T b) {
    { a + b } -> std::convertible_to<T>;
    { a - b } -> std::convertible_to<T>;
    { a * b } -> std::convertible_to<T>;
    { -a } -> std::convertible_to<T>;
    { T(0) };
    { T(1) };
};

// ================================================================
// Polynomial<T> class
// ================================================================

template<PolynomialCoeff T>
class Polynomial {
private:
    // Coefficient array (ascending order): coeffs_[i] is the coefficient of x^i.
    // Empty = zero polynomial.
    std::vector<T> coeffs_;

    // Normalize by trimming trailing zero coefficients
    void normalize() {
        while (!coeffs_.empty() && coeffs_.back() == T(0))
            coeffs_.pop_back();
    }

public:
    // ============================================================
    // Constructors
    // ============================================================

    // Zero polynomial
    Polynomial() = default;

    // Constant polynomial
    Polynomial(const T& c) {
        if (!(c == T(0))) coeffs_ = {c};
    }

    // From a coefficient vector (ascending order)
    explicit Polynomial(std::vector<T> coeffs) : coeffs_(std::move(coeffs)) {
        normalize();
    }

    // From an initializer list (ascending order: coeffs_[i] = coefficient of x^i)
    // Example: {1, 2, 3} -> 1 + 2x + 3x^2
    Polynomial(std::initializer_list<T> il) {
        coeffs_.assign(il.begin(), il.end());
        normalize();
    }

    // ============================================================
    // Accessors
    // ============================================================

    // Degree (-1 means the zero polynomial)
    [[nodiscard]] int degree() const {
        return static_cast<int>(coeffs_.size()) - 1;
    }

    // Whether this is the zero polynomial
    [[nodiscard]] bool isZero() const { return coeffs_.empty(); }

    // Leading coefficient (T(0) if zero polynomial)
    [[nodiscard]] T leadingCoefficient() const {
        return coeffs_.empty() ? T(0) : coeffs_.back();
    }

    // Constant term
    [[nodiscard]] T constantTerm() const {
        return coeffs_.empty() ? T(0) : coeffs_[0];
    }

    // Coefficient access (out of range returns T(0))
    [[nodiscard]] const T& operator[](size_t i) const {
        static const T zero = T(0);
        return (i < coeffs_.size()) ? coeffs_[i] : zero;
    }

    // Coefficient reference access (auto-extends if out of range)
    T& operator[](size_t i) {
        if (i >= coeffs_.size())
            coeffs_.resize(i + 1, T(0));
        return coeffs_[i];
    }

    // Access to the coefficient vector
    [[nodiscard]] const std::vector<T>& coefficients() const { return coeffs_; }

    // ============================================================
    // Comparison operators
    // ============================================================

    [[nodiscard]] bool operator==(const Polynomial& rhs) const {
        return coeffs_ == rhs.coeffs_;
    }
    [[nodiscard]] bool operator!=(const Polynomial& rhs) const {
        return !(*this == rhs);
    }
    [[nodiscard]] bool operator==(const T& rhs) const {
        if (rhs == T(0)) return coeffs_.empty();
        return coeffs_.size() == 1 && coeffs_[0] == rhs;
    }
    [[nodiscard]] bool operator!=(const T& rhs) const {
        return !(*this == rhs);
    }

    // ============================================================
    // Unary operators
    // ============================================================

    [[nodiscard]] Polynomial operator+() const { return *this; }

    [[nodiscard]] Polynomial operator-() const {
        Polynomial result;
        result.coeffs_.resize(coeffs_.size());
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result.coeffs_[i] = -coeffs_[i];
        return result;
    }

    // ============================================================
    // Polynomial + polynomial
    // ============================================================

    [[nodiscard]] Polynomial operator+(const Polynomial& rhs) const {
        size_t n = std::max(coeffs_.size(), rhs.coeffs_.size());
        Polynomial result;
        result.coeffs_.resize(n, T(0));
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result.coeffs_[i] = result.coeffs_[i] + coeffs_[i];
        for (size_t i = 0; i < rhs.coeffs_.size(); ++i)
            result.coeffs_[i] = result.coeffs_[i] + rhs.coeffs_[i];
        result.normalize();
        return result;
    }

    [[nodiscard]] Polynomial operator-(const Polynomial& rhs) const {
        size_t n = std::max(coeffs_.size(), rhs.coeffs_.size());
        Polynomial result;
        result.coeffs_.resize(n, T(0));
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result.coeffs_[i] = result.coeffs_[i] + coeffs_[i];
        for (size_t i = 0; i < rhs.coeffs_.size(); ++i)
            result.coeffs_[i] = result.coeffs_[i] - rhs.coeffs_[i];
        result.normalize();
        return result;
    }

    // Polynomial multiplication O(n*m)
    [[nodiscard]] Polynomial operator*(const Polynomial& rhs) const {
        if (coeffs_.empty() || rhs.coeffs_.empty())
            return Polynomial();

        size_t n = coeffs_.size() + rhs.coeffs_.size() - 1;
        Polynomial result;
        result.coeffs_.resize(n, T(0));
        for (size_t i = 0; i < coeffs_.size(); ++i)
            for (size_t j = 0; j < rhs.coeffs_.size(); ++j)
                result.coeffs_[i + j] = result.coeffs_[i + j] + coeffs_[i] * rhs.coeffs_[j];
        result.normalize();
        return result;
    }

    // ============================================================
    // Polynomial division: this = quotient * rhs + remainder
    // ============================================================

    struct DivResult {
        Polynomial quotient;
        Polynomial remainder;
    };

    [[nodiscard]] DivResult divmod(const Polynomial& divisor) const {
        assert(!divisor.isZero() && "Division by zero polynomial");

        if (degree() < divisor.degree())
            return DivResult{Polynomial(), *this};

        Polynomial rem = *this;
        int qDeg = degree() - divisor.degree();
        Polynomial quot;
        quot.coeffs_.resize(qDeg + 1, T(0));

        T lcDiv = divisor.leadingCoefficient();
        int divDeg = divisor.degree();

        for (int i = qDeg; i >= 0; --i) {
            T coeff = rem.coeffs_[i + divDeg] / lcDiv;
            quot.coeffs_[i] = coeff;
            for (int j = 0; j <= divDeg; ++j)
                rem.coeffs_[i + j] = rem.coeffs_[i + j] - coeff * divisor.coeffs_[j];
        }

        rem.normalize();
        quot.normalize();
        return DivResult{std::move(quot), std::move(rem)};
    }

    // Quotient of polynomial division (Euclidean division).
    // When SANGI_POLY_DIV_RATIONAL is defined, operator/ returns a
    // RationalFunction, so use this function when the quotient is explicitly
    // required.
    [[nodiscard]] Polynomial quo(const Polynomial& rhs) const {
        return divmod(rhs).quotient;
    }

#ifndef SANGI_POLY_DIV_RATIONAL
    // Default: operator/ returns the polynomial quotient (same meaning as int / int).
    // Defining SANGI_POLY_DIV_RATIONAL makes it return a RationalFunction
    // instead (the definition is a free function in RationalFunction.hpp).
    [[nodiscard]] Polynomial operator/(const Polynomial& rhs) const {
        return divmod(rhs).quotient;
    }
#endif

    [[nodiscard]] Polynomial operator%(const Polynomial& rhs) const {
        return divmod(rhs).remainder;
    }

    // ============================================================
    // Polynomial + scalar / scalar + polynomial
    // ============================================================

    [[nodiscard]] Polynomial operator+(const T& c) const {
        Polynomial result = *this;
        if (result.coeffs_.empty())
            result.coeffs_.push_back(c);
        else
            result.coeffs_[0] = result.coeffs_[0] + c;
        result.normalize();
        return result;
    }

    [[nodiscard]] Polynomial operator-(const T& c) const {
        return *this + (-c);
    }

    [[nodiscard]] Polynomial operator*(const T& c) const {
        if (c == T(0)) return Polynomial();
        Polynomial result;
        result.coeffs_.resize(coeffs_.size());
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result.coeffs_[i] = coeffs_[i] * c;
        result.normalize();
        return result;
    }

    [[nodiscard]] Polynomial operator/(const T& c) const {
        Polynomial result;
        result.coeffs_.resize(coeffs_.size());
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result.coeffs_[i] = coeffs_[i] / c;
        result.normalize();
        return result;
    }

    friend Polynomial operator+(const T& lhs, const Polynomial& rhs) {
        return rhs + lhs;
    }
    friend Polynomial operator-(const T& lhs, const Polynomial& rhs) {
        return (-rhs) + lhs;
    }
    friend Polynomial operator*(const T& lhs, const Polynomial& rhs) {
        return rhs * lhs;
    }

    // ============================================================
    // Compound assignment operators
    // ============================================================

    Polynomial& operator+=(const Polynomial& rhs) { *this = *this + rhs; return *this; }
    Polynomial& operator-=(const Polynomial& rhs) { *this = *this - rhs; return *this; }
    Polynomial& operator*=(const Polynomial& rhs) { *this = *this * rhs; return *this; }
    Polynomial& operator/=(const Polynomial& rhs) { *this = *this / rhs; return *this; }
    Polynomial& operator%=(const Polynomial& rhs) { *this = *this % rhs; return *this; }
    Polynomial& operator+=(const T& c) { *this = *this + c; return *this; }
    Polynomial& operator-=(const T& c) { *this = *this - c; return *this; }
    Polynomial& operator*=(const T& c) { *this = *this * c; return *this; }
    Polynomial& operator/=(const T& c) { *this = *this / c; return *this; }

    // ============================================================
    // Power
    // ============================================================

    [[nodiscard]] Polynomial pow(int n) const {
        assert(n >= 0 && "Negative exponent not supported");
        if (n == 0) return Polynomial(T(1));
        Polynomial result(T(1));
        Polynomial base = *this;
        unsigned int m = static_cast<unsigned int>(n);
        while (m > 0) {
            if (m & 1) result *= base;
            base *= base;
            m >>= 1;
        }
        return result;
    }

    // ============================================================
    // Evaluation (Horner's method)
    // ============================================================

    // Compute P(x) by Horner's method: O(n)
    template<typename U = T>
    [[nodiscard]] U operator()(const U& x) const {
        if (coeffs_.empty()) return U(0);

        U result = U(coeffs_.back());
        for (int i = static_cast<int>(coeffs_.size()) - 2; i >= 0; --i)
            result = result * x + U(coeffs_[i]);
        return result;
    }

    // Estrin's scheme (polynomial evaluation suited to SIMD / ILP)
    // Tree-based parallel evaluation: same result as Horner in O(n), but the
    // multiplication dependency chain is reduced to O(log n), giving higher ILP.
    //
    // Example: a0 + a1*x + a2*x^2 + a3*x^3
    //   = (a0 + a1*x) + x^2 * (a2 + a3*x)
    //   level 0: p[i] = a[2i] + a[2i+1]*x       (in parallel)
    //   level 1: p[i] = p[2i] + p[2i+1]*x^2     (in parallel)
    //   ...
    template<typename U = T>
    [[nodiscard]] U evaluateEstrin(const U& x) const {
        size_t n = coeffs_.size();
        if (n == 0) return U(0);
        if (n == 1) return U(coeffs_[0]);

        // Work array: initially a copy of the coefficients
        std::vector<U> buf(n);
        for (size_t i = 0; i < n; ++i)
            buf[i] = U(coeffs_[i]);

        // Square x progressively
        U xpow = x;  // level 0: x, level 1: x^2, level 2: x^4, ...

        while (n > 1) {
            size_t half = n / 2;
            for (size_t i = 0; i < half; ++i)
                buf[i] = buf[2 * i] + buf[2 * i + 1] * xpow;
            // When n is odd, leave the last element as-is
            if (n & 1)
                buf[half] = buf[n - 1];
            n = half + (n & 1);
            xpow = xpow * xpow;
        }

        return buf[0];
    }

    // ============================================================
    // Differentiation / integration
    // ============================================================

    // n-th derivative (n > 0)
    [[nodiscard]] Polynomial derivative(int n = 1) const {
        assert(n >= 0);
        Polynomial result = *this;
        for (int k = 0; k < n; ++k) {
            if (result.coeffs_.size() <= 1) return Polynomial();
            Polynomial diff;
            diff.coeffs_.resize(result.coeffs_.size() - 1);
            for (size_t i = 1; i < result.coeffs_.size(); ++i)
                diff.coeffs_[i - 1] = result.coeffs_[i] * T(static_cast<int>(i));
            result = std::move(diff);
        }
        return result;
    }

    // Indefinite integral (constant of integration = 0)
    // Note: T must support division.
    [[nodiscard]] Polynomial integral() const {
        if (coeffs_.empty()) return Polynomial();
        Polynomial result;
        result.coeffs_.resize(coeffs_.size() + 1, T(0));
        result.coeffs_[0] = T(0);  // Constant of integration = 0
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result.coeffs_[i + 1] = coeffs_[i] / T(static_cast<int>(i + 1));
        return result;
    }

    // ============================================================
    // String conversion / output
    // ============================================================

    [[nodiscard]] std::string toString(const std::string& var = "x") const {
        if (coeffs_.empty()) return "0";

        std::string s;
        bool first = true;

        // Output in descending order
        for (int i = degree(); i >= 0; --i) {
            if (coeffs_[i] == T(0)) continue;

            std::ostringstream oss;
            oss << coeffs_[i];
            std::string coefStr = oss.str();

            if (!first && !coefStr.empty() && coefStr[0] != '-')
                s += "+";

            if (i == 0) {
                s += coefStr;
            } else if (i == 1) {
                if (coefStr == "1") s += var;
                else if (coefStr == "-1") s += "-" + var;
                else s += coefStr + var;
            } else {
                if (coefStr == "1") s += var + "^" + std::to_string(i);
                else if (coefStr == "-1") s += "-" + var + "^" + std::to_string(i);
                else s += coefStr + var + "^" + std::to_string(i);
            }
            first = false;
        }

        return s.empty() ? "0" : s;
    }

    friend std::ostream& operator<<(std::ostream& os, const Polynomial& p) {
        return os << p.toString();
    }
};

// ================================================================
// Free functions
// ================================================================

// --- Helpers for constructing polynomials ---

// Monomial c * x^n
template<PolynomialCoeff T>
[[nodiscard]] Polynomial<T> monomial(const T& c, int n) {
    assert(n >= 0);
    std::vector<T> coeffs(n + 1, T(0));
    coeffs[n] = c;
    return Polynomial<T>(std::move(coeffs));
}

// Variable x (= 1 * x^1)
template<PolynomialCoeff T = double>
[[nodiscard]] Polynomial<T> X() {
    return monomial(T(1), 1);
}

// --- GCD ---
// Integer types: use pseudo-remainder + primitive part to suppress coefficient explosion.
// Floating-point types: use the Euclidean algorithm (divmod computes the quotient exactly).

template<PolynomialCoeff T>
[[nodiscard]] Polynomial<T> gcd(Polynomial<T> a, Polynomial<T> b) {
    if (a.isZero()) return b;
    if (b.isZero()) return a;

    if constexpr (std::is_integral_v<T>) {
        // Integer types: pseudo-remainder approach (avoids the problem that
        // integer divmod fails to reduce the remainder).
        // Use the primitive part to shrink coefficients.
        a = primitivePart(a);
        b = primitivePart(b);
        if (a.degree() < b.degree()) std::swap(a, b);

        while (!b.isZero()) {
            // Pseudo-remainder: remainder of (lc(b)^delta * a) divided by b
            // (delta = deg(a) - deg(b) + 1)
            int delta = a.degree() - b.degree() + 1;
            T lc_b = b.leadingCoefficient();
            T scale = T(1);
            for (int i = 0; i < delta; ++i) scale *= lc_b;
            Polynomial<T> r = (a * scale).divmod(b).remainder;
            if (r.isZero()) break;
            r = primitivePart(r);
            a = std::move(b);
            b = std::move(r);
        }
        auto& result = b.isZero() ? a : b;
        if (result.leadingCoefficient() < T(0)) result = -result;
        return result;
    } else {
        // Floating-point types: ordinary Euclidean algorithm
        while (!b.isZero()) {
            Polynomial<T> r = a % b;
            a = std::move(b);
            b = std::move(r);
        }
        // Make monic (set the leading coefficient to 1)
        if (!a.isZero()) {
            T lc = a.leadingCoefficient();
            if (!(lc == T(1)))
                a = a / lc;
        }
        return a;
    }
}

// --- Approximate GCD ---
// Compute the GCD of polynomials whose coefficients carry numerical error.
// Epsilon-Euclidean method: when the maximum absolute value of the remainder's
// coefficients falls below tolerance, treat the preceding remainder as the GCD.

/**
 * @brief Approximate GCD (epsilon-Euclidean method)
 *
 * Exact GCD's degree can change with tiny coefficient errors. This function
 * stops once the norm of the remainder falls at or below tolerance. Used in
 * signal processing, control theory, and numerical algebra.
 *
 * @param a         input polynomial
 * @param b         input polynomial
 * @param tolerance maximum absolute value of remainder coefficients below
 *                  which the remainder is considered zero
 * @return approximate GCD (monicized)
 */
template<PolynomialCoeff T>
[[nodiscard]] Polynomial<T> approximateGCD(
    Polynomial<T> a, Polynomial<T> b,
    double tolerance = 1e-10) {

    // Ensure |a| >= |b|
    if (a.degree() < b.degree()) std::swap(a, b);

    while (!b.isZero()) {
        Polynomial<T> r = a % b;

        // Compute the maximum absolute value of the remainder's coefficients
        double max_coeff = 0.0;
        for (int i = 0; i <= r.degree(); ++i) {
            double c = std::abs(static_cast<double>(r[i]));
            if (c > max_coeff) max_coeff = c;
        }

        // Stop if it is small enough
        if (max_coeff <= tolerance) break;

        a = std::move(b);
        b = std::move(r);
    }

    // Monicize
    if (!b.isZero()) {
        // If we stopped while b was still non-zero, the GCD candidate would
        // be the one before b (i.e. a), but due to the loop structure, a at
        // the stopping point is the previous b. Here a is the last non-zero
        // remainder.
    }
    if (!a.isZero()) {
        T lc = a.leadingCoefficient();
        if (!(lc == T(1)))
            a = a / lc;
    }
    return a;
}

// --- Content: GCD of all coefficients ---
// Integer polynomials only. Uses std::gcd.

template<PolynomialCoeff T>
[[nodiscard]] T content(const Polynomial<T>& p) {
    static_assert(std::is_integral_v<T>,
                  "content() requires integral coefficient type");
    if (p.isZero()) return T(0);
    T g = T(0);
    for (int i = 0; i <= p.degree(); ++i) {
        T c = p[i];
        if (c < T(0)) c = -c;
        if (g == T(0)) g = c;
        else if (c != T(0)) g = std::gcd(g, c);
    }
    return (g == T(0)) ? T(1) : g;
}

// --- Primitive part: p / content(p) ---

template<PolynomialCoeff T>
[[nodiscard]] Polynomial<T> primitivePart(const Polynomial<T>& p) {
    static_assert(std::is_integral_v<T>,
                  "primitivePart() requires integral coefficient type");
    T c = content(p);
    if (c == T(0) || c == T(1)) return p;
    return p / c;
}

// --- Definite integral ---

template<PolynomialCoeff T>
[[nodiscard]] T definiteIntegral(const Polynomial<T>& p, const T& from, const T& to) {
    Polynomial<T> F = p.integral();
    return F(to) - F(from);
}

// --- Composition p(q(x)) ---

template<PolynomialCoeff T>
[[nodiscard]] Polynomial<T> compose(const Polynomial<T>& p, const Polynomial<T>& q) {
    return p(q);
}

// ================================================================
// Special polynomials
// ================================================================

// --- Chebyshev polynomials of the first kind ---
// T_0 = 1, T_1 = x, T_n = 2x*T_{n-1} - T_{n-2}

template<PolynomialCoeff T = double>
[[nodiscard]] Polynomial<T> chebyshevT(int n) {
    assert(n >= 0);
    if (n == 0) return Polynomial<T>(T(1));
    if (n == 1) return Polynomial<T>({T(0), T(1)});  // x

    Polynomial<T> twoX({T(0), T(2)});  // 2x
    Polynomial<T> prev2(T(1));          // T_0
    Polynomial<T> prev1({T(0), T(1)}); // T_1
    for (int i = 2; i <= n; ++i) {
        Polynomial<T> curr = twoX * prev1 - prev2;
        prev2 = std::move(prev1);
        prev1 = std::move(curr);
    }
    return prev1;
}

// --- Chebyshev polynomials of the second kind ---
// U_0 = 1, U_1 = 2x, U_n = 2x*U_{n-1} - U_{n-2}

template<PolynomialCoeff T = double>
[[nodiscard]] Polynomial<T> chebyshevU(int n) {
    assert(n >= 0);
    if (n == 0) return Polynomial<T>(T(1));
    if (n == 1) return Polynomial<T>({T(0), T(2)});  // 2x

    Polynomial<T> twoX({T(0), T(2)});  // 2x
    Polynomial<T> prev2(T(1));          // U_0
    Polynomial<T> prev1({T(0), T(2)}); // U_1
    for (int i = 2; i <= n; ++i) {
        Polynomial<T> curr = twoX * prev1 - prev2;
        prev2 = std::move(prev1);
        prev1 = std::move(curr);
    }
    return prev1;
}

// --- Legendre polynomials ---
// P_0 = 1, P_1 = x, (n+1)P_{n+1} = (2n+1)x*P_n - n*P_{n-1}

template<PolynomialCoeff T = double>
[[nodiscard]] Polynomial<T> legendre(int n) {
    assert(n >= 0);
    if (n == 0) return Polynomial<T>(T(1));
    if (n == 1) return Polynomial<T>({T(0), T(1)});  // x

    Polynomial<T> xPoly({T(0), T(1)});  // x
    Polynomial<T> prev2(T(1));           // P_0
    Polynomial<T> prev1({T(0), T(1)});  // P_1
    for (int i = 2; i <= n; ++i) {
        // P_i = ((2i-1)*x*P_{i-1} - (i-1)*P_{i-2}) / i
        Polynomial<T> curr = (xPoly * prev1 * T(2 * i - 1) - prev2 * T(i - 1)) / T(i);
        prev2 = std::move(prev1);
        prev1 = std::move(curr);
    }
    return prev1;
}

// --- Hermite polynomials (probabilist's version: He) ---
// H_0 = 1, H_1 = x, H_n = x*H_{n-1} - (n-1)*H_{n-2}
// Note: this differs from the physicist's version (2x*H_{n-1} - 2(n-1)*H_{n-2}).

template<PolynomialCoeff T = double>
[[nodiscard]] Polynomial<T> hermite(int n) {
    assert(n >= 0);
    if (n == 0) return Polynomial<T>(T(1));
    if (n == 1) return Polynomial<T>({T(0), T(2)});  // 2x (physicist's version)

    Polynomial<T> twoX({T(0), T(2)});  // 2x
    Polynomial<T> prev2(T(1));          // H_0
    Polynomial<T> prev1({T(0), T(2)}); // H_1 = 2x
    for (int i = 2; i <= n; ++i) {
        // H_n = 2x*H_{n-1} - 2(n-1)*H_{n-2}  (physicist's version)
        Polynomial<T> curr = twoX * prev1 - prev2 * T(2 * (i - 1));
        prev2 = std::move(prev1);
        prev1 = std::move(curr);
    }
    return prev1;
}

// --- Laguerre polynomials ---
// L_0 = 1, L_1 = 1-x, L_n = ((2n-1-x)*L_{n-1} - (n-1)*L_{n-2}) / n

template<PolynomialCoeff T = double>
[[nodiscard]] Polynomial<T> laguerre(int n) {
    assert(n >= 0);
    if (n == 0) return Polynomial<T>(T(1));
    if (n == 1) return Polynomial<T>({T(1), T(-1)});  // 1-x

    Polynomial<T> xPoly({T(0), T(1)});  // x
    Polynomial<T> prev2(T(1));           // L_0
    Polynomial<T> prev1({T(1), T(-1)}); // L_1 = 1-x
    for (int i = 2; i <= n; ++i) {
        // L_i = ((2i-1-x)*L_{i-1} - (i-1)*L_{i-2}) / i
        Polynomial<T> curr = ((Polynomial<T>(T(2 * i - 1)) - xPoly) * prev1
                              - prev2 * T(i - 1)) / T(i);
        prev2 = std::move(prev1);
        prev1 = std::move(curr);
    }
    return prev1;
}

/// Bessel polynomial y_n(x)
/// y_0(x) = 1, y_1(x) = x + 1
/// y_n(x) = (2n-1) * x * y_{n-1}(x) + y_{n-2}(x)
/// Used in Bessel filter design and signal processing.
template<PolynomialCoeff T = double>
[[nodiscard]] Polynomial<T> bessel(int n) {
    assert(n >= 0);
    if (n == 0) return Polynomial<T>(T(1));
    if (n == 1) return Polynomial<T>({T(1), T(1)});  // 1 + x

    Polynomial<T> xPoly({T(0), T(1)});  // x
    Polynomial<T> prev2(T(1));           // y_0
    Polynomial<T> prev1({T(1), T(1)});   // y_1 = 1 + x
    for (int i = 2; i <= n; ++i) {
        // y_i = (2i-1)*x*y_{i-1} + y_{i-2}
        Polynomial<T> curr = xPoly * prev1 * T(2 * i - 1) + prev2;
        prev2 = std::move(prev1);
        prev1 = std::move(curr);
    }
    return prev1;
}

/// Jacobi polynomial P_n^{(alpha,beta)}(x)
/// P_0 = 1, P_1 = (alpha-beta)/2 + (alpha+beta+2)/2 * x
/// (2n+alpha+beta)(2n+alpha+beta-1)/(2n(n+alpha+beta))
///   * ((2n+alpha+beta-1)*x + (alpha^2-beta^2)/(2n+alpha+beta-2)) * P_{n-1}
///   - (n+alpha-1)(n+beta-1)(2n+alpha+beta) / (n(n+alpha+beta)(2n+alpha+beta-2)) * P_{n-2}
template<PolynomialCoeff T = double>
[[nodiscard]] Polynomial<T> jacobi(int n, T alpha, T beta) {
    assert(n >= 0);
    if (n == 0) return Polynomial<T>(T(1));
    T ab = alpha + beta;
    if (n == 1) return Polynomial<T>({(alpha - beta) / T(2), (ab + T(2)) / T(2)});

    Polynomial<T> xPoly({T(0), T(1)});
    Polynomial<T> prev2(T(1));
    Polynomial<T> prev1({(alpha - beta) / T(2), (ab + T(2)) / T(2)});
    for (int i = 2; i <= n; ++i) {
        T ni = T(i), ab2 = T(2) * ni + ab;
        T c1 = ab2 * (ab2 - T(1)) / (T(2) * ni * (ni + ab));
        T c2 = (alpha * alpha - beta * beta) / ((ab2 - T(2)) * ab2);
        T c3 = (ni + alpha - T(1)) * (ni + beta - T(1)) * ab2
             / (ni * (ni + ab) * (ab2 - T(2)));
        Polynomial<T> curr = (xPoly * c1 + Polynomial<T>(c1 * c2)) * prev1 - prev2 * c3;
        prev2 = std::move(prev1);
        prev1 = std::move(curr);
    }
    return prev1;
}

/// Gegenbauer (ultraspherical) polynomial C_n^{(lambda)}(x)
/// C_0 = 1, C_1 = 2*lambda*x
/// C_n = (2(n+lambda-1)*x*C_{n-1} - (n+2*lambda-2)*C_{n-2}) / n
template<PolynomialCoeff T = double>
[[nodiscard]] Polynomial<T> gegenbauer(int n, T lambda) {
    assert(n >= 0);
    if (n == 0) return Polynomial<T>(T(1));
    if (n == 1) return Polynomial<T>({T(0), T(2) * lambda});

    Polynomial<T> xPoly({T(0), T(1)});
    Polynomial<T> prev2(T(1));
    Polynomial<T> prev1({T(0), T(2) * lambda});
    for (int i = 2; i <= n; ++i) {
        Polynomial<T> curr = (xPoly * prev1 * T(2 * (i + lambda - 1))
                              - prev2 * T(i + 2 * lambda - 2)) / T(i);
        prev2 = std::move(prev1);
        prev1 = std::move(curr);
    }
    return prev1;
}

// ================================================================
// Type traits
// ================================================================

template<typename T>
struct is_polynomial : std::false_type {};

template<PolynomialCoeff T>
struct is_polynomial<Polynomial<T>> : std::true_type {};

template<typename T>
inline constexpr bool is_polynomial_v = is_polynomial<T>::value;

} // namespace sangi