This commit is contained in:
goldenMetteyya 2019-04-27 03:04:52 +03:00
parent a977d93720
commit d4712d31cc
3 changed files with 29 additions and 479 deletions

View File

@ -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<BigUint>,
b_in: Cow<BigUint>,
extended: bool,
) -> (BigInt, Option<BigInt>, Option<BigInt>) {
if a_in.is_zero() && b_in.is_zero() {
if extended {
return (b_in.to_bigint().unwrap(), Some(0.into()), Some(0.into()));
@ -43,7 +40,6 @@ pub fn extended_gcd(
}
}
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<BigUint>, b: Cow<BigUint>) -> (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() {

View File

@ -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::*;

View File

@ -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<BigInt>,
r1_in: Cow<BigInt>,
bound: Cow<BigInt>,
) -> (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);
}
}