|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r49471 - in sandbox/mp_math: boost/mp_math/mp_int boost/mp_math/mp_int/detail libs/mp_math/test
From: baraclese_at_[hidden]
Date: 2008-10-27 20:50:19
Author: baraclese
Date: 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
New Revision: 49471
URL: http://svn.boost.org/trac/boost/changeset/49471
Log:
- large restructuring and overhaul of modular reduction / modpow code
- added file modpow.hpp, moved modpow stuff into it from pow.hpp
- turned most modular reduction functions non-member functions
- added modpow_ctx to hold precalculated variables for modpows against the same modulus
- turned mp_int<>::pow into a non-member function
- moved modpow tests from pow.cpp into modpow.cpp, added some new tests
- renamed primitive_ops<>::add_magnitude to add_smaller_magnitude
- fixed primitive_ops<>::add_smaller_magnitude
- added small optimization to primitive_ops<>::sub_smaller_magnitude
- prime.hpp : use new modpow_ctx
- added push, pop and is_unitialized member functions to mp_int<>
- added include of <boost/cstdint.hpp> to prerequisite.hpp and guarded use of uint64_t with a macro
- misc comment and style fixes
Added:
sandbox/mp_math/boost/mp_math/mp_int/modpow.hpp (contents, props changed)
sandbox/mp_math/libs/mp_math/test/modpow.cpp (contents, props changed)
Text files modified:
sandbox/mp_math/boost/mp_math/mp_int/detail/primitive_ops.hpp | 32 +
sandbox/mp_math/boost/mp_math/mp_int/mod.hpp | 2
sandbox/mp_math/boost/mp_math/mp_int/modinv.hpp | 29 +
sandbox/mp_math/boost/mp_math/mp_int/modular_reduction.hpp | 579 ++++++++++++++-------------------------
sandbox/mp_math/boost/mp_math/mp_int/mp_int.hpp | 34 -
sandbox/mp_math/boost/mp_math/mp_int/pow.hpp | 484 +--------------------------------
sandbox/mp_math/boost/mp_math/mp_int/prime.hpp | 14
sandbox/mp_math/boost/mp_math/mp_int/sub.hpp | 2
sandbox/mp_math/libs/mp_math/test/jamfile.v2 | 1
sandbox/mp_math/libs/mp_math/test/pow.cpp | 78 ----
sandbox/mp_math/libs/mp_math/test/prerequisite.hpp | 3
11 files changed, 290 insertions(+), 968 deletions(-)
Modified: sandbox/mp_math/boost/mp_math/mp_int/detail/primitive_ops.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/detail/primitive_ops.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/detail/primitive_ops.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -47,10 +47,10 @@
digit_type& c);
// z = x + y, where xlen >= ylen
- // z and x may not overlap
- static void add_magnitude(digit_type* z,
- const digit_type* x, size_type xlen,
- const digit_type* y, size_type ylen);
+ // returns last carry
+ static digit_type add_smaller_magnitude(digit_type* z,
+ const digit_type* x, size_type xlen,
+ const digit_type* y, size_type ylen);
// SUB ------------------------------------
@@ -172,21 +172,22 @@
template<typename D, typename W, typename S>
inline
-void basic_primitive_ops<D,W,S>::add_magnitude(
+typename basic_primitive_ops<D,W,S>::digit_type
+basic_primitive_ops<D,W,S>::add_smaller_magnitude(
digit_type* z,
const digit_type* x, size_type xlen,
const digit_type* y, size_type ylen)
{
digit_type carry = add_digits(z, x, y, ylen);
-
+
size_type n = ripple_carry(z + ylen, x + ylen, xlen - ylen, carry);
- if (carry)
- *(z + n++) = carry; // wrong
+ n += ylen;
+
+ if (n < xlen && z != x)
+ std::memcpy(z + n, x + n, sizeof(digit_type) * (xlen - n));
- const size_type cur = ylen + n;
- if (cur < xlen)
- std::memcpy(z + cur, x + cur, (xlen - cur) * sizeof(digit_type));
+ return carry;
}
template<typename D, typename W, typename S>
@@ -257,10 +258,13 @@
{
const digit_type borrow = subtract_digits(z, x, y, ylen);
- const size_type n = ripple_borrow(z + ylen, x + ylen, xlen - ylen, borrow);
+ size_type n = ripple_borrow(z + ylen, x + ylen, xlen - ylen, borrow);
- const size_type cur = ylen + n;
- std::memcpy(z + cur, x + cur, (xlen - cur) * sizeof(digit_type));
+ if (z != x)
+ {
+ n += ylen;
+ std::memcpy(z + n, x + n, (xlen - n) * sizeof(digit_type));
+ }
}
Modified: sandbox/mp_math/boost/mp_math/mp_int/mod.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/mod.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/mod.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -22,5 +22,5 @@
clamp();
if (is_zero())
sign_ = 1;
- }
+}
Modified: sandbox/mp_math/boost/mp_math/mp_int/modinv.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/modinv.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/modinv.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -5,25 +5,26 @@
// hac 14.61, pp608
-// *this = *this**-1 (mod b)
+// Computes the modular multiplicative inverse of *this mod m
+// *this = *this**-1 (mod m)
template<class A, class T>
-void mp_int<A,T>::modinv(const mp_int& b)
+void mp_int<A,T>::modinv(const mp_int& m)
{
- if (b.is_negative() || !b)
+ if (m.is_negative() || !m)
throw std::domain_error("modinv: modulus is negative or zero");
// if the modulus is odd we can use a faster routine
- if (b.is_odd())
- odd_modinv(b);
+ if (m.is_odd())
+ odd_modinv(m);
else
- even_modinv(b);
+ even_modinv(m);
}
-/* hac 14.61, pp608 */
+// hac 14.61, pp608
template<class A1, class T>
void mp_int<A1,T>::even_modinv(const mp_int& y)
{
- assert(y.is_positive() && y);
+ assert(y.is_even());
static const char* const err_msg = "mp_int::modinv: inverse does not exist";
@@ -46,7 +47,7 @@
if (A.is_odd() || B.is_odd())
{
- /* A = (A+y)/2, B = (B-x)/2 */
+ // A = (A+y)/2, B = (B-x)/2
A += y;
B -= x;
}
@@ -60,7 +61,7 @@
if (C.is_odd() || D.is_odd())
{
- /* C = (C+y)/2, D = (D-x)/2 */
+ // C = (C+y)/2, D = (D-x)/2
C += y;
D -= x;
}
@@ -84,9 +85,9 @@
if (u)
goto top;
- /* now a = C, b = D, gcd == g*v */
+ // now a = C, b = D, gcd == g*v
- /* if v != 1 then there is no inverse */
+ // if v != 1 then there is no inverse
if (v != digit_type(1))
throw std::domain_error(err_msg);
@@ -98,7 +99,7 @@
while (C.compare_magnitude(y) != -1)
C -= y;
- swap(C);
+ *this = C;
}
/* computes the modular inverse via binary extended euclidean algorithm,
@@ -168,7 +169,7 @@
while (D.is_negative())
D += x;
- swap(D);
+ *this = D;
}
Added: sandbox/mp_math/boost/mp_math/mp_int/modpow.hpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/boost/mp_math/mp_int/modpow.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -0,0 +1,541 @@
+// Copyright Kevin Sopp 2008.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+
+// r = x mod m given x and m
+template<class A, class T>
+struct modpow_ctx
+{
+ typedef typename mp_int<A,T>::digit_type digit_type;
+ typedef typename mp_int<A,T>::word_type word_type;
+ typedef typename mp_int<A,T>::size_type size_type;
+
+ // dr means diminished radix
+ enum modulus_type_t
+ { // R is our radix, i.e. digit_max
+ mod_restricted_dr, // m = R**k - d; d <= R
+ mod_unrestricted_dr, // m = 2**k - d; d <= R
+ mod_unrestricted_dr_slow, // m = 2**k - d; d < d**(k/2)
+ mod_odd,
+ mod_generic
+ }modulus_type;
+
+ modpow_ctx() : precalculated(false){}
+
+ modulus_type_t do_detect(const mp_int<A,T>& m) const;
+
+ void detect_modulus_type(const mp_int<A,T>& m)
+ {
+ modulus_type = do_detect(m);
+ }
+
+ void precalculate(const mp_int<A,T>& m);
+
+ mp_int<A,T> mu;
+ digit_type rho;
+ bool precalculated;
+};
+
+
+template<class A, class T>
+typename modpow_ctx<A,T>::modulus_type_t
+modpow_ctx<A,T>::do_detect(const mp_int<A,T>& m) const
+{
+ if (m.size() == 1)
+ return mod_unrestricted_dr;
+
+ typename mp_int<A,T>::size_type count = 0;
+
+ const int bits = m.precision() % mp_int<A,T>::valid_bits;
+
+ if (!bits && m[m.size()-1] == mp_int<A,T>::digit_max)
+ ++count;
+
+ for (typename mp_int<A,T>::const_reverse_iterator d = m.rbegin() + 1;
+ d != m.rend(); ++d)
+ {
+ if (*d != mp_int<A,T>::digit_max)
+ break;
+ else
+ ++count;
+ }
+
+ // if all bits are set
+ if (count == m.size() - 1)
+ return mod_restricted_dr;
+
+ // if all bits until the most significant digit are set
+ if (count == m.size() - 2)
+ {
+ bool all_bits_set = true;
+
+ // handle the remaining bits in the most significant digit
+ typename mp_int<A,T>::digit_type mask = 1;
+ for (int i = 0; i < bits; ++i)
+ {
+ if ((m[m.size()-1] & mask) == 0)
+ {
+ all_bits_set = false;
+ break;
+ }
+ mask <<= 1;
+ }
+ if (all_bits_set)
+ return mod_unrestricted_dr;
+ }
+
+ // if more than half of the bits are set
+ if (count >= m.size() / 2)
+ return mod_unrestricted_dr_slow;
+
+ if (m.is_odd())
+ return mod_odd;
+
+ return mod_generic;
+}
+
+template<class A, class T>
+void modpow_ctx<A,T>::precalculate(const mp_int<A,T>& m)
+{
+ typedef typename mp_int<A,T>::digit_type digit_type;
+ typedef typename mp_int<A,T>::word_type word_type;
+
+ switch (modulus_type)
+ {
+ case mod_restricted_dr:
+ {
+ rho = (word_type(1) << static_cast<word_type>(mp_int<A,T>::valid_bits))
+ - static_cast<word_type>(m[0]);
+ break;
+ }
+ case mod_unrestricted_dr:
+ {
+ const size_type p = m.precision();
+
+ mp_int<A,T> tmp;
+ tmp.pow2(p);
+ tmp.sub_smaller_magnitude(m);
+
+ rho = tmp[0];
+ break;
+ }
+ case mod_unrestricted_dr_slow:
+ {
+ mp_int<A,T> tmp;
+
+ tmp.pow2(m.precision());
+ mu = tmp - m;
+ break;
+ }
+ case mod_odd:
+ {
+ assert(m.is_odd());
+
+ /* fast inversion mod 2**k
+ *
+ * Based on the fact that
+ *
+ * XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
+ * => 2*X*A - X*X*A*A = 1
+ * => 2*(1) - (1) = 1
+ */
+ const digit_type b = m[0];
+
+ static const typename mp_int<A,T>::size_type S =
+ sizeof(digit_type) * std::numeric_limits<unsigned char>::digits;
+
+ digit_type x = (((b + 2) & 4) << 1) + b; // here x*a==1 mod 2**4
+ x *= 2 - b * x; // here x*a==1 mod 2**8
+ if (S != 8)
+ x *= 2 - b * x; // here x*a==1 mod 2**16
+ if (S == 64 || !(S == 8 || S == 16))
+ x *= 2 - b * x; // here x*a==1 mod 2**32
+ if (S == 64)
+ x *= 2 - b * x; // here x*a==1 mod 2**64
+
+ // rho = -1/m mod b
+ rho = (word_type(1) << (static_cast<word_type>(mp_int<A,T>::valid_bits))) - x;
+ break;
+ }
+ case mod_generic:
+ {
+ // mu = b**2k/m
+ mu.pow2(m.size() * 2 * mp_int<A,T>::digit_bits);
+ mu /= m;
+ break;
+ }
+ }
+ precalculated = true;
+}
+
+
+template<class A, class T>
+mp_int<A,T> barret_modpow(const mp_int<A,T>& base,
+ const mp_int<A,T>& exp,
+ const mp_int<A,T>& m,
+ int redmode,
+ mp_int<A,T>& mu)
+{
+ typedef typename mp_int<A,T>::size_type size_type;
+ typedef typename mp_int<A,T>::size_type digit_type;
+
+ void (*redux)(mp_int<A,T>&, const mp_int<A,T>&, const mp_int<A,T>&);
+
+ // find window size
+ const size_type x = exp.precision();
+ size_type winsize;
+ if (x <= 7)
+ winsize = 2;
+ else if (x <= 36)
+ winsize = 3;
+ else if (x <= 140)
+ winsize = 4;
+ else if (x <= 450)
+ winsize = 5;
+ else if (x <= 1303)
+ winsize = 6;
+ else if (x <= 3529)
+ winsize = 7;
+ else
+ winsize = 8;
+
+ if (redmode == 0)
+ redux = &barrett_reduce;
+ else
+ redux = &unrestricted_dr_slow_reduce;
+
+ // The M table contains powers of the base, e.g. M[i] = base**i % m
+ // The first half of the table is not computed though except for M[0] and M[1]
+ mp_int<A,T> M[256];
+ M[1] = base % m;
+
+ // compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times
+ M[1 << (winsize - 1)] = M[1];
+
+ for (size_type i = 0; i < (winsize - 1); ++i)
+ {
+ const size_type offset = 1 << (winsize - 1);
+ M[offset].sqr();
+ // reduce modulo m
+ (*redux)(M[offset], m, mu);
+ }
+
+ // create upper table, that is M[i] = M[i-1] * M[1] (mod P)
+ // for i = (2**(winsize - 1) + 1) to (2**winsize - 1)
+ for (size_type i = (1 << (winsize - 1)) + 1;
+ i < (size_type(1) << winsize);
+ ++i)
+ {
+ M[i] = M[i-1] * M[1];
+ (*redux)(M[i], m, mu);
+ }
+
+ // setup result
+ mp_int<A,T> res = digit_type(1);
+
+ int mode = 0;
+ int bitcnt = 1;
+ digit_type buf = 0;
+ typename mp_int<A,T>::const_reverse_iterator exp_digit = exp.rbegin();
+ size_type bitcpy = 0;
+ int bitbuf = 0;
+
+ for (;;) // while (exp_digit != exp.rend())
+ {
+ // grab next digit as required
+ if (--bitcnt == 0)
+ {
+ if (exp_digit == exp.rend())
+ break;
+
+ // read next digit and reset the bitcnt
+ buf = *exp_digit++;
+ bitcnt = mp_int<A,T>::valid_bits;
+ }
+
+ // grab the next msb from the exponent
+ int y = (buf >> static_cast<digit_type>(mp_int<A,T>::valid_bits - 1)) & 1;
+ buf <<= digit_type(1);
+
+ // if the bit is zero and mode == 0 then we ignore it
+ // These represent the leading zero bits before the first 1 bit
+ // in the exponent. Technically this opt is not required but it
+ // does lower the # of trivial squaring/reductions used
+ if (mode == 0 && y == 0)
+ continue;
+
+ // if the bit is zero and mode == 1 then we square
+ if (mode == 1 && y == 0)
+ {
+ res.sqr();
+ (*redux)(res, m, mu);
+ continue;
+ }
+
+ // else we add it to the window
+ bitbuf |= (y << (winsize - ++bitcpy));
+ mode = 2;
+
+ if (bitcpy == winsize)
+ {
+ // ok window is filled so square as required and multiply
+ // square first
+ for (size_type i = 0; i < winsize; ++i)
+ {
+ res.sqr();
+ (*redux)(res, m, mu);
+ }
+
+ // then multiply
+ res *= M[bitbuf];
+ (*redux)(res, m, mu);
+
+ // empty window and reset
+ bitcpy = 0;
+ bitbuf = 0;
+ mode = 1;
+ }
+ }
+
+ // if bits remain then square/multiply
+ if (mode == 2 && bitcpy > 0)
+ {
+ // square then multiply if the bit is set
+ for (size_type i = 0; i < bitcpy; ++i)
+ {
+ res.sqr();
+ (*redux)(res, m, mu);
+
+ bitbuf <<= 1;
+ if ((bitbuf & (1 << winsize)) != 0)
+ {
+ res *= M[1];
+ (*redux)(res, m, mu);
+ }
+ }
+ }
+
+ return res;
+}
+
+/* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
+ *
+ * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
+ * The value of k changes based on the size of the exponent.
+ *
+ * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
+ */
+template<class A, class T>
+mp_int<A,T> montgomery_modpow(const mp_int<A,T>& base,
+ const mp_int<A,T>& exp,
+ const mp_int<A,T>& m,
+ int redmode,
+ typename mp_int<A,T>::digit_type mp)
+{
+ typedef typename mp_int<A,T>::size_type size_type;
+ typedef typename mp_int<A,T>::digit_type digit_type;
+
+ void (*redux)(mp_int<A,T>&, const mp_int<A,T>&, digit_type);
+
+ // find window size
+ const size_type x = exp.precision();
+ size_type winsize;
+ if (x <= 7)
+ winsize = 2;
+ else if (x <= 36)
+ winsize = 3;
+ else if (x <= 140)
+ winsize = 4;
+ else if (x <= 450)
+ winsize = 5;
+ else if (x <= 1303)
+ winsize = 6;
+ else if (x <= 3529)
+ winsize = 7;
+ else
+ winsize = 8;
+
+ //digit_type mp = ctx->mp;
+
+ if (redmode == 0)
+ redux = &montgomery_reduce;
+ else if (redmode == 1)
+ redux = &restricted_dr_reduce;
+ else
+ redux = &unrestricted_dr_reduce;
+
+ mp_int<A,T> res;
+
+ // create M table
+ // The first half of the table is not computed though except for M[0] and M[1]
+ mp_int<A,T> M[256];
+
+ if (redmode == 0)
+ {
+ // now we need R mod m
+ montgomery_normalize(res, m);
+ // now set M[1] to G * R mod m
+ M[1] = (base * res) % m;
+ }
+ else
+ {
+ res = digit_type(1);
+ M[1] = base % m;
+ }
+
+ // compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times
+ M[1 << (winsize - 1)] = M[1];
+
+ for (size_type i = 0; i < (winsize - 1); ++i)
+ {
+ const size_type offset = 1 << (winsize - 1);
+ M[offset].sqr();
+ (*redux)(M[offset], m, mp);
+ }
+
+ // create upper table
+ for (size_type i = (1 << (winsize - 1)) + 1;
+ i < (size_type(1) << winsize); ++i)
+ {
+ M[i] = M[i-1] * M[1];
+ (*redux)(M[i], m, mp);
+ }
+
+ /* set initial mode and bit cnt */
+ int mode = 0;
+ int bitcnt = 1;
+ digit_type buf = 0;
+ typename mp_int<A,T>::const_reverse_iterator exp_digit = exp.rbegin();
+ size_type bitcpy = 0;
+ unsigned int bitbuf = 0;
+
+ for (;;)
+ {
+ // grab next digit as required
+ if (--bitcnt == 0)
+ {
+ if (exp_digit == exp.rend())
+ break;
+ // read next digit and reset bitcnt
+ buf = *exp_digit++;
+ bitcnt = mp_int<A,T>::valid_bits;
+ }
+
+ // grab the next msb from the exponent
+ int y = (buf >> (mp_int<A,T>::valid_bits - 1)) & 1;
+ buf <<= digit_type(1);
+
+ // if the bit is zero and mode == 0 then we ignore it
+ // These represent the leading zero bits before the first 1 bit
+ // in the exponent. Technically this opt is not required but it
+ // does lower the # of trivial squaring/reductions used
+ if (mode == 0 && y == 0)
+ continue;
+
+ // if the bit is zero and mode == 1 then we square
+ if (mode == 1 && y == 0)
+ {
+ res.sqr();
+ (*redux)(res, m, mp);
+ continue;
+ }
+
+ // else we add it to the window
+ bitbuf |= (y << (winsize - ++bitcpy));
+ mode = 2;
+
+ if (bitcpy == winsize)
+ {
+ // ok window is filled so square as required and multiply
+ // square first
+ for (size_type i = 0; i < winsize; ++i)
+ {
+ res.sqr();
+ (*redux)(res, m, mp);
+ }
+
+ // then multiply
+ res *= M[bitbuf];
+ (*redux)(res, m, mp);
+
+ // empty window and reset
+ bitcpy = 0;
+ bitbuf = 0;
+ mode = 1;
+ }
+ }
+
+ // if bits remain then square/multiply
+ if (mode == 2 && bitcpy > 0)
+ {
+ // square then multiply if the bit is set
+ for (size_type i = 0; i < bitcpy; ++i)
+ {
+ res.sqr();
+ (*redux)(res, m, mp);
+
+ // get next bit of the window
+ bitbuf <<= 1;
+ if ((bitbuf & (1 << winsize)) != 0)
+ {
+ // then multiply
+ res *= M[1];
+ (*redux)(res, m, mp);
+ }
+ }
+ }
+
+ // Fixup result if Montgomery reduction is used. Recall that any value in a
+ // Montgomery system is actually multiplied by R mod n. So we have to reduce
+ // one more time to cancel out the factor of R.
+ if (redmode == 0)
+ (*redux)(res, m, mp);
+
+ return res;
+}
+
+
+// z = base^exp % mod
+template<class A, class T>
+mp_int<A,T> modpow(const mp_int<A,T>& base,
+ const mp_int<A,T>& exp,
+ const mp_int<A,T>& mod,
+ modpow_ctx<A,T>* ctx = 0)
+{
+ if (mod.is_negative())
+ throw std::domain_error("modpow: modulus must be positive");
+
+ typedef modpow_ctx<A,T> ctx_t;
+
+ ctx_t tmp_ctx;
+
+ if (!ctx)
+ {
+ tmp_ctx.detect_modulus_type(mod);
+ tmp_ctx.precalculate(mod);
+ ctx = &tmp_ctx;
+ }
+ else
+ {
+ if (!ctx->precalculated)
+ {
+ ctx->detect_modulus_type(mod);
+ ctx->precalculate(mod);
+ }
+ }
+
+ if (exp.is_negative())
+ return modpow(modinv(base, mod), abs(exp), mod, ctx);
+
+ switch (ctx->modulus_type)
+ {
+ case ctx_t::mod_restricted_dr: return montgomery_modpow(base, exp, mod, 1, ctx->rho);
+ case ctx_t::mod_unrestricted_dr: return montgomery_modpow(base, exp, mod, 2, ctx->rho);
+ case ctx_t::mod_unrestricted_dr_slow: return barret_modpow (base, exp, mod, 1, ctx->mu);
+ case ctx_t::mod_odd: return montgomery_modpow(base, exp, mod, 0, ctx->rho);
+ case ctx_t::mod_generic:
+ default: return barret_modpow (base, exp, mod, 0, ctx->mu);
+ }
+}
Modified: sandbox/mp_math/boost/mp_math/mp_int/modular_reduction.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/modular_reduction.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/modular_reduction.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -3,207 +3,147 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
-/* pre-calculate the value required for Barrett reduction
- * For a given modulus "b" it calulates the value required in "a"
- */
-template<class A, class T>
-void mp_int<A,T>::reduce_setup(const mp_int& b)
-{
- pow2(b.used_ * 2 * digit_bits);
- *this /= b;
-}
-// reduces *this mod m, assumes 0 < *this < m**2, mu is
-// precomputed via reduce_setup.
+// reduces x mod m, assumes 0 < x < m**2, mu is precomputed.
// From HAC pp.604 Algorithm 14.42
template<class A, class T>
-void mp_int<A,T>::reduce(const mp_int& m, const mp_int& mu)
+void barrett_reduce(mp_int<A,T>& x, const mp_int<A,T>& m, const mp_int<A,T>& mu)
{
- const size_type k = m.used_;
+ typedef typename mp_int<A,T>::digit_type digit_type;
- mp_int q(*this);
+ const typename mp_int<A,T>::size_type k = m.size();
- /* q1 = x / b**(k-1) */
+ mp_int<A,T> q(x);
+
+ // q = x / base**(k-1)
q.shift_digits_right(k - 1);
- /* according to HAC this optimization is ok */
- if (k > digit_type(1) << (valid_bits - 1))
+ // according to HAC this optimization is ok
+ if (k > digit_type(1) << (mp_int<A,T>::valid_bits - 1))
q *= mu;
else
q.fast_mul_high_digits(mu, k);
- /* q3 = q2 / b**(k+1) */
+ // q = q / base**(k+1)
q.shift_digits_right(k + 1);
- /* x = x mod b**(k+1), quick (no division) */
- modulo_2_to_the_power_of(valid_bits * (k + 1));
+ // r = x mod base**(k+1)
+ x.modulo_2_to_the_power_of(mp_int<A,T>::valid_bits * (k + 1));
- /* q = q * m mod b**(k+1), quick (no division) */
+ // q = q * m mod base**(k+1)
q.mul_digits(m, k + 1);
- /* x = x - q */
- *this -= q;
+ x -= q;
- /* If x < 0, add b**(k+1) to it */
- if (is_negative())
+ // If x < 0, add base**(k+1) to it
+ if (x.is_negative())
{
q = digit_type(1);
q.shift_digits_left(k + 1);
- *this += q;
+ x += q;
}
- /* Back off if it's too big */
- while (compare(m) != -1)
- sub_smaller_magnitude(m);
-}
-
-/* setups the montgomery reduction stuff */
-template<class A, class T>
-typename mp_int<A,T>::digit_type
-mp_int<A,T>::montgomery_setup() const
-{
- /* fast inversion mod 2**k
- *
- * Based on the fact that
- *
- * XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
- * => 2*X*A - X*X*A*A = 1
- * => 2*(1) - (1) = 1
- */
- const digit_type b = digits_[0];
-
- if (is_even())
- throw std::domain_error("montgomery_setup: integer must be odd");
-
- static const size_type S = sizeof(digit_type) * CHAR_BIT;
-
- digit_type x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
- x *= 2 - b * x; /* here x*a==1 mod 2**8 */
- if (S != 8)
- x *= 2 - b * x; /* here x*a==1 mod 2**16 */
- if (S == 64 || !(S == 8 || S == 16))
- x *= 2 - b * x; /* here x*a==1 mod 2**32 */
- if (S == 64)
- x *= 2 - b * x; /* here x*a==1 mod 2**64 */
-
- /* rho = -1/m mod b */
- const digit_type rho = (word_type(1) << (static_cast<word_type>(valid_bits))) - x;
- return rho;
+ while (x >= m)
+ x.sub_smaller_magnitude(m);
}
/* computes xR**-1 == x (mod N) via Montgomery Reduction */
template<class A, class T>
-void mp_int<A,T>::montgomery_reduce(const mp_int& n, digit_type rho)
-{
- /* can the fast reduction [comba] method be used?
- *
- * Note that unlike in mul you're safely allowed *less*
- * than the available columns [255 per default] since carries
- * are fixed up in the inner loop.
- */
- const size_type digs = n.used_ * 2 + 1;
-/* const size_type ws = sizeof(word_type) * CHAR_BIT;
- if (digs < mp_warray && n.used_ < 1 << (ws - (2 * valid_bits)))
- {
- fast_montgomery_reduce(n, rho);
- return;
- }*/
- /* grow the input as required */
- grow_capacity(digs);
- std::memset(digits_ + used_, 0, (capacity_ - used_) * sizeof(digit_type));
- used_ = digs;
-
- for (size_type i = 0; i < n.used_; ++i)
- {
- /* mu = ai * rho mod b
- *
- * The value of rho must be precalculated via
- * montgomery_setup() such that
- * it equals -1/n0 mod b this allows the
- * following inner loop to reduce the
- * input one digit at a time
- */
- const digit_type mu = static_cast<word_type>(digits_[i]) * rho;
+void montgomery_reduce(mp_int<A,T>& x,
+ const mp_int<A,T>& n,
+ typename mp_int<A,T>::digit_type rho)
+{
+ typedef typename mp_int<A,T>::digit_type digit_type;
+ typedef typename mp_int<A,T>::word_type word_type;
+ typedef typename mp_int<A,T>::size_type size_type;
+
+ const size_type digs = n.size() * 2 + 1;
+
+ x.grow_capacity(digs);
+ std::memset(x.digits() + x.size(), 0, (x.capacity() - x.size()) * sizeof(digit_type));
+ x.set_size(digs);
+
+ // TODO rewrite both for loops in terms of const_iterator nd = n.begin();...
+ for (size_type i = 0; i < n.size(); ++i)
+ {
+ // mu = x[i] * rho (mod base)
+ // The value of rho must be precalculated such that it equals -1/n0 mod b
+ // this allows the following inner loop to reduce the input one digit at a
+ // time.
+ const digit_type mu = static_cast<word_type>(x[i]) * rho;
- /* a = a + mu * m * b**i */
+ // x = x + mu * m * base**i
- /* alias for digits of the modulus */
- digit_type* tmpn = n.digits_;
+ // alias for digits of the modulus
+ const digit_type* tmpn = n.digits();
- /* alias for the digits of x [the input] */
- digit_type* tmpx = digits_ + i;
+ // alias for the digits of x
+ digit_type* tmpx = x.digits() + i;
digit_type carry = 0;
- /* Multiply and add in place */
- for (size_type j = 0; j < n.used_; ++j)
+ // Multiply and add in place
+ for (size_type j = 0; j < n.size(); ++j)
{
- /* compute product and sum */
+ // compute product and sum
const word_type r = static_cast<word_type>(mu)
* static_cast<word_type>(*tmpn++)
+ static_cast<word_type>(carry)
+ static_cast<word_type>(*tmpx);
- carry = r >> digit_bits;
+ carry = static_cast<digit_type>(r >> mp_int<A,T>::digit_bits);
- /* fix digit */
*tmpx++ = r;
}
- /* At this point the i'th digit of x should be zero */
+ // At this point the i'th digit of x should be zero
- /* propagate carries upwards as required*/
+ // propagate carries upwards as required
while (carry)
{
const word_type r = static_cast<word_type>(*tmpx) + carry;
*tmpx++ = r;
- carry = r >> digit_bits;
+ carry = r >> mp_int<A,T>::digit_bits;
}
}
- /* at this point the n.used'th least
- * significant digits of x are all zero
- * which means we can shift x to the
- * right by n.used digits and the
- * residue is unchanged.
- */
-
- /* x = x/b**n.used */
- clamp();
- if (is_zero())
- sign_ = 1;
- shift_digits_right(n.used_);
-
- /* if x >= n then x = x - n */
- if (compare_magnitude(n) != -1)
- sub_smaller_magnitude(n);
+ // at this point the n.used'th least significant digits of x are all zero
+ // which means we can shift x to the right by n.used digits and the residue is
+ // unchanged.
+
+ // x = x/base**n.size()
+ x.clamp();
+ if (x.is_zero())
+ x.set_sign(1);
+
+ x.shift_digits_right(n.size());
+
+ if (x.compare_magnitude(n) != -1)
+ x.sub_smaller_magnitude(n);
}
-/*
- * shifts with subtractions when the result is greater than b.
- *
- * The method is slightly modified to shift B unconditionally upto just under
- * the leading bit of b. This saves alot of multiple precision shifting.
- */
+// shifts with subtractions when the result is greater than n.
+// The method is slightly modified to shift B unconditionally upto just under
+// the leading bit of n. This saves alot of multiple precision shifting.
template<class A, class T>
-void mp_int<A,T>::montgomery_calc_normalization(const mp_int& b)
+void montgomery_normalize(mp_int<A,T>& x, const mp_int<A,T>& n)
{
- /* how many bits of last digit does b use */
- size_type bits = b.precision() % valid_bits;
+ // how many bits of last digit does n use
+ typename mp_int<A,T>::size_type bits = n.precision() % mp_int<A,T>::valid_bits;
- if (b.used_ > 1)
- pow2((b.used_ - 1) * valid_bits + bits - 1);
+ if (n.size() > 1)
+ x.pow2((n.size() - 1) * mp_int<A,T>::valid_bits + bits - 1);
else
{
- *this = digit_type(1);
+ x = typename mp_int<A,T>::digit_type(1);
bits = 1;
}
- /* now compute C = A * B mod b */
- for (int x = bits - 1; x < valid_bits; ++x)
+ // now compute C = A * B mod n
+ for (int i = bits - 1; i < mp_int<A,T>::valid_bits; ++i)
{
- multiply_by_2();
- if (compare_magnitude(b) != -1)
- sub_smaller_magnitude(b);
+ x.multiply_by_2();
+ if (x.compare_magnitude(n) != -1)
+ x.sub_smaller_magnitude(n);
}
}
@@ -216,36 +156,35 @@
* Based on Algorithm 14.32 on pp.601 of HAC.
*/
template<class A, class T>
-void mp_int<A,T>::fast_montgomery_reduce(const mp_int& n, digit_type rho)
-{
- word_type W[mp_warray];
-
- /* grow a as required */
- grow_capacity(n.used_ + 1);
-
- /* first we have to get the digits of the input into
- * an array of double precision words W[...]
- */
-
- /* copy the digits of a into W[0..a->used-1] */
- for (size_type i = 0; i < used_; ++i)
- W[i] = digits_[i];
-
- /* zero the high words of W[a->used..m->used*2] */
- std::memset(W + used_, 0, (n.used_ * 2 + 1) * sizeof(word_type));
-
- /* now we proceed to zero successive digits
- * from the least significant upwards
- */
- for (size_type i = 0; i < n.used_; ++i)
- {
- /* mu = ai * m' mod b
- *
- * We avoid a double precision multiplication (which isn't required)
- * by casting the value down to a mp_digit. Note this requires
- * that W[ix-1] have the carry cleared (see after the inner loop)
- */
- const digit_type mu = ((W[i] & mp_mask) * rho) & mp_mask;
+void fast_montgomery_reduce(mp_int<A,T>& x,
+ const mp_int<A,T>& n,
+ typename mp_int<A,T>::digit_type rho)
+{
+ typedef typename mp_int<A,T>::digit_type digit_type;
+ typedef typename mp_int<A,T>::word_type word_type;
+ typedef typename mp_int<A,T>::size_type size_type;
+
+ word_type W[512];
+
+ x.grow_capacity(n.used_ + 1);
+
+ for (size_type i = 0; i < x.size(); ++i)
+ W[i] = x[i];
+
+ // zero the high words of W[a->used..m->used*2]
+ std::memset(W + x.size(), 0, (n.size() * 2 + 1) * sizeof(word_type));
+
+ // now we proceed to zero successive digits
+ // from the least significant upwards
+ for (size_type i = 0; i < n.size(); ++i)
+ {
+ // mu = ai * m' mod b
+ //
+ // We avoid a double precision multiplication (which isn't required) by
+ // casting the value down to a mp_digit. Note this requires that W[ix-1]
+ // have the carry cleared (see after the inner loop)
+ const digit_type mu = static_cast<digit_type>(W[i])
+ * static_cast<word_type>(rho);
/* a = a + mu * m * b**i
*
@@ -262,260 +201,148 @@
* first m->used words of W[] have the carries fixed
*/
/* inner loop */
- for (size_type j = 0; j < n.used_; ++j)
- W[i+j] += static_cast<word_type>(mu) * static_cast<word_type>(n.digits_[j]);
+ for (size_type j = 0; j < n.size(); ++j)
+ W[i+j] += static_cast<word_type>(mu) * static_cast<word_type>(n[j]);
- /* now fix carry for next digit, W[ix+1] */
- W[i + 1] += W[i] >> static_cast<word_type>(valid_bits);
+ // now fix carry for next digit, W[ix+1]
+ W[i + 1] += W[i] >> static_cast<word_type>(mp_int<A,T>::valid_bits);
}
- /* now we have to propagate the carries and
- * shift the words downward [all those least
- * significant digits we zeroed].
- */
-
- /* now fix rest of carries */
-
- for (size_type i = n.used_ + 1; i <= n.used_ * 2 + 1; ++i)
- W[i] += W[i-1] >> static_cast<word_type>(valid_bits);
-
- /* copy out, A = A/b**n
- *
- * The result is A/b**n but instead of converting from an
- * array of mp_word to mp_digit than calling mp_rshd
- * we just copy them in the right order
- */
-
- for (size_type i = 0; i < n.used_ + 1; ++i)
- digits_[i] = W[n.used_ + i] & static_cast<word_type>(mp_mask);
-
- /* set the max used and clamp */
- used_ = n.used_ + 1;
- clamp();
- if (is_zero())
- sign_ = 1;
-
- /* if A >= m then A = A - m */
- if (compare_magnitude(n) != -1)
- sub_smaller_magnitude(n);
+ // now we have to propagate the carries and shift the words downward [all
+ // those least significant digits we zeroed].
+
+ // now fix rest of carries
+ for (size_type i = n.size() + 1; i <= n.size() * 2 + 1; ++i)
+ W[i] += W[i-1] >> static_cast<word_type>(mp_int<A,T>::valid_bits);
+
+ // copy out, A = A/b**n
+ //
+ // The result is A/b**n but instead of converting from an array of word_type
+ // to digit_type than calling mp_rshd we just copy them in the right order
+ for (size_type i = 0; i < n.size() + 1; ++i)
+ x[i] = static_cast<digit_type>(W[n.size() + i]);
+
+ // set the max used and clamp
+ x.set_size(n.size() + 1);
+ x.clamp();
+ if (x.is_zero())
+ x.set_sign(1);
+
+ // if A >= m then A = A - m
+ if (x.compare_magnitude(n) != -1)
+ x.sub_smaller_magnitude(n);
}
-/* determines the setup value */
-template<class A, class T>
-typename mp_int<A,T>::digit_type
-mp_int<A,T>::dr_setup() const
-{
- /* the casts are required if valid_bits is one less than
- * the number of bits in a valid_bits [e.g. valid_bits==31]
- */
- const digit_type d = (word_type(1) << static_cast<word_type>(valid_bits))
- - static_cast<word_type>(digits_[0]);
- return d;
-}
-/* reduce "x" in place modulo "n" using the Diminished Radix algorithm.
- *
- * Based on algorithm from the paper
- *
- * "Generating Efficient Primes for Discrete Log Cryptosystems"
- * Chae Hoon Lim, Pil Joong Lee,
- * POSTECH Information Research Laboratories
- *
- * The modulus must be of a special format [see manual]
- *
- * Has been modified to use algorithm 7.10 from the LTM book instead
- *
- * Input *this must be in the range 0 <= *this <= (n-1)**2
- */
-template<class A, class T>
-void mp_int<A,T>::dr_reduce(const mp_int& n, digit_type k)
-{
- /* m = digits in modulus */
- const size_type m = n.used_;
+// reduce "x" modulo "n" using the Diminished Radix algorithm.
+// Based on algorithm from the paper
+//
+// "Generating Efficient Primes for Discrete Log Cryptosystems"
+// Chae Hoon Lim, Pil Joong Lee,
+// POSTECH Information Research Laboratories
+//
+// The modulus must be of a special format [see manual]
+//
+// Has been modified to use algorithm 7.10 from the LTM book instead
+//
+// Input x must be in the range 0 <= x <= (n-1)**2
+template<class A, class T>
+void restricted_dr_reduce(mp_int<A,T>& x,
+ const mp_int<A,T>& n,
+ typename mp_int<A,T>::digit_type k)
+{
+ typedef typename mp_int<A,T>::digit_type digit_type;
+ typedef typename mp_int<A,T>::word_type word_type;
+ typedef typename mp_int<A,T>::size_type size_type;
+
+ const size_type m = n.size();
+
+ x.grow_capacity(m + m);
+ std::memset(x.digits() + x.size(), 0, (m + m - x.size()) * sizeof(digit_type));
- /* ensure that *this has at least 2m digits */
- grow_capacity(m + m);
- std::memset(digits_ + used_, 0, (m + m - used_) * sizeof(digit_type));
-
-/* top of loop, this is where the code resumes if
- * another reduction pass is required.
- */
top:
- /* set carry to zero */
+ // set carry to zero
digit_type mu = 0;
- /* compute (*this mod B**m) + k * [*this/B**m] inline and inplace */
+ // compute (r mod B**m) + k * [r/B**m] inline and inplace
for (size_type i = 0; i < m; ++i)
{
- const word_type r = static_cast<word_type>(digits_[m+i])
- * static_cast<word_type>(k) + digits_[i] + mu;
- digits_[i] = static_cast<digit_type>(r);
- mu = static_cast<digit_type>(r >> static_cast<word_type>(digit_bits));
+ const word_type r = static_cast<word_type>(x[m+i])
+ * static_cast<word_type>(k) + x[i] + mu;
+ x[i] = static_cast<digit_type>(r);
+ mu = static_cast<digit_type>(r >> static_cast<word_type>(mp_int<A,T>::digit_bits));
}
- /* set final carry */
- digits_[m] = mu;
+ // set final carry
+ x[m] = mu;
- /* zero words above m */
- if (used_ > m + 1) // guard against overflow
- std::memset(digits_ + m + 1, 0, (used_ - (m + 1)) * sizeof(digit_type));
-
- /* clamp, sub and return */
- clamp();
- if (is_zero())
- sign_ = 1;
-
- /* if *this >= n then subtract and reduce again
- * Each successive "recursion" makes the input smaller and smaller.
- */
- if (compare_magnitude(n) != -1)
- {
- sub_smaller_magnitude(n);
- goto top;
- }
-}
+ // zero words above m
+ if (x.size() > m + 1) // guard against overflow
+ std::memset(x.digits() + m + 1, 0, (x.size() - (m + 1)) * sizeof(digit_type));
-/* determines if a number is a valid DR modulus */
-template<class A, class T>
-bool mp_int<A,T>::is_dr_modulus() const
-{
- /* must be at least two digits */
- if (used_ < 2)
- return false;
-
- /* must be of the form b**k - a [a <= b] so all
- * but the first digit must be equal to -1 (mod b).
- */
- for (size_type i = 1; i < used_; ++i)
+ x.clamp();
+ if (x.is_zero())
+ x.set_sign(1);
+
+ if (x.compare_magnitude(n) != -1)
{
- if (digits_[i] != digit_max)
- return false;
+ x.sub_smaller_magnitude(n);
+ goto top;
}
- return true;
-}
-
-/* determines the setup value */
-template<class A, class T>
-typename mp_int<A,T>::digit_type
-mp_int<A,T>::reduce_2k_setup() const
-{
- mp_int tmp;
- const size_type p = precision();
-
- tmp.pow2(p);
- tmp.sub_smaller_magnitude(*this);
-
- return tmp.digits_[0];
}
-/* reduces *this modulo n where n is of the form 2**p - d */
+// reduces x modulo n where n is of the form 2**p - d
template<class A, class T>
-void mp_int<A,T>::reduce_2k(const mp_int& n, digit_type d)
+void unrestricted_dr_reduce(mp_int<A,T>& x,
+ const mp_int<A,T>& n,
+ typename mp_int<A,T>::digit_type d)
{
- const size_type p = n.precision();
+ const typename mp_int<A,T>::size_type p = n.precision();
top:
- mp_int q(*this);
+ mp_int<A,T> q(x);
- /* q = a/2**p, a = a mod 2**p */
- q.shift_right(p, this);
+ /* q = a/2**p, r = r mod 2**p */
+ q.shift_right(p, &x);
if (d != 1)
- /* q = q * d */
q.multiply_by_digit(d);
- /* a = a + q */
- add_magnitude(q);
+
+ x.add_magnitude(q);
- if (compare_magnitude(n) != -1)
+ if (x.compare_magnitude(n) != -1)
{
- sub_smaller_magnitude(n);
+ x.sub_smaller_magnitude(n);
goto top;
}
}
-/* determines if mp_reduce_2k can be used */
+// reduces x modulo n where n is of the form 2**p - d. This differs from
+// unrestricted_dr_reduce since "d" can be larger than a single digit.
template<class A, class T>
-bool mp_int<A,T>::reduce_is_2k() const
+void unrestricted_dr_slow_reduce(mp_int<A,T>& x,
+ const mp_int<A,T>& n,
+ const mp_int<A,T>& d)
{
- if (used_ == 1)
- return true;
- else if (used_ > 1)
- {
- const size_type bits = precision();
- size_type j = 1;
- digit_type k = 1;
-
- /* Test every bit from the second digit up, must be 1 */
- for (size_type i = valid_bits; i < bits; ++i)
- {
- if ((digits_[j] & k) == 0)
- return false;
- k <<= 1;
- if (k > mp_mask) // FIXME this can never happen with the new code
- {
- ++j;
- k = 1;
- }
- }
- }
- return true;
-}
-
-/* determines the setup value */
-template<class A, class T>
-mp_int<A,T> mp_int<A,T>::reduce_2k_l_setup()
-{
- mp_int tmp;
-
- tmp.pow2(precision());
- return tmp - *this;
-}
-
-/* reduces *this modulo n where n is of the form 2**p - d
- This differs from reduce_2k since "d" can be larger
- than a single digit.
-*/
-template<class A, class T>
-void mp_int<A,T>::reduce_2k_l(const mp_int& n, const mp_int& d)
-{
- const size_type p = n.precision();
+ const typename mp_int<A,T>::size_type p = n.precision();
top:
- mp_int q(*this);
+ mp_int<A,T> q(x);
- /* q = a/2**p, a = a mod 2**p */
- q.shift_right(p, this);
+ // q = x/2**p, r = r mod 2**p
+ q.shift_right(p, &x);
- /* q = q * d */
q *= d;
- /* a = a + q */
- add_magnitude(q);
+ x.add_magnitude(q);
- if (compare_magnitude(n) != -1)
+ if (x.compare_magnitude(n) != -1)
{
- sub_smaller_magnitude(n);
+ x.sub_smaller_magnitude(n);
goto top;
}
}
-/* determines if reduce_2k_l can be used */
-template<class A, class T>
-bool mp_int<A,T>::reduce_is_2k_l() const
-{
- if (used_ == 1)
- return true;
- else if (used_ > 1)
- {
- size_type count = 0;
- // if more than half of the digits are -1 we're sold
- for (size_type i = 0; i < used_; ++i)
- if (digits_[i] == digit_max)
- ++count;
- return count >= used_ / 2 ? true : false;
- }
- return false;
-}
Modified: sandbox/mp_math/boost/mp_math/mp_int/mp_int.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/mp_int.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/mp_int.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -33,6 +33,7 @@
namespace boost {
namespace mp_math {
+template<class A, class T> struct modpow_ctx;
// digits are stored in least significant order
@@ -259,12 +260,17 @@
return digits_[i];
}
+ void push(digit_type x) { digits_[used_++] = x; }
+ void pop() { --used_; }
+
void zero();
// debug functionality
void print(bool all=false) const;
bool test_invariants() const;
+ bool is_uninitialized() const { return digits_ == 0; }
+
size_type size() const { return used_; }
size_type capacity() const { return capacity_; }
@@ -316,33 +322,10 @@
size_type precision() const;
size_type count_lsb() const;
void shift_right(size_type b, mp_int* remainder);
-
- void reduce_setup(const mp_int&);
- void reduce(const mp_int& m, const mp_int& mu);
-
- digit_type montgomery_setup() const;
- void montgomery_reduce(const mp_int& n, digit_type rho);
- void montgomery_calc_normalization(const mp_int&);
- void fast_montgomery_reduce(const mp_int& n, digit_type rho);
-
- void dr_reduce(const mp_int& n, digit_type k);
- digit_type dr_setup() const;
- bool is_dr_modulus() const;
-
- digit_type reduce_2k_setup() const;
- void reduce_2k(const mp_int& n, digit_type d);
- bool reduce_is_2k() const;
-
- mp_int reduce_2k_l_setup();
- void reduce_2k_l(const mp_int& n, const mp_int& d);
- bool reduce_is_2k_l() const;
-
+
void pow2(size_type b);
- void pow(digit_type);
- void modpow(const mp_int& exp, const mp_int& m);
- void barret_modpow(const mp_int& exp, const mp_int& m, int reduction_mode);
- void fast_modpow(const mp_int& exp, const mp_int& m, int reduction_mode);
+ void modpow(const mp_int& exp, const mp_int& m, modpow_ctx<Allocator,Traits>* c = 0);
void modinv(const mp_int& modulus);
void even_modinv(const mp_int& modulus);
@@ -875,6 +858,7 @@
#include <boost/mp_math/mp_int/mod.hpp>
#include <boost/mp_math/mp_int/modinv.hpp>
#include <boost/mp_math/mp_int/modular_reduction.hpp>
+#include <boost/mp_math/mp_int/modpow.hpp>
#include <boost/mp_math/mp_int/mul.hpp>
#include <boost/mp_math/mp_int/operators.hpp>
#include <boost/mp_math/mp_int/pow.hpp>
Modified: sandbox/mp_math/boost/mp_math/mp_int/pow.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/pow.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/pow.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -3,494 +3,46 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
-/* computes a = 2**b
- *
- * Simple algorithm which zeroes the int, grows it then just sets one bit
- * as required.
- */
+// computes a = 2**b
+// Simple algorithm which zeroes the int, grows it then just sets one bit
+// as required.
template<class A, class T>
void mp_int<A,T>::pow2(typename mp_int<A,T>::size_type b)
{
- /* grow a to accomodate the single bit */
grow_capacity(b / digit_bits + 1);
- /* set the used count of where the bit will go */
+ // set the used count of where the bit will go
used_ = b / digit_bits + 1;
// set all bits to zero
std::memset(digits_, 0, used_ * sizeof(digit_type));
- /* put the single bit in its place */
+ // put the single bit in its place
digits_[b / digit_bits] = digit_type(1) << (b % digit_bits);
}
-/* calculate c = a**b using a square-multiply algorithm */
+// calculate c = x**y using a square-multiply algorithm
template<class A, class T>
-void mp_int<A,T>::pow(digit_type b)
+mp_int<A,T> pow(const mp_int<A,T>& x, typename mp_int<A,T>::digit_type y)
{
- mp_int g(*this);
+ mp_int<A,T> result;
- /* set initial result */
- *this = digit_type(1);
+ result = typename mp_int<A,T>::digit_type(1);
- for (int x = 0; x < digit_bits; ++x)
- {
- sqr();
-
- /* if the bit is set multiply */
- if (b & (digit_type(1) << (digit_bits - 1)))
- *this *= g;
-
- /* shift to next bit */
- b <<= 1;
- }
-}
-
-// this is a shell function that calls either the normal or Montgomery
-// modpow functions. Originally the call to the montgomery code was
-// embedded in the normal function but that wasted alot of stack space
-// for nothing (since 99% of the time the Montgomery code would be called)
-//
-// y = g^x mod p
-// *this = *this^exp mod m
-template<class A, class T>
-void mp_int<A,T>::modpow(const mp_int& exp, const mp_int& m)
-{
- if (m.is_negative())
- throw std::domain_error("modpow: modulus m must be positive");
-
- if (exp.is_negative())
- {
- // first compute 1 / *this mod m
- modinv(m);
-
- // and now compute (1 / *this)**|exp| instead of *this**exp [exp < 0]
- modpow(abs(exp), m);
- return;
- }
-
- /* modified diminished radix reduction */
- if (m.reduce_is_2k_l())
- {
- //std::cout << "barret_modpow1\n";
- barret_modpow(exp, m, 1);
- return;
- }
-
- /* is it a DR modulus? */
- int dr = m.is_dr_modulus();
-
- /* if not, is it a unrestricted DR modulus? */
- if (!dr)
- dr = m.reduce_is_2k() << 1;
-
- /* if the modulus is odd or dr == true use the montgomery method */
- if (m.is_odd() || dr)
- {
- //std::cout << "fast_modpow dr = " << dr << "\n";
- fast_modpow(exp, m, dr);
- }
- else
- {
- //std::cout << "barret_modpow2\n";
- /* otherwise use the generic Barrett reduction technique */
- barret_modpow(exp, m, 0);
- }
-}
-
-template<class A, class T>
-void mp_int<A,T>::barret_modpow(const mp_int& exp, const mp_int& m, int redmode)
-{
- static const int TAB_SIZE = 256;
- mp_int M[TAB_SIZE];
-
- void (mp_int::*redux)(const mp_int&, const mp_int&);
-
- /* find window size */
- size_type x = exp.precision();
- size_type winsize;
- if (x <= 7)
- winsize = 2;
- else if (x <= 36)
- winsize = 3;
- else if (x <= 140)
- winsize = 4;
- else if (x <= 450)
- winsize = 5;
- else if (x <= 1303)
- winsize = 6;
- else if (x <= 3529)
- winsize = 7;
- else
- winsize = 8;
-
- /* create mu, used for Barrett reduction */
- mp_int mu;
-
- mp_int mm(m);
- if (redmode == 0)
- {
- mu.reduce_setup(m);
- //std::cout << "redux = reduce\n";
- redux = &mp_int::reduce;
- }
- else
- {
- mu = mm.reduce_2k_l_setup();
- //std::cout << "redux = reduce_2k_l\n";
- redux = &mp_int::reduce_2k_l;
- }
-
- /* create M table
- *
- * The M table contains powers of the base,
- * e.g. M[x] = G**x mod P
- *
- * The first half of the table is not
- * computed though accept for M[0] and M[1]
- */
- M[1] = *this % mm;
-
-
- /* compute the value at M[1<<(winsize-1)] by squaring
- * M[1] (winsize-1) times
- */
- M[1 << (winsize - 1)] = M[1];
-
- for (x = 0; x < (winsize - 1); x++)
- {
- /* square it */
- M[1 << (winsize - 1)].sqr();
- /* reduce modulo P */
- (M[1 << (winsize - 1)].*redux)(mm, mu);
- }
-
- /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
- * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
- */
- for (x = (1 << (winsize - 1)) + 1; x < (size_type(1) << winsize); ++x)
- {
- M[x] = M[x-1] * M[1];
- (M[x].*redux)(mm, mu);
- }
-
- /* setup result */
- mp_int res = digit_type(1);
-
- /* set initial mode and bit cnt */
- int mode = 0;
- int bitcnt = 1;
- digit_type buf = 0;
- int digidx = exp.used_ - 1;
- size_type bitcpy = 0;
- int bitbuf = 0;
-
- for (;;)
- {
- /* grab next digit as required */
- if (--bitcnt == 0)
- {
- /* if digidx == -1 we are out of digits */
- if (digidx == -1)
- break;
- /* read next digit and reset the bitcnt */
- buf = exp.digits_[digidx--];
- bitcnt = valid_bits;
- }
-
- /* grab the next msb from the exponent */
- int y = (buf >> static_cast<digit_type>(valid_bits - 1)) & 1;
- buf <<= digit_type(1);
-
- /* if the bit is zero and mode == 0 then we ignore it
- * These represent the leading zero bits before the first 1 bit
- * in the exponent. Technically this opt is not required but it
- * does lower the # of trivial squaring/reductions used
- */
- if (mode == 0 && y == 0)
- continue;
-
- /* if the bit is zero and mode == 1 then we square */
- if (mode == 1 && y == 0)
- {
- res.sqr();
- (res.*redux)(mm, mu);
- continue;
- }
-
- /* else we add it to the window */
- bitbuf |= (y << (winsize - ++bitcpy));
- mode = 2;
-
- if (bitcpy == winsize)
- {
- /* ok window is filled so square as required and multiply */
- /* square first */
- for (x = 0; x < winsize; ++x)
- {
- res.sqr();
- (res.*redux)(mm, mu);
- }
-
- /* then multiply */
- res *= M[bitbuf];
- (res.*redux)(mm, mu);
-
- /* empty window and reset */
- bitcpy = 0;
- bitbuf = 0;
- mode = 1;
- }
- }
-
- /* if bits remain then square/multiply */
- if (mode == 2 && bitcpy > 0)
- {
- /* square then multiply if the bit is set */
- for (x = 0; x < bitcpy; ++x)
- {
- res.sqr();
- (res.*redux)(mm, mu);
- //mu = res * mm;
-
- bitbuf <<= 1;
- if ((bitbuf & (1 << winsize)) != 0)
- {
- /* then multiply */
- res *= M[1];
- (res.*redux)(mm, mu);
- }
- }
- }
-
- swap(res);
-}
-
-/* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
- *
- * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
- * The value of k changes based on the size of the exponent.
- *
- * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
- */
-template<class A, class T>
-void mp_int<A,T>::fast_modpow(const mp_int& exp, const mp_int& m, int redmode)
-{
- static const int TAB_SIZE = 256;
- mp_int M[TAB_SIZE];
-
- /* use a pointer to the reduction algorithm. This allows us to use
- * one of many reduction algorithms without modding the guts of
- * the code with if statements everywhere.
- */
- void (mp_int::*redux)(const mp_int&, digit_type);
+ const typename mp_int<A,T>::digit_type mask = 1 << (mp_int<A,T>::digit_bits - 1);
- /* find window size */
- size_type x = exp.precision();
- size_type winsize;
- if (x <= 7)
- winsize = 2;
- else if (x <= 36)
- winsize = 3;
- else if (x <= 140)
- winsize = 4;
- else if (x <= 450)
- winsize = 5;
- else if (x <= 1303)
- winsize = 6;
- else if (x <= 3529)
- winsize = 7;
- else
- winsize = 8;
-
- digit_type mp;
- /* determine and setup reduction code */
- if (redmode == 0)
- {
- /* now setup montgomery */
- mp = m.montgomery_setup();
-
- /* automatically pick the comba one if available (saves quite a few calls/ifs) */
- /*if (((m.used_ * 2 + 1) < mp_warray) &&
- m.used_ < 1 << ((CHAR_BIT * sizeof(word_type)) - (2 * valid_bits)))
- {
- std::cout << "redux = fast_montgomery_reduce\n";
- redux = &mp_int::fast_montgomery_reduce;
- }
- else
- {*/
- /* use slower baseline Montgomery method */
- //std::cout << "redux = montgomery_reduce\n";
- redux = &mp_int::montgomery_reduce;
- //}
- }
- else if (redmode == 1)
- {
- /* setup DR reduction for moduli of the form B**k - b */
- mp = m.dr_setup();
- redux = &mp_int::dr_reduce;
- }
- else
- {
- /* setup DR reduction for moduli of the form 2**k - b */
- mp = m.reduce_2k_setup();
- redux = &mp_int::reduce_2k;
- }
-
- /* setup result */
- mp_int res;
-
- /* create M table
- *
- * The first half of the table is not computed though accept for M[0] and M[1]
- */
- if (redmode == 0)
- {
- /* now we need R mod m */
- res.montgomery_calc_normalization(m);
- /* now set M[1] to G * R mod m */
- M[1] = (*this * res) % m;
- }
- else
+ for (int i = 0; i < mp_int<A,T>::digit_bits; ++i)
{
- res = digit_type(1);
- M[1] = *this % m;
- }
+ result.sqr();
- /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
- M[1 << (winsize - 1)] = M[1];
+ // if the bit is set multiply
+ if (y & mask)
+ result *= x;
- for (x = 0; x < (winsize - 1); ++x)
- {
- M[1 << (winsize - 1)].sqr();
- (M[1 << (winsize - 1)].*redux)(m, mp);
+ // shift to next bit
+ y <<= 1;
}
- /* create upper table */
- for (x = (1 << (winsize - 1)) + 1; x < (size_type(1) << winsize); ++x)
- {
- M[x] = M[x-1] * M[1];
- (M[x].*redux)(m, mp);
- }
-
- /* set initial mode and bit cnt */
- int mode = 0;
- int bitcnt = 1;
- digit_type buf = 0;
- int digidx = exp.used_ - 1;
- size_type bitcpy = 0;
- unsigned int bitbuf = 0;
-
- for (;;)
- {
- /* grab next digit as required */
- if (--bitcnt == 0)
- {
- /* if digidx == -1 we are out of digits so break */
- if (digidx == -1)
- break;
- /* read next digit and reset bitcnt */
- buf = exp.digits_[digidx--];
- bitcnt = valid_bits;
- }
-
- /* grab the next msb from the exponent */
- int y = (buf >> (valid_bits - 1)) & 1;
- buf <<= digit_type(1);
-
- /* if the bit is zero and mode == 0 then we ignore it
- * These represent the leading zero bits before the first 1 bit
- * in the exponent. Technically this opt is not required but it
- * does lower the # of trivial squaring/reductions used
- */
- if (mode == 0 && y == 0)
- continue;
-
- /* if the bit is zero and mode == 1 then we square */
- if (mode == 1 && y == 0)
- {
- res.sqr();
- (res.*redux)(m, mp);
- continue;
- }
-
- /* else we add it to the window */
- bitbuf |= (y << (winsize - ++bitcpy));
- mode = 2;
-
- if (bitcpy == winsize)
- {
- /* ok window is filled so square as required and multiply */
- /* square first */
- for (x = 0; x < winsize; ++x)
- {
- res.sqr();
- (res.*redux)(m, mp);
- }
-
- /* then multiply */
- res *= M[bitbuf];
- (res.*redux)(m, mp);
-
- /* empty window and reset */
- bitcpy = 0;
- bitbuf = 0;
- mode = 1;
- }
- }
-
- /* if bits remain then square/multiply */
- if (mode == 2 && bitcpy > 0)
- {
- /* square then multiply if the bit is set */
- for (x = 0; x < bitcpy; ++x)
- {
- res.sqr();
- (res.*redux)(m, mp);
-
- /* get next bit of the window */
- bitbuf <<= 1;
- if ((bitbuf & (1 << winsize)) != 0)
- {
- /* then multiply */
- res *= M[1];
- (res.*redux)(m, mp);
- }
- }
- }
-
- /* fixup result if Montgomery reduction is used
- * recall that any value in a Montgomery system is
- * actually multiplied by R mod n. So we have
- * to reduce one more time to cancel out the factor
- * of R.
- */
- if (redmode == 0)
- (res.*redux)(m, mp);
-
- swap(res);
+ return result;
}
-
-
-template<class A, class T>
-inline mp_int<A,T> pow(const mp_int<A,T>& x, typename mp_int<A,T>::digit_type n)
-{
- mp_int<A,T> tmp(x);
- tmp.pow(n);
- return tmp;
-}
-
-template<class A, class T>
-inline mp_int<A,T> modpow(const mp_int<A,T>& base,
- const mp_int<A,T>& exp,
- const mp_int<A,T>& mod)
-{
- mp_int<A,T> tmp(base);
- tmp.modpow(exp, mod);
- return tmp;
-}
-
-
-
-
Modified: sandbox/mp_math/boost/mp_math/mp_int/prime.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/prime.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/prime.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -100,14 +100,15 @@
const mp_int<A,T> p1(p-one);
distribution_type rng(digit_type(2), p1);
-
+
+ modpow_ctx<A,T> ctx;
for (unsigned int i = 0; i < rounds_; ++i)
{
mp_int<A,T> base = rng(eng);
base.set_least_significant_bit(); // test only odd bases
- const mp_int<A,T> tmp = modpow(base, p1, p);
+ const mp_int<A,T> tmp = modpow(base, p1, p, &ctx);
if (tmp != one)
return false; // definitely composite!
@@ -193,11 +194,12 @@
const typename mp_int<A,T>::digit_type one = 1;
const mp_int<A,T> p_minus_one(p-one);
distribution_type rng(one, p_minus_one);
-
+
+ modpow_ctx<A,T> ctx;
for (unsigned int i = 0; i < r; ++i)
{
const mp_int<A,T> base = rng(eng);
- mp_int<A,T> y = modpow(base,n,p);
+ mp_int<A,T> y = modpow(base, n, p, &ctx);
for (typename mp_int<A,T>::size_type j = 0; j < s; ++j)
{
@@ -323,7 +325,9 @@
p.multiply_by_2();
++p;
} while (!is_prime(p, test_));
- } while (p.precision() > bits_); // catch extremely rare case
+ // Catch extremely rare case, this occurs if the carry from ++p ripples
+ // through the whole number thereby adding one more bit to it.
+ } while (p.precision() > bits_);
return p;
}
Modified: sandbox/mp_math/boost/mp_math/mp_int/sub.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/sub.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/sub.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -21,7 +21,7 @@
digits_[0] = b - digits_[0];
sign_ = -1;
}
- else // example 8 - 7 = 1 or 5 - 5 = 0
+ else // example: 8 - 7 = 1 or 5 - 5 = 0
digits_[0] -= b;
}
else
Modified: sandbox/mp_math/libs/mp_math/test/jamfile.v2
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/jamfile.v2 (original)
+++ sandbox/mp_math/libs/mp_math/test/jamfile.v2 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -30,6 +30,7 @@
unit-test jacobi : jacobi.cpp boost_test ;
unit-test lcm : lcm.cpp boost_test ;
unit-test modinv : modinv.cpp boost_test ;
+unit-test modpow : modpow.cpp boost_test ;
unit-test mul : mul.cpp boost_test ;
unit-test pow : pow.cpp boost_test ;
unit-test prime : prime.cpp boost_test ;
Added: sandbox/mp_math/libs/mp_math/test/modpow.cpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/libs/mp_math/test/modpow.cpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -0,0 +1,103 @@
+// Copyright Kevin Sopp 2008.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/test/unit_test.hpp>
+#include "prerequisite.hpp"
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modpow1, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("2");
+ const mp_int_type exp("14");
+ const mp_int_type m("8");
+ const mp_int_type z = modpow(x, exp, m);
+ BOOST_CHECK_EQUAL(z, "0");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modpow2, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("4");
+ const mp_int_type exp("13");
+ const mp_int_type m("497");
+ const mp_int_type z = modpow(x, exp, m);
+ BOOST_CHECK_EQUAL(z, "445");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modpow3, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("2395422");
+ const mp_int_type exp("2424832");
+ const mp_int_type m("2424833");
+ const mp_int_type z = modpow(x, exp, m);
+ BOOST_CHECK_EQUAL(z, "1");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modpow4, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("184");
+ const mp_int_type exp("560");
+ const mp_int_type m("561");
+ const mp_int_type z = modpow(x, exp, m);
+ BOOST_CHECK_EQUAL(z, "1");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modpow5, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("997028168093060821869770104094480850560519901475");
+ const mp_int_type exp("7455602825647884208337395736200454918783366342656");
+ const mp_int_type m("7455602825647884208337395736200454918783366342657");
+ const mp_int_type z = modpow(x, exp, m);
+ BOOST_CHECK_EQUAL(z, "1");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modpow6, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("184");
+ const mp_int_type exp("5600");
+ const mp_int_type m("2668");
+ const mp_int_type z = modpow(x, exp, m);
+ BOOST_CHECK_EQUAL(z, "552");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modpow7, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("0x201abcff00aaffffffffffffffffffffffffffffffffffffffffff"
+ "ffffffffd6d7");
+ const mp_int_type exp("0x123456789abcdef01");
+ // a modulus of type unrestricted diminished radix (2^253 - 41);
+ const mp_int_type m("0x1fffffffffffffffffffffffffffffffffffffffffffffffffffff"
+ "ffffffffd7");
+ const mp_int_type z = modpow(x, exp, m);
+ BOOST_CHECK_EQUAL(z, "0x8f112a89871984cde410bb05621d5a6073557d2da0444b6681699"
+ "80b5ef825a");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modpow8, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("0x201abcff00aaffffffffffffffffffffffffffffffffffffffffff"
+ "ffffffffd6d7123456789abcdef01123456789abcdef01ffffffff");
+ const mp_int_type exp("0x123456789abcdef018978979899");
+ // a modulus of type unrestricted diminished radix slow (2^503 - exp)
+ const mp_int_type m("0x7fffffffffffffffffffffffffffffffffffffffffffffffffffff"
+ "fffffffffffffffffffffffffffffffffffffffffffffedcba987654"
+ "3210fe7687686767");
+ const mp_int_type z = modpow(x, exp, m);
+ BOOST_CHECK_EQUAL(z, "0x1ada35751ae1f5fec7ab0e60f6c924d5f4a4a5d8f6786cdd78838"
+ "5ab16ebd994ee9aaea5faef3f490822ef443fd3e169caa3b608162e"
+ "01d40a593e775c9fa5");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modpow9, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("0x201abcff00aaffffffffffffffffffffffffffffffffffffffffff"
+ "ffffffffd6d7123456789abcdef01123456789abcdef01ffffffff");
+ const mp_int_type exp("0x123456789abcdef018978979899");
+ // a modulus of type restricted diminished radix (2^256 - 53)
+ const mp_int_type m("0xffffffffffffffffffffffffffffffffffffffffffffffffffffff"
+ "ffffffffcb");
+ const mp_int_type z = modpow(x, exp, m);
+ BOOST_CHECK_EQUAL(z, "0x777b5d9b290fbb5f99e4668cf1b0f723d3228fc252da492c54b75"
+ "8379f3024e4");
+}
+
Modified: sandbox/mp_math/libs/mp_math/test/pow.cpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/pow.cpp (original)
+++ sandbox/mp_math/libs/mp_math/test/pow.cpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -27,78 +27,24 @@
BOOST_CHECK_EQUAL(x, "18446744073709551616");
}
-BOOST_AUTO_TEST_CASE_TEMPLATE(pow_1, mp_int_type, mp_int_types)
+BOOST_AUTO_TEST_CASE_TEMPLATE(pow1, mp_int_type, mp_int_types)
{
- mp_int_type x("2");
- x.pow(0);
- BOOST_CHECK_EQUAL(x, "1");
+ const mp_int_type x("2");
+ const mp_int_type z = pow(x, 0);
+ BOOST_CHECK_EQUAL(z, "1");
}
-BOOST_AUTO_TEST_CASE_TEMPLATE(pow_2, mp_int_type, mp_int_types)
+BOOST_AUTO_TEST_CASE_TEMPLATE(pow2, mp_int_type, mp_int_types)
{
- mp_int_type x("2");
- x.pow(1);
- BOOST_CHECK_EQUAL(x, "2");
+ const mp_int_type x("2");
+ const mp_int_type z = pow(x, 1);
+ BOOST_CHECK_EQUAL(z, "2");
}
-BOOST_AUTO_TEST_CASE_TEMPLATE(pow_3, mp_int_type, mp_int_types)
+BOOST_AUTO_TEST_CASE_TEMPLATE(pow3, mp_int_type, mp_int_types)
{
- mp_int_type x("2");
- x.pow(64);
- BOOST_CHECK_EQUAL(x, "18446744073709551616");
-}
-
-BOOST_AUTO_TEST_CASE_TEMPLATE(modpow1, mp_int_type, mp_int_types)
-{
- mp_int_type x("2");
- const mp_int_type exp("14");
- const mp_int_type m("8");
- x.modpow(exp, m);
- BOOST_CHECK_EQUAL(x, "0");
-}
-
-BOOST_AUTO_TEST_CASE_TEMPLATE(modpow2, mp_int_type, mp_int_types)
-{
- mp_int_type x("4");
- const mp_int_type exp("13");
- const mp_int_type m("497");
- x.modpow(exp, m);
- BOOST_CHECK_EQUAL(x, "445");
-}
-
-BOOST_AUTO_TEST_CASE_TEMPLATE(modpow3, mp_int_type, mp_int_types)
-{
- mp_int_type x("2395422");
- const mp_int_type exp("2424832");
- const mp_int_type m("2424833");
- x.modpow(exp, m);
- BOOST_CHECK_EQUAL(x, "1");
-}
-
-BOOST_AUTO_TEST_CASE_TEMPLATE(modpow4, mp_int_type, mp_int_types)
-{
- mp_int_type x("184");
- const mp_int_type exp("560");
- const mp_int_type m("561");
- x.modpow(exp, m);
- BOOST_CHECK_EQUAL(x, "1");
-}
-
-BOOST_AUTO_TEST_CASE_TEMPLATE(modpow5, mp_int_type, mp_int_types)
-{
- mp_int_type x("997028168093060821869770104094480850560519901475");
- const mp_int_type exp("7455602825647884208337395736200454918783366342656");
- const mp_int_type m("7455602825647884208337395736200454918783366342657");
- x.modpow(exp, m);
- BOOST_CHECK_EQUAL(x, "1");
-}
-
-BOOST_AUTO_TEST_CASE_TEMPLATE(modpow6, mp_int_type, mp_int_types)
-{
- mp_int_type x("184");
- const mp_int_type exp("5600");
- const mp_int_type m("2668");
- x.modpow(exp, m);
- BOOST_CHECK_EQUAL(x, "552");
+ const mp_int_type x("2");
+ const mp_int_type z = pow(x, 64);
+ BOOST_CHECK_EQUAL(z, "18446744073709551616");
}
Modified: sandbox/mp_math/libs/mp_math/test/prerequisite.hpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/prerequisite.hpp (original)
+++ sandbox/mp_math/libs/mp_math/test/prerequisite.hpp 2008-10-27 20:50:18 EDT (Mon, 27 Oct 2008)
@@ -3,6 +3,7 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
+#include <boost/cstdint.hpp>
#include <boost/mp_math/mp_int.hpp>
#include <boost/mpl/vector.hpp>
#include <boost/mpl/unique.hpp>
@@ -25,10 +26,12 @@
std::allocator<void>,
boost::mp_math::mp_int_traits<boost::uint16_t, boost::uint32_t>
>,
+#ifndef BOOST_NO_INT64_T
boost::mp_math::mp_int<
std::allocator<void>,
boost::mp_math::mp_int_traits<boost::uint32_t, boost::uint64_t>
>,
+#endif
boost::mp_math::mp_int<>
> some_mp_int_types;
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk