fn half_gcd<T, C>(
p: &FormalPowerSeries<T, C>,
q: &FormalPowerSeries<T, C>,
k: usize,
length: usize,
) -> (FpsMatrix<T, C>, FrequencyMatrix<C>)Examples found in repository?
crates/competitive/src/math/formal_power_series/berlekamp_massey.rs (line 240)
221 pub fn berlekamp_massey(input: &[T]) -> Self {
222 if input.last().is_none_or(|value| value.is_zero())
223 && input.iter().all(|value| value.is_zero())
224 {
225 return Self::one();
226 }
227 let max_work = if input.len() <= 1536 {
228 usize::MAX
229 } else {
230 input.len().saturating_mul(2)
231 };
232 if let Some(recurrence) = berlekamp_massey_naive(input, max_work) {
233 return Self::from_vec(recurrence);
234 }
235 let n = input.len();
236 let leading_zeros = input.iter().take_while(|value| value.is_zero()).count();
237 let sequence = Self::from_vec(input.to_vec()).trimed();
238 let mut modulus = Self::zeros(n + 1);
239 modulus[n] = T::one();
240 let (matrix, _) = half_gcd(&modulus, &sequence, n / 2, n.max(1).next_power_of_two());
241 let (x, y) = matrix.multiply_vector(&modulus, &sequence);
242 let mut recurrence = if y.length() == 0 {
243 matrix.a01.clone()
244 } else {
245 matrix.a11.clone()
246 };
247 let recurrence_leading_zeros = recurrence
248 .iter()
249 .take_while(|value| value.is_zero())
250 .count();
251 if recurrence_leading_zeros > 0 {
252 let (division, _) = x.div_rem(y.clone());
253 recurrence = add(recurrence * division, matrix.a01);
254 }
255 let inverse = T::one() / &recurrence[0];
256 for value in recurrence.iter_mut() {
257 *value *= &inverse;
258 }
259 let minimum_length = (leading_zeros + 2).max(y.length() + 1);
260 if recurrence.length() < minimum_length {
261 recurrence.resize(minimum_length);
262 }
263 recurrence
264 }
265}
266
267fn degree<T, C>(fps: &FormalPowerSeries<T, C>) -> isize {
268 fps.length() as isize - 1
269}
270
271fn add<T, C>(
272 left: FormalPowerSeries<T, C>,
273 right: FormalPowerSeries<T, C>,
274) -> FormalPowerSeries<T, C>
275where
276 T: FormalPowerSeriesCoefficient,
277{
278 (left + right).trimed()
279}
280
281fn tail<T, C>(fps: &FormalPowerSeries<T, C>, start: isize) -> FormalPowerSeries<T, C>
282where
283 T: FormalPowerSeriesCoefficient,
284{
285 let start = start.max(0) as usize;
286 if start >= fps.length() {
287 FormalPowerSeries::zero()
288 } else {
289 FormalPowerSeries::from_vec(fps.data[start..].to_vec())
290 }
291}
292
293fn coefficient<T, C>(fps: &FormalPowerSeries<T, C>, index: isize) -> T
294where
295 T: FormalPowerSeriesCoefficient,
296{
297 if index < 0 {
298 T::zero()
299 } else {
300 fps.coeff(index as usize)
301 }
302}
303
304fn brute_force<T, C>(
305 mut p: FormalPowerSeries<T, C>,
306 mut q: FormalPowerSeries<T, C>,
307 k: usize,
308) -> FpsMatrix<T, C>
309where
310 T: FormalPowerSeriesCoefficient,
311 C: NttReuse<T = Vec<T>>,
312 C::F: Clone,
313{
314 let threshold = degree(&p) - k as isize;
315 let mut matrix = FpsMatrix::identity();
316 while q.length() as isize > threshold {
317 let q_degree = q.length() - 1;
318 let mut negative_quotient = vec![T::zero(); p.length() - q.length() + 1];
319 let inverse = -T::one() / &q[q_degree];
320 for i in (0..negative_quotient.len()).rev() {
321 negative_quotient[i] = p[i + q_degree].clone() * &inverse;
322 p[i + q_degree] = T::zero();
323 for j in 0..q_degree {
324 let value = negative_quotient[i].clone() * &q[j];
325 p[i + j] += &value;
326 }
327 }
328 matrix.left_multiply_step(&negative_quotient);
329 p.truncate(q_degree);
330 p.trim_tail_zeros();
331 swap(&mut p, &mut q);
332 }
333 matrix
334}
335
336fn reduced_transform<T, C>(fps: &FormalPowerSeries<T, C>, length: usize) -> C::F
337where
338 T: FormalPowerSeriesCoefficient,
339 C: NttReuse<T = Vec<T>>,
340{
341 let mut coefficients = vec![T::zero(); length];
342 for (i, value) in fps.iter().enumerate() {
343 coefficients[i & (length - 1)] += value;
344 }
345 C::transform_ntt(coefficients, length)
346}
347
348fn transform_window<T, C>(fps: &FormalPowerSeries<T, C>, end: isize, length: usize) -> C::F
349where
350 T: FormalPowerSeriesCoefficient,
351 C: NttReuse<T = Vec<T>>,
352{
353 let start = end - length as isize;
354 let coefficients = (start..end).map(|index| coefficient(fps, index)).collect();
355 C::transform_ntt(coefficients, length)
356}
357
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}