Trait SemiRing 

pub trait SemiRing {
    type T: Clone;
    type Additive: AbelianMonoid<T = Self::T>;
    type Multiplicative: Monoid<T = Self::T>;

    // Provided methods
    fn zero() -> Self::T { ... }
    fn is_zero(x: &Self::T) -> bool
       where Self::T: PartialEq { ... }
    fn one() -> Self::T { ... }
    fn is_one(x: &Self::T) -> bool
       where Self::T: PartialEq { ... }
    fn add(x: &Self::T, y: &Self::T) -> Self::T { ... }
    fn mul(x: &Self::T, y: &Self::T) -> Self::T { ... }
    fn try_matrix_product(
        _a: &[Vec<Self::T>],
        _b: &[Vec<Self::T>],
    ) -> Option<Vec<Vec<Self::T>>> { ... }
    fn dot_product(x: &[Self::T], y: &[Self::T]) -> Self::T { ... }
    fn add_scaled_assign(x: &mut [Self::T], y: &[Self::T], a: &Self::T) { ... }
    fn add_assign(x: &mut Self::T, y: &Self::T) { ... }
    fn mul_assign(x: &mut Self::T, y: &Self::T) { ... }
}

Required Associated Types

type T: Clone

type Additive: AbelianMonoid<T = Self::T>

type Multiplicative: Monoid<T = Self::T>

Provided Methods

fn zero() -> Self::T

additive identity: $0$

Examples found in repository

crates/competitive/src/graph/shortest_path.rs (line 100)

    fn inf() -> Self::T {
        PartialIgnoredOrd(M::T::maximum(), S::zero())
    }

crates/competitive/src/algebra/ring.rs (line 20)

    fn is_zero(x: &Self::T) -> bool
    where
        Self::T: PartialEq,
    {
        *x == Self::zero()
    }
    /// multiplicative identity: $1$
    fn one() -> Self::T {
        <Self::Multiplicative as Unital>::unit()
    }
    fn is_one(x: &Self::T) -> bool
    where
        Self::T: PartialEq,
    {
        *x == Self::one()
    }
    /// additive operaion: $+$
    fn add(x: &Self::T, y: &Self::T) -> Self::T {
        <Self::Additive as Magma>::operate(x, y)
    }
    /// multiplicative operaion: $+$
    fn mul(x: &Self::T, y: &Self::T) -> Self::T {
        <Self::Multiplicative as Magma>::operate(x, y)
    }

    fn try_matrix_product(_a: &[Vec<Self::T>], _b: &[Vec<Self::T>]) -> Option<Vec<Vec<Self::T>>> {
        None
    }

    fn dot_product(x: &[Self::T], y: &[Self::T]) -> Self::T {
        assert_eq!(x.len(), y.len());
        x.iter().zip(y).fold(Self::zero(), |mut sum, (x, y)| {
            Self::add_assign(&mut sum, &Self::mul(x, y));
            sum
        })
    }

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

    pub fn zeros(shape: (usize, usize)) -> Self {
        Self {
            shape,
            data: vec![vec![R::zero(); shape.1]; shape.0],
            _marker: PhantomData,
        }
    }

    pub fn eye(shape: (usize, usize)) -> Self {
        let mut data = vec![vec![R::zero(); shape.1]; shape.0];
        for (i, d) in data.iter_mut().enumerate().take(shape.1) {
            d[i] = R::one();
        }
        Self {
            shape,
            data,
            _marker: PhantomData,
        }
    }

    pub fn transpose(&self) -> Self {
        Self::new_with((self.shape.1, self.shape.0), |i, j| self[j][i].clone())
    }

    pub fn map<S, F>(&self, mut f: F) -> Matrix<S>
    where
        S: SemiRing,
        F: FnMut(&R::T) -> S::T,
    {
        Matrix::<S>::new_with(self.shape, |i, j| f(&self[i][j]))
    }

    pub fn add_row_with(&mut self, mut f: impl FnMut(usize, usize) -> R::T) {
        self.data
            .push((0..self.shape.1).map(|j| f(self.shape.0, j)).collect());
        self.shape.0 += 1;
    }

    pub fn add_col_with(&mut self, mut f: impl FnMut(usize, usize) -> R::T) {
        for i in 0..self.shape.0 {
            self.data[i].push(f(i, self.shape.1));
        }
        self.shape.1 += 1;
    }

    pub fn pairwise_assign<F>(&mut self, other: &Self, mut f: F)
    where
        F: FnMut(&mut R::T, &R::T),
    {
        assert_eq!(self.shape, other.shape);
        for i in 0..self.shape.0 {
            for j in 0..self.shape.1 {
                f(&mut self[i][j], &other[i][j]);
            }
        }
    }
}

#[derive(Debug)]
pub struct SystemOfLinearEquationsSolution<R>
where
    R: Field<Additive: Invertible, Multiplicative: Invertible>,
{
    pub particular: Vec<R::T>,
    pub basis: Vec<Vec<R::T>>,
}

impl<R> Matrix<R>
where
    R: Field<T: PartialEq, Additive: Invertible, Multiplicative: Invertible>,
{
    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
    }

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

    fn to_x() -> FloorSumData<R, X, Y> {
        let mut dp = array![array![R::zero(); Y]; X];
        dp[0][0] = R::one();
        FloorSumData {
            dp,
            dx: R::one(),
            dy: R::zero(),
            _marker: PhantomData,
        }
    }
    fn to_y() -> FloorSumData<R, X, Y> {
        FloorSumData {
            dp: array![array![R::zero(); Y]; X],
            dx: R::zero(),
            dy: R::one(),
            _marker: PhantomData,
        }
    }
}

impl<R, const X: usize, const Y: usize> FloorSum<R, X, Y>
where
    R: Ring<Additive: Invertible>,
{
    fn offset(x: i64, y: i64) -> FloorSumData<R, X, Y> {
        FloorSumData {
            dp: array![array![R::zero(); Y]; X],
            dx: R::Additive::signed_pow(R::one(), x),
            dy: R::Additive::signed_pow(R::one(), y),
            _marker: PhantomData,
        }
    }
}

impl<R, const X: usize, const Y: usize> Magma for FloorSum<R, X, Y>
where
    R: SemiRing,
{
    type T = FloorSumData<R, X, Y>;

    fn operate(a: &Self::T, b: &Self::T) -> Self::T {
        let mut a = a.clone();
        let mut b = b.clone();
        let mut pow_x = array![R::zero(); X];
        let mut pow_y = array![R::zero(); Y];
        pow_x[0] = R::one();
        pow_y[0] = R::one();
        for i in 1..X {
            pow_x[i] = R::mul(&pow_x[i - 1], &a.dx);
        }
        for j in 1..Y {
            pow_y[j] = R::mul(&pow_y[j - 1], &a.dy);
        }
        macro_rules! go {
            ($N:ident) => {
                let mut comb = array![array![R::zero(); $N]; $N];
                comb[0][0] = R::one();
                let mut i = 0;
                while i + 1 < $N {
                    let mut j = 0;
                    while j <= i {
                        let x = comb[i][j].clone();
                        R::add_assign(&mut comb[i + 1][j], &x);
                        R::add_assign(&mut comb[i + 1][j + 1], &x);
                        j += 1;
                    }
                    i += 1;
                }
                for i in 0..X {
                    for j in (0..Y).rev() {
                        for k in j + 1..Y {
                            let mut x = b.dp[i][j].clone();
                            R::mul_assign(&mut x, &comb[k][j]);
                            R::mul_assign(&mut x, &pow_y[k - j]);
                            R::add_assign(&mut b.dp[i][k], &x);
                        }
                    }
                }
                for j in 0..Y {
                    for i in (0..X).rev() {
                        for k in i..X {
                            let mut x = b.dp[i][j].clone();
                            R::mul_assign(&mut x, &comb[k][i]);
                            R::mul_assign(&mut x, &pow_x[k - i]);
                            R::add_assign(&mut a.dp[k][j], &x);
                        }
                    }
                }
            };
        }
        if X <= Y {
            go!(Y);
        } else {
            go!(X);
        }
        R::add_assign(&mut a.dx, &b.dx);
        R::add_assign(&mut a.dy, &b.dy);
        a
    }
}

impl<R, const X: usize, const Y: usize> Unital for FloorSum<R, X, Y>
where
    R: SemiRing,
{
    fn unit() -> Self::T {
        FloorSumData {
            dp: array![array![R::zero(); Y]; X],
            dx: R::zero(),
            dy: R::zero(),
            _marker: PhantomData,
        }
    }
}

impl<R, const X: usize, const Y: usize> Associative for FloorSum<R, X, Y> where R: SemiRing {}

fn floor_monoid_product<M>(
    mut x: M::T,
    mut y: M::T,
    mut n: u64,
    mut a: u64,
    mut b: u64,
    mut m: u64,
) -> M::T
where
    M: Monoid,
{
    let mut c = (a * n + b) / m;
    let mut pre = M::unit();
    let mut suf = M::unit();
    loop {
        let (p, q) = (a / m, b / m);
        a %= m;
        b %= m;
        x = M::operate(&x, &M::pow(y.clone(), p));
        pre = M::operate(&pre, &M::pow(y.clone(), q));
        c -= p * n + q;
        if c == 0 {
            break;
        }
        let d = (m * c - b - 1) / a + 1;
        suf = M::operate(&y, &M::operate(&M::pow(x.clone(), n - d), &suf));
        b = m - b - 1 + a;
        n = c - 1;
        c = d;
        swap(&mut m, &mut a);
        swap(&mut x, &mut y);
    }
    x = M::pow(x.clone(), n);
    M::operate(&M::operate(&pre, &x), &suf)
}

/// $$\sum_{i=0}^{n-1}i^X\left\lfloor\frac{a\times i+b}{m}\right\rfloor^Y$$
pub fn floor_sum_polynomial<T, const X: usize, const Y: usize>(
    n: u64,
    a: u64,
    b: u64,
    m: u64,
) -> [[T; Y]; X]
where
    T: DotProduct + One,
{
    debug_assert!(a == 0 || n < (u64::MAX - b) / a);
    floor_monoid_product::<FloorSum<AddMulOperation<T>, X, Y>>(
        FloorSum::<AddMulOperation<T>, X, Y>::to_x(),
        FloorSum::<AddMulOperation<T>, X, Y>::to_y(),
        n,
        a,
        b,
        m,
    )
    .dp
}

/// $$\sum_{i=l}^{r-1}i^X\left\lfloor\frac{a\times i+b}{m}\right\rfloor^Y$$
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
}

#[derive(Debug)]
struct FloorPowerSum<R>
where
    R: SemiRing,
{
    x: R::T,
    sum: R::T,
}

impl<R> Clone for FloorPowerSum<R>
where
    R: SemiRing,
{
    fn clone(&self) -> Self {
        Self {
            x: self.x.clone(),
            sum: self.sum.clone(),
        }
    }
}

impl<R> FloorPowerSum<R>
where
    R: SemiRing,
{
    fn to_x(x: R::T) -> Self {
        Self { x, sum: R::one() }
    }
    fn to_y(y: R::T) -> Self {
        Self {
            x: y,
            sum: R::zero(),
        }
    }
}

impl<R> Magma for FloorPowerSum<R>
where
    R: SemiRing,
{
    type T = Self;

    fn operate(a: &Self::T, b: &Self::T) -> Self::T {
        Self {
            x: R::mul(&a.x, &b.x),
            sum: R::add(&a.sum, &R::mul(&a.x, &b.sum)),
        }
    }
}

impl<R> Unital for FloorPowerSum<R>
where
    R: SemiRing,
{
    fn unit() -> Self::T {
        Self {
            x: R::one(),
            sum: R::zero(),
        }
    }

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

    fn ranked(t: Vec<R::T>, len: usize) -> (Vec<R::T>, usize) {
        let width = len.trailing_zeros() as usize + 1;
        let mut ranked = Vec::with_capacity(len * width);
        advise_huge_pages(&mut ranked);
        ranked.resize(len * width, R::zero());
        for (i, value) in t.into_iter().enumerate() {
            ranked[i * width + i.count_ones() as usize] = value;
        }
        (ranked, width)
    }

    fn diagonal(ranked: Vec<R::T>, width: usize) -> Vec<R::T> {
        ranked
            .chunks_exact(width)
            .enumerate()
            .map(|(i, row)| row[i.count_ones() as usize].clone())
            .collect()
    }

    #[inline]
    fn multiply_row(
        x: &[R::T],
        y: &[R::T],
        right: &mut [R::T],
        output: &mut [R::T],
        rank: usize,
    ) -> usize {
        for (right, y) in right[..=rank].iter_mut().zip(y[..=rank].iter().rev()) {
            right.clone_from(y);
        }
        let end = (rank * 2).min(x.len() - 1);
        for (degree, output) in output.iter_mut().enumerate().take(end + 1).skip(rank) {
            let first = degree - rank;
            *output = R::dot_product(&x[first..=rank], &right[..=rank - first]);
        }
        end
    }
}

impl<R> ConvolveSteps for SubsetConvolve<R>
where
    R: Ring<T: PartialEq, Additive: Invertible>,
{
    type T = Vec<R::T>;
    type F = (Vec<R::T>, usize);

    fn length(t: &Self::T) -> usize {
        t.len()
    }

    fn transform(t: Self::T, len: usize) -> Self::F {
        let (mut f, width) = Self::ranked(t, len);
        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 ranks = lower.count_ones() as usize + 1;
                    let (lower_rows, upper_rows) = f.split_at_mut(upper * width);
                    let lower_row = &lower_rows[lower * width..lower * width + ranks];
                    let upper_row = &mut upper_rows[..ranks];
                    for (upper, lower) in upper_row.iter_mut().zip(lower_row) {
                        R::add_assign(upper, lower);
                    }
                }
            }
        }
        (f, width)
    }

    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/algorithm/automata_learning.rs (line 222)

    fn behavior<I>(&self, input: I) -> Self::Output
    where
        I: IntoIterator<Item = usize>,
    {
        let mut weights = self.initial_weights.clone();
        for x in input {
            weights = &weights * &self.transitions[x];
        }
        let result = &weights * &self.final_weights;
        if result.shape != (0, 0) {
            result[0][0].clone()
        } else {
            F::zero()
        }
    }
}

impl<F> SerdeByteStr for WeightedFiniteAutomaton<F>
where
    F: Field<T: SerdeByteStr, Additive: Invertible, Multiplicative: Invertible>,
{
    fn serialize(&self, buf: &mut Vec<u8>) {
        self.initial_weights.serialize(buf);
        self.transitions.serialize(buf);
        self.final_weights.serialize(buf);
    }

    fn deserialize<I>(iter: &mut I) -> Self
    where
        I: Iterator<Item = u8>,
    {
        let initial_weights = Matrix::deserialize(iter);
        let transitions = Vec::deserialize(iter);
        let final_weights = Matrix::deserialize(iter);
        Self {
            initial_weights,
            transitions,
            final_weights,
        }
    }
}

pub fn dense_sampling(sigma: usize, max_len: usize) -> impl Iterator<Item = Vec<usize>> {
    assert_ne!(sigma, 0, "Sigma must be greater than 0");
    let mut current = vec![];
    once_with(Vec::new).chain(from_fn(move || {
        let mut carry = true;
        for i in (0..current.len()).rev() {
            current[i] += 1;
            if current[i] == sigma {
                current[i] = 0;
            } else {
                carry = false;
                break;
            }
        }
        if carry {
            current.push(0);
        }
        if current.len() <= max_len {
            Some(current.to_vec())
        } else {
            None
        }
    }))
}

pub fn random_sampling(
    sigma: usize,
    len_spec: impl RandomSpec<usize>,
    seconds: f64,
) -> impl Iterator<Item = Vec<usize>> {
    assert_ne!(sigma, 0, "Sigma must be greater than 0");
    let now = Instant::now();
    let mut rng = Xorshift::new();
    from_fn(move || {
        if now.elapsed().as_secs_f64() > seconds {
            None
        } else {
            let n = rng.random(&len_spec);
            Some(rng.random_iter(0..sigma).take(n).collect())
        }
    })
}

#[derive(Debug, Clone)]
pub struct DfaLearning<A>
where
    A: BlackBoxAutomaton<Output = bool>,
{
    automaton: A,
    prefixes: Vec<Vec<usize>>,
    suffixes: Vec<Vec<usize>>,
    table: Vec<BitSet>,
    row_map: HashMap<BitSet, usize>,
}

impl<A> DfaLearning<A>
where
    A: BlackBoxAutomaton<Output = bool>,
{
    pub fn new(automaton: A) -> Self {
        let mut this = Self {
            automaton,
            prefixes: vec![],
            suffixes: vec![],
            table: vec![],
            row_map: HashMap::new(),
        };
        this.add_suffix(vec![]);
        this.add_prefix(vec![]);
        this
    }
    fn add_prefix(&mut self, prefix: Vec<usize>) -> usize {
        let row: BitSet = self
            .suffixes
            .iter()
            .map(|s| {
                self.automaton
                    .behavior(prefix.iter().cloned().chain(s.iter().cloned()))
            })
            .collect();
        *self.row_map.entry(row.clone()).or_insert_with(|| {
            let idx = self.table.len();
            self.table.push(row);
            self.prefixes.push(prefix);
            idx
        })
    }
    fn add_suffix(&mut self, suffix: Vec<usize>) {
        if self.suffixes.contains(&suffix) {
            return;
        }
        for (prefix, table) in self.prefixes.iter_mut().zip(&mut self.table) {
            table.push(
                self.automaton
                    .behavior(prefix.iter().cloned().chain(suffix.iter().cloned())),
            );
        }
        self.suffixes.push(suffix);
        self.row_map.clear();
        for (i_prefix, row) in self.table.iter().enumerate() {
            self.row_map.insert(row.clone(), i_prefix);
        }
    }
    pub fn construct_dfa(&mut self) -> DeterministicFiniteAutomaton {
        let sigma = self.automaton.sigma();
        let mut dfa = DeterministicFiniteAutomaton {
            states: vec![],
            initial_state: 0,
        };
        let mut i_prefix = 0;
        while i_prefix < self.prefixes.len() {
            let mut delta = vec![];
            for x in 0..sigma {
                let prefix: Vec<usize> =
                    self.prefixes[i_prefix].iter().cloned().chain([x]).collect();
                let index = self.add_prefix(prefix);
                delta.push(index);
            }
            dfa.states.push(DfaState {
                delta,
                accept: self.table[i_prefix].get(0),
            });
            i_prefix += 1;
        }
        dfa
    }
    pub fn train_sample(&mut self, dfa: &DeterministicFiniteAutomaton, sample: &[usize]) -> bool {
        let expected = self.automaton.behavior(sample.iter().cloned());
        if expected == dfa.behavior(sample.iter().cloned()) {
            return false;
        }
        let n = sample.len();
        let mut states: Vec<(usize, usize)> = Vec::with_capacity(n + 1);
        let mut s = 0usize;
        states.push((s, 0));
        for (k, &x) in sample.iter().enumerate() {
            s = dfa.states[s].delta[x];
            states.push((s, k + 1));
        }
        let split = states.partition_point(|&(state, k)| {
            self.automaton.behavior(
                self.prefixes[state]
                    .iter()
                    .cloned()
                    .chain(sample[k..].iter().cloned()),
            ) == expected
        });
        let new_prefix = sample[..split].to_vec();
        let new_suffix = sample[split..].to_vec();
        self.add_suffix(new_suffix);
        self.add_prefix(new_prefix);
        true
    }
    pub fn train(
        &mut self,
        samples: impl IntoIterator<Item = Vec<usize>>,
    ) -> DeterministicFiniteAutomaton {
        let mut dfa = self.construct_dfa();
        for sample in samples {
            if self.train_sample(&dfa, &sample) {
                dfa = self.construct_dfa();
            }
        }
        dfa
    }
}

pub struct WfaLearning<F, A>
where
    F: Field<Additive: Invertible, Multiplicative: Invertible>,
    A: BlackBoxAutomaton<Output = F::T>,
{
    automaton: A,
    prefixes: Vec<Vec<usize>>,
    suffixes: Vec<Vec<usize>>,
    inv_h: Matrix<F>,
    nh: Vec<Matrix<F>>,
    wfa: WeightedFiniteAutomaton<F>,
    _marker: PhantomData<fn() -> F>,
}

impl<F, A> Debug for WfaLearning<F, A>
where
    F: Field<T: Debug, Additive: Invertible, Multiplicative: Invertible>,
    A: BlackBoxAutomaton<Output = F::T> + Debug,
{
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.debug_struct("WfaLearning")
            .field("automaton", &self.automaton)
            .field("prefixes", &self.prefixes)
            .field("suffixes", &self.suffixes)
            .field("inv_h", &self.inv_h)
            .field("nh", &self.nh)
            .field("wfa", &self.wfa)
            .finish()
    }
}

impl<F, A> Clone for WfaLearning<F, A>
where
    F: Field<Additive: Invertible, Multiplicative: Invertible>,
    A: BlackBoxAutomaton<Output = F::T> + Clone,
{
    fn clone(&self) -> Self {
        Self {
            automaton: self.automaton.clone(),
            prefixes: self.prefixes.clone(),
            suffixes: self.suffixes.clone(),
            inv_h: self.inv_h.clone(),
            nh: self.nh.clone(),
            wfa: self.wfa.clone(),
            _marker: self._marker,
        }
    }
}

impl<F, A> WfaLearning<F, A>
where
    F: Field<T: PartialEq, Additive: Invertible, Multiplicative: Invertible>,
    A: BlackBoxAutomaton<Output = F::T>,
{
    pub fn new(automaton: A) -> Self {
        let sigma = automaton.sigma();
        Self {
            automaton,
            prefixes: vec![],
            suffixes: vec![],
            inv_h: Matrix::zeros((0, 0)),
            nh: vec![Matrix::zeros((0, 0)); sigma],
            wfa: WeightedFiniteAutomaton {
                initial_weights: Matrix::zeros((1, 0)),
                transitions: vec![Matrix::zeros((0, 0)); sigma],
                final_weights: Matrix::zeros((0, 1)),
            },
            _marker: PhantomData,
        }
    }
    pub fn wfa(&self) -> &WeightedFiniteAutomaton<F> {
        &self.wfa
    }
    fn split_sample(&mut self, sample: &[usize]) -> Option<(Vec<usize>, Vec<usize>)> {
        if self.prefixes.is_empty() && !F::is_zero(&self.automaton.behavior(sample.iter().cloned()))
        {
            return Some((vec![], sample.to_vec()));
        }
        let expected = self.automaton.behavior(sample.iter().cloned());
        if expected == self.wfa.behavior(sample.iter().cloned()) {
            return None;
        }
        let n = sample.len();
        let dim = self.wfa.final_weights.shape.0;
        let mut states: Vec<(Matrix<F>, usize)> = Vec::with_capacity(n + 1);
        let mut v = self.wfa.final_weights.clone();
        states.push((v.clone(), n));
        for k in (0..n).rev() {
            v = &self.wfa.transitions[sample[k]] * &v;
            states.push((v.clone(), k));
        }
        states.reverse();
        let split = states.partition_point(|(state, k)| {
            (0..dim).any(|j| {
                self.automaton.behavior(
                    self.prefixes[j]
                        .iter()
                        .cloned()
                        .chain(sample[*k..].iter().cloned()),
                ) != state[j][0]
            })
        });
        Some((sample[..split].to_vec(), sample[split..].to_vec()))
    }
    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
    }
    pub fn train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
        for sample in samples {
            self.train_sample(&sample);
        }
    }
    pub fn batch_train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
        let mut prefix_set: HashSet<_> = self.prefixes.iter().cloned().collect();
        let mut suffix_set: HashSet<_> = self.suffixes.iter().cloned().collect();
        for sample in samples {
            if prefix_set.insert(sample.to_vec()) {
                self.prefixes.push(sample.to_vec());
            }
            if suffix_set.insert(sample.to_vec()) {
                self.suffixes.push(sample);
            }
        }
        let mut h = Matrix::<F>::new_with((self.prefixes.len(), self.suffixes.len()), |i, j| {
            self.automaton.behavior(
                self.prefixes[i]
                    .iter()
                    .cloned()
                    .chain(self.suffixes[j].iter().cloned()),
            )
        });
        if !self.prefixes.is_empty() && !self.suffixes.is_empty() && F::is_zero(&h[0][0]) {
            for j in 1..self.suffixes.len() {
                if !F::is_zero(&h[0][j]) {
                    self.suffixes.swap(0, j);
                    for row in &mut h.data {
                        row.swap(0, j);
                    }
                    break;
                }
            }
        }
        let mut row_id: Vec<usize> = (0..h.shape.0).collect();
        let mut pivots = vec![];
        h.row_reduction_with(false, |r, p, c| {
            row_id.swap(r, p);
            pivots.push((row_id[r], c));
        });
        let mut new_prefixes = vec![];
        let mut new_suffixes = vec![];
        for (i, j) in pivots {
            new_prefixes.push(self.prefixes[i].clone());
            new_suffixes.push(self.suffixes[j].clone());
        }
        self.prefixes = new_prefixes;
        self.suffixes = new_suffixes;
        assert_eq!(self.prefixes.len(), self.suffixes.len());
        let n = self.prefixes.len();
        let h = Matrix::<F>::new_with((n, n), |i, j| {
            self.automaton.behavior(
                self.prefixes[i]
                    .iter()
                    .cloned()
                    .chain(self.suffixes[j].iter().cloned()),
            )
        });
        self.inv_h = h.inverse().expect("Hankel matrix must be invertible");
        self.wfa = WeightedFiniteAutomaton::<F> {
            initial_weights: Matrix::new_with((1, n), |_, j| {
                if self.prefixes[j].is_empty() {
                    F::one()
                } else {
                    F::zero()
                }
            }),
            transitions: (0..self.automaton.sigma())
                .map(|x| {
                    &Matrix::new_with((n, n), |i, j| {
                        self.automaton.behavior(
                            self.prefixes[i]
                                .iter()
                                .cloned()
                                .chain([x])
                                .chain(self.suffixes[j].iter().cloned()),
                        )
                    }) * &self.inv_h
                })
                .collect(),
            final_weights: Matrix::new_with((n, 1), |i, _| {
                self.automaton.behavior(self.prefixes[i].iter().cloned())
            }),
        };
    }

Additional examples can be found in:

• crates/competitive/src/algorithm/esper.rs
• crates/competitive/src/math/black_box_matrix.rs

fn is_zero(x: &Self::T) -> bool

where Self::T: PartialEq,

Examples found in repository

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

    pub fn new_with<F>(shape: (usize, usize), f: F) -> Self
    where
        R: SemiRing<T: PartialEq>,
        F: Fn(usize, usize) -> R::T,
    {
        let mut nonzero = vec![];
        for i in 0..shape.0 {
            for j in 0..shape.1 {
                let v = f(i, j);
                if !R::is_zero(&v) {
                    nonzero.push((i, j, v));
                }
            }
        }
        Self { shape, nonzero }
    }
    pub fn from_nonzero(shape: (usize, usize), nonzero: Vec<(usize, usize, R::T)>) -> Self {
        Self { shape, nonzero }
    }
}

impl<R> SparseMatrix<R>
where
    R: Field<T: PartialEq, Additive: Invertible, Multiplicative: Invertible>,
{
    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
    }
}

impl<R> From<Matrix<R>> for SparseMatrix<R>
where
    R: SemiRing<T: PartialEq>,
{
    fn from(mat: Matrix<R>) -> Self {
        let mut nonzero = vec![];
        for i in 0..mat.shape.0 {
            for j in 0..mat.shape.1 {
                let v = mat[(i, j)].clone();
                if !R::is_zero(&v) {
                    nonzero.push((i, j, v));
                }
            }
        }
        Self {
            shape: mat.shape,
            nonzero,
        }
    }

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

    fn split_sample(&mut self, sample: &[usize]) -> Option<(Vec<usize>, Vec<usize>)> {
        if self.prefixes.is_empty() && !F::is_zero(&self.automaton.behavior(sample.iter().cloned()))
        {
            return Some((vec![], sample.to_vec()));
        }
        let expected = self.automaton.behavior(sample.iter().cloned());
        if expected == self.wfa.behavior(sample.iter().cloned()) {
            return None;
        }
        let n = sample.len();
        let dim = self.wfa.final_weights.shape.0;
        let mut states: Vec<(Matrix<F>, usize)> = Vec::with_capacity(n + 1);
        let mut v = self.wfa.final_weights.clone();
        states.push((v.clone(), n));
        for k in (0..n).rev() {
            v = &self.wfa.transitions[sample[k]] * &v;
            states.push((v.clone(), k));
        }
        states.reverse();
        let split = states.partition_point(|(state, k)| {
            (0..dim).any(|j| {
                self.automaton.behavior(
                    self.prefixes[j]
                        .iter()
                        .cloned()
                        .chain(sample[*k..].iter().cloned()),
                ) != state[j][0]
            })
        });
        Some((sample[..split].to_vec(), sample[split..].to_vec()))
    }
    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
    }
    pub fn train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
        for sample in samples {
            self.train_sample(&sample);
        }
    }
    pub fn batch_train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
        let mut prefix_set: HashSet<_> = self.prefixes.iter().cloned().collect();
        let mut suffix_set: HashSet<_> = self.suffixes.iter().cloned().collect();
        for sample in samples {
            if prefix_set.insert(sample.to_vec()) {
                self.prefixes.push(sample.to_vec());
            }
            if suffix_set.insert(sample.to_vec()) {
                self.suffixes.push(sample);
            }
        }
        let mut h = Matrix::<F>::new_with((self.prefixes.len(), self.suffixes.len()), |i, j| {
            self.automaton.behavior(
                self.prefixes[i]
                    .iter()
                    .cloned()
                    .chain(self.suffixes[j].iter().cloned()),
            )
        });
        if !self.prefixes.is_empty() && !self.suffixes.is_empty() && F::is_zero(&h[0][0]) {
            for j in 1..self.suffixes.len() {
                if !F::is_zero(&h[0][j]) {
                    self.suffixes.swap(0, j);
                    for row in &mut h.data {
                        row.swap(0, j);
                    }
                    break;
                }
            }
        }
        let mut row_id: Vec<usize> = (0..h.shape.0).collect();
        let mut pivots = vec![];
        h.row_reduction_with(false, |r, p, c| {
            row_id.swap(r, p);
            pivots.push((row_id[r], c));
        });
        let mut new_prefixes = vec![];
        let mut new_suffixes = vec![];
        for (i, j) in pivots {
            new_prefixes.push(self.prefixes[i].clone());
            new_suffixes.push(self.suffixes[j].clone());
        }
        self.prefixes = new_prefixes;
        self.suffixes = new_suffixes;
        assert_eq!(self.prefixes.len(), self.suffixes.len());
        let n = self.prefixes.len();
        let h = Matrix::<F>::new_with((n, n), |i, j| {
            self.automaton.behavior(
                self.prefixes[i]
                    .iter()
                    .cloned()
                    .chain(self.suffixes[j].iter().cloned()),
            )
        });
        self.inv_h = h.inverse().expect("Hankel matrix must be invertible");
        self.wfa = WeightedFiniteAutomaton::<F> {
            initial_weights: Matrix::new_with((1, n), |_, j| {
                if self.prefixes[j].is_empty() {
                    F::one()
                } else {
                    F::zero()
                }
            }),
            transitions: (0..self.automaton.sigma())
                .map(|x| {
                    &Matrix::new_with((n, n), |i, j| {
                        self.automaton.behavior(
                            self.prefixes[i]
                                .iter()
                                .cloned()
                                .chain([x])
                                .chain(self.suffixes[j].iter().cloned()),
                        )
                    }) * &self.inv_h
                })
                .collect(),
            final_weights: Matrix::new_with((n, 1), |i, _| {
                self.automaton.behavior(self.prefixes[i].iter().cloned())
            }),
        };
    }

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

    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()
    }

fn one() -> Self::T

multiplicative identity: $1$

Examples found in repository

crates/competitive/src/graph/shortest_path.rs (line 96)

    fn source() -> Self::T {
        PartialIgnoredOrd(M::unit(), S::one())
    }

crates/competitive/src/algebra/ring.rs (line 30)

    fn is_one(x: &Self::T) -> bool
    where
        Self::T: PartialEq,
    {
        *x == Self::one()
    }

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

    fn new(base: R::T) -> Self {
        Self {
            base,
            pow: vec![R::one()],
        }
    }
    fn ensure_pow(&mut self, len: usize) {
        if self.pow.len() <= len {
            self.pow.reserve(len - self.pow.len() + 1);
            if self.pow.is_empty() {
                self.pow.push(R::one());
            }
            for _ in 0..=len - self.pow.len() {
                self.pow.push(R::mul(self.pow.last().unwrap(), &self.base));
            }
        }
    }

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

    fn to_x() -> FloorSumData<R, X, Y> {
        let mut dp = array![array![R::zero(); Y]; X];
        dp[0][0] = R::one();
        FloorSumData {
            dp,
            dx: R::one(),
            dy: R::zero(),
            _marker: PhantomData,
        }
    }
    fn to_y() -> FloorSumData<R, X, Y> {
        FloorSumData {
            dp: array![array![R::zero(); Y]; X],
            dx: R::zero(),
            dy: R::one(),
            _marker: PhantomData,
        }
    }
}

impl<R, const X: usize, const Y: usize> FloorSum<R, X, Y>
where
    R: Ring<Additive: Invertible>,
{
    fn offset(x: i64, y: i64) -> FloorSumData<R, X, Y> {
        FloorSumData {
            dp: array![array![R::zero(); Y]; X],
            dx: R::Additive::signed_pow(R::one(), x),
            dy: R::Additive::signed_pow(R::one(), y),
            _marker: PhantomData,
        }
    }
}

impl<R, const X: usize, const Y: usize> Magma for FloorSum<R, X, Y>
where
    R: SemiRing,
{
    type T = FloorSumData<R, X, Y>;

    fn operate(a: &Self::T, b: &Self::T) -> Self::T {
        let mut a = a.clone();
        let mut b = b.clone();
        let mut pow_x = array![R::zero(); X];
        let mut pow_y = array![R::zero(); Y];
        pow_x[0] = R::one();
        pow_y[0] = R::one();
        for i in 1..X {
            pow_x[i] = R::mul(&pow_x[i - 1], &a.dx);
        }
        for j in 1..Y {
            pow_y[j] = R::mul(&pow_y[j - 1], &a.dy);
        }
        macro_rules! go {
            ($N:ident) => {
                let mut comb = array![array![R::zero(); $N]; $N];
                comb[0][0] = R::one();
                let mut i = 0;
                while i + 1 < $N {
                    let mut j = 0;
                    while j <= i {
                        let x = comb[i][j].clone();
                        R::add_assign(&mut comb[i + 1][j], &x);
                        R::add_assign(&mut comb[i + 1][j + 1], &x);
                        j += 1;
                    }
                    i += 1;
                }
                for i in 0..X {
                    for j in (0..Y).rev() {
                        for k in j + 1..Y {
                            let mut x = b.dp[i][j].clone();
                            R::mul_assign(&mut x, &comb[k][j]);
                            R::mul_assign(&mut x, &pow_y[k - j]);
                            R::add_assign(&mut b.dp[i][k], &x);
                        }
                    }
                }
                for j in 0..Y {
                    for i in (0..X).rev() {
                        for k in i..X {
                            let mut x = b.dp[i][j].clone();
                            R::mul_assign(&mut x, &comb[k][i]);
                            R::mul_assign(&mut x, &pow_x[k - i]);
                            R::add_assign(&mut a.dp[k][j], &x);
                        }
                    }
                }
            };
        }
        if X <= Y {
            go!(Y);
        } else {
            go!(X);
        }
        R::add_assign(&mut a.dx, &b.dx);
        R::add_assign(&mut a.dy, &b.dy);
        a
    }
}

impl<R, const X: usize, const Y: usize> Unital for FloorSum<R, X, Y>
where
    R: SemiRing,
{
    fn unit() -> Self::T {
        FloorSumData {
            dp: array![array![R::zero(); Y]; X],
            dx: R::zero(),
            dy: R::zero(),
            _marker: PhantomData,
        }
    }
}

impl<R, const X: usize, const Y: usize> Associative for FloorSum<R, X, Y> where R: SemiRing {}

fn floor_monoid_product<M>(
    mut x: M::T,
    mut y: M::T,
    mut n: u64,
    mut a: u64,
    mut b: u64,
    mut m: u64,
) -> M::T
where
    M: Monoid,
{
    let mut c = (a * n + b) / m;
    let mut pre = M::unit();
    let mut suf = M::unit();
    loop {
        let (p, q) = (a / m, b / m);
        a %= m;
        b %= m;
        x = M::operate(&x, &M::pow(y.clone(), p));
        pre = M::operate(&pre, &M::pow(y.clone(), q));
        c -= p * n + q;
        if c == 0 {
            break;
        }
        let d = (m * c - b - 1) / a + 1;
        suf = M::operate(&y, &M::operate(&M::pow(x.clone(), n - d), &suf));
        b = m - b - 1 + a;
        n = c - 1;
        c = d;
        swap(&mut m, &mut a);
        swap(&mut x, &mut y);
    }
    x = M::pow(x.clone(), n);
    M::operate(&M::operate(&pre, &x), &suf)
}

/// $$\sum_{i=0}^{n-1}i^X\left\lfloor\frac{a\times i+b}{m}\right\rfloor^Y$$
pub fn floor_sum_polynomial<T, const X: usize, const Y: usize>(
    n: u64,
    a: u64,
    b: u64,
    m: u64,
) -> [[T; Y]; X]
where
    T: DotProduct + One,
{
    debug_assert!(a == 0 || n < (u64::MAX - b) / a);
    floor_monoid_product::<FloorSum<AddMulOperation<T>, X, Y>>(
        FloorSum::<AddMulOperation<T>, X, Y>::to_x(),
        FloorSum::<AddMulOperation<T>, X, Y>::to_y(),
        n,
        a,
        b,
        m,
    )
    .dp
}

/// $$\sum_{i=l}^{r-1}i^X\left\lfloor\frac{a\times i+b}{m}\right\rfloor^Y$$
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
}

#[derive(Debug)]
struct FloorPowerSum<R>
where
    R: SemiRing,
{
    x: R::T,
    sum: R::T,
}

impl<R> Clone for FloorPowerSum<R>
where
    R: SemiRing,
{
    fn clone(&self) -> Self {
        Self {
            x: self.x.clone(),
            sum: self.sum.clone(),
        }
    }
}

impl<R> FloorPowerSum<R>
where
    R: SemiRing,
{
    fn to_x(x: R::T) -> Self {
        Self { x, sum: R::one() }
    }
    fn to_y(y: R::T) -> Self {
        Self {
            x: y,
            sum: R::zero(),
        }
    }
}

impl<R> Magma for FloorPowerSum<R>
where
    R: SemiRing,
{
    type T = Self;

    fn operate(a: &Self::T, b: &Self::T) -> Self::T {
        Self {
            x: R::mul(&a.x, &b.x),
            sum: R::add(&a.sum, &R::mul(&a.x, &b.sum)),
        }
    }
}

impl<R> Unital for FloorPowerSum<R>
where
    R: SemiRing,
{
    fn unit() -> Self::T {
        Self {
            x: R::one(),
            sum: R::zero(),
        }
    }

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

    pub fn eye(shape: (usize, usize)) -> Self {
        let mut data = vec![vec![R::zero(); shape.1]; shape.0];
        for (i, d) in data.iter_mut().enumerate().take(shape.1) {
            d[i] = R::one();
        }
        Self {
            shape,
            data,
            _marker: PhantomData,
        }
    }

    pub fn transpose(&self) -> Self {
        Self::new_with((self.shape.1, self.shape.0), |i, j| self[j][i].clone())
    }

    pub fn map<S, F>(&self, mut f: F) -> Matrix<S>
    where
        S: SemiRing,
        F: FnMut(&R::T) -> S::T,
    {
        Matrix::<S>::new_with(self.shape, |i, j| f(&self[i][j]))
    }

    pub fn add_row_with(&mut self, mut f: impl FnMut(usize, usize) -> R::T) {
        self.data
            .push((0..self.shape.1).map(|j| f(self.shape.0, j)).collect());
        self.shape.0 += 1;
    }

    pub fn add_col_with(&mut self, mut f: impl FnMut(usize, usize) -> R::T) {
        for i in 0..self.shape.0 {
            self.data[i].push(f(i, self.shape.1));
        }
        self.shape.1 += 1;
    }

    pub fn pairwise_assign<F>(&mut self, other: &Self, mut f: F)
    where
        F: FnMut(&mut R::T, &R::T),
    {
        assert_eq!(self.shape, other.shape);
        for i in 0..self.shape.0 {
            for j in 0..self.shape.1 {
                f(&mut self[i][j], &other[i][j]);
            }
        }
    }
}

#[derive(Debug)]
pub struct SystemOfLinearEquationsSolution<R>
where
    R: Field<Additive: Invertible, Multiplicative: Invertible>,
{
    pub particular: Vec<R::T>,
    pub basis: Vec<Vec<R::T>>,
}

impl<R> Matrix<R>
where
    R: Field<T: PartialEq, Additive: Invertible, Multiplicative: Invertible>,
{
    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()
    }

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

    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
    }

Additional examples can be found in:

• crates/competitive/src/algorithm/automata_learning.rs

fn is_one(x: &Self::T) -> bool

where Self::T: PartialEq,

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

additive operaion: $+$

Examples found in repository

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

    fn mul1_add(&self, x: &R::T, y: &R::T) -> R::T {
        R::add(&R::mul(x, &self.base), y)
    }
    fn muln_add(&mut self, x: &R::T, y: &R::T, n: usize) -> R::T {
        if let Some(pow) = self.pow.get(n) {
            R::add(&R::mul(x, pow), y)
        } else {
            let pow = <R::Multiplicative as Monoid>::pow(self.base.clone(), n);
            R::add(&R::mul(x, &pow), y)
        }
    }

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

    fn operate(a: &Self::T, b: &Self::T) -> Self::T {
        Self {
            x: R::mul(&a.x, &b.x),
            sum: R::add(&a.sum, &R::mul(&a.x, &b.sum)),
        }
    }

crates/competitive/src/graph/shortest_path.rs (line 110)

    fn add_assign(x: &mut Self::T, y: &Self::T) -> bool {
        match x.0.cmp(&y.0) {
            Ordering::Equal => {
                x.1 = S::add(&x.1, &y.1);
                false
            }
            Ordering::Greater => {
                *x = y.clone();
                true
            }
            _ => false,
        }
    }

crates/competitive/src/algorithm/esper.rs (line 152)

    pub fn solve(&self, input: Input) -> R::T {
        let coeff = self
            .data
            .get(&(self.class)(&input))
            .expect("unrecognized class")
            .as_ref()
            .expect("failed to solve");
        let feature = (self.feature)(&input);
        feature
            .into_iter()
            .zip(coeff)
            .map(|(x, y)| R::mul(&x, y))
            .fold(R::zero(), |x, y| R::add(&x, &y))
    }

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

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 103)

    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 mul(x: &Self::T, y: &Self::T) -> Self::T

multiplicative operaion: $+$

Examples found in repository

crates/competitive/src/graph/shortest_path.rs (line 104)

    fn mul(x: &Self::T, y: &Self::T) -> Self::T {
        PartialIgnoredOrd(M::operate(&x.0, &y.0), S::mul(&x.1, &y.1))
    }

crates/competitive/src/math/bitwiseand_convolve.rs (line 51)

    fn multiply(f: &mut Self::F, g: &Self::F) {
        for (f, g) in f.iter_mut().zip(g) {
            *f = R::mul(f, g);
        }
    }

    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 a = Self::transform(a, len);
        if same {
            for a in a.iter_mut() {
                *a = R::mul(a, a);
            }
        } else {
            let b = Self::transform(b, len);
            Self::multiply(&mut a, &b);
        }
        Self::inverse_transform(a, len)
    }

crates/competitive/src/math/bitwiseor_convolve.rs (line 51)

    fn multiply(f: &mut Self::F, g: &Self::F) {
        for (f, g) in f.iter_mut().zip(g) {
            *f = R::mul(f, g);
        }
    }

    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 a = Self::transform(a, len);
        if same {
            for a in a.iter_mut() {
                *a = R::mul(a, a);
            }
        } else {
            let b = Self::transform(b, len);
            Self::multiply(&mut a, &b);
        }
        Self::inverse_transform(a, len)
    }

crates/competitive/src/math/gcd_convolve.rs (line 65)

    fn multiply(f: &mut Self::F, g: &Self::F) {
        for (f, g) in f.iter_mut().zip(g) {
            *f = R::mul(f, g);
        }
    }

crates/competitive/src/math/lcm_convolve.rs (line 65)

    fn multiply(f: &mut Self::F, g: &Self::F) {
        for (f, g) in f.iter_mut().zip(g) {
            *f = R::mul(f, g);
        }
    }

crates/competitive/src/algebra/ring.rs (line 48)

    fn dot_product(x: &[Self::T], y: &[Self::T]) -> Self::T {
        assert_eq!(x.len(), y.len());
        x.iter().zip(y).fold(Self::zero(), |mut sum, (x, y)| {
            Self::add_assign(&mut sum, &Self::mul(x, y));
            sum
        })
    }

    fn add_scaled_assign(x: &mut [Self::T], y: &[Self::T], a: &Self::T) {
        assert_eq!(x.len(), y.len());
        for (x, y) in x.iter_mut().zip(y) {
            Self::add_assign(x, &Self::mul(a, y));
        }
    }

Additional examples can be found in:

• crates/competitive/src/string/rolling_hash.rs
• crates/competitive/src/algorithm/esper.rs
• crates/competitive/src/math/bitwisexor_convolve.rs
• crates/competitive/src/math/quotient_array.rs
• crates/competitive/src/math/matrix.rs
• crates/competitive/src/math/floor_sum.rs
• crates/competitive/src/math/black_box_matrix.rs
• crates/competitive/src/algorithm/automata_learning.rs

fn try_matrix_product( _a: &[Vec<Self::T>], _b: &[Vec<Self::T>], ) -> Option<Vec<Vec<Self::T>>>

Examples found in repository

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

    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
    }

fn dot_product(x: &[Self::T], y: &[Self::T]) -> Self::T

Examples found in repository

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

    fn apply(&self, v: &[R::T]) -> Vec<R::T> {
        assert_eq!(self.shape.1, v.len());
        self.data.iter().map(|row| R::dot_product(row, v)).collect()
    }

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

    fn multiply_row(
        x: &[R::T],
        y: &[R::T],
        right: &mut [R::T],
        output: &mut [R::T],
        rank: usize,
    ) -> usize {
        for (right, y) in right[..=rank].iter_mut().zip(y[..=rank].iter().rev()) {
            right.clone_from(y);
        }
        let end = (rank * 2).min(x.len() - 1);
        for (degree, output) in output.iter_mut().enumerate().take(end + 1).skip(rank) {
            let first = degree - rank;
            *output = R::dot_product(&x[first..=rank], &right[..=rank - first]);
        }
        end
    }

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

    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);
    }
}

fn add_scaled_assign(x: &mut [Self::T], y: &[Self::T], a: &Self::T)

Examples found in repository

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

    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()
    }

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

Examples found in repository

crates/competitive/src/algebra/ring.rs (line 48)

    fn dot_product(x: &[Self::T], y: &[Self::T]) -> Self::T {
        assert_eq!(x.len(), y.len());
        x.iter().zip(y).fold(Self::zero(), |mut sum, (x, y)| {
            Self::add_assign(&mut sum, &Self::mul(x, y));
            sum
        })
    }

    fn add_scaled_assign(x: &mut [Self::T], y: &[Self::T], a: &Self::T) {
        assert_eq!(x.len(), y.len());
        for (x, y) in x.iter_mut().zip(y) {
            Self::add_assign(x, &Self::mul(a, y));
        }
    }

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

    fn transform(t: Self::T, len: usize) -> Self::F {
        let (mut f, width) = Self::ranked(t, len);
        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 ranks = lower.count_ones() as usize + 1;
                    let (lower_rows, upper_rows) = f.split_at_mut(upper * width);
                    let lower_row = &lower_rows[lower * width..lower * width + ranks];
                    let upper_row = &mut upper_rows[..ranks];
                    for (upper, lower) in upper_row.iter_mut().zip(lower_row) {
                        R::add_assign(upper, lower);
                    }
                }
            }
        }
        (f, width)
    }

    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/floor_sum.rs (line 268)

    fn operate(a: &Self::T, b: &Self::T) -> Self::T {
        let mut a = a.clone();
        let mut b = b.clone();
        let mut pow_x = array![R::zero(); X];
        let mut pow_y = array![R::zero(); Y];
        pow_x[0] = R::one();
        pow_y[0] = R::one();
        for i in 1..X {
            pow_x[i] = R::mul(&pow_x[i - 1], &a.dx);
        }
        for j in 1..Y {
            pow_y[j] = R::mul(&pow_y[j - 1], &a.dy);
        }
        macro_rules! go {
            ($N:ident) => {
                let mut comb = array![array![R::zero(); $N]; $N];
                comb[0][0] = R::one();
                let mut i = 0;
                while i + 1 < $N {
                    let mut j = 0;
                    while j <= i {
                        let x = comb[i][j].clone();
                        R::add_assign(&mut comb[i + 1][j], &x);
                        R::add_assign(&mut comb[i + 1][j + 1], &x);
                        j += 1;
                    }
                    i += 1;
                }
                for i in 0..X {
                    for j in (0..Y).rev() {
                        for k in j + 1..Y {
                            let mut x = b.dp[i][j].clone();
                            R::mul_assign(&mut x, &comb[k][j]);
                            R::mul_assign(&mut x, &pow_y[k - j]);
                            R::add_assign(&mut b.dp[i][k], &x);
                        }
                    }
                }
                for j in 0..Y {
                    for i in (0..X).rev() {
                        for k in i..X {
                            let mut x = b.dp[i][j].clone();
                            R::mul_assign(&mut x, &comb[k][i]);
                            R::mul_assign(&mut x, &pow_x[k - i]);
                            R::add_assign(&mut a.dp[k][j], &x);
                        }
                    }
                }
            };
        }
        if X <= Y {
            go!(Y);
        } else {
            go!(X);
        }
        R::add_assign(&mut a.dx, &b.dx);
        R::add_assign(&mut a.dy, &b.dy);
        a
    }

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

    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);
    }
}

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

    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
    }
}

impl<R> From<Matrix<R>> for SparseMatrix<R>
where
    R: SemiRing<T: PartialEq>,
{
    fn from(mat: Matrix<R>) -> Self {
        let mut nonzero = vec![];
        for i in 0..mat.shape.0 {
            for j in 0..mat.shape.1 {
                let v = mat[(i, j)].clone();
                if !R::is_zero(&v) {
                    nonzero.push((i, j, v));
                }
            }
        }
        Self {
            shape: mat.shape,
            nonzero,
        }
    }
}

impl<R> From<SparseMatrix<R>> for Matrix<R>
where
    R: SemiRing,
{
    fn from(smat: SparseMatrix<R>) -> Self {
        let mut mat = Matrix::zeros(smat.shape);
        for &(i, j, ref v) in &smat.nonzero {
            R::add_assign(&mut mat[(i, j)], v);
        }
        mat
    }
}

impl<R> BlackBoxMatrix<R> for SparseMatrix<R>
where
    R: SemiRing,
{
    fn apply(&self, v: &[R::T]) -> Vec<R::T> {
        assert_eq!(self.shape.1, v.len());
        let mut res = vec![R::zero(); self.shape.0];
        for &(i, j, ref val) in &self.nonzero {
            R::add_assign(&mut res[i], &R::mul(val, &v[j]));
        }
        res
    }

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

    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 mul_assign(x: &mut Self::T, y: &Self::T)

Examples found in repository

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

    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);
            }
        }
    }

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

    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
    }

Dyn Compatibility

This trait is not dyn compatible.

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

Implementors

impl SemiRing for Gf2_63

type T = u64

type Additive = BitXorOperation<u64>

type Multiplicative = Gf2_63

impl SemiRing for Mersenne61

type T = u64

type Additive = Mersenne61Add

type Multiplicative = Mersenne61

impl<T> SemiRing for AddMulOperation<T>

where T: DotProduct + One,

type T = T

type Additive = AdditiveOperation<T>

type Multiplicative = MultiplicativeOperation<T>