alkahest-cas 3.12.0

//! Dense `fmpz_mat_t` wrapper for Hermite / Smith normal forms.

use super::ffi;
use super::integer::FlintInteger;

/// Owned `m × n` integer matrix backed by FLINT.
pub(crate) struct FlintMat {
    inner: ffi::FmpzMatStruct,
}

impl FlintMat {
    pub(crate) fn new(rows: usize, cols: usize) -> Self {
        assert!(rows <= i64::MAX as usize && cols <= i64::MAX as usize);
        #[cfg(not(flint3_stride))]
        let mut inner = ffi::FmpzMatStruct {
            entries: std::ptr::null_mut(),
            r: 0,
            c: 0,
            rows: std::ptr::null_mut(),
        };
        #[cfg(flint3_stride)]
        let mut inner = ffi::FmpzMatStruct {
            entries: std::ptr::null_mut(),
            r: 0,
            c: 0,
            stride: 0,
        };
        crate::flint::note_thread_uses_flint();
        unsafe { ffi::fmpz_mat_init(&mut inner, rows as ffi::slong, cols as ffi::slong) };
        Self { inner }
    }

    pub(crate) fn rows(&self) -> usize {
        self.inner.r as usize
    }

    pub(crate) fn cols(&self) -> usize {
        self.inner.c as usize
    }

    unsafe fn entry_ptr(&mut self, i: usize, j: usize) -> *mut ffi::fmpz {
        debug_assert!(i < self.rows() && j < self.cols());
        #[cfg(not(flint3_stride))]
        {
            let row_start = *self.inner.rows.add(i);
            row_start.add(j)
        }
        #[cfg(flint3_stride)]
        {
            let stride = self.inner.stride as usize;
            self.inner.entries.add(i * stride + j)
        }
    }

    pub(crate) unsafe fn entry_const(&self, i: usize, j: usize) -> *const ffi::fmpz {
        debug_assert!(i < self.rows() && j < self.cols());
        #[cfg(not(flint3_stride))]
        {
            let row_start = *self.inner.rows.add(i);
            row_start.add(j)
        }
        #[cfg(flint3_stride)]
        {
            let stride = self.inner.stride as usize;
            self.inner.entries.add(i * stride + j)
        }
    }

    pub(crate) fn set_entry(&mut self, i: usize, j: usize, v: &FlintInteger) {
        unsafe {
            let e = self.entry_ptr(i, j);
            ffi::fmpz_set(e, v.inner_ptr());
        }
    }

    pub(crate) fn get_flint(&self, i: usize, j: usize) -> FlintInteger {
        let mut out = FlintInteger::new();
        unsafe {
            ffi::fmpz_set(out.inner_mut_ptr(), self.entry_const(i, j));
        }
        out
    }

    /// Exact determinant (`fmpz_mat_det`). The matrix must be square.
    pub(crate) fn det(&self) -> FlintInteger {
        assert_eq!(self.rows(), self.cols(), "det of a non-square matrix");
        let mut d = FlintInteger::new();
        unsafe { ffi::fmpz_mat_det(d.inner_mut_ptr(), &self.inner) };
        d
    }

    /// Fraction-free reduced row echelon form (`fmpz_mat_rref`): returns
    /// `(rank, den, B)` with `B / den` the (unique) reduced row echelon form of
    /// `self` over ℚ, `den > 0` when the rank is positive.
    pub(crate) fn rref(&self) -> (usize, FlintInteger, FlintMat) {
        let mut b = FlintMat::new(self.rows(), self.cols());
        let mut den = FlintInteger::new();
        let rank = unsafe { ffi::fmpz_mat_rref(&mut b.inner, den.inner_mut_ptr(), &self.inner) };
        (rank as usize, den, b)
    }

    pub(crate) fn hnf_transform(&self, h: &mut FlintMat, u: &mut FlintMat) {
        unsafe {
            ffi::fmpz_mat_hnf_transform(&mut h.inner, &mut u.inner, &self.inner);
        }
    }

    /// Hermite normal form **without** a transform matrix.
    ///
    /// [`FlintMat::hnf_transform`] also produces the unimodular `U` with
    /// `U·A = H`, and on a rank-deficient input `U`'s entries can be orders of
    /// magnitude larger than `H`'s — the cost of computing it dominates.
    /// Callers that only want the echelon basis should use this.
    pub(crate) fn hnf(&self, h: &mut FlintMat) {
        unsafe {
            ffi::fmpz_mat_hnf(&mut h.inner, &self.inner);
        }
    }

    /// Smith normal form, diagonal only — no transform matrices.
    ///
    /// `matrix::normal_form` uses it two ways: as the oracle its own
    /// transform-tracking port is tested against, and as the implementation of
    /// `smith_invariants` for callers that never look at `U` and `V`.
    pub(crate) fn snf_diagonal(&self, s: &mut FlintMat) {
        unsafe {
            ffi::fmpz_mat_snf(&mut s.inner, &self.inner);
        }
    }
}

/// Shape checks used only from `matrix::normal_form` unit tests.
#[cfg(test)]
impl FlintMat {
    pub(crate) fn is_in_hnf(&self) -> bool {
        unsafe { ffi::fmpz_mat_is_in_hnf(&self.inner) != 0 }
    }

    pub(crate) fn is_in_snf(&self) -> bool {
        unsafe { ffi::fmpz_mat_is_in_snf(&self.inner) != 0 }
    }

    pub(crate) fn equals(&self, other: &FlintMat) -> bool {
        unsafe { ffi::fmpz_mat_equal(&self.inner, &other.inner) != 0 }
    }
}

impl Drop for FlintMat {
    fn drop(&mut self) {
        unsafe { ffi::fmpz_mat_clear(&mut self.inner) };
    }
}