Trait Ring 

pub trait Ring: SemiRing<Additive: Invertible> {
    // Provided methods
    fn neg(x: &Self::T) -> Self::T { ... }
    fn sub(x: &Self::T, y: &Self::T) -> Self::T { ... }
    fn sub_assign(x: &mut Self::T, y: &Self::T) { ... }
}

Provided Methods

fn neg(x: &Self::T) -> Self::T

additive inverse: $-$

Examples found in repository

crates/competitive/src/math/floor_sum.rs (line 367)

pub fn floor_sum_polynomial_i64<T, const X: usize, const Y: usize>(
    l: i64,
    r: i64,
    a: i64,
    b: i64,
    m: u64,
) -> [[T; Y]; X]
where
    T: DotProduct + One,
    AddMulOperation<T>: SemiRing<T = T, Additive: Invertible>,
{
    assert!(l <= r);
    assert!(m > 0);

    if a < 0 {
        let mut ans = floor_sum_polynomial_i64::<T, X, Y>(-r + 1, -l + 1, -a, b, m);
        for ans in ans.iter_mut().skip(1).step_by(2) {
            for ans in ans.iter_mut() {
                *ans = AddMulOperation::<T>::neg(ans);
            }
        }
        return ans;
    }

    let add_x = l;
    let n = (r - l) as u64;
    let b = a * add_x + b;

    let add_y = b.div_euclid(m as i64);
    let b = b.rem_euclid(m as i64);
    assert!(a >= 0);
    assert!(b >= 0);
    let data = floor_monoid_product::<FloorSum<AddMulOperation<T>, X, Y>>(
        FloorSum::<AddMulOperation<T>, X, Y>::to_x(),
        FloorSum::<AddMulOperation<T>, X, Y>::to_y(),
        n,
        a as u64,
        b as u64,
        m,
    );

    let offset = FloorSum::<AddMulOperation<T>, X, Y>::offset(add_x, add_y);
    FloorSum::<AddMulOperation<T>, X, Y>::operate(&offset, &data).dp
}

crates/competitive/src/math/matrix.rs (line 215)

    fn eliminate<const DETERMINANT: bool>(&mut self) -> (usize, R::T) {
        let (n, m) = self.shape;
        let mut rank = 0;
        let mut determinant = R::one();
        let mut negative = false;
        for first in (0..m).step_by(64) {
            if rank == n {
                break;
            }
            let end = (first + 64).min(m);
            let start = rank;
            let mut pivots = Vec::new();
            let mut panel = vec![vec![R::zero(); end - first]; end - first];
            for col in first..end {
                if panel[col - first][..col - first]
                    .iter()
                    .any(|x| !R::is_zero(x))
                {
                    for row in &mut self.data[rank..] {
                        let value =
                            R::dot_product(&row[first..col], &panel[col - first][..col - first]);
                        R::sub_assign(&mut row[col], &value);
                    }
                }
                let Some(pivot) = (rank..n).find(|&i| !R::is_zero(&self[i][col])) else {
                    continue;
                };
                if pivot != rank {
                    self.data.swap(rank, pivot);
                    negative = !negative;
                }
                if DETERMINANT {
                    R::mul_assign(&mut determinant, &self[rank][col]);
                }
                let inv = R::inv(&self[rank][col]);
                let row = &mut self.data[rank];
                for c in col + 1..end {
                    let value = R::dot_product(&row[first..col], &panel[c - first][..col - first]);
                    R::sub_assign(&mut row[c], &value);
                    panel[c - first][col - first] = row[c].clone();
                }
                for row in &mut self.data[rank + 1..] {
                    R::mul_assign(&mut row[col], &inv);
                }
                pivots.push(col);
                rank += 1;
                if rank == n {
                    break;
                }
            }
            for i in start..if rank - start < 32 { n } else { rank } {
                let (upper, lower) = self.data.split_at_mut(i);
                let row = &mut lower[0];
                for (j, &col) in pivots[..(i - start).min(pivots.len())].iter().enumerate() {
                    if R::is_zero(&row[col]) {
                        continue;
                    }
                    let factor = R::neg(&row[col]);
                    R::add_scaled_assign(&mut row[end..], &upper[start + j][end..], &factor);
                }
            }
            if rank < n
                && end < m
                && rank - start >= 32
                && self.data[rank..]
                    .iter()
                    .any(|row| pivots.iter().any(|&col| !R::is_zero(&row[col])))
            {
                let lower = Self::new_with((n - rank, rank - start), |i, j| {
                    R::neg(&self[rank + i][pivots[j]])
                });
                let upper = Self::new_with((rank - start, m - end), |i, j| {
                    self[start + i][end + j].clone()
                });
                let update = &lower * &upper;
                for (row, update) in self.data[rank..].iter_mut().zip(update.data) {
                    for (x, y) in row[end..].iter_mut().zip(update) {
                        R::add_assign(x, &y);
                    }
                }
            }
            for (i, &col) in pivots.iter().enumerate() {
                for row in &mut self.data[start + i + 1..] {
                    row[col] = R::zero();
                }
            }
            if DETERMINANT && rank < end {
                return (rank, R::zero());
            }
        }
        if DETERMINANT && negative {
            determinant = R::neg(&determinant);
        }
        (rank, determinant)
    }

    /// f: (row, pivot_row, col)
    pub fn row_reduction_with<F>(&mut self, normalize: bool, mut f: F)
    where
        F: FnMut(usize, usize, usize),
    {
        let (n, m) = self.shape;
        let mut c = 0;
        for r in 0..n {
            loop {
                if c >= m {
                    return;
                }
                if let Some(pivot) = (r..n).find(|&p| !R::is_zero(&self[p][c])) {
                    f(r, pivot, c);
                    self.data.swap(r, pivot);
                    break;
                };
                c += 1;
            }
            let d = R::inv(&self[r][c]);
            if normalize {
                for value in &mut self[r][c..m] {
                    R::mul_assign(value, &d);
                }
            }
            for i in (0..n).filter(|&i| i != r) {
                let mut e = self[i][c].clone();
                if !normalize {
                    R::mul_assign(&mut e, &d);
                }
                for j in c..m {
                    let e = R::mul(&e, &self[r][j]);
                    R::sub_assign(&mut self[i][j], &e);
                }
            }
            c += 1;
        }
    }

    pub fn row_reduction(&mut self, normalize: bool) {
        self.row_reduction_with(normalize, |_, _, _| {});
    }

    pub fn rank(&mut self) -> usize {
        self.eliminate::<false>().0
    }

    pub fn determinant(&mut self) -> R::T {
        assert_eq!(self.shape.0, self.shape.1);
        self.eliminate::<true>().1
    }

    pub fn solve_system_of_linear_equations(
        &self,
        b: &[R::T],
    ) -> Option<SystemOfLinearEquationsSolution<R>> {
        assert_eq!(self.shape.0, b.len());
        let m = self.shape.1;
        let mut a = Self::new_with((self.shape.0, m + 1), |i, j| {
            if j == m {
                b[i].clone()
            } else {
                self[i][j].clone()
            }
        });
        let rank = a.eliminate::<false>().0;
        let mut pivots = Vec::with_capacity(rank);
        let mut b = Vec::with_capacity(rank);
        for row in &a.data[..rank] {
            let c = row.iter().position(|x| !R::is_zero(x)).unwrap();
            if c == m {
                return None;
            }
            pivots.push(c);
            b.push(row[m].clone());
        }

        let mut free = Vec::with_capacity(m - rank);
        let mut pivot = 0;
        for c in 0..m {
            if pivot < rank && pivots[pivot] == c {
                pivot += 1;
            } else {
                free.push(c);
            }
        }
        let mut coefficients: Vec<Vec<_>> = (0..rank)
            .map(|i| free.iter().map(|&c| a[i][c].clone()).collect())
            .collect();
        for k in (0..rank).rev() {
            let c = pivots[k];
            let inv = R::inv(&a[k][c]);
            R::mul_assign(&mut b[k], &inv);
            let pivot_b = b[k].clone();
            let (upper, lower) = coefficients.split_at_mut(k);
            let pivot_coefficients = &mut lower[0];
            for x in pivot_coefficients.iter_mut() {
                R::mul_assign(x, &inv);
            }
            for ((row, value), coefficients) in a.data[..k].iter_mut().zip(&mut b[..k]).zip(upper) {
                if R::is_zero(&row[c]) {
                    continue;
                }
                let factor = row[c].clone();
                row[c] = R::zero();
                R::sub_assign(value, &R::mul(&factor, &pivot_b));
                R::add_scaled_assign(coefficients, pivot_coefficients, &R::neg(&factor));
            }
        }

        let mut particular = vec![R::zero(); m];
        for i in 0..rank {
            particular[pivots[i]] = b[i].clone();
        }
        let mut basis = Vec::with_capacity(free.len());
        for (j, &c) in free.iter().enumerate() {
            let mut vector = vec![R::zero(); m];
            vector[c] = R::one();
            for i in 0..rank {
                vector[pivots[i]] = R::neg(&coefficients[i][j]);
            }
            basis.push(vector);
        }
        Some(SystemOfLinearEquationsSolution { particular, basis })
    }

    pub fn inverse(&self) -> Option<Matrix<R>> {
        assert_eq!(self.shape.0, self.shape.1);
        let n = self.shape.0;
        if n >= 64 {
            let m = n / 2;
            let a = Self::new_with((m, m), |i, j| self[i][j].clone());
            if let Some(mut ai) = a.inverse() {
                let b = Self::new_with((m, n - m), |i, j| self[i][j + m].clone());
                let c = Self::new_with((n - m, m), |i, j| self[i + m][j].clone());
                let mut d = Self::new_with((n - m, n - m), |i, j| self[i + m][j + m].clone());
                let u = &ai * &b;
                let v = &c * &ai;
                d -= &v * &b;
                let di = d.inverse()?;
                let r = &u * &di;
                let t = &di * &v;
                ai += &r * &v;
                let mut inverse = Self::zeros((n, n));
                for i in 0..m {
                    inverse[i][..m].clone_from_slice(&ai[i]);
                    for (x, y) in inverse[i][m..].iter_mut().zip(&r[i]) {
                        *x = R::neg(y);
                    }
                }
                for i in m..n {
                    for (x, y) in inverse[i][..m].iter_mut().zip(&t[i - m]) {
                        *x = R::neg(y);
                    }
                    inverse[i][m..].clone_from_slice(&di[i - m]);
                }
                return Some(inverse);
            }
        }
        let mut a = self.clone();
        let mut inverse = Self::eye((n, n));
        let mut ranges: Vec<_> = (0..n).map(|i| (i, i + 1)).collect();
        for r in 0..n {
            let pivot = (r..n).find(|&i| !R::is_zero(&a[i][r]))?;
            a.data.swap(r, pivot);
            inverse.data.swap(r, pivot);
            ranges.swap(r, pivot);

            let d = R::inv(&a[r][r]);
            for x in &mut a[r][r..] {
                R::mul_assign(x, &d);
            }
            let (left, right) = ranges[r];
            for x in &mut inverse[r][left..right] {
                R::mul_assign(x, &d);
            }

            let (a_upper, a_lower) = a.data.split_at_mut(r + 1);
            let pivot_a = &a_upper[r];
            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r + 1);
            let pivot_inverse = &inverse_upper[r];
            let (ranges_upper, ranges_lower) = ranges.split_at_mut(r + 1);
            let (left, right) = ranges_upper[r];
            for ((a, inverse), range) in a_lower.iter_mut().zip(inverse_lower).zip(ranges_lower) {
                if R::is_zero(&a[r]) {
                    continue;
                }
                let e = a[r].clone();
                a[r] = R::zero();
                R::add_scaled_assign(&mut a[(r + 1)..], &pivot_a[(r + 1)..], &R::neg(&e));
                R::add_scaled_assign(
                    &mut inverse[left..right],
                    &pivot_inverse[left..right],
                    &R::neg(&e),
                );
                range.0 = range.0.min(left);
                range.1 = range.1.max(right);
            }
        }
        for r in (0..n).rev() {
            let (left, right) = ranges[r];
            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r);
            let pivot_inverse = &inverse_lower[0];
            let (ranges_upper, _) = ranges.split_at_mut(r);
            for ((a, inverse), range) in a.data[..r].iter_mut().zip(inverse_upper).zip(ranges_upper)
            {
                if R::is_zero(&a[r]) {
                    continue;
                }
                let e = a[r].clone();
                a[r] = R::zero();
                R::add_scaled_assign(
                    &mut inverse[left..right],
                    &pivot_inverse[left..right],
                    &R::neg(&e),
                );
                range.0 = range.0.min(left);
                range.1 = range.1.max(right);
            }
        }
        Some(inverse)
    }

    pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
        let n = self.shape.0;
        if n == 0 {
            return vec![R::one()];
        }
        assert!(self.data.iter().all(|a| a.len() == n));
        for j in 0..(n - 1) {
            if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
                self.data.swap(j + 1, x);
                self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
                let inv = R::inv(&self[j + 1][j]);
                let mut v = vec![];
                let src = std::mem::take(&mut self[j + 1]);
                for a in self.data[(j + 2)..].iter_mut() {
                    let mul = R::mul(&a[j], &inv);
                    R::add_scaled_assign(&mut a[j..], &src[j..], &R::neg(&mul));
                    v.push(mul);
                }
                self[j + 1] = src;
                for a in self.data.iter_mut() {
                    let v = R::dot_product(&a[(j + 2)..], &v);
                    R::add_assign(&mut a[j + 1], &v);
                }
            }
        }
        // dp[k][j - k] stores [x^k] det(xI - A[..j, ..j]).
        let mut dp: Vec<Vec<R::T>> = (0..=n).map(|i| Vec::with_capacity(n + 1 - i)).collect();
        dp[0].push(R::one());
        for i in 0..n {
            let mut c = vec![R::zero(); i + 1];
            c[i] = R::neg(&self[i][i]);
            let mut mul = R::one();
            for j in (0..i).rev() {
                mul = R::mul(&mul, &self[j + 1][j]);
                c[j] = R::neg(&R::mul(&mul, &self[j][i]));
            }
            for k in (0..=i).rev() {
                let mut value = R::dot_product(&dp[k], &c[k..]);
                if k > 0 {
                    R::add_assign(&mut value, dp[k - 1].last().unwrap());
                }
                dp[k].push(value);
            }
            dp[i + 1].push(R::one());
        }
        dp.into_iter().map(|mut c| c.pop().unwrap()).collect()
    }
}

impl<R> Index<usize> for Matrix<R>
where
    R: SemiRing,
{
    type Output = Vec<R::T>;
    fn index(&self, index: usize) -> &Self::Output {
        &self.data[index]
    }
}

impl<R> IndexMut<usize> for Matrix<R>
where
    R: SemiRing,
{
    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
        &mut self.data[index]
    }
}

impl<R> Index<(usize, usize)> for Matrix<R>
where
    R: SemiRing,
{
    type Output = R::T;
    fn index(&self, index: (usize, usize)) -> &Self::Output {
        &self.data[index.0][index.1]
    }
}

impl<R> IndexMut<(usize, usize)> for Matrix<R>
where
    R: SemiRing,
{
    fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
        &mut self.data[index.0][index.1]
    }
}

macro_rules! impl_matrix_pairwise_binop {
    ($imp:ident, $method:ident, $imp_assign:ident, $method_assign:ident $(where [$($clauses:tt)*])?) => {
        impl<R> $imp_assign for Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            fn $method_assign(&mut self, rhs: Self) {
                self.pairwise_assign(&rhs, |a, b| R::$method_assign(a, b));
            }
        }
        impl<R> $imp_assign<&Matrix<R>> for Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            fn $method_assign(&mut self, rhs: &Self) {
                self.pairwise_assign(rhs, |a, b| R::$method_assign(a, b));
            }
        }
        impl<R> $imp for Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            type Output = Matrix<R>;
            fn $method(mut self, rhs: Self) -> Self::Output {
                self.$method_assign(rhs);
                self
            }
        }
        impl<R> $imp<&Matrix<R>> for Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            type Output = Matrix<R>;
            fn $method(mut self, rhs: &Self) -> Self::Output {
                self.$method_assign(rhs);
                self
            }
        }
        impl<R> $imp<Matrix<R>> for &Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            type Output = Matrix<R>;
            fn $method(self, mut rhs: Matrix<R>) -> Self::Output {
                rhs.pairwise_assign(self, |a, b| *a = R::$method(b, a));
                rhs
            }
        }
        impl<R> $imp<&Matrix<R>> for &Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            type Output = Matrix<R>;
            fn $method(self, rhs: &Matrix<R>) -> Self::Output {
                let mut this = self.clone();
                this.$method_assign(rhs);
                this
            }
        }
    };
}

impl_matrix_pairwise_binop!(Add, add, AddAssign, add_assign);
impl_matrix_pairwise_binop!(Sub, sub, SubAssign, sub_assign where [R: SemiRing<Additive: Invertible>]);

impl<R> Mul for Matrix<R>
where
    R: SemiRing,
{
    type Output = Matrix<R>;
    fn mul(self, rhs: Self) -> Self::Output {
        (&self).mul(&rhs)
    }
}
impl<R> Mul<&Matrix<R>> for Matrix<R>
where
    R: SemiRing,
{
    type Output = Matrix<R>;
    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
        (&self).mul(rhs)
    }
}
impl<R> Mul<Matrix<R>> for &Matrix<R>
where
    R: SemiRing,
{
    type Output = Matrix<R>;
    fn mul(self, rhs: Matrix<R>) -> Self::Output {
        self.mul(&rhs)
    }
}
impl<R> Mul<&Matrix<R>> for &Matrix<R>
where
    R: SemiRing,
{
    type Output = Matrix<R>;
    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
        assert_eq!(self.shape.1, rhs.shape.0);
        if let Some(data) = R::try_matrix_product(&self.data, &rhs.data) {
            return Matrix::from_vec(data);
        }
        let rhs = rhs.transpose();
        Matrix::new_with((self.shape.0, rhs.shape.0), |i, j| {
            R::dot_product(&self[i], &rhs[j])
        })
    }
}

fn strassen_rec<R: Ring>(
    a: &[R::T],
    b: &[R::T],
    c: &mut [R::T],
    shape: (usize, usize, usize),
    stride_a: usize,
    stride_b: usize,
) {
    let (n, m, p) = shape;
    fn add_block<R: Ring>(
        a: &[R::T],
        b: &[R::T],
        out: &mut [R::T],
        n: usize,
        stride_a: usize,
        stride_b: usize,
    ) {
        for ((a, b), c) in a
            .chunks(stride_a)
            .zip(b.chunks(stride_b))
            .zip(out.chunks_exact_mut(n))
        {
            for ((a, b), c) in a.iter().zip(b.iter()).zip(c.iter_mut()) {
                *c = R::add(a, b);
            }
        }
    }

    fn sub_block<R: Ring>(
        a: &[R::T],
        b: &[R::T],
        out: &mut [R::T],
        n: usize,
        stride_a: usize,
        stride_b: usize,
    ) {
        for ((a, b), c) in a
            .chunks(stride_a)
            .zip(b.chunks(stride_b))
            .zip(out.chunks_exact_mut(n))
        {
            for ((a, b), c) in a.iter().zip(b.iter()).zip(c.iter_mut()) {
                *c = R::sub(a, b);
            }
        }
    }

    if n.min(m).min(p) <= 128 {
        let transposed: Vec<_> = (0..p)
            .flat_map(|j| (0..m).map(move |i| b[i * stride_b + j].clone()))
            .collect();
        for (a, c) in a.chunks(stride_a).zip(c.chunks_exact_mut(p)) {
            for (b, c) in transposed.chunks_exact(m).zip(c) {
                *c = R::dot_product(&a[..m], b);
            }
        }
        return;
    }
    let (h, k, w) = (n / 2, m / 2, p / 2);
    let a11 = 0;
    let a12 = k;
    let a21 = h * stride_a;
    let a22 = a21 + k;
    let b11 = 0;
    let b12 = w;
    let b21 = k * stride_b;
    let b22 = b21 + w;

    let block = h * w;
    let mut buf = vec![R::zero(); h * k + k * w + block * 7];
    let (s1, rest) = buf.split_at_mut(h * k);
    let (s2, m_buf) = rest.split_at_mut(k * w);
    let (m1, rest) = m_buf.split_at_mut(block);
    let (m2, rest) = rest.split_at_mut(block);
    let (m3, rest) = rest.split_at_mut(block);
    let (m4, rest) = rest.split_at_mut(block);
    let (m5, rest) = rest.split_at_mut(block);
    let (m6, m7) = rest.split_at_mut(block);

    // (A11 + A22)(B11 + B22)
    add_block::<R>(&a[a11..], &a[a22..], s1, k, stride_a, stride_a);
    add_block::<R>(&b[b11..], &b[b22..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(s1, s2, m1, (h, k, w), k, w);

    // (A21 + A22) B11
    add_block::<R>(&a[a21..], &a[a22..], s1, k, stride_a, stride_a);
    strassen_rec::<R>(s1, &b[b11..], m2, (h, k, w), k, stride_b);

    // A11 (B12 - B22)
    sub_block::<R>(&b[b12..], &b[b22..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(&a[a11..], s2, m3, (h, k, w), stride_a, w);

    // A22 (B21 - B11)
    sub_block::<R>(&b[b21..], &b[b11..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(&a[a22..], s2, m4, (h, k, w), stride_a, w);

    // (A11 + A12) B22
    add_block::<R>(&a[a11..], &a[a12..], s1, k, stride_a, stride_a);
    strassen_rec::<R>(s1, &b[b22..], m5, (h, k, w), k, stride_b);

    // (A21 - A11)(B11 + B12)
    sub_block::<R>(&a[a21..], &a[a11..], s1, k, stride_a, stride_a);
    add_block::<R>(&b[b11..], &b[b12..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(s1, s2, m6, (h, k, w), k, w);

    // (A12 - A22)(B21 + B22)
    sub_block::<R>(&a[a12..], &a[a22..], s1, k, stride_a, stride_a);
    add_block::<R>(&b[b21..], &b[b22..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(s1, s2, m7, (h, k, w), k, w);

    let c11 = 0;
    let c12 = w;
    let c21 = h * p;
    let c22 = c21 + w;
    for ((((m1, m4), m5), m7), c) in m1
        .iter()
        .zip(m4.iter())
        .zip(m5.iter())
        .zip(m7.iter())
        .zip(c[c11..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::add(m1, m4);
        R::sub_assign(c, m5);
        R::add_assign(c, m7);
    }
    for ((m3, m5), c) in m3
        .iter()
        .zip(m5.iter())
        .zip(c[c12..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::add(m3, m5);
    }
    for ((m2, m4), c) in m2
        .iter()
        .zip(m4.iter())
        .zip(c[c21..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::add(m2, m4);
    }
    for ((((m1, m2), m3), m6), c) in m1
        .iter()
        .zip(m2.iter())
        .zip(m3.iter())
        .zip(m6.iter())
        .zip(c[c22..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::sub(m1, m2);
        R::add_assign(c, m3);
        R::add_assign(c, m6);
    }
}

impl<R> Matrix<R>
where
    R: Ring,
{
    pub fn mul_strassen(&self, rhs: &Matrix<R>) -> Matrix<R> {
        assert_eq!(self.shape.1, rhs.shape.0);
        if let Some(data) = R::try_matrix_product(&self.data, &rhs.data) {
            return Matrix::from_vec(data);
        }
        let (n, m) = self.shape;
        let p = rhs.shape.1;
        if n == 0 || m == 0 || p == 0 {
            return Matrix::zeros((n, p));
        }
        let split = n.min(m).min(p).div_ceil(128).next_power_of_two();
        if split <= 2 {
            return self * rhs;
        }
        let rows = n.div_ceil(split) * split;
        let inner = m.div_ceil(split) * split;
        let cols = p.div_ceil(split) * split;
        let mut a = vec![R::zero(); rows * inner];
        for (a, data) in a.chunks_exact_mut(inner).zip(&self.data) {
            a[..m].clone_from_slice(data);
        }
        let mut b = vec![R::zero(); inner * cols];
        for (b, data) in b.chunks_exact_mut(cols).zip(&rhs.data) {
            b[..p].clone_from_slice(data);
        }
        let mut c = vec![R::zero(); rows * cols];
        strassen_rec::<R>(&a, &b, &mut c, (rows, inner, cols), inner, cols);
        let mut res = Matrix::zeros((n, p));
        for (data, c) in res.data.iter_mut().zip(c.chunks_exact(cols)) {
            data.clone_from_slice(&c[..p]);
        }
        res
    }
}

impl<R> MulAssign<&R::T> for Matrix<R>
where
    R: SemiRing,
{
    fn mul_assign(&mut self, rhs: &R::T) {
        for i in 0..self.shape.0 {
            for j in 0..self.shape.1 {
                R::mul_assign(&mut self[(i, j)], rhs);
            }
        }
    }
}

impl<R> Neg for Matrix<R>
where
    R: SemiRing<Additive: Invertible>,
{
    type Output = Self;

    fn neg(self) -> Self::Output {
        self.map(|x| R::neg(x))
    }
}

impl<R> Neg for &Matrix<R>
where
    R: SemiRing<Additive: Invertible>,
{
    type Output = Matrix<R>;

    fn neg(self) -> Self::Output {
        self.map(|x| R::neg(x))
    }

crates/competitive/src/math/black_box_matrix.rs (line 166)

    pub fn determinant(&self) -> R::T {
        assert_eq!(self.shape.0, self.shape.1);
        let n = self.shape.0;
        let mut columns = vec![Vec::<(usize, R::T)>::new(); n];
        for &(i, j, ref value) in &self.nonzero {
            columns[j].push((i, value.clone()));
        }
        let mut degrees = vec![0; n];
        for column in &mut columns {
            column.sort_unstable_by_key(|&(i, _)| i);
            let mut merged: Vec<(usize, R::T)> = Vec::with_capacity(column.len());
            for (i, value) in column.drain(..) {
                if let Some((_, x)) = merged.last_mut().filter(|(last, _)| *last == i) {
                    R::add_assign(x, &value);
                } else {
                    merged.push((i, value));
                }
            }
            merged.retain(|(i, value)| {
                if R::is_zero(value) {
                    false
                } else {
                    degrees[*i] += 1;
                    true
                }
            });
            *column = merged;
        }
        let mut order: Vec<_> = (0..n).collect();
        order.sort_unstable_by_key(|&j| columns[j].len());
        let mut lower: Vec<Vec<(usize, R::T)>> = Vec::with_capacity(n);
        let mut pivots: Vec<Option<usize>> = vec![None; n];
        let mut x = vec![R::zero(); n];
        let mut seen = vec![0; n];
        let mut stack = Vec::new();
        let mut support = Vec::new();
        let mut determinant = R::one();
        for (k, &j) in order.iter().enumerate() {
            support.clear();
            for &(i, _) in &columns[j] {
                if seen[i] == k + 1 {
                    continue;
                }
                seen[i] = k + 1;
                x[i] = R::zero();
                stack.push((i, 0));
                while let Some((i, next)) = stack.last_mut() {
                    if let Some(pivot) = pivots[*i]
                        && *next < lower[pivot].len()
                    {
                        let row = lower[pivot][*next].0;
                        *next += 1;
                        if seen[row] != k + 1 {
                            seen[row] = k + 1;
                            x[row] = R::zero();
                            stack.push((row, 0));
                        }
                        continue;
                    }
                    support.push(*i);
                    stack.pop();
                }
            }
            for &(i, ref value) in &columns[j] {
                x[i] = value.clone();
            }
            let mut pivot = None;
            for &i in support.iter().rev() {
                if let Some(p) = pivots[i] {
                    let factor = R::neg(&x[i]);
                    for &(row, ref value) in &lower[p] {
                        R::add_assign(&mut x[row], &R::mul(&factor, value));
                    }
                } else if !R::is_zero(&x[i]) && pivot.is_none_or(|p| degrees[i] < degrees[p]) {
                    pivot = Some(i);
                }
            }
            let Some(pivot) = pivot else { return R::zero() };
            R::mul_assign(&mut determinant, &x[pivot]);
            let inv = R::inv(&x[pivot]);
            pivots[pivot] = Some(k);
            lower.push(
                support
                    .iter()
                    .filter(|&&i| pivots[i].is_none() && !R::is_zero(&x[i]))
                    .map(|&i| (i, R::mul(&x[i], &inv)))
                    .collect(),
            );
        }
        for mut permutation in [order, pivots.into_iter().map(Option::unwrap).collect()] {
            for i in 0..n {
                while permutation[i] != i {
                    let j = permutation[i];
                    permutation.swap(i, j);
                    determinant = R::neg(&determinant);
                }
            }
        }
        determinant
    }

crates/competitive/src/algorithm/automata_learning.rs (line 562)

    pub fn train_sample(&mut self, sample: &[usize]) -> bool {
        let Some((prefix, suffix)) = self.split_sample(sample) else {
            return false;
        };
        self.prefixes.push(prefix);
        self.suffixes.push(suffix);
        let n = self.inv_h.shape.0;
        let prefix = &self.prefixes[n];
        let suffix = &self.suffixes[n];
        let u = Matrix::<F>::new_with((n, 1), |i, _| {
            self.automaton.behavior(
                self.prefixes[i]
                    .iter()
                    .cloned()
                    .chain(suffix.iter().cloned()),
            )
        });
        let v = Matrix::<F>::new_with((1, n), |_, j| {
            self.automaton.behavior(
                prefix
                    .iter()
                    .cloned()
                    .chain(self.suffixes[j].iter().cloned()),
            )
        });
        let w = Matrix::<F>::new_with((1, 1), |_, _| {
            self.automaton
                .behavior(prefix.iter().cloned().chain(suffix.iter().cloned()))
        });
        let t = &self.inv_h * &u;
        let s = &v * &self.inv_h;
        let d = F::inv(&(&w - &(&v * &t))[0][0]);
        let dh = &t * &s;
        for i in 0..n {
            for j in 0..n {
                F::add_assign(&mut self.inv_h[i][j], &F::mul(&dh[i][j], &d));
            }
        }
        self.inv_h
            .add_col_with(|i, _| F::neg(&F::mul(&t[i][0], &d)));
        self.inv_h.add_row_with(|_, j| {
            if j != n {
                F::neg(&F::mul(&s[0][j], &d))
            } else {
                d.clone()
            }
        });

        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
            let b = &(&self.nh[x] * &t) * &s;
            for i in 0..n {
                for j in 0..n {
                    F::add_assign(&mut transition[i][j], &F::mul(&b[i][j], &d));
                }
            }
        }
        for (x, nh) in self.nh.iter_mut().enumerate() {
            nh.add_col_with(|i, j| {
                self.automaton.behavior(
                    self.prefixes[i]
                        .iter()
                        .cloned()
                        .chain([x])
                        .chain(self.suffixes[j].iter().cloned()),
                )
            });
            nh.add_row_with(|i, j| {
                self.automaton.behavior(
                    self.prefixes[i]
                        .iter()
                        .cloned()
                        .chain([x])
                        .chain(self.suffixes[j].iter().cloned()),
                )
            });
        }
        self.wfa
            .initial_weights
            .add_col_with(|_, _| if n == 0 { F::one() } else { F::zero() });
        self.wfa
            .final_weights
            .add_row_with(|_, _| self.automaton.behavior(prefix.iter().cloned()));
        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
            transition.add_col_with(|_, _| F::zero());
            transition.add_row_with(|_, _| F::zero());
            for i in 0..=n {
                for j in 0..=n {
                    if i == n || j == n {
                        for k in 0..=n {
                            if i != n && j != n && k != n {
                                continue;
                            }
                            F::add_assign(
                                &mut transition[i][k],
                                &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
                            );
                        }
                    } else {
                        let k = n;
                        F::add_assign(
                            &mut transition[i][k],
                            &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
                        );
                    }
                }
            }
        }
        true
    }

fn sub(x: &Self::T, y: &Self::T) -> Self::T

additive right inversed operaion: $-$

Examples found in repository

crates/competitive/src/string/rolling_hash.rs (line 527)

    fn muln_sub(&mut self, l: &R::T, r: &R::T, n: usize) -> R::T {
        if let Some(pow) = self.pow.get(n) {
            R::sub(r, &R::mul(l, pow))
        } else {
            let pow = <R::Multiplicative as Monoid>::pow(self.base.clone(), n);
            R::sub(r, &R::mul(l, &pow))
        }
    }

crates/competitive/src/math/matrix.rs (line 721)

fn strassen_rec<R: Ring>(
    a: &[R::T],
    b: &[R::T],
    c: &mut [R::T],
    shape: (usize, usize, usize),
    stride_a: usize,
    stride_b: usize,
) {
    let (n, m, p) = shape;
    fn add_block<R: Ring>(
        a: &[R::T],
        b: &[R::T],
        out: &mut [R::T],
        n: usize,
        stride_a: usize,
        stride_b: usize,
    ) {
        for ((a, b), c) in a
            .chunks(stride_a)
            .zip(b.chunks(stride_b))
            .zip(out.chunks_exact_mut(n))
        {
            for ((a, b), c) in a.iter().zip(b.iter()).zip(c.iter_mut()) {
                *c = R::add(a, b);
            }
        }
    }

    fn sub_block<R: Ring>(
        a: &[R::T],
        b: &[R::T],
        out: &mut [R::T],
        n: usize,
        stride_a: usize,
        stride_b: usize,
    ) {
        for ((a, b), c) in a
            .chunks(stride_a)
            .zip(b.chunks(stride_b))
            .zip(out.chunks_exact_mut(n))
        {
            for ((a, b), c) in a.iter().zip(b.iter()).zip(c.iter_mut()) {
                *c = R::sub(a, b);
            }
        }
    }

    if n.min(m).min(p) <= 128 {
        let transposed: Vec<_> = (0..p)
            .flat_map(|j| (0..m).map(move |i| b[i * stride_b + j].clone()))
            .collect();
        for (a, c) in a.chunks(stride_a).zip(c.chunks_exact_mut(p)) {
            for (b, c) in transposed.chunks_exact(m).zip(c) {
                *c = R::dot_product(&a[..m], b);
            }
        }
        return;
    }
    let (h, k, w) = (n / 2, m / 2, p / 2);
    let a11 = 0;
    let a12 = k;
    let a21 = h * stride_a;
    let a22 = a21 + k;
    let b11 = 0;
    let b12 = w;
    let b21 = k * stride_b;
    let b22 = b21 + w;

    let block = h * w;
    let mut buf = vec![R::zero(); h * k + k * w + block * 7];
    let (s1, rest) = buf.split_at_mut(h * k);
    let (s2, m_buf) = rest.split_at_mut(k * w);
    let (m1, rest) = m_buf.split_at_mut(block);
    let (m2, rest) = rest.split_at_mut(block);
    let (m3, rest) = rest.split_at_mut(block);
    let (m4, rest) = rest.split_at_mut(block);
    let (m5, rest) = rest.split_at_mut(block);
    let (m6, m7) = rest.split_at_mut(block);

    // (A11 + A22)(B11 + B22)
    add_block::<R>(&a[a11..], &a[a22..], s1, k, stride_a, stride_a);
    add_block::<R>(&b[b11..], &b[b22..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(s1, s2, m1, (h, k, w), k, w);

    // (A21 + A22) B11
    add_block::<R>(&a[a21..], &a[a22..], s1, k, stride_a, stride_a);
    strassen_rec::<R>(s1, &b[b11..], m2, (h, k, w), k, stride_b);

    // A11 (B12 - B22)
    sub_block::<R>(&b[b12..], &b[b22..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(&a[a11..], s2, m3, (h, k, w), stride_a, w);

    // A22 (B21 - B11)
    sub_block::<R>(&b[b21..], &b[b11..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(&a[a22..], s2, m4, (h, k, w), stride_a, w);

    // (A11 + A12) B22
    add_block::<R>(&a[a11..], &a[a12..], s1, k, stride_a, stride_a);
    strassen_rec::<R>(s1, &b[b22..], m5, (h, k, w), k, stride_b);

    // (A21 - A11)(B11 + B12)
    sub_block::<R>(&a[a21..], &a[a11..], s1, k, stride_a, stride_a);
    add_block::<R>(&b[b11..], &b[b12..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(s1, s2, m6, (h, k, w), k, w);

    // (A12 - A22)(B21 + B22)
    sub_block::<R>(&a[a12..], &a[a22..], s1, k, stride_a, stride_a);
    add_block::<R>(&b[b21..], &b[b22..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(s1, s2, m7, (h, k, w), k, w);

    let c11 = 0;
    let c12 = w;
    let c21 = h * p;
    let c22 = c21 + w;
    for ((((m1, m4), m5), m7), c) in m1
        .iter()
        .zip(m4.iter())
        .zip(m5.iter())
        .zip(m7.iter())
        .zip(c[c11..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::add(m1, m4);
        R::sub_assign(c, m5);
        R::add_assign(c, m7);
    }
    for ((m3, m5), c) in m3
        .iter()
        .zip(m5.iter())
        .zip(c[c12..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::add(m3, m5);
    }
    for ((m2, m4), c) in m2
        .iter()
        .zip(m4.iter())
        .zip(c[c21..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::add(m2, m4);
    }
    for ((((m1, m2), m3), m6), c) in m1
        .iter()
        .zip(m2.iter())
        .zip(m3.iter())
        .zip(m6.iter())
        .zip(c[c22..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::sub(m1, m2);
        R::add_assign(c, m3);
        R::add_assign(c, m6);
    }
}

crates/competitive/src/math/quotient_array.rs (line 102)

    pub fn min_25_sieve<R>(&self, mut f: impl FnMut(u64, u32) -> T) -> Self
    where
        T: Clone + One,
        R: Ring<T = T, Additive: Invertible>,
    {
        let mut dp = self.clone();
        let max_n = self.isqrtn as u32;
        with_prime_list(max_n, |pl| {
            for p in pl.primes_lte(max_n).rev() {
                let p = u64::from(p);
                let k = self.quotient_index(p);
                for (i, q) in Self::index_iter(self.n, self.isqrtn).enumerate() {
                    let mut pc = p;
                    if pc * p > q {
                        break;
                    }
                    let mut c = 1;
                    while q / p >= pc {
                        let x = R::mul(&f(p, c), &(R::sub(&dp[q / pc], &self.data[k])));
                        let x = R::add(&x, &f(p, c + 1));
                        dp.data[i] = R::add(&dp.data[i], &x);
                        c += 1;
                        pc *= p;
                    }
                }
            }
        });
        for x in &mut dp.data {
            *x = R::add(x, &T::one());
        }
        dp
    }

fn sub_assign(x: &mut Self::T, y: &Self::T)

Examples found in repository

crates/competitive/src/math/subset_convolve.rs (line 95)

    fn inverse_transform((mut f, width): Self::F, len: usize) -> Self::T {
        let k = width - 1;
        for bit in 0..k {
            let half = 1 << bit;
            for base in (0..len).step_by(half * 2) {
                for lower in base..base + half {
                    let upper = lower + half;
                    let rank = lower.count_ones() as usize;
                    let (lower_rows, upper_rows) = f.split_at_mut(upper * width);
                    let lower_row = &lower_rows[lower * width + rank..lower * width + width];
                    let upper_row = &mut upper_rows[rank..width];
                    for (upper, lower) in upper_row.iter_mut().zip(lower_row) {
                        R::sub_assign(upper, lower);
                    }
                }
            }
        }
        Self::diagonal(f, width)
    }

    fn multiply(f: &mut Self::F, g: &Self::F) {
        let (f, width) = f;
        let (g, _) = g;
        let mut right = vec![R::zero(); *width];
        let mut output = vec![R::zero(); *width];
        for (i, f) in f.chunks_exact_mut(*width).enumerate() {
            let rank = i.count_ones() as usize;
            let g = &g[i * *width..(i + 1) * *width];
            let end = Self::multiply_row(f, g, &mut right, &mut output, rank);
            f[rank..=end].clone_from_slice(&output[rank..=end]);
        }
    }

    fn convolve(a: Self::T, b: Self::T) -> Self::T {
        assert_eq!(a.len(), b.len());
        let len = a.len();
        let same = a == b;
        let (mut x, width) = Self::ranked(a, len);
        let (mut y, _) = if same {
            (x.clone(), width)
        } else {
            Self::ranked(b, len)
        };
        let mut right = vec![R::zero(); width];
        let mut output = vec![R::zero(); width];
        for i in 0..len {
            for bit in (0..(i | len).trailing_zeros() as usize).rev() {
                let half = width << bit;
                let start = i * width;
                let (lower, upper) = x[start..start + half * 2].split_at_mut(half);
                for (upper, lower) in upper.iter_mut().zip(lower) {
                    R::add_assign(upper, lower);
                }
                let (lower, upper) = y[start..start + half * 2].split_at_mut(half);
                for (upper, lower) in upper.iter_mut().zip(lower) {
                    R::add_assign(upper, lower);
                }
            }

            let rank = i.count_ones() as usize;
            let start = i * width;
            let x_row = &x[start..start + width];
            let y_row = &y[start..start + width];
            output.fill(R::zero());
            Self::multiply_row(x_row, y_row, &mut right, &mut output, rank);
            x[start..start + width].clone_from_slice(&output);

            for bit in 0..i.trailing_ones() as usize {
                let end = (i + 1) * width;
                let half = width << bit;
                let (lower, upper) = x[end - half * 2..end].split_at_mut(half);
                for (upper, lower) in upper.iter_mut().zip(lower) {
                    R::sub_assign(upper, lower);
                }
            }
        }
        Self::diagonal(x, width)
    }

crates/competitive/src/math/matrix.rs (line 179)

    fn eliminate<const DETERMINANT: bool>(&mut self) -> (usize, R::T) {
        let (n, m) = self.shape;
        let mut rank = 0;
        let mut determinant = R::one();
        let mut negative = false;
        for first in (0..m).step_by(64) {
            if rank == n {
                break;
            }
            let end = (first + 64).min(m);
            let start = rank;
            let mut pivots = Vec::new();
            let mut panel = vec![vec![R::zero(); end - first]; end - first];
            for col in first..end {
                if panel[col - first][..col - first]
                    .iter()
                    .any(|x| !R::is_zero(x))
                {
                    for row in &mut self.data[rank..] {
                        let value =
                            R::dot_product(&row[first..col], &panel[col - first][..col - first]);
                        R::sub_assign(&mut row[col], &value);
                    }
                }
                let Some(pivot) = (rank..n).find(|&i| !R::is_zero(&self[i][col])) else {
                    continue;
                };
                if pivot != rank {
                    self.data.swap(rank, pivot);
                    negative = !negative;
                }
                if DETERMINANT {
                    R::mul_assign(&mut determinant, &self[rank][col]);
                }
                let inv = R::inv(&self[rank][col]);
                let row = &mut self.data[rank];
                for c in col + 1..end {
                    let value = R::dot_product(&row[first..col], &panel[c - first][..col - first]);
                    R::sub_assign(&mut row[c], &value);
                    panel[c - first][col - first] = row[c].clone();
                }
                for row in &mut self.data[rank + 1..] {
                    R::mul_assign(&mut row[col], &inv);
                }
                pivots.push(col);
                rank += 1;
                if rank == n {
                    break;
                }
            }
            for i in start..if rank - start < 32 { n } else { rank } {
                let (upper, lower) = self.data.split_at_mut(i);
                let row = &mut lower[0];
                for (j, &col) in pivots[..(i - start).min(pivots.len())].iter().enumerate() {
                    if R::is_zero(&row[col]) {
                        continue;
                    }
                    let factor = R::neg(&row[col]);
                    R::add_scaled_assign(&mut row[end..], &upper[start + j][end..], &factor);
                }
            }
            if rank < n
                && end < m
                && rank - start >= 32
                && self.data[rank..]
                    .iter()
                    .any(|row| pivots.iter().any(|&col| !R::is_zero(&row[col])))
            {
                let lower = Self::new_with((n - rank, rank - start), |i, j| {
                    R::neg(&self[rank + i][pivots[j]])
                });
                let upper = Self::new_with((rank - start, m - end), |i, j| {
                    self[start + i][end + j].clone()
                });
                let update = &lower * &upper;
                for (row, update) in self.data[rank..].iter_mut().zip(update.data) {
                    for (x, y) in row[end..].iter_mut().zip(update) {
                        R::add_assign(x, &y);
                    }
                }
            }
            for (i, &col) in pivots.iter().enumerate() {
                for row in &mut self.data[start + i + 1..] {
                    row[col] = R::zero();
                }
            }
            if DETERMINANT && rank < end {
                return (rank, R::zero());
            }
        }
        if DETERMINANT && negative {
            determinant = R::neg(&determinant);
        }
        (rank, determinant)
    }

    /// f: (row, pivot_row, col)
    pub fn row_reduction_with<F>(&mut self, normalize: bool, mut f: F)
    where
        F: FnMut(usize, usize, usize),
    {
        let (n, m) = self.shape;
        let mut c = 0;
        for r in 0..n {
            loop {
                if c >= m {
                    return;
                }
                if let Some(pivot) = (r..n).find(|&p| !R::is_zero(&self[p][c])) {
                    f(r, pivot, c);
                    self.data.swap(r, pivot);
                    break;
                };
                c += 1;
            }
            let d = R::inv(&self[r][c]);
            if normalize {
                for value in &mut self[r][c..m] {
                    R::mul_assign(value, &d);
                }
            }
            for i in (0..n).filter(|&i| i != r) {
                let mut e = self[i][c].clone();
                if !normalize {
                    R::mul_assign(&mut e, &d);
                }
                for j in c..m {
                    let e = R::mul(&e, &self[r][j]);
                    R::sub_assign(&mut self[i][j], &e);
                }
            }
            c += 1;
        }
    }

    pub fn row_reduction(&mut self, normalize: bool) {
        self.row_reduction_with(normalize, |_, _, _| {});
    }

    pub fn rank(&mut self) -> usize {
        self.eliminate::<false>().0
    }

    pub fn determinant(&mut self) -> R::T {
        assert_eq!(self.shape.0, self.shape.1);
        self.eliminate::<true>().1
    }

    pub fn solve_system_of_linear_equations(
        &self,
        b: &[R::T],
    ) -> Option<SystemOfLinearEquationsSolution<R>> {
        assert_eq!(self.shape.0, b.len());
        let m = self.shape.1;
        let mut a = Self::new_with((self.shape.0, m + 1), |i, j| {
            if j == m {
                b[i].clone()
            } else {
                self[i][j].clone()
            }
        });
        let rank = a.eliminate::<false>().0;
        let mut pivots = Vec::with_capacity(rank);
        let mut b = Vec::with_capacity(rank);
        for row in &a.data[..rank] {
            let c = row.iter().position(|x| !R::is_zero(x)).unwrap();
            if c == m {
                return None;
            }
            pivots.push(c);
            b.push(row[m].clone());
        }

        let mut free = Vec::with_capacity(m - rank);
        let mut pivot = 0;
        for c in 0..m {
            if pivot < rank && pivots[pivot] == c {
                pivot += 1;
            } else {
                free.push(c);
            }
        }
        let mut coefficients: Vec<Vec<_>> = (0..rank)
            .map(|i| free.iter().map(|&c| a[i][c].clone()).collect())
            .collect();
        for k in (0..rank).rev() {
            let c = pivots[k];
            let inv = R::inv(&a[k][c]);
            R::mul_assign(&mut b[k], &inv);
            let pivot_b = b[k].clone();
            let (upper, lower) = coefficients.split_at_mut(k);
            let pivot_coefficients = &mut lower[0];
            for x in pivot_coefficients.iter_mut() {
                R::mul_assign(x, &inv);
            }
            for ((row, value), coefficients) in a.data[..k].iter_mut().zip(&mut b[..k]).zip(upper) {
                if R::is_zero(&row[c]) {
                    continue;
                }
                let factor = row[c].clone();
                row[c] = R::zero();
                R::sub_assign(value, &R::mul(&factor, &pivot_b));
                R::add_scaled_assign(coefficients, pivot_coefficients, &R::neg(&factor));
            }
        }

        let mut particular = vec![R::zero(); m];
        for i in 0..rank {
            particular[pivots[i]] = b[i].clone();
        }
        let mut basis = Vec::with_capacity(free.len());
        for (j, &c) in free.iter().enumerate() {
            let mut vector = vec![R::zero(); m];
            vector[c] = R::one();
            for i in 0..rank {
                vector[pivots[i]] = R::neg(&coefficients[i][j]);
            }
            basis.push(vector);
        }
        Some(SystemOfLinearEquationsSolution { particular, basis })
    }

    pub fn inverse(&self) -> Option<Matrix<R>> {
        assert_eq!(self.shape.0, self.shape.1);
        let n = self.shape.0;
        if n >= 64 {
            let m = n / 2;
            let a = Self::new_with((m, m), |i, j| self[i][j].clone());
            if let Some(mut ai) = a.inverse() {
                let b = Self::new_with((m, n - m), |i, j| self[i][j + m].clone());
                let c = Self::new_with((n - m, m), |i, j| self[i + m][j].clone());
                let mut d = Self::new_with((n - m, n - m), |i, j| self[i + m][j + m].clone());
                let u = &ai * &b;
                let v = &c * &ai;
                d -= &v * &b;
                let di = d.inverse()?;
                let r = &u * &di;
                let t = &di * &v;
                ai += &r * &v;
                let mut inverse = Self::zeros((n, n));
                for i in 0..m {
                    inverse[i][..m].clone_from_slice(&ai[i]);
                    for (x, y) in inverse[i][m..].iter_mut().zip(&r[i]) {
                        *x = R::neg(y);
                    }
                }
                for i in m..n {
                    for (x, y) in inverse[i][..m].iter_mut().zip(&t[i - m]) {
                        *x = R::neg(y);
                    }
                    inverse[i][m..].clone_from_slice(&di[i - m]);
                }
                return Some(inverse);
            }
        }
        let mut a = self.clone();
        let mut inverse = Self::eye((n, n));
        let mut ranges: Vec<_> = (0..n).map(|i| (i, i + 1)).collect();
        for r in 0..n {
            let pivot = (r..n).find(|&i| !R::is_zero(&a[i][r]))?;
            a.data.swap(r, pivot);
            inverse.data.swap(r, pivot);
            ranges.swap(r, pivot);

            let d = R::inv(&a[r][r]);
            for x in &mut a[r][r..] {
                R::mul_assign(x, &d);
            }
            let (left, right) = ranges[r];
            for x in &mut inverse[r][left..right] {
                R::mul_assign(x, &d);
            }

            let (a_upper, a_lower) = a.data.split_at_mut(r + 1);
            let pivot_a = &a_upper[r];
            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r + 1);
            let pivot_inverse = &inverse_upper[r];
            let (ranges_upper, ranges_lower) = ranges.split_at_mut(r + 1);
            let (left, right) = ranges_upper[r];
            for ((a, inverse), range) in a_lower.iter_mut().zip(inverse_lower).zip(ranges_lower) {
                if R::is_zero(&a[r]) {
                    continue;
                }
                let e = a[r].clone();
                a[r] = R::zero();
                R::add_scaled_assign(&mut a[(r + 1)..], &pivot_a[(r + 1)..], &R::neg(&e));
                R::add_scaled_assign(
                    &mut inverse[left..right],
                    &pivot_inverse[left..right],
                    &R::neg(&e),
                );
                range.0 = range.0.min(left);
                range.1 = range.1.max(right);
            }
        }
        for r in (0..n).rev() {
            let (left, right) = ranges[r];
            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r);
            let pivot_inverse = &inverse_lower[0];
            let (ranges_upper, _) = ranges.split_at_mut(r);
            for ((a, inverse), range) in a.data[..r].iter_mut().zip(inverse_upper).zip(ranges_upper)
            {
                if R::is_zero(&a[r]) {
                    continue;
                }
                let e = a[r].clone();
                a[r] = R::zero();
                R::add_scaled_assign(
                    &mut inverse[left..right],
                    &pivot_inverse[left..right],
                    &R::neg(&e),
                );
                range.0 = range.0.min(left);
                range.1 = range.1.max(right);
            }
        }
        Some(inverse)
    }

    pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
        let n = self.shape.0;
        if n == 0 {
            return vec![R::one()];
        }
        assert!(self.data.iter().all(|a| a.len() == n));
        for j in 0..(n - 1) {
            if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
                self.data.swap(j + 1, x);
                self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
                let inv = R::inv(&self[j + 1][j]);
                let mut v = vec![];
                let src = std::mem::take(&mut self[j + 1]);
                for a in self.data[(j + 2)..].iter_mut() {
                    let mul = R::mul(&a[j], &inv);
                    R::add_scaled_assign(&mut a[j..], &src[j..], &R::neg(&mul));
                    v.push(mul);
                }
                self[j + 1] = src;
                for a in self.data.iter_mut() {
                    let v = R::dot_product(&a[(j + 2)..], &v);
                    R::add_assign(&mut a[j + 1], &v);
                }
            }
        }
        // dp[k][j - k] stores [x^k] det(xI - A[..j, ..j]).
        let mut dp: Vec<Vec<R::T>> = (0..=n).map(|i| Vec::with_capacity(n + 1 - i)).collect();
        dp[0].push(R::one());
        for i in 0..n {
            let mut c = vec![R::zero(); i + 1];
            c[i] = R::neg(&self[i][i]);
            let mut mul = R::one();
            for j in (0..i).rev() {
                mul = R::mul(&mul, &self[j + 1][j]);
                c[j] = R::neg(&R::mul(&mul, &self[j][i]));
            }
            for k in (0..=i).rev() {
                let mut value = R::dot_product(&dp[k], &c[k..]);
                if k > 0 {
                    R::add_assign(&mut value, dp[k - 1].last().unwrap());
                }
                dp[k].push(value);
            }
            dp[i + 1].push(R::one());
        }
        dp.into_iter().map(|mut c| c.pop().unwrap()).collect()
    }
}

impl<R> Index<usize> for Matrix<R>
where
    R: SemiRing,
{
    type Output = Vec<R::T>;
    fn index(&self, index: usize) -> &Self::Output {
        &self.data[index]
    }
}

impl<R> IndexMut<usize> for Matrix<R>
where
    R: SemiRing,
{
    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
        &mut self.data[index]
    }
}

impl<R> Index<(usize, usize)> for Matrix<R>
where
    R: SemiRing,
{
    type Output = R::T;
    fn index(&self, index: (usize, usize)) -> &Self::Output {
        &self.data[index.0][index.1]
    }
}

impl<R> IndexMut<(usize, usize)> for Matrix<R>
where
    R: SemiRing,
{
    fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
        &mut self.data[index.0][index.1]
    }
}

macro_rules! impl_matrix_pairwise_binop {
    ($imp:ident, $method:ident, $imp_assign:ident, $method_assign:ident $(where [$($clauses:tt)*])?) => {
        impl<R> $imp_assign for Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            fn $method_assign(&mut self, rhs: Self) {
                self.pairwise_assign(&rhs, |a, b| R::$method_assign(a, b));
            }
        }
        impl<R> $imp_assign<&Matrix<R>> for Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            fn $method_assign(&mut self, rhs: &Self) {
                self.pairwise_assign(rhs, |a, b| R::$method_assign(a, b));
            }
        }
        impl<R> $imp for Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            type Output = Matrix<R>;
            fn $method(mut self, rhs: Self) -> Self::Output {
                self.$method_assign(rhs);
                self
            }
        }
        impl<R> $imp<&Matrix<R>> for Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            type Output = Matrix<R>;
            fn $method(mut self, rhs: &Self) -> Self::Output {
                self.$method_assign(rhs);
                self
            }
        }
        impl<R> $imp<Matrix<R>> for &Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            type Output = Matrix<R>;
            fn $method(self, mut rhs: Matrix<R>) -> Self::Output {
                rhs.pairwise_assign(self, |a, b| *a = R::$method(b, a));
                rhs
            }
        }
        impl<R> $imp<&Matrix<R>> for &Matrix<R>
        where
            R: SemiRing,
            $($($clauses)*)?
        {
            type Output = Matrix<R>;
            fn $method(self, rhs: &Matrix<R>) -> Self::Output {
                let mut this = self.clone();
                this.$method_assign(rhs);
                this
            }
        }
    };
}

impl_matrix_pairwise_binop!(Add, add, AddAssign, add_assign);
impl_matrix_pairwise_binop!(Sub, sub, SubAssign, sub_assign where [R: SemiRing<Additive: Invertible>]);

impl<R> Mul for Matrix<R>
where
    R: SemiRing,
{
    type Output = Matrix<R>;
    fn mul(self, rhs: Self) -> Self::Output {
        (&self).mul(&rhs)
    }
}
impl<R> Mul<&Matrix<R>> for Matrix<R>
where
    R: SemiRing,
{
    type Output = Matrix<R>;
    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
        (&self).mul(rhs)
    }
}
impl<R> Mul<Matrix<R>> for &Matrix<R>
where
    R: SemiRing,
{
    type Output = Matrix<R>;
    fn mul(self, rhs: Matrix<R>) -> Self::Output {
        self.mul(&rhs)
    }
}
impl<R> Mul<&Matrix<R>> for &Matrix<R>
where
    R: SemiRing,
{
    type Output = Matrix<R>;
    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
        assert_eq!(self.shape.1, rhs.shape.0);
        if let Some(data) = R::try_matrix_product(&self.data, &rhs.data) {
            return Matrix::from_vec(data);
        }
        let rhs = rhs.transpose();
        Matrix::new_with((self.shape.0, rhs.shape.0), |i, j| {
            R::dot_product(&self[i], &rhs[j])
        })
    }
}

fn strassen_rec<R: Ring>(
    a: &[R::T],
    b: &[R::T],
    c: &mut [R::T],
    shape: (usize, usize, usize),
    stride_a: usize,
    stride_b: usize,
) {
    let (n, m, p) = shape;
    fn add_block<R: Ring>(
        a: &[R::T],
        b: &[R::T],
        out: &mut [R::T],
        n: usize,
        stride_a: usize,
        stride_b: usize,
    ) {
        for ((a, b), c) in a
            .chunks(stride_a)
            .zip(b.chunks(stride_b))
            .zip(out.chunks_exact_mut(n))
        {
            for ((a, b), c) in a.iter().zip(b.iter()).zip(c.iter_mut()) {
                *c = R::add(a, b);
            }
        }
    }

    fn sub_block<R: Ring>(
        a: &[R::T],
        b: &[R::T],
        out: &mut [R::T],
        n: usize,
        stride_a: usize,
        stride_b: usize,
    ) {
        for ((a, b), c) in a
            .chunks(stride_a)
            .zip(b.chunks(stride_b))
            .zip(out.chunks_exact_mut(n))
        {
            for ((a, b), c) in a.iter().zip(b.iter()).zip(c.iter_mut()) {
                *c = R::sub(a, b);
            }
        }
    }

    if n.min(m).min(p) <= 128 {
        let transposed: Vec<_> = (0..p)
            .flat_map(|j| (0..m).map(move |i| b[i * stride_b + j].clone()))
            .collect();
        for (a, c) in a.chunks(stride_a).zip(c.chunks_exact_mut(p)) {
            for (b, c) in transposed.chunks_exact(m).zip(c) {
                *c = R::dot_product(&a[..m], b);
            }
        }
        return;
    }
    let (h, k, w) = (n / 2, m / 2, p / 2);
    let a11 = 0;
    let a12 = k;
    let a21 = h * stride_a;
    let a22 = a21 + k;
    let b11 = 0;
    let b12 = w;
    let b21 = k * stride_b;
    let b22 = b21 + w;

    let block = h * w;
    let mut buf = vec![R::zero(); h * k + k * w + block * 7];
    let (s1, rest) = buf.split_at_mut(h * k);
    let (s2, m_buf) = rest.split_at_mut(k * w);
    let (m1, rest) = m_buf.split_at_mut(block);
    let (m2, rest) = rest.split_at_mut(block);
    let (m3, rest) = rest.split_at_mut(block);
    let (m4, rest) = rest.split_at_mut(block);
    let (m5, rest) = rest.split_at_mut(block);
    let (m6, m7) = rest.split_at_mut(block);

    // (A11 + A22)(B11 + B22)
    add_block::<R>(&a[a11..], &a[a22..], s1, k, stride_a, stride_a);
    add_block::<R>(&b[b11..], &b[b22..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(s1, s2, m1, (h, k, w), k, w);

    // (A21 + A22) B11
    add_block::<R>(&a[a21..], &a[a22..], s1, k, stride_a, stride_a);
    strassen_rec::<R>(s1, &b[b11..], m2, (h, k, w), k, stride_b);

    // A11 (B12 - B22)
    sub_block::<R>(&b[b12..], &b[b22..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(&a[a11..], s2, m3, (h, k, w), stride_a, w);

    // A22 (B21 - B11)
    sub_block::<R>(&b[b21..], &b[b11..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(&a[a22..], s2, m4, (h, k, w), stride_a, w);

    // (A11 + A12) B22
    add_block::<R>(&a[a11..], &a[a12..], s1, k, stride_a, stride_a);
    strassen_rec::<R>(s1, &b[b22..], m5, (h, k, w), k, stride_b);

    // (A21 - A11)(B11 + B12)
    sub_block::<R>(&a[a21..], &a[a11..], s1, k, stride_a, stride_a);
    add_block::<R>(&b[b11..], &b[b12..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(s1, s2, m6, (h, k, w), k, w);

    // (A12 - A22)(B21 + B22)
    sub_block::<R>(&a[a12..], &a[a22..], s1, k, stride_a, stride_a);
    add_block::<R>(&b[b21..], &b[b22..], s2, w, stride_b, stride_b);
    strassen_rec::<R>(s1, s2, m7, (h, k, w), k, w);

    let c11 = 0;
    let c12 = w;
    let c21 = h * p;
    let c22 = c21 + w;
    for ((((m1, m4), m5), m7), c) in m1
        .iter()
        .zip(m4.iter())
        .zip(m5.iter())
        .zip(m7.iter())
        .zip(c[c11..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::add(m1, m4);
        R::sub_assign(c, m5);
        R::add_assign(c, m7);
    }
    for ((m3, m5), c) in m3
        .iter()
        .zip(m5.iter())
        .zip(c[c12..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::add(m3, m5);
    }
    for ((m2, m4), c) in m2
        .iter()
        .zip(m4.iter())
        .zip(c[c21..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::add(m2, m4);
    }
    for ((((m1, m2), m3), m6), c) in m1
        .iter()
        .zip(m2.iter())
        .zip(m3.iter())
        .zip(m6.iter())
        .zip(c[c22..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
    {
        *c = R::sub(m1, m2);
        R::add_assign(c, m3);
        R::add_assign(c, m6);
    }
}

Dyn Compatibility

This trait is not dyn compatible.

In older versions of Rust, dyn compatibility was called "object safety".

Implementors

impl<R> Ring for R

where R: SemiRing<Additive: Invertible>,