392 lines
12 KiB
C++
392 lines
12 KiB
C++
#ifndef CLICK_BIGINT_HH
|
|
#define CLICK_BIGINT_HH 1
|
|
/*
|
|
* bigint.hh -- multiple-precision arithmetic
|
|
* Eddie Kohler and the authors of GMP
|
|
*
|
|
* Copyright (c) 2008 Meraki, Inc.
|
|
* Derived from the GNU Multiple Precision Arithmetic Library, which is
|
|
* Copyright (c) 2001-2008 Free Software Foundation, Inc.
|
|
*
|
|
* This source code is free software; you can redistribute it and/or modify it
|
|
* under the terms of the GNU Lesser General Public License as published by
|
|
* the Free Software Foundation, version 3.
|
|
*
|
|
* This program is distributed in the hope that it will be useful, but WITHOUT
|
|
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
|
|
* License for more details.
|
|
*
|
|
* Using the interfaces provided by this file does not create a derived work.
|
|
*/
|
|
#include <click/integers.hh>
|
|
#include <click/type_traits.hh>
|
|
#include <click/string.hh>
|
|
CLICK_DECLS
|
|
|
|
/** @file <click/bigint.hh>
|
|
* @brief Simple multiple-precision integer arithmetic.
|
|
*/
|
|
|
|
/** @class Bigint
|
|
* @brief Template for multiple-precision integer arithmetic.
|
|
*
|
|
* The Bigint template provides a couple simple, but useful, functions for
|
|
* multiple-precision integer arithmetic. It is a port of a very limited
|
|
* subset of the GNU Multiple Precision Arithmetic library's mpn functionality
|
|
* (architecture-generic functions that work on unsigned integers of known
|
|
* size).
|
|
*
|
|
* A multiple-precision integer is an array of "limbs," or machine integers.
|
|
* The multiply_add() and divide() functions take arrays of limbs as
|
|
* arguments. In such arrays, the least significant limb is stored at array
|
|
* index 0.
|
|
*
|
|
* For example, here's Bigint speak for "x = 0x123_4567_89AB_CDEF * 0x45":
|
|
* @code
|
|
* // Initialize result
|
|
* uint32_t x[2] = { 0, 0 };
|
|
* // Define multiplicand; note that least significant comes first
|
|
* uint32_t a[2] = { 0x89ABCDEF, 0x1234567 };
|
|
* Bigint<uint32_t>::multiply_add(x, a, 2, 0x45);
|
|
* @endcode
|
|
*
|
|
* Bigint is not optimized for speed.
|
|
*
|
|
* Bigint's template parameter is the type of a limb. Most users will use
|
|
* bigint, which is equivalent to Bigint<uint32_t>.
|
|
*/
|
|
template <typename L, typename Lhalf = L>
|
|
class Bigint { public:
|
|
|
|
/** @brief Type of a limb, the unit of single arithmetic. */
|
|
typedef L limb_type;
|
|
|
|
/** @brief Type of an arithmetic unit no smaller than one-half limb_type. */
|
|
typedef Lhalf half_limb_type;
|
|
|
|
enum {
|
|
limb_bits = sizeof(limb_type) * 8,
|
|
half_limb_bits = sizeof(limb_type) * 4
|
|
};
|
|
enum {
|
|
limb_high_bit = limb_type(1) << (limb_bits - 1),
|
|
limb_half = limb_type(1) << half_limb_bits,
|
|
limb_low_mask = limb_half - 1,
|
|
limb_zero = limb_type(0)
|
|
};
|
|
|
|
/** @brief Return the less significant half of @a x.
|
|
* @invariant low(@a x) + ((limb_type) high(@a x) << half_limb_bits) ==
|
|
* @a x */
|
|
static inline half_limb_type low(limb_type x) {
|
|
return x & limb_low_mask;
|
|
}
|
|
|
|
/** @brief Return the more significant half of @a x. */
|
|
static inline half_limb_type high(limb_type x) {
|
|
return x >> half_limb_bits;
|
|
}
|
|
|
|
/** @brief Add two-limb integers [@a a1,@a a0] and [@a b1,@a b0], storing
|
|
* the result in [@a x1,@a x0].
|
|
* @param[out] x1 most significant limb of result
|
|
* @param[out] x0 least significant limb of result
|
|
* @param a1 most significant limb of addend
|
|
* @param a0 least significant limb of addend
|
|
* @param b1 most significant limb of addend
|
|
* @param b0 least significant limb of addend
|
|
*
|
|
* Handles carry in the low order limbs. */
|
|
static void add(limb_type &x1, limb_type &x0,
|
|
limb_type a1, limb_type a0, limb_type b1, limb_type b0) {
|
|
x0 = a0 + b0;
|
|
x1 = a1 + b1 + (x0 < a0);
|
|
}
|
|
|
|
/** @brief Multiply one-limb integers @a a and @a b, storing the result in
|
|
* [@a x1,@a x0].
|
|
* @param[out] x1 most significant limb of result
|
|
* @param[out] x0 least significant limb of result
|
|
* @param a multiplicand
|
|
* @param b multiplicand */
|
|
static void multiply(limb_type &x1, limb_type &x0,
|
|
limb_type a, limb_type b) {
|
|
int_multiply(a, b, x0, x1);
|
|
}
|
|
|
|
/** @brief Multiply one-limb integer @a a by half-limb integer @a b,
|
|
* storing the result in [@a x1,@a x0].
|
|
* @param[out] x1 most significant limb of result
|
|
* @param[out] x0 least significant limb of result
|
|
* @param a multiplicand
|
|
* @param b multiplicand */
|
|
static void multiply_half(limb_type &x1, limb_type &x0,
|
|
limb_type a, half_limb_type b) {
|
|
if (has_fast_int_multiply<limb_type>::value)
|
|
multiply(x1, x0, a, b);
|
|
else {
|
|
half_limb_type al = low(a), ah = high(a);
|
|
|
|
limb_type r0 = (limb_type) al * b;
|
|
|
|
limb_type r1 = (limb_type) ah * b;
|
|
limb_type r2 = high(r0) + r1;
|
|
|
|
x1 = (r2 < r1 ? limb_half : limb_zero) + high(r2);
|
|
x0 = (r2 << half_limb_bits) + low(r0);
|
|
}
|
|
}
|
|
|
|
/** @brief Return the inverse of limb @a x.
|
|
*
|
|
* Returns the limb @a y that is the largest limb not larger than
|
|
* (2**(2*limb_bits))/@a x - (2**limb_bits). If this would yield
|
|
* overflow, @a y is the largest possible number (i.e., only ones). */
|
|
static limb_type inverse(limb_type x) {
|
|
limb_type y1 = ~x, y0 = ~limb_zero;
|
|
|
|
// now calculate [y1y0] / x
|
|
limb_type x1 = high(x), x0 = low(x);
|
|
limb_type q1 = y1 / x1;
|
|
limb_type r1 = y1 - q1 * x1;
|
|
limb_type m = q1 * x0;
|
|
r1 = r1 * limb_half | high(y0);
|
|
if (r1 < m) {
|
|
--q1, r1 += x;
|
|
if (r1 >= x && r1 < m)
|
|
--q1, r1 += x;
|
|
}
|
|
r1 -= m;
|
|
|
|
limb_type q0 = r1 / x1;
|
|
limb_type r0 = r1 - q0 * x1;
|
|
m = q0 * x0;
|
|
r0 = r0 * limb_half | low(y0);
|
|
if (r0 < m) {
|
|
--q0, r0 += x;
|
|
if (r0 >= x && r0 < m)
|
|
--q0, r0 += x;
|
|
}
|
|
|
|
return q1 * limb_half | q0;
|
|
}
|
|
|
|
/** @brief Set @a n-limb integer @a a to the value @a b.
|
|
* @param[out] a points to @a n-limb result
|
|
* @param n number of limbs in @a a
|
|
* @param b input value
|
|
* @return carry */
|
|
static click_uintmax_t set(limb_type *a, int n, click_uintmax_t b) {
|
|
while (n > 0) {
|
|
*a++ = b;
|
|
n--;
|
|
b >>= limb_bits;
|
|
}
|
|
return b;
|
|
}
|
|
|
|
/** @brief Multiply @a n-limb integer @a a by 1-limb integer @a b and add
|
|
* the result to @a n-limb integer @a x.
|
|
* @param[in,out] x points to @a n-limb addend and result
|
|
* @param a points to @a n-limb multiplicand
|
|
* @param n number of limbs in @a x and @a a
|
|
* @param b 1-limb multiplicand
|
|
* @return overflow
|
|
*
|
|
* Like @a x += @a a * @a b. Both @a x and @a a must have @a n limbs. It
|
|
* is safe for @a x and @a a to point to exactly the same memory, but they
|
|
* must not otherwise overlap. */
|
|
static limb_type multiply_add(limb_type *x, const limb_type *a, int n,
|
|
limb_type b) {
|
|
limb_type carry = 0;
|
|
do {
|
|
limb_type x0, x1;
|
|
multiply(x1, x0, *a++, b);
|
|
x0 += carry;
|
|
carry = (x0 < carry) + x1;
|
|
x0 += *x;
|
|
carry += (x0 < *x);
|
|
*x++ = x0;
|
|
} while (--n != 0);
|
|
return carry;
|
|
}
|
|
|
|
/** @brief Multiply @a n-limb integer @a a by 1-limb integer @a b, add
|
|
* 1-limb integer @a carry, and store the result in @a n-limb integer @a x.
|
|
* @param[in,out] x points to @a n-limb result
|
|
* @param a points to @a n-limb multiplicand
|
|
* @param n number of limbs in @a x and @a a
|
|
* @param b 1-limb multiplicand
|
|
* @param carry 1-limb initial carry
|
|
* @return overflow
|
|
*
|
|
* Like @a x = (@a a * @a b) + @a carry. Both @a x and @a a must have @a
|
|
* n limbs. It is safe for @a x and @a a to point to exactly the same
|
|
* memory, but they must not otherwise overlap. */
|
|
static limb_type multiply(limb_type *x, const limb_type *a, int n,
|
|
limb_type b, limb_type carry = 0) {
|
|
do {
|
|
limb_type x0, x1;
|
|
multiply(x1, x0, *a++, b);
|
|
x0 += carry;
|
|
carry = (x0 < carry) + x1;
|
|
*x++ = x0;
|
|
} while (--n != 0);
|
|
return carry;
|
|
}
|
|
|
|
/** @brief Multiply @a n-limb integer @a a by 1/2-limb integer @a b, add
|
|
* 1-limb integer @a carry, and store the result in @a n-limb integer @a x.
|
|
* @param[in,out] x points to @a n-limb result
|
|
* @param a points to @a n-limb multiplicand
|
|
* @param n number of limbs in @a x and @a a
|
|
* @param b 1/2-limb multiplicand
|
|
* @param carry 1-limb initial carry
|
|
* @return overflow
|
|
*
|
|
* Like @a x = (@a a * @a b) + @a carry. Both @a x and @a a must have @a
|
|
* n limbs. It is safe for @a x and @a a to point to exactly the same
|
|
* memory, but they must not otherwise overlap. */
|
|
static limb_type multiply_half(limb_type *x, const limb_type *a, int n,
|
|
half_limb_type b, limb_type carry = 0) {
|
|
do {
|
|
limb_type x0, x1;
|
|
multiply_half(x1, x0, *a++, b);
|
|
x0 += carry;
|
|
carry = (x0 < carry) + x1;
|
|
*x++ = x0;
|
|
} while (--n != 0);
|
|
return carry;
|
|
}
|
|
|
|
/** @brief Divide @a n-limb integer @a a by 1-limb integer @a b and store
|
|
* the result in @a n-limb integer @a x.
|
|
* @param[out] x points to @a n-limb result
|
|
* @param a points to @a n-limb dividend
|
|
* @param n number of limbs in @a x and @a a
|
|
* @param b 1-limb divisor
|
|
* @return the remainder
|
|
*
|
|
* Like @a x = @a a / @a b. Both @a x and @a a must have @a n limbs. It
|
|
* is safe for @a x and @a a to point to exactly the same memory, but they
|
|
* must not otherwise overlap. */
|
|
static limb_type divide(limb_type *x, const limb_type *a, int n,
|
|
limb_type b) {
|
|
x += n - 1; /* Make x point at most significant quotient limb */
|
|
a += n - 1; /* Make a point at most significant dividend limb */
|
|
limb_type r = 0;
|
|
|
|
if (b & limb_high_bit) {
|
|
/* High quotient limb is 0 or 1, skip a divide step. */
|
|
r = *a;
|
|
*x = (r >= b);
|
|
r -= (b & -*x);
|
|
--x, --a, --n;
|
|
|
|
/* Multiply-by-inverse, divisor already normalized. */
|
|
limb_type b_inverse = inverse(b);
|
|
while (n > 0) {
|
|
preinverted_divide(*x, r, r, *a, b, b_inverse);
|
|
--x, --a, --n;
|
|
}
|
|
|
|
} else {
|
|
/* Most significant bit of divisor == 0. */
|
|
/* Skip a division if high < divisor (high quotient 0). Testing
|
|
here before normalizing will still skip as often as
|
|
possible. */
|
|
if (*a < b) {
|
|
r = *a;
|
|
*x = 0;
|
|
--x, --a, --n;
|
|
}
|
|
|
|
if (n != 0) {
|
|
int norm = ffs_msb(b) - 1; // number of most significant 0 bits
|
|
b <<= norm;
|
|
r <<= norm;
|
|
limb_type b_inverse = inverse(b);
|
|
|
|
limb_type a1 = *a;
|
|
r |= (a1 >> (limb_bits - norm));
|
|
--a, --n;
|
|
while (n > 0) {
|
|
limb_type a0 = *a;
|
|
preinverted_divide(*x, r, r, ((a1 << norm) | (a0 >> (limb_bits - norm))), b, b_inverse);
|
|
--x, --a, --n;
|
|
a1 = a0;
|
|
}
|
|
preinverted_divide(*x, r, r, a1 << norm, b, b_inverse);
|
|
r >>= norm;
|
|
}
|
|
}
|
|
|
|
return r;
|
|
}
|
|
|
|
/** @brief Return a string representation of @a n-limb integer @a x and
|
|
* set @a x to 0.
|
|
* @param[in,out] x @a n-limb input, set to zero on output
|
|
* @param n number of limbs in @a x
|
|
* @param base base (between 2 and 36)
|
|
* @param uppercase if true, use uppercase letters for digits >= 10 */
|
|
static String unparse_clear(limb_type *x, int n, int base = 10, bool uppercase = false) {
|
|
// need d chars, min d s.t. 10^d >= 2^(sizeof(limb_type) * 8 * n)
|
|
// == min d s.t. d >= sizeof(limb_type) * 8 * n / lg 10
|
|
int div = (base >= 16 ? 4 : (base >= 8 ? 3 : 1));
|
|
String s = String::make_uninitialized((n * limb_bits) / div + 1);
|
|
char *q = const_cast<char *>(s.end());
|
|
assert(base >= 2 && base <= 36);
|
|
while (1) {
|
|
while (n > 0 && x[n - 1] == 0)
|
|
--n;
|
|
if (n == 0)
|
|
break;
|
|
int r = divide(x, x, n, base);
|
|
if (r <= 9)
|
|
*--q = '0' + r;
|
|
else if (uppercase)
|
|
*--q = 'A' + r - 10;
|
|
else
|
|
*--q = 'a' + r - 10;
|
|
}
|
|
assert(q >= s.begin());
|
|
if (q == s.end())
|
|
*--q = '0';
|
|
return s.substring(q, s.end());
|
|
}
|
|
|
|
private:
|
|
|
|
/* Divide the two-limb number in (a1,a0) by b, with b_inverse being the
|
|
largest limb not larger than (2**(2*limb_bits))/b - (2**limb_bits). If
|
|
this would yield overflow, b_inverse should be the largest possible
|
|
number (i.e., only ones). For correct operation, the most significant
|
|
bit of b has to be set. Put the quotient in q and the remainder in
|
|
r. */
|
|
/* Like udiv_qrnnd_preinv, but branch-free. */
|
|
static void preinverted_divide(limb_type &q, limb_type &r,
|
|
limb_type a1, limb_type a0, limb_type b,
|
|
limb_type b_inverse) {
|
|
limb_type a0_mask = (a0 & limb_high_bit ? ~limb_zero : limb_zero);
|
|
limb_type a0_adjusted = a0 + (b & a0_mask);
|
|
limb_type x1, x0;
|
|
multiply(x1, x0, b_inverse, a1 - a0_mask);
|
|
add(x1, x0, x1, x0, a1, a0_adjusted);
|
|
limb_type q1 = ~x1;
|
|
multiply(x1, x0, q1, b);
|
|
add(x1, x0, x1, x0, a1, a0);
|
|
x1 -= b;
|
|
r = x0 + (b & x1);
|
|
q = x1 - q1;
|
|
}
|
|
|
|
};
|
|
|
|
/** @brief Typical Bigint usage with uint32_t limb_type. */
|
|
typedef Bigint<uint32_t> bigint;
|
|
|
|
CLICK_ENDDECLS
|
|
#endif
|