fn generate_frobenius_block<M>(
a: &Matrix<AddMulOperation<MInt<M>>>,
v: Vec<MInt<M>>,
rows: &mut Vec<EchelonRow<M>>,
t: &mut Vec<Vec<MInt<M>>>,
) -> Polynomial<M>where
M: MIntDotProduct,Examples found in repository?
crates/competitive/src/math/mint_matrix.rs (line 269)
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}