pub trait SemiRing {
type T: Clone;
type Additive: AbelianMonoid<T = Self::T>;
type Multiplicative: Monoid<T = Self::T>;
// Provided methods
fn zero() -> Self::T { ... }
fn is_zero(x: &Self::T) -> bool
where Self::T: PartialEq { ... }
fn one() -> Self::T { ... }
fn is_one(x: &Self::T) -> bool
where Self::T: PartialEq { ... }
fn add(x: &Self::T, y: &Self::T) -> Self::T { ... }
fn mul(x: &Self::T, y: &Self::T) -> Self::T { ... }
fn add_assign(x: &mut Self::T, y: &Self::T) { ... }
fn mul_assign(x: &mut Self::T, y: &Self::T) { ... }
}Required Associated Types§
type T: Clone
type Additive: AbelianMonoid<T = Self::T>
type Multiplicative: Monoid<T = Self::T>
Provided Methods§
Sourcefn zero() -> Self::T
fn zero() -> Self::T
additive identity: $0$
Examples found in repository?
More examples
crates/competitive/src/math/matrix.rs (line 90)
87 pub fn zeros(shape: (usize, usize)) -> Self {
88 Self {
89 shape,
90 data: vec![vec![R::zero(); shape.1]; shape.0],
91 _marker: PhantomData,
92 }
93 }
94
95 pub fn eye(shape: (usize, usize)) -> Self {
96 let mut data = vec![vec![R::zero(); shape.1]; shape.0];
97 for (i, d) in data.iter_mut().enumerate().take(shape.1) {
98 d[i] = R::one();
99 }
100 Self {
101 shape,
102 data,
103 _marker: PhantomData,
104 }
105 }
106
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 /// f: (row, pivot_row, col)
159 pub fn row_reduction_with<F>(&mut self, normalize: bool, mut f: F)
160 where
161 F: FnMut(usize, usize, usize),
162 {
163 let (n, m) = self.shape;
164 let mut c = 0;
165 for r in 0..n {
166 loop {
167 if c >= m {
168 return;
169 }
170 if let Some(pivot) = (r..n).find(|&p| !R::is_zero(&self[p][c])) {
171 f(r, pivot, c);
172 self.data.swap(r, pivot);
173 break;
174 };
175 c += 1;
176 }
177 let d = R::inv(&self[r][c]);
178 if normalize {
179 for j in c..m {
180 R::mul_assign(&mut self[r][j], &d);
181 }
182 }
183 for i in (0..n).filter(|&i| i != r) {
184 let mut e = self[i][c].clone();
185 if !normalize {
186 R::mul_assign(&mut e, &d);
187 }
188 for j in c..m {
189 let e = R::mul(&e, &self[r][j]);
190 R::sub_assign(&mut self[i][j], &e);
191 }
192 }
193 c += 1;
194 }
195 }
196
197 pub fn row_reduction(&mut self, normalize: bool) {
198 self.row_reduction_with(normalize, |_, _, _| {});
199 }
200
201 pub fn rank(&mut self) -> usize {
202 let n = self.shape.0;
203 self.row_reduction(false);
204 (0..n)
205 .filter(|&i| !self.data[i].iter().all(|x| R::is_zero(x)))
206 .count()
207 }
208
209 pub fn determinant(&mut self) -> R::T {
210 assert_eq!(self.shape.0, self.shape.1);
211 let mut neg = false;
212 self.row_reduction_with(false, |r, p, _| neg ^= r != p);
213 let mut d = R::one();
214 if neg {
215 d = R::neg(&d);
216 }
217 for i in 0..self.shape.0 {
218 R::mul_assign(&mut d, &self[i][i]);
219 }
220 d
221 }
222
223 pub fn solve_system_of_linear_equations(
224 &self,
225 b: &[R::T],
226 ) -> Option<SystemOfLinearEquationsSolution<R>> {
227 assert_eq!(self.shape.0, b.len());
228 let (n, m) = self.shape;
229 let mut c = Matrix::<R>::zeros((n, m + 1));
230 for i in 0..n {
231 c[i][..m].clone_from_slice(&self[i]);
232 c[i][m] = b[i].clone();
233 }
234 let mut reduced = vec![!0; m + 1];
235 c.row_reduction_with(true, |r, _, c| reduced[c] = r);
236 if reduced[m] != !0 {
237 return None;
238 }
239 let mut particular = vec![R::zero(); m];
240 let mut basis = vec![];
241 for j in 0..m {
242 if reduced[j] != !0 {
243 particular[j] = c[reduced[j]][m].clone();
244 } else {
245 let mut v = vec![R::zero(); m];
246 v[j] = R::one();
247 for i in 0..m {
248 if reduced[i] != !0 {
249 R::sub_assign(&mut v[i], &c[reduced[i]][j]);
250 }
251 }
252 basis.push(v);
253 }
254 }
255 Some(SystemOfLinearEquationsSolution { particular, basis })
256 }
257
258 pub fn inverse(&self) -> Option<Matrix<R>> {
259 assert_eq!(self.shape.0, self.shape.1);
260 let n = self.shape.0;
261 let mut c = Matrix::<R>::zeros((n, n * 2));
262 for i in 0..n {
263 c[i][..n].clone_from_slice(&self[i]);
264 c[i][n + i] = R::one();
265 }
266 c.row_reduction(true);
267 if (0..n).any(|i| R::is_zero(&c[i][i])) {
268 None
269 } else {
270 Some(Self::from_vec(
271 c.data.into_iter().map(|r| r[n..].to_vec()).collect(),
272 ))
273 }
274 }
275
276 pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
277 let n = self.shape.0;
278 if n == 0 {
279 return vec![R::one()];
280 }
281 assert!(self.data.iter().all(|a| a.len() == n));
282 for j in 0..(n - 1) {
283 if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
284 self.data.swap(j + 1, x);
285 self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
286 let inv = R::inv(&self[j + 1][j]);
287 let mut v = vec![];
288 let src = std::mem::take(&mut self[j + 1]);
289 for a in self.data[(j + 2)..].iter_mut() {
290 let mul = R::mul(&a[j], &inv);
291 for (a, src) in a[j..].iter_mut().zip(src[j..].iter()) {
292 R::sub_assign(a, &R::mul(&mul, src));
293 }
294 v.push(mul);
295 }
296 self[j + 1] = src;
297 for a in self.data.iter_mut() {
298 let v = a[(j + 2)..]
299 .iter()
300 .zip(v.iter())
301 .fold(R::zero(), |s, a| R::add(&s, &R::mul(a.0, a.1)));
302 R::add_assign(&mut a[j + 1], &v);
303 }
304 }
305 }
306 let mut dp = vec![vec![R::one()]];
307 for i in 0..n {
308 let mut next = vec![R::zero(); i + 2];
309 for (j, dp) in dp[i].iter().enumerate() {
310 R::sub_assign(&mut next[j], &R::mul(dp, &self[i][i]));
311 R::add_assign(&mut next[j + 1], dp);
312 }
313 let mut mul = R::one();
314 for j in (0..i).rev() {
315 mul = R::mul(&mul, &self[j + 1][j]);
316 let c = R::mul(&mul, &self[j][i]);
317 for (next, dp) in next.iter_mut().zip(dp[j].iter()) {
318 R::sub_assign(next, &R::mul(&c, dp));
319 }
320 }
321 dp.push(next);
322 }
323 dp.pop().unwrap()
324 }crates/competitive/src/math/floor_sum.rs (line 125)
124 fn to_x() -> FloorSumData<R, X, Y> {
125 let mut dp = array![array![R::zero(); Y]; X];
126 dp[0][0] = R::one();
127 FloorSumData {
128 dp,
129 dx: R::one(),
130 dy: R::zero(),
131 _marker: PhantomData,
132 }
133 }
134 fn to_y() -> FloorSumData<R, X, Y> {
135 FloorSumData {
136 dp: array![array![R::zero(); Y]; X],
137 dx: R::zero(),
138 dy: R::one(),
139 _marker: PhantomData,
140 }
141 }
142}
143
144impl<R, const X: usize, const Y: usize> FloorSum<R, X, Y>
145where
146 R: Ring<Additive: Invertible>,
147{
148 fn offset(x: i64, y: i64) -> FloorSumData<R, X, Y> {
149 FloorSumData {
150 dp: array![array![R::zero(); Y]; X],
151 dx: R::Additive::signed_pow(R::one(), x),
152 dy: R::Additive::signed_pow(R::one(), y),
153 _marker: PhantomData,
154 }
155 }
156}
157
158impl<R, const X: usize, const Y: usize> Magma for FloorSum<R, X, Y>
159where
160 R: SemiRing,
161{
162 type T = FloorSumData<R, X, Y>;
163
164 fn operate(a: &Self::T, b: &Self::T) -> Self::T {
165 let mut a = a.clone();
166 let mut b = b.clone();
167 let mut pow_x = array![R::zero(); X];
168 let mut pow_y = array![R::zero(); Y];
169 pow_x[0] = R::one();
170 pow_y[0] = R::one();
171 for i in 1..X {
172 pow_x[i] = R::mul(&pow_x[i - 1], &a.dx);
173 }
174 for j in 1..Y {
175 pow_y[j] = R::mul(&pow_y[j - 1], &a.dy);
176 }
177 macro_rules! go {
178 ($N:ident) => {
179 let mut comb = array![array![R::zero(); $N]; $N];
180 comb[0][0] = R::one();
181 let mut i = 0;
182 while i + 1 < $N {
183 let mut j = 0;
184 while j <= i {
185 let x = comb[i][j].clone();
186 R::add_assign(&mut comb[i + 1][j], &x);
187 R::add_assign(&mut comb[i + 1][j + 1], &x);
188 j += 1;
189 }
190 i += 1;
191 }
192 for i in 0..X {
193 for j in (0..Y).rev() {
194 for k in j + 1..Y {
195 let mut x = b.dp[i][j].clone();
196 R::mul_assign(&mut x, &comb[k][j]);
197 R::mul_assign(&mut x, &pow_y[k - j]);
198 R::add_assign(&mut b.dp[i][k], &x);
199 }
200 }
201 }
202 for j in 0..Y {
203 for i in (0..X).rev() {
204 for k in i..X {
205 let mut x = b.dp[i][j].clone();
206 R::mul_assign(&mut x, &comb[k][i]);
207 R::mul_assign(&mut x, &pow_x[k - i]);
208 R::add_assign(&mut a.dp[k][j], &x);
209 }
210 }
211 }
212 };
213 }
214 if X <= Y {
215 go!(Y);
216 } else {
217 go!(X);
218 }
219 R::add_assign(&mut a.dx, &b.dx);
220 R::add_assign(&mut a.dy, &b.dy);
221 a
222 }
223}
224
225impl<R, const X: usize, const Y: usize> Unital for FloorSum<R, X, Y>
226where
227 R: SemiRing,
228{
229 fn unit() -> Self::T {
230 FloorSumData {
231 dp: array![array![R::zero(); Y]; X],
232 dx: R::zero(),
233 dy: R::zero(),
234 _marker: PhantomData,
235 }
236 }
237}
238
239impl<R, const X: usize, const Y: usize> Associative for FloorSum<R, X, Y> where R: SemiRing {}
240
241fn floor_monoid_product<M>(
242 mut x: M::T,
243 mut y: M::T,
244 mut n: u64,
245 mut a: u64,
246 mut b: u64,
247 mut m: u64,
248) -> M::T
249where
250 M: Monoid,
251{
252 let mut c = (a * n + b) / m;
253 let mut pre = M::unit();
254 let mut suf = M::unit();
255 loop {
256 let (p, q) = (a / m, b / m);
257 a %= m;
258 b %= m;
259 x = M::operate(&x, &M::pow(y.clone(), p));
260 pre = M::operate(&pre, &M::pow(y.clone(), q));
261 c -= p * n + q;
262 if c == 0 {
263 break;
264 }
265 let d = (m * c - b - 1) / a + 1;
266 suf = M::operate(&y, &M::operate(&M::pow(x.clone(), n - d), &suf));
267 b = m - b - 1 + a;
268 n = c - 1;
269 c = d;
270 swap(&mut m, &mut a);
271 swap(&mut x, &mut y);
272 }
273 x = M::pow(x.clone(), n);
274 M::operate(&M::operate(&pre, &x), &suf)
275}
276
277/// $$\sum_{i=0}^{n-1}i^X\left\lfloor\frac{a\times i+b}{m}\right\rfloor^Y$$
278pub fn floor_sum_polynomial<T, const X: usize, const Y: usize>(
279 n: u64,
280 a: u64,
281 b: u64,
282 m: u64,
283) -> [[T; Y]; X]
284where
285 T: Clone + Zero + One + Add<Output = T> + Mul<Output = T>,
286{
287 debug_assert!(a == 0 || n < (u64::MAX - b) / a);
288 floor_monoid_product::<FloorSum<AddMulOperation<T>, X, Y>>(
289 FloorSum::<AddMulOperation<T>, X, Y>::to_x(),
290 FloorSum::<AddMulOperation<T>, X, Y>::to_y(),
291 n,
292 a,
293 b,
294 m,
295 )
296 .dp
297}
298
299/// $$\sum_{i=l}^{r-1}i^X\left\lfloor\frac{a\times i+b}{m}\right\rfloor^Y$$
300pub fn floor_sum_polynomial_i64<T, const X: usize, const Y: usize>(
301 l: i64,
302 r: i64,
303 a: i64,
304 b: i64,
305 m: u64,
306) -> [[T; Y]; X]
307where
308 T: Clone + Zero + One + Add<Output = T> + Mul<Output = T>,
309 AddMulOperation<T>: SemiRing<T = T, Additive: Invertible>,
310{
311 assert!(l <= r);
312 assert!(m > 0);
313
314 if a < 0 {
315 let mut ans = floor_sum_polynomial_i64::<T, X, Y>(-r + 1, -l + 1, -a, b, m);
316 for ans in ans.iter_mut().skip(1).step_by(2) {
317 for ans in ans.iter_mut() {
318 *ans = AddMulOperation::<T>::neg(ans);
319 }
320 }
321 return ans;
322 }
323
324 let add_x = l;
325 let n = (r - l) as u64;
326 let b = a * add_x + b;
327
328 let add_y = b.div_euclid(m as i64);
329 let b = b.rem_euclid(m as i64);
330 assert!(a >= 0);
331 assert!(b >= 0);
332 let data = floor_monoid_product::<FloorSum<AddMulOperation<T>, X, Y>>(
333 FloorSum::<AddMulOperation<T>, X, Y>::to_x(),
334 FloorSum::<AddMulOperation<T>, X, Y>::to_y(),
335 n,
336 a as u64,
337 b as u64,
338 m,
339 );
340
341 let offset = FloorSum::<AddMulOperation<T>, X, Y>::offset(add_x, add_y);
342 FloorSum::<AddMulOperation<T>, X, Y>::operate(&offset, &data).dp
343}
344
345#[derive(Debug)]
346struct FloorPowerSum<R>
347where
348 R: SemiRing,
349{
350 x: R::T,
351 sum: R::T,
352}
353
354impl<R> Clone for FloorPowerSum<R>
355where
356 R: SemiRing,
357{
358 fn clone(&self) -> Self {
359 Self {
360 x: self.x.clone(),
361 sum: self.sum.clone(),
362 }
363 }
364}
365
366impl<R> FloorPowerSum<R>
367where
368 R: SemiRing,
369{
370 fn to_x(x: R::T) -> Self {
371 Self { x, sum: R::one() }
372 }
373 fn to_y(y: R::T) -> Self {
374 Self {
375 x: y,
376 sum: R::zero(),
377 }
378 }
379}
380
381impl<R> Magma for FloorPowerSum<R>
382where
383 R: SemiRing,
384{
385 type T = Self;
386
387 fn operate(a: &Self::T, b: &Self::T) -> Self::T {
388 Self {
389 x: R::mul(&a.x, &b.x),
390 sum: R::add(&a.sum, &R::mul(&a.x, &b.sum)),
391 }
392 }
393}
394
395impl<R> Unital for FloorPowerSum<R>
396where
397 R: SemiRing,
398{
399 fn unit() -> Self::T {
400 Self {
401 x: R::one(),
402 sum: R::zero(),
403 }
404 }crates/competitive/src/math/black_box_matrix.rs (line 25)
23 fn apply(&self, v: &[R::T]) -> Vec<R::T> {
24 assert_eq!(self.shape.1, v.len());
25 let mut res = vec![R::zero(); self.shape.0];
26 for i in 0..self.shape.0 {
27 for j in 0..self.shape.1 {
28 R::add_assign(&mut res[i], &R::mul(&self[(i, j)], &v[j]));
29 }
30 }
31 res
32 }
33
34 fn shape(&self) -> (usize, usize) {
35 self.shape
36 }
37}
38
39pub struct SparseMatrix<R>
40where
41 R: SemiRing,
42{
43 shape: (usize, usize),
44 nonzero: Vec<(usize, usize, R::T)>,
45}
46
47impl<R> Debug for SparseMatrix<R>
48where
49 R: SemiRing<T: Debug>,
50{
51 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
52 f.debug_struct("SparseMatrix")
53 .field("shape", &self.shape)
54 .field("nonzero", &self.nonzero)
55 .finish()
56 }
57}
58
59impl<R> Clone for SparseMatrix<R>
60where
61 R: SemiRing,
62{
63 fn clone(&self) -> Self {
64 Self {
65 shape: self.shape,
66 nonzero: self.nonzero.clone(),
67 }
68 }
69}
70
71impl<R> SparseMatrix<R>
72where
73 R: SemiRing,
74{
75 pub fn new(shape: (usize, usize)) -> Self {
76 Self {
77 shape,
78 nonzero: vec![],
79 }
80 }
81 pub fn new_with<F>(shape: (usize, usize), f: F) -> Self
82 where
83 R: SemiRing<T: PartialEq>,
84 F: Fn(usize, usize) -> R::T,
85 {
86 let mut nonzero = vec![];
87 for i in 0..shape.0 {
88 for j in 0..shape.1 {
89 let v = f(i, j);
90 if !R::is_zero(&v) {
91 nonzero.push((i, j, v));
92 }
93 }
94 }
95 Self { shape, nonzero }
96 }
97 pub fn from_nonzero(shape: (usize, usize), nonzero: Vec<(usize, usize, R::T)>) -> Self {
98 Self { shape, nonzero }
99 }
100}
101
102impl<R> From<Matrix<R>> for SparseMatrix<R>
103where
104 R: SemiRing<T: PartialEq>,
105{
106 fn from(mat: Matrix<R>) -> Self {
107 let mut nonzero = vec![];
108 for i in 0..mat.shape.0 {
109 for j in 0..mat.shape.1 {
110 let v = mat[(i, j)].clone();
111 if !R::is_zero(&v) {
112 nonzero.push((i, j, v));
113 }
114 }
115 }
116 Self {
117 shape: mat.shape,
118 nonzero,
119 }
120 }
121}
122
123impl<R> From<SparseMatrix<R>> for Matrix<R>
124where
125 R: SemiRing,
126{
127 fn from(smat: SparseMatrix<R>) -> Self {
128 let mut mat = Matrix::zeros(smat.shape);
129 for &(i, j, ref v) in &smat.nonzero {
130 R::add_assign(&mut mat[(i, j)], v);
131 }
132 mat
133 }
134}
135
136impl<R> BlackBoxMatrix<R> for SparseMatrix<R>
137where
138 R: SemiRing,
139{
140 fn apply(&self, v: &[R::T]) -> Vec<R::T> {
141 assert_eq!(self.shape.1, v.len());
142 let mut res = vec![R::zero(); self.shape.0];
143 for &(i, j, ref val) in &self.nonzero {
144 R::add_assign(&mut res[i], &R::mul(val, &v[j]));
145 }
146 res
147 }crates/competitive/src/math/subset_convolve.rs (line 21)
19 fn transform(t: Self::T, len: usize) -> Self::F {
20 let k = len.trailing_zeros() as usize;
21 let mut f = vec![vec![R::zero(); len]; k + 1];
22 for (i, t) in t.iter().enumerate() {
23 let f = &mut f[i.count_ones() as usize][i];
24 *f = R::add(f, t);
25 }
26 for f in f.iter_mut() {
27 BitwiseorConvolve::<R::Additive>::zeta_transform(f);
28 }
29 f
30 }
31
32 fn inverse_transform(mut f: Self::F, len: usize) -> Self::T {
33 for f in f.iter_mut() {
34 BitwiseorConvolve::<R::Additive>::mobius_transform(f);
35 }
36 let mut t = vec![R::zero(); len];
37 for (i, t) in t.iter_mut().enumerate() {
38 *t = R::add(t, &f[i.count_ones() as usize][i]);
39 }
40 t
41 }Additional examples can be found in:
Sourcefn is_zero(x: &Self::T) -> bool
fn is_zero(x: &Self::T) -> bool
checks if the element is zero
Examples found in repository?
crates/competitive/src/math/black_box_matrix.rs (line 90)
81 pub fn new_with<F>(shape: (usize, usize), f: F) -> Self
82 where
83 R: SemiRing<T: PartialEq>,
84 F: Fn(usize, usize) -> R::T,
85 {
86 let mut nonzero = vec![];
87 for i in 0..shape.0 {
88 for j in 0..shape.1 {
89 let v = f(i, j);
90 if !R::is_zero(&v) {
91 nonzero.push((i, j, v));
92 }
93 }
94 }
95 Self { shape, nonzero }
96 }
97 pub fn from_nonzero(shape: (usize, usize), nonzero: Vec<(usize, usize, R::T)>) -> Self {
98 Self { shape, nonzero }
99 }
100}
101
102impl<R> From<Matrix<R>> for SparseMatrix<R>
103where
104 R: SemiRing<T: PartialEq>,
105{
106 fn from(mat: Matrix<R>) -> Self {
107 let mut nonzero = vec![];
108 for i in 0..mat.shape.0 {
109 for j in 0..mat.shape.1 {
110 let v = mat[(i, j)].clone();
111 if !R::is_zero(&v) {
112 nonzero.push((i, j, v));
113 }
114 }
115 }
116 Self {
117 shape: mat.shape,
118 nonzero,
119 }
120 }More examples
crates/competitive/src/math/matrix.rs (line 170)
159 pub fn row_reduction_with<F>(&mut self, normalize: bool, mut f: F)
160 where
161 F: FnMut(usize, usize, usize),
162 {
163 let (n, m) = self.shape;
164 let mut c = 0;
165 for r in 0..n {
166 loop {
167 if c >= m {
168 return;
169 }
170 if let Some(pivot) = (r..n).find(|&p| !R::is_zero(&self[p][c])) {
171 f(r, pivot, c);
172 self.data.swap(r, pivot);
173 break;
174 };
175 c += 1;
176 }
177 let d = R::inv(&self[r][c]);
178 if normalize {
179 for j in c..m {
180 R::mul_assign(&mut self[r][j], &d);
181 }
182 }
183 for i in (0..n).filter(|&i| i != r) {
184 let mut e = self[i][c].clone();
185 if !normalize {
186 R::mul_assign(&mut e, &d);
187 }
188 for j in c..m {
189 let e = R::mul(&e, &self[r][j]);
190 R::sub_assign(&mut self[i][j], &e);
191 }
192 }
193 c += 1;
194 }
195 }
196
197 pub fn row_reduction(&mut self, normalize: bool) {
198 self.row_reduction_with(normalize, |_, _, _| {});
199 }
200
201 pub fn rank(&mut self) -> usize {
202 let n = self.shape.0;
203 self.row_reduction(false);
204 (0..n)
205 .filter(|&i| !self.data[i].iter().all(|x| R::is_zero(x)))
206 .count()
207 }
208
209 pub fn determinant(&mut self) -> R::T {
210 assert_eq!(self.shape.0, self.shape.1);
211 let mut neg = false;
212 self.row_reduction_with(false, |r, p, _| neg ^= r != p);
213 let mut d = R::one();
214 if neg {
215 d = R::neg(&d);
216 }
217 for i in 0..self.shape.0 {
218 R::mul_assign(&mut d, &self[i][i]);
219 }
220 d
221 }
222
223 pub fn solve_system_of_linear_equations(
224 &self,
225 b: &[R::T],
226 ) -> Option<SystemOfLinearEquationsSolution<R>> {
227 assert_eq!(self.shape.0, b.len());
228 let (n, m) = self.shape;
229 let mut c = Matrix::<R>::zeros((n, m + 1));
230 for i in 0..n {
231 c[i][..m].clone_from_slice(&self[i]);
232 c[i][m] = b[i].clone();
233 }
234 let mut reduced = vec![!0; m + 1];
235 c.row_reduction_with(true, |r, _, c| reduced[c] = r);
236 if reduced[m] != !0 {
237 return None;
238 }
239 let mut particular = vec![R::zero(); m];
240 let mut basis = vec![];
241 for j in 0..m {
242 if reduced[j] != !0 {
243 particular[j] = c[reduced[j]][m].clone();
244 } else {
245 let mut v = vec![R::zero(); m];
246 v[j] = R::one();
247 for i in 0..m {
248 if reduced[i] != !0 {
249 R::sub_assign(&mut v[i], &c[reduced[i]][j]);
250 }
251 }
252 basis.push(v);
253 }
254 }
255 Some(SystemOfLinearEquationsSolution { particular, basis })
256 }
257
258 pub fn inverse(&self) -> Option<Matrix<R>> {
259 assert_eq!(self.shape.0, self.shape.1);
260 let n = self.shape.0;
261 let mut c = Matrix::<R>::zeros((n, n * 2));
262 for i in 0..n {
263 c[i][..n].clone_from_slice(&self[i]);
264 c[i][n + i] = R::one();
265 }
266 c.row_reduction(true);
267 if (0..n).any(|i| R::is_zero(&c[i][i])) {
268 None
269 } else {
270 Some(Self::from_vec(
271 c.data.into_iter().map(|r| r[n..].to_vec()).collect(),
272 ))
273 }
274 }
275
276 pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
277 let n = self.shape.0;
278 if n == 0 {
279 return vec![R::one()];
280 }
281 assert!(self.data.iter().all(|a| a.len() == n));
282 for j in 0..(n - 1) {
283 if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
284 self.data.swap(j + 1, x);
285 self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
286 let inv = R::inv(&self[j + 1][j]);
287 let mut v = vec![];
288 let src = std::mem::take(&mut self[j + 1]);
289 for a in self.data[(j + 2)..].iter_mut() {
290 let mul = R::mul(&a[j], &inv);
291 for (a, src) in a[j..].iter_mut().zip(src[j..].iter()) {
292 R::sub_assign(a, &R::mul(&mul, src));
293 }
294 v.push(mul);
295 }
296 self[j + 1] = src;
297 for a in self.data.iter_mut() {
298 let v = a[(j + 2)..]
299 .iter()
300 .zip(v.iter())
301 .fold(R::zero(), |s, a| R::add(&s, &R::mul(a.0, a.1)));
302 R::add_assign(&mut a[j + 1], &v);
303 }
304 }
305 }
306 let mut dp = vec![vec![R::one()]];
307 for i in 0..n {
308 let mut next = vec![R::zero(); i + 2];
309 for (j, dp) in dp[i].iter().enumerate() {
310 R::sub_assign(&mut next[j], &R::mul(dp, &self[i][i]));
311 R::add_assign(&mut next[j + 1], dp);
312 }
313 let mut mul = R::one();
314 for j in (0..i).rev() {
315 mul = R::mul(&mul, &self[j + 1][j]);
316 let c = R::mul(&mul, &self[j][i]);
317 for (next, dp) in next.iter_mut().zip(dp[j].iter()) {
318 R::sub_assign(next, &R::mul(&c, dp));
319 }
320 }
321 dp.push(next);
322 }
323 dp.pop().unwrap()
324 }crates/competitive/src/algorithm/automata_learning.rs (line 493)
492 fn split_sample(&mut self, sample: &[usize]) -> Option<(Vec<usize>, Vec<usize>)> {
493 if self.prefixes.is_empty() && !F::is_zero(&self.automaton.behavior(sample.iter().cloned()))
494 {
495 return Some((vec![], sample.to_vec()));
496 }
497 let expected = self.automaton.behavior(sample.iter().cloned());
498 if expected == self.wfa.behavior(sample.iter().cloned()) {
499 return None;
500 }
501 let n = sample.len();
502 let dim = self.wfa.final_weights.shape.0;
503 let mut states: Vec<(Matrix<F>, usize)> = Vec::with_capacity(n + 1);
504 let mut v = self.wfa.final_weights.clone();
505 states.push((v.clone(), n));
506 for k in (0..n).rev() {
507 v = &self.wfa.transitions[sample[k]] * &v;
508 states.push((v.clone(), k));
509 }
510 states.reverse();
511 let split = states.partition_point(|(state, k)| {
512 (0..dim).any(|j| {
513 self.automaton.behavior(
514 self.prefixes[j]
515 .iter()
516 .cloned()
517 .chain(sample[*k..].iter().cloned()),
518 ) != state[j][0]
519 })
520 });
521 Some((sample[..split].to_vec(), sample[split..].to_vec()))
522 }
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 i in 0..self.prefixes.len() {
661 h.data[i].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 }Sourcefn one() -> Self::T
fn one() -> Self::T
multiplicative identity: $1$
Examples found in repository?
More examples
crates/competitive/src/string/rolling_hash.rs (line 500)
497 fn new(base: R::T) -> Self {
498 Self {
499 base,
500 pow: vec![R::one()],
501 }
502 }
503 fn ensure_pow(&mut self, len: usize) {
504 if self.pow.len() <= len {
505 self.pow.reserve(len - self.pow.len() + 1);
506 if self.pow.is_empty() {
507 self.pow.push(R::one());
508 }
509 for _ in 0..=len - self.pow.len() {
510 self.pow.push(R::mul(self.pow.last().unwrap(), &self.base));
511 }
512 }
513 }crates/competitive/src/math/floor_sum.rs (line 126)
124 fn to_x() -> FloorSumData<R, X, Y> {
125 let mut dp = array![array![R::zero(); Y]; X];
126 dp[0][0] = R::one();
127 FloorSumData {
128 dp,
129 dx: R::one(),
130 dy: R::zero(),
131 _marker: PhantomData,
132 }
133 }
134 fn to_y() -> FloorSumData<R, X, Y> {
135 FloorSumData {
136 dp: array![array![R::zero(); Y]; X],
137 dx: R::zero(),
138 dy: R::one(),
139 _marker: PhantomData,
140 }
141 }
142}
143
144impl<R, const X: usize, const Y: usize> FloorSum<R, X, Y>
145where
146 R: Ring<Additive: Invertible>,
147{
148 fn offset(x: i64, y: i64) -> FloorSumData<R, X, Y> {
149 FloorSumData {
150 dp: array![array![R::zero(); Y]; X],
151 dx: R::Additive::signed_pow(R::one(), x),
152 dy: R::Additive::signed_pow(R::one(), y),
153 _marker: PhantomData,
154 }
155 }
156}
157
158impl<R, const X: usize, const Y: usize> Magma for FloorSum<R, X, Y>
159where
160 R: SemiRing,
161{
162 type T = FloorSumData<R, X, Y>;
163
164 fn operate(a: &Self::T, b: &Self::T) -> Self::T {
165 let mut a = a.clone();
166 let mut b = b.clone();
167 let mut pow_x = array![R::zero(); X];
168 let mut pow_y = array![R::zero(); Y];
169 pow_x[0] = R::one();
170 pow_y[0] = R::one();
171 for i in 1..X {
172 pow_x[i] = R::mul(&pow_x[i - 1], &a.dx);
173 }
174 for j in 1..Y {
175 pow_y[j] = R::mul(&pow_y[j - 1], &a.dy);
176 }
177 macro_rules! go {
178 ($N:ident) => {
179 let mut comb = array![array![R::zero(); $N]; $N];
180 comb[0][0] = R::one();
181 let mut i = 0;
182 while i + 1 < $N {
183 let mut j = 0;
184 while j <= i {
185 let x = comb[i][j].clone();
186 R::add_assign(&mut comb[i + 1][j], &x);
187 R::add_assign(&mut comb[i + 1][j + 1], &x);
188 j += 1;
189 }
190 i += 1;
191 }
192 for i in 0..X {
193 for j in (0..Y).rev() {
194 for k in j + 1..Y {
195 let mut x = b.dp[i][j].clone();
196 R::mul_assign(&mut x, &comb[k][j]);
197 R::mul_assign(&mut x, &pow_y[k - j]);
198 R::add_assign(&mut b.dp[i][k], &x);
199 }
200 }
201 }
202 for j in 0..Y {
203 for i in (0..X).rev() {
204 for k in i..X {
205 let mut x = b.dp[i][j].clone();
206 R::mul_assign(&mut x, &comb[k][i]);
207 R::mul_assign(&mut x, &pow_x[k - i]);
208 R::add_assign(&mut a.dp[k][j], &x);
209 }
210 }
211 }
212 };
213 }
214 if X <= Y {
215 go!(Y);
216 } else {
217 go!(X);
218 }
219 R::add_assign(&mut a.dx, &b.dx);
220 R::add_assign(&mut a.dy, &b.dy);
221 a
222 }
223}
224
225impl<R, const X: usize, const Y: usize> Unital for FloorSum<R, X, Y>
226where
227 R: SemiRing,
228{
229 fn unit() -> Self::T {
230 FloorSumData {
231 dp: array![array![R::zero(); Y]; X],
232 dx: R::zero(),
233 dy: R::zero(),
234 _marker: PhantomData,
235 }
236 }
237}
238
239impl<R, const X: usize, const Y: usize> Associative for FloorSum<R, X, Y> where R: SemiRing {}
240
241fn floor_monoid_product<M>(
242 mut x: M::T,
243 mut y: M::T,
244 mut n: u64,
245 mut a: u64,
246 mut b: u64,
247 mut m: u64,
248) -> M::T
249where
250 M: Monoid,
251{
252 let mut c = (a * n + b) / m;
253 let mut pre = M::unit();
254 let mut suf = M::unit();
255 loop {
256 let (p, q) = (a / m, b / m);
257 a %= m;
258 b %= m;
259 x = M::operate(&x, &M::pow(y.clone(), p));
260 pre = M::operate(&pre, &M::pow(y.clone(), q));
261 c -= p * n + q;
262 if c == 0 {
263 break;
264 }
265 let d = (m * c - b - 1) / a + 1;
266 suf = M::operate(&y, &M::operate(&M::pow(x.clone(), n - d), &suf));
267 b = m - b - 1 + a;
268 n = c - 1;
269 c = d;
270 swap(&mut m, &mut a);
271 swap(&mut x, &mut y);
272 }
273 x = M::pow(x.clone(), n);
274 M::operate(&M::operate(&pre, &x), &suf)
275}
276
277/// $$\sum_{i=0}^{n-1}i^X\left\lfloor\frac{a\times i+b}{m}\right\rfloor^Y$$
278pub fn floor_sum_polynomial<T, const X: usize, const Y: usize>(
279 n: u64,
280 a: u64,
281 b: u64,
282 m: u64,
283) -> [[T; Y]; X]
284where
285 T: Clone + Zero + One + Add<Output = T> + Mul<Output = T>,
286{
287 debug_assert!(a == 0 || n < (u64::MAX - b) / a);
288 floor_monoid_product::<FloorSum<AddMulOperation<T>, X, Y>>(
289 FloorSum::<AddMulOperation<T>, X, Y>::to_x(),
290 FloorSum::<AddMulOperation<T>, X, Y>::to_y(),
291 n,
292 a,
293 b,
294 m,
295 )
296 .dp
297}
298
299/// $$\sum_{i=l}^{r-1}i^X\left\lfloor\frac{a\times i+b}{m}\right\rfloor^Y$$
300pub fn floor_sum_polynomial_i64<T, const X: usize, const Y: usize>(
301 l: i64,
302 r: i64,
303 a: i64,
304 b: i64,
305 m: u64,
306) -> [[T; Y]; X]
307where
308 T: Clone + Zero + One + Add<Output = T> + Mul<Output = T>,
309 AddMulOperation<T>: SemiRing<T = T, Additive: Invertible>,
310{
311 assert!(l <= r);
312 assert!(m > 0);
313
314 if a < 0 {
315 let mut ans = floor_sum_polynomial_i64::<T, X, Y>(-r + 1, -l + 1, -a, b, m);
316 for ans in ans.iter_mut().skip(1).step_by(2) {
317 for ans in ans.iter_mut() {
318 *ans = AddMulOperation::<T>::neg(ans);
319 }
320 }
321 return ans;
322 }
323
324 let add_x = l;
325 let n = (r - l) as u64;
326 let b = a * add_x + b;
327
328 let add_y = b.div_euclid(m as i64);
329 let b = b.rem_euclid(m as i64);
330 assert!(a >= 0);
331 assert!(b >= 0);
332 let data = floor_monoid_product::<FloorSum<AddMulOperation<T>, X, Y>>(
333 FloorSum::<AddMulOperation<T>, X, Y>::to_x(),
334 FloorSum::<AddMulOperation<T>, X, Y>::to_y(),
335 n,
336 a as u64,
337 b as u64,
338 m,
339 );
340
341 let offset = FloorSum::<AddMulOperation<T>, X, Y>::offset(add_x, add_y);
342 FloorSum::<AddMulOperation<T>, X, Y>::operate(&offset, &data).dp
343}
344
345#[derive(Debug)]
346struct FloorPowerSum<R>
347where
348 R: SemiRing,
349{
350 x: R::T,
351 sum: R::T,
352}
353
354impl<R> Clone for FloorPowerSum<R>
355where
356 R: SemiRing,
357{
358 fn clone(&self) -> Self {
359 Self {
360 x: self.x.clone(),
361 sum: self.sum.clone(),
362 }
363 }
364}
365
366impl<R> FloorPowerSum<R>
367where
368 R: SemiRing,
369{
370 fn to_x(x: R::T) -> Self {
371 Self { x, sum: R::one() }
372 }
373 fn to_y(y: R::T) -> Self {
374 Self {
375 x: y,
376 sum: R::zero(),
377 }
378 }
379}
380
381impl<R> Magma for FloorPowerSum<R>
382where
383 R: SemiRing,
384{
385 type T = Self;
386
387 fn operate(a: &Self::T, b: &Self::T) -> Self::T {
388 Self {
389 x: R::mul(&a.x, &b.x),
390 sum: R::add(&a.sum, &R::mul(&a.x, &b.sum)),
391 }
392 }
393}
394
395impl<R> Unital for FloorPowerSum<R>
396where
397 R: SemiRing,
398{
399 fn unit() -> Self::T {
400 Self {
401 x: R::one(),
402 sum: R::zero(),
403 }
404 }crates/competitive/src/math/matrix.rs (line 98)
95 pub fn eye(shape: (usize, usize)) -> Self {
96 let mut data = vec![vec![R::zero(); shape.1]; shape.0];
97 for (i, d) in data.iter_mut().enumerate().take(shape.1) {
98 d[i] = R::one();
99 }
100 Self {
101 shape,
102 data,
103 _marker: PhantomData,
104 }
105 }
106
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 /// f: (row, pivot_row, col)
159 pub fn row_reduction_with<F>(&mut self, normalize: bool, mut f: F)
160 where
161 F: FnMut(usize, usize, usize),
162 {
163 let (n, m) = self.shape;
164 let mut c = 0;
165 for r in 0..n {
166 loop {
167 if c >= m {
168 return;
169 }
170 if let Some(pivot) = (r..n).find(|&p| !R::is_zero(&self[p][c])) {
171 f(r, pivot, c);
172 self.data.swap(r, pivot);
173 break;
174 };
175 c += 1;
176 }
177 let d = R::inv(&self[r][c]);
178 if normalize {
179 for j in c..m {
180 R::mul_assign(&mut self[r][j], &d);
181 }
182 }
183 for i in (0..n).filter(|&i| i != r) {
184 let mut e = self[i][c].clone();
185 if !normalize {
186 R::mul_assign(&mut e, &d);
187 }
188 for j in c..m {
189 let e = R::mul(&e, &self[r][j]);
190 R::sub_assign(&mut self[i][j], &e);
191 }
192 }
193 c += 1;
194 }
195 }
196
197 pub fn row_reduction(&mut self, normalize: bool) {
198 self.row_reduction_with(normalize, |_, _, _| {});
199 }
200
201 pub fn rank(&mut self) -> usize {
202 let n = self.shape.0;
203 self.row_reduction(false);
204 (0..n)
205 .filter(|&i| !self.data[i].iter().all(|x| R::is_zero(x)))
206 .count()
207 }
208
209 pub fn determinant(&mut self) -> R::T {
210 assert_eq!(self.shape.0, self.shape.1);
211 let mut neg = false;
212 self.row_reduction_with(false, |r, p, _| neg ^= r != p);
213 let mut d = R::one();
214 if neg {
215 d = R::neg(&d);
216 }
217 for i in 0..self.shape.0 {
218 R::mul_assign(&mut d, &self[i][i]);
219 }
220 d
221 }
222
223 pub fn solve_system_of_linear_equations(
224 &self,
225 b: &[R::T],
226 ) -> Option<SystemOfLinearEquationsSolution<R>> {
227 assert_eq!(self.shape.0, b.len());
228 let (n, m) = self.shape;
229 let mut c = Matrix::<R>::zeros((n, m + 1));
230 for i in 0..n {
231 c[i][..m].clone_from_slice(&self[i]);
232 c[i][m] = b[i].clone();
233 }
234 let mut reduced = vec![!0; m + 1];
235 c.row_reduction_with(true, |r, _, c| reduced[c] = r);
236 if reduced[m] != !0 {
237 return None;
238 }
239 let mut particular = vec![R::zero(); m];
240 let mut basis = vec![];
241 for j in 0..m {
242 if reduced[j] != !0 {
243 particular[j] = c[reduced[j]][m].clone();
244 } else {
245 let mut v = vec![R::zero(); m];
246 v[j] = R::one();
247 for i in 0..m {
248 if reduced[i] != !0 {
249 R::sub_assign(&mut v[i], &c[reduced[i]][j]);
250 }
251 }
252 basis.push(v);
253 }
254 }
255 Some(SystemOfLinearEquationsSolution { particular, basis })
256 }
257
258 pub fn inverse(&self) -> Option<Matrix<R>> {
259 assert_eq!(self.shape.0, self.shape.1);
260 let n = self.shape.0;
261 let mut c = Matrix::<R>::zeros((n, n * 2));
262 for i in 0..n {
263 c[i][..n].clone_from_slice(&self[i]);
264 c[i][n + i] = R::one();
265 }
266 c.row_reduction(true);
267 if (0..n).any(|i| R::is_zero(&c[i][i])) {
268 None
269 } else {
270 Some(Self::from_vec(
271 c.data.into_iter().map(|r| r[n..].to_vec()).collect(),
272 ))
273 }
274 }
275
276 pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
277 let n = self.shape.0;
278 if n == 0 {
279 return vec![R::one()];
280 }
281 assert!(self.data.iter().all(|a| a.len() == n));
282 for j in 0..(n - 1) {
283 if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
284 self.data.swap(j + 1, x);
285 self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
286 let inv = R::inv(&self[j + 1][j]);
287 let mut v = vec![];
288 let src = std::mem::take(&mut self[j + 1]);
289 for a in self.data[(j + 2)..].iter_mut() {
290 let mul = R::mul(&a[j], &inv);
291 for (a, src) in a[j..].iter_mut().zip(src[j..].iter()) {
292 R::sub_assign(a, &R::mul(&mul, src));
293 }
294 v.push(mul);
295 }
296 self[j + 1] = src;
297 for a in self.data.iter_mut() {
298 let v = a[(j + 2)..]
299 .iter()
300 .zip(v.iter())
301 .fold(R::zero(), |s, a| R::add(&s, &R::mul(a.0, a.1)));
302 R::add_assign(&mut a[j + 1], &v);
303 }
304 }
305 }
306 let mut dp = vec![vec![R::one()]];
307 for i in 0..n {
308 let mut next = vec![R::zero(); i + 2];
309 for (j, dp) in dp[i].iter().enumerate() {
310 R::sub_assign(&mut next[j], &R::mul(dp, &self[i][i]));
311 R::add_assign(&mut next[j + 1], dp);
312 }
313 let mut mul = R::one();
314 for j in (0..i).rev() {
315 mul = R::mul(&mul, &self[j + 1][j]);
316 let c = R::mul(&mul, &self[j][i]);
317 for (next, dp) in next.iter_mut().zip(dp[j].iter()) {
318 R::sub_assign(next, &R::mul(&c, dp));
319 }
320 }
321 dp.push(next);
322 }
323 dp.pop().unwrap()
324 }crates/competitive/src/algorithm/automata_learning.rs (line 601)
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 i in 0..self.prefixes.len() {
661 h.data[i].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 }Sourcefn add(x: &Self::T, y: &Self::T) -> Self::T
fn add(x: &Self::T, y: &Self::T) -> Self::T
additive operaion: $+$
Examples found in repository?
More examples
crates/competitive/src/graph/shortest_path.rs (line 98)
95 fn add_assign(x: &mut Self::T, y: &Self::T) -> bool {
96 match x.0.cmp(&y.0) {
97 Ordering::Equal => {
98 x.1 = S::add(&x.1, &y.1);
99 false
100 }
101 Ordering::Greater => {
102 *x = y.clone();
103 true
104 }
105 _ => false,
106 }
107 }crates/competitive/src/math/subset_convolve.rs (line 24)
19 fn transform(t: Self::T, len: usize) -> Self::F {
20 let k = len.trailing_zeros() as usize;
21 let mut f = vec![vec![R::zero(); len]; k + 1];
22 for (i, t) in t.iter().enumerate() {
23 let f = &mut f[i.count_ones() as usize][i];
24 *f = R::add(f, t);
25 }
26 for f in f.iter_mut() {
27 BitwiseorConvolve::<R::Additive>::zeta_transform(f);
28 }
29 f
30 }
31
32 fn inverse_transform(mut f: Self::F, len: usize) -> Self::T {
33 for f in f.iter_mut() {
34 BitwiseorConvolve::<R::Additive>::mobius_transform(f);
35 }
36 let mut t = vec![R::zero(); len];
37 for (i, t) in t.iter_mut().enumerate() {
38 *t = R::add(t, &f[i.count_ones() as usize][i]);
39 }
40 t
41 }
42
43 fn multiply(f: &mut Self::F, g: &Self::F) {
44 for i in (0..f.len()).rev() {
45 let (lf, rf) = f.split_at_mut(i);
46 for (x, y) in rf[0].iter_mut().zip(&g[0]) {
47 *x = R::mul(x, y);
48 }
49 for (x, y) in lf.iter().rev().zip(g.iter().skip(1)) {
50 for ((x, y), z) in x.iter().zip(y).zip(&mut rf[0]) {
51 *z = R::add(z, &R::mul(x, y));
52 }
53 }
54 }
55 }crates/competitive/src/algorithm/esper.rs (line 152)
140 pub fn solve(&self, input: Input) -> R::T {
141 let coeff = self
142 .data
143 .get(&(self.class)(&input))
144 .expect("unrecognized class")
145 .as_ref()
146 .expect("failed to solve");
147 let feature = (self.feature)(&input);
148 feature
149 .into_iter()
150 .zip(coeff)
151 .map(|(x, y)| R::mul(&x, y))
152 .fold(R::zero(), |x, y| R::add(&x, &y))
153 }crates/competitive/src/math/quotient_array.rs (line 99)
82 pub fn min_25_sieve<R>(&self, mut f: impl FnMut(u64, u32) -> T) -> Self
83 where
84 T: Clone + One,
85 R: Ring<T = T, Additive: Invertible>,
86 {
87 let mut dp = self.clone();
88 with_prime_list(self.isqrtn, |pl| {
89 for &p in pl.primes_lte(self.isqrtn).iter().rev() {
90 let k = self.quotient_index(p);
91 for (i, q) in Self::index_iter(self.n, self.isqrtn).enumerate() {
92 let mut pc = p;
93 if pc * p > q {
94 break;
95 }
96 let mut c = 1;
97 while q / p >= pc {
98 let x = R::mul(&f(p, c), &(R::sub(&dp[q / pc], &self.data[k])));
99 let x = R::add(&x, &f(p, c + 1));
100 dp.data[i] = R::add(&dp.data[i], &x);
101 c += 1;
102 pc *= p;
103 }
104 }
105 }
106 });
107 for x in &mut dp.data {
108 *x = R::add(x, &T::one());
109 }
110 dp
111 }Additional examples can be found in:
Sourcefn mul(x: &Self::T, y: &Self::T) -> Self::T
fn mul(x: &Self::T, y: &Self::T) -> Self::T
multiplicative operaion: $+$
Examples found in repository?
More examples
crates/competitive/src/math/bitwiseand_convolve.rs (line 51)
49 fn multiply(f: &mut Self::F, g: &Self::F) {
50 for (f, g) in f.iter_mut().zip(g) {
51 *f = R::mul(f, g);
52 }
53 }
54
55 fn convolve(a: Self::T, b: Self::T) -> Self::T {
56 assert_eq!(a.len(), b.len());
57 let len = a.len();
58 let same = a == b;
59 let mut a = Self::transform(a, len);
60 if same {
61 for a in a.iter_mut() {
62 *a = R::mul(a, a);
63 }
64 } else {
65 let b = Self::transform(b, len);
66 Self::multiply(&mut a, &b);
67 }
68 Self::inverse_transform(a, len)
69 }crates/competitive/src/math/bitwiseor_convolve.rs (line 51)
49 fn multiply(f: &mut Self::F, g: &Self::F) {
50 for (f, g) in f.iter_mut().zip(g) {
51 *f = R::mul(f, g);
52 }
53 }
54
55 fn convolve(a: Self::T, b: Self::T) -> Self::T {
56 assert_eq!(a.len(), b.len());
57 let len = a.len();
58 let same = a == b;
59 let mut a = Self::transform(a, len);
60 if same {
61 for a in a.iter_mut() {
62 *a = R::mul(a, a);
63 }
64 } else {
65 let b = Self::transform(b, len);
66 Self::multiply(&mut a, &b);
67 }
68 Self::inverse_transform(a, len)
69 }crates/competitive/src/math/bitwisexor_convolve.rs (line 48)
46 fn multiply(f: &mut Self::F, g: &Self::F) {
47 for (f, g) in f.iter_mut().zip(g) {
48 *f = R::mul(f, g);
49 }
50 }
51
52 fn convolve(a: Self::T, b: Self::T) -> Self::T {
53 assert_eq!(a.len(), b.len());
54 let len = a.len();
55 let same = a == b;
56 let mut a = Self::transform(a, len);
57 if same {
58 for a in a.iter_mut() {
59 *a = R::mul(a, a);
60 }
61 } else {
62 let b = Self::transform(b, len);
63 Self::multiply(&mut a, &b);
64 }
65 Self::inverse_transform(a, len)
66 }
67}
68
69impl<R> ConvolveSteps for BitwisexorConvolve<R, true>
70where
71 R: Field<T: PartialEq + TryFrom<usize>, Additive: Invertible, Multiplicative: Invertible>,
72 <R::T as TryFrom<usize>>::Error: Debug,
73{
74 type T = Vec<R::T>;
75 type F = Vec<R::T>;
76
77 fn length(t: &Self::T) -> usize {
78 t.len()
79 }
80
81 fn transform(mut t: Self::T, _len: usize) -> Self::F {
82 BitwisexorConvolve::<R::Additive, true>::hadamard_transform(&mut t);
83 t
84 }
85
86 fn inverse_transform(mut f: Self::F, len: usize) -> Self::T {
87 BitwisexorConvolve::<R::Additive, true>::hadamard_transform(&mut f);
88 let len = R::T::try_from(len).unwrap();
89 for f in f.iter_mut() {
90 *f = R::div(f, &len);
91 }
92 f
93 }
94
95 fn multiply(f: &mut Self::F, g: &Self::F) {
96 for (f, g) in f.iter_mut().zip(g) {
97 *f = R::mul(f, g);
98 }
99 }
100
101 fn convolve(a: Self::T, b: Self::T) -> Self::T {
102 assert_eq!(a.len(), b.len());
103 let len = a.len();
104 let same = a == b;
105 let mut a = Self::transform(a, len);
106 if same {
107 for a in a.iter_mut() {
108 *a = R::mul(a, a);
109 }
110 } else {
111 let b = Self::transform(b, len);
112 Self::multiply(&mut a, &b);
113 }
114 Self::inverse_transform(a, len)
115 }Additional examples can be found in:
- crates/competitive/src/math/black_box_matrix.rs
- crates/competitive/src/string/rolling_hash.rs
- crates/competitive/src/algorithm/esper.rs
- crates/competitive/src/math/subset_convolve.rs
- crates/competitive/src/math/quotient_array.rs
- crates/competitive/src/math/matrix.rs
- crates/competitive/src/math/floor_sum.rs
- crates/competitive/src/algorithm/automata_learning.rs
Sourcefn add_assign(x: &mut Self::T, y: &Self::T)
fn add_assign(x: &mut Self::T, y: &Self::T)
Examples found in repository?
crates/competitive/src/math/black_box_matrix.rs (line 28)
23 fn apply(&self, v: &[R::T]) -> Vec<R::T> {
24 assert_eq!(self.shape.1, v.len());
25 let mut res = vec![R::zero(); self.shape.0];
26 for i in 0..self.shape.0 {
27 for j in 0..self.shape.1 {
28 R::add_assign(&mut res[i], &R::mul(&self[(i, j)], &v[j]));
29 }
30 }
31 res
32 }
33
34 fn shape(&self) -> (usize, usize) {
35 self.shape
36 }
37}
38
39pub struct SparseMatrix<R>
40where
41 R: SemiRing,
42{
43 shape: (usize, usize),
44 nonzero: Vec<(usize, usize, R::T)>,
45}
46
47impl<R> Debug for SparseMatrix<R>
48where
49 R: SemiRing<T: Debug>,
50{
51 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
52 f.debug_struct("SparseMatrix")
53 .field("shape", &self.shape)
54 .field("nonzero", &self.nonzero)
55 .finish()
56 }
57}
58
59impl<R> Clone for SparseMatrix<R>
60where
61 R: SemiRing,
62{
63 fn clone(&self) -> Self {
64 Self {
65 shape: self.shape,
66 nonzero: self.nonzero.clone(),
67 }
68 }
69}
70
71impl<R> SparseMatrix<R>
72where
73 R: SemiRing,
74{
75 pub fn new(shape: (usize, usize)) -> Self {
76 Self {
77 shape,
78 nonzero: vec![],
79 }
80 }
81 pub fn new_with<F>(shape: (usize, usize), f: F) -> Self
82 where
83 R: SemiRing<T: PartialEq>,
84 F: Fn(usize, usize) -> R::T,
85 {
86 let mut nonzero = vec![];
87 for i in 0..shape.0 {
88 for j in 0..shape.1 {
89 let v = f(i, j);
90 if !R::is_zero(&v) {
91 nonzero.push((i, j, v));
92 }
93 }
94 }
95 Self { shape, nonzero }
96 }
97 pub fn from_nonzero(shape: (usize, usize), nonzero: Vec<(usize, usize, R::T)>) -> Self {
98 Self { shape, nonzero }
99 }
100}
101
102impl<R> From<Matrix<R>> for SparseMatrix<R>
103where
104 R: SemiRing<T: PartialEq>,
105{
106 fn from(mat: Matrix<R>) -> Self {
107 let mut nonzero = vec![];
108 for i in 0..mat.shape.0 {
109 for j in 0..mat.shape.1 {
110 let v = mat[(i, j)].clone();
111 if !R::is_zero(&v) {
112 nonzero.push((i, j, v));
113 }
114 }
115 }
116 Self {
117 shape: mat.shape,
118 nonzero,
119 }
120 }
121}
122
123impl<R> From<SparseMatrix<R>> for Matrix<R>
124where
125 R: SemiRing,
126{
127 fn from(smat: SparseMatrix<R>) -> Self {
128 let mut mat = Matrix::zeros(smat.shape);
129 for &(i, j, ref v) in &smat.nonzero {
130 R::add_assign(&mut mat[(i, j)], v);
131 }
132 mat
133 }
134}
135
136impl<R> BlackBoxMatrix<R> for SparseMatrix<R>
137where
138 R: SemiRing,
139{
140 fn apply(&self, v: &[R::T]) -> Vec<R::T> {
141 assert_eq!(self.shape.1, v.len());
142 let mut res = vec![R::zero(); self.shape.0];
143 for &(i, j, ref val) in &self.nonzero {
144 R::add_assign(&mut res[i], &R::mul(val, &v[j]));
145 }
146 res
147 }More examples
crates/competitive/src/math/matrix.rs (line 302)
276 pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
277 let n = self.shape.0;
278 if n == 0 {
279 return vec![R::one()];
280 }
281 assert!(self.data.iter().all(|a| a.len() == n));
282 for j in 0..(n - 1) {
283 if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
284 self.data.swap(j + 1, x);
285 self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
286 let inv = R::inv(&self[j + 1][j]);
287 let mut v = vec![];
288 let src = std::mem::take(&mut self[j + 1]);
289 for a in self.data[(j + 2)..].iter_mut() {
290 let mul = R::mul(&a[j], &inv);
291 for (a, src) in a[j..].iter_mut().zip(src[j..].iter()) {
292 R::sub_assign(a, &R::mul(&mul, src));
293 }
294 v.push(mul);
295 }
296 self[j + 1] = src;
297 for a in self.data.iter_mut() {
298 let v = a[(j + 2)..]
299 .iter()
300 .zip(v.iter())
301 .fold(R::zero(), |s, a| R::add(&s, &R::mul(a.0, a.1)));
302 R::add_assign(&mut a[j + 1], &v);
303 }
304 }
305 }
306 let mut dp = vec![vec![R::one()]];
307 for i in 0..n {
308 let mut next = vec![R::zero(); i + 2];
309 for (j, dp) in dp[i].iter().enumerate() {
310 R::sub_assign(&mut next[j], &R::mul(dp, &self[i][i]));
311 R::add_assign(&mut next[j + 1], dp);
312 }
313 let mut mul = R::one();
314 for j in (0..i).rev() {
315 mul = R::mul(&mul, &self[j + 1][j]);
316 let c = R::mul(&mul, &self[j][i]);
317 for (next, dp) in next.iter_mut().zip(dp[j].iter()) {
318 R::sub_assign(next, &R::mul(&c, dp));
319 }
320 }
321 dp.push(next);
322 }
323 dp.pop().unwrap()
324 }
325}
326
327impl<R> Index<usize> for Matrix<R>
328where
329 R: SemiRing,
330{
331 type Output = Vec<R::T>;
332 fn index(&self, index: usize) -> &Self::Output {
333 &self.data[index]
334 }
335}
336
337impl<R> IndexMut<usize> for Matrix<R>
338where
339 R: SemiRing,
340{
341 fn index_mut(&mut self, index: usize) -> &mut Self::Output {
342 &mut self.data[index]
343 }
344}
345
346impl<R> Index<(usize, usize)> for Matrix<R>
347where
348 R: SemiRing,
349{
350 type Output = R::T;
351 fn index(&self, index: (usize, usize)) -> &Self::Output {
352 &self.data[index.0][index.1]
353 }
354}
355
356impl<R> IndexMut<(usize, usize)> for Matrix<R>
357where
358 R: SemiRing,
359{
360 fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
361 &mut self.data[index.0][index.1]
362 }
363}
364
365macro_rules! impl_matrix_pairwise_binop {
366 ($imp:ident, $method:ident, $imp_assign:ident, $method_assign:ident $(where [$($clauses:tt)*])?) => {
367 impl<R> $imp_assign for Matrix<R>
368 where
369 R: SemiRing,
370 $($($clauses)*)?
371 {
372 fn $method_assign(&mut self, rhs: Self) {
373 self.pairwise_assign(&rhs, |a, b| R::$method_assign(a, b));
374 }
375 }
376 impl<R> $imp_assign<&Matrix<R>> for Matrix<R>
377 where
378 R: SemiRing,
379 $($($clauses)*)?
380 {
381 fn $method_assign(&mut self, rhs: &Self) {
382 self.pairwise_assign(rhs, |a, b| R::$method_assign(a, b));
383 }
384 }
385 impl<R> $imp for Matrix<R>
386 where
387 R: SemiRing,
388 $($($clauses)*)?
389 {
390 type Output = Matrix<R>;
391 fn $method(mut self, rhs: Self) -> Self::Output {
392 self.$method_assign(rhs);
393 self
394 }
395 }
396 impl<R> $imp<&Matrix<R>> for Matrix<R>
397 where
398 R: SemiRing,
399 $($($clauses)*)?
400 {
401 type Output = Matrix<R>;
402 fn $method(mut self, rhs: &Self) -> Self::Output {
403 self.$method_assign(rhs);
404 self
405 }
406 }
407 impl<R> $imp<Matrix<R>> for &Matrix<R>
408 where
409 R: SemiRing,
410 $($($clauses)*)?
411 {
412 type Output = Matrix<R>;
413 fn $method(self, mut rhs: Matrix<R>) -> Self::Output {
414 rhs.pairwise_assign(self, |a, b| *a = R::$method(b, a));
415 rhs
416 }
417 }
418 impl<R> $imp<&Matrix<R>> for &Matrix<R>
419 where
420 R: SemiRing,
421 $($($clauses)*)?
422 {
423 type Output = Matrix<R>;
424 fn $method(self, rhs: &Matrix<R>) -> Self::Output {
425 let mut this = self.clone();
426 this.$method_assign(rhs);
427 this
428 }
429 }
430 };
431}
432
433impl_matrix_pairwise_binop!(Add, add, AddAssign, add_assign);
434impl_matrix_pairwise_binop!(Sub, sub, SubAssign, sub_assign where [R: SemiRing<Additive: Invertible>]);
435
436impl<R> Mul for Matrix<R>
437where
438 R: SemiRing,
439{
440 type Output = Matrix<R>;
441 fn mul(self, rhs: Self) -> Self::Output {
442 (&self).mul(&rhs)
443 }
444}
445impl<R> Mul<&Matrix<R>> for Matrix<R>
446where
447 R: SemiRing,
448{
449 type Output = Matrix<R>;
450 fn mul(self, rhs: &Matrix<R>) -> Self::Output {
451 (&self).mul(rhs)
452 }
453}
454impl<R> Mul<Matrix<R>> for &Matrix<R>
455where
456 R: SemiRing,
457{
458 type Output = Matrix<R>;
459 fn mul(self, rhs: Matrix<R>) -> Self::Output {
460 self.mul(&rhs)
461 }
462}
463impl<R> Mul<&Matrix<R>> for &Matrix<R>
464where
465 R: SemiRing,
466{
467 type Output = Matrix<R>;
468 fn mul(self, rhs: &Matrix<R>) -> Self::Output {
469 assert_eq!(self.shape.1, rhs.shape.0);
470 let mut res = Matrix::zeros((self.shape.0, rhs.shape.1));
471 for i in 0..self.shape.0 {
472 for k in 0..self.shape.1 {
473 for j in 0..rhs.shape.1 {
474 R::add_assign(&mut res[i][j], &R::mul(&self[i][k], &rhs[k][j]));
475 }
476 }
477 }
478 res
479 }crates/competitive/src/math/floor_sum.rs (line 219)
164 fn operate(a: &Self::T, b: &Self::T) -> Self::T {
165 let mut a = a.clone();
166 let mut b = b.clone();
167 let mut pow_x = array![R::zero(); X];
168 let mut pow_y = array![R::zero(); Y];
169 pow_x[0] = R::one();
170 pow_y[0] = R::one();
171 for i in 1..X {
172 pow_x[i] = R::mul(&pow_x[i - 1], &a.dx);
173 }
174 for j in 1..Y {
175 pow_y[j] = R::mul(&pow_y[j - 1], &a.dy);
176 }
177 macro_rules! go {
178 ($N:ident) => {
179 let mut comb = array![array![R::zero(); $N]; $N];
180 comb[0][0] = R::one();
181 let mut i = 0;
182 while i + 1 < $N {
183 let mut j = 0;
184 while j <= i {
185 let x = comb[i][j].clone();
186 R::add_assign(&mut comb[i + 1][j], &x);
187 R::add_assign(&mut comb[i + 1][j + 1], &x);
188 j += 1;
189 }
190 i += 1;
191 }
192 for i in 0..X {
193 for j in (0..Y).rev() {
194 for k in j + 1..Y {
195 let mut x = b.dp[i][j].clone();
196 R::mul_assign(&mut x, &comb[k][j]);
197 R::mul_assign(&mut x, &pow_y[k - j]);
198 R::add_assign(&mut b.dp[i][k], &x);
199 }
200 }
201 }
202 for j in 0..Y {
203 for i in (0..X).rev() {
204 for k in i..X {
205 let mut x = b.dp[i][j].clone();
206 R::mul_assign(&mut x, &comb[k][i]);
207 R::mul_assign(&mut x, &pow_x[k - i]);
208 R::add_assign(&mut a.dp[k][j], &x);
209 }
210 }
211 }
212 };
213 }
214 if X <= Y {
215 go!(Y);
216 } else {
217 go!(X);
218 }
219 R::add_assign(&mut a.dx, &b.dx);
220 R::add_assign(&mut a.dy, &b.dy);
221 a
222 }crates/competitive/src/algorithm/automata_learning.rs (line 558)
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 }Sourcefn mul_assign(x: &mut Self::T, y: &Self::T)
fn mul_assign(x: &mut Self::T, y: &Self::T)
Examples found in repository?
crates/competitive/src/math/matrix.rs (line 180)
159 pub fn row_reduction_with<F>(&mut self, normalize: bool, mut f: F)
160 where
161 F: FnMut(usize, usize, usize),
162 {
163 let (n, m) = self.shape;
164 let mut c = 0;
165 for r in 0..n {
166 loop {
167 if c >= m {
168 return;
169 }
170 if let Some(pivot) = (r..n).find(|&p| !R::is_zero(&self[p][c])) {
171 f(r, pivot, c);
172 self.data.swap(r, pivot);
173 break;
174 };
175 c += 1;
176 }
177 let d = R::inv(&self[r][c]);
178 if normalize {
179 for j in c..m {
180 R::mul_assign(&mut self[r][j], &d);
181 }
182 }
183 for i in (0..n).filter(|&i| i != r) {
184 let mut e = self[i][c].clone();
185 if !normalize {
186 R::mul_assign(&mut e, &d);
187 }
188 for j in c..m {
189 let e = R::mul(&e, &self[r][j]);
190 R::sub_assign(&mut self[i][j], &e);
191 }
192 }
193 c += 1;
194 }
195 }
196
197 pub fn row_reduction(&mut self, normalize: bool) {
198 self.row_reduction_with(normalize, |_, _, _| {});
199 }
200
201 pub fn rank(&mut self) -> usize {
202 let n = self.shape.0;
203 self.row_reduction(false);
204 (0..n)
205 .filter(|&i| !self.data[i].iter().all(|x| R::is_zero(x)))
206 .count()
207 }
208
209 pub fn determinant(&mut self) -> R::T {
210 assert_eq!(self.shape.0, self.shape.1);
211 let mut neg = false;
212 self.row_reduction_with(false, |r, p, _| neg ^= r != p);
213 let mut d = R::one();
214 if neg {
215 d = R::neg(&d);
216 }
217 for i in 0..self.shape.0 {
218 R::mul_assign(&mut d, &self[i][i]);
219 }
220 d
221 }
222
223 pub fn solve_system_of_linear_equations(
224 &self,
225 b: &[R::T],
226 ) -> Option<SystemOfLinearEquationsSolution<R>> {
227 assert_eq!(self.shape.0, b.len());
228 let (n, m) = self.shape;
229 let mut c = Matrix::<R>::zeros((n, m + 1));
230 for i in 0..n {
231 c[i][..m].clone_from_slice(&self[i]);
232 c[i][m] = b[i].clone();
233 }
234 let mut reduced = vec![!0; m + 1];
235 c.row_reduction_with(true, |r, _, c| reduced[c] = r);
236 if reduced[m] != !0 {
237 return None;
238 }
239 let mut particular = vec![R::zero(); m];
240 let mut basis = vec![];
241 for j in 0..m {
242 if reduced[j] != !0 {
243 particular[j] = c[reduced[j]][m].clone();
244 } else {
245 let mut v = vec![R::zero(); m];
246 v[j] = R::one();
247 for i in 0..m {
248 if reduced[i] != !0 {
249 R::sub_assign(&mut v[i], &c[reduced[i]][j]);
250 }
251 }
252 basis.push(v);
253 }
254 }
255 Some(SystemOfLinearEquationsSolution { particular, basis })
256 }
257
258 pub fn inverse(&self) -> Option<Matrix<R>> {
259 assert_eq!(self.shape.0, self.shape.1);
260 let n = self.shape.0;
261 let mut c = Matrix::<R>::zeros((n, n * 2));
262 for i in 0..n {
263 c[i][..n].clone_from_slice(&self[i]);
264 c[i][n + i] = R::one();
265 }
266 c.row_reduction(true);
267 if (0..n).any(|i| R::is_zero(&c[i][i])) {
268 None
269 } else {
270 Some(Self::from_vec(
271 c.data.into_iter().map(|r| r[n..].to_vec()).collect(),
272 ))
273 }
274 }
275
276 pub fn characteristic_polynomial(&mut self) -> Vec<R::T> {
277 let n = self.shape.0;
278 if n == 0 {
279 return vec![R::one()];
280 }
281 assert!(self.data.iter().all(|a| a.len() == n));
282 for j in 0..(n - 1) {
283 if let Some(x) = ((j + 1)..n).find(|&x| !R::is_zero(&self[x][j])) {
284 self.data.swap(j + 1, x);
285 self.data.iter_mut().for_each(|a| a.swap(j + 1, x));
286 let inv = R::inv(&self[j + 1][j]);
287 let mut v = vec![];
288 let src = std::mem::take(&mut self[j + 1]);
289 for a in self.data[(j + 2)..].iter_mut() {
290 let mul = R::mul(&a[j], &inv);
291 for (a, src) in a[j..].iter_mut().zip(src[j..].iter()) {
292 R::sub_assign(a, &R::mul(&mul, src));
293 }
294 v.push(mul);
295 }
296 self[j + 1] = src;
297 for a in self.data.iter_mut() {
298 let v = a[(j + 2)..]
299 .iter()
300 .zip(v.iter())
301 .fold(R::zero(), |s, a| R::add(&s, &R::mul(a.0, a.1)));
302 R::add_assign(&mut a[j + 1], &v);
303 }
304 }
305 }
306 let mut dp = vec![vec![R::one()]];
307 for i in 0..n {
308 let mut next = vec![R::zero(); i + 2];
309 for (j, dp) in dp[i].iter().enumerate() {
310 R::sub_assign(&mut next[j], &R::mul(dp, &self[i][i]));
311 R::add_assign(&mut next[j + 1], dp);
312 }
313 let mut mul = R::one();
314 for j in (0..i).rev() {
315 mul = R::mul(&mul, &self[j + 1][j]);
316 let c = R::mul(&mul, &self[j][i]);
317 for (next, dp) in next.iter_mut().zip(dp[j].iter()) {
318 R::sub_assign(next, &R::mul(&c, dp));
319 }
320 }
321 dp.push(next);
322 }
323 dp.pop().unwrap()
324 }
325}
326
327impl<R> Index<usize> for Matrix<R>
328where
329 R: SemiRing,
330{
331 type Output = Vec<R::T>;
332 fn index(&self, index: usize) -> &Self::Output {
333 &self.data[index]
334 }
335}
336
337impl<R> IndexMut<usize> for Matrix<R>
338where
339 R: SemiRing,
340{
341 fn index_mut(&mut self, index: usize) -> &mut Self::Output {
342 &mut self.data[index]
343 }
344}
345
346impl<R> Index<(usize, usize)> for Matrix<R>
347where
348 R: SemiRing,
349{
350 type Output = R::T;
351 fn index(&self, index: (usize, usize)) -> &Self::Output {
352 &self.data[index.0][index.1]
353 }
354}
355
356impl<R> IndexMut<(usize, usize)> for Matrix<R>
357where
358 R: SemiRing,
359{
360 fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
361 &mut self.data[index.0][index.1]
362 }
363}
364
365macro_rules! impl_matrix_pairwise_binop {
366 ($imp:ident, $method:ident, $imp_assign:ident, $method_assign:ident $(where [$($clauses:tt)*])?) => {
367 impl<R> $imp_assign for Matrix<R>
368 where
369 R: SemiRing,
370 $($($clauses)*)?
371 {
372 fn $method_assign(&mut self, rhs: Self) {
373 self.pairwise_assign(&rhs, |a, b| R::$method_assign(a, b));
374 }
375 }
376 impl<R> $imp_assign<&Matrix<R>> for Matrix<R>
377 where
378 R: SemiRing,
379 $($($clauses)*)?
380 {
381 fn $method_assign(&mut self, rhs: &Self) {
382 self.pairwise_assign(rhs, |a, b| R::$method_assign(a, b));
383 }
384 }
385 impl<R> $imp for Matrix<R>
386 where
387 R: SemiRing,
388 $($($clauses)*)?
389 {
390 type Output = Matrix<R>;
391 fn $method(mut self, rhs: Self) -> Self::Output {
392 self.$method_assign(rhs);
393 self
394 }
395 }
396 impl<R> $imp<&Matrix<R>> for Matrix<R>
397 where
398 R: SemiRing,
399 $($($clauses)*)?
400 {
401 type Output = Matrix<R>;
402 fn $method(mut self, rhs: &Self) -> Self::Output {
403 self.$method_assign(rhs);
404 self
405 }
406 }
407 impl<R> $imp<Matrix<R>> for &Matrix<R>
408 where
409 R: SemiRing,
410 $($($clauses)*)?
411 {
412 type Output = Matrix<R>;
413 fn $method(self, mut rhs: Matrix<R>) -> Self::Output {
414 rhs.pairwise_assign(self, |a, b| *a = R::$method(b, a));
415 rhs
416 }
417 }
418 impl<R> $imp<&Matrix<R>> for &Matrix<R>
419 where
420 R: SemiRing,
421 $($($clauses)*)?
422 {
423 type Output = Matrix<R>;
424 fn $method(self, rhs: &Matrix<R>) -> Self::Output {
425 let mut this = self.clone();
426 this.$method_assign(rhs);
427 this
428 }
429 }
430 };
431}
432
433impl_matrix_pairwise_binop!(Add, add, AddAssign, add_assign);
434impl_matrix_pairwise_binop!(Sub, sub, SubAssign, sub_assign where [R: SemiRing<Additive: Invertible>]);
435
436impl<R> Mul for Matrix<R>
437where
438 R: SemiRing,
439{
440 type Output = Matrix<R>;
441 fn mul(self, rhs: Self) -> Self::Output {
442 (&self).mul(&rhs)
443 }
444}
445impl<R> Mul<&Matrix<R>> for Matrix<R>
446where
447 R: SemiRing,
448{
449 type Output = Matrix<R>;
450 fn mul(self, rhs: &Matrix<R>) -> Self::Output {
451 (&self).mul(rhs)
452 }
453}
454impl<R> Mul<Matrix<R>> for &Matrix<R>
455where
456 R: SemiRing,
457{
458 type Output = Matrix<R>;
459 fn mul(self, rhs: Matrix<R>) -> Self::Output {
460 self.mul(&rhs)
461 }
462}
463impl<R> Mul<&Matrix<R>> for &Matrix<R>
464where
465 R: SemiRing,
466{
467 type Output = Matrix<R>;
468 fn mul(self, rhs: &Matrix<R>) -> Self::Output {
469 assert_eq!(self.shape.1, rhs.shape.0);
470 let mut res = Matrix::zeros((self.shape.0, rhs.shape.1));
471 for i in 0..self.shape.0 {
472 for k in 0..self.shape.1 {
473 for j in 0..rhs.shape.1 {
474 R::add_assign(&mut res[i][j], &R::mul(&self[i][k], &rhs[k][j]));
475 }
476 }
477 }
478 res
479 }
480}
481
482impl<R> MulAssign<&R::T> for Matrix<R>
483where
484 R: SemiRing,
485{
486 fn mul_assign(&mut self, rhs: &R::T) {
487 for i in 0..self.shape.0 {
488 for j in 0..self.shape.1 {
489 R::mul_assign(&mut self[(i, j)], rhs);
490 }
491 }
492 }Dyn Compatibility§
This trait is not dyn compatible.
In older versions of Rust, dyn compatibility was called "object safety", so this trait is not object safe.