Struct Matrix

Source

pub struct Matrix<R>where
    R: SemiRing,{
    pub shape: (usize, usize),
    pub data: Vec<Vec<R::T>>,
    /* private fields */
}

Fields§

§shape: (usize, usize)§data: Vec<Vec<R::T>>

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/library_checker/src/linear_algebra/characteristic_polynomial.rs (line 10)

6pub fn characteristic_polynomial(reader: impl Read, mut writer: impl Write) {
7    let s = read_all_unchecked(reader);
8    let mut scanner = Scanner::new(&s);
9    scan!(scanner, n, a: [[MInt998244353; n]; n]);
10    let p = Matrix::<AddMulOperation<_>>::from_vec(a).characteristic_polynomial();
11    iter_print!(writer, @it p);
12}

More examples

Hide additional examples

crates/competitive/src/math/matrix.rs (lines 282-284)

270    pub fn inverse(&self) -> Option<Matrix<R>> {
271        assert_eq!(self.shape.0, self.shape.1);
272        let n = self.shape.0;
273        let mut c = Matrix::<R>::zeros((n, n * 2));
274        for i in 0..n {
275            c[i][..n].clone_from_slice(&self[i]);
276            c[i][n + i] = R::one();
277        }
278        c.row_reduction(true);
279        if (0..n).any(|i| R::is_zero(&c[i][i])) {
280            None
281        } else {
282            Some(Self::from_vec(
283                c.data.into_iter().map(|r| r[n..].to_vec()).collect(),
284            ))
285        }
286    }
287
288    pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
289        let n = self.shape.0;
290        if n == 0 {
291            return vec![R::one()];
292        }
293        assert!(self.data.iter().all(|a| a.len() == n));
294        for j in 0..(n - 1) {
295            if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
296                self.data.swap(j + 1, x);
297                self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
298                let inv = R::inv(&self[j + 1][j]);
299                let mut v = vec![];
300                let src = std::mem::take(&mut self[j + 1]);
301                for a in self.data[(j + 2)..].iter_mut() {
302                    let mul = R::mul(&a[j], &inv);
303                    for (a, src) in a[j..].iter_mut().zip(src[j..].iter()) {
304                        R::sub_assign(a, &R::mul(&mul, src));
305                    }
306                    v.push(mul);
307                }
308                self[j + 1] = src;
309                for a in self.data.iter_mut() {
310                    let v = a[(j + 2)..]
311                        .iter()
312                        .zip(v.iter())
313                        .fold(R::zero(), |s, a| R::add(&s, &R::mul(a.0, a.1)));
314                    R::add_assign(&mut a[j + 1], &v);
315                }
316            }
317        }
318        let mut dp = vec![vec![R::one()]];
319        for i in 0..n {
320            let mut next = vec![R::zero(); i + 2];
321            for (j, dp) in dp[i].iter().enumerate() {
322                R::sub_assign(&mut next[j], &R::mul(dp, &self[i][i]));
323                R::add_assign(&mut next[j + 1], dp);
324            }
325            let mut mul = R::one();
326            for j in (0..i).rev() {
327                mul = R::mul(&mul, &self[j + 1][j]);
328                let c = R::mul(&mul, &self[j][i]);
329                for (next, dp) in next.iter_mut().zip(dp[j].iter()) {
330                    R::sub_assign(next, &R::mul(&c, dp));
331                }
332            }
333            dp.push(next);
334        }
335        dp.pop().unwrap()
336    }
337}
338
339impl<R> Index<usize> for Matrix<R>
340where
341    R: SemiRing,
342{
343    type Output = Vec<R::T>;
344    fn index(&self, index: usize) -> &Self::Output {
345        &self.data[index]
346    }
347}
348
349impl<R> IndexMut<usize> for Matrix<R>
350where
351    R: SemiRing,
352{
353    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
354        &mut self.data[index]
355    }
356}
357
358impl<R> Index<(usize, usize)> for Matrix<R>
359where
360    R: SemiRing,
361{
362    type Output = R::T;
363    fn index(&self, index: (usize, usize)) -> &Self::Output {
364        &self.data[index.0][index.1]
365    }
366}
367
368impl<R> IndexMut<(usize, usize)> for Matrix<R>
369where
370    R: SemiRing,
371{
372    fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
373        &mut self.data[index.0][index.1]
374    }
375}
376
377macro_rules! impl_matrix_pairwise_binop {
378    ($imp:ident, $method:ident, $imp_assign:ident, $method_assign:ident $(where [$($clauses:tt)*])?) => {
379        impl<R> $imp_assign for Matrix<R>
380        where
381            R: SemiRing,
382            $($($clauses)*)?
383        {
384            fn $method_assign(&mut self, rhs: Self) {
385                self.pairwise_assign(&rhs, |a, b| R::$method_assign(a, b));
386            }
387        }
388        impl<R> $imp_assign<&Matrix<R>> for Matrix<R>
389        where
390            R: SemiRing,
391            $($($clauses)*)?
392        {
393            fn $method_assign(&mut self, rhs: &Self) {
394                self.pairwise_assign(rhs, |a, b| R::$method_assign(a, b));
395            }
396        }
397        impl<R> $imp for Matrix<R>
398        where
399            R: SemiRing,
400            $($($clauses)*)?
401        {
402            type Output = Matrix<R>;
403            fn $method(mut self, rhs: Self) -> Self::Output {
404                self.$method_assign(rhs);
405                self
406            }
407        }
408        impl<R> $imp<&Matrix<R>> for Matrix<R>
409        where
410            R: SemiRing,
411            $($($clauses)*)?
412        {
413            type Output = Matrix<R>;
414            fn $method(mut self, rhs: &Self) -> Self::Output {
415                self.$method_assign(rhs);
416                self
417            }
418        }
419        impl<R> $imp<Matrix<R>> for &Matrix<R>
420        where
421            R: SemiRing,
422            $($($clauses)*)?
423        {
424            type Output = Matrix<R>;
425            fn $method(self, mut rhs: Matrix<R>) -> Self::Output {
426                rhs.pairwise_assign(self, |a, b| *a = R::$method(b, a));
427                rhs
428            }
429        }
430        impl<R> $imp<&Matrix<R>> for &Matrix<R>
431        where
432            R: SemiRing,
433            $($($clauses)*)?
434        {
435            type Output = Matrix<R>;
436            fn $method(self, rhs: &Matrix<R>) -> Self::Output {
437                let mut this = self.clone();
438                this.$method_assign(rhs);
439                this
440            }
441        }
442    };
443}
444
445impl_matrix_pairwise_binop!(Add, add, AddAssign, add_assign);
446impl_matrix_pairwise_binop!(Sub, sub, SubAssign, sub_assign where [R::Additive: Invertible]);
447
448impl<R> Mul for Matrix<R>
449where
450    R: SemiRing,
451{
452    type Output = Matrix<R>;
453    fn mul(self, rhs: Self) -> Self::Output {
454        (&self).mul(&rhs)
455    }
456}
457impl<R> Mul<&Matrix<R>> for Matrix<R>
458where
459    R: SemiRing,
460{
461    type Output = Matrix<R>;
462    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
463        (&self).mul(rhs)
464    }
465}
466impl<R> Mul<Matrix<R>> for &Matrix<R>
467where
468    R: SemiRing,
469{
470    type Output = Matrix<R>;
471    fn mul(self, rhs: Matrix<R>) -> Self::Output {
472        self.mul(&rhs)
473    }
474}
475impl<R> Mul<&Matrix<R>> for &Matrix<R>
476where
477    R: SemiRing,
478{
479    type Output = Matrix<R>;
480    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
481        assert_eq!(self.shape.1, rhs.shape.0);
482        let mut res = Matrix::zeros((self.shape.0, rhs.shape.1));
483        for i in 0..self.shape.0 {
484            for k in 0..self.shape.1 {
485                for j in 0..rhs.shape.1 {
486                    R::add_assign(&mut res[i][j], &R::mul(&self[i][k], &rhs[k][j]));
487                }
488            }
489        }
490        res
491    }
492}
493
494impl<R> MulAssign<&R::T> for Matrix<R>
495where
496    R: SemiRing,
497{
498    fn mul_assign(&mut self, rhs: &R::T) {
499        for i in 0..self.shape.0 {
500            for j in 0..self.shape.1 {
501                R::mul_assign(&mut self[(i, j)], rhs);
502            }
503        }
504    }
505}
506
507impl<R> Matrix<R>
508where
509    R: SemiRing,
510{
511    pub fn pow(self, mut n: usize) -> Self {
512        assert_eq!(self.shape.0, self.shape.1);
513        let mut res = Matrix::eye(self.shape);
514        let mut x = self;
515        while n > 0 {
516            if n & 1 == 1 {
517                res = &res * &x;
518            }
519            x = &x * &x;
520            n >>= 1;
521        }
522        res
523    }
524}
525
526impl<R> SerdeByteStr for Matrix<R>
527where
528    R: SemiRing,
529    R::T: SerdeByteStr,
530{
531    fn serialize(&self, buf: &mut Vec<u8>) {
532        self.data.serialize(buf);
533    }
534
535    fn deserialize<I>(iter: &mut I) -> Self
536    where
537        I: Iterator<Item = u8>,
538    {
539        Self::from_vec(Vec::deserialize(iter))
540    }

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

6pub fn system_of_linear_equations(reader: impl Read, mut writer: impl Write) {
7    let s = read_all_unchecked(reader);
8    let mut scanner = Scanner::new(&s);
9    scan!(scanner, n, m, a: [[MInt998244353; m]; n], b: [MInt998244353; n]);
10    let a = Matrix::<AddMulOperation<MInt998244353>>::from_vec(a);
11    if let Some(sol) = a.solve_system_of_linear_equations(&b) {
12        iter_print!(writer, sol.basis.len(); @it &sol.particular);
13        for b in sol.basis {
14            iter_print!(writer, @it &b);
15        }
16    } else {
17        iter_print!(writer, -1);
18    }
19}

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

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

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

114    pub fn transpose(&self) -> Self {
115        Self::new_with((self.shape.1, self.shape.0), |i, j| self[j][i].clone())
116    }
117
118    pub fn map<S, F>(&self, mut f: F) -> Matrix<S>
119    where
120        S: SemiRing,
121        F: FnMut(&R::T) -> S::T,
122    {
123        Matrix::<S>::new_with(self.shape, |i, j| f(&self[i][j]))
124    }

More examples

Hide additional examples

crates/competitive/src/algorithm/automata_learning.rs (lines 554-561)

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

Source

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

Examples found in repository ?

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

129    fn from(smat: SparseMatrix<R>) -> Self {
130        let mut mat = Matrix::zeros(smat.shape);
131        for &(i, j, ref v) in &smat.nonzero {
132            R::add_assign(&mut mat[(i, j)], v);
133        }
134        mat
135    }

More examples

Hide additional examples

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

495    pub fn new(automaton: A) -> Self {
496        let sigma = automaton.sigma();
497        Self {
498            automaton,
499            prefixes: vec![],
500            suffixes: vec![],
501            inv_h: Matrix::zeros((0, 0)),
502            nh: vec![Matrix::zeros((0, 0)); sigma],
503            wfa: WeightedFiniteAutomaton {
504                initial_weights: Matrix::zeros((1, 0)),
505                transitions: vec![Matrix::zeros((0, 0)); sigma],
506                final_weights: Matrix::zeros((0, 1)),
507            },
508            _marker: PhantomData,
509        }
510    }

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

235    pub fn solve_system_of_linear_equations(
236        &self,
237        b: &[R::T],
238    ) -> Option<SystemOfLinearEquationsSolution<R>> {
239        assert_eq!(self.shape.0, b.len());
240        let (n, m) = self.shape;
241        let mut c = Matrix::<R>::zeros((n, m + 1));
242        for i in 0..n {
243            c[i][..m].clone_from_slice(&self[i]);
244            c[i][m] = b[i].clone();
245        }
246        let mut reduced = vec![!0; m + 1];
247        c.row_reduction_with(true, |r, _, c| reduced[c] = r);
248        if reduced[m] != !0 {
249            return None;
250        }
251        let mut particular = vec![R::zero(); m];
252        let mut basis = vec![];
253        for j in 0..m {
254            if reduced[j] != !0 {
255                particular[j] = c[reduced[j]][m].clone();
256            } else {
257                let mut v = vec![R::zero(); m];
258                v[j] = R::one();
259                for i in 0..m {
260                    if reduced[i] != !0 {
261                        R::sub_assign(&mut v[i], &c[reduced[i]][j]);
262                    }
263                }
264                basis.push(v);
265            }
266        }
267        Some(SystemOfLinearEquationsSolution { particular, basis })
268    }
269
270    pub fn inverse(&self) -> Option<Matrix<R>> {
271        assert_eq!(self.shape.0, self.shape.1);
272        let n = self.shape.0;
273        let mut c = Matrix::<R>::zeros((n, n * 2));
274        for i in 0..n {
275            c[i][..n].clone_from_slice(&self[i]);
276            c[i][n + i] = R::one();
277        }
278        c.row_reduction(true);
279        if (0..n).any(|i| R::is_zero(&c[i][i])) {
280            None
281        } else {
282            Some(Self::from_vec(
283                c.data.into_iter().map(|r| r[n..].to_vec()).collect(),
284            ))
285        }
286    }
287
288    pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
289        let n = self.shape.0;
290        if n == 0 {
291            return vec![R::one()];
292        }
293        assert!(self.data.iter().all(|a| a.len() == n));
294        for j in 0..(n - 1) {
295            if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
296                self.data.swap(j + 1, x);
297                self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
298                let inv = R::inv(&self[j + 1][j]);
299                let mut v = vec![];
300                let src = std::mem::take(&mut self[j + 1]);
301                for a in self.data[(j + 2)..].iter_mut() {
302                    let mul = R::mul(&a[j], &inv);
303                    for (a, src) in a[j..].iter_mut().zip(src[j..].iter()) {
304                        R::sub_assign(a, &R::mul(&mul, src));
305                    }
306                    v.push(mul);
307                }
308                self[j + 1] = src;
309                for a in self.data.iter_mut() {
310                    let v = a[(j + 2)..]
311                        .iter()
312                        .zip(v.iter())
313                        .fold(R::zero(), |s, a| R::add(&s, &R::mul(a.0, a.1)));
314                    R::add_assign(&mut a[j + 1], &v);
315                }
316            }
317        }
318        let mut dp = vec![vec![R::one()]];
319        for i in 0..n {
320            let mut next = vec![R::zero(); i + 2];
321            for (j, dp) in dp[i].iter().enumerate() {
322                R::sub_assign(&mut next[j], &R::mul(dp, &self[i][i]));
323                R::add_assign(&mut next[j + 1], dp);
324            }
325            let mut mul = R::one();
326            for j in (0..i).rev() {
327                mul = R::mul(&mul, &self[j + 1][j]);
328                let c = R::mul(&mul, &self[j][i]);
329                for (next, dp) in next.iter_mut().zip(dp[j].iter()) {
330                    R::sub_assign(next, &R::mul(&c, dp));
331                }
332            }
333            dp.push(next);
334        }
335        dp.pop().unwrap()
336    }
337}
338
339impl<R> Index<usize> for Matrix<R>
340where
341    R: SemiRing,
342{
343    type Output = Vec<R::T>;
344    fn index(&self, index: usize) -> &Self::Output {
345        &self.data[index]
346    }
347}
348
349impl<R> IndexMut<usize> for Matrix<R>
350where
351    R: SemiRing,
352{
353    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
354        &mut self.data[index]
355    }
356}
357
358impl<R> Index<(usize, usize)> for Matrix<R>
359where
360    R: SemiRing,
361{
362    type Output = R::T;
363    fn index(&self, index: (usize, usize)) -> &Self::Output {
364        &self.data[index.0][index.1]
365    }
366}
367
368impl<R> IndexMut<(usize, usize)> for Matrix<R>
369where
370    R: SemiRing,
371{
372    fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
373        &mut self.data[index.0][index.1]
374    }
375}
376
377macro_rules! impl_matrix_pairwise_binop {
378    ($imp:ident, $method:ident, $imp_assign:ident, $method_assign:ident $(where [$($clauses:tt)*])?) => {
379        impl<R> $imp_assign for Matrix<R>
380        where
381            R: SemiRing,
382            $($($clauses)*)?
383        {
384            fn $method_assign(&mut self, rhs: Self) {
385                self.pairwise_assign(&rhs, |a, b| R::$method_assign(a, b));
386            }
387        }
388        impl<R> $imp_assign<&Matrix<R>> for Matrix<R>
389        where
390            R: SemiRing,
391            $($($clauses)*)?
392        {
393            fn $method_assign(&mut self, rhs: &Self) {
394                self.pairwise_assign(rhs, |a, b| R::$method_assign(a, b));
395            }
396        }
397        impl<R> $imp for Matrix<R>
398        where
399            R: SemiRing,
400            $($($clauses)*)?
401        {
402            type Output = Matrix<R>;
403            fn $method(mut self, rhs: Self) -> Self::Output {
404                self.$method_assign(rhs);
405                self
406            }
407        }
408        impl<R> $imp<&Matrix<R>> for Matrix<R>
409        where
410            R: SemiRing,
411            $($($clauses)*)?
412        {
413            type Output = Matrix<R>;
414            fn $method(mut self, rhs: &Self) -> Self::Output {
415                self.$method_assign(rhs);
416                self
417            }
418        }
419        impl<R> $imp<Matrix<R>> for &Matrix<R>
420        where
421            R: SemiRing,
422            $($($clauses)*)?
423        {
424            type Output = Matrix<R>;
425            fn $method(self, mut rhs: Matrix<R>) -> Self::Output {
426                rhs.pairwise_assign(self, |a, b| *a = R::$method(b, a));
427                rhs
428            }
429        }
430        impl<R> $imp<&Matrix<R>> for &Matrix<R>
431        where
432            R: SemiRing,
433            $($($clauses)*)?
434        {
435            type Output = Matrix<R>;
436            fn $method(self, rhs: &Matrix<R>) -> Self::Output {
437                let mut this = self.clone();
438                this.$method_assign(rhs);
439                this
440            }
441        }
442    };
443}
444
445impl_matrix_pairwise_binop!(Add, add, AddAssign, add_assign);
446impl_matrix_pairwise_binop!(Sub, sub, SubAssign, sub_assign where [R::Additive: Invertible]);
447
448impl<R> Mul for Matrix<R>
449where
450    R: SemiRing,
451{
452    type Output = Matrix<R>;
453    fn mul(self, rhs: Self) -> Self::Output {
454        (&self).mul(&rhs)
455    }
456}
457impl<R> Mul<&Matrix<R>> for Matrix<R>
458where
459    R: SemiRing,
460{
461    type Output = Matrix<R>;
462    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
463        (&self).mul(rhs)
464    }
465}
466impl<R> Mul<Matrix<R>> for &Matrix<R>
467where
468    R: SemiRing,
469{
470    type Output = Matrix<R>;
471    fn mul(self, rhs: Matrix<R>) -> Self::Output {
472        self.mul(&rhs)
473    }
474}
475impl<R> Mul<&Matrix<R>> for &Matrix<R>
476where
477    R: SemiRing,
478{
479    type Output = Matrix<R>;
480    fn mul(self, rhs: &Matrix<R>) -> Self::Output {
481        assert_eq!(self.shape.1, rhs.shape.0);
482        let mut res = Matrix::zeros((self.shape.0, rhs.shape.1));
483        for i in 0..self.shape.0 {
484            for k in 0..self.shape.1 {
485                for j in 0..rhs.shape.1 {
486                    R::add_assign(&mut res[i][j], &R::mul(&self[i][k], &rhs[k][j]));
487                }
488            }
489        }
490        res
491    }

Source

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

Examples found in repository ?

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

511    pub fn pow(self, mut n: usize) -> Self {
512        assert_eq!(self.shape.0, self.shape.1);
513        let mut res = Matrix::eye(self.shape);
514        let mut x = self;
515        while n > 0 {
516            if n & 1 == 1 {
517                res = &res * &x;
518            }
519            x = &x * &x;
520            n >>= 1;
521        }
522        res
523    }

Source

pub fn transpose(&self) -> Self

Source

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

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 585-591)

545    pub fn train_sample(&mut self, sample: &[usize]) -> bool {
546        let Some((prefix, suffix)) = self.split_sample(sample) else {
547            return false;
548        };
549        self.prefixes.push(prefix);
550        self.suffixes.push(suffix);
551        let n = self.inv_h.shape.0;
552        let prefix = &self.prefixes[n];
553        let suffix = &self.suffixes[n];
554        let u = Matrix::<F>::new_with((n, 1), |i, _| {
555            self.automaton.behavior(
556                self.prefixes[i]
557                    .iter()
558                    .cloned()
559                    .chain(suffix.iter().cloned()),
560            )
561        });
562        let v = Matrix::<F>::new_with((1, n), |_, j| {
563            self.automaton.behavior(
564                prefix
565                    .iter()
566                    .cloned()
567                    .chain(self.suffixes[j].iter().cloned()),
568            )
569        });
570        let w = Matrix::<F>::new_with((1, 1), |_, _| {
571            self.automaton
572                .behavior(prefix.iter().cloned().chain(suffix.iter().cloned()))
573        });
574        let t = &self.inv_h * &u;
575        let s = &v * &self.inv_h;
576        let d = F::inv(&(&w - &(&v * &t))[0][0]);
577        let dh = &t * &s;
578        for i in 0..n {
579            for j in 0..n {
580                F::add_assign(&mut self.inv_h[i][j], &F::mul(&dh[i][j], &d));
581            }
582        }
583        self.inv_h
584            .add_col_with(|i, _| F::neg(&F::mul(&t[i][0], &d)));
585        self.inv_h.add_row_with(|_, j| {
586            if j != n {
587                F::neg(&F::mul(&s[0][j], &d))
588            } else {
589                d.clone()
590            }
591        });
592
593        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
594            let b = &(&self.nh[x] * &t) * &s;
595            for i in 0..n {
596                for j in 0..n {
597                    F::add_assign(&mut transition[i][j], &F::mul(&b[i][j], &d));
598                }
599            }
600        }
601        for (x, nh) in self.nh.iter_mut().enumerate() {
602            nh.add_col_with(|i, j| {
603                self.automaton.behavior(
604                    self.prefixes[i]
605                        .iter()
606                        .cloned()
607                        .chain([x])
608                        .chain(self.suffixes[j].iter().cloned()),
609                )
610            });
611            nh.add_row_with(|i, j| {
612                self.automaton.behavior(
613                    self.prefixes[i]
614                        .iter()
615                        .cloned()
616                        .chain([x])
617                        .chain(self.suffixes[j].iter().cloned()),
618                )
619            });
620        }
621        self.wfa
622            .initial_weights
623            .add_col_with(|_, _| if n == 0 { F::one() } else { F::zero() });
624        self.wfa
625            .final_weights
626            .add_row_with(|_, _| self.automaton.behavior(prefix.iter().cloned()));
627        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
628            transition.add_col_with(|_, _| F::zero());
629            transition.add_row_with(|_, _| F::zero());
630            for i in 0..=n {
631                for j in 0..=n {
632                    if i == n || j == n {
633                        for k in 0..=n {
634                            if i != n && j != n && k != n {
635                                continue;
636                            }
637                            F::add_assign(
638                                &mut transition[i][k],
639                                &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
640                            );
641                        }
642                    } else {
643                        let k = n;
644                        F::add_assign(
645                            &mut transition[i][k],
646                            &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
647                        );
648                    }
649                }
650            }
651        }
652        true
653    }

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

545    pub fn train_sample(&mut self, sample: &[usize]) -> bool {
546        let Some((prefix, suffix)) = self.split_sample(sample) else {
547            return false;
548        };
549        self.prefixes.push(prefix);
550        self.suffixes.push(suffix);
551        let n = self.inv_h.shape.0;
552        let prefix = &self.prefixes[n];
553        let suffix = &self.suffixes[n];
554        let u = Matrix::<F>::new_with((n, 1), |i, _| {
555            self.automaton.behavior(
556                self.prefixes[i]
557                    .iter()
558                    .cloned()
559                    .chain(suffix.iter().cloned()),
560            )
561        });
562        let v = Matrix::<F>::new_with((1, n), |_, j| {
563            self.automaton.behavior(
564                prefix
565                    .iter()
566                    .cloned()
567                    .chain(self.suffixes[j].iter().cloned()),
568            )
569        });
570        let w = Matrix::<F>::new_with((1, 1), |_, _| {
571            self.automaton
572                .behavior(prefix.iter().cloned().chain(suffix.iter().cloned()))
573        });
574        let t = &self.inv_h * &u;
575        let s = &v * &self.inv_h;
576        let d = F::inv(&(&w - &(&v * &t))[0][0]);
577        let dh = &t * &s;
578        for i in 0..n {
579            for j in 0..n {
580                F::add_assign(&mut self.inv_h[i][j], &F::mul(&dh[i][j], &d));
581            }
582        }
583        self.inv_h
584            .add_col_with(|i, _| F::neg(&F::mul(&t[i][0], &d)));
585        self.inv_h.add_row_with(|_, j| {
586            if j != n {
587                F::neg(&F::mul(&s[0][j], &d))
588            } else {
589                d.clone()
590            }
591        });
592
593        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
594            let b = &(&self.nh[x] * &t) * &s;
595            for i in 0..n {
596                for j in 0..n {
597                    F::add_assign(&mut transition[i][j], &F::mul(&b[i][j], &d));
598                }
599            }
600        }
601        for (x, nh) in self.nh.iter_mut().enumerate() {
602            nh.add_col_with(|i, j| {
603                self.automaton.behavior(
604                    self.prefixes[i]
605                        .iter()
606                        .cloned()
607                        .chain([x])
608                        .chain(self.suffixes[j].iter().cloned()),
609                )
610            });
611            nh.add_row_with(|i, j| {
612                self.automaton.behavior(
613                    self.prefixes[i]
614                        .iter()
615                        .cloned()
616                        .chain([x])
617                        .chain(self.suffixes[j].iter().cloned()),
618                )
619            });
620        }
621        self.wfa
622            .initial_weights
623            .add_col_with(|_, _| if n == 0 { F::one() } else { F::zero() });
624        self.wfa
625            .final_weights
626            .add_row_with(|_, _| self.automaton.behavior(prefix.iter().cloned()));
627        for (x, transition) in self.wfa.transitions.iter_mut().enumerate() {
628            transition.add_col_with(|_, _| F::zero());
629            transition.add_row_with(|_, _| F::zero());
630            for i in 0..=n {
631                for j in 0..=n {
632                    if i == n || j == n {
633                        for k in 0..=n {
634                            if i != n && j != n && k != n {
635                                continue;
636                            }
637                            F::add_assign(
638                                &mut transition[i][k],
639                                &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
640                            );
641                        }
642                    } else {
643                        let k = n;
644                        F::add_assign(
645                            &mut transition[i][k],
646                            &F::mul(&self.nh[x][i][j], &self.inv_h[j][k]),
647                        );
648                    }
649                }
650            }
651        }
652        true
653    }

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, R::Additive: Invertible, R::Multiplicative: Invertible, R::T: PartialEq,

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

209    pub fn row_reduction(&mut self, normalize: bool) {
210        self.row_reduction_with(normalize, |_, _, _| {});
211    }
212
213    pub fn rank(&mut self) -> usize {
214        let n = self.shape.0;
215        self.row_reduction(false);
216        (0..n)
217            .filter(|&i| !self.data[i].iter().all(|x| R::is_zero(x)))
218            .count()
219    }
220
221    pub fn determinant(&mut self) -> R::T {
222        assert_eq!(self.shape.0, self.shape.1);
223        let mut neg = false;
224        self.row_reduction_with(false, |r, p, _| neg ^= r != p);
225        let mut d = R::one();
226        if neg {
227            d = R::neg(&d);
228        }
229        for i in 0..self.shape.0 {
230            R::mul_assign(&mut d, &self[i][i]);
231        }
232        d
233    }
234
235    pub fn solve_system_of_linear_equations(
236        &self,
237        b: &[R::T],
238    ) -> Option<SystemOfLinearEquationsSolution<R>> {
239        assert_eq!(self.shape.0, b.len());
240        let (n, m) = self.shape;
241        let mut c = Matrix::<R>::zeros((n, m + 1));
242        for i in 0..n {
243            c[i][..m].clone_from_slice(&self[i]);
244            c[i][m] = b[i].clone();
245        }
246        let mut reduced = vec![!0; m + 1];
247        c.row_reduction_with(true, |r, _, c| reduced[c] = r);
248        if reduced[m] != !0 {
249            return None;
250        }
251        let mut particular = vec![R::zero(); m];
252        let mut basis = vec![];
253        for j in 0..m {
254            if reduced[j] != !0 {
255                particular[j] = c[reduced[j]][m].clone();
256            } else {
257                let mut v = vec![R::zero(); m];
258                v[j] = R::one();
259                for i in 0..m {
260                    if reduced[i] != !0 {
261                        R::sub_assign(&mut v[i], &c[reduced[i]][j]);
262                    }
263                }
264                basis.push(v);
265            }
266        }
267        Some(SystemOfLinearEquationsSolution { particular, basis })
268    }

More examples

Hide additional examples

crates/competitive/src/algorithm/automata_learning.rs (lines 691-694)

659    pub fn batch_train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
660        let mut prefix_set: HashSet<_> = self.prefixes.iter().cloned().collect();
661        let mut suffix_set: HashSet<_> = self.suffixes.iter().cloned().collect();
662        for sample in samples {
663            if prefix_set.insert(sample.to_vec()) {
664                self.prefixes.push(sample.to_vec());
665            }
666            if suffix_set.insert(sample.to_vec()) {
667                self.suffixes.push(sample);
668            }
669        }
670        let mut h = Matrix::<F>::new_with((self.prefixes.len(), self.suffixes.len()), |i, j| {
671            self.automaton.behavior(
672                self.prefixes[i]
673                    .iter()
674                    .cloned()
675                    .chain(self.suffixes[j].iter().cloned()),
676            )
677        });
678        if !self.prefixes.is_empty() && !self.suffixes.is_empty() && F::is_zero(&h[0][0]) {
679            for j in 1..self.suffixes.len() {
680                if !F::is_zero(&h[0][j]) {
681                    self.suffixes.swap(0, j);
682                    for i in 0..self.prefixes.len() {
683                        h.data[i].swap(0, j);
684                    }
685                    break;
686                }
687            }
688        }
689        let mut row_id: Vec<usize> = (0..h.shape.0).collect();
690        let mut pivots = vec![];
691        h.row_reduction_with(false, |r, p, c| {
692            row_id.swap(r, p);
693            pivots.push((row_id[r], c));
694        });
695        let mut new_prefixes = vec![];
696        let mut new_suffixes = vec![];
697        for (i, j) in pivots {
698            new_prefixes.push(self.prefixes[i].clone());
699            new_suffixes.push(self.suffixes[j].clone());
700        }
701        self.prefixes = new_prefixes;
702        self.suffixes = new_suffixes;
703        assert_eq!(self.prefixes.len(), self.suffixes.len());
704        let n = self.prefixes.len();
705        let h = Matrix::<F>::new_with((n, n), |i, j| {
706            self.automaton.behavior(
707                self.prefixes[i]
708                    .iter()
709                    .cloned()
710                    .chain(self.suffixes[j].iter().cloned()),
711            )
712        });
713        self.inv_h = h.inverse().expect("Hankel matrix must be invertible");
714        self.wfa = WeightedFiniteAutomaton::<F> {
715            initial_weights: Matrix::new_with((1, n), |_, j| {
716                if self.prefixes[j].is_empty() {
717                    F::one()
718                } else {
719                    F::zero()
720                }
721            }),
722            transitions: (0..self.automaton.sigma())
723                .map(|x| {
724                    &Matrix::new_with((n, n), |i, j| {
725                        self.automaton.behavior(
726                            self.prefixes[i]
727                                .iter()
728                                .cloned()
729                                .chain([x])
730                                .chain(self.suffixes[j].iter().cloned()),
731                        )
732                    }) * &self.inv_h
733                })
734                .collect(),
735            final_weights: Matrix::new_with((n, 1), |i, _| {
736                self.automaton.behavior(self.prefixes[i].iter().cloned())
737            }),
738        };
739    }

Source

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

Examples found in repository ?

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

213    pub fn rank(&mut self) -> usize {
214        let n = self.shape.0;
215        self.row_reduction(false);
216        (0..n)
217            .filter(|&i| !self.data[i].iter().all(|x| R::is_zero(x)))
218            .count()
219    }
220
221    pub fn determinant(&mut self) -> R::T {
222        assert_eq!(self.shape.0, self.shape.1);
223        let mut neg = false;
224        self.row_reduction_with(false, |r, p, _| neg ^= r != p);
225        let mut d = R::one();
226        if neg {
227            d = R::neg(&d);
228        }
229        for i in 0..self.shape.0 {
230            R::mul_assign(&mut d, &self[i][i]);
231        }
232        d
233    }
234
235    pub fn solve_system_of_linear_equations(
236        &self,
237        b: &[R::T],
238    ) -> Option<SystemOfLinearEquationsSolution<R>> {
239        assert_eq!(self.shape.0, b.len());
240        let (n, m) = self.shape;
241        let mut c = Matrix::<R>::zeros((n, m + 1));
242        for i in 0..n {
243            c[i][..m].clone_from_slice(&self[i]);
244            c[i][m] = b[i].clone();
245        }
246        let mut reduced = vec![!0; m + 1];
247        c.row_reduction_with(true, |r, _, c| reduced[c] = r);
248        if reduced[m] != !0 {
249            return None;
250        }
251        let mut particular = vec![R::zero(); m];
252        let mut basis = vec![];
253        for j in 0..m {
254            if reduced[j] != !0 {
255                particular[j] = c[reduced[j]][m].clone();
256            } else {
257                let mut v = vec![R::zero(); m];
258                v[j] = R::one();
259                for i in 0..m {
260                    if reduced[i] != !0 {
261                        R::sub_assign(&mut v[i], &c[reduced[i]][j]);
262                    }
263                }
264                basis.push(v);
265            }
266        }
267        Some(SystemOfLinearEquationsSolution { particular, basis })
268    }
269
270    pub fn inverse(&self) -> Option<Matrix<R>> {
271        assert_eq!(self.shape.0, self.shape.1);
272        let n = self.shape.0;
273        let mut c = Matrix::<R>::zeros((n, n * 2));
274        for i in 0..n {
275            c[i][..n].clone_from_slice(&self[i]);
276            c[i][n + i] = R::one();
277        }
278        c.row_reduction(true);
279        if (0..n).any(|i| R::is_zero(&c[i][i])) {
280            None
281        } else {
282            Some(Self::from_vec(
283                c.data.into_iter().map(|r| r[n..].to_vec()).collect(),
284            ))
285        }
286    }

Source

pub fn rank(&mut self) -> usize

Source

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

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

6pub fn system_of_linear_equations(reader: impl Read, mut writer: impl Write) {
7    let s = read_all_unchecked(reader);
8    let mut scanner = Scanner::new(&s);
9    scan!(scanner, n, m, a: [[MInt998244353; m]; n], b: [MInt998244353; n]);
10    let a = Matrix::<AddMulOperation<MInt998244353>>::from_vec(a);
11    if let Some(sol) = a.solve_system_of_linear_equations(&b) {
12        iter_print!(writer, sol.basis.len(); @it &sol.particular);
13        for b in sol.basis {
14            iter_print!(writer, @it &b);
15        }
16    } else {
17        iter_print!(writer, -1);
18    }
19}

More examples

Hide additional examples

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

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

Source

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

Examples found in repository ?

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

659    pub fn batch_train(&mut self, samples: impl IntoIterator<Item = Vec<usize>>) {
660        let mut prefix_set: HashSet<_> = self.prefixes.iter().cloned().collect();
661        let mut suffix_set: HashSet<_> = self.suffixes.iter().cloned().collect();
662        for sample in samples {
663            if prefix_set.insert(sample.to_vec()) {
664                self.prefixes.push(sample.to_vec());
665            }
666            if suffix_set.insert(sample.to_vec()) {
667                self.suffixes.push(sample);
668            }
669        }
670        let mut h = Matrix::<F>::new_with((self.prefixes.len(), self.suffixes.len()), |i, j| {
671            self.automaton.behavior(
672                self.prefixes[i]
673                    .iter()
674                    .cloned()
675                    .chain(self.suffixes[j].iter().cloned()),
676            )
677        });
678        if !self.prefixes.is_empty() && !self.suffixes.is_empty() && F::is_zero(&h[0][0]) {
679            for j in 1..self.suffixes.len() {
680                if !F::is_zero(&h[0][j]) {
681                    self.suffixes.swap(0, j);
682                    for i in 0..self.prefixes.len() {
683                        h.data[i].swap(0, j);
684                    }
685                    break;
686                }
687            }
688        }
689        let mut row_id: Vec<usize> = (0..h.shape.0).collect();
690        let mut pivots = vec![];
691        h.row_reduction_with(false, |r, p, c| {
692            row_id.swap(r, p);
693            pivots.push((row_id[r], c));
694        });
695        let mut new_prefixes = vec![];
696        let mut new_suffixes = vec![];
697        for (i, j) in pivots {
698            new_prefixes.push(self.prefixes[i].clone());
699            new_suffixes.push(self.suffixes[j].clone());
700        }
701        self.prefixes = new_prefixes;
702        self.suffixes = new_suffixes;
703        assert_eq!(self.prefixes.len(), self.suffixes.len());
704        let n = self.prefixes.len();
705        let h = Matrix::<F>::new_with((n, n), |i, j| {
706            self.automaton.behavior(
707                self.prefixes[i]
708                    .iter()
709                    .cloned()
710                    .chain(self.suffixes[j].iter().cloned()),
711            )
712        });
713        self.inv_h = h.inverse().expect("Hankel matrix must be invertible");
714        self.wfa = WeightedFiniteAutomaton::<F> {
715            initial_weights: Matrix::new_with((1, n), |_, j| {
716                if self.prefixes[j].is_empty() {
717                    F::one()
718                } else {
719                    F::zero()
720                }
721            }),
722            transitions: (0..self.automaton.sigma())
723                .map(|x| {
724                    &Matrix::new_with((n, n), |i, j| {
725                        self.automaton.behavior(
726                            self.prefixes[i]
727                                .iter()
728                                .cloned()
729                                .chain([x])
730                                .chain(self.suffixes[j].iter().cloned()),
731                        )
732                    }) * &self.inv_h
733                })
734                .collect(),
735            final_weights: Matrix::new_with((n, 1), |i, _| {
736                self.automaton.behavior(self.prefixes[i].iter().cloned())
737            }),
738        };
739    }

Source

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

Examples found in repository ?

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

6pub fn characteristic_polynomial(reader: impl Read, mut writer: impl Write) {
7    let s = read_all_unchecked(reader);
8    let mut scanner = Scanner::new(&s);
9    scan!(scanner, n, a: [[MInt998244353; n]; n]);
10    let p = Matrix::<AddMulOperation<_>>::from_vec(a).characteristic_polynomial();
11    iter_print!(writer, @it p);
12}

More examples

Hide additional examples

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

41    fn determinant_linear_non_singular(mut self, mut other: Self) -> Option<Vec<MInt<M>>>
42    where
43        M: MIntBase,
44    {
45        let n = self.data.len();
46        let mut f = MInt::one();
47        for d in 0..n {
48            let i = other.data.iter().position(|other| !other[d].is_zero())?;
49            if i != d {
50                self.data.swap(i, d);
51                other.data.swap(i, d);
52                f = -f;
53            }
54            f *= other[d][d];
55            let r = other[d][d].inv();
56            for j in 0..n {
57                self[d][j] *= r;
58                other[d][j] *= r;
59            }
60            assert!(other[d][d].is_one());
61            for i in d + 1..n {
62                let a = other[i][d];
63                for k in 0..n {
64                    self[i][k] = self[i][k] - a * self[d][k];
65                    other[i][k] = other[i][k] - a * other[d][k];
66                }
67            }
68            for j in d + 1..n {
69                let a = other[d][j];
70                for k in 0..n {
71                    self[k][j] = self[k][j] - a * self[k][d];
72                    other[k][j] = other[k][j] - a * other[k][d];
73                }
74            }
75        }
76        for s in self.data.iter_mut() {
77            for s in s.iter_mut() {
78                *s = -*s;
79            }
80        }
81        let mut p = self.characteristic_polynomial();
82        for p in p.iter_mut() {
83            *p *= f;
84        }
85        Some(p)
86    }