#![allow(clippy::module_name_repetitions)]
mod matrix;
pub use self::matrix::Matrix as LehmerMatrix;
use crate::Uint;
use core::mem::swap;
#[must_use]
pub fn gcd<const BITS: usize, const LIMBS: usize>(
mut a: Uint<BITS, LIMBS>,
mut b: Uint<BITS, LIMBS>,
) -> Uint<BITS, LIMBS> {
if b > a {
swap(&mut a, &mut b);
}
while b != Uint::ZERO {
debug_assert!(a >= b);
let m = LehmerMatrix::from(a, b);
if m == LehmerMatrix::IDENTITY {
a %= b;
swap(&mut a, &mut b);
} else {
m.apply(&mut a, &mut b);
}
}
a
}
#[must_use]
pub fn gcd_extended<const BITS: usize, const LIMBS: usize>(
mut a: Uint<BITS, LIMBS>,
mut b: Uint<BITS, LIMBS>,
) -> (
Uint<BITS, LIMBS>,
Uint<BITS, LIMBS>,
Uint<BITS, LIMBS>,
bool,
) {
if BITS == 0 {
return (Uint::ZERO, Uint::ZERO, Uint::ZERO, false);
}
let swapped = a < b;
if swapped {
swap(&mut a, &mut b);
}
let mut s0 = Uint::from(1);
let mut s1 = Uint::ZERO;
let mut t0 = Uint::ZERO;
let mut t1 = Uint::from(1);
let mut even = true;
while b != Uint::ZERO {
debug_assert!(a >= b);
let m = LehmerMatrix::from(a, b);
if m == LehmerMatrix::IDENTITY {
let q = a / b;
a -= q * b;
swap(&mut a, &mut b);
s0 -= q * s1;
swap(&mut s0, &mut s1);
t0 -= q * t1;
swap(&mut t0, &mut t1);
even = !even;
} else {
m.apply(&mut a, &mut b);
m.apply(&mut s0, &mut s1);
m.apply(&mut t0, &mut t1);
even ^= !m.4;
}
}
if even {
t0 = Uint::ZERO - t0;
} else {
s0 = Uint::ZERO - s0;
}
if swapped {
swap(&mut s0, &mut t0);
even = !even;
}
(a, s0, t0, even)
}
#[must_use]
pub fn inv_mod<const BITS: usize, const LIMBS: usize>(
num: Uint<BITS, LIMBS>,
modulus: Uint<BITS, LIMBS>,
) -> Option<Uint<BITS, LIMBS>> {
if BITS == 0 || modulus == Uint::ZERO {
return None;
}
let mut a = modulus;
let mut b = num;
if b >= a {
b %= a;
}
if b == Uint::ZERO {
return None;
}
let mut t0 = Uint::ZERO;
let mut t1 = Uint::from(1);
let mut even = true;
while b != Uint::ZERO {
debug_assert!(a >= b);
let m = LehmerMatrix::from(a, b);
if m == LehmerMatrix::IDENTITY {
let q = a / b;
a -= q * b;
swap(&mut a, &mut b);
t0 -= q * t1;
swap(&mut t0, &mut t1);
even = !even;
} else {
m.apply(&mut a, &mut b);
m.apply(&mut t0, &mut t1);
even ^= !m.4;
}
}
if a == Uint::from(1) {
Some(if even { modulus + t0 } else { t0 })
} else {
None
}
}
#[cfg(test)]
#[allow(clippy::cast_lossless)]
mod tests {
use super::*;
use crate::{const_for, nlimbs};
use core::cmp::min;
use proptest::{proptest, test_runner::Config};
#[test]
fn test_gcd_one() {
use std::str::FromStr;
const BITS: usize = 129;
const LIMBS: usize = nlimbs(BITS);
type U = Uint<BITS, LIMBS>;
let a = U::from_str("0x006d7c4641f88b729a97889164dd8d07db").unwrap();
let b = U::from_str("0x01de6ef6f3caa963a548d7a411b05b9988").unwrap();
assert_eq!(gcd(a, b), gcd_ref(a, b));
}
fn gcd_ref<const BITS: usize, const LIMBS: usize>(
mut a: Uint<BITS, LIMBS>,
mut b: Uint<BITS, LIMBS>,
) -> Uint<BITS, LIMBS> {
while b != Uint::ZERO {
a %= b;
swap(&mut a, &mut b);
}
a
}
#[test]
#[allow(clippy::absurd_extreme_comparisons)] fn test_gcd() {
const_for!(BITS in SIZES {
const LIMBS: usize = nlimbs(BITS);
type U = Uint<BITS, LIMBS>;
let mut config = Config::default();
config.cases = min(config.cases, if BITS > 500 { 9 } else { 30 });
proptest!(config, |(a: U, b: U)| {
assert_eq!(gcd(a, b), gcd_ref(a, b));
});
});
}
#[test]
#[allow(clippy::absurd_extreme_comparisons)] fn test_gcd_extended() {
const_for!(BITS in SIZES {
const LIMBS: usize = nlimbs(BITS);
type U = Uint<BITS, LIMBS>;
let mut config = Config::default();
config.cases = min(config.cases, if BITS > 500 { 3 } else { 10 });
proptest!(config, |(a: U, b: U)| {
let (g, x, y, sign) = gcd_extended(a, b);
assert_eq!(g, gcd_ref(a, b));
if sign {
assert_eq!(a * x - b * y, g);
} else {
assert_eq!(b * y - a * x, g);
}
});
});
}
}
#[cfg(feature = "bench")]
#[doc(hidden)]
pub mod bench {
use super::*;
use criterion::Criterion;
pub fn group(criterion: &mut Criterion) {
matrix::bench::group(criterion);
}
}