Skip to content

Commit

Permalink
Improve division performance
Browse files Browse the repository at this point in the history 
Before:
test divide_0      ... bench:       4,058 ns/iter (+/- 255)
test divide_1      ... bench:     304,507 ns/iter (+/- 28,063)
test divide_2      ... bench: 668,293,969 ns/iter (+/- 25,383,239)

After:
test divide_0      ... bench:         874 ns/iter (+/- 71)
test divide_1      ... bench:      16,641 ns/iter (+/- 1,205)
test divide_2      ... bench:   1,336,888 ns/iter (+/- 77,450)
  • Loading branch information
@koverstreet
koverstreet committed Nov 26, 2015
1 parent 7c1bf88 commit 7976257
Showing 1 changed file with 118 additions and 107 deletions.
225 changes: 118 additions & 107 deletions src/bigint.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,12 @@ fn mac_with_carry(a: BigDigit, b: BigDigit, c: BigDigit, carry: &mut BigDigit) -
lo
}

#[inline]
fn div_wide(hi: BigDigit, lo: BigDigit, divisor: BigDigit) -> (BigDigit, BigDigit) {
let lhs = big_digit::to_doublebigdigit(hi, lo);
let rhs = divisor as DoubleBigDigit;
((lhs / rhs) as BigDigit, (lhs % rhs) as BigDigit)
}
/// A big unsigned integer type.
///
/// A `BigUint`-typed value `BigUint { data: vec!(a, b, c) }` represents a number
Expand All @@ -182,18 +188,25 @@ impl PartialOrd for BigUint {
}
}

fn cmp_slice(a: &[BigDigit], b: &[BigDigit]) -> Ordering {
debug_assert!(a.last() != Some(&0));
debug_assert!(b.last() != Some(&0));

let (s_len, o_len) = (a.len(), b.len());
if s_len < o_len { return Less; }
if s_len > o_len { return Greater; }

for (&self_i, &other_i) in a.iter().rev().zip(b.iter().rev()) {
if self_i < other_i { return Less; }
if self_i > other_i { return Greater; }
}
return Equal;
}

impl Ord for BigUint {
#[inline]
fn cmp(&self, other: &BigUint) -> Ordering {
let (s_len, o_len) = (self.data.len(), other.data.len());
if s_len < o_len { return Less; }
if s_len > o_len { return Greater; }

for (&self_i, &other_i) in self.data.iter().rev().zip(other.data.iter().rev()) {
if self_i < other_i { return Less; }
if self_i > other_i { return Greater; }
}
return Equal;
cmp_slice(&self.data[..], &other.data[..])
}
}

Expand Down Expand Up @@ -523,7 +536,7 @@ fn add2(a: &mut [BigDigit], b: &[BigDigit]) {
}
}

debug_assert!(b_iter.next() == None && carry == 0);
assert!(b_iter.next() == None && carry == 0);
}

impl<'a> Add<&'a BigUint> for BigUint {
Expand Down Expand Up @@ -787,6 +800,18 @@ impl<'a, 'b> Mul<&'b BigUint> for &'a BigUint {
}
}

fn div_rem_digit(mut a: BigUint, b: BigDigit) -> (BigUint, BigDigit) {
let mut rem = 0;

for d in a.data.iter_mut().rev() {
let (q, r) = div_wide(rem, *d, b);
*d = q;
rem = r;
}

(a.normalize(), rem)
}

forward_all_binop_to_ref_ref!(impl Div for BigUint, div);

impl<'a, 'b> Div<&'b BigUint> for &'a BigUint {
Expand Down Expand Up @@ -888,77 +913,88 @@ impl Integer for BigUint {
Greater => {} // Do nothing
}

let mut shift = 0;
let mut n = *other.data.last().unwrap();
while n < (1 << big_digit::BITS - 2) {
n <<= 1;
shift += 1;
}
assert!(shift < big_digit::BITS);
let (d, m) = div_mod_floor_inner(self << shift, other << shift);
return (d, m >> shift);

/*
* This algorithm is from Knuth, TAOCP vol 2 section 4.3, algorithm D:
*
* First, normalize the arguments so the highest bit in the highest digit of the divisor is
* set: the main loop uses the highest digit of the divisor for generating guesses, so we
* want it to be the largest number we can efficiently divide by.
*/

fn div_mod_floor_inner(a: BigUint, b: BigUint) -> (BigUint, BigUint) {
let mut m = a;
let mut d: BigUint = Zero::zero();
let mut n = 1;
while m >= b {
let (d0, d_unit, b_unit) = div_estimate(&m, &b, n);
let mut d0 = d0;
let mut prod = &b * &d0;
while prod > m {
// FIXME(#5992): assignment operator overloads
// d0 -= &d_unit
d0 = d0 - &d_unit;
// FIXME(#5992): assignment operator overloads
// prod -= &b_unit;
prod = prod - &b_unit
}
if d0.is_zero() {
n = 2;
let shift = other.data.last().unwrap().leading_zeros() as usize;
let (q, m) = div_mod_floor_inner(self << shift, other << shift);
return (q, m >> shift);


fn div_mod_floor_inner(mut a: BigUint, b: BigUint) -> (BigUint, BigUint) {
/*
* The algorithm works by incrementally calculating "guesses", q0, for part of the
* remainder. Once we have any number q0 such that q0 * b <= a, we can set
*
* q += q0
* a -= q0 * b
*
* and then iterate until a < b. Then, (q, a) will be our desired quotient and
* remainder.
*
* q0, our guess, is calculated by dividing the last few digits of a by the last digit
* of b - this should give us a guess that is "close" to the actual quotient, but is
* possibly greater than the actual quotient. If q0 * b > a, we simply use iterated
* subtraction until we have a guess such that q8 & b <= a.
*/

let bn = *b.data.last().unwrap();
let q_len = a.data.len() - b.data.len() + 1;
let mut q: BigUint = BigUint { data: Vec::with_capacity(q_len) };

q.data.extend(repeat(0).take(q_len));

/* We reuse the same temporary to avoid hitting the allocator in our inner loop: */
let mut tmp: BigUint = BigUint { data: Vec::with_capacity(2) };

for j in (0..q_len).rev() {
/*
* When calculating our next guess q0, we don't need to consider the digits below j
* + b.data.len() - 1: we're guessing digit j of the quotient (i.e. q0 << j) from
* digit bn of the divisor (i.e. bn << (b.data.len() - 1) - so the product of those
* two numbers will be zero in all digits up to (j + b.data.len() - 1).
*/
let offset = j + b.data.len() - 1;
if offset >= a.data.len() {
continue;
}
n = 1;
// FIXME(#5992): assignment operator overloads
// d += d0;
d = d + d0;
// FIXME(#5992): assignment operator overloads
// m -= prod;
m = m - prod;

/* just avoiding a heap allocation: */
let mut a0 = tmp;
a0.data.truncate(0);
a0.data.extend(&a.data[offset..]);

/*
* q0 << j * big_digit::BITS is our actual quotient estimate - we do the shifts
* implicitly at the end, when adding and subtracting to a and q. Not only do we
* save the cost of the shifts, the rest of the arithmetic gets to work with
* smaller numbers.
*/
let (mut q0, _) = div_rem_digit(a0, bn);
let mut prod = &b * &q0;

while cmp_slice(&prod.data[..], &a.data[j..]) == Greater {
let one: BigUint = One::one();
q0 = q0 - one;
prod = prod - &b;
}

add2(&mut q.data[j..], &q0.data[..]);
sub2(&mut a.data[j..], &prod.data[..]);
a = a.normalize();

tmp = q0;
}
return (d, m);
}


fn div_estimate(a: &BigUint, b: &BigUint, n: usize)
-> (BigUint, BigUint, BigUint) {
if a.data.len() < n {
return (Zero::zero(), Zero::zero(), (*a).clone());
}

let an = &a.data[a.data.len() - n ..];
let bn = *b.data.last().unwrap();
let mut d = Vec::with_capacity(an.len());
let mut carry = 0;
for elt in an.iter().rev() {
let ai = big_digit::to_doublebigdigit(carry, *elt);
let di = ai / (bn as DoubleBigDigit);
assert!(di < big_digit::BASE);
carry = (ai % (bn as DoubleBigDigit)) as BigDigit;
d.push(di as BigDigit)
}
d.reverse();

let shift = (a.data.len() - an.len()) - (b.data.len() - 1);
if shift == 0 {
return (BigUint::new(d), One::one(), (*b).clone());
}
let one: BigUint = One::one();
return (BigUint::new(d).shl_unit(shift),
one.shl_unit(shift),
b.shl_unit(shift));
}

debug_assert!(a < b);

(q.normalize(), a)
}
}

/// Calculates the Greatest Common Divisor (GCD) of the number and `other`.
Expand Down Expand Up @@ -1111,43 +1147,18 @@ fn to_str_radix_reversed(u: &BigUint, radix: u32) -> Vec<u8> {
vec![b'0']
} else {
let mut res = Vec::new();
let mut digits = u.data.to_vec();

while !digits.is_empty() {
let rem = div_rem_in_place(&mut digits, radix);
res.push(to_digit(rem as u8));
let mut digits = u.clone();

// If we finished the most significant digit, drop it
if let Some(&0) = digits.last() {
digits.pop();
}
while digits != Zero::zero() {
let (q, r) = div_rem_digit(digits, radix as BigDigit);
res.push(to_digit(r as u8));
digits = q;
}

res
}
}

fn div_rem_in_place(digits: &mut [BigDigit], divisor: BigDigit) -> BigDigit {
let mut rem = 0;

for d in digits.iter_mut().rev() {
let (q, r) = full_div_rem(*d, divisor, rem);
*d = q;
rem = r;
}

rem
}

fn full_div_rem(a: BigDigit, b: BigDigit, borrow: BigDigit) -> (BigDigit, BigDigit) {
let lo = a as DoubleBigDigit;
let hi = borrow as DoubleBigDigit;

let lhs = lo | (hi << big_digit::BITS);
let rhs = b as DoubleBigDigit;
((lhs / rhs) as BigDigit, (lhs % rhs) as BigDigit)
}

fn to_digit(b: u8) -> u8 {
match b {
0 ... 9 => b'0' + b,
Expand Down

0 comments on commit 7976257

Please sign in to comment.