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

// HNF.hpp
//
// Hermite 標準形 (Hermite Normal Form)
//
// 整数行列 A (m×n) を上三角行列 H に変換する。
// H = U·A (U は m×m のユニモジュラー行列、det(U) = ±1)。
// 格子理論の基盤: LLL 簡約、格子の同値判定、整数線形方程式に使用。
//
// 定義 (行 HNF):
//   H は上三角 (または各行の先頭非零要素が列で最左)
//   各行のピボットは正
//   ピボットより上の同列要素は 0 ≤ h_{i,j} < h_{j,j}
//
// アルゴリズム: 列操作ベースの古典的 HNF (GCD を用いた列掃き出し)
//
// 主な API:
//   hermiteNormalForm(A)            — 行 HNF を返す
//   hermiteNormalFormWithTransform(A) — {H, U} のペア (H = U·A)

#ifndef CALX_LINALG_HNF_HPP
#define CALX_LINALG_HNF_HPP

#include <math/core/matrix.hpp>
#include <math/core/mp/Int.hpp>
#include <cstddef>
#include <stdexcept>
#include <utility>
#include <cmath>

namespace calx {

/// HNF の結果 (変換行列付き)
template<typename T>
struct HNFResult {
    Matrix<T> H;   ///< Hermite 標準形
    Matrix<T> U;   ///< ユニモジュラー変換行列 (H = U·A)
};

namespace detail_hnf {

// 整数の GCD + Bezout 係数 (ax + by = gcd)
// T = int, int64_t, Int 等に対応
template<typename T>
struct ExtGcdResult { T g, x, y; };

template<typename T>
ExtGcdResult<T> extgcd(T a, T b) {
    if (b == T(0)) return {a, T(1), T(0)};
    auto [g, x1, y1] = extgcd(b, a % b);
    return {g, y1, x1 - (a / b) * y1};
}

// Int 特殊化: IntGCD::extendedGcd を使用
template<>
inline ExtGcdResult<Int> extgcd(Int a, Int b) {
    if (b.isZero()) return {a, Int(1), Int(0)};
    Int x, y;
    Int g = IntGCD::extendedGcd(a, b, x, y);
    return {g, x, y};
}

// T の符号 (-1, 0, 1)
template<typename T>
int signOf(const T& v) {
    if (v > T(0)) return 1;
    if (v < T(0)) return -1;
    return 0;
}

template<>
inline int signOf(const Int& v) {
    return v.getSign();
}

// T の絶対値
template<typename T>
T absVal(const T& v) { return (v < T(0)) ? -v : v; }

template<>
inline Int absVal(const Int& v) { return abs(v); }

// mod (C++ の % は符号付き、0 ≤ result < |b| にする)
template<typename T>
T posMod(const T& a, const T& b) {
    T r = a % b;
    if (r < T(0)) r += absVal(b);
    return r;
}

} // namespace detail_hnf

/// 行 Hermite 標準形を計算
/// @param A  整数行列 (m×n)
/// @return   HNF 行列 H (m×n)
template<typename T>
Matrix<T> hermiteNormalForm(const Matrix<T>& A) {
    return hermiteNormalFormWithTransform(A).H;
}

/// 行 Hermite 標準形 + 変換行列を計算
/// @param A  整数行列 (m×n)
/// @return   {H, U} where H = U·A, U はユニモジュラー
template<typename T>
HNFResult<T> hermiteNormalFormWithTransform(const Matrix<T>& A) {
    using std::size_t;
    size_t m = A.rows(), n = A.cols();
    if (m == 0 || n == 0) return {A, Matrix<T>(m, m, T(0))};

    // H = A のコピー、U = 単位行列
    Matrix<T> H = A;
    Matrix<T> U(m, m, T(0));
    for (size_t i = 0; i < m; ++i) U(i, i) = T(1);

    // 列ごとに掃き出し
    size_t pivotRow = 0;
    for (size_t col = 0; col < n && pivotRow < m; ++col) {
        // この列で絶対値最小の非ゼロ行を選択 (係数膨張を抑制)
        size_t nonzeroRow = m;
        T bestAbs{};
        for (size_t i = pivotRow; i < m; ++i) {
            if (H(i, col) == T(0)) continue;
            T av = detail_hnf::absVal(H(i, col));
            if (nonzeroRow == m || av < bestAbs) {
                bestAbs = av;
                nonzeroRow = i;
            }
        }
        if (nonzeroRow == m) continue;  // この列は全てゼロ

        // pivotRow と nonzeroRow を交換
        if (nonzeroRow != pivotRow) {
            for (size_t j = 0; j < n; ++j) std::swap(H(pivotRow, j), H(nonzeroRow, j));
            for (size_t j = 0; j < m; ++j) std::swap(U(pivotRow, j), U(nonzeroRow, j));
        }

        // pivotRow より下の行を GCD で消去
        for (size_t i = pivotRow + 1; i < m; ++i) {
            if (H(i, col) == T(0)) continue;

            auto [g, x, y] = detail_hnf::extgcd(H(pivotRow, col), H(i, col));
            T a = H(pivotRow, col) / g;
            T b = H(i, col) / g;

            // 行変換: [row_pivot, row_i] = [[x, y], [-b, a]] * [row_pivot, row_i]
            for (size_t j = 0; j < n; ++j) {
                T hp = H(pivotRow, j), hi = H(i, j);
                H(pivotRow, j) = x * hp + y * hi;
                H(i, j) = -b * hp + a * hi;
            }
            for (size_t j = 0; j < m; ++j) {
                T up = U(pivotRow, j), ui = U(i, j);
                U(pivotRow, j) = x * up + y * ui;
                U(i, j) = -b * up + a * ui;
            }
        }

        // ピボットを正にする
        if (detail_hnf::signOf(H(pivotRow, col)) < 0) {
            for (size_t j = 0; j < n; ++j) H(pivotRow, j) = -H(pivotRow, j);
            for (size_t j = 0; j < m; ++j) U(pivotRow, j) = -U(pivotRow, j);
        }

        // ピボットより上の要素を 0 ≤ h < pivot に縮約
        T pivot = H(pivotRow, col);
        if (pivot != T(0)) {
            for (size_t i = 0; i < pivotRow; ++i) {
                if (H(i, col) == T(0)) continue;
                T q = H(i, col) / pivot;
                // H(i, col) - q*pivot が負にならないよう調整
                if (H(i, col) - q * pivot < T(0)) q -= T(1);
                if (q == T(0)) continue;
                for (size_t j = 0; j < n; ++j) H(i, j) -= q * H(pivotRow, j);
                for (size_t j = 0; j < m; ++j) U(i, j) -= q * U(pivotRow, j);
            }
        }

        ++pivotRow;
    }

    return {std::move(H), std::move(U)};
}

/// HNF が正しいか検証: H = U·A かつ H が HNF 条件を満たすか
template<typename T>
bool verifyHNF(const Matrix<T>& A, const Matrix<T>& H, const Matrix<T>& U) {
    // H = U·A を検証
    size_t m = A.rows(), n = A.cols();
    for (size_t i = 0; i < m; ++i) {
        for (size_t j = 0; j < n; ++j) {
            T sum = T(0);
            for (size_t k = 0; k < m; ++k) sum += U(i, k) * A(k, j);
            if (sum != H(i, j)) return false;
        }
    }
    return true;
}

} // namespace calx

#endif // CALX_LINALG_HNF_HPP