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

// Polynomial_factorization.hpp
// 整数多項式の厳密な因数分解
// lib++20 多項式.因数分解.cpp からの移植 (2026-03-13)
//
// アルゴリズム:
//   1. content 抽出 (全係数の GCD)
//   2. Yun の無平方分解 (gcd(f, f') による重複因子分離)
//   3. 有理根定理による 1 次因子分離
//   4. Kronecker 法による高次因子分離 (Lagrange 補間)
//
// 使用法:
//   #include <math/core/Polynomial.hpp>
//   #include <math/core/Polynomial_factorization.hpp>
//   #include <math/core/FactoredInteger.hpp>  // 結果の型

#pragma once

#include <math/core/Polynomial.hpp>
#include <math/core/FactoredInteger.hpp>
#include <numeric>
#include <cmath>
#include <algorithm>
#include <set>

namespace calx {

// ================================================================
// 内部ヘルパー
// ================================================================

namespace detail {

// 整数の約数を列挙 (正の約数のみ、ソート済み)
inline std::vector<int> divisors(int n) {
    if (n < 0) n = -n;
    if (n == 0) return {};
    std::vector<int> result;
    for (int i = 1; i * i <= n; ++i) {
        if (n % i == 0) {
            result.push_back(i);
            if (i != n / i)
                result.push_back(n / i);
        }
    }
    std::sort(result.begin(), result.end());
    return result;
}

// 整数の約数を列挙 (int64_t 版)
inline std::vector<int64_t> divisors(int64_t n) {
    if (n < 0) n = -n;
    if (n == 0) return {};
    std::vector<int64_t> result;
    for (int64_t i = 1; i * i <= n; ++i) {
        if (n % i == 0) {
            result.push_back(i);
            if (i != n / i)
                result.push_back(n / i);
        }
    }
    std::sort(result.begin(), result.end());
    return result;
}

// Lagrange 補間で多項式の係数を復元する
// xs, ys: 評価点と値 (n+1 個 → 最大 n 次多項式)
template<typename T>
Polynomial<T> lagrangeInterpolationPoly(
    const std::vector<T>& xs,
    const std::vector<T>& ys)
{
    int n = static_cast<int>(xs.size());
    Polynomial<T> result(T(0));
    Polynomial<T> one(T(1));

    for (int i = 0; i < n; ++i) {
        if (ys[i] == T(0)) continue;

        // L_i(x) = Π_{j≠i} (x - x_j) / (x_i - x_j)
        Polynomial<T> basis = one;
        T denom = T(1);
        for (int j = 0; j < n; ++j) {
            if (j == i) continue;
            // (x - x_j)
            basis = basis * Polynomial<T>({T(0) - xs[j], T(1)});
            denom = denom * (xs[i] - xs[j]);
        }
        // ys[i] / denom * basis
        result = result + basis * (ys[i] / denom);
    }
    return result;
}

// 多項式が整数係数かチェック (T が整数型なら常に true)
template<typename T>
bool hasIntegerCoeffs(const Polynomial<T>& p) {
    if constexpr (std::is_integral_v<T>) {
        return true;
    } else {
        // 浮動小数点の場合: 各係数が整数に十分近いかチェック
        for (int i = 0; i <= p.degree(); ++i) {
            double v = static_cast<double>(p[i]);
            if (std::abs(v - std::round(v)) > 1e-9) return false;
        }
        return true;
    }
}

// 多項式を整数化して返す (有理数 Lagrange 結果 → 整数係数チェック)
template<typename T>
bool tryRoundToInt(const Polynomial<double>& src, Polynomial<T>& dst) {
    std::vector<T> coeffs(src.degree() + 1);
    for (int i = 0; i <= src.degree(); ++i) {
        double v = src[i];
        double r = std::round(v);
        if (std::abs(v - r) > 1e-6) return false;
        coeffs[i] = static_cast<T>(static_cast<int64_t>(r));
    }
    dst = Polynomial<T>(std::move(coeffs));
    return true;
}

} // namespace detail

// ================================================================
// 整数多項式の GCD (擬剰余ベース)
// ================================================================

/// 整数多項式の擬剰余 (pseudo-remainder)
/// prem(f, g) = lc(g)^(deg(f)-deg(g)+1) * f mod g
/// 結果は必ず整数係数を保つ
template<PolynomialCoeff T>
[[nodiscard]] Polynomial<T> pseudoRemainder(
    const Polynomial<T>& f, const Polynomial<T>& g)
{
    static_assert(std::is_integral_v<T>);
    if (g.isZero()) return f;
    if (f.degree() < g.degree()) return f;

    int delta = f.degree() - g.degree() + 1;
    T lc_g = g.leadingCoefficient();

    // lc(g)^delta * f
    T scale = T(1);
    for (int i = 0; i < delta; ++i) scale *= lc_g;

    Polynomial<T> r = f * scale;
    auto dr = r.divmod(g);
    return dr.remainder;
}

/// 整数多項式の GCD (subresultant PRS)
/// 各ステップで擬剰余を取り、primitivePart で係数を縮小する
template<PolynomialCoeff T>
[[nodiscard]] Polynomial<T> intPolyGcd(Polynomial<T> a, Polynomial<T> b) {
    static_assert(std::is_integral_v<T>);

    // 零チェック
    if (a.isZero()) return b;
    if (b.isZero()) return a;

    // 両方 primitive にする
    a = primitivePart(a);
    b = primitivePart(b);

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

    while (!b.isZero()) {
        Polynomial<T> r = pseudoRemainder(a, b);
        if (r.isZero()) break;
        r = primitivePart(r);
        a = std::move(b);
        b = std::move(r);
    }

    // b が 0 なら a が GCD、そうでなければ b が GCD
    Polynomial<T>& result = b.isZero() ? a : b;

    // 主係数を正にする
    if (result.leadingCoefficient() < T(0))
        result = -result;

    return primitivePart(result);
}

// ================================================================
// 無平方分解 (Yun's algorithm)
// ================================================================

/// Yun のアルゴリズムで多項式 f を無平方因子に分解する
/// 結果: f = c * Π f_i^i  (f_i は無平方で互いに素)
/// 返値は FactoredInteger<Polynomial<T>> で、各因子と冪を保持
///
/// T は整数型 (int, int64_t 等)
template<PolynomialCoeff T>
[[nodiscard]] FactoredInteger<Polynomial<T>> squareFreeFactorization(
    const Polynomial<T>& f)
{
    static_assert(std::is_integral_v<T>,
                  "squareFreeFactorization requires integral coefficient type");

    FactoredInteger<Polynomial<T>> result;

    if (f.isZero() || f.degree() <= 0) {
        if (!f.isZero()) result.addFactor(f, 1);
        return result;
    }

    Polynomial<T> fd = f.derivative();
    if (fd.isZero()) {
        // f' = 0 → f は定数の可能性 (or char p で全冪が p の倍数)
        result.addFactor(f, 1);
        return result;
    }

    // 擬剰余ベースの整数多項式 GCD を使用
    Polynomial<T> pp = primitivePart(f);
    Polynomial<T> ppd = pp.derivative();

    Polynomial<T> g = intPolyGcd(pp, ppd);
    Polynomial<T> v = pp / g;
    Polynomial<T> w = ppd / g;

    for (int n = 1; !v.isZero() && v.degree() > 0; ++n) {
        Polynomial<T> vd = v.derivative();
        Polynomial<T> diff = w - vd;

        Polynomial<T> h = intPolyGcd(v, diff);

        if (h.degree() > 0) {
            // 主係数を正にする
            if (h.leadingCoefficient() < T(0))
                h = -h;
            result.addFactor(h, n);
        }

        v = v / h;
        if (diff.isZero()) break;
        w = diff / h;
    }

    // 残りの v が定数でなければ追加
    if (!v.isZero() && v.degree() > 0) {
        if (v.leadingCoefficient() < T(0)) v = -v;
        result.addFactor(v, 1);
    }

    return result;
}

// ================================================================
// 有理根定理による 1 次因子分離
// ================================================================

namespace detail {

/// 有理根定理: f(b/a) = 0 となる有理数 b/a を探し、
/// (ax - b) の因子を全て分離する
/// f は primitive (content = 1) であること
template<PolynomialCoeff T>
void extractLinearFactors(
    Polynomial<T>& f,
    FactoredInteger<Polynomial<T>>& factors,
    int multiplicity = 1)
{
    static_assert(std::is_integral_v<T>);

    while (f.degree() >= 1) {
        T f0 = f[0];
        T fn = f.leadingCoefficient();

        if (f0 < T(0)) f0 = -f0;
        if (fn < T(0)) fn = -fn;

        auto divs_f0 = divisors(static_cast<int64_t>(f0));
        auto divs_fn = divisors(static_cast<int64_t>(fn));

        // f0=0 なら x=0 が根 → (x) で割り切れる
        if (f[0] == T(0)) {
            divs_f0 = {0};
        }

        bool found = false;
        // 候補: x = ±b/a (b | f0, a | fn)
        // 整数多項式のまま判定: f(b/a) = 0 ⟺ f を (ax - b) で割った余りが 0
        for (auto a : divs_fn) {
            if (a == 0) continue;
            for (auto b : divs_f0) {
                for (int sign : {1, -1}) {
                    T bv = static_cast<T>(b * sign);
                    T av = static_cast<T>(a);
                    // 候補因子: (av*x - bv) = av * x + (-bv)
                    Polynomial<T> candidate({T(0) - bv, av});

                    auto dr = f.divmod(candidate);
                    if (dr.remainder.isZero()) {
                        if (candidate.leadingCoefficient() < T(0))
                            candidate = -candidate;
                        factors.addFactor(candidate, multiplicity);
                        f = std::move(dr.quotient);
                        found = true;
                        goto next_round;
                    }
                }
            }
        }
    next_round:
        if (!found) break;
    }
}

} // namespace detail

// ================================================================
// Kronecker 法による因数分解
// ================================================================

namespace detail {

/// Kronecker のアルゴリズム:
/// deg(f) = n の因子は最大 n/2 次 → n/2+1 個の点で評価し、
/// 値の約数の全組み合わせを Lagrange 補間して因子を試す
template<PolynomialCoeff T>
void kroneckerFactor(
    Polynomial<T>& f,
    FactoredInteger<Polynomial<T>>& factors,
    int multiplicity = 1)
{
    static_assert(std::is_integral_v<T>);

    while (f.degree() >= 2) {
        int deg = f.degree();
        int half = deg / 2;
        int npts = half + 1;

        // 評価点を選択 (0 中心)
        std::vector<int64_t> xs(npts);
        int offset = half / 2;
        for (int i = 0; i < npts; ++i)
            xs[i] = i - offset;

        // 各点での値を計算し、約数を列挙
        std::vector<std::vector<int64_t>> allDivs(npts);
        for (int i = 0; i < npts; ++i) {
            // f(xs[i]) を計算
            T val = f(static_cast<T>(xs[i]));
            int64_t v = static_cast<int64_t>(val);
            auto d = divisors(v);
            // 正と負の約数
            std::vector<int64_t> fullDivs;
            fullDivs.reserve(d.size() * 2);
            for (auto dv : d) {
                fullDivs.push_back(dv);
                fullDivs.push_back(-dv);
            }
            // v=0 の場合は特別: 任意の整数が約数 → 小さい範囲で列挙
            if (v == 0) {
                fullDivs.clear();
                for (int64_t k = -10; k <= 10; ++k) fullDivs.push_back(k);
            }
            allDivs[i] = std::move(fullDivs);
        }

        // 約数の全組み合わせを試す
        // 組み合わせ数 = Π |allDivs[i]|
        size_t totalCombs = 1;
        for (int i = 0; i < npts; ++i) {
            if (allDivs[i].empty()) { totalCombs = 0; break; }
            totalCombs *= allDivs[i].size();
            if (totalCombs > 100000) { totalCombs = 100000; break; }  // 上限
        }

        bool found = false;
        for (size_t idx = 0; idx < totalCombs; ++idx) {
            // idx をインデックスの組に分解
            std::vector<double> xsd(npts), ysd(npts);
            size_t rem = idx;
            for (int i = 0; i < npts; ++i) {
                size_t j = rem % allDivs[i].size();
                rem /= allDivs[i].size();
                xsd[i] = static_cast<double>(xs[i]);
                ysd[i] = static_cast<double>(allDivs[i][j]);
            }

            // Lagrange 補間で候補多項式を得る
            Polynomial<double> cand_d = lagrangeInterpolationPoly(xsd, ysd);

            // 整数係数かチェック
            Polynomial<T> candidate;
            if (!tryRoundToInt(cand_d, candidate)) continue;
            if (candidate.degree() <= 0) continue;
            if (candidate.degree() > half) continue;

            // 主係数が f の主係数を割り切るか
            T candLC = candidate.leadingCoefficient();
            if (candLC < T(0)) { candidate = -candidate; candLC = -candLC; }
            T fLC = f.leadingCoefficient();
            if (fLC < T(0)) fLC = -fLC;
            if (candLC == T(0) || fLC % candLC != T(0)) continue;

            // f を candidate で割ってみる
            auto dr = f.divmod(candidate);
            if (dr.remainder.isZero()) {
                factors.addFactor(candidate, multiplicity);
                f = std::move(dr.quotient);
                found = true;
                break;
            }
        }

        if (!found) break;  // これ以上分解できない
    }
}

} // namespace detail

// ================================================================
// 因数分解 (メインエントリ)
// ================================================================

/// 整数多項式 f を既約因子に分解する
/// 結果: FactoredInteger<Polynomial<T>> で f = content * Π factor_i ^ exp_i
///
/// アルゴリズム:
///   1. content (全係数の GCD) を定数因子として抽出
///   2. Yun の無平方分解で重複因子を分離
///   3. 各無平方因子に対して:
///      a. 有理根定理で 1 次因子を分離
///      b. Kronecker 法で高次因子を分離
template<PolynomialCoeff T>
[[nodiscard]] FactoredInteger<Polynomial<T>> factorize(const Polynomial<T>& f) {
    static_assert(std::is_integral_v<T>,
                  "factorize requires integral coefficient type");

    FactoredInteger<Polynomial<T>> result;

    if (f.isZero()) return result;

    // 1. content 抽出
    T c = content(f);
    Polynomial<T> prim = f;
    if (c != T(1) && c != T(0)) {
        prim = primitivePart(f);
        if (c != T(1)) {
            result.addFactor(Polynomial<T>(c), 1);
        }
    }

    // 主係数を正にする
    if (!prim.isZero() && prim.leadingCoefficient() < T(0)) {
        prim = -prim;
        // 符号を content に吸収 (result に -1 因子を追加)
        if (result.size() > 0 && result[0].degree() == 0) {
            // 既に定数因子があるなら符号反転
            result[0] = Polynomial<T>(T(0) - result[0][0]);
        } else {
            result.addFactor(Polynomial<T>(T(-1)), 1);
        }
    }

    if (prim.degree() <= 0) {
        if (!prim.isZero() && !(prim[0] == T(1))) {
            result.addFactor(prim, 1);
        }
        return result;
    }

    if (prim.degree() == 1) {
        result.addFactor(prim, 1);
        return result;
    }

    // 2. 無平方分解
    auto sqfree = squareFreeFactorization(prim);

    if (sqfree.empty()) {
        // 無平方分解が空 → prim 自体を因子として追加
        result.addFactor(prim, 1);
        return result;
    }

    // 3. 各無平方因子をさらに分解
    for (size_t i = 0; i < sqfree.size(); ++i) {
        Polynomial<T> fi = sqfree[i];
        int mult = sqfree.exponent(i);

        if (fi.degree() <= 0) continue;

        if (fi.degree() == 1) {
            if (fi.leadingCoefficient() < T(0)) fi = -fi;
            result.addFactor(fi, mult);
            continue;
        }

        // 3a. 有理根定理で 1 次因子を分離
        detail::extractLinearFactors(fi, result, mult);

        if (fi.degree() <= 0) continue;

        if (fi.degree() == 1) {
            if (fi.leadingCoefficient() < T(0)) fi = -fi;
            result.addFactor(fi, mult);
            continue;
        }

        // 3b. Kronecker 法で高次因子を分離
        detail::kroneckerFactor(fi, result, mult);

        // 残りがあれば既約因子として追加
        if (fi.degree() > 0) {
            if (fi.leadingCoefficient() < T(0)) fi = -fi;
            result.addFactor(fi, mult);
        }
    }

    return result;
}

// ========================================================================
// CAS-4: Resultant / Discriminant (Sylvester 行列ベース)
// ========================================================================

/// Sylvester 行列を構築
/// f(x) の次数 m, g(x) の次数 n → (m+n)×(m+n) 行列
template<PolynomialCoeff T>
Matrix<T> sylvesterMatrix(const Polynomial<T>& f, const Polynomial<T>& g) {
    int m = f.degree(), n = g.degree();
    int sz = m + n;
    Matrix<T> S(sz, sz, T(0));

    // f の係数を n 行分シフトして配置
    for (int i = 0; i < n; ++i)
        for (int j = 0; j <= m; ++j)
            S(i, i + j) = f[m - j];  // 降順

    // g の係数を m 行分シフトして配置
    for (int i = 0; i < m; ++i)
        for (int j = 0; j <= n; ++j)
            S(n + i, i + j) = g[n - j];  // 降順

    return S;
}

/// 終結式 Res(f, g) = det(Sylvester(f, g))
template<PolynomialCoeff T>
[[nodiscard]] T resultant(const Polynomial<T>& f, const Polynomial<T>& g) {
    if (f.degree() == 0 || g.degree() == 0) return T(0);
    auto S = sylvesterMatrix(f, g);
    // LU 分解で行列式を計算 (簡易: ガウス消去)
    int n = S.rows();
    T det = T(1);
    for (int i = 0; i < n; ++i) {
        // ピボット選択
        int pivot = -1;
        for (int k = i; k < n; ++k)
            if (S(k, i) != T(0)) { pivot = k; break; }
        if (pivot < 0) return T(0);
        if (pivot != i) {
            for (int j = 0; j < n; ++j) std::swap(S(i, j), S(pivot, j));
            det = -det;
        }
        det *= S(i, i);
        T inv = S(i, i);
        for (int k = i + 1; k < n; ++k) {
            T factor = S(k, i);
            if (factor == T(0)) continue;
            for (int j = i; j < n; ++j)
                S(k, j) = S(k, j) * inv - S(i, j) * factor;
            // 整数除算でスケール調整
            if (i > 0) {
                for (int j = i; j < n; ++j)
                    S(k, j) /= det / S(i, i);  // 仮 (浮動小数点の場合)
            }
        }
    }
    return det;
}

/// 判別式 Disc(f) = (-1)^(n(n-1)/2) * Res(f, f') / a_n
/// a_n = 最高次係数
template<PolynomialCoeff T>
[[nodiscard]] T discriminant(const Polynomial<T>& f) {
    int n = f.degree();
    if (n < 1) return T(0);
    auto fp = f.derivative();
    T res = resultant(f, fp);
    T an = f[n];
    T sign = ((n * (n - 1) / 2) % 2 == 0) ? T(1) : T(-1);
    if (an != T(0)) return sign * res / an;
    return sign * res;
}

// ========================================================================
// CAS-10: 整数分割
// ========================================================================

/// n の全分割を列挙 (降順パーティション)
[[nodiscard]] inline std::vector<std::vector<int>> integerPartitions(int n) {
    std::vector<std::vector<int>> result;
    if (n <= 0) { result.push_back({}); return result; }

    std::function<void(int, int, std::vector<int>&)> gen;
    gen = [&](int remaining, int maxPart, std::vector<int>& current) {
        if (remaining == 0) {
            result.push_back(current);
            return;
        }
        for (int k = std::min(remaining, maxPart); k >= 1; --k) {
            current.push_back(k);
            gen(remaining - k, k, current);
            current.pop_back();
        }
    };
    std::vector<int> buf;
    gen(n, n, buf);
    return result;
}

/// 分割数 p(n) (動的計画法)
[[nodiscard]] inline uint64_t partitionCount(int n) {
    if (n < 0) return 0;
    std::vector<uint64_t> dp(n + 1, 0);
    dp[0] = 1;
    for (int k = 1; k <= n; ++k)
        for (int j = k; j <= n; ++j)
            dp[j] += dp[j - k];
    return dp[n];
}

// ========================================================================
// CAS-3b-5: 絶対因数分解 (Trager)
// ========================================================================
//
// Z[x] (or Q[x]) の多項式を代数的数体 Q(α) 上で因数分解する。
// 例: x^4+1 は Z[x] で既約だが Q(√2)[x] で (x²+√2x+1)(x²-√2x+1)
//
// Trager のアルゴリズム:
//   1. f(x) の splitting field の最小多項式を求める
//   2. Norm_{Q(α)/Q}(f(x) - t·g(x)) を t について計算 (resultant)
//   3. 結果の Z[t] 多項式を因数分解 → α の最小多項式の候補
//   4. 各候補で GCD を取って因子を復元
//
// 簡易版: f(x) が可分 (separable) で deg ≤ 小 の場合に対応。
// 出力: (既約因子, 最小多項式) のペア。最小多項式が x なら Q 上で既約。

/// 絶対因数分解の結果
template<typename T>
struct AbsoluteFactorResult {
    Polynomial<T> factor;       ///< Q(α)[x] 上の既約因子 (係数は α の多項式で表現)
    Polynomial<T> minPoly;      ///< α の最小多項式 (Q[t] 上)
    int degree;                 ///< 因子の次数
};

/// 絶対因数分解: f(x) を代数閉包上の既約因子に分解
/// @param f  Z[x] の多項式 (square-free を仮定)
/// @return   分解情報のリスト
///           minPoly が {0, 1} (= t) なら Q 上で既約
///           minPoly が非自明なら Q(α)[x] 上の因子
template<typename T>
std::vector<AbsoluteFactorResult<T>> absoluteFactorize(const Polynomial<T>& f) {
    static_assert(std::is_integral_v<T>, "absoluteFactorize requires integral type");
    std::vector<AbsoluteFactorResult<T>> results;

    int n = f.degree();
    if (n <= 0) return results;

    // 1次は Q 上で既約
    if (n == 1) {
        results.push_back({f, Polynomial<T>({T(0), T(1)}), 1});
        return results;
    }

    // Z[x] 上で因数分解を試みる
    auto zFactors = factorize(f);
    if (zFactors.size() > 1 || (zFactors.size() == 1 && zFactors.exponent(0) > 1)) {
        // Z[x] 上で可約 → 各因子を再帰的に絶対因数分解
        for (size_t i = 0; i < zFactors.size(); ++i) {
            if (zFactors[i].degree() <= 0) continue;
            auto sub = absoluteFactorize(zFactors[i]);
            for (auto& r : sub) {
                r.degree = zFactors[i].degree();
                for (int j = 0; j < zFactors.exponent(i); ++j)
                    results.push_back(r);
            }
        }
        return results;
    }

    // Z[x] 上で既約。代数拡大上で分解できるか調べる。
    // Trager の方法: h(t) = Res_x(f(x), x - t·f'(x)/n) を計算
    // h(t) の因数分解から分裂体の構造がわかる。

    // 簡易判定: 次数 2 の場合、判別式が完全平方数でなければ Q(√D) に分解
    if (n == 2) {
        T a = f[2], b = f[1], c = f[0];
        T D = b * b - T(4) * a * c;  // 判別式
        // D が完全平方数か?
        if (D >= T(0)) {
            T sqrtD = static_cast<T>(std::sqrt(static_cast<double>(D)));
            if (sqrtD * sqrtD == D) {
                // Q 上で分解可能 (既に factorize で処理済みのはず)
                results.push_back({f, Polynomial<T>({T(0), T(1)}), n});
                return results;
            }
        }
        // Q(√D) 上で 1 次因子に分解
        // f(x) = a(x - (-b+√D)/(2a))(x - (-b-√D)/(2a))
        // 最小多項式: t² - D
        results.push_back({f, Polynomial<T>({-D, T(0), T(1)}), n});
        return results;
    }

    // 次数 4: x^4 + bx^2 + c 型 (biquadratic)
    if (n == 4 && f[3] == T(0) && f[1] == T(0)) {
        // y = x² として y² + by + c を因数分解
        Polynomial<T> g({f[0], f[2], f[4]});  // f[4]y² + f[2]y + f[0]
        T D = f[2] * f[2] - T(4) * f[4] * f[0];
        if (D >= T(0)) {
            T sqrtD = static_cast<T>(std::sqrt(static_cast<double>(D)));
            if (sqrtD * sqrtD == D) {
                // Q 上で 2 つの 2 次因子に分解
                // (これも factorize で見つかるはず → Q上既約として返す)
                results.push_back({f, Polynomial<T>({T(0), T(1)}), n});
                return results;
            }
        }
        // Q(√D) 上で 2 次因子に分解
        results.push_back({f, Polynomial<T>({-D, T(0), T(1)}), n});
        return results;
    }

    // 一般の場合: Q 上で既約として返す (完全な Trager は将来課題)
    results.push_back({f, Polynomial<T>({T(0), T(1)}), n});
    return results;
}

} // namespace calx
