fn tail<T, C>(
fps: &FormalPowerSeries<T, C>,
start: isize,
) -> FormalPowerSeries<T, C>where
T: FormalPowerSeriesCoefficient,Examples found in repository?
crates/competitive/src/math/formal_power_series/berlekamp_massey.rs (line 380)
358fn half_gcd<T, C>(
359 p: &FormalPowerSeries<T, C>,
360 q: &FormalPowerSeries<T, C>,
361 k: usize,
362 length: usize,
363) -> (FpsMatrix<T, C>, FrequencyMatrix<C>)
364where
365 T: FormalPowerSeriesCoefficient,
366 C: NttReuse<T = Vec<T>>,
367 C::F: Clone,
368{
369 let d = degree(p);
370 if degree(q) < d - k as isize {
371 let matrix = FpsMatrix::identity();
372 let frequency = matrix.transform(length);
373 return (matrix, frequency);
374 }
375 if k == 1 {
376 let matrix = FpsMatrix {
377 a00: FormalPowerSeries::zero(),
378 a01: FormalPowerSeries::one(),
379 a10: FormalPowerSeries::one(),
380 a11: -(tail(p, d - 2) / tail(q, d - 2)),
381 };
382 let frequency = matrix.transform(length);
383 return (matrix, frequency);
384 }
385 if p.length().min(q.length()) <= 32 {
386 let matrix = brute_force(p.clone(), q.clone(), k);
387 let frequency = matrix.transform(length);
388 return (matrix, frequency);
389 }
390
391 let half = length / 2;
392 if k <= half {
393 let (matrix, frequency) = half_gcd(p, q, k, half);
394 let frequency = matrix.extend_transform(frequency, length);
395 return (matrix, frequency);
396 }
397
398 let (matrix, mut matrix_frequency) = half_gcd(
399 &tail(p, d - 2 * half as isize),
400 &tail(q, d - 2 * half as isize),
401 half,
402 length,
403 );
404 let degeneracy = half as isize - degree(&matrix.a11);
405
406 let (p0, q0) = matrix_frequency.apply(
407 &transform_window(p, d - half as isize + degeneracy, length),
408 &transform_window(q, d - half as isize + degeneracy, length),
409 length,
410 );
411 let (p1, q1) = matrix_frequency.apply(
412 &transform_window(p, d - 2 * half as isize, length),
413 &transform_window(q, d - 2 * half as isize, length),
414 length,
415 );
416 let part_length = (half as isize + degeneracy) as usize;
417 let mut p_reduced = p1[length - part_length..].to_vec();
418 p_reduced.extend_from_slice(&p0[length - part_length..]);
419 let mut q_reduced = q1[length - part_length..].to_vec();
420 q_reduced.extend_from_slice(&q0[length - part_length..]);
421 let mut q_reduced: FormalPowerSeries<T, C> = FormalPowerSeries::from_vec(q_reduced).trimed();
422
423 let position = d - half as isize + degeneracy;
424 let mut leading = T::zero();
425 for i in 0..=position {
426 leading += coefficient(p, i) * coefficient(&matrix.a00, position - i)
427 + coefficient(q, i) * coefficient(&matrix.a01, position - i);
428 }
429 p_reduced.push(leading);
430 let mut p_reduced: FormalPowerSeries<T, C> = FormalPowerSeries::from_vec(p_reduced);
431 if degree(&q_reduced) < 3 * half as isize + degeneracy - k as isize {
432 return (matrix, matrix_frequency);
433 }
434
435 let mut remaining = k as isize - degree(&matrix.a11);
436 let mut top_product = matrix.a11.data.last().unwrap().clone();
437 let mut product_degree = degree(&matrix.a11);
438 if degeneracy > 0 {
439 let skip = (2 * half as isize + 2 * degeneracy - (d - half as isize + degeneracy)).max(0);
440 let (division, remainder) = tail(&p_reduced, skip).div_rem(tail(&q_reduced, skip));
441 remaining -= degree(&division);
442 top_product *= -division.data.last().unwrap().clone();
443 product_degree += degree(&division);
444 matrix_frequency = matrix_frequency.left_multiply_step(&division, length);
445 swap(&mut p_reduced, &mut q_reduced);
446 q_reduced = FormalPowerSeries::zeros(skip as usize);
447 q_reduced.data.extend(remainder.data);
448 }
449
450 let start = 3 * half as isize + degeneracy - k as isize - remaining;
451 let (right_matrix, right_frequency) = half_gcd(
452 &tail(&p_reduced, start),
453 &tail(&q_reduced, start),
454 remaining as usize,
455 length,
456 );
457 let product_frequency = right_frequency.multiply(&matrix_frequency);
458 let mut product = product_frequency.clone().inverse_transform(length);
459 product.a00.truncate(k);
460 product.a00.trim_tail_zeros();
461 product.a01.truncate(k);
462 product.a01.trim_tail_zeros();
463 product.a10.truncate(k);
464 product.a10.trim_tail_zeros();
465 product_degree += degree(&right_matrix.a11);
466 if product_degree == length as isize {
467 product.a11.resize(k + 1);
468 let highest = top_product * right_matrix.a11.data.last().unwrap();
469 product.a11[k] = highest.clone();
470 product.a11[0] -= highest;
471 }
472 product.a11.trim_tail_zeros();
473 let product_frequency = if C::MULTIPLE {
474 product.transform(length)
475 } else {
476 product_frequency
477 };
478 (product, product_frequency)
479}