feanor-math 3.5.18

use std::alloc::Allocator;

use crate::algorithms::lll::float::lll_quadratic_form;
use crate::algorithms::matmul::*;
use crate::field::*;
use crate::homomorphism::*;
use crate::integer::*;
use crate::matrix::transform::{TransformCols, TransformTarget};
use crate::matrix::*;
use crate::ordered::{OrderedRing, OrderedRingStore};
use crate::ring::*;
use crate::rings::approx_real::float::Real64;
use crate::rings::fraction::FractionFieldStore;
use crate::rings::rational::*;
use crate::seq::*;

/// Size-reduces `target` w.r.t. the GSO matrix, and also sends the performed
/// operations to `col_ops`.
fn size_reduce<R, I, V, T>(
    ring: R,
    mut target: SubmatrixMut<V, El<R>>,
    target_j: usize,
    gso_part: Submatrix<V, El<R>>,
    col_ops: &mut T,
) where
    R: RingStore<Type = RationalFieldBase<I>> + Copy,
    I: RingStore,
    I::Type: IntegerRing,
    V: AsPointerToSlice<El<RationalField<I>>>,
    T: TransformTarget<R::Type>,
{
    assert!(!ring.get_ring().is_approximate());
    for j in (0..gso_part.col_count()).rev() {
        let target_const = target.as_const();
        let mu = target_const.at(j, 0);
        let factor = ring.base_ring().rounded_div(
            ring.base_ring().clone_el(ring.get_ring().num(mu)),
            ring.get_ring().den(mu),
        );
        let factor = ring.inclusion().map(factor);
        col_ops.subtract(ring, j, target_j, &factor);
        ring.sub_assign_ref(target.at_mut(j, 0), &factor);
        for k in 0..j {
            ring.sub_assign(target.at_mut(k, 0), ring.mul_ref(gso_part.at(k, j), &factor));
        }
    }
}

/// Updates the "gso"-matrix resulting from swapping cols i and i + 1.
///
/// This uses explicit formulas, in particular if `b_i` and `b_(i + 1)`
/// represent the vectors, `b_i*` and `b_(i + 1)*` represent their GS-orthogonalizations
/// and the corresponding `b'` values are the ones after swapping, then we find
/// then this gives
/// ```text
///   b'_i = b_(i + 1)
///   b'_(i + 1) = b_i
///
///   b'_i* = b_(i + 1)* + mu b_i*
///   b'_(i + 1) = (1 - gamma^2 mu^2) b_i* - mu * gamma^2 b_(i + 1)*
///     where gamma^2 = |b_i*|^2 / |b'_i*|^2
///   mu' = gamma^2 mu
/// ```
fn swap_gso_cols<R, V>(ring: R, mut gso: SubmatrixMut<V, El<R>>, i: usize, j: usize)
where
    R: RingStore,
    R::Type: OrderedRing + Field,
    V: AsPointerToSlice<El<R>>,
{
    // numerically very unstable
    assert!(!ring.get_ring().is_approximate());
    assert!(j == i + 1);

    let col_count = gso.col_count();

    // swap the columns
    let (mut col_i, mut col_i1) = gso.reborrow().restrict_cols(i..(i + 2)).split_cols(0..1, 1..2);
    for k in 0..i {
        std::mem::swap(col_i.at_mut(k, 0), col_i1.at_mut(k, 0));
    }

    // re-orthogonalize the triangle `i..(i + 2) x i..(i + 2)`

    // | b_i* |^2
    let bi_star_norm_sqr = ring.clone_el(gso.at(i, i));
    // | b_(i + 1)* |^2
    let bi1_star_norm_sqr = ring.clone_el(gso.at(i + 1, i + 1));
    // mu_(i + 1)i = <b_(i + 1), bi*> / <bi*, bi*>
    let mu = ring.clone_el(gso.at(i, i + 1));
    let mu_sqr = ring.pow(ring.clone_el(&mu), 2);

    let new_bi_star_norm_sqr = ring.add_ref_fst(&bi1_star_norm_sqr, ring.mul_ref(&mu_sqr, &bi_star_norm_sqr));
    // `gamma := |b_(i + 1)*|^2 / |bnew_i*|^2`
    let gamma = if ring.is_leq(&new_bi_star_norm_sqr, &bi1_star_norm_sqr) {
        // in this case, we must have `|bnew_i*|^2 = |b_(i + 1)*|^2`
        ring.one()
    } else {
        ring.div(&bi1_star_norm_sqr, &new_bi_star_norm_sqr)
    };
    let new_bi1_star_norm_sqr = ring.mul_ref(&gamma, &bi_star_norm_sqr);

    // we now update the `mu_ki` resp. `mu_k(i + 1)` by a linear transform;
    // Concretely it is given by the matrix multiplication
    // mu_(j, i)' = mu_(i + 1, i)' mu_(j, i) + gamma mu_(j, i + 1)
    // mu_(j, i + 1)' = mu_(j, i) - mu_(i, i + 1) mu_(j, i + 1)

    // we use `mu_(i + 1, i)' = mu_(i + 1, i) * |b_i*|^2 / |bnew_i*|^2`
    let new_mu = if ring.is_zero(&new_bi_star_norm_sqr) {
        ring.zero()
    } else {
        // this is why we cannot work with floats here, since this might cause a huge error
        ring.div(&ring.mul_ref(&mu, &bi_star_norm_sqr), &new_bi_star_norm_sqr)
    };
    let (row_i, row_i1) = gso.reborrow().restrict_rows(i..(i + 2)).split_rows(0..1, 1..2);
    let row_i = row_i.into_row_mut_at(0);
    let row_i1 = row_i1.into_row_mut_at(0);
    for k in (i + 2)..col_count {
        let new_mu_j_i = ring.add(ring.mul_ref(&new_mu, row_i.at(k)), ring.mul_ref(&gamma, row_i1.at(k)));
        let new_mu_j_i1 = ring.sub_ref_fst(row_i.at(k), ring.mul_ref(&mu, row_i1.at(k)));
        *row_i.at_mut(k) = new_mu_j_i;
        *row_i1.at_mut(k) = new_mu_j_i1;
    }

    *gso.at_mut(i, i + 1) = new_mu;
    *gso.at_mut(i, i) = new_bi_star_norm_sqr;
    *gso.at_mut(i + 1, i + 1) = new_bi1_star_norm_sqr;
}

/// Computes the LDL-decomposition of the given matrix, i.e. writes it as
/// a product `L * D * L^T`, where `D` is diagonal and `L` is lower triangle.
///
/// If the matrix is not invertible, this function cannot proceed, and will return
/// the index of the  column in which on to small values have been detected as `Err()`.
/// The top left square until this index will contain a valid LDL-decomposition of the
/// corresponding square of the input.
///
/// `D` is returned on the diagonal of the matrix, and `L^T` is returned in
/// the upper triangle of the matrix.
fn ldl<R, V>(ring: R, mut matrix: SubmatrixMut<V, El<R>>) -> Result<(), usize>
where
    R: RingStore,
    R::Type: Field,
    V: AsPointerToSlice<El<R>>,
{
    // only the upper triangle part of matrix is used
    assert_eq!(matrix.row_count(), matrix.col_count());
    let n = matrix.row_count();
    for i in 0..n {
        let pivot = ring.clone_el(matrix.at(i, i));
        if ring.is_zero(&pivot) {
            return Err(i);
        }
        let pivot_inv = ring.div(&ring.one(), matrix.at(i, i));
        for j in (i + 1)..n {
            ring.mul_assign_ref(matrix.at_mut(i, j), &pivot_inv);
        }
        for k in (i + 1)..n {
            for l in k..n {
                let subtract = ring.mul_ref_snd(
                    ring.mul_ref(matrix.as_const().at(i, k), matrix.as_const().at(i, l)),
                    &pivot,
                );
                ring.sub_assign(matrix.at_mut(k, l), subtract);
            }
        }
    }
    return Ok(());
}

/// LLL-reduces the lattice basis given by the columns of the given matrix, w.r.t.
/// the norm induced by the given quadratic form.
///
/// The exact restrictions imposed on `B'` are that its columns `b1, ..., bn`
/// are `delta`-LLL-reduced. This means
///  - (size-reduced) `|<bi,bj*>| < |bj*|^2 / 2` whenever `i > j`
///  - (Lovasz-condition) `|bk*|^2 >= delta |b(k - 1)*|^2 - <bk, b(k - 1)*>^2 / |b(k - 1)*|^2`
///
/// Here the `bi*` refer to the Gram-Schmidt orthogonalization of the `bi`.
///
/// The given quadratic form must be positive definite.
///
/// # Internal computations with floating point numbers
///
/// If `disable_float_lll` is not set, this function will first heuristically reduce
/// the matrix using [`crate::algorithms::lll::float::lll_quadratic_form()`].
/// This can significantly speed up the whole computation, as it means we have to do
/// less rational arithmetic, which can be very slow.
#[stability::unstable(feature = "enable")]
pub fn lll<R, I, V1, V2, A>(
    ring: R,
    quadratic_form: Submatrix<V1, El<R>>,
    mut matrix: SubmatrixMut<V2, El<R>>,
    delta: &El<R>,
    allocator: A,
    disable_float_lll: bool,
) where
    R: RingStore<Type = RationalFieldBase<I>> + Copy,
    I: RingStore,
    I::Type: IntegerRing,
    V1: AsPointerToSlice<El<R>>,
    V2: AsPointerToSlice<El<R>>,
    A: Allocator,
{
    let n = matrix.col_count();
    assert_eq!(matrix.row_count(), quadratic_form.col_count());
    assert_eq!(matrix.row_count(), quadratic_form.col_count());
    assert!(ring.is_lt(delta, &ring.one()));
    assert!(ring.is_gt(
        delta,
        &ring.from_fraction(ring.base_ring().one(), ring.base_ring().int_hom().map(4))
    ));

    let mut tmp = OwnedMatrix::zero_in(n, matrix.row_count(), ring, &allocator);
    let mut compute_gram_matrix = |matrix: Submatrix<_, _>, gram_matrix: SubmatrixMut<_, _>| {
        let mut tmp_data = tmp.data_mut().submatrix(0..matrix.col_count(), 0..matrix.row_count());
        STANDARD_MATMUL.matmul(
            TransposableSubmatrix::from(matrix).transpose(),
            TransposableSubmatrix::from(quadratic_form),
            TransposableSubmatrixMut::from(tmp_data.reborrow()),
            ring,
        );
        STANDARD_MATMUL.matmul(
            TransposableSubmatrix::from(tmp_data.as_const()),
            TransposableSubmatrix::from(matrix),
            TransposableSubmatrixMut::from(gram_matrix),
            &ring,
        );
    };

    let mut gram_matrix = OwnedMatrix::zero_in(n, n, ring, &allocator);
    let mut gram_matrix = gram_matrix.data_mut();
    compute_gram_matrix(matrix.as_const(), gram_matrix.reborrow());

    if !disable_float_lll {
        let RR = Real64::RING;
        let QQ_to_RR = RR.can_hom(&ring).unwrap();
        _ = lll_quadratic_form(
            gram_matrix.reborrow(),
            QQ_to_RR,
            &0.999,
            &0.51,
            &mut TransformCols(matrix.reborrow(), ring.get_ring()),
        );
    }

    // we have an outer loop which removes all generated zero vectors from
    // `matrix` and `gram_matrix`, thus shrinking these matrices; It will
    // then perform the actual LLL on the top left part of the `gram_matrix`;
    // we might not be able to run it on the whole `gram_matrix` (with zero
    // vectors removed), since LDL can handle at most a single vector that is
    // in the Q-span of the previous ones
    'remove_zero_vectors: loop {
        compute_gram_matrix(matrix.as_const(), gram_matrix.reborrow());
        while gram_matrix.row_count() > 0 && ring.is_zero(gram_matrix.at(0, 0)) {
            let n = matrix.col_count();
            gram_matrix = gram_matrix.submatrix(1..n, 1..n);
            matrix = matrix.restrict_cols(1..n);
        }
        let n = matrix.col_count();
        let mut gso = match ldl(ring, gram_matrix.reborrow()) {
            Ok(()) => gram_matrix.reborrow(),
            Err(valid_cols) => gram_matrix
                .reborrow()
                .submatrix(0..(valid_cols + 1), 0..(valid_cols + 1)),
        };

        let mut col_ops = TransformCols(matrix.reborrow(), ring.get_ring());
        let mut i = 1;

        #[allow(unused_labels)]
        'lll_main_loop: while i < n {
            assert!(i > 0);
            let (target, gso_part) = gso.reborrow().split_cols(i..(i + 1), 0..i);
            size_reduce(ring, target, i, gso_part.as_const(), &mut col_ops);

            if ring.is_gt(
                &ring.mul_ref_snd(
                    ring.sub_ref_fst(delta, ring.mul_ref(gso.at(i - 1, i), gso.at(i - 1, i))),
                    gso.at(i - 1, i - 1),
                ),
                gso.at(i, i),
            ) {
                col_ops.swap(ring, i - 1, i);
                swap_gso_cols(&ring, gso.reborrow(), i - 1, i);
                i -= 1;
                if i == 0 {
                    // we might have generated a new zero vector, continue
                    // with outer loop where it will be removed
                    continue 'remove_zero_vectors;
                }
            } else {
                i += 1;
            }
        }

        // if we ran the lll main loop on the whole matrix, we are done now
        let just_lll_reduced_dimension = gso.row_count();
        if just_lll_reduced_dimension == gram_matrix.row_count() {
            return;
        }
    }
}

#[cfg(test)]
use std::alloc::Global;

#[cfg(test)]
use test::Bencher;

#[cfg(test)]
use crate::algorithms::lll::{assert_rational_lattice_isomorphic, norm_squared};
#[cfg(test)]
use crate::assert_matrix_eq;
#[cfg(test)]
use crate::primitive_int::StaticRing;

#[cfg(test)]
macro_rules! matrix {
    (# $hom:expr; $num:literal) => {
        ($hom).map($num)
    };
    (# $hom:expr; $num:literal, $den:literal) => {
        ($hom).codomain().div(&($hom).map($num), &($hom).map($den))
    };
    ($hom:expr,$([$($num:literal $(/ $den:literal)?),*]),*) => {
        {
            let ZZ_to_ring = $hom;
            [
                $([$(
                    matrix!(# ZZ_to_ring; $num $(, $den)?)
                ),*]),*
            ]
        }
    };
    ($hom:expr,$(DerefArray::from([$($num:literal $(/ $den:literal)?),*])),*) => {
        {
            let ZZ_to_ring = $hom;
            [
                $(DerefArray::from([$(
                    matrix!(# ZZ_to_ring; $num $(, $den)?)
                ),*])),*
            ]
        }
    };
}

#[test]
fn test_ldl() {
    let ZZ = StaticRing::<i64>::RING;
    let QQ = RationalField::new(ZZ);
    let mut data = matrix!(
        QQ.inclusion(),
        DerefArray::from([1, 2, 1]),
        DerefArray::from([2, 5, 0]),
        DerefArray::from([1, 0, 7])
    );
    let mut matrix = SubmatrixMut::from_2d(&mut data);
    let mut expected = matrix!(QQ.inclusion(), [1, 2, 1], [0, 1, -2], [0, 0, 2]);
    ldl(QQ, matrix.reborrow()).unwrap();

    // only the upper triangle is filled
    expected[1][0] = *matrix.at(1, 0);
    expected[2][0] = *matrix.at(2, 0);
    expected[2][1] = *matrix.at(2, 1);

    assert_matrix_eq!(&QQ, &expected, &matrix);
}

#[test]
fn test_swap_gso_cols() {
    let ZZ = StaticRing::<i64>::RING;
    let QQ = RationalField::new(ZZ);
    let mut matrix = matrix!(
        QQ.inclusion(),
        DerefArray::from([2, 1 / 2, 2 / 5]),
        DerefArray::from([0, 3 / 2, 1 / 4]),
        DerefArray::from([0, 0, 1])
    );
    let expected = matrix!(QQ.inclusion(), [2, 1 / 2, 31 / 80], [0, 3 / 2, 11 / 40], [0, 0, 1]);
    let matrix_view = SubmatrixMut::from_2d(&mut matrix);

    swap_gso_cols(&QQ, matrix_view, 0, 1);

    assert_matrix_eq!(&QQ, &expected, &matrix);
}

#[test]
fn test_lll_2d() {
    let ZZ = StaticRing::<i64>::RING;
    let QQ = RationalField::new(ZZ);
    let original = matrix!(QQ.inclusion(), DerefArray::from([5, 9]), DerefArray::from([11, 20]));
    let mut reduced = original;
    let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 2>, _>::from_2d(&mut reduced);
    lll(
        QQ,
        OwnedMatrix::identity(2, 2, QQ).data(),
        reduced_matrix.reborrow(),
        &QQ.from_fraction(9, 10),
        Global,
        true,
    );

    assert_rational_lattice_isomorphic(
        QQ,
        RationalField::new(BigIntRing::RING),
        Submatrix::from_2d(&original),
        reduced_matrix.as_const(),
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(1),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(0))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(1),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(1))
    );

    let original = matrix!(QQ.inclusion(), DerefArray::from([10, 8]), DerefArray::from([27, 22]));
    let mut reduced = original;
    let mut reduced_matrix = SubmatrixMut::<DerefArray<_, 2>, _>::from_2d(&mut reduced);
    lll(
        QQ,
        OwnedMatrix::identity(2, 2, QQ).data(),
        reduced_matrix.reborrow(),
        &QQ.from_fraction(9, 10),
        Global,
        true,
    );

    assert_rational_lattice_isomorphic(
        QQ,
        RationalField::new(BigIntRing::RING),
        Submatrix::from_2d(&original),
        reduced_matrix.as_const(),
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(4),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(0))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(5),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(1))
    );
}

#[test]
fn test_lll_3d() {
    let ZZ = StaticRing::<i128>::RING;
    let QQ = RationalField::new(ZZ);
    // in this case, the shortest vector is shorter than half the second successive minimum,
    // so LLL will find it (for delta = 0.9 > 0.75)
    let original = matrix!(
        QQ.inclusion(),
        DerefArray::from([72, 0, 0]),
        DerefArray::from([0, 9, 0]),
        DerefArray::from([8432, 7344, 16864])
    );

    let mut reduced = original;
    let mut reduced_matrix = SubmatrixMut::from_2d(&mut reduced);
    lll(
        QQ,
        OwnedMatrix::identity(3, 3, QQ).data(),
        reduced_matrix.reborrow(),
        &QQ.from_fraction(999, 1000),
        Global,
        true,
    );

    assert_rational_lattice_isomorphic(
        QQ,
        RationalField::new(BigIntRing::RING),
        Submatrix::from_2d(&original),
        reduced_matrix.as_const(),
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(144 * 144),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(0))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(72 * 72 + 279 * 279),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(1))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(72 * 72 * 2 + 272 * 272),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(2))
    );
}

#[test]
fn test_lll_generating_set() {
    let ZZ = StaticRing::<i64>::RING;
    let QQ = RationalField::new(ZZ);

    let original = matrix!(
        QQ.inclusion(),
        DerefArray::from([-6, -1, 6, 116, -2]),
        DerefArray::from([-14, -12, 8, 232, -2]),
        DerefArray::from([-10, 2, 12, 0, 2])
    );

    let mut reduced = original;
    let mut reduced_matrix = SubmatrixMut::from_2d(&mut reduced);
    lll(
        QQ,
        OwnedMatrix::identity(3, 3, QQ).data(),
        reduced_matrix.reborrow(),
        &QQ.from_fraction(999, 1000),
        Global,
        true,
    );

    assert_rational_lattice_isomorphic(
        QQ,
        RationalField::new(BigIntRing::RING),
        Submatrix::from_2d(&original),
        reduced_matrix.as_const(),
    );
    assert_el_eq!(QQ, QQ.zero(), norm_squared(QQ, &reduced_matrix.as_const().col_at(0)));
    assert_el_eq!(QQ, QQ.zero(), norm_squared(QQ, &reduced_matrix.as_const().col_at(1)));
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(5),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(2))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(5),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(3))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(12),
        norm_squared(QQ, &reduced_matrix.as_const().col_at(4))
    );

    let ZZ = BigIntRing::RING;
    let QQ = RationalField::new(ZZ);
    let original = matrix!(
        QQ.inclusion().compose(ZZ.can_hom(&StaticRing::<i128>::RING).unwrap()),
        DerefArray::from([-4, 8, -54, -1, 42, 15, -23, -259]),
        DerefArray::from([-3, 10, -36, 18, -48, -473, -1200, -6493]),
        DerefArray::from([5, -13, 62, -15, 17, 398, 1043, 5721]),
        DerefArray::from([-8, 10, -68, 18, -18, -434, -1118, -6126]),
        DerefArray::from([11, -5, 90, 26, -215, -910, -2227, -11637])
    );
    let mut reduced = original.clone();
    let mut reduced_matrix = SubmatrixMut::from_2d(&mut reduced);
    lll(
        &QQ,
        OwnedMatrix::identity(5, 5, &QQ).data(),
        reduced_matrix.reborrow(),
        &QQ.from_fraction(
            int_cast(999, ZZ, StaticRing::<i64>::RING),
            int_cast(1000, ZZ, StaticRing::<i64>::RING),
        ),
        Global,
        true,
    );

    assert_rational_lattice_isomorphic(&QQ, &QQ, Submatrix::from_2d(&original), reduced_matrix.as_const());
    assert_el_eq!(QQ, QQ.zero(), norm_squared(&QQ, &reduced_matrix.as_const().col_at(0)));
    assert_el_eq!(QQ, QQ.zero(), norm_squared(&QQ, &reduced_matrix.as_const().col_at(1)));
    assert_el_eq!(QQ, QQ.zero(), norm_squared(&QQ, &reduced_matrix.as_const().col_at(2)));
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(1),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(3))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(1),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(4))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(1),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(5))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(5),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(6))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(40),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(7))
    );

    let original = matrix!(
        QQ.inclusion().compose(ZZ.can_hom(&StaticRing::<i128>::RING).unwrap()),
        DerefArray::from([
            -60725263117,
            -448122081513,
            -218368759847,
            2100701846793,
            216156377534,
            -3137996709827,
            14835704835919,
            67504381450573
        ]),
        DerefArray::from([
            -1310716961940,
            -9682451257943,
            -4729935920987,
            45413204073392,
            4667627712725,
            -67791459966817,
            320528485599331,
            1458334256347773
        ]),
        DerefArray::from([
            1159398893231,
            8564380015666,
            4183444050825,
            -40168532351902,
            -4128711773154,
            59963582382418,
            -283516375460965,
            -1289940218804617
        ]),
        DerefArray::from([
            -1236320093452,
            -9132639612642,
            -4461079566742,
            42833893239948,
            4402644221490,
            -63942205977848,
            302328006373120,
            1375528654002990
        ]),
        DerefArray::from([
            -2344979577397,
            -17323890604545,
            -8464219959177,
            81256383857324,
            8351008895542,
            -121291584595649,
            573488494361461,
            2609233431319737
        ])
    );
    let mut reduced = original.clone();
    let mut reduced_matrix = SubmatrixMut::from_2d(&mut reduced);
    lll(
        &QQ,
        OwnedMatrix::identity(5, 5, &QQ).data(),
        reduced_matrix.reborrow(),
        &QQ.from_fraction(
            int_cast(999, ZZ, StaticRing::<i64>::RING),
            int_cast(1000, ZZ, StaticRing::<i64>::RING),
        ),
        Global,
        true,
    );

    assert_rational_lattice_isomorphic(&QQ, &QQ, Submatrix::from_2d(&original), reduced_matrix.as_const());
    assert_el_eq!(QQ, QQ.zero(), norm_squared(&QQ, &reduced_matrix.as_const().col_at(0)));
    assert_el_eq!(QQ, QQ.zero(), norm_squared(&QQ, &reduced_matrix.as_const().col_at(1)));
    assert_el_eq!(QQ, QQ.zero(), norm_squared(&QQ, &reduced_matrix.as_const().col_at(2)));
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(1),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(3))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(1),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(4))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(1),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(5))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(5),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(6))
    );
    assert_el_eq!(
        QQ,
        QQ.int_hom().map(40),
        norm_squared(&QQ, &reduced_matrix.as_const().col_at(7))
    );
}

#[bench]
fn bench_lll_10d(bencher: &mut Bencher) {
    let ZZ = BigIntRing::RING;
    let QQ = RationalField::new(ZZ);

    bencher.iter(|| {
        let original = matrix!(
            QQ.inclusion().compose(ZZ.can_hom(&StaticRing::<i128>::RING).unwrap()),
            DerefArray::from([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
            DerefArray::from([0, 1, 0, 0, 0, 0, 0, 0, 0, 0]),
            DerefArray::from([0, 0, 1, 0, 0, 0, 0, 0, 0, 0]),
            DerefArray::from([0, 0, 0, 1, 0, 0, 0, 0, 0, 0]),
            DerefArray::from([0, 0, 0, 0, 1, 0, 0, 0, 0, 0]),
            DerefArray::from([0, 0, 0, 0, 0, 1, 0, 0, 0, 0]),
            DerefArray::from([0, 0, 0, 0, 0, 0, 1, 0, 0, 0]),
            DerefArray::from([2, 2, 2, 2, 0, 0, 1, 4, 0, 0]),
            DerefArray::from([4, 3, 3, 3, 1, 2, 1, 0, 5, 0]),
            DerefArray::from([
                3433883, 14315221, 24549008, 6570781, 32725387, 33674813, 27390657, 15726308, 43003827, 43364304
            ])
        );
        let mut reduced = original.clone();
        let mut reduced_matrix = SubmatrixMut::from_2d(&mut reduced);
        let delta = QQ.from_fraction(
            int_cast(9, ZZ, StaticRing::<i64>::RING),
            int_cast(10, ZZ, StaticRing::<i64>::RING),
        );
        lll(
            &QQ,
            OwnedMatrix::identity(10, 10, &QQ).data(),
            reduced_matrix.reborrow(),
            &delta,
            Global,
            true,
        );

        assert_rational_lattice_isomorphic(&QQ, &QQ, Submatrix::from_2d(&original), reduced_matrix.as_const());
        assert!(QQ.is_geq(
            &QQ.int_hom().map(16 * 16),
            &norm_squared(&QQ, &reduced_matrix.as_const().col_at(0))
        ));
    });
}