// Copyright (C) 2026 Kiyotsugu Arai // SPDX-License-Identifier: LGPL-3.0-or-later // Polynomial_factorization.hpp // Exact factorization of integer polynomials // // Algorithm: // 1. Extract content (GCD of all coefficients) // 2. Separate linear factors via the rational-root theorem // 3. Yun's square-free decomposition (separate repeated factors via gcd(f, f')) // 4. Try higher-degree factors via Kronecker's method // // Usage: // #include

//   #include //   #include // result type

#pragma once

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

// ================================================================
// Internal helpers
// ================================================================

namespace detail {

// Enumerate divisors of an integer (positive divisors only, sorted)
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;
}

// Enumerate divisors of an integer (int64_t version)
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;
}

// Recover polynomial coefficients via Lagrange interpolation
// xs, ys: evaluation points and values (n+1 of them -> polynomial of degree at most 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;
}

// Check whether the polynomial has integer coefficients (always true when T is an integer type)
template bool hasIntegerCoeffs(const Polynomial & p) {
    if constexpr (std::is_integral_v) {
        return true;
    } else {
        // Floating-point case: check whether each coefficient is sufficiently close to an integer
        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;
    }
}

// Convert the polynomial to integer coefficients and return (Lagrange rational result -> integer-coefficient check)
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 of integer polynomials (pseudo-remainder based)
// ================================================================

/// Pseudo-remainder of integer polynomials
/// prem(f, g) = lc(g)^(deg(f)-deg(g)+1) * f mod g
/// The result always keeps integer coefficients
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 of integer polynomials (subresultant PRS)
/// At each step take the pseudo-remainder and shrink the coefficients via primitivePart
template [[nodiscard]] Polynomial intPolyGcd(Polynomial a, Polynomial b) {
    static_assert(std::is_integral_v);

    // Zero check
    if (a.isZero()) return b;
    if (b.isZero()) return a;

    // Make both primitive
    a = primitivePart(a);
    b = primitivePart(b);

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

    // If b is 0, a is the GCD; otherwise b is the GCD
    Polynomial & result = b.isZero() ? a : b;

    // Make the leading coefficient positive
    if (result.leadingCoefficient() < T(0))
        result = -result;

    return primitivePart(result);
}

// ================================================================
// Square-free decomposition (Yun's algorithm)
// ================================================================

/// Decompose polynomial f into square-free factors via Yun's algorithm
/// Result: f = c * Product f_i^i  (f_i are square-free and mutually coprime)
/// Returns FactoredInteger > holding each factor and its exponent
///
/// T is an integer type (int, int64_t, etc.)
template [[nodiscard]] FactoredInteger