Struct BarrettReduction 

pub struct BarrettReduction<T> {
    m: T,
    im: T,
}

Fields

m: T

im: T

Implementations

impl<T> BarrettReduction<T>

where T: Barrettable,

pub fn new(m: T) -> Self

Examples found in repository

crates/competitive/src/math/miller_rabin.rs (line 121)

pub fn miller_rabin(n: u64) -> bool {
    miller_rabin_with_br(n, &BarrettReduction::<u128>::new(n as u128))
}

crates/competitive/src/num/mint/mint_basic.rs (line 132)

    pub fn set_mod(m: u32) {
        DYN_MODULUS_U32
            .with(|cell| unsafe { *cell.get() = BarrettReduction::<u64>::new(m as u64) });
    }
    fn rem(x: u64) -> u64 {
        DYN_MODULUS_U32.with(|cell| unsafe { (*cell.get()).rem(x) })
    }
}
impl DynMIntU32 {
    pub fn set_mod(m: u32) {
        DynModuloU32::set_mod(m)
    }
}

thread_local!(static DYN_MODULUS_U64: UnsafeCell<BarrettReduction<u128>> = const { UnsafeCell::new(BarrettReduction::<u128>::new_with_im(1_000_000_007, !0 / 1_000_000_007)) });
impl DynModuloU64 {
    pub fn set_mod(m: u64) {
        DYN_MODULUS_U64
            .with(|cell| unsafe { *cell.get() = BarrettReduction::<u128>::new(m as u128) })
    }

crates/competitive/src/math/primitive_root.rs (line 9)

pub fn primitive_root(p: u64) -> u64 {
    if p == 2 {
        return 1;
    }
    let phi = p - 1;
    let pf = prime_factors(phi);
    let br = BarrettReduction::<u128>::new(p as _);
    (2..)
        .find(|&g| check_primitive_root(g, phi, &br, &pf))
        .unwrap()
}

crates/competitive/src/math/floor_sum.rs (line 26)

pub fn floor_sum(n: u64, a: u64, b: u64, m: u64) -> u64 {
    let mut ans = Wrapping(0u64);
    let (mut n, mut m, mut a, mut b) = (Wrapping(n), m, a, b);
    loop {
        let br = BarrettReduction::<u64>::new(m);
        if a >= m {
            let (q, r) = br.div_rem(a);
            ans += choose2(n) * q;
            a = r;
        }
        if b >= m {
            let (q, r) = br.div_rem(b);
            ans += n * q;
            b = r;
        }
        let y_max = (n * a + b).0;
        if y_max < m {
            break;
        }
        let (q, r) = br.div_rem(y_max);
        n = Wrapping(q);
        b = r;
        swap(&mut m, &mut a);
    }
    ans.0
}

crates/competitive/src/math/discrete_logarithm.rs (line 73)

    fn discrete_logarithm(&mut self, a: u64, b: u64, p: u64) -> Option<(u64, u64)> {
        let lim = ((((p as f64).log2() * (p as f64).log2().log2()).sqrt() / 2.0 + 1.).exp2() * 0.9)
            as u64;
        self.primes.reserve(lim);
        let primes = self.primes.primes_lte(lim);
        while self.br_primes.len() < primes.len() {
            let br = BarrettReduction::<u64>::new(primes[self.br_primes.len()]);
            self.br_primes.push(br);
        }
        let br_primes = &self.br_primes[..primes.len()];
        self.ic
            .entry(p)
            .or_insert_with(|| IndexCalculusWithPrimitiveRoot::new(p, br_primes))
            .discrete_logarithm(a, b, br_primes)
    }
}

const A: [u32; 150] = [
    62, 61, 60, 60, 59, 58, 58, 58, 57, 56, 56, 56, 56, 55, 55, 55, 54, 54, 54, 53, 53, 53, 53, 52,
    52, 52, 52, 52, 52, 51, 50, 50, 50, 50, 49, 49, 49, 48, 48, 48, 48, 48, 47, 47, 47, 47, 47, 47,
    47, 47, 47, 47, 47, 47, 47, 47, 45, 42, 42, 41, 41, 41, 41, 41, 41, 41, 40, 40, 40, 40, 40, 40,
    40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 38, 38, 38, 38, 38, 32, 32, 32, 32,
    32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 31, 31, 31, 31, 31,
    31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 22, 22, 22, 22,
    22, 22, 22, 22, 22, 22,
];

fn factorize_smooth(mut x: u64, row: &mut [u64], br_primes: &[BarrettReduction<u64>]) -> bool {
    for (j, (&br, r)) in br_primes.iter().zip(row).enumerate() {
        *r = 0;
        loop {
            let (div, rem) = br.div_rem(x);
            if rem != 0 {
                break;
            }
            *r += 1;
            x = div;
        }
        if j < 150 && x >= (1u64 << A[j]) {
            break;
        }
    }
    x == 1
}

#[derive(Debug)]
struct QdrtPowPrec {
    br_qdrt: BarrettReduction<u64>,
    p0: Vec<u64>,
    p1: Vec<u64>,
    p2: Vec<u64>,
    p3: Vec<u64>,
}

impl QdrtPowPrec {
    fn new(a: u64, ord: u64, br: &BarrettReduction<u128>) -> Self {
        let qdrt = (ord as f64).powf(0.25).ceil() as u64;
        let br_qdrt = BarrettReduction::<u64>::new(qdrt);
        let mut p0 = Vec::with_capacity(qdrt as usize);
        let mut p1 = Vec::with_capacity(qdrt as usize);
        let mut p2 = Vec::with_capacity(qdrt as usize);
        let mut p3 = Vec::with_capacity(qdrt as usize);
        let mut acc = 1u64;
        for _ in 0..qdrt {
            p0.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        let a = acc;
        acc = 1;
        for _ in 0..qdrt {
            p1.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        let a = acc;
        acc = 1;
        for _ in 0..qdrt {
            p2.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        let a = acc;
        acc = 1;
        for _ in 0..qdrt {
            p3.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        Self {
            br_qdrt,
            p0,
            p1,
            p2,
            p3,
        }
    }
    fn pow(&self, mut k: u64, br: &BarrettReduction<u128>) -> u64 {
        let (a, b) = self.br_qdrt.div_rem(k);
        let mut x = self.p0[b as usize];
        k = a;
        if k > 0 {
            let (a, b) = self.br_qdrt.div_rem(k);
            x = br.rem(x as u128 * self.p1[b as usize] as u128) as u64;
            k = a;
        }
        if k > 0 {
            let (a, b) = self.br_qdrt.div_rem(k);
            x = br.rem(x as u128 * self.p2[b as usize] as u128) as u64;
            k = a;
        }
        if k > 0 {
            let (_, b) = self.br_qdrt.div_rem(k);
            x = br.rem(x as u128 * self.p3[b as usize] as u128) as u64;
        }
        x
    }
}

fn index_calculus_for_primitive_root(
    p: u64,
    ord: u64,
    br_primes: &[BarrettReduction<u64>],
    prec: &QdrtPowPrec,
) -> Vec<u64> {
    let br_ord = BarrettReduction::<u128>::new(ord as u128);
    let mul = |x: u64, y: u64| br_ord.rem(x as u128 * y as u128) as u64;
    let sub = |x: u64, y: u64| if x < y { x + ord - y } else { x - y };

    let pc = br_primes.len();
    let mut mat: Vec<Vec<u64>> = vec![];
    let mut rows: Vec<Vec<u64>> = vec![];

    let mut rng = Xorshift::default();
    let br = BarrettReduction::<u128>::new(p as u128);

    for i in 0..pc {
        for ri in 0usize.. {
            let mut row = vec![0u64; pc + 1];
            let mut kk = rng.rand(ord - 1) + 1;
            let mut gkk = prec.pow(kk, &br);
            let mut k = kk;
            let mut gk = gkk;
            while ri >= rows.len() {
                row[pc] = k;
                if factorize_smooth(gk, &mut row, br_primes) {
                    rows.push(row);
                    break;
                }
                if k + kk < ord {
                    k += kk;
                    gk = br.rem(gk as u128 * gkk as u128) as u64;
                } else {
                    kk = rng.rand(ord - 1) + 1;
                    gkk = prec.pow(kk, &br);
                    k = kk;
                    gk = gkk;
                }
            }
            let row = &mut rows[ri];
            for j in 0..i {
                if row[j] != 0 {
                    let b = mul(modinv(mat[j][j], ord), row[j]);
                    for (r, a) in row[j..].iter_mut().zip(&mat[j][j..]) {
                        *r = sub(*r, mul(*a, b));
                    }
                }
                assert_eq!(row[j], 0);
            }
            if gcd(row[i], ord) == 1 {
                let last = rows.len() - 1;
                rows.swap(ri, last);
                mat.push(rows.pop().unwrap());
                break;
            }
        }
    }
    for i in (0..pc).rev() {
        for j in i + 1..pc {
            mat[i][pc] = sub(mat[i][pc], mul(mat[i][j], mat[j][pc]));
        }
        mat[i][pc] = mul(mat[i][pc], modinv(mat[i][i], ord));
    }
    (0..pc).map(|i| mat[i][pc]).collect()
}

#[derive(Debug)]
struct IndexCalculusWithPrimitiveRoot {
    p: u64,
    ord: u64,
    prec: QdrtPowPrec,
    coeff: Vec<u64>,
}

impl IndexCalculusWithPrimitiveRoot {
    fn new(p: u64, br_primes: &[BarrettReduction<u64>]) -> Self {
        let ord = p - 1;
        let g = primitive_root(p);
        let br = BarrettReduction::<u128>::new(p as u128);
        let prec = QdrtPowPrec::new(g, ord, &br);
        let coeff = index_calculus_for_primitive_root(p, ord, br_primes, &prec);
        Self {
            p,
            ord,
            prec,
            coeff,
        }
    }
    fn index_calculus(&self, a: u64, br_primes: &[BarrettReduction<u64>]) -> Option<u64> {
        let p = self.p;
        let ord = self.ord;
        let br = BarrettReduction::<u128>::new(p as u128);
        let a = br.rem(a as _) as u64;
        if a == 1 {
            return Some(0);
        }
        if p == 2 {
            return None;
        }

        let mut rng = Xorshift::new();
        let mut row = vec![0u64; br_primes.len()];
        let mut kk = rng.rand(ord - 1) + 1;
        let mut gkk = self.prec.pow(kk, &br);
        let mut k = kk;
        let mut gk = br.rem(gkk as u128 * a as u128) as u64;
        loop {
            if factorize_smooth(gk, &mut row, br_primes) {
                let mut res = ord - k;
                for (&c, &r) in self.coeff.iter().zip(&row) {
                    for _ in 0..r {
                        res += c;
                        if res >= ord {
                            res -= ord;
                        }
                    }
                }
                return Some(res);
            }
            if k + kk < ord {
                k += kk;
                gk = br.rem(gk as u128 * gkk as u128) as u64;
            } else {
                kk = rng.rand(ord - 1) + 1;
                gkk = self.prec.pow(kk, &br);
                k = kk;
                gk = br.rem(gkk as u128 * a as u128) as u64;
            }
        }
    }
    fn discrete_logarithm(
        &self,
        a: u64,
        b: u64,
        br_primes: &[BarrettReduction<u64>],
    ) -> Option<(u64, u64)> {
        let p = self.p;
        let ord = self.ord;
        let br = BarrettReduction::<u128>::new(p as u128);
        let a = br.rem(a as _) as u64;
        let b = br.rem(b as _) as u64;
        if a == 0 {
            return if b == 0 { Some((1, 1)) } else { None };
        }
        if b == 0 {
            return None;
        }

        let x = self.index_calculus(a, br_primes)?;
        let y = self.index_calculus(b, br_primes)?;
        solve_linear_congruence(x, y, ord)
    }
}

thread_local!(
    static IC: UnsafeCell<IndexCalculus> = UnsafeCell::new(IndexCalculus::new());
);

pub fn discrete_logarithm_prime_mod(a: u64, b: u64, p: u64) -> Option<u64> {
    IC.with(|ic| unsafe { &mut *ic.get() }.discrete_logarithm(a, b, p))
        .map(|t| t.0)
}

/// a^x ≡ b (mod n), a has order p^e
fn pohlig_hellman_prime_power_order(a: u64, b: u64, n: u64, p: u64, e: u32) -> Option<u64> {
    let br = BarrettReduction::<u128>::new(n as u128);
    let mul = |x: u64, y: u64| br.rem(x as u128 * y as u128) as u64;
    let block_size = (p as f64).sqrt().ceil() as u64;
    let mut baby = HashMap::<u64, u64>::new();
    let g = pow(a, p.pow(e - 1), &br);
    let mut xj = 1;
    for j in 0..block_size {
        baby.entry(xj).or_insert(j);
        xj = mul(xj, g);
    }
    let xi = modinv(xj, n);
    let mut t = 0u64;
    for k in 0..e {
        let mut h = pow(mul(modinv(pow(a, t, &br), n), b), p.pow(e - 1 - k), &br);
        let mut ok = false;
        for i in (0..block_size * block_size).step_by(block_size as usize) {
            if let Some(j) = baby.get(&h) {
                t += (i + j) * p.pow(k);
                ok = true;
                break;
            }
            h = mul(h, xi);
        }
        if !ok {
            return None;
        }
    }
    Some(t)
}

/// a^x ≡ b (mod p^e)
fn discrete_logarithm_prime_power(a: u64, b: u64, p: u64, e: u32) -> Option<(u64, u64)> {
    assert_ne!(p, 0);
    assert_ne!(e, 0);
    let n = p.pow(e);
    assert!(a < n);
    assert!(b < n);
    assert_eq!(gcd(a, p), 1);
    if p == 1 {
        return Some((0, 1));
    }
    if a == 0 {
        return if b == 0 { Some((1, 1)) } else { None };
    }
    if b == 0 {
        return None;
    }
    if e == 1 {
        return IC.with(|ic| unsafe { &mut *ic.get() }.discrete_logarithm(a, b, p));
    }
    let br = BarrettReduction::<u128>::new(n as _);
    if p == 2 {
        if e >= 3 {
            if a % 4 == 1 && b % 4 != 1 {
                return None;
            }
            let aa = if a % 4 == 1 { a } else { n - a };
            let bb = if b % 4 == 1 { b } else { n - b };
            let g = 5;
            let ord = n / 4;
            let x = pohlig_hellman_prime_power_order(g, aa, n, p, e - 2)?;
            let y = pohlig_hellman_prime_power_order(g, bb, n, p, e - 2)?;
            let t = solve_linear_congruence(x, y, ord)?;
            match (a % 4 == 1, b % 4 == 1) {
                (true, true) => Some(t),
                (false, true) if t.0 % 2 == 0 => Some((t.0, lcm(t.1, 2))),
                (false, false) if t.0 % 2 == 1 => Some((t.0, lcm(t.1, 2))),
                (false, false) if a == b => Some((1, lcm(t.1, 2))),
                _ => None,
            }
        } else if a == 1 {
            if b == 1 { Some((0, 1)) } else { None }
        } else {
            assert_eq!(a, 3);
            if b == 1 {
                Some((0, 2))
            } else if b == 3 {
                Some((1, 2))
            } else {
                None
            }
        }
    } else {
        let ord = n - n / p;
        let pf_ord = prime_factors(ord);
        let g = (2..)
            .find(|&g| check_primitive_root(g, ord, &br, &pf_ord))
            .unwrap();
        let mut pf_p = prime_factors(p - 1);
        pf_p.push((p, e - 1));
        let mut abm = vec![];
        for (q, c) in pf_p {
            let m = q.pow(c);
            let d = ord / m;
            let gg = pow(g, d, &br);
            let aa = pow(a, d, &br);
            let bb = pow(b, d, &br);
            let x = pohlig_hellman_prime_power_order(gg, aa, n, q, c)?;
            let y = pohlig_hellman_prime_power_order(gg, bb, n, q, c)?;
            abm.push((x, y, m));
        }
        solve_linear_congruences(abm)
    }
}

crates/competitive/src/math/prime_factors.rs (line 4)

fn find_factor(n: u64) -> Option<u64> {
    let br = BarrettReduction::<u128>::new(n as u128);
    if miller_rabin_with_br(n, &br) {
        return None;
    }
    let m = 1u64 << ((63 - n.leading_zeros()) / 5);
    let sub = |x: u64, y: u64| x.abs_diff(y);
    let mul = |x: u64, y: u64| br.rem(x as u128 * y as u128) as u64;
    for c in 12.. {
        let f = |x: u64| (br.rem(x as u128 * x as u128) + c) as u64;
        let (mut x, mut y, mut r, mut g, mut k, mut ys) = (0, 2, 1, 1, 0, 0);
        while g == 1 {
            x = y;
            for _ in 0..r {
                y = f(y);
            }
            while r > k && g == 1 {
                ys = y;
                let mut q = 1;
                for _ in 0..m.min(r - k) {
                    y = f(y);
                    q = mul(q, sub(x, y));
                }
                g = gcd(q, n);
                k += m;
            }
            r <<= 1;
        }
        if g == n {
            g = 1;
            while g == 1 {
                ys = f(ys);
                g = gcd(sub(x, ys), n);
            }
        }
        if g < n {
            return Some(g);
        }
    }
    unreachable!();
}

Additional examples can be found in:

• crates/competitive/src/math/arbitrary_mod_binomial.rs

pub const fn new_with_im(m: T, im: T) -> Self

pub const fn get_mod(&self) -> T

pub fn div_rem(&self, a: T) -> (T, T)

Examples found in repository

crates/competitive/src/num/barrett_reduction.rs (line 30)

    pub fn div(&self, a: T) -> T {
        self.div_rem(a).0
    }
    pub fn rem(&self, a: T) -> T {
        self.div_rem(a).1
    }

crates/competitive/src/math/discrete_logarithm.rs (line 98)

fn factorize_smooth(mut x: u64, row: &mut [u64], br_primes: &[BarrettReduction<u64>]) -> bool {
    for (j, (&br, r)) in br_primes.iter().zip(row).enumerate() {
        *r = 0;
        loop {
            let (div, rem) = br.div_rem(x);
            if rem != 0 {
                break;
            }
            *r += 1;
            x = div;
        }
        if j < 150 && x >= (1u64 << A[j]) {
            break;
        }
    }
    x == 1
}

#[derive(Debug)]
struct QdrtPowPrec {
    br_qdrt: BarrettReduction<u64>,
    p0: Vec<u64>,
    p1: Vec<u64>,
    p2: Vec<u64>,
    p3: Vec<u64>,
}

impl QdrtPowPrec {
    fn new(a: u64, ord: u64, br: &BarrettReduction<u128>) -> Self {
        let qdrt = (ord as f64).powf(0.25).ceil() as u64;
        let br_qdrt = BarrettReduction::<u64>::new(qdrt);
        let mut p0 = Vec::with_capacity(qdrt as usize);
        let mut p1 = Vec::with_capacity(qdrt as usize);
        let mut p2 = Vec::with_capacity(qdrt as usize);
        let mut p3 = Vec::with_capacity(qdrt as usize);
        let mut acc = 1u64;
        for _ in 0..qdrt {
            p0.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        let a = acc;
        acc = 1;
        for _ in 0..qdrt {
            p1.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        let a = acc;
        acc = 1;
        for _ in 0..qdrt {
            p2.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        let a = acc;
        acc = 1;
        for _ in 0..qdrt {
            p3.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        Self {
            br_qdrt,
            p0,
            p1,
            p2,
            p3,
        }
    }
    fn pow(&self, mut k: u64, br: &BarrettReduction<u128>) -> u64 {
        let (a, b) = self.br_qdrt.div_rem(k);
        let mut x = self.p0[b as usize];
        k = a;
        if k > 0 {
            let (a, b) = self.br_qdrt.div_rem(k);
            x = br.rem(x as u128 * self.p1[b as usize] as u128) as u64;
            k = a;
        }
        if k > 0 {
            let (a, b) = self.br_qdrt.div_rem(k);
            x = br.rem(x as u128 * self.p2[b as usize] as u128) as u64;
            k = a;
        }
        if k > 0 {
            let (_, b) = self.br_qdrt.div_rem(k);
            x = br.rem(x as u128 * self.p3[b as usize] as u128) as u64;
        }
        x
    }

crates/competitive/src/math/floor_sum.rs (line 28)

pub fn floor_sum(n: u64, a: u64, b: u64, m: u64) -> u64 {
    let mut ans = Wrapping(0u64);
    let (mut n, mut m, mut a, mut b) = (Wrapping(n), m, a, b);
    loop {
        let br = BarrettReduction::<u64>::new(m);
        if a >= m {
            let (q, r) = br.div_rem(a);
            ans += choose2(n) * q;
            a = r;
        }
        if b >= m {
            let (q, r) = br.div_rem(b);
            ans += n * q;
            b = r;
        }
        let y_max = (n * a + b).0;
        if y_max < m {
            break;
        }
        let (q, r) = br.div_rem(y_max);
        n = Wrapping(q);
        b = r;
        swap(&mut m, &mut a);
    }
    ans.0
}

crates/competitive/src/math/arbitrary_mod_binomial.rs (line 111)

    fn combination(&self, mut n: u64, mut k: u64) -> u64 {
        if k > n {
            return 0;
        }
        assert!(self.size as u64 >= n.min(self.m - 1));
        if self.m < 1 << 31 {
            let mut res = 1u64;
            if self.e == 1 {
                while n > 0 {
                    let (nn, n0) = self.bp.div_rem(n);
                    let (nk, k0) = self.bp.div_rem(k);
                    if n0 < k0 {
                        return 0;
                    }
                    res = self.bm.rem(res * self.fact[n0 as usize]);
                    res = self.bm.rem(res * self.inv_fact[k0 as usize]);
                    res = self.bm.rem(res * self.inv_fact[(n0 - k0) as usize]);
                    n = nn;
                    k = nk;
                }
            } else {
                let mut r = n - k;
                let mut e0 = 0;
                let mut eq = 0;
                let mut i = 0;
                while n > 0 {
                    res = self.bm.rem(res * self.fact[self.bm.rem(n) as usize]);
                    res = self.bm.rem(res * self.inv_fact[self.bm.rem(k) as usize]);
                    res = self.bm.rem(res * self.inv_fact[self.bm.rem(r) as usize]);
                    n = self.bp.div(n);
                    k = self.bp.div(k);
                    r = self.bp.div(r);
                    let eps = n - k - r;
                    e0 += eps;
                    if e0 >= self.e as u64 {
                        return 0;
                    }
                    i += 1;
                    if i >= self.e {
                        eq ^= eps & 1;
                    }
                }
                if eq == 1 {
                    res = self.bm.rem(res * self.delta);
                }
                res = self
                    .bm
                    .rem(res * pow32(self.p as _, e0 as _, &self.bm) as u64);
            }
            res
        } else {
            let mut res = 1u128;
            if self.e == 1 {
                while n > 0 {
                    let (nn, n0) = self.bp.div_rem(n);
                    let (nk, k0) = self.bp.div_rem(k);
                    if n0 < k0 {
                        return 0;
                    }
                    res = self.bm128.rem(res * self.fact[n0 as usize] as u128);
                    res = self.bm128.rem(res * self.inv_fact[k0 as usize] as u128);
                    res = self
                        .bm128
                        .rem(res * self.inv_fact[(n0 - k0) as usize] as u128);
                    n = nn;
                    k = nk;
                }
            } else {
                let mut r = n - k;
                let mut e0 = 0;
                let mut eq = 0;
                let mut i = 0;
                while n > 0 {
                    res = self
                        .bm128
                        .rem(res * self.fact[self.bm.rem(n) as usize] as u128);
                    res = self
                        .bm128
                        .rem(res * self.inv_fact[self.bm.rem(k) as usize] as u128);
                    res = self
                        .bm128
                        .rem(res * self.inv_fact[self.bm.rem(r) as usize] as u128);
                    n = self.bp.div(n);
                    k = self.bp.div(k);
                    r = self.bp.div(r);
                    let eps = n - k - r;
                    e0 += eps;
                    if e0 >= self.e as u64 {
                        return 0;
                    }
                    i += 1;
                    if i >= self.e {
                        eq ^= eps & 1;
                    }
                }
                if eq == 1 {
                    res = self.bm128.rem(res * self.delta as u128);
                }
                res = self.bm128.rem(res * pow64(self.p, e0, &self.bm128) as u128);
            }
            res as u64
        }
    }

pub fn div(&self, a: T) -> T

Examples found in repository

crates/competitive/src/math/arbitrary_mod_binomial.rs (line 131)

    fn combination(&self, mut n: u64, mut k: u64) -> u64 {
        if k > n {
            return 0;
        }
        assert!(self.size as u64 >= n.min(self.m - 1));
        if self.m < 1 << 31 {
            let mut res = 1u64;
            if self.e == 1 {
                while n > 0 {
                    let (nn, n0) = self.bp.div_rem(n);
                    let (nk, k0) = self.bp.div_rem(k);
                    if n0 < k0 {
                        return 0;
                    }
                    res = self.bm.rem(res * self.fact[n0 as usize]);
                    res = self.bm.rem(res * self.inv_fact[k0 as usize]);
                    res = self.bm.rem(res * self.inv_fact[(n0 - k0) as usize]);
                    n = nn;
                    k = nk;
                }
            } else {
                let mut r = n - k;
                let mut e0 = 0;
                let mut eq = 0;
                let mut i = 0;
                while n > 0 {
                    res = self.bm.rem(res * self.fact[self.bm.rem(n) as usize]);
                    res = self.bm.rem(res * self.inv_fact[self.bm.rem(k) as usize]);
                    res = self.bm.rem(res * self.inv_fact[self.bm.rem(r) as usize]);
                    n = self.bp.div(n);
                    k = self.bp.div(k);
                    r = self.bp.div(r);
                    let eps = n - k - r;
                    e0 += eps;
                    if e0 >= self.e as u64 {
                        return 0;
                    }
                    i += 1;
                    if i >= self.e {
                        eq ^= eps & 1;
                    }
                }
                if eq == 1 {
                    res = self.bm.rem(res * self.delta);
                }
                res = self
                    .bm
                    .rem(res * pow32(self.p as _, e0 as _, &self.bm) as u64);
            }
            res
        } else {
            let mut res = 1u128;
            if self.e == 1 {
                while n > 0 {
                    let (nn, n0) = self.bp.div_rem(n);
                    let (nk, k0) = self.bp.div_rem(k);
                    if n0 < k0 {
                        return 0;
                    }
                    res = self.bm128.rem(res * self.fact[n0 as usize] as u128);
                    res = self.bm128.rem(res * self.inv_fact[k0 as usize] as u128);
                    res = self
                        .bm128
                        .rem(res * self.inv_fact[(n0 - k0) as usize] as u128);
                    n = nn;
                    k = nk;
                }
            } else {
                let mut r = n - k;
                let mut e0 = 0;
                let mut eq = 0;
                let mut i = 0;
                while n > 0 {
                    res = self
                        .bm128
                        .rem(res * self.fact[self.bm.rem(n) as usize] as u128);
                    res = self
                        .bm128
                        .rem(res * self.inv_fact[self.bm.rem(k) as usize] as u128);
                    res = self
                        .bm128
                        .rem(res * self.inv_fact[self.bm.rem(r) as usize] as u128);
                    n = self.bp.div(n);
                    k = self.bp.div(k);
                    r = self.bp.div(r);
                    let eps = n - k - r;
                    e0 += eps;
                    if e0 >= self.e as u64 {
                        return 0;
                    }
                    i += 1;
                    if i >= self.e {
                        eq ^= eps & 1;
                    }
                }
                if eq == 1 {
                    res = self.bm128.rem(res * self.delta as u128);
                }
                res = self.bm128.rem(res * pow64(self.p, e0, &self.bm128) as u128);
            }
            res as u64
        }
    }

pub fn rem(&self, a: T) -> T

Examples found in repository

crates/competitive/src/num/mint/mint_basic.rs (line 135)

    fn rem(x: u64) -> u64 {
        DYN_MODULUS_U32.with(|cell| unsafe { (*cell.get()).rem(x) })
    }
}
impl DynMIntU32 {
    pub fn set_mod(m: u32) {
        DynModuloU32::set_mod(m)
    }
}

thread_local!(static DYN_MODULUS_U64: UnsafeCell<BarrettReduction<u128>> = const { UnsafeCell::new(BarrettReduction::<u128>::new_with_im(1_000_000_007, !0 / 1_000_000_007)) });
impl DynModuloU64 {
    pub fn set_mod(m: u64) {
        DYN_MODULUS_U64
            .with(|cell| unsafe { *cell.get() = BarrettReduction::<u128>::new(m as u128) })
    }
    fn rem(x: u128) -> u128 {
        DYN_MODULUS_U64.with(|cell| unsafe { (*cell.get()).rem(x) })
    }

crates/competitive/src/math/discrete_logarithm.rs (line 12)

fn pow(x: u64, mut y: u64, br: &BarrettReduction<u128>) -> u64 {
    let mut x = x as u128;
    let mut z: u128 = 1;
    while y > 0 {
        if y & 1 == 1 {
            z = br.rem(z * x);
        }
        x = br.rem(x * x);
        y >>= 1;
    }
    z as u64
}

fn solve_linear_congruence(a: u64, b: u64, m: u64) -> Option<(u64, u64)> {
    let g = gcd(a, m);
    if !b.is_multiple_of(g) {
        return None;
    }
    let (a, b, m) = (a / g, b / g, m / g);
    Some(((b as u128 * modinv(a, m) as u128 % m as u128) as _, m))
}

fn solve_linear_congruences<I>(abm: I) -> Option<(u64, u64)>
where
    I: IntoIterator<Item = (u64, u64, u64)>,
{
    let mut x = 0u64;
    let mut m0 = 1u64;
    for (a, b, m) in abm {
        let mut b = b + m - a * x % m;
        if b >= m {
            b -= m;
        }
        let a = a * m0;
        let g = gcd(a, m);
        if !b.is_multiple_of(g) {
            return None;
        }
        let (a, b, m) = (a / g, b / g, m / g);
        x += (b as u128 * modinv(a, m) as u128 % m as u128 * m0 as u128) as u64;
        m0 *= m;
    }
    Some((x, m0))
}

#[derive(Debug)]
struct IndexCalculus {
    primes: PrimeList,
    br_primes: Vec<BarrettReduction<u64>>,
    ic: HashMap<u64, IndexCalculusWithPrimitiveRoot>,
}

impl IndexCalculus {
    fn new() -> Self {
        Self {
            primes: PrimeList::new(2),
            br_primes: Default::default(),
            ic: Default::default(),
        }
    }
    fn discrete_logarithm(&mut self, a: u64, b: u64, p: u64) -> Option<(u64, u64)> {
        let lim = ((((p as f64).log2() * (p as f64).log2().log2()).sqrt() / 2.0 + 1.).exp2() * 0.9)
            as u64;
        self.primes.reserve(lim);
        let primes = self.primes.primes_lte(lim);
        while self.br_primes.len() < primes.len() {
            let br = BarrettReduction::<u64>::new(primes[self.br_primes.len()]);
            self.br_primes.push(br);
        }
        let br_primes = &self.br_primes[..primes.len()];
        self.ic
            .entry(p)
            .or_insert_with(|| IndexCalculusWithPrimitiveRoot::new(p, br_primes))
            .discrete_logarithm(a, b, br_primes)
    }
}

const A: [u32; 150] = [
    62, 61, 60, 60, 59, 58, 58, 58, 57, 56, 56, 56, 56, 55, 55, 55, 54, 54, 54, 53, 53, 53, 53, 52,
    52, 52, 52, 52, 52, 51, 50, 50, 50, 50, 49, 49, 49, 48, 48, 48, 48, 48, 47, 47, 47, 47, 47, 47,
    47, 47, 47, 47, 47, 47, 47, 47, 45, 42, 42, 41, 41, 41, 41, 41, 41, 41, 40, 40, 40, 40, 40, 40,
    40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 38, 38, 38, 38, 38, 32, 32, 32, 32,
    32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 31, 31, 31, 31, 31,
    31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 22, 22, 22, 22,
    22, 22, 22, 22, 22, 22,
];

fn factorize_smooth(mut x: u64, row: &mut [u64], br_primes: &[BarrettReduction<u64>]) -> bool {
    for (j, (&br, r)) in br_primes.iter().zip(row).enumerate() {
        *r = 0;
        loop {
            let (div, rem) = br.div_rem(x);
            if rem != 0 {
                break;
            }
            *r += 1;
            x = div;
        }
        if j < 150 && x >= (1u64 << A[j]) {
            break;
        }
    }
    x == 1
}

#[derive(Debug)]
struct QdrtPowPrec {
    br_qdrt: BarrettReduction<u64>,
    p0: Vec<u64>,
    p1: Vec<u64>,
    p2: Vec<u64>,
    p3: Vec<u64>,
}

impl QdrtPowPrec {
    fn new(a: u64, ord: u64, br: &BarrettReduction<u128>) -> Self {
        let qdrt = (ord as f64).powf(0.25).ceil() as u64;
        let br_qdrt = BarrettReduction::<u64>::new(qdrt);
        let mut p0 = Vec::with_capacity(qdrt as usize);
        let mut p1 = Vec::with_capacity(qdrt as usize);
        let mut p2 = Vec::with_capacity(qdrt as usize);
        let mut p3 = Vec::with_capacity(qdrt as usize);
        let mut acc = 1u64;
        for _ in 0..qdrt {
            p0.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        let a = acc;
        acc = 1;
        for _ in 0..qdrt {
            p1.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        let a = acc;
        acc = 1;
        for _ in 0..qdrt {
            p2.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        let a = acc;
        acc = 1;
        for _ in 0..qdrt {
            p3.push(acc);
            acc = br.rem(acc as u128 * a as u128) as u64;
        }
        Self {
            br_qdrt,
            p0,
            p1,
            p2,
            p3,
        }
    }
    fn pow(&self, mut k: u64, br: &BarrettReduction<u128>) -> u64 {
        let (a, b) = self.br_qdrt.div_rem(k);
        let mut x = self.p0[b as usize];
        k = a;
        if k > 0 {
            let (a, b) = self.br_qdrt.div_rem(k);
            x = br.rem(x as u128 * self.p1[b as usize] as u128) as u64;
            k = a;
        }
        if k > 0 {
            let (a, b) = self.br_qdrt.div_rem(k);
            x = br.rem(x as u128 * self.p2[b as usize] as u128) as u64;
            k = a;
        }
        if k > 0 {
            let (_, b) = self.br_qdrt.div_rem(k);
            x = br.rem(x as u128 * self.p3[b as usize] as u128) as u64;
        }
        x
    }
}

fn index_calculus_for_primitive_root(
    p: u64,
    ord: u64,
    br_primes: &[BarrettReduction<u64>],
    prec: &QdrtPowPrec,
) -> Vec<u64> {
    let br_ord = BarrettReduction::<u128>::new(ord as u128);
    let mul = |x: u64, y: u64| br_ord.rem(x as u128 * y as u128) as u64;
    let sub = |x: u64, y: u64| if x < y { x + ord - y } else { x - y };

    let pc = br_primes.len();
    let mut mat: Vec<Vec<u64>> = vec![];
    let mut rows: Vec<Vec<u64>> = vec![];

    let mut rng = Xorshift::default();
    let br = BarrettReduction::<u128>::new(p as u128);

    for i in 0..pc {
        for ri in 0usize.. {
            let mut row = vec![0u64; pc + 1];
            let mut kk = rng.rand(ord - 1) + 1;
            let mut gkk = prec.pow(kk, &br);
            let mut k = kk;
            let mut gk = gkk;
            while ri >= rows.len() {
                row[pc] = k;
                if factorize_smooth(gk, &mut row, br_primes) {
                    rows.push(row);
                    break;
                }
                if k + kk < ord {
                    k += kk;
                    gk = br.rem(gk as u128 * gkk as u128) as u64;
                } else {
                    kk = rng.rand(ord - 1) + 1;
                    gkk = prec.pow(kk, &br);
                    k = kk;
                    gk = gkk;
                }
            }
            let row = &mut rows[ri];
            for j in 0..i {
                if row[j] != 0 {
                    let b = mul(modinv(mat[j][j], ord), row[j]);
                    for (r, a) in row[j..].iter_mut().zip(&mat[j][j..]) {
                        *r = sub(*r, mul(*a, b));
                    }
                }
                assert_eq!(row[j], 0);
            }
            if gcd(row[i], ord) == 1 {
                let last = rows.len() - 1;
                rows.swap(ri, last);
                mat.push(rows.pop().unwrap());
                break;
            }
        }
    }
    for i in (0..pc).rev() {
        for j in i + 1..pc {
            mat[i][pc] = sub(mat[i][pc], mul(mat[i][j], mat[j][pc]));
        }
        mat[i][pc] = mul(mat[i][pc], modinv(mat[i][i], ord));
    }
    (0..pc).map(|i| mat[i][pc]).collect()
}

#[derive(Debug)]
struct IndexCalculusWithPrimitiveRoot {
    p: u64,
    ord: u64,
    prec: QdrtPowPrec,
    coeff: Vec<u64>,
}

impl IndexCalculusWithPrimitiveRoot {
    fn new(p: u64, br_primes: &[BarrettReduction<u64>]) -> Self {
        let ord = p - 1;
        let g = primitive_root(p);
        let br = BarrettReduction::<u128>::new(p as u128);
        let prec = QdrtPowPrec::new(g, ord, &br);
        let coeff = index_calculus_for_primitive_root(p, ord, br_primes, &prec);
        Self {
            p,
            ord,
            prec,
            coeff,
        }
    }
    fn index_calculus(&self, a: u64, br_primes: &[BarrettReduction<u64>]) -> Option<u64> {
        let p = self.p;
        let ord = self.ord;
        let br = BarrettReduction::<u128>::new(p as u128);
        let a = br.rem(a as _) as u64;
        if a == 1 {
            return Some(0);
        }
        if p == 2 {
            return None;
        }

        let mut rng = Xorshift::new();
        let mut row = vec![0u64; br_primes.len()];
        let mut kk = rng.rand(ord - 1) + 1;
        let mut gkk = self.prec.pow(kk, &br);
        let mut k = kk;
        let mut gk = br.rem(gkk as u128 * a as u128) as u64;
        loop {
            if factorize_smooth(gk, &mut row, br_primes) {
                let mut res = ord - k;
                for (&c, &r) in self.coeff.iter().zip(&row) {
                    for _ in 0..r {
                        res += c;
                        if res >= ord {
                            res -= ord;
                        }
                    }
                }
                return Some(res);
            }
            if k + kk < ord {
                k += kk;
                gk = br.rem(gk as u128 * gkk as u128) as u64;
            } else {
                kk = rng.rand(ord - 1) + 1;
                gkk = self.prec.pow(kk, &br);
                k = kk;
                gk = br.rem(gkk as u128 * a as u128) as u64;
            }
        }
    }
    fn discrete_logarithm(
        &self,
        a: u64,
        b: u64,
        br_primes: &[BarrettReduction<u64>],
    ) -> Option<(u64, u64)> {
        let p = self.p;
        let ord = self.ord;
        let br = BarrettReduction::<u128>::new(p as u128);
        let a = br.rem(a as _) as u64;
        let b = br.rem(b as _) as u64;
        if a == 0 {
            return if b == 0 { Some((1, 1)) } else { None };
        }
        if b == 0 {
            return None;
        }

        let x = self.index_calculus(a, br_primes)?;
        let y = self.index_calculus(b, br_primes)?;
        solve_linear_congruence(x, y, ord)
    }
}

thread_local!(
    static IC: UnsafeCell<IndexCalculus> = UnsafeCell::new(IndexCalculus::new());
);

pub fn discrete_logarithm_prime_mod(a: u64, b: u64, p: u64) -> Option<u64> {
    IC.with(|ic| unsafe { &mut *ic.get() }.discrete_logarithm(a, b, p))
        .map(|t| t.0)
}

/// a^x ≡ b (mod n), a has order p^e
fn pohlig_hellman_prime_power_order(a: u64, b: u64, n: u64, p: u64, e: u32) -> Option<u64> {
    let br = BarrettReduction::<u128>::new(n as u128);
    let mul = |x: u64, y: u64| br.rem(x as u128 * y as u128) as u64;
    let block_size = (p as f64).sqrt().ceil() as u64;
    let mut baby = HashMap::<u64, u64>::new();
    let g = pow(a, p.pow(e - 1), &br);
    let mut xj = 1;
    for j in 0..block_size {
        baby.entry(xj).or_insert(j);
        xj = mul(xj, g);
    }
    let xi = modinv(xj, n);
    let mut t = 0u64;
    for k in 0..e {
        let mut h = pow(mul(modinv(pow(a, t, &br), n), b), p.pow(e - 1 - k), &br);
        let mut ok = false;
        for i in (0..block_size * block_size).step_by(block_size as usize) {
            if let Some(j) = baby.get(&h) {
                t += (i + j) * p.pow(k);
                ok = true;
                break;
            }
            h = mul(h, xi);
        }
        if !ok {
            return None;
        }
    }
    Some(t)
}

crates/competitive/src/math/arbitrary_mod_binomial.rs (line 8)

fn pow64(x: u64, mut y: u64, br: &BarrettReduction<u128>) -> u64 {
    let mut x = x as u128;
    let mut z: u128 = 1;
    while y > 0 {
        if y & 1 == 1 {
            z = br.rem(z * x);
        }
        x = br.rem(x * x);
        y >>= 1;
    }
    z as u64
}

fn pow32(x: u32, mut y: u64, br: &BarrettReduction<u64>) -> u32 {
    let mut x = x as u64;
    let mut z: u64 = 1;
    while y > 0 {
        if y & 1 == 1 {
            z = br.rem(z * x);
        }
        x = br.rem(x * x);
        y >>= 1;
    }
    z as u32
}

#[derive(Debug)]
struct PrimePowerBinomial {
    p: u64,
    e: u32,
    m: u64,
    size: usize,
    fact: Vec<u64>,
    inv_fact: Vec<u64>,
    delta: u64,
    bp: BarrettReduction<u64>,
    bm: BarrettReduction<u64>,
    bm128: BarrettReduction<u128>,
}

impl PrimePowerBinomial {
    fn new(p: u64, e: u32, max_n: u64) -> Self {
        let m = p.checked_pow(e).expect("prime power overflow");
        let bp = BarrettReduction::new(p);
        let bm = BarrettReduction::new(m);
        let bm128 = BarrettReduction::new(m as u128);
        let size = max_n.min(m - 1);
        assert!(size < usize::MAX as u64);
        let size = size as usize;
        let mut fact = vec![1u64; size + 1];
        let mut inv_fact = vec![1u64; size + 1];
        if m < 1 << 31 {
            for i in 2..=size {
                fact[i] = if bp.rem(i as u64) == 0 {
                    fact[i - 1]
                } else {
                    bm.rem(fact[i - 1] * i as u64)
                };
            }
            inv_fact[size] = fact[size].mod_inv(m);
            for i in (3..=size).rev() {
                inv_fact[i - 1] = if bp.rem(i as u64) == 0 {
                    inv_fact[i]
                } else {
                    bm.rem(inv_fact[i] * i as u64)
                };
            }
        } else {
            for i in 2..=size {
                fact[i] = if bp.rem(i as u64) == 0 {
                    fact[i - 1]
                } else {
                    bm128.rem(fact[i - 1] as u128 * i as u128) as u64
                };
            }
            inv_fact[size] = fact[size].mod_inv(m);
            for i in (3..=size).rev() {
                inv_fact[i - 1] = if bp.rem(i as u64) == 0 {
                    inv_fact[i]
                } else {
                    bm128.rem(inv_fact[i] as u128 * i as u128) as u64
                };
            }
        }
        let delta = if p == 2 && e >= 3 { 1 } else { m - 1 };
        Self {
            p,
            e,
            m,
            size,
            fact,
            inv_fact,
            delta,
            bp,
            bm,
            bm128,
        }
    }

    fn combination(&self, mut n: u64, mut k: u64) -> u64 {
        if k > n {
            return 0;
        }
        assert!(self.size as u64 >= n.min(self.m - 1));
        if self.m < 1 << 31 {
            let mut res = 1u64;
            if self.e == 1 {
                while n > 0 {
                    let (nn, n0) = self.bp.div_rem(n);
                    let (nk, k0) = self.bp.div_rem(k);
                    if n0 < k0 {
                        return 0;
                    }
                    res = self.bm.rem(res * self.fact[n0 as usize]);
                    res = self.bm.rem(res * self.inv_fact[k0 as usize]);
                    res = self.bm.rem(res * self.inv_fact[(n0 - k0) as usize]);
                    n = nn;
                    k = nk;
                }
            } else {
                let mut r = n - k;
                let mut e0 = 0;
                let mut eq = 0;
                let mut i = 0;
                while n > 0 {
                    res = self.bm.rem(res * self.fact[self.bm.rem(n) as usize]);
                    res = self.bm.rem(res * self.inv_fact[self.bm.rem(k) as usize]);
                    res = self.bm.rem(res * self.inv_fact[self.bm.rem(r) as usize]);
                    n = self.bp.div(n);
                    k = self.bp.div(k);
                    r = self.bp.div(r);
                    let eps = n - k - r;
                    e0 += eps;
                    if e0 >= self.e as u64 {
                        return 0;
                    }
                    i += 1;
                    if i >= self.e {
                        eq ^= eps & 1;
                    }
                }
                if eq == 1 {
                    res = self.bm.rem(res * self.delta);
                }
                res = self
                    .bm
                    .rem(res * pow32(self.p as _, e0 as _, &self.bm) as u64);
            }
            res
        } else {
            let mut res = 1u128;
            if self.e == 1 {
                while n > 0 {
                    let (nn, n0) = self.bp.div_rem(n);
                    let (nk, k0) = self.bp.div_rem(k);
                    if n0 < k0 {
                        return 0;
                    }
                    res = self.bm128.rem(res * self.fact[n0 as usize] as u128);
                    res = self.bm128.rem(res * self.inv_fact[k0 as usize] as u128);
                    res = self
                        .bm128
                        .rem(res * self.inv_fact[(n0 - k0) as usize] as u128);
                    n = nn;
                    k = nk;
                }
            } else {
                let mut r = n - k;
                let mut e0 = 0;
                let mut eq = 0;
                let mut i = 0;
                while n > 0 {
                    res = self
                        .bm128
                        .rem(res * self.fact[self.bm.rem(n) as usize] as u128);
                    res = self
                        .bm128
                        .rem(res * self.inv_fact[self.bm.rem(k) as usize] as u128);
                    res = self
                        .bm128
                        .rem(res * self.inv_fact[self.bm.rem(r) as usize] as u128);
                    n = self.bp.div(n);
                    k = self.bp.div(k);
                    r = self.bp.div(r);
                    let eps = n - k - r;
                    e0 += eps;
                    if e0 >= self.e as u64 {
                        return 0;
                    }
                    i += 1;
                    if i >= self.e {
                        eq ^= eps & 1;
                    }
                }
                if eq == 1 {
                    res = self.bm128.rem(res * self.delta as u128);
                }
                res = self.bm128.rem(res * pow64(self.p, e0, &self.bm128) as u128);
            }
            res as u64
        }
    }

crates/competitive/src/math/primitive_root.rs (line 27)

pub fn check_primitive_root(
    g: u64,
    phi: u64,
    br: &BarrettReduction<u128>,
    pf: &[(u64, u32)],
) -> bool {
    pf.iter().all(|&(q, _)| {
        let mut g = g as u128;
        let mut k = phi / q;
        let mut r: u128 = 1;
        while k > 0 {
            if k & 1 == 1 {
                r = br.rem(r * g);
            }
            g = br.rem(g * g);
            k >>= 1;
        }
        r > 1
    })
}

crates/competitive/src/math/prime_factors.rs (line 10)

fn find_factor(n: u64) -> Option<u64> {
    let br = BarrettReduction::<u128>::new(n as u128);
    if miller_rabin_with_br(n, &br) {
        return None;
    }
    let m = 1u64 << ((63 - n.leading_zeros()) / 5);
    let sub = |x: u64, y: u64| x.abs_diff(y);
    let mul = |x: u64, y: u64| br.rem(x as u128 * y as u128) as u64;
    for c in 12.. {
        let f = |x: u64| (br.rem(x as u128 * x as u128) + c) as u64;
        let (mut x, mut y, mut r, mut g, mut k, mut ys) = (0, 2, 1, 1, 0, 0);
        while g == 1 {
            x = y;
            for _ in 0..r {
                y = f(y);
            }
            while r > k && g == 1 {
                ys = y;
                let mut q = 1;
                for _ in 0..m.min(r - k) {
                    y = f(y);
                    q = mul(q, sub(x, y));
                }
                g = gcd(q, n);
                k += m;
            }
            r <<= 1;
        }
        if g == n {
            g = 1;
            while g == 1 {
                ys = f(ys);
                g = gcd(sub(x, ys), n);
            }
        }
        if g < n {
            return Some(g);
        }
    }
    unreachable!();
}

Trait Implementations

impl<T: Clone> Clone for BarrettReduction<T>

fn clone(&self) -> BarrettReduction<T>

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<T: Debug> Debug for BarrettReduction<T>

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<T: Copy> Copy for BarrettReduction<T>