|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r49479 - in sandbox/mp_math: boost/mp_math boost/mp_math/mp_int boost/mp_math/mp_int/detail libs/mp_math/test
From: baraclese_at_[hidden]
Date: 2008-10-28 16:42:10
Author: baraclese
Date: 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
New Revision: 49479
URL: http://svn.boost.org/trac/boost/changeset/49479
Log:
* General restructuring of code
- Made some member functions non-member functions and put them into namespace
detail.
- Also deleted a few friend declarations in mp_int<>.
- Created mp_int_fwd.hpp to facilitate independence of header files.
- Use include guards in gcd.hpp, jacobi.hpp, lcm.hpp, modinv.hpp,
modpow_ctx.hpp, modpow.hpp, prime.hpp, root.hpp
* A few more comment fixes.
Added:
sandbox/mp_math/boost/mp_math/mp_int/detail/modinv.hpp (contents, props changed)
sandbox/mp_math/boost/mp_math/mp_int/detail/modpow.hpp (contents, props changed)
sandbox/mp_math/boost/mp_math/mp_int/detail/modular_reduction.hpp (contents, props changed)
- copied, changed from r49471, /sandbox/mp_math/boost/mp_math/mp_int/modular_reduction.hpp
sandbox/mp_math/boost/mp_math/mp_int/modpow_ctx.hpp (contents, props changed)
sandbox/mp_math/boost/mp_math/mp_int/mp_int_fwd.hpp (contents, props changed)
Removed:
sandbox/mp_math/boost/mp_math/mp_int/modular_reduction.hpp
Text files modified:
sandbox/mp_math/boost/mp_math/mp_int.hpp | 10
sandbox/mp_math/boost/mp_math/mp_int/abs.hpp | 2
sandbox/mp_math/boost/mp_math/mp_int/detail/modular_reduction.hpp | 14 +
sandbox/mp_math/boost/mp_math/mp_int/gcd.hpp | 16 +
sandbox/mp_math/boost/mp_math/mp_int/jacobi.hpp | 19 +
sandbox/mp_math/boost/mp_math/mp_int/lcm.hpp | 15 +
sandbox/mp_math/boost/mp_math/mp_int/modinv.hpp | 189 +------------
sandbox/mp_math/boost/mp_math/mp_int/modpow.hpp | 521 +--------------------------------------
sandbox/mp_math/boost/mp_math/mp_int/mp_int.hpp | 32 --
sandbox/mp_math/boost/mp_math/mp_int/prime.hpp | 14
sandbox/mp_math/boost/mp_math/mp_int/root.hpp | 75 +++--
sandbox/mp_math/libs/mp_math/test/prerequisite.hpp | 2
12 files changed, 178 insertions(+), 731 deletions(-)
Modified: sandbox/mp_math/boost/mp_math/mp_int.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -6,6 +6,16 @@
#ifndef BOOST_MP_MATH_MP_INT_HPP
#define BOOST_MP_MATH_MP_INT_HPP
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
#include <boost/mp_math/mp_int/mp_int.hpp>
+#include <boost/mp_math/mp_int/gcd.hpp>
+#include <boost/mp_math/mp_int/jacobi.hpp>
+#include <boost/mp_math/mp_int/lcm.hpp>
+#include <boost/mp_math/mp_int/modinv.hpp>
+#include <boost/mp_math/mp_int/modpow.hpp>
+#include <boost/mp_math/mp_int/numeric_limits.hpp>
+#include <boost/mp_math/mp_int/prime.hpp>
+#include <boost/mp_math/mp_int/root.hpp>
+
#endif
Modified: sandbox/mp_math/boost/mp_math/mp_int/abs.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/abs.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/abs.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -7,6 +7,6 @@
mp_int<A,T> abs(const mp_int<A,T>& x)
{
mp_int<A,T> tmp(x);
- tmp.sign_ = 1;
+ tmp.set_sign(1);
return tmp;
}
Added: sandbox/mp_math/boost/mp_math/mp_int/detail/modinv.hpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/boost/mp_math/mp_int/detail/modinv.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -0,0 +1,177 @@
+// 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)
+
+#ifndef BOOST_MP_MATH_MP_INT_DETAIL_MODINV_HPP
+#define BOOST_MP_MATH_MP_INT_DETAIL_MODINV_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+
+
+namespace boost {
+namespace mp_math {
+namespace detail {
+
+// hac 14.61, pp608
+template<class A1, class T>
+mp_int<A1,T> even_modinv(const mp_int<A1,T>& z, const mp_int<A1,T>& m)
+{
+ typedef typename mp_int<A1,T>::digit_type digit_type;
+
+ assert(m.is_even());
+
+ static const char* const err_msg = "mp_int::modinv: inverse does not exist";
+
+ const mp_int<A1,T> x = z % m;
+
+ if (x.is_even())
+ throw std::domain_error(err_msg);
+
+ mp_int<A1,T> u(x);
+ mp_int<A1,T> v(m);
+ mp_int<A1,T> A = digit_type(1);
+ mp_int<A1,T> B = digit_type(0);
+ mp_int<A1,T> C = digit_type(0);
+ mp_int<A1,T> D = digit_type(1);
+
+top:
+ while (u.is_even())
+ {
+ u.divide_by_2();
+
+ if (A.is_odd() || B.is_odd())
+ {
+ // A = (A+m)/2, B = (B-x)/2
+ A += m;
+ B -= x;
+ }
+ A.divide_by_2();
+ B.divide_by_2();
+ }
+
+ while (v.is_even())
+ {
+ v.divide_by_2();
+
+ if (C.is_odd() || D.is_odd())
+ {
+ // C = (C+m)/2, D = (D-x)/2
+ C += m;
+ D -= x;
+ }
+ C.divide_by_2();
+ D.divide_by_2();
+ }
+
+ if (u >= v)
+ {
+ u -= v;
+ A -= C;
+ B -= D;
+ }
+ else
+ {
+ v -= u;
+ C -= A;
+ D -= B;
+ }
+
+ if (u)
+ goto top;
+
+ // now a = C, b = D, gcd == g*v
+
+ // if v != 1 then there is no inverse
+ if (v != digit_type(1))
+ throw std::domain_error(err_msg);
+
+ // if it's too low
+ while (C.compare_to_digit(0) == -1)
+ C += m;
+
+ // too big
+ while (C.compare_magnitude(m) != -1)
+ C -= m;
+
+ return C;
+}
+
+/* computes the modular inverse via binary extended euclidean algorithm,
+ * that is z = 1 / z mod m
+ *
+ * Based on even modinv except this is optimized for the case where m is
+ * odd as per HAC Note 14.64 on pp. 610
+ */
+template<class A1, class T>
+mp_int<A1,T> odd_modinv(const mp_int<A1,T>& z, const mp_int<A1,T>& m)
+{
+ typedef typename mp_int<A1,T>::digit_type digit_type;
+
+ assert(m.is_odd());
+
+ // m == modulus, y == value to invert
+ // we need y = |a|
+ const mp_int<A1,T> y = z % m;
+
+ mp_int<A1,T> u(m);
+ mp_int<A1,T> v(y);
+ mp_int<A1,T> A = digit_type(1);
+ mp_int<A1,T> B = digit_type(0);
+ mp_int<A1,T> C = digit_type(0);
+ mp_int<A1,T> D = digit_type(1);
+
+top:
+ while (u.is_even())
+ {
+ u.divide_by_2();
+
+ if (B.is_odd())
+ B -= m;
+
+ B.divide_by_2();
+ }
+
+ while (v.is_even())
+ {
+ v.divide_by_2();
+
+ if (D.is_odd())
+ D -= m;
+
+ D.divide_by_2();
+ }
+
+ if (u >= v)
+ {
+ /* u = u - v, B = B - D */
+ u -= v;
+ B -=D;
+ }
+ else
+ {
+ v -= u;
+ D -= B;
+ }
+
+ if (u)
+ goto top;
+
+ /* now a = C, m = D, gcd == g*v */
+
+ if (v != digit_type(1))
+ throw std::domain_error("mp_int::modinv: inverse does not exist");
+
+ while (D.is_negative())
+ D += m;
+
+ return D;
+}
+
+
+} // namespace detail
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
Added: sandbox/mp_math/boost/mp_math/mp_int/detail/modpow.hpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/boost/mp_math/mp_int/detail/modpow.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -0,0 +1,349 @@
+// 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)
+
+
+#ifndef BOOST_MP_MATH_MP_INT_DETAIL_MODPOW_HPP
+#define BOOST_MP_MATH_MP_INT_DETAIL_MODPOW_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+#include <boost/mp_math/mp_int/detail/modular_reduction.hpp>
+
+namespace boost {
+namespace mp_math {
+namespace detail {
+
+
+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;
+}
+
+
+} // namespace detail
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
Copied: sandbox/mp_math/boost/mp_math/mp_int/detail/modular_reduction.hpp (from r49471, /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/detail/modular_reduction.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -3,6 +3,14 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MP_MATH_MP_INT_DETAIL_MODULAR_REDUCTION_HPP
+#define BOOST_MP_MATH_MP_INT_DETAIL_MODULAR_REDUCTION_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+
+namespace boost {
+namespace mp_math {
+namespace detail {
// reduces x mod m, assumes 0 < x < m**2, mu is precomputed.
// From HAC pp.604 Algorithm 14.42
@@ -346,3 +354,9 @@
}
+} // namespace detail
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
Modified: sandbox/mp_math/boost/mp_math/mp_int/gcd.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/gcd.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/gcd.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -3,6 +3,15 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MP_MATH_MP_INT_GCD_HPP
+#define BOOST_MP_MATH_MP_INT_GCD_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+
+
+namespace boost {
+namespace mp_math {
+
// Greatest Common Divisor using the binary method
template<class A, class T>
mp_int<A,T> gcd(const mp_int<A,T>& a, const mp_int<A,T>& b)
@@ -52,3 +61,10 @@
return gcd(gcd(a, b), args...);
}
#endif
+
+
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
Modified: sandbox/mp_math/boost/mp_math/mp_int/jacobi.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/jacobi.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/jacobi.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -3,6 +3,15 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MP_MATH_MP_INT_JACOBI_HPP
+#define BOOST_MP_MATH_MP_INT_JACOBI_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+
+
+namespace boost {
+namespace mp_math {
+
// computes the jacobi c = (a | p) (or Legendre if p is prime)
// HAC pp. 73 Algorithm 2.149
template<class A, class T>
@@ -37,7 +46,7 @@
else
{
// calculate p.digits_[0] mod 8
- const digit_type residue = p.digits_[0] & 7;
+ const digit_type residue = p[0] & 7;
if (residue == 1 || residue == 7)
s = 1;
@@ -46,7 +55,7 @@
}
/* if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
- if (((p.digits_[0] & 3) == 3) && ((a1.digits_[0] & 3) == 3))
+ if (((p[0] & 3) == 3) && ((a1[0] & 3) == 3))
s = -s;
if (a1 == digit_type(1))
@@ -58,3 +67,9 @@
}
}
+
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
Modified: sandbox/mp_math/boost/mp_math/mp_int/lcm.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/lcm.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/lcm.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -3,6 +3,15 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MP_MATH_MP_INT_LCM_HPP
+#define BOOST_MP_MATH_MP_INT_LCM_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+
+
+namespace boost {
+namespace mp_math {
+
// computes least common multiple as |a*b|/gcd(a,b)
template<class A, class T>
mp_int<A,T> lcm(const mp_int<A,T>& a, const mp_int<A,T>& b)
@@ -30,3 +39,9 @@
}
#endif
+
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
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-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -3,176 +3,18 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MP_MATH_MP_INT_MODINV_HPP
+#define BOOST_MP_MATH_MP_INT_MODINV_HPP
-// hac 14.61, pp608
-// 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& m)
-{
- 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 (m.is_odd())
- odd_modinv(m);
- else
- even_modinv(m);
-}
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+#include <boost/mp_math/mp_int/detail/modinv.hpp>
-// hac 14.61, pp608
-template<class A1, class T>
-void mp_int<A1,T>::even_modinv(const mp_int& y)
-{
- assert(y.is_even());
- static const char* const err_msg = "mp_int::modinv: inverse does not exist";
-
- const mp_int x = *this % y;
-
- if (x.is_even())
- throw std::domain_error(err_msg);
-
- mp_int u(x);
- mp_int v(y);
- mp_int A = digit_type(1);
- mp_int B = digit_type(0);
- mp_int C = digit_type(0);
- mp_int D = digit_type(1);
-
-top:
- while (u.is_even())
- {
- u.divide_by_2();
-
- if (A.is_odd() || B.is_odd())
- {
- // A = (A+y)/2, B = (B-x)/2
- A += y;
- B -= x;
- }
- A.divide_by_2();
- B.divide_by_2();
- }
-
- while (v.is_even())
- {
- v.divide_by_2();
-
- if (C.is_odd() || D.is_odd())
- {
- // C = (C+y)/2, D = (D-x)/2
- C += y;
- D -= x;
- }
- C.divide_by_2();
- D.divide_by_2();
- }
-
- if (u >= v)
- {
- u -= v;
- A -= C;
- B -= D;
- }
- else
- {
- v -= u;
- C -= A;
- D -= B;
- }
-
- if (u)
- goto top;
-
- // now a = C, b = D, gcd == g*v
-
- // if v != 1 then there is no inverse
- if (v != digit_type(1))
- throw std::domain_error(err_msg);
-
- // if it's too low
- while (C.compare_to_digit(0) == -1)
- C += y;
-
- // too big
- while (C.compare_magnitude(y) != -1)
- C -= y;
-
- *this = C;
-}
-
-/* computes the modular inverse via binary extended euclidean algorithm,
- * that is *this = 1 / *this mod x
- *
- * Based on even modinv except this is optimized for the case where x is
- * odd as per HAC Note 14.64 on pp. 610
- */
-template<class A1, class T>
-void mp_int<A1,T>::odd_modinv(const mp_int& x)
-{
- assert(x.is_odd());
-
- // x == modulus, y == value to invert
- // we need y = |a|
- const mp_int y = *this % x;
-
- // 3. u=x, v=y, A=1, B=0, C=0, D=1
- mp_int u(x);
- mp_int v(y);
- mp_int A = digit_type(1);
- mp_int B = digit_type(0);
- mp_int C = digit_type(0);
- mp_int D = digit_type(1);
-
-top:
- while (u.is_even())
- {
- u.divide_by_2();
-
- if (B.is_odd())
- B -= x;
-
- B.divide_by_2();
- }
-
- while (v.is_even())
- {
- v.divide_by_2();
-
- if (D.is_odd())
- D -= x;
-
- D.divide_by_2();
- }
-
- if (u >= v)
- {
- /* u = u - v, B = B - D */
- u -= v;
- B -=D;
- }
- else
- {
- v -= u;
- D -= B;
- }
-
- if (u)
- goto top;
-
- /* now a = C, x = D, gcd == g*v */
-
- if (v != digit_type(1))
- throw std::domain_error("mp_int::modinv: inverse does not exist");
-
- while (D.is_negative())
- D += x;
-
- *this = D;
-}
+namespace boost {
+namespace mp_math {
+// hac 14.61, pp608
// returns the modular multiplicative inverse x of a (mod m) such that
// x*a = 1 (mod m) =>
// a^-1 = x (mod m)
@@ -181,8 +23,19 @@
template<class A, class T>
mp_int<A,T> modinv(const mp_int<A,T>& a, const mp_int<A,T>& m)
{
- mp_int<A,T> x(a);
- x.modinv(m);
- return x;
+ 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 (m.is_odd())
+ return detail::odd_modinv(a, m);
+ else
+ return detail::even_modinv(a, m);
}
+
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
Modified: sandbox/mp_math/boost/mp_math/mp_int/modpow.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/modpow.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/modpow.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -3,499 +3,16 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MP_MATH_MP_INT_MODPOW_HPP
+#define BOOST_MP_MATH_MP_INT_MODPOW_HPP
-// 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;
-}
+#include <boost/mp_math/mp_int/modpow_ctx.hpp>
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+#include <boost/mp_math/mp_int/detail/modpow.hpp>
+#include <boost/mp_math/mp_int/detail/modular_reduction.hpp>
+namespace boost {
+namespace mp_math {
// z = base^exp % mod
template<class A, class T>
@@ -531,11 +48,23 @@
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_restricted_dr:
+ return detail::montgomery_modpow(base, exp, mod, 1, ctx->rho);
+ case ctx_t::mod_unrestricted_dr:
+ return detail::montgomery_modpow(base, exp, mod, 2, ctx->rho);
+ case ctx_t::mod_unrestricted_dr_slow:
+ return detail::barret_modpow (base, exp, mod, 1, ctx->mu);
+ case ctx_t::mod_odd:
+ return detail::montgomery_modpow(base, exp, mod, 0, ctx->rho);
case ctx_t::mod_generic:
- default: return barret_modpow (base, exp, mod, 0, ctx->mu);
+ default:
+ return detail::barret_modpow (base, exp, mod, 0, ctx->mu);
}
}
+
+
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
Added: sandbox/mp_math/boost/mp_math/mp_int/modpow_ctx.hpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/boost/mp_math/mp_int/modpow_ctx.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -0,0 +1,186 @@
+// 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)
+
+#ifndef BOOST_MP_MATH_MP_INT_MODPOW_CTX_HPP
+#define BOOST_MP_MATH_MP_INT_MODPOW_CTX_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+
+namespace boost {
+namespace mp_math {
+
+
+// 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;
+}
+
+
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
Deleted: sandbox/mp_math/boost/mp_math/mp_int/modular_reduction.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/modular_reduction.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
+++ (empty file)
@@ -1,348 +0,0 @@
-// 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)
-
-
-// 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 barrett_reduce(mp_int<A,T>& x, const mp_int<A,T>& m, const mp_int<A,T>& mu)
-{
- typedef typename mp_int<A,T>::digit_type digit_type;
-
- const typename mp_int<A,T>::size_type k = m.size();
-
- 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) << (mp_int<A,T>::valid_bits - 1))
- q *= mu;
- else
- q.fast_mul_high_digits(mu, k);
-
- // q = q / base**(k+1)
- q.shift_digits_right(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 base**(k+1)
- q.mul_digits(m, k + 1);
-
- x -= q;
-
- // If x < 0, add base**(k+1) to it
- if (x.is_negative())
- {
- q = digit_type(1);
- q.shift_digits_left(k + 1);
- x += q;
- }
-
- while (x >= m)
- x.sub_smaller_magnitude(m);
-}
-
-/* computes xR**-1 == x (mod N) via Montgomery Reduction */
-template<class A, class T>
-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;
-
- // x = x + mu * m * base**i
-
- // alias for digits of the modulus
- const digit_type* tmpn = n.digits();
-
- // 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.size(); ++j)
- {
- // 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 = static_cast<digit_type>(r >> mp_int<A,T>::digit_bits);
-
- *tmpx++ = r;
- }
- // At this point the i'th digit of x should be zero
-
- // propagate carries upwards as required
- while (carry)
- {
- const word_type r = static_cast<word_type>(*tmpx) + carry;
- *tmpx++ = r;
- 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/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 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 montgomery_normalize(mp_int<A,T>& x, const mp_int<A,T>& n)
-{
- // 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 (n.size() > 1)
- x.pow2((n.size() - 1) * mp_int<A,T>::valid_bits + bits - 1);
- else
- {
- x = typename mp_int<A,T>::digit_type(1);
- bits = 1;
- }
-
- // now compute C = A * B mod n
- for (int i = bits - 1; i < mp_int<A,T>::valid_bits; ++i)
- {
- x.multiply_by_2();
- if (x.compare_magnitude(n) != -1)
- x.sub_smaller_magnitude(n);
- }
-}
-
-/* computes xR**-1 == x (mod N) via Montgomery Reduction
- *
- * This is an optimized implementation of montgomery_reduce
- * which uses the comba method to quickly calculate the columns of the
- * reduction.
- *
- * Based on Algorithm 14.32 on pp.601 of HAC.
-*/
-template<class A, class T>
-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
- *
- * This is computed in place and on the fly. The multiplication
- * by b**i is handled by offseting which columns the results
- * are added to.
- *
- * Note the comba method normally doesn't handle carries in the
- * inner loop In this case we fix the carry from the previous
- * column since the Montgomery reduction requires digits of the
- * result (so far) [see above] to work. This is
- * handled by fixing up one carry after the inner loop. The
- * carry fixups are done in order so after these loops the
- * first m->used words of W[] have the carries fixed
- */
- /* inner loop */
- 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>(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.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);
-}
-
-
-// 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));
-
-top:
- // set carry to zero
- digit_type mu = 0;
-
- // 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>(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
- x[m] = mu;
-
- // 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));
-
- x.clamp();
- if (x.is_zero())
- x.set_sign(1);
-
- if (x.compare_magnitude(n) != -1)
- {
- x.sub_smaller_magnitude(n);
- goto top;
- }
-}
-
-// reduces x modulo n where n is of the form 2**p - d
-template<class A, class T>
-void unrestricted_dr_reduce(mp_int<A,T>& x,
- const mp_int<A,T>& n,
- typename mp_int<A,T>::digit_type d)
-{
- const typename mp_int<A,T>::size_type p = n.precision();
-
-top:
-
- mp_int<A,T> q(x);
-
- /* q = a/2**p, r = r mod 2**p */
- q.shift_right(p, &x);
-
- if (d != 1)
- q.multiply_by_digit(d);
-
- x.add_magnitude(q);
-
- if (x.compare_magnitude(n) != -1)
- {
- x.sub_smaller_magnitude(n);
- goto top;
- }
-}
-
-// 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>
-void unrestricted_dr_slow_reduce(mp_int<A,T>& x,
- const mp_int<A,T>& n,
- const mp_int<A,T>& d)
-{
- const typename mp_int<A,T>::size_type p = n.precision();
-
-top:
-
- mp_int<A,T> q(x);
-
- // q = x/2**p, r = r mod 2**p
- q.shift_right(p, &x);
-
- q *= d;
-
- x.add_magnitude(q);
-
- if (x.compare_magnitude(n) != -1)
- {
- x.sub_smaller_magnitude(n);
- goto top;
- }
-}
-
-
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-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -22,24 +22,20 @@
#include <boost/serialization/string.hpp>
#include <boost/type_traits/is_integral.hpp>
#include <boost/utility/enable_if.hpp>
-#include <boost/mp_math/mp_int/traits.hpp>
#include <boost/mp_math/mp_int/detail/string_conversion_constants.hpp>
#include <boost/mp_math/mp_int/detail/integral_ops.hpp>
#include <boost/mp_math/mp_int/detail/meta_math.hpp>
-#include <boost/mp_math/mp_int/detail/prime_tab.hpp>
#include <boost/mp_math/mp_int/detail/primitive_ops.hpp>
namespace boost {
namespace mp_math {
-template<class A, class T> struct modpow_ctx;
-
// digits are stored in least significant order
template<
- class Allocator = std::allocator<void>,
- class Traits = mp_int_traits<>
+ class Allocator,
+ class Traits
>
struct mp_int
:
@@ -325,12 +321,6 @@
void pow2(size_type b);
- 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);
- void odd_modinv(const mp_int& modulus);
-
void set_least_significant_bit()
{
digits_[0] |= digit_type(1);
@@ -348,13 +338,6 @@
template<class A, class T>
friend bool operator < (const mp_int<A,T>&, const mp_int<A,T>&);
- template<class A, class T>
- friend mp_int<A,T> abs(const mp_int<A,T>& x);
- template<class A, class T>
- friend mp_int<A,T> gcd(const mp_int<A,T>& a, const mp_int<A,T>& b);
- template<class A, class T>
- friend int jacobi(const mp_int<A,T>& a, const mp_int<A,T>& b);
-
template<typename Iter>
void from_string(Iter first, Iter last, unsigned radix);
@@ -852,19 +835,11 @@
#include <boost/mp_math/mp_int/add.hpp>
#include <boost/mp_math/mp_int/ctors.hpp>
#include <boost/mp_math/mp_int/div.hpp>
-#include <boost/mp_math/mp_int/gcd.hpp>
-#include <boost/mp_math/mp_int/jacobi.hpp>
-#include <boost/mp_math/mp_int/lcm.hpp>
#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>
#include <boost/mp_math/mp_int/random.hpp>
-#include <boost/mp_math/mp_int/prime.hpp>
-#include <boost/mp_math/mp_int/root.hpp>
#include <boost/mp_math/mp_int/sqr.hpp>
#include <boost/mp_math/mp_int/sub.hpp>
#include <boost/mp_math/mp_int/string_conversion.hpp>
@@ -872,7 +847,4 @@
} // namespace mp_math
} // namespace boost
-#include <boost/mp_math/mp_int/numeric_limits.hpp>
-
-
#endif
Added: sandbox/mp_math/boost/mp_math/mp_int/mp_int_fwd.hpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/boost/mp_math/mp_int/mp_int_fwd.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -0,0 +1,26 @@
+// 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)
+
+#ifndef BOOST_MP_MATH_MP_INT_MP_INT_FWD_HPP
+#define BOOST_MP_MATH_MP_INT_MP_INT_FWD_HPP
+
+#include <memory>
+#include <boost/mp_math/mp_int/traits.hpp>
+
+
+namespace boost {
+namespace mp_math {
+
+template<
+ class Allocator = std::allocator<void>,
+ class Traits = mp_int_traits<>
+>
+struct mp_int;
+
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
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-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -3,6 +3,16 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MP_MATH_MP_INT_PRIME_HPP
+#define BOOST_MP_MATH_MP_INT_PRIME_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+#include <boost/mp_math/mp_int/detail/prime_tab.hpp>
+
+
+namespace boost {
+namespace mp_math {
+
// this one simply divides x by some small primes found in prime_tab
// returns false if definitely composite, returns true if probably prime
// this test is to be used to quickly filter out composite numbers
@@ -333,6 +343,8 @@
}
+} // namespace mp_math
+} // namespace boost
-
+#endif
Modified: sandbox/mp_math/boost/mp_math/mp_int/root.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/root.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/root.hpp 2008-10-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -3,16 +3,23 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MP_MATH_MP_INT_ROOT_HPP
+#define BOOST_MP_MATH_MP_INT_ROOT_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+
+
+namespace boost {
+namespace mp_math {
+
template<class A, class T>
mp_int<A,T> sqrt(const mp_int<A,T>& x)
{
- /* must be positive */
- if (x.sign() == -1)
- throw std::domain_error("sqrt: argument is negative");
+ if (x.is_negative())
+ throw std::domain_error("sqrt: argument must be positive");
mp_int<A,T> t1;
- /* easy out */
if (x.is_zero())
{
t1.zero();
@@ -21,80 +28,78 @@
t1 = x;
- /* First approx. (not very bad for large arg) */
+ // First approx. (not very bad for large arg)
t1.shift_digits_right(t1.size()/2);
- /* t1 > 0 */
+ // t1 > 0
mp_int<A,T> t2 = x / t1;
t1 += t2;
t1.divide_by_2();
- /* And now t1 > sqrt(arg) */
+ // And now t1 > sqrt(arg)
do
{
t2 = x / t1;
t1 += t2;
t1.divide_by_2();
- /* t1 >= sqrt(arg) >= t2 at this point */
+ // t1 >= sqrt(arg) >= t2 at this point
} while (t1.compare_magnitude(t2) == 1);
return t1;
}
-/* find the n'th root of an integer
- *
- * Result found such that (c)**b <= a and (c+1)**b > a
- *
- * This algorithm uses Newton's approximation
- * x[i+1] = x[i] - f(x[i])/f'(x[i])
- * which will find the root in log(N) time where
- * each step involves a fair bit. This is not meant to
- * find huge roots [square and cube, etc].
- */
+// find the n'th root of an integer
+//
+// Result found such that (c)**b <= a and (c+1)**b > a
+//
+// This algorithm uses Newton's approximation
+// x[i+1] = x[i] - f(x[i])/f'(x[i])
+// which will find the root in log(N) time where
+// each step involves a fair bit. This is not meant to
+// find huge roots [square and cube, etc].
template<class A, class T>
mp_int<A,T> nth_root(const mp_int<A,T>& x, typename mp_int<A,T>::digit_type n)
{
mp_int<A,T> t1, t2, t3;
- /* input must be positive if n is even */
if ((n & 1) == 0 && x.is_negative())
throw std::domain_error("nth_root: argument must be positive if n is even");
- /* if x is negative fudge the sign but keep track */
+ // if x is negative fudge the sign but keep track
const int neg = x.sign();
const_cast<mp_int<A,T>*>(&x)->set_sign(1);
- /* t2 = 2 */
+ // t2 = 2
t2 = typename mp_int<A,T>::digit_type(2);
do
{
t1 = t2;
- /* t2 = t1 - ((t1**n - x) / (n * t1**(n-1))) */
+ // t2 = t1 - ((t1**n - x) / (n * t1**(n-1)))
- /* t3 = t1**(n-1) */
+ // t3 = t1**(n-1)
t3 = pow(t1, n-1);
- /* numerator */
- /* t2 = t1**n */
+ // numerator
+ // t2 = t1**n
t2 = t3 * t1;
- /* t2 = t1**n - x */
+ // t2 = t1**n - x
t2 -= x;
- /* denominator */
- /* t3 = t1**(n-1) * n */
+ // denominator
+ // t3 = t1**(n-1) * n
t3.multiply_by_digit(n);
- /* t3 = (t1**n - x)/(n * t1**(n-1)) */
+ // t3 = (t1**n - x)/(n * t1**(n-1))
t3 = t2 / t3;
t2 = t1 - t3;
} while (t1 != t2);
- /* result can be off by a few so check */
+ // result can be off by a few so check
for (;;)
{
t2 = pow(t1, n);
@@ -105,12 +110,18 @@
break;
}
- /* reset the sign of x first */
+ // reset the sign of x first
const_cast<mp_int<A,T>*>(&x)->set_sign(neg);
- /* set the sign of the result */
+ // set the sign of the result
t1.set_sign(neg);
return t1;
}
+
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
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-28 16:42:08 EDT (Tue, 28 Oct 2008)
@@ -5,8 +5,8 @@
#include <boost/cstdint.hpp>
#include <boost/mp_math/mp_int.hpp>
-#include <boost/mpl/vector.hpp>
#include <boost/mpl/unique.hpp>
+#include <boost/mpl/vector.hpp>
#include <boost/type_traits/is_same.hpp>
//typedef boost::mp_math::mp_int_traits<boost::uint8_t, boost::uint16_t> traits_type;
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