Struct Matrix

Source

pub struct Matrix<R>where
    R: SemiRing,{
    pub shape: (usize, usize),
    pub data: Vec<Vec<R::T>>,
    _marker: PhantomData<fn() -> R>,
}

Fields§

§shape: (usize, usize)§data: Vec<Vec<R::T>>§_marker: PhantomData<fn() -> R>

Implementations§

Source §

impl<R> Matrix<R>
where R: SemiRing,

Source

pub fn new(shape: (usize, usize), z: R::T) -> Self

Source

pub fn from_vec(data: Vec<Vec<R::T>>) -> Self

Examples found in repository ?

crates/competitive/src/math/mod.rs (line 325)

321    fn deserialize<I>(iter: &mut I) -> Self
322    where
323        I: Iterator<Item = u8>,
324    {
325        Self::from_vec(Vec::deserialize(iter))
326    }

More examples

Hide additional examples

crates/library_checker/src/linear_algebra/matrix_det.rs (line 8)

5pub fn matrix_det(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, a: [[M; n]; n]);
8    let mut a = Matrix::<AddMulOperation<_>>::from_vec(a);
9    let det = a.determinant();
10    pp!(det);
11}

crates/library_checker/src/linear_algebra/characteristic_polynomial.rs (line 8)

5pub fn characteristic_polynomial(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, a: [[M; n]; n]);
8    let p = Matrix::<AddMulOperation<_>>::from_vec(a).characteristic_polynomial();
9    pp!(@it p);
10}

crates/library_checker/src/linear_algebra/pow_of_matrix.rs (line 12)

9pub fn pow_of_matrix(reader: impl Read, writer: impl Write) {
10    prepare_io!(reader, writer);
11    sc!(n, k, a: [[M; n]; n]);
12    let a = Matrix::<AddMulOperation<_>>::from_vec(a);
13    let b = a.pow_frobenius(k);
14    pp!(@it2d b.data);
15}
16
17#[verify::library_checker("pow_of_matrix")]
18pub fn pow_of_matrix_strassen(reader: impl Read, writer: impl Write) {
19    prepare_io!(reader, writer);
20    sc!(n, k, a: [[M; n]; n]);
21    let a = Matrix::<AddMulOperation<_>>::from_vec(a);
22    let b = a.pow_strassen(k);
23    pp!(@it2d b.data);
24}

crates/library_checker/src/linear_algebra/inverse_matrix.rs (line 8)

5pub fn inverse_matrix(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, a: [[M; n]; n]);
8    let a = Matrix::<AddMulOperation<_>>::from_vec(a);
9    if let Some(b) = a.inverse() {
10        pp!(@it2d b.data);
11    } else {
12        pp!("-1");
13    }
14}

crates/library_checker/src/linear_algebra/matrix_product.rs (line 8)

5pub fn matrix_product(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, m, k, a: [[M; m]; n], b: [[M; k]; m]);
8    let a = Matrix::<AddMulOperation<_>>::from_vec(a);
9    let b = Matrix::<AddMulOperation<_>>::from_vec(b);
10    let c = a * b;
11    pp!(@it2d c.data);
12}
13
14#[verify::library_checker("matrix_product")]
15pub fn matrix_product_strassen(reader: impl Read, writer: impl Write) {
16    prepare_io!(reader, writer);
17    sc!(n, m, k, a: [[M; m]; n], b: [[M; k]; m]);
18    let a = Matrix::<AddMulOperation<_>>::from_vec(a);
19    let b = Matrix::<AddMulOperation<_>>::from_vec(b);
20    let c = a.mul_strassen(&b);
21    pp!(@it2d c.data);
22}

Additional examples can be found in:

Source

pub fn new_with( shape: (usize, usize), f: impl FnMut(usize, usize) -> R::T, ) -> Self

Examples found in repository ?

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

107    pub fn transpose(&self) -> Self {
108        Self::new_with((self.shape.1, self.shape.0), |i, j| self[j][i].clone())
109    }
110
111    pub fn map<S, F>(&self, mut f: F) -> Matrix<S>
112    where
113        S: SemiRing,
114        F: FnMut(&R::T) -> S::T,
115    {
116        Matrix::<S>::new_with(self.shape, |i, j| f(&self[i][j]))
117    }
118
119    pub fn add_row_with(&mut self, mut f: impl FnMut(usize, usize) -> R::T) {
120        self.data
121            .push((0..self.shape.1).map(|j| f(self.shape.0, j)).collect());
122        self.shape.0 += 1;
123    }
124
125    pub fn add_col_with(&mut self, mut f: impl FnMut(usize, usize) -> R::T) {
126        for i in 0..self.shape.0 {
127            self.data[i].push(f(i, self.shape.1));
128        }
129        self.shape.1 += 1;
130    }
131
132    pub fn pairwise_assign<F>(&mut self, other: &Self, mut f: F)
133    where
134        F: FnMut(&mut R::T, &R::T),
135    {
136        assert_eq!(self.shape, other.shape);
137        for i in 0..self.shape.0 {
138            for j in 0..self.shape.1 {
139                f(&mut self[i][j], &other[i][j]);
140            }
141        }
142    }
143}
144
145#[derive(Debug)]
146pub struct SystemOfLinearEquationsSolution<R>
147where
148    R: Field<Additive: Invertible, Multiplicative: Invertible>,
149{
150    pub particular: Vec<R::T>,
151    pub basis: Vec<Vec<R::T>>,
152}
153
154impl<R> Matrix<R>
155where
156    R: Field<T: PartialEq, Additive: Invertible, Multiplicative: Invertible>,
157{
158    fn eliminate<const DETERMINANT: bool>(&mut self) -> (usize, R::T) {
159        let (n, m) = self.shape;
160        let mut rank = 0;
161        let mut determinant = R::one();
162        let mut negative = false;
163        for first in (0..m).step_by(64) {
164            if rank == n {
165                break;
166            }
167            let end = (first + 64).min(m);
168            let start = rank;
169            let mut pivots = Vec::new();
170            let mut panel = vec![vec![R::zero(); end - first]; end - first];
171            for col in first..end {
172                if panel[col - first][..col - first]
173                    .iter()
174                    .any(|x| !R::is_zero(x))
175                {
176                    for row in &mut self.data[rank..] {
177                        let value =
178                            R::dot_product(&row[first..col], &panel[col - first][..col - first]);
179                        R::sub_assign(&mut row[col], &value);
180                    }
181                }
182                let Some(pivot) = (rank..n).find(|&i| !R::is_zero(&self[i][col])) else {
183                    continue;
184                };
185                if pivot != rank {
186                    self.data.swap(rank, pivot);
187                    negative = !negative;
188                }
189                if DETERMINANT {
190                    R::mul_assign(&mut determinant, &self[rank][col]);
191                }
192                let inv = R::inv(&self[rank][col]);
193                let row = &mut self.data[rank];
194                for c in col + 1..end {
195                    let value = R::dot_product(&row[first..col], &panel[c - first][..col - first]);
196                    R::sub_assign(&mut row[c], &value);
197                    panel[c - first][col - first] = row[c].clone();
198                }
199                for row in &mut self.data[rank + 1..] {
200                    R::mul_assign(&mut row[col], &inv);
201                }
202                pivots.push(col);
203                rank += 1;
204                if rank == n {
205                    break;
206                }
207            }
208            for i in start..if rank - start < 32 { n } else { rank } {
209                let (upper, lower) = self.data.split_at_mut(i);
210                let row = &mut lower[0];
211                for (j, &col) in pivots[..(i - start).min(pivots.len())].iter().enumerate() {
212                    if R::is_zero(&row[col]) {
213                        continue;
214                    }
215                    let factor = R::neg(&row[col]);
216                    R::add_scaled_assign(&mut row[end..], &upper[start + j][end..], &factor);
217                }
218            }
219            if rank < n
220                && end < m
221                && rank - start >= 32
222                && self.data[rank..]
223                    .iter()
224                    .any(|row| pivots.iter().any(|&col| !R::is_zero(&row[col])))
225            {
226                let lower = Self::new_with((n - rank, rank - start), |i, j| {
227                    R::neg(&self[rank + i][pivots[j]])
228                });
229                let upper = Self::new_with((rank - start, m - end), |i, j| {
230                    self[start + i][end + j].clone()
231                });
232                let update = &lower * &upper;
233                for (row, update) in self.data[rank..].iter_mut().zip(update.data) {
234                    for (x, y) in row[end..].iter_mut().zip(update) {
235                        R::add_assign(x, &y);
236                    }
237                }
238            }
239            for (i, &col) in pivots.iter().enumerate() {
240                for row in &mut self.data[start + i + 1..] {
241                    row[col] = R::zero();
242                }
243            }
244            if DETERMINANT && rank < end {
245                return (rank, R::zero());
246            }
247        }
248        if DETERMINANT && negative {
249            determinant = R::neg(&determinant);
250        }
251        (rank, determinant)
252    }
253
254    /// f: (row, pivot_row, col)
255    pub fn row_reduction_with<F>(&mut self, normalize: bool, mut f: F)
256    where
257        F: FnMut(usize, usize, usize),
258    {
259        let (n, m) = self.shape;
260        let mut c = 0;
261        for r in 0..n {
262            loop {
263                if c >= m {
264                    return;
265                }
266                if let Some(pivot) = (r..n).find(|&p| !R::is_zero(&self[p][c])) {
267                    f(r, pivot, c);
268                    self.data.swap(r, pivot);
269                    break;
270                };
271                c += 1;
272            }
273            let d = R::inv(&self[r][c]);
274            if normalize {
275                for value in &mut self[r][c..m] {
276                    R::mul_assign(value, &d);
277                }
278            }
279            for i in (0..n).filter(|&i| i != r) {
280                let mut e = self[i][c].clone();
281                if !normalize {
282                    R::mul_assign(&mut e, &d);
283                }
284                for j in c..m {
285                    let e = R::mul(&e, &self[r][j]);
286                    R::sub_assign(&mut self[i][j], &e);
287                }
288            }
289            c += 1;
290        }
291    }
292
293    pub fn row_reduction(&mut self, normalize: bool) {
294        self.row_reduction_with(normalize, |_, _, _| {});
295    }
296
297    pub fn rank(&mut self) -> usize {
298        self.eliminate::<false>().0
299    }
300
301    pub fn determinant(&mut self) -> R::T {
302        assert_eq!(self.shape.0, self.shape.1);
303        self.eliminate::<true>().1
304    }
305
306    pub fn solve_system_of_linear_equations(
307        &self,
308        b: &[R::T],
309    ) -> Option<SystemOfLinearEquationsSolution<R>> {
310        assert_eq!(self.shape.0, b.len());
311        let m = self.shape.1;
312        let mut a = Self::new_with((self.shape.0, m + 1), |i, j| {
313            if j == m {
314                b[i].clone()
315            } else {
316                self[i][j].clone()
317            }
318        });
319        let rank = a.eliminate::<false>().0;
320        let mut pivots = Vec::with_capacity(rank);
321        let mut b = Vec::with_capacity(rank);
322        for row in &a.data[..rank] {
323            let c = row.iter().position(|x| !R::is_zero(x)).unwrap();
324            if c == m {
325                return None;
326            }
327            pivots.push(c);
328            b.push(row[m].clone());
329        }
330
331        let mut free = Vec::with_capacity(m - rank);
332        let mut pivot = 0;
333        for c in 0..m {
334            if pivot < rank && pivots[pivot] == c {
335                pivot += 1;
336            } else {
337                free.push(c);
338            }
339        }
340        let mut coefficients: Vec<Vec<_>> = (0..rank)
341            .map(|i| free.iter().map(|&c| a[i][c].clone()).collect())
342            .collect();
343        for k in (0..rank).rev() {
344            let c = pivots[k];
345            let inv = R::inv(&a[k][c]);
346            R::mul_assign(&mut b[k], &inv);
347            let pivot_b = b[k].clone();
348            let (upper, lower) = coefficients.split_at_mut(k);
349            let pivot_coefficients = &mut lower[0];
350            for x in pivot_coefficients.iter_mut() {
351                R::mul_assign(x, &inv);
352            }
353            for ((row, value), coefficients) in a.data[..k].iter_mut().zip(&mut b[..k]).zip(upper) {
354                if R::is_zero(&row[c]) {
355                    continue;
356                }
357                let factor = row[c].clone();
358                row[c] = R::zero();
359                R::sub_assign(value, &R::mul(&factor, &pivot_b));
360                R::add_scaled_assign(coefficients, pivot_coefficients, &R::neg(&factor));
361            }
362        }
363
364        let mut particular = vec![R::zero(); m];
365        for i in 0..rank {
366            particular[pivots[i]] = b[i].clone();
367        }
368        let mut basis = Vec::with_capacity(free.len());
369        for (j, &c) in free.iter().enumerate() {
370            let mut vector = vec![R::zero(); m];
371            vector[c] = R::one();
372            for i in 0..rank {
373                vector[pivots[i]] = R::neg(&coefficients[i][j]);
374            }
375            basis.push(vector);
376        }
377        Some(SystemOfLinearEquationsSolution { particular, basis })
378    }
379
380    pub fn inverse(&self) -> Option<Matrix<R>> {
381        assert_eq!(self.shape.0, self.shape.1);
382        let n = self.shape.0;
383        if n >= 64 {
384            let m = n / 2;
385            let a = Self::new_with((m, m), |i, j| self[i][j].clone());
386            if let Some(mut ai) = a.inverse() {
387                let b = Self::new_with((m, n - m), |i, j| self[i][j + m].clone());
388                let c = Self::new_with((n - m, m), |i, j| self[i + m][j].clone());
389                let mut d = Self::new_with((n - m, n - m), |i, j| self[i + m][j + m].clone());
390                let u = &ai * &b;
391                let v = &c * &ai;
392                d -= &v * &b;
393                let di = d.inverse()?;
394                let r = &u * &di;
395                let t = &di * &v;
396                ai += &r * &v;
397                let mut inverse = Self::zeros((n, n));
398                for i in 0..m {
399                    inverse[i][..m].clone_from_slice(&ai[i]);
400                    for (x, y) in inverse[i][m..].iter_mut().zip(&r[i]) {
401                        *x = R::neg(y);
402                    }
403                }
404                for i in m..n {
405                    for (x, y) in inverse[i][..m].iter_mut().zip(&t[i - m]) {
406                        *x = R::neg(y);
407                    }
408                    inverse[i][m..].clone_from_slice(&di[i - m]);
409                }
410                return Some(inverse);
411            }
412        }
413        let mut a = self.clone();
414        let mut inverse = Self::eye((n, n));
415        let mut ranges: Vec<_> = (0..n).map(|i| (i, i + 1)).collect();
416        for r in 0..n {
417            let pivot = (r..n).find(|&i| !R::is_zero(&a[i][r]))?;
418            a.data.swap(r, pivot);
419            inverse.data.swap(r, pivot);
420            ranges.swap(r, pivot);
421
422            let d = R::inv(&a[r][r]);
423            for x in &mut a[r][r..] {
424                R::mul_assign(x, &d);
425            }
426            let (left, right) = ranges[r];
427            for x in &mut inverse[r][left..right] {
428                R::mul_assign(x, &d);
429            }
430
431            let (a_upper, a_lower) = a.data.split_at_mut(r + 1);
432            let pivot_a = &a_upper[r];
433            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r + 1);
434            let pivot_inverse = &inverse_upper[r];
435            let (ranges_upper, ranges_lower) = ranges.split_at_mut(r + 1);
436            let (left, right) = ranges_upper[r];
437            for ((a, inverse), range) in a_lower.iter_mut().zip(inverse_lower).zip(ranges_lower) {
438                if R::is_zero(&a[r]) {
439                    continue;
440                }
441                let e = a[r].clone();
442                a[r] = R::zero();
443                R::add_scaled_assign(&mut a[(r + 1)..], &pivot_a[(r + 1)..], &R::neg(&e));
444                R::add_scaled_assign(
445                    &mut inverse[left..right],
446                    &pivot_inverse[left..right],
447                    &R::neg(&e),
448                );
449                range.0 = range.0.min(left);
450                range.1 = range.1.max(right);
451            }
452        }
453        for r in (0..n).rev() {
454            let (left, right) = ranges[r];
455            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r);
456            let pivot_inverse = &inverse_lower[0];
457            let (ranges_upper, _) = ranges.split_at_mut(r);
458            for ((a, inverse), range) in a.data[..r].iter_mut().zip(inverse_upper).zip(ranges_upper)
459            {
460                if R::is_zero(&a[r]) {
461                    continue;
462                }
463                let e = a[r].clone();
464                a[r] = R::zero();
465                R::add_scaled_assign(
466                    &mut inverse[left..right],
467                    &pivot_inverse[left..right],
468                    &R::neg(&e),
469                );
470                range.0 = range.0.min(left);
471                range.1 = range.1.max(right);
472            }
473        }
474        Some(inverse)
475    }
476
477    pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
478        let n = self.shape.0;
479        if n == 0 {
480            return vec![R::one()];
481        }
482        assert!(self.data.iter().all(|a| a.len() == n));
483        for j in 0..(n - 1) {
484            if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
485                self.data.swap(j + 1, x);
486                self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
487                let inv = R::inv(&self[j + 1][j]);
488                let mut v = vec![];
489                let src = std::mem::take(&mut self[j + 1]);
490                for a in self.data[(j + 2)..].iter_mut() {
491                    let mul = R::mul(&a[j], &inv);
492                    R::add_scaled_assign(&mut a[j..], &src[j..], &R::neg(&mul));
493                    v.push(mul);
494                }
495                self[j + 1] = src;
496                for a in self.data.iter_mut() {
497                    let v = R::dot_product(&a[(j + 2)..], &v);
498                    R::add_assign(&mut a[j + 1], &v);
499                }
500            }
501        }
502        // dp[k][j - k] stores [x^k] det(xI - A[..j, ..j]).
503        let mut dp: Vec<Vec<R::T>> = (0..=n).map(|i| Vec::with_capacity(n + 1 - i)).collect();
504        dp[0].push(R::one());
505        for i in 0..n {
506            let mut c = vec![R::zero(); i + 1];
507            c[i] = R::neg(&self[i][i]);
508            let mut mul = R::one();
509            for j in (0..i).rev() {
510                mul = R::mul(&mul, &self[j + 1][j]);
511                c[j] = R::neg(&R::mul(&mul, &self[j][i]));
512            }
513            for k in (0..=i).rev() {
514                let mut value = R::dot_product(&dp[k], &c[k..]);
515                if k > 0 {
516                    R::add_assign(&mut value, dp[k - 1].last().unwrap());
517                }
518                dp[k].push(value);
519            }
520            dp[i + 1].push(R::one());
521        }
522        dp.into_iter().map(|mut c| c.pop().unwrap()).collect()
523    }
524}
525
526impl<R> Index<usize> for Matrix<R>
527where
528    R: SemiRing,
529{
530    type Output = Vec<R::T>;
531    fn index(&self, index: usize) -> &Self::Output {
532        &self.data[index]
533    }
534}
535
536impl<R> IndexMut<usize> for Matrix<R>
537where
538    R: SemiRing,
539{
540    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
541        &mut self.data[index]
542    }
543}
544
545impl<R> Index<(usize, usize)> for Matrix<R>
546where
547    R: SemiRing,
548{
549    type Output = R::T;
550    fn index(&self, index: (usize, usize)) -> &Self::Output {
551        &self.data[index.0][index.1]
552    }
553}
554
555impl<R> IndexMut<(usize, usize)> for Matrix<R>
556where
557    R: SemiRing,
558{
559    fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
560        &mut self.data[index.0][index.1]
561    }
562}
563
564macro_rules! impl_matrix_pairwise_binop {
565    ($imp:ident, $method:ident, $imp_assign:ident, $method_assign:ident $(where [$($clauses:tt)*])?) => {
566        impl<R> $imp_assign for Matrix<R>
567        where
568            R: SemiRing,
569            $($($clauses)*)?
570        {
571            fn $method_assign(&mut self, rhs: Self) {
572                self.pairwise_assign(&rhs, |a, b| R::$method_assign(a, b));
573            }
574        }
575        impl<R> $imp_assign<&Matrix<R>> for Matrix<R>
576        where
577            R: SemiRing,
578            $($($clauses)*)?
579        {
580            fn $method_assign(&mut self, rhs: &Self) {
581                self.pairwise_assign(rhs, |a, b| R::$method_assign(a, b));
582            }
583        }
584        impl<R> $imp for Matrix<R>
585        where
586            R: SemiRing,
587            $($($clauses)*)?
588        {
589            type Output = Matrix<R>;
590            fn $method(mut self, rhs: Self) -> Self::Output {
591                self.$method_assign(rhs);
592                self
593            }
594        }
595        impl<R> $imp<&Matrix<R>> for Matrix<R>
596        where
597            R: SemiRing,
598            $($($clauses)*)?
599        {
600            type Output = Matrix<R>;
601            fn $method(mut self, rhs: &Self) -> Self::Output {
602                self.$method_assign(rhs);
603                self
604            }
605        }
606        impl<R> $imp<Matrix<R>> for &Matrix<R>
607        where
608            R: SemiRing,
609            $($($clauses)*)?
610        {
611            type Output = Matrix<R>;
612            fn $method(self, mut rhs: Matrix<R>) -> Self::Output {
613                rhs.pairwise_assign(self, |a, b| *a = R::$method(b, a));
614                rhs
615            }
616        }
617        impl<R> $imp<&Matrix<R>> for &Matrix<R>
618        where
619            R: SemiRing,
620            $($($clauses)*)?
621        {
622            type Output = Matrix<R>;
623            fn $method(self, rhs: &Matrix<R>) -> Self::Output {
624                let mut this = self.clone();
625                this.$method_assign(rhs);
626                this
627            }
628        }
629    };
630}
631
632impl_matrix_pairwise_binop!(Add, add, AddAssign, add_assign);
633impl_matrix_pairwise_binop!(Sub, sub, SubAssign, sub_assign where [R: SemiRing<Additive: Invertible>]);
634
635impl<R> Mul for Matrix<R>
636where
637    R: SemiRing,
638{
639    type Output = Matrix<R>;
640    fn mul(self, rhs: Self) -> Self::Output {
641        (&self).mul(&rhs)
642    }
643}
644impl<R> Mul<&Matrix<R>> for Matrix<R>
645where
646    R: SemiRing,
647{
648    type Output = Matrix<R>;
649    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
650        (&self).mul(rhs)
651    }
652}
653impl<R> Mul<Matrix<R>> for &Matrix<R>
654where
655    R: SemiRing,
656{
657    type Output = Matrix<R>;
658    fn mul(self, rhs: Matrix<R>) -> Self::Output {
659        self.mul(&rhs)
660    }
661}
662impl<R> Mul<&Matrix<R>> for &Matrix<R>
663where
664    R: SemiRing,
665{
666    type Output = Matrix<R>;
667    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
668        assert_eq!(self.shape.1, rhs.shape.0);
669        if let Some(data) = R::try_matrix_product(&self.data, &rhs.data) {
670            return Matrix::from_vec(data);
671        }
672        let rhs = rhs.transpose();
673        Matrix::new_with((self.shape.0, rhs.shape.0), |i, j| {
674            R::dot_product(&self[i], &rhs[j])
675        })
676    }

More examples

Hide additional examples

crates/competitive/src/algorithm/automata_learning.rs (lines 532-539)

523    pub fn train_sample(&mut self, sample: &[usize]) -> bool {
524        let Some((prefix, suffix)) = self.split_sample(sample) else {
525            return false;
526        };
527        self.prefixes.push(prefix);
528        self.suffixes.push(suffix);
529        let n = self.inv_h.shape.0;
530        let prefix = &self.prefixes[n];
531        let suffix = &self.suffixes[n];
532        let u = Matrix::<F>::new_with((n, 1), |i, _| {
533            self.automaton.behavior(
534                self.prefixes[i]
535                    .iter()
536                    .cloned()
537                    .chain(suffix.iter().cloned()),
538            )
539        });
540        let v = Matrix::<F>::new_with((1, n), |_, j| {
541            self.automaton.behavior(
542                prefix
543                    .iter()
544                    .cloned()
545                    .chain(self.suffixes[j].iter().cloned()),
546            )
547        });
548        let w = Matrix::<F>::new_with((1, 1), |_, _| {
549            self.automaton
550                .behavior(prefix.iter().cloned().chain(suffix.iter().cloned()))
551        });
552        let t = &self.inv_h * &u;
553        let s = &v * &self.inv_h;
554        let d = F::inv(&(&w - &(&v * &t))[0][0]);
555        let dh = &t * &s;
556        for i in 0..n {
557            for j in 0..n {
558                F::add_assign(&mut self.inv_h[i][j], &F::mul(&dh[i][j], &d));
559            }
560        }
561        self.inv_h
562            .add_col_with(|i, _| F::neg(&F::mul(&t[i][0], &d)));
563        self.inv_h.add_row_with(|_, j| {
564            if j != n {
565                F::neg(&F::mul(&s[0][j], &d))
566            } else {
567                d.clone()
568            }
569        });
570
571        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
572            let b = &(&self.nh[x] * &t) * &s;
573            for i in 0..n {
574                for j in 0..n {
575                    F::add_assign(&mut transition[i][j], &F::mul(&b[i][j], &d));
576                }
577            }
578        }
579        for (x, nh) in self.nh.iter_mut().enumerate() {
580            nh.add_col_with(|i, j| {
581                self.automaton.behavior(
582                    self.prefixes[i]
583                        .iter()
584                        .cloned()
585                        .chain([x])
586                        .chain(self.suffixes[j].iter().cloned()),
587                )
588            });
589            nh.add_row_with(|i, j| {
590                self.automaton.behavior(
591                    self.prefixes[i]
592                        .iter()
593                        .cloned()
594                        .chain([x])
595                        .chain(self.suffixes[j].iter().cloned()),
596                )
597            });
598        }
599        self.wfa
600            .initial_weights
601            .add_col_with(|_, _| if n == 0 { F::one() } else { F::zero() });
602        self.wfa
603            .final_weights
604            .add_row_with(|_, _| self.automaton.behavior(prefix.iter().cloned()));
605        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
606            transition.add_col_with(|_, _| F::zero());
607            transition.add_row_with(|_, _| F::zero());
608            for i in 0..=n {
609                for j in 0..=n {
610                    if i == n || j == n {
611                        for k in 0..=n {
612                            if i != n && j != n && k != n {
613                                continue;
614                            }
615                            F::add_assign(
616                                &mut transition[i][k],
617                                &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
618                            );
619                        }
620                    } else {
621                        let k = n;
622                        F::add_assign(
623                            &mut transition[i][k],
624                            &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
625                        );
626                    }
627                }
628            }
629        }
630        true
631    }
632    pub fn train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
633        for sample in samples {
634            self.train_sample(&sample);
635        }
636    }
637    pub fn batch_train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
638        let mut prefix_set: HashSet<_> = self.prefixes.iter().cloned().collect();
639        let mut suffix_set: HashSet<_> = self.suffixes.iter().cloned().collect();
640        for sample in samples {
641            if prefix_set.insert(sample.to_vec()) {
642                self.prefixes.push(sample.to_vec());
643            }
644            if suffix_set.insert(sample.to_vec()) {
645                self.suffixes.push(sample);
646            }
647        }
648        let mut h = Matrix::<F>::new_with((self.prefixes.len(), self.suffixes.len()), |i, j| {
649            self.automaton.behavior(
650                self.prefixes[i]
651                    .iter()
652                    .cloned()
653                    .chain(self.suffixes[j].iter().cloned()),
654            )
655        });
656        if !self.prefixes.is_empty() && !self.suffixes.is_empty() && F::is_zero(&h[0][0]) {
657            for j in 1..self.suffixes.len() {
658                if !F::is_zero(&h[0][j]) {
659                    self.suffixes.swap(0, j);
660                    for row in &mut h.data {
661                        row.swap(0, j);
662                    }
663                    break;
664                }
665            }
666        }
667        let mut row_id: Vec<usize> = (0..h.shape.0).collect();
668        let mut pivots = vec![];
669        h.row_reduction_with(false, |r, p, c| {
670            row_id.swap(r, p);
671            pivots.push((row_id[r], c));
672        });
673        let mut new_prefixes = vec![];
674        let mut new_suffixes = vec![];
675        for (i, j) in pivots {
676            new_prefixes.push(self.prefixes[i].clone());
677            new_suffixes.push(self.suffixes[j].clone());
678        }
679        self.prefixes = new_prefixes;
680        self.suffixes = new_suffixes;
681        assert_eq!(self.prefixes.len(), self.suffixes.len());
682        let n = self.prefixes.len();
683        let h = Matrix::<F>::new_with((n, n), |i, j| {
684            self.automaton.behavior(
685                self.prefixes[i]
686                    .iter()
687                    .cloned()
688                    .chain(self.suffixes[j].iter().cloned()),
689            )
690        });
691        self.inv_h = h.inverse().expect("Hankel matrix must be invertible");
692        self.wfa = WeightedFiniteAutomaton::<F> {
693            initial_weights: Matrix::new_with((1, n), |_, j| {
694                if self.prefixes[j].is_empty() {
695                    F::one()
696                } else {
697                    F::zero()
698                }
699            }),
700            transitions: (0..self.automaton.sigma())
701                .map(|x| {
702                    &Matrix::new_with((n, n), |i, j| {
703                        self.automaton.behavior(
704                            self.prefixes[i]
705                                .iter()
706                                .cloned()
707                                .chain([x])
708                                .chain(self.suffixes[j].iter().cloned()),
709                        )
710                    }) * &self.inv_h
711                })
712                .collect(),
713            final_weights: Matrix::new_with((n, 1), |i, _| {
714                self.automaton.behavior(self.prefixes[i].iter().cloned())
715            }),
716        };
717    }

Source

pub fn zeros(shape: (usize, usize)) -> Self

Examples found in repository ?

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

224    fn from(smat: SparseMatrix<R>) -> Self {
225        let mut mat = Matrix::zeros(smat.shape);
226        for &(i, j, ref v) in &smat.nonzero {
227            R::add_assign(&mut mat[(i, j)], v);
228        }
229        mat
230    }

More examples

Hide additional examples

crates/library_checker/src/linear_algebra/matrix_rank.rs (line 12)

5pub fn matrix_rank(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, m);
8    let mut a = if n <= m {
9        sc!(a: [[M; m]; n]);
10        Matrix::<AddMulOperation<_>>::from_vec(a)
11    } else {
12        let mut a = Matrix::<AddMulOperation<_>>::zeros((m, n));
13        for j in 0..n {
14            for row in &mut a.data {
15                sc!(x: M);
16                row[j] = x;
17            }
18        }
19        a
20    };
21    let rank = a.rank();
22    pp!(rank);
23}

crates/library_checker/src/graph/counting_spanning_tree_directed.rs (line 8)

5pub fn counting_spanning_tree_directed(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, m, r: usize, edges: [(usize, usize); iter m]);
8    let mut a = Matrix::<AddMulOperation<M>>::zeros((n - 1, n - 1));
9    for (u, v) in edges {
10        if v != r {
11            let v = v - usize::from(v > r);
12            a[v][v] += M::from(1);
13            if u != r {
14                a[u - usize::from(u > r)][v] -= M::from(1);
15            }
16        }
17    }
18    pp!(a.determinant());
19}

crates/library_checker/src/graph/counting_spanning_tree_undirected.rs (line 8)

5pub fn counting_spanning_tree_undirected(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, m, edges: [(usize, usize); iter m]);
8    let mut a = Matrix::<AddMulOperation<M>>::zeros((n - 1, n - 1));
9    for (u, v) in edges {
10        if u < n - 1 {
11            a[u][u] += M::from(1);
12        }
13        if v < n - 1 {
14            a[v][v] += M::from(1);
15        }
16        if u < n - 1 && v < n - 1 {
17            a[u][v] -= M::from(1);
18            a[v][u] -= M::from(1);
19        }
20    }
21    pp!(a.determinant());
22}

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

473    pub fn new(automaton: A) -> Self {
474        let sigma = automaton.sigma();
475        Self {
476            automaton,
477            prefixes: vec![],
478            suffixes: vec![],
479            inv_h: Matrix::zeros((0, 0)),
480            nh: vec![Matrix::zeros((0, 0)); sigma],
481            wfa: WeightedFiniteAutomaton {
482                initial_weights: Matrix::zeros((1, 0)),
483                transitions: vec![Matrix::zeros((0, 0)); sigma],
484                final_weights: Matrix::zeros((0, 1)),
485            },
486            _marker: PhantomData,
487        }
488    }

crates/library_checker/src/graph/counting_eulerian_circuits.rs (line 12)

9pub fn counting_eulerian_circuits(reader: impl Read, writer: impl Write) {
10    prepare_io!(reader, writer);
11    sc!(n, m, edges: [(usize, usize); iter m]);
12    let mut a = Matrix::<AddMulOperation<M>>::zeros((n, n));
13    let mut indegree = vec![0; n];
14    let mut outdegree = vec![0; n];
15    for (u, v) in edges {
16        a[u][v] -= M::from(1);
17        a[v][v] += M::from(1);
18        outdegree[u] += 1;
19        indegree[v] += 1;
20    }
21    if indegree != outdegree {
22        pp!(0);
23        return;
24    }
25    let root = outdegree.iter().position(|&d| d != 0).unwrap();
26    for i in 0..n {
27        a[root][i] = M::from(0);
28        a[i][root] = M::from(0);
29        if outdegree[i] == 0 {
30            a[i][i] = M::one();
31        }
32    }
33    a[root][root] = M::one();
34    let factorial = MemorizedFactorial::new(*outdegree.iter().max().unwrap() - 1);
35    let mut ans = a.determinant();
36    for d in outdegree {
37        if d != 0 {
38            ans *= factorial.fact[d - 1];
39        }
40    }
41    pp!(ans);
42}

Additional examples can be found in:

Source

pub fn eye(shape: (usize, usize)) -> Self

Examples found in repository ?

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

380    pub fn inverse(&self) -> Option<Matrix<R>> {
381        assert_eq!(self.shape.0, self.shape.1);
382        let n = self.shape.0;
383        if n >= 64 {
384            let m = n / 2;
385            let a = Self::new_with((m, m), |i, j| self[i][j].clone());
386            if let Some(mut ai) = a.inverse() {
387                let b = Self::new_with((m, n - m), |i, j| self[i][j + m].clone());
388                let c = Self::new_with((n - m, m), |i, j| self[i + m][j].clone());
389                let mut d = Self::new_with((n - m, n - m), |i, j| self[i + m][j + m].clone());
390                let u = &ai * &b;
391                let v = &c * &ai;
392                d -= &v * &b;
393                let di = d.inverse()?;
394                let r = &u * &di;
395                let t = &di * &v;
396                ai += &r * &v;
397                let mut inverse = Self::zeros((n, n));
398                for i in 0..m {
399                    inverse[i][..m].clone_from_slice(&ai[i]);
400                    for (x, y) in inverse[i][m..].iter_mut().zip(&r[i]) {
401                        *x = R::neg(y);
402                    }
403                }
404                for i in m..n {
405                    for (x, y) in inverse[i][..m].iter_mut().zip(&t[i - m]) {
406                        *x = R::neg(y);
407                    }
408                    inverse[i][m..].clone_from_slice(&di[i - m]);
409                }
410                return Some(inverse);
411            }
412        }
413        let mut a = self.clone();
414        let mut inverse = Self::eye((n, n));
415        let mut ranges: Vec<_> = (0..n).map(|i| (i, i + 1)).collect();
416        for r in 0..n {
417            let pivot = (r..n).find(|&i| !R::is_zero(&a[i][r]))?;
418            a.data.swap(r, pivot);
419            inverse.data.swap(r, pivot);
420            ranges.swap(r, pivot);
421
422            let d = R::inv(&a[r][r]);
423            for x in &mut a[r][r..] {
424                R::mul_assign(x, &d);
425            }
426            let (left, right) = ranges[r];
427            for x in &mut inverse[r][left..right] {
428                R::mul_assign(x, &d);
429            }
430
431            let (a_upper, a_lower) = a.data.split_at_mut(r + 1);
432            let pivot_a = &a_upper[r];
433            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r + 1);
434            let pivot_inverse = &inverse_upper[r];
435            let (ranges_upper, ranges_lower) = ranges.split_at_mut(r + 1);
436            let (left, right) = ranges_upper[r];
437            for ((a, inverse), range) in a_lower.iter_mut().zip(inverse_lower).zip(ranges_lower) {
438                if R::is_zero(&a[r]) {
439                    continue;
440                }
441                let e = a[r].clone();
442                a[r] = R::zero();
443                R::add_scaled_assign(&mut a[(r + 1)..], &pivot_a[(r + 1)..], &R::neg(&e));
444                R::add_scaled_assign(
445                    &mut inverse[left..right],
446                    &pivot_inverse[left..right],
447                    &R::neg(&e),
448                );
449                range.0 = range.0.min(left);
450                range.1 = range.1.max(right);
451            }
452        }
453        for r in (0..n).rev() {
454            let (left, right) = ranges[r];
455            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r);
456            let pivot_inverse = &inverse_lower[0];
457            let (ranges_upper, _) = ranges.split_at_mut(r);
458            for ((a, inverse), range) in a.data[..r].iter_mut().zip(inverse_upper).zip(ranges_upper)
459            {
460                if R::is_zero(&a[r]) {
461                    continue;
462                }
463                let e = a[r].clone();
464                a[r] = R::zero();
465                R::add_scaled_assign(
466                    &mut inverse[left..right],
467                    &pivot_inverse[left..right],
468                    &R::neg(&e),
469                );
470                range.0 = range.0.min(left);
471                range.1 = range.1.max(right);
472            }
473        }
474        Some(inverse)
475    }
476
477    pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
478        let n = self.shape.0;
479        if n == 0 {
480            return vec![R::one()];
481        }
482        assert!(self.data.iter().all(|a| a.len() == n));
483        for j in 0..(n - 1) {
484            if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
485                self.data.swap(j + 1, x);
486                self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
487                let inv = R::inv(&self[j + 1][j]);
488                let mut v = vec![];
489                let src = std::mem::take(&mut self[j + 1]);
490                for a in self.data[(j + 2)..].iter_mut() {
491                    let mul = R::mul(&a[j], &inv);
492                    R::add_scaled_assign(&mut a[j..], &src[j..], &R::neg(&mul));
493                    v.push(mul);
494                }
495                self[j + 1] = src;
496                for a in self.data.iter_mut() {
497                    let v = R::dot_product(&a[(j + 2)..], &v);
498                    R::add_assign(&mut a[j + 1], &v);
499                }
500            }
501        }
502        // dp[k][j - k] stores [x^k] det(xI - A[..j, ..j]).
503        let mut dp: Vec<Vec<R::T>> = (0..=n).map(|i| Vec::with_capacity(n + 1 - i)).collect();
504        dp[0].push(R::one());
505        for i in 0..n {
506            let mut c = vec![R::zero(); i + 1];
507            c[i] = R::neg(&self[i][i]);
508            let mut mul = R::one();
509            for j in (0..i).rev() {
510                mul = R::mul(&mul, &self[j + 1][j]);
511                c[j] = R::neg(&R::mul(&mul, &self[j][i]));
512            }
513            for k in (0..=i).rev() {
514                let mut value = R::dot_product(&dp[k], &c[k..]);
515                if k > 0 {
516                    R::add_assign(&mut value, dp[k - 1].last().unwrap());
517                }
518                dp[k].push(value);
519            }
520            dp[i + 1].push(R::one());
521        }
522        dp.into_iter().map(|mut c| c.pop().unwrap()).collect()
523    }
524}
525
526impl<R> Index<usize> for Matrix<R>
527where
528    R: SemiRing,
529{
530    type Output = Vec<R::T>;
531    fn index(&self, index: usize) -> &Self::Output {
532        &self.data[index]
533    }
534}
535
536impl<R> IndexMut<usize> for Matrix<R>
537where
538    R: SemiRing,
539{
540    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
541        &mut self.data[index]
542    }
543}
544
545impl<R> Index<(usize, usize)> for Matrix<R>
546where
547    R: SemiRing,
548{
549    type Output = R::T;
550    fn index(&self, index: (usize, usize)) -> &Self::Output {
551        &self.data[index.0][index.1]
552    }
553}
554
555impl<R> IndexMut<(usize, usize)> for Matrix<R>
556where
557    R: SemiRing,
558{
559    fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
560        &mut self.data[index.0][index.1]
561    }
562}
563
564macro_rules! impl_matrix_pairwise_binop {
565    ($imp:ident, $method:ident, $imp_assign:ident, $method_assign:ident $(where [$($clauses:tt)*])?) => {
566        impl<R> $imp_assign for Matrix<R>
567        where
568            R: SemiRing,
569            $($($clauses)*)?
570        {
571            fn $method_assign(&mut self, rhs: Self) {
572                self.pairwise_assign(&rhs, |a, b| R::$method_assign(a, b));
573            }
574        }
575        impl<R> $imp_assign<&Matrix<R>> for Matrix<R>
576        where
577            R: SemiRing,
578            $($($clauses)*)?
579        {
580            fn $method_assign(&mut self, rhs: &Self) {
581                self.pairwise_assign(rhs, |a, b| R::$method_assign(a, b));
582            }
583        }
584        impl<R> $imp for Matrix<R>
585        where
586            R: SemiRing,
587            $($($clauses)*)?
588        {
589            type Output = Matrix<R>;
590            fn $method(mut self, rhs: Self) -> Self::Output {
591                self.$method_assign(rhs);
592                self
593            }
594        }
595        impl<R> $imp<&Matrix<R>> for Matrix<R>
596        where
597            R: SemiRing,
598            $($($clauses)*)?
599        {
600            type Output = Matrix<R>;
601            fn $method(mut self, rhs: &Self) -> Self::Output {
602                self.$method_assign(rhs);
603                self
604            }
605        }
606        impl<R> $imp<Matrix<R>> for &Matrix<R>
607        where
608            R: SemiRing,
609            $($($clauses)*)?
610        {
611            type Output = Matrix<R>;
612            fn $method(self, mut rhs: Matrix<R>) -> Self::Output {
613                rhs.pairwise_assign(self, |a, b| *a = R::$method(b, a));
614                rhs
615            }
616        }
617        impl<R> $imp<&Matrix<R>> for &Matrix<R>
618        where
619            R: SemiRing,
620            $($($clauses)*)?
621        {
622            type Output = Matrix<R>;
623            fn $method(self, rhs: &Matrix<R>) -> Self::Output {
624                let mut this = self.clone();
625                this.$method_assign(rhs);
626                this
627            }
628        }
629    };
630}
631
632impl_matrix_pairwise_binop!(Add, add, AddAssign, add_assign);
633impl_matrix_pairwise_binop!(Sub, sub, SubAssign, sub_assign where [R: SemiRing<Additive: Invertible>]);
634
635impl<R> Mul for Matrix<R>
636where
637    R: SemiRing,
638{
639    type Output = Matrix<R>;
640    fn mul(self, rhs: Self) -> Self::Output {
641        (&self).mul(&rhs)
642    }
643}
644impl<R> Mul<&Matrix<R>> for Matrix<R>
645where
646    R: SemiRing,
647{
648    type Output = Matrix<R>;
649    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
650        (&self).mul(rhs)
651    }
652}
653impl<R> Mul<Matrix<R>> for &Matrix<R>
654where
655    R: SemiRing,
656{
657    type Output = Matrix<R>;
658    fn mul(self, rhs: Matrix<R>) -> Self::Output {
659        self.mul(&rhs)
660    }
661}
662impl<R> Mul<&Matrix<R>> for &Matrix<R>
663where
664    R: SemiRing,
665{
666    type Output = Matrix<R>;
667    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
668        assert_eq!(self.shape.1, rhs.shape.0);
669        if let Some(data) = R::try_matrix_product(&self.data, &rhs.data) {
670            return Matrix::from_vec(data);
671        }
672        let rhs = rhs.transpose();
673        Matrix::new_with((self.shape.0, rhs.shape.0), |i, j| {
674            R::dot_product(&self[i], &rhs[j])
675        })
676    }
677}
678
679fn strassen_rec<R: Ring>(
680    a: &[R::T],
681    b: &[R::T],
682    c: &mut [R::T],
683    shape: (usize, usize, usize),
684    stride_a: usize,
685    stride_b: usize,
686) {
687    let (n, m, p) = shape;
688    fn add_block<R: Ring>(
689        a: &[R::T],
690        b: &[R::T],
691        out: &mut [R::T],
692        n: usize,
693        stride_a: usize,
694        stride_b: usize,
695    ) {
696        for ((a, b), c) in a
697            .chunks(stride_a)
698            .zip(b.chunks(stride_b))
699            .zip(out.chunks_exact_mut(n))
700        {
701            for ((a, b), c) in a.iter().zip(b.iter()).zip(c.iter_mut()) {
702                *c = R::add(a, b);
703            }
704        }
705    }
706
707    fn sub_block<R: Ring>(
708        a: &[R::T],
709        b: &[R::T],
710        out: &mut [R::T],
711        n: usize,
712        stride_a: usize,
713        stride_b: usize,
714    ) {
715        for ((a, b), c) in a
716            .chunks(stride_a)
717            .zip(b.chunks(stride_b))
718            .zip(out.chunks_exact_mut(n))
719        {
720            for ((a, b), c) in a.iter().zip(b.iter()).zip(c.iter_mut()) {
721                *c = R::sub(a, b);
722            }
723        }
724    }
725
726    if n.min(m).min(p) <= 128 {
727        let transposed: Vec<_> = (0..p)
728            .flat_map(|j| (0..m).map(move |i| b[i * stride_b + j].clone()))
729            .collect();
730        for (a, c) in a.chunks(stride_a).zip(c.chunks_exact_mut(p)) {
731            for (b, c) in transposed.chunks_exact(m).zip(c) {
732                *c = R::dot_product(&a[..m], b);
733            }
734        }
735        return;
736    }
737    let (h, k, w) = (n / 2, m / 2, p / 2);
738    let a11 = 0;
739    let a12 = k;
740    let a21 = h * stride_a;
741    let a22 = a21 + k;
742    let b11 = 0;
743    let b12 = w;
744    let b21 = k * stride_b;
745    let b22 = b21 + w;
746
747    let block = h * w;
748    let mut buf = vec![R::zero(); h * k + k * w + block * 7];
749    let (s1, rest) = buf.split_at_mut(h * k);
750    let (s2, m_buf) = rest.split_at_mut(k * w);
751    let (m1, rest) = m_buf.split_at_mut(block);
752    let (m2, rest) = rest.split_at_mut(block);
753    let (m3, rest) = rest.split_at_mut(block);
754    let (m4, rest) = rest.split_at_mut(block);
755    let (m5, rest) = rest.split_at_mut(block);
756    let (m6, m7) = rest.split_at_mut(block);
757
758    // (A11 + A22)(B11 + B22)
759    add_block::<R>(&a[a11..], &a[a22..], s1, k, stride_a, stride_a);
760    add_block::<R>(&b[b11..], &b[b22..], s2, w, stride_b, stride_b);
761    strassen_rec::<R>(s1, s2, m1, (h, k, w), k, w);
762
763    // (A21 + A22) B11
764    add_block::<R>(&a[a21..], &a[a22..], s1, k, stride_a, stride_a);
765    strassen_rec::<R>(s1, &b[b11..], m2, (h, k, w), k, stride_b);
766
767    // A11 (B12 - B22)
768    sub_block::<R>(&b[b12..], &b[b22..], s2, w, stride_b, stride_b);
769    strassen_rec::<R>(&a[a11..], s2, m3, (h, k, w), stride_a, w);
770
771    // A22 (B21 - B11)
772    sub_block::<R>(&b[b21..], &b[b11..], s2, w, stride_b, stride_b);
773    strassen_rec::<R>(&a[a22..], s2, m4, (h, k, w), stride_a, w);
774
775    // (A11 + A12) B22
776    add_block::<R>(&a[a11..], &a[a12..], s1, k, stride_a, stride_a);
777    strassen_rec::<R>(s1, &b[b22..], m5, (h, k, w), k, stride_b);
778
779    // (A21 - A11)(B11 + B12)
780    sub_block::<R>(&a[a21..], &a[a11..], s1, k, stride_a, stride_a);
781    add_block::<R>(&b[b11..], &b[b12..], s2, w, stride_b, stride_b);
782    strassen_rec::<R>(s1, s2, m6, (h, k, w), k, w);
783
784    // (A12 - A22)(B21 + B22)
785    sub_block::<R>(&a[a12..], &a[a22..], s1, k, stride_a, stride_a);
786    add_block::<R>(&b[b21..], &b[b22..], s2, w, stride_b, stride_b);
787    strassen_rec::<R>(s1, s2, m7, (h, k, w), k, w);
788
789    let c11 = 0;
790    let c12 = w;
791    let c21 = h * p;
792    let c22 = c21 + w;
793    for ((((m1, m4), m5), m7), c) in m1
794        .iter()
795        .zip(m4.iter())
796        .zip(m5.iter())
797        .zip(m7.iter())
798        .zip(c[c11..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
799    {
800        *c = R::add(m1, m4);
801        R::sub_assign(c, m5);
802        R::add_assign(c, m7);
803    }
804    for ((m3, m5), c) in m3
805        .iter()
806        .zip(m5.iter())
807        .zip(c[c12..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
808    {
809        *c = R::add(m3, m5);
810    }
811    for ((m2, m4), c) in m2
812        .iter()
813        .zip(m4.iter())
814        .zip(c[c21..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
815    {
816        *c = R::add(m2, m4);
817    }
818    for ((((m1, m2), m3), m6), c) in m1
819        .iter()
820        .zip(m2.iter())
821        .zip(m3.iter())
822        .zip(m6.iter())
823        .zip(c[c22..].chunks_mut(p).flat_map(|c| c.iter_mut().take(w)))
824    {
825        *c = R::sub(m1, m2);
826        R::add_assign(c, m3);
827        R::add_assign(c, m6);
828    }
829}
830
831impl<R> Matrix<R>
832where
833    R: Ring,
834{
835    pub fn mul_strassen(&self, rhs: &Matrix<R>) -> Matrix<R> {
836        assert_eq!(self.shape.1, rhs.shape.0);
837        if let Some(data) = R::try_matrix_product(&self.data, &rhs.data) {
838            return Matrix::from_vec(data);
839        }
840        let (n, m) = self.shape;
841        let p = rhs.shape.1;
842        if n == 0 || m == 0 || p == 0 {
843            return Matrix::zeros((n, p));
844        }
845        let split = n.min(m).min(p).div_ceil(128).next_power_of_two();
846        if split <= 2 {
847            return self * rhs;
848        }
849        let rows = n.div_ceil(split) * split;
850        let inner = m.div_ceil(split) * split;
851        let cols = p.div_ceil(split) * split;
852        let mut a = vec![R::zero(); rows * inner];
853        for (a, data) in a.chunks_exact_mut(inner).zip(&self.data) {
854            a[..m].clone_from_slice(data);
855        }
856        let mut b = vec![R::zero(); inner * cols];
857        for (b, data) in b.chunks_exact_mut(cols).zip(&rhs.data) {
858            b[..p].clone_from_slice(data);
859        }
860        let mut c = vec![R::zero(); rows * cols];
861        strassen_rec::<R>(&a, &b, &mut c, (rows, inner, cols), inner, cols);
862        let mut res = Matrix::zeros((n, p));
863        for (data, c) in res.data.iter_mut().zip(c.chunks_exact(cols)) {
864            data.clone_from_slice(&c[..p]);
865        }
866        res
867    }
868}
869
870impl<R> MulAssign<&R::T> for Matrix<R>
871where
872    R: SemiRing,
873{
874    fn mul_assign(&mut self, rhs: &R::T) {
875        for i in 0..self.shape.0 {
876            for j in 0..self.shape.1 {
877                R::mul_assign(&mut self[(i, j)], rhs);
878            }
879        }
880    }
881}
882
883impl<R> Neg for Matrix<R>
884where
885    R: SemiRing<Additive: Invertible>,
886{
887    type Output = Self;
888
889    fn neg(self) -> Self::Output {
890        self.map(|x| R::neg(x))
891    }
892}
893
894impl<R> Neg for &Matrix<R>
895where
896    R: SemiRing<Additive: Invertible>,
897{
898    type Output = Matrix<R>;
899
900    fn neg(self) -> Self::Output {
901        self.map(|x| R::neg(x))
902    }
903}
904
905impl<R> Matrix<R>
906where
907    R: SemiRing,
908{
909    pub fn pow(self, mut n: usize) -> Self {
910        assert_eq!(self.shape.0, self.shape.1);
911        let mut res = Matrix::eye(self.shape);
912        let mut x = self;
913        while n > 0 {
914            if n & 1 == 1 {
915                res = &res * &x;
916            }
917            x = &x * &x;
918            n >>= 1;
919        }
920        res
921    }
922}
923
924impl<R> Matrix<R>
925where
926    R: Ring,
927{
928    pub fn pow_strassen(self, mut n: usize) -> Self {
929        assert_eq!(self.shape.0, self.shape.1);
930        let mut res = Matrix::eye(self.shape);
931        let mut x = self;
932        while n > 0 {
933            if n & 1 == 1 {
934                res = res.mul_strassen(&x);
935            }
936            x = x.mul_strassen(&x);
937            n >>= 1;
938        }
939        res
940    }

Source

pub fn transpose(&self) -> Self

Examples found in repository ?

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

667    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
668        assert_eq!(self.shape.1, rhs.shape.0);
669        if let Some(data) = R::try_matrix_product(&self.data, &rhs.data) {
670            return Matrix::from_vec(data);
671        }
672        let rhs = rhs.transpose();
673        Matrix::new_with((self.shape.0, rhs.shape.0), |i, j| {
674            R::dot_product(&self[i], &rhs[j])
675        })
676    }

More examples

Hide additional examples

crates/competitive/src/math/mint_matrix.rs (line 46)

41    fn pow_frobenius(self, k: usize) -> Self
42    where
43        M: MIntConvert<u64>,
44    {
45        assert_eq!(self.shape.0, self.shape.1);
46        let a = self.transpose();
47        let mut rng = Xorshift::new();
48        let f = loop {
49            if let Some(f) = frobenius_decomposition(&a, &mut rng) {
50                break f;
51            }
52        };
53        let fk = f.pow(k);
54        let n = f.t.shape.0;
55        if f.blocks
56            .iter()
57            .map(|p| (p.0.len() - 1).pow(2))
58            .sum::<usize>()
59            * 4
60            <= n * n
61        {
62            let mut ft = Matrix::zeros((n, n));
63            let mut first = 0;
64            for p in &f.blocks {
65                let d = p.0.len() - 1;
66                for i in first..first + d {
67                    for j in first..first + d {
68                        MInt::add_scaled_assign(&mut ft[i], &f.t[j], &fk[i][j]);
69                    }
70                }
71                first += d;
72            }
73            &f.t_inv * &ft
74        } else {
75            &(&f.t_inv * &fk) * &f.t
76        }
77    }

Source

pub fn map<S, F>(&self, f: F) -> Matrix<S>
where S: SemiRing, F: FnMut(&R::T) -> S::T,

Examples found in repository ?

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

889    fn neg(self) -> Self::Output {
890        self.map(|x| R::neg(x))
891    }
892}
893
894impl<R> Neg for &Matrix<R>
895where
896    R: SemiRing<Additive: Invertible>,
897{
898    type Output = Matrix<R>;
899
900    fn neg(self) -> Self::Output {
901        self.map(|x| R::neg(x))
902    }

Source

pub fn add_row_with(&mut self, f: impl FnMut(usize, usize) -> R::T)

Examples found in repository ?

crates/competitive/src/algorithm/automata_learning.rs (lines 563-569)

523    pub fn train_sample(&mut self, sample: &[usize]) -> bool {
524        let Some((prefix, suffix)) = self.split_sample(sample) else {
525            return false;
526        };
527        self.prefixes.push(prefix);
528        self.suffixes.push(suffix);
529        let n = self.inv_h.shape.0;
530        let prefix = &self.prefixes[n];
531        let suffix = &self.suffixes[n];
532        let u = Matrix::<F>::new_with((n, 1), |i, _| {
533            self.automaton.behavior(
534                self.prefixes[i]
535                    .iter()
536                    .cloned()
537                    .chain(suffix.iter().cloned()),
538            )
539        });
540        let v = Matrix::<F>::new_with((1, n), |_, j| {
541            self.automaton.behavior(
542                prefix
543                    .iter()
544                    .cloned()
545                    .chain(self.suffixes[j].iter().cloned()),
546            )
547        });
548        let w = Matrix::<F>::new_with((1, 1), |_, _| {
549            self.automaton
550                .behavior(prefix.iter().cloned().chain(suffix.iter().cloned()))
551        });
552        let t = &self.inv_h * &u;
553        let s = &v * &self.inv_h;
554        let d = F::inv(&(&w - &(&v * &t))[0][0]);
555        let dh = &t * &s;
556        for i in 0..n {
557            for j in 0..n {
558                F::add_assign(&mut self.inv_h[i][j], &F::mul(&dh[i][j], &d));
559            }
560        }
561        self.inv_h
562            .add_col_with(|i, _| F::neg(&F::mul(&t[i][0], &d)));
563        self.inv_h.add_row_with(|_, j| {
564            if j != n {
565                F::neg(&F::mul(&s[0][j], &d))
566            } else {
567                d.clone()
568            }
569        });
570
571        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
572            let b = &(&self.nh[x] * &t) * &s;
573            for i in 0..n {
574                for j in 0..n {
575                    F::add_assign(&mut transition[i][j], &F::mul(&b[i][j], &d));
576                }
577            }
578        }
579        for (x, nh) in self.nh.iter_mut().enumerate() {
580            nh.add_col_with(|i, j| {
581                self.automaton.behavior(
582                    self.prefixes[i]
583                        .iter()
584                        .cloned()
585                        .chain([x])
586                        .chain(self.suffixes[j].iter().cloned()),
587                )
588            });
589            nh.add_row_with(|i, j| {
590                self.automaton.behavior(
591                    self.prefixes[i]
592                        .iter()
593                        .cloned()
594                        .chain([x])
595                        .chain(self.suffixes[j].iter().cloned()),
596                )
597            });
598        }
599        self.wfa
600            .initial_weights
601            .add_col_with(|_, _| if n == 0 { F::one() } else { F::zero() });
602        self.wfa
603            .final_weights
604            .add_row_with(|_, _| self.automaton.behavior(prefix.iter().cloned()));
605        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
606            transition.add_col_with(|_, _| F::zero());
607            transition.add_row_with(|_, _| F::zero());
608            for i in 0..=n {
609                for j in 0..=n {
610                    if i == n || j == n {
611                        for k in 0..=n {
612                            if i != n && j != n && k != n {
613                                continue;
614                            }
615                            F::add_assign(
616                                &mut transition[i][k],
617                                &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
618                            );
619                        }
620                    } else {
621                        let k = n;
622                        F::add_assign(
623                            &mut transition[i][k],
624                            &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
625                        );
626                    }
627                }
628            }
629        }
630        true
631    }

Source

pub fn add_col_with(&mut self, f: impl FnMut(usize, usize) -> R::T)

Examples found in repository ?

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

523    pub fn train_sample(&mut self, sample: &[usize]) -> bool {
524        let Some((prefix, suffix)) = self.split_sample(sample) else {
525            return false;
526        };
527        self.prefixes.push(prefix);
528        self.suffixes.push(suffix);
529        let n = self.inv_h.shape.0;
530        let prefix = &self.prefixes[n];
531        let suffix = &self.suffixes[n];
532        let u = Matrix::<F>::new_with((n, 1), |i, _| {
533            self.automaton.behavior(
534                self.prefixes[i]
535                    .iter()
536                    .cloned()
537                    .chain(suffix.iter().cloned()),
538            )
539        });
540        let v = Matrix::<F>::new_with((1, n), |_, j| {
541            self.automaton.behavior(
542                prefix
543                    .iter()
544                    .cloned()
545                    .chain(self.suffixes[j].iter().cloned()),
546            )
547        });
548        let w = Matrix::<F>::new_with((1, 1), |_, _| {
549            self.automaton
550                .behavior(prefix.iter().cloned().chain(suffix.iter().cloned()))
551        });
552        let t = &self.inv_h * &u;
553        let s = &v * &self.inv_h;
554        let d = F::inv(&(&w - &(&v * &t))[0][0]);
555        let dh = &t * &s;
556        for i in 0..n {
557            for j in 0..n {
558                F::add_assign(&mut self.inv_h[i][j], &F::mul(&dh[i][j], &d));
559            }
560        }
561        self.inv_h
562            .add_col_with(|i, _| F::neg(&F::mul(&t[i][0], &d)));
563        self.inv_h.add_row_with(|_, j| {
564            if j != n {
565                F::neg(&F::mul(&s[0][j], &d))
566            } else {
567                d.clone()
568            }
569        });
570
571        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
572            let b = &(&self.nh[x] * &t) * &s;
573            for i in 0..n {
574                for j in 0..n {
575                    F::add_assign(&mut transition[i][j], &F::mul(&b[i][j], &d));
576                }
577            }
578        }
579        for (x, nh) in self.nh.iter_mut().enumerate() {
580            nh.add_col_with(|i, j| {
581                self.automaton.behavior(
582                    self.prefixes[i]
583                        .iter()
584                        .cloned()
585                        .chain([x])
586                        .chain(self.suffixes[j].iter().cloned()),
587                )
588            });
589            nh.add_row_with(|i, j| {
590                self.automaton.behavior(
591                    self.prefixes[i]
592                        .iter()
593                        .cloned()
594                        .chain([x])
595                        .chain(self.suffixes[j].iter().cloned()),
596                )
597            });
598        }
599        self.wfa
600            .initial_weights
601            .add_col_with(|_, _| if n == 0 { F::one() } else { F::zero() });
602        self.wfa
603            .final_weights
604            .add_row_with(|_, _| self.automaton.behavior(prefix.iter().cloned()));
605        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
606            transition.add_col_with(|_, _| F::zero());
607            transition.add_row_with(|_, _| F::zero());
608            for i in 0..=n {
609                for j in 0..=n {
610                    if i == n || j == n {
611                        for k in 0..=n {
612                            if i != n && j != n && k != n {
613                                continue;
614                            }
615                            F::add_assign(
616                                &mut transition[i][k],
617                                &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
618                            );
619                        }
620                    } else {
621                        let k = n;
622                        F::add_assign(
623                            &mut transition[i][k],
624                            &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
625                        );
626                    }
627                }
628            }
629        }
630        true
631    }

Source

pub fn pairwise_assign<F>(&mut self, other: &Self, f: F)
where F: FnMut(&mut R::T, &R::T),

Source §

impl<R> Matrix<R>
where R: Field<T: PartialEq, Additive: Invertible, Multiplicative: Invertible>,

Source

fn eliminate<const DETERMINANT: bool>(&mut self) -> (usize, R::T)

Examples found in repository ?

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

297    pub fn rank(&mut self) -> usize {
298        self.eliminate::<false>().0
299    }
300
301    pub fn determinant(&mut self) -> R::T {
302        assert_eq!(self.shape.0, self.shape.1);
303        self.eliminate::<true>().1
304    }
305
306    pub fn solve_system_of_linear_equations(
307        &self,
308        b: &[R::T],
309    ) -> Option<SystemOfLinearEquationsSolution<R>> {
310        assert_eq!(self.shape.0, b.len());
311        let m = self.shape.1;
312        let mut a = Self::new_with((self.shape.0, m + 1), |i, j| {
313            if j == m {
314                b[i].clone()
315            } else {
316                self[i][j].clone()
317            }
318        });
319        let rank = a.eliminate::<false>().0;
320        let mut pivots = Vec::with_capacity(rank);
321        let mut b = Vec::with_capacity(rank);
322        for row in &a.data[..rank] {
323            let c = row.iter().position(|x| !R::is_zero(x)).unwrap();
324            if c == m {
325                return None;
326            }
327            pivots.push(c);
328            b.push(row[m].clone());
329        }
330
331        let mut free = Vec::with_capacity(m - rank);
332        let mut pivot = 0;
333        for c in 0..m {
334            if pivot < rank && pivots[pivot] == c {
335                pivot += 1;
336            } else {
337                free.push(c);
338            }
339        }
340        let mut coefficients: Vec<Vec<_>> = (0..rank)
341            .map(|i| free.iter().map(|&c| a[i][c].clone()).collect())
342            .collect();
343        for k in (0..rank).rev() {
344            let c = pivots[k];
345            let inv = R::inv(&a[k][c]);
346            R::mul_assign(&mut b[k], &inv);
347            let pivot_b = b[k].clone();
348            let (upper, lower) = coefficients.split_at_mut(k);
349            let pivot_coefficients = &mut lower[0];
350            for x in pivot_coefficients.iter_mut() {
351                R::mul_assign(x, &inv);
352            }
353            for ((row, value), coefficients) in a.data[..k].iter_mut().zip(&mut b[..k]).zip(upper) {
354                if R::is_zero(&row[c]) {
355                    continue;
356                }
357                let factor = row[c].clone();
358                row[c] = R::zero();
359                R::sub_assign(value, &R::mul(&factor, &pivot_b));
360                R::add_scaled_assign(coefficients, pivot_coefficients, &R::neg(&factor));
361            }
362        }
363
364        let mut particular = vec![R::zero(); m];
365        for i in 0..rank {
366            particular[pivots[i]] = b[i].clone();
367        }
368        let mut basis = Vec::with_capacity(free.len());
369        for (j, &c) in free.iter().enumerate() {
370            let mut vector = vec![R::zero(); m];
371            vector[c] = R::one();
372            for i in 0..rank {
373                vector[pivots[i]] = R::neg(&coefficients[i][j]);
374            }
375            basis.push(vector);
376        }
377        Some(SystemOfLinearEquationsSolution { particular, basis })
378    }

Source

pub fn row_reduction_with<F>(&mut self, normalize: bool, f: F)
where F: FnMut(usize, usize, usize),

f: (row, pivot_row, col)

Examples found in repository ?

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

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

More examples

Hide additional examples

crates/competitive/src/algorithm/automata_learning.rs (lines 669-672)

637    pub fn batch_train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
638        let mut prefix_set: HashSet<_> = self.prefixes.iter().cloned().collect();
639        let mut suffix_set: HashSet<_> = self.suffixes.iter().cloned().collect();
640        for sample in samples {
641            if prefix_set.insert(sample.to_vec()) {
642                self.prefixes.push(sample.to_vec());
643            }
644            if suffix_set.insert(sample.to_vec()) {
645                self.suffixes.push(sample);
646            }
647        }
648        let mut h = Matrix::<F>::new_with((self.prefixes.len(), self.suffixes.len()), |i, j| {
649            self.automaton.behavior(
650                self.prefixes[i]
651                    .iter()
652                    .cloned()
653                    .chain(self.suffixes[j].iter().cloned()),
654            )
655        });
656        if !self.prefixes.is_empty() && !self.suffixes.is_empty() && F::is_zero(&h[0][0]) {
657            for j in 1..self.suffixes.len() {
658                if !F::is_zero(&h[0][j]) {
659                    self.suffixes.swap(0, j);
660                    for row in &mut h.data {
661                        row.swap(0, j);
662                    }
663                    break;
664                }
665            }
666        }
667        let mut row_id: Vec<usize> = (0..h.shape.0).collect();
668        let mut pivots = vec![];
669        h.row_reduction_with(false, |r, p, c| {
670            row_id.swap(r, p);
671            pivots.push((row_id[r], c));
672        });
673        let mut new_prefixes = vec![];
674        let mut new_suffixes = vec![];
675        for (i, j) in pivots {
676            new_prefixes.push(self.prefixes[i].clone());
677            new_suffixes.push(self.suffixes[j].clone());
678        }
679        self.prefixes = new_prefixes;
680        self.suffixes = new_suffixes;
681        assert_eq!(self.prefixes.len(), self.suffixes.len());
682        let n = self.prefixes.len();
683        let h = Matrix::<F>::new_with((n, n), |i, j| {
684            self.automaton.behavior(
685                self.prefixes[i]
686                    .iter()
687                    .cloned()
688                    .chain(self.suffixes[j].iter().cloned()),
689            )
690        });
691        self.inv_h = h.inverse().expect("Hankel matrix must be invertible");
692        self.wfa = WeightedFiniteAutomaton::<F> {
693            initial_weights: Matrix::new_with((1, n), |_, j| {
694                if self.prefixes[j].is_empty() {
695                    F::one()
696                } else {
697                    F::zero()
698                }
699            }),
700            transitions: (0..self.automaton.sigma())
701                .map(|x| {
702                    &Matrix::new_with((n, n), |i, j| {
703                        self.automaton.behavior(
704                            self.prefixes[i]
705                                .iter()
706                                .cloned()
707                                .chain([x])
708                                .chain(self.suffixes[j].iter().cloned()),
709                        )
710                    }) * &self.inv_h
711                })
712                .collect(),
713            final_weights: Matrix::new_with((n, 1), |i, _| {
714                self.automaton.behavior(self.prefixes[i].iter().cloned())
715            }),
716        };
717    }

Source

pub fn row_reduction(&mut self, normalize: bool)

Source

pub fn rank(&mut self) -> usize

Examples found in repository ?

crates/library_checker/src/linear_algebra/matrix_rank.rs (line 21)

5pub fn matrix_rank(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, m);
8    let mut a = if n <= m {
9        sc!(a: [[M; m]; n]);
10        Matrix::<AddMulOperation<_>>::from_vec(a)
11    } else {
12        let mut a = Matrix::<AddMulOperation<_>>::zeros((m, n));
13        for j in 0..n {
14            for row in &mut a.data {
15                sc!(x: M);
16                row[j] = x;
17            }
18        }
19        a
20    };
21    let rank = a.rank();
22    pp!(rank);
23}

Source

pub fn determinant(&mut self) -> R::T

Examples found in repository ?

crates/library_checker/src/linear_algebra/matrix_det.rs (line 9)

5pub fn matrix_det(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, a: [[M; n]; n]);
8    let mut a = Matrix::<AddMulOperation<_>>::from_vec(a);
9    let det = a.determinant();
10    pp!(det);
11}

More examples

Hide additional examples

crates/library_checker/src/graph/counting_spanning_tree_directed.rs (line 18)

5pub fn counting_spanning_tree_directed(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, m, r: usize, edges: [(usize, usize); iter m]);
8    let mut a = Matrix::<AddMulOperation<M>>::zeros((n - 1, n - 1));
9    for (u, v) in edges {
10        if v != r {
11            let v = v - usize::from(v > r);
12            a[v][v] += M::from(1);
13            if u != r {
14                a[u - usize::from(u > r)][v] -= M::from(1);
15            }
16        }
17    }
18    pp!(a.determinant());
19}

crates/library_checker/src/graph/counting_spanning_tree_undirected.rs (line 21)

5pub fn counting_spanning_tree_undirected(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, m, edges: [(usize, usize); iter m]);
8    let mut a = Matrix::<AddMulOperation<M>>::zeros((n - 1, n - 1));
9    for (u, v) in edges {
10        if u < n - 1 {
11            a[u][u] += M::from(1);
12        }
13        if v < n - 1 {
14            a[v][v] += M::from(1);
15        }
16        if u < n - 1 && v < n - 1 {
17            a[u][v] -= M::from(1);
18            a[v][u] -= M::from(1);
19        }
20    }
21    pp!(a.determinant());
22}

crates/library_checker/src/graph/counting_eulerian_circuits.rs (line 35)

9pub fn counting_eulerian_circuits(reader: impl Read, writer: impl Write) {
10    prepare_io!(reader, writer);
11    sc!(n, m, edges: [(usize, usize); iter m]);
12    let mut a = Matrix::<AddMulOperation<M>>::zeros((n, n));
13    let mut indegree = vec![0; n];
14    let mut outdegree = vec![0; n];
15    for (u, v) in edges {
16        a[u][v] -= M::from(1);
17        a[v][v] += M::from(1);
18        outdegree[u] += 1;
19        indegree[v] += 1;
20    }
21    if indegree != outdegree {
22        pp!(0);
23        return;
24    }
25    let root = outdegree.iter().position(|&d| d != 0).unwrap();
26    for i in 0..n {
27        a[root][i] = M::from(0);
28        a[i][root] = M::from(0);
29        if outdegree[i] == 0 {
30            a[i][i] = M::one();
31        }
32    }
33    a[root][root] = M::one();
34    let factorial = MemorizedFactorial::new(*outdegree.iter().max().unwrap() - 1);
35    let mut ans = a.determinant();
36    for d in outdegree {
37        if d != 0 {
38            ans *= factorial.fact[d - 1];
39        }
40    }
41    pp!(ans);
42}

Source

pub fn solve_system_of_linear_equations( &self, b: &[R::T], ) -> Option<SystemOfLinearEquationsSolution<R>>

Examples found in repository ?

crates/library_checker/src/linear_algebra/system_of_linear_equations.rs (line 9)

5pub fn system_of_linear_equations(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, m, a: [[M; m]; n], b: [M; n]);
8    let a = Matrix::<AddMulOperation<M>>::from_vec(a);
9    if let Some(sol) = a.solve_system_of_linear_equations(&b) {
10        pp!(sol.basis.len(); @it2d std::iter::once(sol.particular).chain(sol.basis));
11    } else {
12        pp!(-1);
13    }
14}

More examples

Hide additional examples

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

81    pub fn solve(self) -> EsperSolver<R, Input, Class, FC, FF> {
82        let data: HashMap<_, _> = self
83            .data
84            .into_iter()
85            .map(|(key, SystemOfLinearEquation { a, b })| {
86                (
87                    key,
88                    Matrix::<R>::from_vec(a)
89                        .solve_system_of_linear_equations(&b)
90                        .map(|sol| sol.particular),
91                )
92            })
93            .collect();
94        EsperSolver {
95            class: self.class,
96            feature: self.feature,
97            data,
98            _marker: PhantomData,
99        }
100    }
101
102    pub fn solve_checked(self) -> EsperSolver<R, Input, Class, FC, FF>
103    where
104        Class: Debug,
105        R: Field<T: Debug, Additive: Invertible, Multiplicative: Invertible>,
106    {
107        let data: HashMap<_, _> = self
108            .data
109            .into_iter()
110            .map(|(key, SystemOfLinearEquation { a, b })| {
111                let mat = Matrix::<R>::from_vec(a);
112                let coeff = mat
113                    .solve_system_of_linear_equations(&b)
114                    .map(|sol| sol.particular);
115                if coeff.is_none() {
116                    eprintln!(
117                        "failed to solve linear equations: key={:?} A={:?} b={:?}",
118                        key, mat.data, b
119                    );
120                }
121                (key, coeff)
122            })
123            .collect();
124        EsperSolver {
125            class: self.class,
126            feature: self.feature,
127            data,
128            _marker: PhantomData,
129        }
130    }

Source

pub fn inverse(&self) -> Option<Matrix<R>>

Examples found in repository ?

crates/library_checker/src/linear_algebra/inverse_matrix.rs (line 9)

5pub fn inverse_matrix(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, a: [[M; n]; n]);
8    let a = Matrix::<AddMulOperation<_>>::from_vec(a);
9    if let Some(b) = a.inverse() {
10        pp!(@it2d b.data);
11    } else {
12        pp!("-1");
13    }
14}

More examples

Hide additional examples

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

637    pub fn batch_train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
638        let mut prefix_set: HashSet<_> = self.prefixes.iter().cloned().collect();
639        let mut suffix_set: HashSet<_> = self.suffixes.iter().cloned().collect();
640        for sample in samples {
641            if prefix_set.insert(sample.to_vec()) {
642                self.prefixes.push(sample.to_vec());
643            }
644            if suffix_set.insert(sample.to_vec()) {
645                self.suffixes.push(sample);
646            }
647        }
648        let mut h = Matrix::<F>::new_with((self.prefixes.len(), self.suffixes.len()), |i, j| {
649            self.automaton.behavior(
650                self.prefixes[i]
651                    .iter()
652                    .cloned()
653                    .chain(self.suffixes[j].iter().cloned()),
654            )
655        });
656        if !self.prefixes.is_empty() && !self.suffixes.is_empty() && F::is_zero(&h[0][0]) {
657            for j in 1..self.suffixes.len() {
658                if !F::is_zero(&h[0][j]) {
659                    self.suffixes.swap(0, j);
660                    for row in &mut h.data {
661                        row.swap(0, j);
662                    }
663                    break;
664                }
665            }
666        }
667        let mut row_id: Vec<usize> = (0..h.shape.0).collect();
668        let mut pivots = vec![];
669        h.row_reduction_with(false, |r, p, c| {
670            row_id.swap(r, p);
671            pivots.push((row_id[r], c));
672        });
673        let mut new_prefixes = vec![];
674        let mut new_suffixes = vec![];
675        for (i, j) in pivots {
676            new_prefixes.push(self.prefixes[i].clone());
677            new_suffixes.push(self.suffixes[j].clone());
678        }
679        self.prefixes = new_prefixes;
680        self.suffixes = new_suffixes;
681        assert_eq!(self.prefixes.len(), self.suffixes.len());
682        let n = self.prefixes.len();
683        let h = Matrix::<F>::new_with((n, n), |i, j| {
684            self.automaton.behavior(
685                self.prefixes[i]
686                    .iter()
687                    .cloned()
688                    .chain(self.suffixes[j].iter().cloned()),
689            )
690        });
691        self.inv_h = h.inverse().expect("Hankel matrix must be invertible");
692        self.wfa = WeightedFiniteAutomaton::<F> {
693            initial_weights: Matrix::new_with((1, n), |_, j| {
694                if self.prefixes[j].is_empty() {
695                    F::one()
696                } else {
697                    F::zero()
698                }
699            }),
700            transitions: (0..self.automaton.sigma())
701                .map(|x| {
702                    &Matrix::new_with((n, n), |i, j| {
703                        self.automaton.behavior(
704                            self.prefixes[i]
705                                .iter()
706                                .cloned()
707                                .chain([x])
708                                .chain(self.suffixes[j].iter().cloned()),
709                        )
710                    }) * &self.inv_h
711                })
712                .collect(),
713            final_weights: Matrix::new_with((n, 1), |i, _| {
714                self.automaton.behavior(self.prefixes[i].iter().cloned())
715            }),
716        };
717    }

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

380    pub fn inverse(&self) -> Option<Matrix<R>> {
381        assert_eq!(self.shape.0, self.shape.1);
382        let n = self.shape.0;
383        if n >= 64 {
384            let m = n / 2;
385            let a = Self::new_with((m, m), |i, j| self[i][j].clone());
386            if let Some(mut ai) = a.inverse() {
387                let b = Self::new_with((m, n - m), |i, j| self[i][j + m].clone());
388                let c = Self::new_with((n - m, m), |i, j| self[i + m][j].clone());
389                let mut d = Self::new_with((n - m, n - m), |i, j| self[i + m][j + m].clone());
390                let u = &ai * &b;
391                let v = &c * &ai;
392                d -= &v * &b;
393                let di = d.inverse()?;
394                let r = &u * &di;
395                let t = &di * &v;
396                ai += &r * &v;
397                let mut inverse = Self::zeros((n, n));
398                for i in 0..m {
399                    inverse[i][..m].clone_from_slice(&ai[i]);
400                    for (x, y) in inverse[i][m..].iter_mut().zip(&r[i]) {
401                        *x = R::neg(y);
402                    }
403                }
404                for i in m..n {
405                    for (x, y) in inverse[i][..m].iter_mut().zip(&t[i - m]) {
406                        *x = R::neg(y);
407                    }
408                    inverse[i][m..].clone_from_slice(&di[i - m]);
409                }
410                return Some(inverse);
411            }
412        }
413        let mut a = self.clone();
414        let mut inverse = Self::eye((n, n));
415        let mut ranges: Vec<_> = (0..n).map(|i| (i, i + 1)).collect();
416        for r in 0..n {
417            let pivot = (r..n).find(|&i| !R::is_zero(&a[i][r]))?;
418            a.data.swap(r, pivot);
419            inverse.data.swap(r, pivot);
420            ranges.swap(r, pivot);
421
422            let d = R::inv(&a[r][r]);
423            for x in &mut a[r][r..] {
424                R::mul_assign(x, &d);
425            }
426            let (left, right) = ranges[r];
427            for x in &mut inverse[r][left..right] {
428                R::mul_assign(x, &d);
429            }
430
431            let (a_upper, a_lower) = a.data.split_at_mut(r + 1);
432            let pivot_a = &a_upper[r];
433            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r + 1);
434            let pivot_inverse = &inverse_upper[r];
435            let (ranges_upper, ranges_lower) = ranges.split_at_mut(r + 1);
436            let (left, right) = ranges_upper[r];
437            for ((a, inverse), range) in a_lower.iter_mut().zip(inverse_lower).zip(ranges_lower) {
438                if R::is_zero(&a[r]) {
439                    continue;
440                }
441                let e = a[r].clone();
442                a[r] = R::zero();
443                R::add_scaled_assign(&mut a[(r + 1)..], &pivot_a[(r + 1)..], &R::neg(&e));
444                R::add_scaled_assign(
445                    &mut inverse[left..right],
446                    &pivot_inverse[left..right],
447                    &R::neg(&e),
448                );
449                range.0 = range.0.min(left);
450                range.1 = range.1.max(right);
451            }
452        }
453        for r in (0..n).rev() {
454            let (left, right) = ranges[r];
455            let (inverse_upper, inverse_lower) = inverse.data.split_at_mut(r);
456            let pivot_inverse = &inverse_lower[0];
457            let (ranges_upper, _) = ranges.split_at_mut(r);
458            for ((a, inverse), range) in a.data[..r].iter_mut().zip(inverse_upper).zip(ranges_upper)
459            {
460                if R::is_zero(&a[r]) {
461                    continue;
462                }
463                let e = a[r].clone();
464                a[r] = R::zero();
465                R::add_scaled_assign(
466                    &mut inverse[left..right],
467                    &pivot_inverse[left..right],
468                    &R::neg(&e),
469                );
470                range.0 = range.0.min(left);
471                range.1 = range.1.max(right);
472            }
473        }
474        Some(inverse)
475    }

Source

pub fn characteristic_polynomial(&mut self) -> Vec<R::T>

Examples found in repository ?

crates/library_checker/src/linear_algebra/characteristic_polynomial.rs (line 8)

5pub fn characteristic_polynomial(reader: impl Read, writer: impl Write) {
6    prepare_io!(reader, writer);
7    sc!(n, a: [[M; n]; n]);
8    let p = Matrix::<AddMulOperation<_>>::from_vec(a).characteristic_polynomial();
9    pp!(@it p);
10}

More examples

Hide additional examples

crates/competitive/src/math/mint_matrix.rs (line 124)

84    fn determinant_linear_non_singular(mut self, mut other: Self) -> Option<Vec<MInt<M>>>
85    where
86        M: MIntDotProduct,
87    {
88        let n = self.data.len();
89        let mut f = MInt::one();
90        for d in 0..n {
91            let i = other.data.iter().position(|other| !other[d].is_zero())?;
92            if i != d {
93                self.data.swap(i, d);
94                other.data.swap(i, d);
95                f = -f;
96            }
97            f *= other[d][d];
98            let r = other[d][d].inv();
99            for j in 0..n {
100                self[d][j] *= r;
101                other[d][j] *= r;
102            }
103            assert!(other[d][d].is_one());
104            for i in d + 1..n {
105                let a = other[i][d];
106                for k in 0..n {
107                    self[i][k] = self[i][k] - a * self[d][k];
108                    other[i][k] = other[i][k] - a * other[d][k];
109                }
110            }
111            for j in d + 1..n {
112                let a = other[d][j];
113                for k in 0..n {
114                    self[k][j] = self[k][j] - a * self[k][d];
115                    other[k][j] = other[k][j] - a * other[k][d];
116                }
117            }
118        }
119        for s in self.data.iter_mut() {
120            for s in s.iter_mut() {
121                *s = -*s;
122            }
123        }
124        let mut p = self.characteristic_polynomial();
125        for p in p.iter_mut() {
126            *p *= f;
127        }
128        Some(p)
129    }

Source §

impl<R> Matrix<R>
where R: Ring,

Source

pub fn mul_strassen(&self, rhs: &Matrix<R>) -> Matrix<R>

Examples found in repository ?

crates/library_checker/src/linear_algebra/matrix_product.rs (line 20)

15pub fn matrix_product_strassen(reader: impl Read, writer: impl Write) {
16    prepare_io!(reader, writer);
17    sc!(n, m, k, a: [[M; m]; n], b: [[M; k]; m]);
18    let a = Matrix::<AddMulOperation<_>>::from_vec(a);
19    let b = Matrix::<AddMulOperation<_>>::from_vec(b);
20    let c = a.mul_strassen(&b);
21    pp!(@it2d c.data);
22}

More examples

Hide additional examples

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

928    pub fn pow_strassen(self, mut n: usize) -> Self {
929        assert_eq!(self.shape.0, self.shape.1);
930        let mut res = Matrix::eye(self.shape);
931        let mut x = self;
932        while n > 0 {
933            if n & 1 == 1 {
934                res = res.mul_strassen(&x);
935            }
936            x = x.mul_strassen(&x);
937            n >>= 1;
938        }
939        res
940    }

Source §

pub fn pow_strassen(self, n: usize) -> Self

Examples found in repository ?

crates/library_checker/src/linear_algebra/pow_of_matrix.rs (line 22)

18pub fn pow_of_matrix_strassen(reader: impl Read, writer: impl Write) {
19    prepare_io!(reader, writer);
20    sc!(n, k, a: [[M; n]; n]);
21    let a = Matrix::<AddMulOperation<_>>::from_vec(a);
22    let b = a.pow_strassen(k);
23    pp!(@it2d b.data);
24}

Source §

impl<M> Matrix<AddMulOperation<MInt<M>>>
where M: MIntDotProduct,

Source

fn determinant_linear_non_singular(self, other: Self) -> Option<Vec<MInt<M>>>
where M: MIntDotProduct,

Examples found in repository ?

crates/competitive/src/math/mint_matrix.rs (line 36)

24    fn determinant_linear(mut self, other: Self) -> Option<Vec<MInt<M>>>
25    where
26        M: MIntConvert<usize> + MIntConvert<u64>,
27    {
28        let mut rng = Xorshift::new();
29        let a = MInt::from(rng.rand64());
30        let n = self.data.len();
31        for i in 0..n {
32            for j in 0..n {
33                self[i][j] += other[i][j] * a;
34            }
35        }
36        let mut f = other.determinant_linear_non_singular(self)?;
37        f.reverse();
38        Some(taylor_shift::<M>(f, -a))
39    }