Skip to main content

transform_window

Function transform_window 

Source
fn transform_window<T, C>(
    fps: &FormalPowerSeries<T, C>,
    end: isize,
    length: usize,
) -> C::F
Examples found in repository?
crates/competitive/src/math/formal_power_series/berlekamp_massey.rs (line 407)
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}