zippel-gcd 0.2.0

Sparse multivariate polynomial GCD by Zippel's modular interpolation.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98

//! Sparse multivariate polynomial GCD by Zippel's modular interpolation.
#![allow(
    clippy::multiple_crate_versions,
    clippy::redundant_pub_crate,
    clippy::cast_possible_truncation,
    clippy::missing_panics_doc,
    clippy::many_single_char_names,
    clippy::similar_names,
    clippy::module_name_repetitions
)]

mod frac;
mod gcd;
mod linzip;
mod pgcd;
mod poly;

use num_bigint::BigInt;
use num_integer::Integer;
use num_rational::BigRational;
use num_traits::One;
use poly::IntPoly;
use polycore::Monomial;

pub use frac::Frac;

type Poly = polycore::Poly<BigRational>;

/// `f = F / d` with `F` integral.
fn integral(f: &Poly) -> (IntPoly, BigInt) {
    let d = f
        .terms
        .iter()
        .fold(BigInt::one(), |d, (_, c)| d.lcm(c.denom()));
    let terms = f
        .terms
        .iter()
        .map(|(m, c)| (m.exps().to_vec(), c.numer() * (&d / c.denom())));
    (IntPoly::new(f.nvars, terms), d)
}

fn rational(h: &IntPoly, d: &BigInt, like: &Poly) -> Poly {
    let terms = h.terms.iter().map(|(e, c)| {
        (
            Monomial::new(e.clone()),
            BigRational::new(c.clone(), d.clone()),
        )
    });
    Poly::new(terms.collect(), like.nvars, like.order.clone())
}

/// Primitive over Z with a positive leading coefficient in `like`'s monomial order.
fn normalized(h: &IntPoly, like: &Poly) -> Poly {
    rational(h, &BigInt::one(), like).primitive()
}

fn gcd_int(f: &Poly, g: &Poly) -> IntPoly {
    assert_eq!(f.nvars, g.nvars, "polynomials must share a ring");
    gcd::gcd_z(&integral(f).0, &integral(g).0)
}

/// The greatest common divisor, as a primitive integer polynomial with positive leading
/// coefficient. `gcd(0, 0) = 0`.
#[must_use]
pub fn gcd(f: &Poly, g: &Poly) -> Poly {
    normalized(&gcd_int(f, g), f)
}

/// `(h, f / h, g / h)` where `h = gcd(f, g)`.
#[must_use]
pub fn cofactors(f: &Poly, g: &Poly) -> (Poly, Poly, Poly) {
    let h = gcd_int(f, g);
    let zero = Poly::zero(f.nvars, f.order.clone());
    if h.is_zero() {
        return (zero.clone(), zero.clone(), zero);
    }
    let h = normalized(&h, f);
    let hi = integral(&h).0;
    let (hf, df) = integral(f);
    let (hg, dg) = integral(g);
    let qf = hf.div_exact(&hi).expect("gcd divides f");
    let qg = hg.div_exact(&hi).expect("gcd divides g");
    (h, rational(&qf, &df, f), rational(&qg, &dg, g))
}

/// The least common multiple, normalized like [`gcd`]. `lcm(f, 0) = 0`.
#[must_use]
pub fn lcm(f: &Poly, g: &Poly) -> Poly {
    let h = gcd_int(f, g);
    if f.is_zero() || g.is_zero() {
        return Poly::zero(f.nvars, f.order.clone());
    }
    let (hg, _) = integral(g);
    normalized(
        &integral(f).0.mul(&hg.div_exact(&h).expect("gcd divides g")),
        f,
    )
}