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

// HNF.hpp
//
// Hermite Normal Form
//
// Transforms an integer matrix A (m×n) into an upper triangular matrix H.
// H = U·A (U is an m×m unimodular matrix, det(U) = ±1).
// Foundation of lattice theory: used in LLL reduction, lattice equivalence
// testing, and integer linear equations.
//
// Definition (row HNF):
//   H is upper triangular (or the leading nonzero entry of each row is leftmost in its column)
//   The pivot of each row is positive
//   Entries above a pivot in the same column satisfy 0 ≤ h_{i,j} < h_{j,j}
//
// Algorithm: classical column-operation based HNF (column elimination using GCD)
//
// Main API:
//   hermiteNormalForm(A)            — returns the row HNF
//   hermiteNormalFormWithTransform(A) — the pair {H, U} (H = U·A)

#ifndef SANGI_LINALG_HNF_HPP
#define SANGI_LINALG_HNF_HPP

#include <math/core/matrix.hpp>
#include <math/core/mp/Int.hpp>
#include <math/core/mp/Int/IntGCD.hpp>   // IntGCD::extendedGcd (used by extended_gcd below)
#include <cstddef>
#include <stdexcept>
#include <utility>
#include <cmath>

namespace sangi {

/// HNF result (with transformation matrix)
template<typename T>
struct HNFResult {
    Matrix<T> H;   ///< Hermite Normal Form
    Matrix<T> U;   ///< unimodular transformation matrix (H = U·A)
};

namespace detail_hnf {

// Integer GCD + Bezout coefficients (ax + by = gcd)
// Supports T = int, int64_t, Int, etc.
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 specialization: uses 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};
}

// Sign of 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();
}

// Absolute value of 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++ % is signed; make 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

/// Computes the row Hermite Normal Form
/// @param A  integer matrix (m×n)
/// @return   HNF matrix H (m×n)
template<typename T>
Matrix<T> hermiteNormalForm(const BaseMatrix<T>& A) {
    return hermiteNormalFormWithTransform(A).H;
}

/// Computes the row Hermite Normal Form + transformation matrix
/// @param A  integer matrix (m×n)
/// @return   {H, U} where H = U·A, U is unimodular
template<typename T>
HNFResult<T> hermiteNormalFormWithTransform(const BaseMatrix<T>& A) {
    using std::size_t;
    size_t m = A.rows(), n = A.cols();
    if (m == 0 || n == 0) return {Matrix<T>(m, n, T(0)), Matrix<T>(m, m, T(0))};

    // H = copy of A, U = identity matrix
    Matrix<T> H(A);
    Matrix<T> U(m, m, T(0));
    for (size_t i = 0; i < m; ++i) U(i, i) = T(1);

    // Eliminate column by column
    size_t pivotRow = 0;
    for (size_t col = 0; col < n && pivotRow < m; ++col) {
        // Choose the nonzero row with smallest absolute value in this column (limits coefficient growth)
        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;  // this column is all zeros

        // Swap pivotRow and 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));
        }

        // Eliminate rows below pivotRow using 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 transform: [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;
            }
        }

        // Make the pivot positive
        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);
        }

        // Reduce entries above the pivot to 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;
                // Adjust so that H(i, col) - q*pivot is not negative
                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)};
}

/// Verifies HNF correctness: whether H = U·A and H satisfies the HNF conditions
template<typename T>
bool verifyHNF(const BaseMatrix<T>& A, const BaseMatrix<T>& H, const BaseMatrix<T>& U) {
    // Verify 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 sangi

#endif // SANGI_LINALG_HNF_HPP