// Copyright (C) 2026 Kiyotsugu Arai // SPDX-License-Identifier: LGPL-3.0-or-later // Polynomial_factorization.hpp // 整数多項式の厳密な因数分解 // // アルゴリズム: // 1. content 抽出 (全係数の GCD) // 2. 有理根定理による 1 次因子分離 // 3. Yun の無平方分解 (gcd(f, f') による重複因子分離) // 4. Kronecker 法による高次因子の試行 // // 使用法: // #include

//   #include //   #include // 結果の型

#pragma once

#include #include #include #include #include #include #include #include #include namespace sangi {

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

namespace detail {

// 整数の約数を列挙 (正の約数のみ、ソート済み)
inline std::vector divisors(int n) {
    if (n < 0) n = -n;
    if (n == 0) return {};
    std::vector 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 divisors(int64_t n) {
    if (n < 0) n = -n;
    if (n == 0) return {};
    std::vector 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 Polynomial lagrangeInterpolationPoly(
    const std::vector & xs,
    const std::vector & ys)
{
    int n = static_cast (xs.size());
    Polynomial result(T(0));
    Polynomial 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 basis = one;
        T denom = T(1);
        for (int j = 0; j < n; ++j) {
            if (j == i) continue;
            // (x - x_j)
            basis = basis * Polynomial ({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 bool hasIntegerCoeffs(const Polynomial & p) {
    if constexpr (std::is_integral_v) {
        return true;
    } else {
        // 浮動小数点の場合: 各係数が整数に十分近いかチェック
        for (int i = 0; i <= p.degree(); ++i) {
            double v = static_cast (p[i]);
            if (std::abs(v - std::round(v)) > 1e-9) return false;
        }
        return true;
    }
}

// 多項式を整数化して返す (有理数 Lagrange 結果 → 整数係数チェック)
template bool tryRoundToInt(const Polynomial & src, Polynomial & dst) {
    std::vector 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 (static_cast (r));
    }
    dst = Polynomial (std::move(coeffs));
    return true;
}

} // namespace detail

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

/// 整数多項式の擬剰余 (pseudo-remainder)
/// prem(f, g) = lc(g)^(deg(f)-deg(g)+1) * f mod g
/// 結果は必ず整数係数を保つ
template [[nodiscard]] Polynomial pseudoRemainder(
    const Polynomial & f, const Polynomial & g)
{
    static_assert(std::is_integral_v);
    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 r = f * scale;
    auto dr = r.divmod(g);
    return dr.remainder;
}

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

    // 零チェック
    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 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 & 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 > で、各因子と冪を保持
///
/// T は整数型 (int, int64_t 等)
template [[nodiscard]] FactoredInteger