pub trait DotProduct:
Sized
+ Clone
+ Zero
+ Add<Output = Self>
+ Mul<Output = Self> {
// Provided methods
fn try_matrix_product(
_a: &[Vec<Self>],
_b: &[Vec<Self>],
) -> Option<Vec<Vec<Self>>> { ... }
fn add_scaled_assign(x: &mut [Self], y: &[Self], a: &Self) { ... }
fn dot_product(x: &[Self], y: &[Self]) -> Self { ... }
}Provided Methods§
Sourcefn add_scaled_assign(x: &mut [Self], y: &[Self], a: &Self)
fn add_scaled_assign(x: &mut [Self], y: &[Self], a: &Self)
Examples found in repository?
More examples
crates/competitive/src/math/mint_matrix.rs (line 68)
41 fn pow_frobenius(self, k: usize) -> Self
42 where
43 M: MIntConvert<u64>,
44 {
45 assert_eq!(self.shape.0, self.shape.1);
46 let a = self.transpose();
47 let mut rng = Xorshift::new();
48 let f = loop {
49 if let Some(f) = frobenius_decomposition(&a, &mut rng) {
50 break f;
51 }
52 };
53 let fk = f.pow(k);
54 let n = f.t.shape.0;
55 if f.blocks
56 .iter()
57 .map(|p| (p.0.len() - 1).pow(2))
58 .sum::<usize>()
59 * 4
60 <= n * n
61 {
62 let mut ft = Matrix::zeros((n, n));
63 let mut first = 0;
64 for p in &f.blocks {
65 let d = p.0.len() - 1;
66 for i in first..first + d {
67 for j in first..first + d {
68 MInt::add_scaled_assign(&mut ft[i], &f.t[j], &fk[i][j]);
69 }
70 }
71 first += d;
72 }
73 &f.t_inv * &ft
74 } else {
75 &(&f.t_inv * &fk) * &f.t
76 }
77 }
78}
79
80impl<M> Matrix<AddMulOperation<MInt<M>>>
81where
82 M: MIntDotProduct,
83{
84 fn determinant_linear_non_singular(mut self, mut other: Self) -> Option<Vec<MInt<M>>>
85 where
86 M: MIntDotProduct,
87 {
88 let n = self.data.len();
89 let mut f = MInt::one();
90 for d in 0..n {
91 let i = other.data.iter().position(|other| !other[d].is_zero())?;
92 if i != d {
93 self.data.swap(i, d);
94 other.data.swap(i, d);
95 f = -f;
96 }
97 f *= other[d][d];
98 let r = other[d][d].inv();
99 for j in 0..n {
100 self[d][j] *= r;
101 other[d][j] *= r;
102 }
103 assert!(other[d][d].is_one());
104 for i in d + 1..n {
105 let a = other[i][d];
106 for k in 0..n {
107 self[i][k] = self[i][k] - a * self[d][k];
108 other[i][k] = other[i][k] - a * other[d][k];
109 }
110 }
111 for j in d + 1..n {
112 let a = other[d][j];
113 for k in 0..n {
114 self[k][j] = self[k][j] - a * self[k][d];
115 other[k][j] = other[k][j] - a * other[k][d];
116 }
117 }
118 }
119 for s in self.data.iter_mut() {
120 for s in s.iter_mut() {
121 *s = -*s;
122 }
123 }
124 let mut p = self.characteristic_polynomial();
125 for p in p.iter_mut() {
126 *p *= f;
127 }
128 Some(p)
129 }
130}
131
132struct EchelonRow<M>
133where
134 M: MIntDotProduct,
135{
136 pivot: usize,
137 inv: MInt<M>,
138 row: Vec<MInt<M>>,
139}
140
141struct Polynomial<M>(Vec<MInt<M>>)
142where
143 M: MIntDotProduct;
144
145struct FrobeniusDecomposition<M>
146where
147 M: MIntDotProduct,
148{
149 t: Matrix<AddMulOperation<MInt<M>>>,
150 t_inv: Matrix<AddMulOperation<MInt<M>>>,
151 blocks: Vec<Polynomial<M>>,
152}
153
154impl<M> EchelonRow<M>
155where
156 M: MIntDotProduct,
157{
158 fn reduce(&self, row: &mut [MInt<M>]) {
159 let a = -row[self.pivot] * self.inv;
160 if a.is_zero() {
161 return;
162 }
163 let end = self.row.len();
164 MInt::add_scaled_assign(&mut row[self.pivot..end], &self.row[self.pivot..], &a);
165 }
166}
167
168fn generate_frobenius_block<M>(
169 a: &Matrix<AddMulOperation<MInt<M>>>,
170 mut v: Vec<MInt<M>>,
171 rows: &mut Vec<EchelonRow<M>>,
172 t: &mut Vec<Vec<MInt<M>>>,
173) -> Polynomial<M>
174where
175 M: MIntDotProduct,
176{
177 let n = a.shape.0;
178 loop {
179 let mut row = vec![MInt::zero(); n + rows.len() + 1];
180 let (x, c) = row.split_at_mut(n);
181 x.copy_from_slice(&v);
182 c[rows.len()] = MInt::one();
183 for r in rows.iter() {
184 r.reduce(&mut row);
185 }
186 if let Some(pivot) = row[..n].iter().position(|x| !x.is_zero()) {
187 t.push(v);
188 let u = t.last().unwrap();
189 v = a.data.iter().map(|row| MInt::dot_product(u, row)).collect();
190 rows.push(EchelonRow {
191 pivot,
192 inv: row[pivot].inv(),
193 row,
194 });
195 } else {
196 let mut p = row.split_off(n);
197 while p.last().is_some_and(|x| x.is_zero()) {
198 p.pop();
199 }
200 return Polynomial(p);
201 }
202 }
203}
204
205impl<M> Polynomial<M>
206where
207 M: MIntDotProduct,
208{
209 fn exact_div(mut self, rhs: &Self) -> Option<Self> {
210 let mut q = vec![MInt::zero(); self.0.len() - rhs.0.len() + 1];
211 let inv = rhs.0.last().unwrap().inv();
212 for i in (0..q.len()).rev() {
213 q[i] = self.0[i + rhs.0.len() - 1] * inv;
214 MInt::add_scaled_assign(&mut self.0[i..i + rhs.0.len()], &rhs.0, &-q[i]);
215 }
216 self.0.iter().all(|x| x.is_zero()).then_some(Self(q))
217 }
218
219 fn square_mod(&self, p: &Self) -> Self {
220 let d = p.0.len() - 1;
221 let mut c = vec![MInt::zero(); 2 * d - 1];
222 for (i, &x) in self.0.iter().enumerate() {
223 MInt::add_scaled_assign(&mut c[i..2 * i], &self.0[..i], &(x + x));
224 c[2 * i] += x * x;
225 }
226 for i in (d..c.len()).rev() {
227 let x = c[i];
228 MInt::add_scaled_assign(&mut c[i - d..=i], &p.0, &-x);
229 }
230 c.truncate(d);
231 Self(c)
232 }
233
234 fn x_pow_mod(&self, k: usize) -> Self {
235 let d = self.0.len() - 1;
236 if d == 1 {
237 return Self(vec![(-self.0[0]).pow(k)]);
238 }
239 let mut r = Self(vec![MInt::zero(); d]);
240 r.0[0] = MInt::one();
241 for bit in (0..usize::BITS - k.leading_zeros()).rev() {
242 r = r.square_mod(self);
243 if k >> bit & 1 != 0 {
244 let x = r.0[d - 1];
245 for i in (1..d).rev() {
246 r.0[i] = r.0[i - 1] - x * self.0[i];
247 }
248 r.0[0] = -x * self.0[0];
249 }
250 }
251 r
252 }
253}
254
255fn frobenius_decomposition<M>(
256 a: &Matrix<AddMulOperation<MInt<M>>>,
257 rng: &mut Xorshift,
258) -> Option<FrobeniusDecomposition<M>>
259where
260 M: MIntDotProduct + MIntConvert<u64>,
261{
262 let n = a.shape.0;
263 let mut rows = Vec::with_capacity(n);
264 let mut t = Vec::with_capacity(n);
265 let mut blocks: Vec<Polynomial<M>> = Vec::new();
266 while rows.len() < n {
267 let s = rows.len();
268 let v = (0..n).map(|_| MInt::from(rng.rand64())).collect();
269 let c = generate_frobenius_block(a, v, &mut rows, &mut t);
270 if rows.len() == s {
271 continue;
272 }
273 let p = Polynomial(c.0[s..].to_vec());
274 if c.0[..s].iter().any(|x| !x.is_zero()) {
275 let q = c.exact_div(&p)?;
276 let d = rows.len() - s;
277 let mut coefficients = q.0[..s].to_vec();
278 let mut shifts = Vec::with_capacity(d);
279 for _ in 0..d {
280 shifts.push(coefficients.clone());
281 let mut first = 0;
282 for block in &blocks {
283 let len = block.0.len() - 1;
284 let c = &mut coefficients[first..first + len];
285 let last = c[len - 1];
286 for j in (1..len).rev() {
287 c[j] = c[j - 1] - last * block.0[j];
288 }
289 c[0] = -last * block.0[0];
290 first += len;
291 }
292 }
293 let shifts: Matrix<AddMulOperation<MInt<M>>> = Matrix::from_vec(shifts);
294 if d < 32 {
295 let (previous, current) = t.split_at_mut(s);
296 for (shift, row) in shifts.data.iter().zip(current) {
297 for (factor, source) in shift.iter().zip(previous.iter()) {
298 if !factor.is_zero() {
299 MInt::add_scaled_assign(row, source, factor);
300 }
301 }
302 }
303 } else {
304 let previous = Matrix::from_vec(t[..s].to_vec());
305 let correction = &shifts * &previous;
306 for (row, correction) in t[s..].iter_mut().zip(&correction.data) {
307 for (x, &y) in row.iter_mut().zip(correction) {
308 *x += y;
309 }
310 }
311 }
312 for row in &mut rows[s..] {
313 // Keep the reduced vector fixed: T_new += S*T_old gives C_old -= C_new*S.
314 let (previous, current) = row.row[n..].split_at_mut(s);
315 for (&x, shift) in current.iter().zip(&shifts.data) {
316 MInt::add_scaled_assign(previous, shift, &-x);
317 }
318 }
319 }
320 blocks.push(p);
321 }
322
323 let mut t_inv = vec![vec![MInt::zero(); n]; n];
324 for i in (0..n).rev() {
325 let row = &rows[i];
326 let mut c = row.row[n..].to_vec();
327 c.resize(n, MInt::zero());
328 for x in &mut c {
329 *x *= row.inv;
330 }
331 for next in &rows[i + 1..] {
332 let factor = -row.row[next.pivot] * row.inv;
333 if !factor.is_zero() {
334 MInt::add_scaled_assign(&mut c, &t_inv[next.pivot], &factor);
335 }
336 }
337 t_inv[row.pivot] = c;
338 }
339 Some(FrobeniusDecomposition {
340 t: Matrix::from_vec(t),
341 t_inv: Matrix::from_vec(t_inv),
342 blocks,
343 })
344}Sourcefn dot_product(x: &[Self], y: &[Self]) -> Self
fn dot_product(x: &[Self], y: &[Self]) -> Self
Examples found in repository?
More examples
crates/competitive/src/math/black_box_mint_matrix.rs (line 22)
14 fn minimal_polynomial(&self) -> Vec<MInt<M>> {
15 assert_eq!(self.shape().0, self.shape().1);
16 let n = self.shape().0;
17 let mut rng = Xorshift::new();
18 let b: Vec<MInt<M>> = (0..n).map(|_| MInt::from(rng.rand64())).collect();
19 let u: Vec<MInt<M>> = (0..n).map(|_| MInt::from(rng.rand64())).collect();
20 let a: Vec<MInt<M>> = (0..2 * n)
21 .scan(b, |b, _| {
22 let a = MInt::dot_product(b, &u);
23 *b = self.apply(b);
24 Some(a)
25 })
26 .collect();
27 let polynomial: Fps<M> = FormalPowerSeries::berlekamp_massey(&a);
28 let mut p = polynomial.data;
29 p.reverse();
30 p
31 }crates/competitive/src/math/mint_matrix.rs (line 189)
168fn generate_frobenius_block<M>(
169 a: &Matrix<AddMulOperation<MInt<M>>>,
170 mut v: Vec<MInt<M>>,
171 rows: &mut Vec<EchelonRow<M>>,
172 t: &mut Vec<Vec<MInt<M>>>,
173) -> Polynomial<M>
174where
175 M: MIntDotProduct,
176{
177 let n = a.shape.0;
178 loop {
179 let mut row = vec![MInt::zero(); n + rows.len() + 1];
180 let (x, c) = row.split_at_mut(n);
181 x.copy_from_slice(&v);
182 c[rows.len()] = MInt::one();
183 for r in rows.iter() {
184 r.reduce(&mut row);
185 }
186 if let Some(pivot) = row[..n].iter().position(|x| !x.is_zero()) {
187 t.push(v);
188 let u = t.last().unwrap();
189 v = a.data.iter().map(|row| MInt::dot_product(u, row)).collect();
190 rows.push(EchelonRow {
191 pivot,
192 inv: row[pivot].inv(),
193 row,
194 });
195 } else {
196 let mut p = row.split_off(n);
197 while p.last().is_some_and(|x| x.is_zero()) {
198 p.pop();
199 }
200 return Polynomial(p);
201 }
202 }
203}Dyn Compatibility§
This trait is not dyn compatible.
In older versions of Rust, dyn compatibility was called "object safety".