diff --git a/src/algorithms/gcd.rs b/src/algorithms/gcd.rs index 4a03bb4..bc3834e 100644 --- a/src/algorithms/gcd.rs +++ b/src/algorithms/gcd.rs @@ -11,14 +11,11 @@ use crate::biguint::{BigUint, IntDigits}; /// Uses the lehemer algorithm. /// Based on https://github.com/golang/go/blob/master/src/math/big/int.go#L612 /// If `extended` is set, the Bezout coefficients are calculated, otherwise they are `None`. -/// If x or y are not nil, GCD sets their value such that z = a*x + b*y. -/// If either a or b is <= 0, GCD sets z = x = y = 0. pub fn extended_gcd( a_in: Cow, b_in: Cow, extended: bool, ) -> (BigInt, Option, Option) { - if a_in.is_zero() && b_in.is_zero() { if extended { return (b_in.to_bigint().unwrap(), Some(0.into()), Some(0.into())); @@ -42,7 +39,6 @@ pub fn extended_gcd( return (a_in.to_bigint().unwrap(), None, None); } } - let a_in = a_in.to_bigint().unwrap(); let b_in = b_in.to_bigint().unwrap(); @@ -60,17 +56,14 @@ pub fn extended_gcd( std::mem::swap(&mut ua, &mut ub); } - // temp variables for multiprecision update let mut q: BigInt = 0.into(); let mut r: BigInt = 0.into(); let mut s: BigInt = 0.into(); let mut t: BigInt = 0.into(); - // loop invariant A >= B while b.len() > 1 { // Attempt to calculate in single-precision using leading words of a and b. let (u0, u1, v0, v1, even) = lehmer_simulate(&a, &b); - // multiprecision step if v0 != 0 { // Simulate the effect of the single-precision steps using cosequences. @@ -144,19 +137,14 @@ pub fn extended_gcd( t.data.set_digit(ua_word); s.data.set_digit(va); - t.sign = if even { Plus } else { Minus }; s.sign = if even { Minus } else { Plus }; - if let Some(ua) = ua.as_mut() { - t *= &*ua; - s *= ub.unwrap(); *ua = &t + &s; - } } else { while b_word != 0 { @@ -177,12 +165,10 @@ pub fn extended_gcd( } else { None }; - //println!("exteneded gcd: {:?}", (a, ua, y)); + (a, ua, y) } - - /// Attempts to simulate several Euclidean update steps using leading digits of `a` and `b`. /// It returns `u0`, `u1`, `v0`, `v1` such that `a` and `b` can be updated as: /// a = u0 * a + v0 * b @@ -219,14 +205,9 @@ fn lehmer_simulate(a: &BigInt, b: &BigInt) -> (BigDigit, BigDigit, BigDigit, Big 0 }; - // Since we are calculating with full words to avoid overflow, - // we use 'even' to track the sign of the cosequences. - // For even iterations: u0, v1 >= 0 && u1, v0 <= 0 - // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0 - // The first iteration starts with k=1 (odd). + // odd, even tracking let mut even = false; - // variables to track the cosequences let mut u0 = 0; let mut u1 = 1; let mut u2 = 0; @@ -235,10 +216,7 @@ fn lehmer_simulate(a: &BigInt, b: &BigInt) -> (BigDigit, BigDigit, BigDigit, Big let mut v1 = 0; let mut v2 = 1; - // Calculate the quotient and cosequences using Collins' stopping condition. - // Note that overflow of a Word is not possible when computing the remainder - // sequence and cosequences since the cosequence size is bounded by the input size. - // See section 4.2 of Jebelean for details. + // Calculate the quotient and cosequences using Collins' stoppting condition. while a2 >= v2 && a1.wrapping_sub(a2) >= v1 + v2 { let q = a1 / a2; let r = a1 % a2; @@ -262,13 +240,6 @@ fn lehmer_simulate(a: &BigInt, b: &BigInt) -> (BigDigit, BigDigit, BigDigit, Big (u0, u1, v0, v1, even) } -// lehmerUpdate updates the inputs A and B such that: -// A = u0*A + v0*B -// B = u1*A + v1*B -// where the signs of u0, u1, v0, v1 are given by even -// For even == true: u0, v1 >= 0 && u1, v0 <= 0 -// For even == false: u0, v1 <= 0 && u1, v0 >= 0 -// q, r, s, t are temporary variables to avoid allocations in the multiplication fn lehmer_update( a: &mut BigInt, b: &mut BigInt, @@ -282,26 +253,13 @@ fn lehmer_update( v1: BigDigit, even: bool, ) { - t.data.set_digit(u0); s.data.set_digit(v0); - - //We handle to edge case that s or t is a zero and make sure that we dont have negative zero. if even { t.sign = Plus; - if s.data.is_zero() { - s.sign = Plus; - } else { - s.sign = Minus; - } - + s.sign = Minus } else { - // t.sign = Minus; - if t.data.is_zero() { - t.sign = Plus; - } else { - t.sign = Minus; - } + t.sign = Minus; s.sign = Plus; } @@ -325,8 +283,6 @@ fn lehmer_update( *b = r + q; } -// euclidUpdate performs a single step of the Euclidean GCD algorithm -// if extended is true, it also updates the cosequence Ua, Ub fn euclid_udpate( a: &mut BigInt, b: &mut BigInt, @@ -369,10 +325,14 @@ mod tests { #[cfg(feature = "rand")] use num_traits::{One, Zero}; #[cfg(feature = "rand")] - use rand::{SeedableRng, XorShiftRng}; + use rand::SeedableRng; + #[cfg(feature = "rand")] + use rand_xorshift::XorShiftRng; #[cfg(feature = "rand")] fn extended_gcd_euclid(a: Cow, b: Cow) -> (BigInt, BigInt, BigInt) { + use crate::bigint::ToBigInt; + if a.is_zero() && b.is_zero() { return (0.into(), 0.into(), 0.into()); } @@ -438,6 +398,25 @@ mod tests { assert_eq!(t_k, None); } + #[test] + fn test_extended_gcd_example_wolfram() { + // https://www.wolframalpha.com/input/?i=ExtendedGCD%5B-565721958+,+4486780496%5D + // https://github.com/Chia-Network/oldvdf-competition/blob/master/tests/test_classgroup.py#L109 + + let a = BigInt::from(-565721958); + let b = BigInt::from(4486780496u64); + + let (q, _s_k, _t_k) = extended_gcd( + Cow::Owned(a.to_biguint().unwrap()), + Cow::Owned(b.to_biguint().unwrap()), + true, + ); + + assert_eq!(q, BigInt::from(2)); + assert_eq!(_s_k, Some(BigInt::from(-1090996795))); + assert_eq!(_t_k, Some(BigInt::from(-137559848))); + } + #[test] #[cfg(feature = "rand")] fn test_gcd_lehmer_euclid_extended() { diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs index dec5f2a..38014ec 100644 --- a/src/algorithms/mod.rs +++ b/src/algorithms/mod.rs @@ -9,7 +9,6 @@ mod jacobi; mod mac; mod mod_inverse; mod mul; -mod partial; mod shl; mod shr; mod sub; @@ -23,7 +22,6 @@ pub use self::jacobi::*; pub use self::mac::*; pub use self::mod_inverse::*; pub use self::mul::*; -pub use self::partial::*; pub use self::shl::*; pub use self::shr::*; pub use self::sub::*; diff --git a/src/algorithms/partial.rs b/src/algorithms/partial.rs deleted file mode 100644 index a9fd188..0000000 --- a/src/algorithms/partial.rs +++ /dev/null @@ -1,427 +0,0 @@ -// Partial Euclidean algorithm. -// Lehmer's version for computing GCD -// Ported from Flint-2.4.1, fmpz_xgcd_partial - -use crate::big_digit::BITS as LIMB_BITS; -use crate::bigint::{BigInt, Sign}; -use bigint::ToBigInt; -use integer::Integer; -use num_traits::{One, Pow, ToPrimitive, Zero}; -use std::borrow::Cow; -use std::ops::Neg; -use biguint::IntDigits; - - - -#[inline] -fn signed_shift(op: u64, shift: i32) -> u64 { - //let ushift = shift as u32; - - if shift > 0 { - return op << shift; - } - - if shift <= -64 { - return 0; - } - - return op >> (-shift); -} - -// //https://github.com/thomaso-mirodin/intmath/blob/master/u64/log2.go#L5 -// // Log2 returns log base 2 of n. It's the same as index of the highest -// // bit set in n. n == 0 returns 0 -// fn log2_64(n: u64) -> u64 { -// // Using uint instead of uint64 is about 25% faster -// // on x86 systems with the default Go compiler. -// let mut r: usize = 0; -// let mut _v: usize = 0; - -// if n < 1<<32 { -// _v = n as usize; -// } else { -// r = 32; -// _v = (n >> 32) as usize; -// } -// if _v >= 1<<16 { -// r += 16; -// _v >>= 16; -// } -// if _v >= 1<<8 { -// r += 8; -// _v >>= 8; -// } -// if _v >= 1<<4 { -// r += 4; -// _v >>= 4; -// } -// if _v >= 1<<2 { -// r += 2; -// _v >>= 2; -// } -// r += _v >> 1; -// return r as u64; -// } - - - -fn log2_64(mut value: u64) -> u32 { - let table = [ - 63, 0, 58, 1, 59, 47, 53, 2, - 60, 39, 48, 27, 54, 33, 42, 3, - 61, 51, 37, 40, 49, 18, 28, 20, - 55, 30, 34, 11, 43, 14, 22, 4, - 62, 57, 46, 52, 38, 26, 32, 41, - 50, 36, 17, 19, 29, 10, 13, 21, - 56, 45, 25, 31, 35, 16, 9, 12, - 44, 24, 15, 8, 23, 7, 6, 5 - ]; - - value |= value >> 1; - value |= value >> 2; - value |= value >> 4; - value |= value >> 8; - value |= value >> 16; - value |= value >> 32; - println!("-------lg2 value------: {:?}", value); - - return table[((value * 0x03f6eaf2cd271461) >> 58) as usize]; - //return tab64[((((value - (value >> 1))*0x07EDD5E59A4E28C2)) >> 58) as usize]; -} - -// const int tab32[32] = { -// 0, 9, 1, 10, 13, 21, 2, 29, -// 11, 14, 16, 18, 22, 25, 3, 30, -// 8, 12, 20, 28, 15, 17, 24, 7, -// 19, 27, 23, 6, 26, 5, 4, 31}; - -// int log2_32 (uint32_t value) -// { -// value |= value >> 1; -// value |= value >> 2; -// value |= value >> 4; -// value |= value >> 8; -// value |= value >> 16; -// return tab32[(uint32_t)(value*0x07C4ACDD) >> 27]; -// } - - -//Port to rust from c++ -//https://github.com/Chia-Network/vdftrack1results/blob/b171a420ecfffb2d956979da7b1b183c4c88d2a0/akashnil/entry/vdf4.cpp#L82 -// Return an approximation x of the large mpz_t op by an int64_t and the exponent e adjustment. -// We must have (x * 2^e) / op = constant approximately. -#[inline] -pub fn partial_bigint(op: &BigInt) -> (i64, i32) { - //uint64_t size = mpz_size(op); - //number if limbs used to represent this number - let digies = op.digits(); - let size = digies.len(); - // - let last: u64 = digies[size-1]; - - //uint64_t last = mpz_getlimbn(op, size - 1); - //let last: f64 = (size % LIMB_BITS) as f64; - //uint64_t ret; - let mut _ret = 0u64; - - //let lg2: f64 = (last as f64).log2() + 1.0; - // // extract the top word of bits from a and b - // let h = a.digits()[n - 1].leading_zeros(); - //let lg2 = last.leading_zeros(); - let lg2 = log2_64(last) + 1; - println!("-------lg2 ------: {:?}", lg2); - let mut exp = lg2 as i32; - - _ret = signed_shift(last, 63 - exp); - - if size > 1 { - let side = size - 1; - println!("-------lside ------: {:?}", side); - exp += (side * 64) as i32; - // uint64_t prev = mpz_getlimbn(op, size - 2); - let prev: u64 = digies[size-2]; - _ret += signed_shift(prev, -1i32 - (lg2 as i32)); - } - - //if (mpz_sgn(op) < 0) return - ((int64_t)ret); - if op.sign() == Sign::Minus { - return (-(_ret as i64), exp); - } - - return ((_ret as i64), exp); -} - - - - -/// This function is an implementation of Lehmer extended GCD with early termination. -/// It terminates early when remainders fall below the specified bound. -/// The initial values r1 and r2 are treated as successive remainders in the Euclidean algorithm -/// and are replaced with the last two remainders computed. The values _co1 and _co2 are the last two -/// cofactors and satisfy the identity _co2*r1 - _co1*r2 == +/- r2_orig upon termination, where -/// r2_orig is the starting value of r2 supplied, and r1 and r2 are the final values. -pub fn partial_extended_gcd( - r2_in: Cow, - r1_in: Cow, - bound: Cow, -) -> (BigInt, BigInt, BigInt, BigInt) { - //return (R2, R1, C1, C2) - - let mut _r: BigInt = Zero::zero(); - - let mut _co1 = BigInt::one(); - _co1.sign = _co1.sign.neg(); - let mut _co2 = BigInt::zero(); - - let r1_in = r1_in.to_bigint().unwrap(); - let r2_in = r2_in.to_bigint().unwrap(); - - let mut r1 = r1_in.clone(); - let mut r2 = r2_in.clone(); - - //loop index - let mut _index = 0; - - while !r1.is_zero() && r1 > *bound { - //get bits length - let r2_bits = r2.bits(); - let r1_bits = r1.bits(); - let one_under_limb = LIMB_BITS - 1; - - let mut _t = (r2_bits as isize) - (one_under_limb as isize); - let mut _t1 = (r1_bits as isize) - (one_under_limb as isize); - - //Bits - if _t < _t1 { - _t = _t1 - } - if _t < 0 { - _t = 0 - } - - let mut t: isize = _t; - - //truncate for a positive number is same as floor - //truncate for a negative number is same as ceil - - //r = R2 / (2 ^ T); truncate r - let (d, m) = r2.div_mod_floor(&BigInt::from(2).pow(t as usize)); - - //positive sign or no sign - if d.sign() == Sign::Minus { - //negative numbers we all to ceil - - if m.is_zero() { - _r = d; - } else { - _r = d + BigInt::one(); - } - } else { - _r = d; - } - - let mut rr2: isize = _r.to_isize().unwrap(); - - //r = R1 / (2 ^ T); truncate r - let (d, m) = r1.div_mod_floor(&BigInt::from(2).pow(t as usize)); - - //positive sign or no sign - if d.sign() == Sign::Minus { - //negative numbers we all to ceil - println!("suckerssss 2"); - if m.is_zero() { - _r = d; - } else { - _r = d + BigInt::one(); - } - } else { - _r = d; - } - - //positive sign or no sign - let mut rr1: isize = _r.to_isize().unwrap(); - - //r = R1 / (2 ^ T); truncate r - _r = bound.div_floor(&BigInt::from(2).pow(t as usize)); - - //r = bound / (2 ^ T); truncate r - let (d, m) = bound.div_mod_floor(&BigInt::from(2).pow(t as usize)); - - //positive sign or no sign - if d.sign() == Sign::Minus { - //negative numbers we all to ceil - - if m.is_zero() { - _r = d; - } else { - _r = d + BigInt::one(); - } - } else { - _r = d; - } - - let bb: isize = _r.to_isize().unwrap(); //might need tobe isize - - //reset values - let mut a1: isize = 1; - let mut a2: isize = 0; - let mut b1: isize = 0; - let mut b2: isize = 1; - - //reset inner loop index - _index = 0; - - // Euclidean Step - while rr1 != 0 && rr1 > bb { - let qq: isize = rr2 / rr1; - - //t1 - t = rr2 - qq * rr1; - rr2 = rr1; - rr1 = t; - - //t2 - t = a2 - qq * a1; - a2 = a1; - a1 = t; - - //t3 - t = b2 - qq * b1; - b2 = b1; - b1 = t; - - //check if it is even or odd - if _index % 2 != 0 { - //index & 1 - //its odd - if rr1 < -b1 || rr2 - rr1 < a1 - a2 { - break; - } - } else { - //its even - if rr1 < -a1 || rr2 - rr1 < b1 - b2 { - break; - } - } - - //increment counter - _index += 1; - } - - if _index == 0 { - // multiprecision step - let q: BigInt = r2.div_floor(&r1); - r2 = &r2 % &r1; - std::mem::swap(&mut r2, &mut r1); - _co2 = &_co2 - (&q * &_co1); - std::mem::swap(&mut _co2, &mut _co1); - } else { - // recombination - _r = &r2 * &b2; - - if a2 >= 0 { - _r = &_r + &r1 * &a2; - } else { - _r = &_r - (&r1 * &-a2); - } - - r1 = &r1 * &a1; - if b1 >= 0 { - r1 = &r1 + &r2 * &b1; - } else { - r1 = &r1 - (&r2 * &-b1); - } - - r2 = _r.clone(); - - _r = &_co2 * &b2; - - if a2 >= 0 { - _r = &_r + &_co1 * &a2; - } else { - _r = &_r - (&_co1 * &-a2); - } - - _co1 = &_co1 * &a1; - if b1 >= 0 { - _co1 = &_co1 + &_co2 * &b1; - } else { - _co1 = &_co1 - (&_co2 * &-b1); - } - - // C2 = r; - _co2 = _r.clone(); - - // make sure R1 is positive - if r1.sign() == Sign::Minus { - _co1.sign = _co1.sign.neg(); - r1.sign = r1.sign.neg(); - } - // make sure R2 is positive - if r2.sign() == Sign::Minus { - _co2.sign = _co2.sign.neg(); - r2.sign = r2.sign.neg(); - } - } - } - - // make sure R2 is positive - if r2.sign() == Sign::Minus { - _co1.sign = _co1.sign.neg(); - _co2.sign = _co2.sign.neg(); - r2.sign = r2.sign.neg(); - } - - //return back - (_co2, _co1, r2, r1) -} - -#[cfg(test)] -mod test { - - #[cfg(feature = "rand")] - use crate::bigrand::RandBigInt; - #[cfg(feature = "rand")] - use num_traits::Signed; - #[cfg(feature = "rand")] - use rand::{SeedableRng, XorShiftRng}; - - #[test] - #[cfg(feature = "rand")] - fn test_partial_extended_gcd() { - use super::*; - - let mut rng = XorShiftRng::from_seed([1u8; 16]); - - /* Test _co2*r1 - _co1*r2 = r2_orig */ - let mut _co1 = BigInt::zero(); - let mut _co2 = BigInt::zero(); - let mut _f = BigInt::zero(); - let mut g = rng.gen_bigint(2000); - //let mut g = BigInt::zero(); - let mut _t1 = BigInt::zero(); - let mut _t2 = BigInt::zero(); - let mut _l = BigInt::zero(); - - g += BigInt::one(); - _f = BigInt::from_biguint(Sign::Plus, rng.gen_biguint_below(&g.to_biguint().unwrap())); - _l = rng.gen_bigint(1000); - //println!("L: {:?}", L); - - _t2 = g.clone(); - _t2 = _t2.abs(); - - let (_co2, _co1, r2, r1) = - partial_extended_gcd(Cow::Borrowed(&g), Cow::Borrowed(&_f), Cow::Borrowed(&_l)); - - _t1 = &_co2 * &r1; - _t1 -= &_co1 * &r2; - _t1 = _t1.abs(); - // println!("------------------test_partial_extended_gcd"); - // println!("t1: {:?}", t1); - // println!("t2: {:?}", t2); - - assert_eq!(&_t1, &_t2); - } - -}