Skip to main content

solve_linear_congruence

Function solve_linear_congruence 

Source
fn solve_linear_congruence(a: u64, b: u64, m: u64) -> Option<(u64, u64)>
Examples found in repository?
crates/competitive/src/math/discrete_logarithm.rs (line 335)
315    fn discrete_logarithm(
316        &self,
317        a: u64,
318        b: u64,
319        br_primes: &[BarrettReduction<u64>],
320    ) -> Option<(u64, u64)> {
321        let p = self.p;
322        let ord = self.ord;
323        let br = BarrettReduction::<u128>::new(p as u128);
324        let a = br.rem(a as _) as u64;
325        let b = br.rem(b as _) as u64;
326        if a == 0 {
327            return if b == 0 { Some((1, 1)) } else { None };
328        }
329        if b == 0 {
330            return None;
331        }
332
333        let x = self.index_calculus(a, br_primes)?;
334        let y = self.index_calculus(b, br_primes)?;
335        solve_linear_congruence(x, y, ord)
336    }
337}
338
339thread_local!(
340    static IC: UnsafeCell<IndexCalculus> = UnsafeCell::new(IndexCalculus::new());
341);
342
343pub fn discrete_logarithm_prime_mod(a: u64, b: u64, p: u64) -> Option<u64> {
344    IC.with(|ic| unsafe { &mut *ic.get() }.discrete_logarithm(a, b, p))
345        .map(|t| t.0)
346}
347
348/// a^x ≡ b (mod n), a has order p^e
349fn pohlig_hellman_prime_power_order(a: u64, b: u64, n: u64, p: u64, e: u32) -> Option<u64> {
350    let br = BarrettReduction::<u128>::new(n as u128);
351    let mul = |x: u64, y: u64| br.rem(x as u128 * y as u128) as u64;
352    let block_size = (p as f64).sqrt().ceil() as u64;
353    let mut baby = HashMap::<u64, u64>::new();
354    let g = pow(a, p.pow(e - 1), &br);
355    let mut xj = 1;
356    for j in 0..block_size {
357        baby.entry(xj).or_insert(j);
358        xj = mul(xj, g);
359    }
360    let xi = modinv(xj, n);
361    let mut t = 0u64;
362    for k in 0..e {
363        let mut h = pow(mul(modinv(pow(a, t, &br), n), b), p.pow(e - 1 - k), &br);
364        let mut ok = false;
365        for i in (0..block_size * block_size).step_by(block_size as usize) {
366            if let Some(j) = baby.get(&h) {
367                t += (i + j) * p.pow(k);
368                ok = true;
369                break;
370            }
371            h = mul(h, xi);
372        }
373        if !ok {
374            return None;
375        }
376    }
377    Some(t)
378}
379
380/// a^x ≡ b (mod p^e)
381fn discrete_logarithm_prime_power(a: u64, b: u64, p: u64, e: u32) -> Option<(u64, u64)> {
382    assert_ne!(p, 0);
383    assert_ne!(e, 0);
384    let n = p.pow(e);
385    assert!(a < n);
386    assert!(b < n);
387    assert_eq!(gcd(a, p), 1);
388    if p == 1 {
389        return Some((0, 1));
390    }
391    if a == 0 {
392        return if b == 0 { Some((1, 1)) } else { None };
393    }
394    if b == 0 {
395        return None;
396    }
397    if e == 1 {
398        return IC.with(|ic| unsafe { &mut *ic.get() }.discrete_logarithm(a, b, p));
399    }
400    let br = BarrettReduction::<u128>::new(n as _);
401    if p == 2 {
402        if e >= 3 {
403            if a % 4 == 1 && b % 4 != 1 {
404                return None;
405            }
406            let aa = if a % 4 == 1 { a } else { n - a };
407            let bb = if b % 4 == 1 { b } else { n - b };
408            let g = 5;
409            let ord = n / 4;
410            let x = pohlig_hellman_prime_power_order(g, aa, n, p, e - 2)?;
411            let y = pohlig_hellman_prime_power_order(g, bb, n, p, e - 2)?;
412            let t = solve_linear_congruence(x, y, ord)?;
413            match (a % 4 == 1, b % 4 == 1) {
414                (true, true) => Some(t),
415                (false, true) if t.0 % 2 == 0 => Some((t.0, lcm(t.1, 2))),
416                (false, false) if t.0 % 2 == 1 => Some((t.0, lcm(t.1, 2))),
417                (false, false) if a == b => Some((1, lcm(t.1, 2))),
418                _ => None,
419            }
420        } else if a == 1 {
421            if b == 1 { Some((0, 1)) } else { None }
422        } else {
423            assert_eq!(a, 3);
424            if b == 1 {
425                Some((0, 2))
426            } else if b == 3 {
427                Some((1, 2))
428            } else {
429                None
430            }
431        }
432    } else {
433        let ord = n - n / p;
434        let pf_ord = prime_factors(ord);
435        let g = (2..)
436            .find(|&g| check_primitive_root(g, ord, &br, &pf_ord))
437            .unwrap();
438        let mut pf_p = prime_factors(p - 1);
439        pf_p.push((p, e - 1));
440        let mut abm = vec![];
441        for (q, c) in pf_p {
442            let m = q.pow(c);
443            let d = ord / m;
444            let gg = pow(g, d, &br);
445            let aa = pow(a, d, &br);
446            let bb = pow(b, d, &br);
447            let x = pohlig_hellman_prime_power_order(gg, aa, n, q, c)?;
448            let y = pohlig_hellman_prime_power_order(gg, bb, n, q, c)?;
449            abm.push((x, y, m));
450        }
451        solve_linear_congruences(abm)
452    }
453}