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

// Polynomial.hpp
// 多項式テンプレートクラス
// lib++20 多項式.h からの移植 (2026-02-15)
//
// 内部表現: 昇べきの順 (c[0]=定数項, c[n]=主係数)
// 設計:
//   - T = double, float, Int, Float, Rational, Complex<T> 等で使用可能
//   - Horner 法による O(n) 評価
//   - 多項式除算 (商と剰余)
//   - ユークリッド GCD
//   - 微分・積分
//   - 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 {

// ================================================================
// コンセプト
// ================================================================

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> クラス
// ================================================================

template<PolynomialCoeff T>
class Polynomial {
private:
    // 係数配列 (昇べきの順): coeffs_[i] は x^i の係数
    // 空 = 零多項式
    std::vector<T> coeffs_;

    // 先頭のゼロ係数を除去して正規化
    void normalize() {
        while (!coeffs_.empty() && coeffs_.back() == T(0))
            coeffs_.pop_back();
    }

public:
    // ============================================================
    // コンストラクタ
    // ============================================================

    // 零多項式
    Polynomial() = default;

    // 定数多項式
    Polynomial(const T& c) {
        if (!(c == T(0))) coeffs_ = {c};
    }

    // 係数ベクトルから (昇べきの順)
    explicit Polynomial(std::vector<T> coeffs) : coeffs_(std::move(coeffs)) {
        normalize();
    }

    // 初期化子リストから (昇べきの順: coeffs_[i] = x^i の係数)
    // 例: {1, 2, 3} → 1 + 2x + 3x^2
    Polynomial(std::initializer_list<T> il) {
        coeffs_.assign(il.begin(), il.end());
        normalize();
    }

    // ============================================================
    // アクセサ
    // ============================================================

    // 次数 (-1 = 零多項式)
    [[nodiscard]] int degree() const {
        return static_cast<int>(coeffs_.size()) - 1;
    }

    // 零多項式か
    [[nodiscard]] bool isZero() const { return coeffs_.empty(); }

    // 主係数 (零多項式なら T(0))
    [[nodiscard]] T leadingCoefficient() const {
        return coeffs_.empty() ? T(0) : coeffs_.back();
    }

    // 定数項
    [[nodiscard]] T constantTerm() const {
        return coeffs_.empty() ? T(0) : coeffs_[0];
    }

    // 係数アクセス (範囲外は T(0))
    [[nodiscard]] const T& operator[](size_t i) const {
        static const T zero = T(0);
        return (i < coeffs_.size()) ? coeffs_[i] : zero;
    }

    // 係数への参照アクセス (範囲外は自動拡張)
    T& operator[](size_t i) {
        if (i >= coeffs_.size())
            coeffs_.resize(i + 1, T(0));
        return coeffs_[i];
    }

    // 係数ベクトルへのアクセス
    [[nodiscard]] const std::vector<T>& coefficients() const { return coeffs_; }

    // ============================================================
    // 比較演算子
    // ============================================================

    [[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);
    }

    // ============================================================
    // 単項演算子
    // ============================================================

    [[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;
    }

    // ============================================================
    // 多項式 + 多項式
    // ============================================================

    [[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;
    }

    // 多項式乗算 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;
    }

    // ============================================================
    // 多項式除算: 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)};
    }

    // 多項式除算の商 (ユークリッド除算)
    // SANGI_POLY_DIV_RATIONAL 定義時は operator/ が RationalFunction を返すため、
    // 明示的に商が必要な場合はこちらを使用する。
    [[nodiscard]] Polynomial quo(const Polynomial& rhs) const {
        return divmod(rhs).quotient;
    }

#ifndef SANGI_POLY_DIV_RATIONAL
    // デフォルト: operator/ は多項式除算の商を返す (int / int と同じ意味)
    // SANGI_POLY_DIV_RATIONAL を定義すると RationalFunction を返すようになる
    // (定義は 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;
    }

    // ============================================================
    // 多項式 + スカラー / スカラー + 多項式
    // ============================================================

    [[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;
    }

    // ============================================================
    // 複合代入演算子
    // ============================================================

    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; }

    // ============================================================
    // 冪乗
    // ============================================================

    [[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;
    }

    // ============================================================
    // 評価 (Horner 法)
    // ============================================================

    // P(x) を Horner 法で計算: 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 法 (SIMD / ILP 向け多項式評価)
    // ツリーベースの並列評価: Horner と同じ結果を O(n) で計算するが
    // 乗算の依存チェーンが O(log n) に短縮されるため ILP が高い。
    //
    // 例: a0 + a1·x + a2·x² + a3·x³
    //   = (a0 + a1·x) + x²·(a2 + a3·x)
    //   level 0: p[i] = a[2i] + a[2i+1]·x     (並列)
    //   level 1: p[i] = p[2i] + p[2i+1]·x²     (並列)
    //   ...
    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]);

        // 作業配列: 初期値は係数のコピー
        std::vector<U> buf(n);
        for (size_t i = 0; i < n; ++i)
            buf[i] = U(coeffs_[i]);

        // x のべき乗を段階的に二乗
        U xpow = x;  // level 0: x, level 1: x², level 2: x⁴, ...

        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;
            // n が奇数の場合、最後の要素をそのまま残す
            if (n & 1)
                buf[half] = buf[n - 1];
            n = half + (n & 1);
            xpow = xpow * xpow;
        }

        return buf[0];
    }

    // ============================================================
    // 微分・積分
    // ============================================================

    // n 階微分 (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;
    }

    // 不定積分 (定数項 = 0)
    // 注: T が除算をサポートしている必要がある
    [[nodiscard]] Polynomial integral() const {
        if (coeffs_.empty()) return Polynomial();
        Polynomial result;
        result.coeffs_.resize(coeffs_.size() + 1, T(0));
        result.coeffs_[0] = T(0);  // 積分定数 = 0
        for (size_t i = 0; i < coeffs_.size(); ++i)
            result.coeffs_[i + 1] = coeffs_[i] / T(static_cast<int>(i + 1));
        return result;
    }

    // ============================================================
    // 文字列変換・出力
    // ============================================================

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

        std::string s;
        bool first = true;

        // 降べきの順で出力
        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();
    }
};

// ================================================================
// 自由関数
// ================================================================

// --- 多項式の構築ヘルパー ---

// 単項式 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));
}

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

// --- GCD ---
// 整数型: 擬剰余 (pseudo-remainder) + primitive part で係数爆発を抑制
// 浮動小数点型: ユークリッド互除法 (divmod で商が正確に計算できる)

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>) {
        // 整数型: 擬剰余方式 (divmod の整数除算で剰余が減らない問題を回避)
        // primitive part で係数を縮小
        a = primitivePart(a);
        b = primitivePart(b);
        if (a.degree() < b.degree()) std::swap(a, b);

        while (!b.isZero()) {
            // 擬剰余: lc(b)^δ · a を b で割った剰余 (δ = 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 {
        // 浮動小数点型: 通常のユークリッド互除法
        while (!b.isZero()) {
            Polynomial<T> r = a % b;
            a = std::move(b);
            b = std::move(r);
        }
        // モニック化 (主係数を 1 にする)
        if (!a.isZero()) {
            T lc = a.leadingCoefficient();
            if (!(lc == T(1)))
                a = a / lc;
        }
        return a;
    }
}

// --- 近似 GCD (Approximate GCD) ---
// 係数に数値誤差がある多項式の GCD を求める。
// ε-ユークリッド法: 剰余の係数の最大絶対値が tolerance 以下になったら
// その手前の剰余を GCD とみなす。

/**
 * @brief 近似 GCD (ε-ユークリッド法)
 *
 * 厳密な GCD は係数の微小な誤差で次数が変わるが、
 * この関数は剰余のノルムが tolerance 以下になった時点で打ち切る。
 * 信号処理、制御理論、数値代数で使用される。
 *
 * @param a         入力多項式
 * @param b         入力多項式
 * @param tolerance 剰余の係数最大絶対値がこれ以下なら零とみなす
 * @return 近似的な GCD (モニック化済み)
 */
template<PolynomialCoeff T>
[[nodiscard]] Polynomial<T> approximateGCD(
    Polynomial<T> a, Polynomial<T> b,
    double tolerance = 1e-10) {

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

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

        // 剰余の係数の最大絶対値を計算
        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;
        }

        // 十分小さければ打ち切り
        if (max_coeff <= tolerance) break;

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

    // モニック化
    if (!b.isZero()) {
        // b が非零のまま打ち切った場合は b の手前 = a が GCD 候補だが
        // ループ構造上、打ち切り時の a は直前の b
        // ここでは a が最後の非零剰余
    }
    if (!a.isZero()) {
        T lc = a.leadingCoefficient();
        if (!(lc == T(1)))
            a = a / lc;
    }
    return a;
}

// --- 内容 (content): 全係数の GCD ---
// 整数型多項式のみ。GCD は 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;
}

// --- 定積分 ---

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);
}

// --- 合成 p(q(x)) ---

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

// ================================================================
// 特殊多項式
// ================================================================

// --- Chebyshev 多項式 (第 1 種) ---
// 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 多項式 (第 2 種) ---
// 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 多項式 ---
// 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 多項式 (確率論版: He) ---
// H_0 = 1, H_1 = x, H_n = x*H_{n-1} - (n-1)*H_{n-2}
// 注: 物理版 (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 (物理版)

    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}  (物理版)
        Polynomial<T> curr = twoX * prev1 - prev2 * T(2 * (i - 1));
        prev2 = std::move(prev1);
        prev1 = std::move(curr);
    }
    return prev1;
}

// --- Laguerre 多項式 ---
// 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 多項式 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)
/// Bessel フィルタ設計、信号処理に使用
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 多項式 P_n^{(α,β)}(x)
/// P_0 = 1, P_1 = (α-β)/2 + (α+β+2)/2 · x
/// (2n+α+β)(2n+α+β-1)/(2n(n+α+β)) · ((2n+α+β-1)·x + (α²-β²)/(2n+α+β-2)) · P_{n-1}
///   - (n+α-1)(n+β-1)(2n+α+β) / (n(n+α+β)(2n+α+β-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 (超球面) 多項式 C_n^{(λ)}(x)
/// C_0 = 1, C_1 = 2λx
/// C_n = (2(n+λ-1)·x·C_{n-1} - (n+2λ-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;
}

// ================================================================
// 型特性
// ================================================================

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