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

// Polynomial_numerical.hpp
// 多項式の数値アルゴリズム (線形代数依存)
//
// SVD ベースの近似 GCD など、Matrix / decomposition に依存する
// 数値的多項式アルゴリズムをここに配置する。

#pragma once

#include <math/core/Polynomial.hpp>
#include <math/core/matrix.hpp>
#include <math/core/vector.hpp>
#include <math/linalg/decomposition.hpp>

namespace calx {

/**
 * @brief 近似 GCD (Sylvester 行列 + SVD)
 *
 * Sylvester 部分行列 S_k の特異値を分析し、GCD の次数を推定する。
 * S_k の最小特異値が tolerance 以下なら rank-deficient → GCD 次数 >= k。
 * ε-ユークリッド法より高精度で、ノイズに対してロバスト。
 *
 * アルゴリズム:
 *   1. k = min(deg a, deg b) から 1 まで降順に探索
 *   2. 各 k に対し Sylvester 部分行列 S_k を構築
 *   3. SVD で最小特異値を取得し、tolerance 以下なら GCD 次数 = k
 *   4. 特定した次数まで ε-ユークリッド法で実際の GCD を計算
 *
 * @param a         入力多項式
 * @param b         入力多項式
 * @param tolerance 特異値がこれ以下なら零とみなす (デフォルト 1e-10)
 * @return 近似的な GCD (モニック化済み)
 */
template<PolynomialCoeff T>
[[nodiscard]] Polynomial<T> approximateGCD_SVD(
    const Polynomial<T>& a, const Polynomial<T>& b,
    double tolerance = 1e-10) {

    int m = a.degree();
    int n = b.degree();
    if (m < 0 || n < 0) return Polynomial<T>({T(1)});
    if (m < n) return approximateGCD_SVD(b, a, tolerance);

    // GCD の次数 k を min(m,n) から 0 まで探索
    int max_k = n;
    int gcd_degree = 0;

    for (int k = max_k; k >= 1; --k) {
        // S_k: (m+n-2k+2) 行 × (n-k+1 + m-k+1) 列
        int num_rows = m + n - 2 * k + 2;
        int cols_a = n - k + 1;
        int cols_b = m - k + 1;
        int num_cols = cols_a + cols_b;

        if (num_rows <= 0 || num_cols <= 0) continue;

        // Sylvester 部分行列を構築
        Matrix<double> S(num_rows, num_cols);
        for (int r = 0; r < num_rows; ++r)
            for (int c = 0; c < num_cols; ++c)
                S(r, c) = 0.0;

        // a の列: x^j * a(x) の係数 (j = 0..cols_a-1)
        for (int j = 0; j < cols_a; ++j) {
            for (int i = 0; i <= m; ++i) {
                int row = i + j;
                if (row < num_rows)
                    S(row, j) = static_cast<double>(a[i]);
            }
        }

        // b の列: x^j * b(x) の係数 (j = 0..cols_b-1)
        for (int j = 0; j < cols_b; ++j) {
            for (int i = 0; i <= n; ++i) {
                int row = i + j;
                if (row < num_rows)
                    S(row, cols_a + j) = static_cast<double>(b[i]);
            }
        }

        // SVD で最小特異値を取得
        auto [U, sigma, Vt] = svd_decomposition(S);

        double min_sv = sigma[sigma.size() - 1];
        if (min_sv <= tolerance) {
            gcd_degree = k;
            break;
        }
    }

    if (gcd_degree == 0) {
        return Polynomial<T>({T(1)});  // 互いに素
    }

    // GCD を実際に計算: ε-ユークリッド法で gcd_degree 次まで
    Polynomial<T> r0 = a;
    Polynomial<T> r1 = b;
    if (r0.degree() < r1.degree()) std::swap(r0, r1);

    while (r1.degree() > gcd_degree) {
        Polynomial<T> r = r0 % r1;
        r0 = std::move(r1);
        r1 = std::move(r);
    }

    Polynomial<T>& result = (r1.degree() == gcd_degree) ? r1 : r0;

    // モニック化
    if (!result.isZero()) {
        T lc = result.leadingCoefficient();
        if (!(lc == T(1)))
            result = result / lc;
    }
    return result;
}

} // namespace calx
