Subject: [Boost-commit] svn:boost r50541 - in sandbox/mp_math: . boost/mp_math/mp_int boost/mp_math/mp_int/detail libs/mp_math/test libs/mp_math/tools/benchmark
From: baraclese_at_[hidden]
Date: 2009-01-10 19:30:29
Author: baraclese
Date: 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
New Revision: 50541
URL: http://svn.boost.org/trac/boost/changeset/50541
Log:
* mp_int/ctors.hpp
Ctors now have an allocator reference parameter.
Fixed a bug in dtor reported by Nathan Kitchen on Boost devel mailing list on 10th January.
* mp_int/sqr.hpp
Added small optimization to comba_sqr().
* mp_int/detail/primitive_ops.hpp
Added multiply_by_two, divide_by_two and divide_by_digit functions.
* mp_int/detail/integral_ops.hpp
Added parenthesis around expression to calculate q constant, reported by Nathan Kitchen.
* mp_int/div.hpp
Moved division algo to detail/div.hpp and renamed function to classic_divide.
* mp_int/detail/div.hpp
New.
* mp_int/root.hpp
Fixed nth_root. Added better initial approximation for Newton method.
* mp_int/mp_int.hpp
Added clamp_high_digit().
Added clear_bit, set_bits, clear_bits, truncate functions.
Removed is_zero().
Removed divide().
Fixed deallocate() null ptr bug as reported by Nathan Kitchen.
* test/bitmanipulation.cpp
New file to test set_bit, clear_bit functions.
Added:
sandbox/mp_math/boost/mp_math/mp_int/detail/div.hpp (contents, props changed)
sandbox/mp_math/libs/mp_math/test/bitmanipulation.cpp (contents, props changed)
Text files modified:
sandbox/mp_math/boost/mp_math/mp_int/add.hpp | 7
sandbox/mp_math/boost/mp_math/mp_int/ctors.hpp | 46 +++++-
sandbox/mp_math/boost/mp_math/mp_int/detail/integral_ops.hpp | 45 +++---
sandbox/mp_math/boost/mp_math/mp_int/detail/modular_reduction.hpp | 6
sandbox/mp_math/boost/mp_math/mp_int/detail/primitive_ops.hpp | 262 +++++++++++++++++++++++++++++----------
sandbox/mp_math/boost/mp_math/mp_int/div.hpp | 223 ++++-----------------------------
sandbox/mp_math/boost/mp_math/mp_int/jacobi.hpp | 2
sandbox/mp_math/boost/mp_math/mp_int/mod.hpp | 3
sandbox/mp_math/boost/mp_math/mp_int/mp_int.hpp | 173 ++++++++++++++++++++------
sandbox/mp_math/boost/mp_math/mp_int/mul.hpp | 33 +---
sandbox/mp_math/boost/mp_math/mp_int/operators.hpp | 24 +-
sandbox/mp_math/boost/mp_math/mp_int/root.hpp | 53 +++++--
sandbox/mp_math/boost/mp_math/mp_int/sqr.hpp | 34 ++++-
sandbox/mp_math/boost/mp_math/mp_int/string_conversion.hpp | 39 +++--
sandbox/mp_math/libs/mp_math/test/cmp.cpp | 6
sandbox/mp_math/libs/mp_math/test/ctors.cpp | 4
sandbox/mp_math/libs/mp_math/test/div.cpp | 16 ++
sandbox/mp_math/libs/mp_math/test/integral_ops.cpp | 17 ++
sandbox/mp_math/libs/mp_math/test/jamfile.v2 | 19 +-
sandbox/mp_math/libs/mp_math/test/root.cpp | 14 ++
sandbox/mp_math/libs/mp_math/test/shift.cpp | 7 +
sandbox/mp_math/libs/mp_math/test/string_ops.cpp | 2
sandbox/mp_math/libs/mp_math/tools/benchmark/benchmark_gmp.hpp | 2
sandbox/mp_math/project-root.jam | 2
24 files changed, 603 insertions(+), 436 deletions(-)
Modified: sandbox/mp_math/boost/mp_math/mp_int/add.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/add.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/add.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -19,7 +19,7 @@
if (digits_[0] >= b) // example: -16 + 5 = -11; or -5 + 5 = 0
{
digits_[0] -= b;
- if (is_zero())
+ if (!*this)
set_sign(1);
}
else
@@ -77,10 +77,7 @@
return;
}
else if (carry) // at this point n equals x->size_
- {
- digits_[n] = carry;
- ++n;
- }
+ digits_[n++] = carry;
size_ = n;
}
Modified: sandbox/mp_math/boost/mp_math/mp_int/ctors.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/ctors.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/ctors.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -1,4 +1,4 @@
-// Copyright Kevin Sopp 2008.
+// Copyright Kevin Sopp 2008 - 2009.
// 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)
@@ -13,11 +13,23 @@
}
template<class A, class T>
+mp_int<A,T>::mp_int(const allocator_type& a)
+:
+ base_allocator_type(a),
+ digits_(0),
+ size_(0),
+ capacity_(0)
+{
+}
+
+template<class A, class T>
template<typename IntegralT>
mp_int<A,T>::mp_int(
IntegralT b,
+ const allocator_type& a,
typename enable_if<is_integral<IntegralT> >::type*)
:
+ base_allocator_type(a),
digits_(0),
size_(0),
capacity_(0)
@@ -50,10 +62,10 @@
++c;
sign = 1;
}
-
+
// detect the radix
unsigned int radix;
-
+
if (c != last)
{
if (*c == '0') // octal
@@ -111,7 +123,7 @@
}
set_sign(1);
}
-
+
const bool uppercase = f & std::ios_base::uppercase;
const bool showbase = f & std::ios_base::showbase;
@@ -158,8 +170,11 @@
template<class A, class T>
template<typename RandomAccessIterator>
-mp_int<A,T>::mp_int(RandomAccessIterator first, RandomAccessIterator last)
+mp_int<A,T>::mp_int(RandomAccessIterator first,
+ RandomAccessIterator last,
+ const allocator_type& a)
:
+ base_allocator_type(a),
digits_(0),
size_(0),
capacity_(0)
@@ -169,8 +184,9 @@
template<class A, class T>
template<typename charT>
-mp_int<A,T>::mp_int(const charT* s)
+mp_int<A,T>::mp_int(const charT* s, const allocator_type& a)
:
+ base_allocator_type(a),
digits_(0),
size_(0),
capacity_(0)
@@ -180,8 +196,11 @@
template<class A, class T>
template<typename charT>
-mp_int<A,T>::mp_int(const charT* s, std::ios_base::fmtflags f)
+mp_int<A,T>::mp_int(const charT* s,
+ std::ios_base::fmtflags f,
+ const allocator_type& a)
:
+ base_allocator_type(a),
digits_(0),
size_(0),
capacity_(0)
@@ -191,8 +210,10 @@
template<class A, class T>
template<typename charT, class traits, class Alloc>
-mp_int<A,T>::mp_int(const std::basic_string<charT,traits,Alloc>& s)
+mp_int<A,T>::mp_int(const std::basic_string<charT,traits,Alloc>& s,
+ const allocator_type& a)
:
+ base_allocator_type(a),
digits_(0),
size_(0),
capacity_(0)
@@ -203,8 +224,10 @@
template<class A, class T>
template<typename charT, class traits, class Alloc>
mp_int<A,T>::mp_int(const std::basic_string<charT,traits,Alloc>& s,
- std::ios_base::fmtflags f)
+ std::ios_base::fmtflags f,
+ const allocator_type& a)
:
+ base_allocator_type(a),
digits_(0),
size_(0),
capacity_(0)
@@ -215,6 +238,8 @@
template<class A, class T>
mp_int<A,T>::mp_int(const mp_int& copy)
+:
+ base_allocator_type(copy.get_allocator())
{
digits_ = this->allocate(copy.size_);
std::memcpy(digits_, copy.digits_, copy.size_ * sizeof(digit_type));
@@ -241,6 +266,7 @@
template<class A, class T>
mp_int<A,T>::~mp_int()
{
- this->deallocate(digits_, capacity());
+ if (digits_)
+ this->deallocate(digits_, capacity());
}
Added: sandbox/mp_math/boost/mp_math/mp_int/detail/div.hpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/boost/mp_math/mp_int/detail/div.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -0,0 +1,174 @@
+// 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_DIV_HPP
+#define BOOST_MP_MATH_MP_INT_DETAIL_DIV_HPP
+
+#include <boost/mp_math/mp_int/mp_int_fwd.hpp>
+
+namespace boost {
+namespace mp_math {
+namespace detail {
+
+// integer signed division.
+// q*b + r == a
+// HAC pp.598 Algorithm 14.20
+//
+// Note that the description in HAC is horribly incomplete. For example, it
+// doesn't consider the case where digits are removed from 'x' in the inner
+// loop. It also doesn't consider the case that y has fewer than three digits,
+// etc..
+// The overall algorithm is as described as 14.20 from HAC but fixed to treat
+// these cases.
+
+// divide a by b, optionally store remainder
+template<class A, class T>
+void classic_divide(const mp_int<A,T>& a, const mp_int<A,T>& b,
+ mp_int<A,T>& q, mp_int<A,T>* remainder = 0)
+{
+ 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;
+
+ if (!b)
+ throw std::domain_error("mp_int::divide: division by zero");
+
+ // if *this < b then q=0, r = *this
+ if (a.compare_magnitude(b) == -1)
+ {
+ if (remainder)
+ *remainder = a;
+ q.zero();
+ return;
+ }
+
+ q.grow_capacity(a.size() + 2);
+ q.set_size(a.size() + 2);
+ std::memset(q.digits(), 0, q.size() * sizeof(digit_type));
+
+ mp_int<A,T> x(a);
+ mp_int<A,T> y(b);
+
+ // fix the sign
+ const int neg = (a.sign() == b.sign()) ? 1 : -1;
+ x.set_sign(1);
+ y.set_sign(1);
+
+ // normalize both x and y, ensure that y >= beta/2, [beta == 2**valid_bits]
+ size_type norm = y.precision() % mp_int<A,T>::valid_bits;
+ if (norm < mp_int<A,T>::valid_bits - 1)
+ {
+ norm = mp_int<A,T>::valid_bits - 1 - norm;
+ x <<= norm;
+ y <<= norm;
+ }
+ else
+ norm = 0;
+
+ // note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4
+ const size_type n = x.size() - 1;
+ const size_type t = y.size() - 1;
+
+ // find leading digit of the quotient
+ // while (x >= y*beta**(n-t)) do { q[n-t] += 1; x -= y*beta**(n-t) }
+ y.shift_digits_left(n - t); // y = y*beta**(n-t)
+
+ while (x.compare(y) != -1)
+ {
+ ++q[n - t];
+ x -= y;
+ }
+
+ // reset y by shifting it back down
+ y.shift_digits_right(n - t);
+
+ // find the remainder of the digits
+ // step 3. for i from n down to (t + 1)
+ for (size_type i = n; i >= (t + 1); --i)
+ {
+ if (i > x.size())
+ continue;
+
+ // step 3.1 if xi == yt then set q{i-t-1} to beta-1,
+ // otherwise set q{i-t-1} to (xi*beta + x{i-1})/yt
+ if (x[i] == y[t])
+ q[i - t - 1] = mp_int<A,T>::digit_max;
+ else
+ {
+ word_type tmp = static_cast<word_type>(x[i])
+ << static_cast<word_type>(mp_int<A,T>::valid_bits);
+ tmp |= x[i - 1];
+ tmp /= y[t];
+ q[i - t - 1] = static_cast<digit_type>(tmp);
+ }
+
+ // now fixup quotient estimation
+ // while (q{i-t-1} * (yt * beta + y{t-1})) >
+ // xi * beta**2 + xi-1 * beta + xi-2
+ //
+ // do q{i-t-1} -= 1;
+
+ mp_int<A,T> t1, t2;
+ t1.grow_capacity(3);
+ t2.grow_capacity(3);
+
+ ++q[i - t - 1];
+ do
+ {
+ --q[i - t - 1];
+
+ // find left hand
+ t1.zero();
+ t1[0] = (t == 0) ? 0 : y[t - 1];
+ t1[1] = y[t];
+ t1.set_size(2);
+ t1.multiply_by_digit(q[i - t - 1]);
+
+ // find right hand
+ t2[0] = (i < 2) ? 0 : x[i - 2];
+ t2[1] = (i == 0) ? 0 : x[i - 1];
+ t2[2] = x[i];
+ t2.set_size(3);
+ } while (t1.compare_magnitude(t2) == 1);
+
+ // step 3.3 x = x - q{i-t-1} * y * beta**{i-t-1}
+ t1 = y;
+ t1.multiply_by_digit(q[i - t -1]);
+ t1.shift_digits_left(i - t - 1);
+ x -= t1;
+
+ // if x < 0 then { x = x + y*beta**{i-t-1}; q{i-t-1} -= 1; }
+ if (x.is_negative())
+ {
+ t1 = y;
+ t1.shift_digits_left(i - t -1);
+ x += t1;
+
+ --q[i - t - 1] = q[i - t - 1];
+ }
+ }
+
+ // now q is the quotient and x is the remainder [which we have to normalize]
+
+ // get sign before writing to q
+ x.set_sign(!x ? 1 : a.sign());
+
+ q.clamp();
+ q.set_sign(neg);
+
+ if (remainder)
+ {
+ x >>= norm;
+ remainder->swap(x);
+ }
+}
+
+
+} // namespace detail
+} // namespace mp_math
+} // namespace boost
+
+#endif
+
Modified: sandbox/mp_math/boost/mp_math/mp_int/detail/integral_ops.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/detail/integral_ops.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/detail/integral_ops.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -1,4 +1,4 @@
-// Copyright Kevin Sopp 2008.
+// Copyright Kevin Sopp 2008 - 2009.
// 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)
@@ -42,7 +42,7 @@
static bool equal (const MpInt& lhs, IntegralT rhs);
static bool less (const MpInt& lhs, IntegralT rhs);
-
+
static void add (MpInt& lhs, IntegralT rhs);
static void subtract(MpInt& lhs, IntegralT rhs);
static void multiply(MpInt& lhs, IntegralT rhs);
@@ -82,7 +82,7 @@
static bool equal (const MpInt& lhs, IntegralT rhs);
static bool less (const MpInt& lhs, IntegralT rhs);
-
+
static void add (MpInt& lhs, IntegralT rhs);
static void subtract(MpInt& lhs, IntegralT rhs);
static void multiply(MpInt& lhs, IntegralT rhs);
@@ -98,7 +98,7 @@
typedef typename make_unsigned<IntegralT>::type unsigned_type;
static const typename MpInt::size_type q =
- std::numeric_limits<IntegralT>::digits + (MpInt::valid_bits - 1)
+ (std::numeric_limits<IntegralT>::digits + (MpInt::valid_bits - 1))
/ MpInt::valid_bits;
static bool equal_to_integral_type(const MpInt& x, IntegralT y)
@@ -121,12 +121,13 @@
integral_ops_impl<IntegralT, MpInt, false>::assign(MpInt& lhs, IntegralT rhs)
{
lhs.zero();
-
+
+ typedef typename MpInt::digit_type digit_type;
+
if (rhs <= MpInt::digit_max)
- lhs[0] = rhs;
+ lhs[0] = static_cast<digit_type>(rhs);
else
{
- typedef typename MpInt::digit_type digit_type;
static const int ud = std::numeric_limits<IntegralT>::digits;
static const int dd = MpInt::valid_bits;
static const int m = dd < ud ? dd : ud;
@@ -181,7 +182,7 @@
{
if (lhs.size() > q)
return false;
-
+
return equal_to_integral_type(lhs, rhs);
}
@@ -190,7 +191,7 @@
integral_ops_impl<IntegralT, MpInt, true>::equal(const MpInt& lhs, IntegralT rhs)
{
const int r_sign = rhs < 0 ? -1 : 1;
-
+
if (lhs.sign() != r_sign)
return false;
@@ -219,7 +220,7 @@
{
if (lhs.sign() == -1 && rhs >= 0)
return true;
-
+
if (lhs.sign() == 1 && rhs < 0)
return false;
@@ -414,7 +415,7 @@
integral_ops_impl<IntegralT, MpInt, false>::bitwise_or(MpInt& lhs, IntegralT rhs)
{
if (rhs <= MpInt::digit_max)
- lhs[0] |= rhs;
+ lhs[0] |= static_cast<typename MpInt::digit_type>(rhs);
else
lhs |= MpInt(rhs);
}
@@ -441,7 +442,7 @@
integral_ops_impl<IntegralT, MpInt, false>::bitwise_and(MpInt& lhs, IntegralT rhs)
{
if (rhs <= MpInt::digit_max)
- lhs[0] &= rhs;
+ lhs[0] &= static_cast<typename MpInt::digit_type>(rhs);
else
lhs &= MpInt(rhs);
}
@@ -467,7 +468,7 @@
integral_ops_impl<IntegralT, MpInt, false>::bitwise_xor(MpInt& lhs, IntegralT rhs)
{
if (rhs <= MpInt::digit_max)
- lhs[0] ^= rhs;
+ lhs[0] ^= static_cast<typename MpInt::digit_type>(rhs);
else
lhs ^= MpInt(rhs);
}
@@ -507,7 +508,7 @@
{
typedef std::numeric_limits<IntegralT> integral_type_limits;
typedef typename make_unsigned<IntegralT>::type unsigned_type;
-
+
static IntegralT to_integral(const MpInt& x);
};
@@ -516,7 +517,7 @@
{
typedef std::numeric_limits<IntegralT> integral_type_limits;
typedef typename make_unsigned<IntegralT>::type unsigned_type;
-
+
static IntegralT to_integral(const MpInt& x);
};
@@ -524,7 +525,7 @@
struct integral_ops2<IntegralT, MpInt, false, true>
{
typedef std::numeric_limits<IntegralT> integral_type_limits;
-
+
static IntegralT apply_shift(const MpInt& x);
static IntegralT to_integral(const MpInt& x);
};
@@ -533,7 +534,7 @@
struct integral_ops2<IntegralT, MpInt, false, false>
{
typedef std::numeric_limits<IntegralT> integral_type_limits;
-
+
static IntegralT apply_shift(const MpInt& x);
static IntegralT to_integral(const MpInt& x);
};
@@ -572,7 +573,7 @@
{
int shift_count = 0;
IntegralT tmp = 0;
-
+
for (typename MpInt::const_reverse_iterator digit = x.rbegin();
digit != x.rend(); ++digit)
{
@@ -597,7 +598,7 @@
else
throw std::overflow_error(
"to_integral: cannot convert negative number to unsigned integral type");
-
+
return apply_shift(x);
}
@@ -631,13 +632,13 @@
{
return integral_ops2<IntegralT, mp_int<A,T> >::to_integral(x);
}
-
+
template<class A, class T>
static void assign(mp_int<A,T>& lhs, IntegralT rhs)
{
integral_ops_impl<IntegralT, mp_int<A,T> >::assign(lhs, rhs);
}
-
+
template<class A, class T>
static bool equal(const mp_int<A,T>& lhs, IntegralT rhs)
{
@@ -667,7 +668,7 @@
{
integral_ops_impl<IntegralT, mp_int<A,T> >::multiply(lhs, rhs);
}
-
+
template<class A, class T>
static void divide(mp_int<A,T>& lhs, IntegralT rhs)
{
Modified: sandbox/mp_math/boost/mp_math/mp_int/detail/modular_reduction.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/detail/modular_reduction.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/detail/modular_reduction.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -103,7 +103,8 @@
// x = x/base**m.size()
x.clamp();
- if (x.is_zero())
+
+ if (!x)
x.set_sign(1);
x.shift_digits_right(m.size());
@@ -186,7 +187,8 @@
std::memset(x.digits() + m + 1, 0, (x.size() - (m + 1)) * sizeof(digit_type));
x.clamp();
- if (x.is_zero())
+
+ if (!x)
x.set_sign(1);
if (x.compare_magnitude(n) != -1)
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 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -20,17 +20,17 @@
typedef DigitT digit_type;
typedef WordT word_type;
typedef SizeT size_type;
-
+
static const word_type digit_bits = std::numeric_limits<digit_type>::digits;
// ADD ------------------------------------
- // add x to the digits in y and store result in z
- // ylen must be > 0
+ // add y to the digits in x and store result in z
+ // xlen must be > 0
// returns: the last carry (it will not get stored in z)
static digit_type add_single_digit(digit_type* z,
- const digit_type* y, size_type ylen,
- digit_type x);
+ const digit_type* x, size_type xlen,
+ digit_type y);
// z = x + y, returns last carry
static digit_type add_digits(digit_type* z,
@@ -56,8 +56,8 @@
// subtracts x from the digits in y and store result in z
static void subtract_single_digit(digit_type* z,
- const digit_type* y, size_type ylen,
- digit_type x);
+ const digit_type* x, size_type xlen,
+ digit_type y);
// z = x - y, returns last borrow
static digit_type subtract_digits(digit_type* z,
@@ -82,8 +82,12 @@
// multiply y of length ylen with x and store result in z
// returns: the last carry (it will not get stored in z)
static digit_type multiply_by_digit(digit_type* z,
- const digit_type* y, size_type ylen,
- digit_type x);
+ const digit_type* x, size_type xlen,
+ digit_type y);
+
+ // z = x * 2
+ static digit_type multiply_by_two(digit_type* z,
+ const digit_type* x, size_type len);
// z = x * y
static void classic_mul(digit_type* z, const digit_type* x, size_type xlen,
@@ -111,24 +115,35 @@
digit_type x,
const digit_type* y,
size_type n);
+
+ // DIV -------------------------------------
+
+ // z = x / 2
+ static void divide_by_two(digit_type* z, const digit_type* x, size_type len);
+
+ // z = x / y
+ // returns remainder
+ static digit_type divide_by_digit(digit_type* z,
+ const digit_type* x, size_type xlen,
+ digit_type y);
};
template<typename D, typename W, typename S>
inline
-typename basic_primitive_ops<D,W,S>::digit_type
-basic_primitive_ops<D,W,S>::add_single_digit(digit_type* z,
- const digit_type* y, size_type ylen,
- digit_type x)
+D basic_primitive_ops<D,W,S>::add_single_digit(digit_type* z,
+ const digit_type* x,
+ size_type xlen,
+ digit_type y)
{
- word_type carry = static_cast<word_type>(*y++) + x;
+ word_type carry = static_cast<word_type>(*x++) + y;
*z++ = static_cast<digit_type>(carry);
carry >>= digit_bits;
- while (carry && --ylen)
+ while (carry && --xlen)
{
- carry += static_cast<word_type>(*y++);
+ carry += static_cast<word_type>(*x++);
*z++ = static_cast<digit_type>(carry);
carry >>= digit_bits;
}
@@ -139,10 +154,9 @@
template<typename D, typename W, typename S>
inline
-typename basic_primitive_ops<D,W,S>::digit_type
-basic_primitive_ops<D,W,S>::add_digits(digit_type* z,
- const digit_type* x,
- const digit_type* y, size_type n)
+D basic_primitive_ops<D,W,S>::add_digits(digit_type* z,
+ const digit_type* x,
+ const digit_type* y, size_type n)
{
word_type carry = 0;
@@ -158,11 +172,10 @@
template<typename D, typename W, typename S>
inline
-typename basic_primitive_ops<D,W,S>::size_type
-basic_primitive_ops<D,W,S>::ripple_carry(digit_type* z,
- const digit_type* x,
- size_type n,
- digit_type& carry)
+S basic_primitive_ops<D,W,S>::ripple_carry(digit_type* z,
+ const digit_type* x,
+ size_type n,
+ digit_type& carry)
{
word_type c = carry;
size_type i = 0;
@@ -175,27 +188,27 @@
}
carry = static_cast<digit_type>(c);
-
+
return i;
}
template<typename D, typename W, typename S>
inline
-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)
+D 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);
n += ylen;
-
+
if (n < xlen && z != x)
std::memcpy(z + n, x + n, sizeof(digit_type) * (xlen - n));
-
+
return carry;
}
@@ -203,7 +216,8 @@
inline
void
basic_primitive_ops<D,W,S>::subtract_single_digit(digit_type* z,
- const digit_type* y, size_type ylen,
+ const digit_type* y,
+ size_type ylen,
digit_type x)
{
word_type borrow = static_cast<word_type>(*y++) - x;
@@ -220,11 +234,10 @@
template<typename D, typename W, typename S>
inline
-typename basic_primitive_ops<D,W,S>::digit_type
-basic_primitive_ops<D,W,S>::subtract_digits(digit_type* z,
- const digit_type* x,
- const digit_type* y,
- size_type n)
+D basic_primitive_ops<D,W,S>::subtract_digits(digit_type* z,
+ const digit_type* x,
+ const digit_type* y,
+ size_type n)
{
word_type borrow = 0;
@@ -240,11 +253,10 @@
template<typename D, typename W, typename S>
inline
-typename basic_primitive_ops<D,W,S>::size_type
-basic_primitive_ops<D,W,S>::ripple_borrow(digit_type* z,
- const digit_type* x,
- size_type n,
- digit_type borrow)
+S basic_primitive_ops<D,W,S>::ripple_borrow(digit_type* z,
+ const digit_type* x,
+ size_type n,
+ digit_type borrow)
{
word_type b = borrow;
size_type i = 0;
@@ -254,7 +266,7 @@
*z++ = static_cast<digit_type>(b);
b >>= std::numeric_limits<word_type>::digits - 1;
}
-
+
return i;
}
@@ -276,13 +288,12 @@
}
}
-
template<typename D, typename W, typename S>
inline
-typename basic_primitive_ops<D,W,S>::digit_type
-basic_primitive_ops<D,W,S>::multiply_by_digit(digit_type* z,
- const digit_type* y, size_type ylen,
- digit_type x)
+D basic_primitive_ops<D,W,S>::multiply_by_digit(digit_type* z,
+ const digit_type* y,
+ size_type ylen,
+ digit_type x)
{
digit_type carry = 0;
@@ -294,6 +305,29 @@
*z++ = static_cast<digit_type>(tmp);
carry = static_cast<digit_type>(tmp >> digit_bits);
}
+
+ return carry;
+}
+
+template<typename D, typename W, typename S>
+inline
+D basic_primitive_ops<D,W,S>::multiply_by_two(digit_type* z,
+ const digit_type* x, size_type n)
+{
+ static const digit_type one = 1U;
+
+ digit_type carry = 0;
+
+ while (n--)
+ {
+ // get carry bit for next iteration
+ const digit_type r = *x >> (static_cast<digit_type>(digit_bits) - one);
+
+ *z++ = (*x++ << one) | carry;
+
+ carry = r;
+ }
+
return carry;
}
@@ -306,7 +340,7 @@
// phase 1
word_type tmp = static_cast<word_type>(x[0]) * static_cast<word_type>(y[0]);
z[0] = static_cast<digit_type>(tmp);
-
+
for (size_type i = 1; i < xlen; ++i)
{
tmp = (tmp >> digit_bits)
@@ -314,10 +348,10 @@
* static_cast<word_type>(y[0]);
z[i] = static_cast<digit_type>(tmp);
}
-
+
tmp >>= digit_bits;
z[xlen] = static_cast<digit_type>(tmp);
-
+
// phase 2
for (size_type i = 1; i < ylen; ++i)
{
@@ -325,7 +359,7 @@
* static_cast<word_type>(x[0])
+ static_cast<word_type>(z[i]);
z[i] = static_cast<digit_type>(tmp);
-
+
for (size_type j = 1; j < xlen; ++j)
{
tmp = (tmp >> digit_bits)
@@ -333,7 +367,7 @@
+ static_cast<word_type>(z[i+j]);
z[i+j] = static_cast<digit_type>(tmp);
}
-
+
tmp >>= digit_bits;
z[i + xlen] = static_cast<digit_type>(tmp);
}
@@ -348,8 +382,8 @@
//----------------------------------
// | 18 15 | 12 9 6
// 12 | 10 8 | 6 4
-// 12 10 | 8 6 | 4
-// | |
+// 12 10 | 8 6 | 4
+// | |
// phase: 3 | 2 | 1
//
// TODO handle too long columns => carry may overflow. This can happen if
@@ -365,7 +399,7 @@
word_type acc = 0; // accumulator for each column
word_type carry = 0;
-
+
// phase 1
for (size_type i = 0; i < ylen; ++i)
{
@@ -401,7 +435,7 @@
acc = static_cast<digit_type>(carry);
carry >>= digit_bits;
}
-
+
// phase 3
for (size_type i = 0; i < ylen - 1; ++i)
{
@@ -477,6 +511,54 @@
const digit_type* x,
size_type xlen)
{
+/* word_type acc = 0;
+ word_type carry = 0;
+ word_type acc2;
+
+ for (size_type i = 0; i < xlen; ++i)
+ {
+ // even colum
+ acc += static_cast<word_type>(x[i]) * static_cast<word_type>(x[i]);
+
+ const digit_type* a = x + i;
+ const digit_type* b = x + i;
+
+ acc2 = 0;
+ for (size_type j = 0; j < (i - n&1) >> 1; ++j)
+ {
+ acc2 += static_cast<word_type>(*(--a)) * static_cast<word_type>(*(--b));
+ carry += acc2 >> digit_bits;
+ acc2 = static_cast<digit_type>(acc2);
+ }
+
+ acc += acc2 + acc2;
+
+ *z++ = static_cast<digit_type>(acc);
+ acc = static_cast<digit_type>(carry);
+ carry >>= digit_bits;
+
+ // odd column
+ const digit_type* a = x + i;
+ const digit_type* b = x + i + 1;
+
+ acc2 = 0;
+ for (size_type j = 0; j <= i; ++j)
+ {
+ acc2 += static_cast<word_type>(*a++) * static_cast<word_type>(*b--);
+ carry += acc2 >> digit_bits;
+ acc2 = static_cast<digit_type>(acc2);
+ }
+
+ acc += acc2 + acc2;
+
+ *z++ = static_cast<digit_type>(acc);
+ acc = static_cast<digit_type>(carry);
+ carry >>= digit_bits;
+ }
+
+ *z = static_cast<digit_type>(acc);*/
+
+
word_type acc = 0; // accumulator for each column
word_type carry = 0;
@@ -520,12 +602,11 @@
}
template<typename D, typename W, typename S>
-typename basic_primitive_ops<D,W,S>::digit_type
-basic_primitive_ops<D,W,S>::multiply_add_digits(digit_type* z,
- const digit_type* w,
- digit_type x,
- const digit_type* y,
- size_type n)
+D basic_primitive_ops<D,W,S>::multiply_add_digits(digit_type* z,
+ const digit_type* w,
+ digit_type x,
+ const digit_type* y,
+ size_type n)
{
word_type carry = 0;
while (n--)
@@ -543,6 +624,55 @@
}
+template<typename D, typename W, typename S>
+inline
+void basic_primitive_ops<D,W,S>::divide_by_two(digit_type* z,
+ const digit_type* x, size_type n)
+{
+ z += n - 1;
+ x += n - 1;
+
+ digit_type carry = 0;
+
+ while (n--)
+ {
+ // get carry for next iteration
+ const digit_type r = *x & 1;
+
+ *z-- = (*x-- >> 1) | (carry << (digit_bits - 1));
+
+ carry = r;
+ }
+}
+
+template<typename D, typename W, typename S>
+D basic_primitive_ops<D,W,S>::divide_by_digit(digit_type* z,
+ const digit_type* x, size_type n,
+ digit_type y)
+{
+ z += n - 1;
+ x += n - 1;
+
+ word_type w = 0;
+
+ while (n--)
+ {
+ w = (w << digit_bits) | static_cast<word_type>(*x--);
+ digit_type tmp;
+ if (w >= y)
+ {
+ tmp = static_cast<digit_type>(w / y);
+ w -= tmp * y;
+ }
+ else
+ tmp = 0;
+ *z-- = tmp;
+ }
+
+ return static_cast<digit_type>(w);
+}
+
+
// This exists to ease development of primitive_ops specializations. If a
@@ -558,7 +688,7 @@
// Here we include primitive_ops specializations that use assembler
#if defined(BOOST_MP_MATH_MP_INT_USE_ASM)
-
+
#if defined(__GNU__)
#if defined(__386__)
#include <boost/mp_math/mp_int/detail/asm/x86/gnu_386_primitive_ops.hpp>
Modified: sandbox/mp_math/boost/mp_math/mp_int/div.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/div.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/div.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -10,61 +10,42 @@
if (b == 0)
throw std::domain_error("mp_int::divide_by_digit: division by zero");
- if (b == 1 || is_zero())
+ if (b == 1 || !*this)
return 0;
const bool is_power_of_two = (b & (b-1)) == 0;
- if (is_power_of_two)
- {
- for (int i = 0; i < valid_bits; ++i)
- {
- if (b == digit_type(1) << i)
- {
- const digit_type remainder = digits_[0] & ((digit_type(1) << i) - 1);
- *this >>= i;
- return remainder;
- }
- }
- }
- word_type w = 0;
- for (reverse_iterator d = rbegin(); d != rend(); ++d)
+ if (!is_power_of_two)
{
- w = (w << static_cast<word_type>(valid_bits)) | static_cast<word_type>(*d);
- digit_type tmp;
- if (w >= b)
- {
- tmp = static_cast<digit_type>(w / b);
- w -= tmp * b;
- }
- else
- tmp = 0;
- *d = tmp;
+ const digit_type remainder =
+ ops_type::divide_by_digit(digits(), digits(), size(), b);
+
+ clamp_high_digit();
+
+ if (!*this)
+ set_sign(1);
+
+ return remainder;
}
- clamp();
- if (is_zero())
- set_sign(1);
-
- return static_cast<digit_type>(w);
+ int i = 0;
+ while ((i < valid_bits) && (b != digit_type(1) << i))
+ ++i;
+
+ const digit_type remainder = digits_[0] & ((digit_type(1) << i) - 1);
+ *this >>= i;
+ return remainder;
}
// *this /= 2
template<class A, class T>
void mp_int<A,T>::divide_by_2()
{
- digit_type carry = 0;
- for (reverse_iterator d = rbegin(); d != rend(); ++d)
- {
- // get the carry for the next iteration
- const digit_type rr = *d & 1;
- // shift the current digit, add in carry and store
- *d = (*d >> 1) | (carry << (valid_bits - 1));
- // forward carry to next iteration
- carry = rr;
- }
- clamp();
- if (is_zero())
+ ops_type::divide_by_two(digits(), digits(), size());
+
+ clamp_high_digit();
+
+ if (!*this)
set_sign(1);
}
@@ -154,165 +135,15 @@
carry = rr;
}
}
- clamp();
- if (is_zero())
- set_sign(1);
-}
-
-// integer signed division.
-// c*b + d == a [e.g. a/b, c=quotient, d=remainder]
-// HAC pp.598 Algorithm 14.20
-//
-// Note that the description in HAC is horribly incomplete. For example, it
-// doesn't consider the case where digits are removed from 'x' in the inner
-// loop. It also doesn't consider the case that y has fewer than three digits,
-// etc..
-// The overall algorithm is as described as 14.20 from HAC but fixed to treat
-// these cases.
-// divide *this by rhs, optionally store remainder
-template<class A, class T>
-void mp_int<A,T>::divide(const mp_int& rhs, mp_int* remainder)
-{
- if (rhs.is_zero())
- throw std::domain_error("mp_int::divide: division by zero");
-
- // if *this < rhs then q=0, r = *this
- if (compare_magnitude(rhs) == -1)
- {
- if (remainder)
- *remainder = *this;
- zero();
- return;
- }
-
- mp_int q;
- q.grow_capacity(size_ + 2);
- q.size_ = size_ + 2;
- std::memset(q.digits_, 0, q.size_ * sizeof(digit_type));
-
- mp_int x(*this);
- mp_int y(rhs);
-
- // fix the sign
- const int neg = (sign() == rhs.sign()) ? 1 : -1;
- x.set_sign(1);
- y.set_sign(1);
-
- // normalize both x and y, ensure that y >= beta/2, [beta == 2**valid_bits]
- size_type norm = y.precision() % valid_bits;
- if (norm < valid_bits-1)
- {
- norm = valid_bits - 1 - norm;
- x <<= norm;
- y <<= norm;
- }
- else
- norm = 0;
-
- // note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4
- const size_type n = x.size_ - 1;
- const size_type t = y.size_ - 1;
-
- // find leading digit of the quotient
- // while (x >= y*beta**(n-t)) do { q[n-t] += 1; x -= y*beta**(n-t) }
- y.shift_digits_left(n - t); // y = y*beta**(n-t)
-
- while (x.compare(y) != -1)
- {
- ++q[n - t];
- x -= y;
- }
-
- // reset y by shifting it back down
- y.shift_digits_right(n - t);
-
- // find the remainder of the digits
- // step 3. for i from n down to (t + 1)
- for (size_type i = n; i >= (t + 1); i--)
- {
- if (i > x.size_)
- continue;
-
- // step 3.1 if xi == yt then set q{i-t-1} to beta-1,
- // otherwise set q{i-t-1} to (xi*beta + x{i-1})/yt
- if (x[i] == y[t])
- q[i - t - 1] = std::numeric_limits<digit_type>::max();
- else
- {
- word_type tmp = static_cast<word_type>(x[i])
- << static_cast<word_type>(valid_bits);
- tmp |= x[i - 1];
- tmp /= y[t];
- q[i - t - 1] = static_cast<digit_type>(tmp);
- }
-
- // now fixup quotient estimation
- // while (q{i-t-1} * (yt * beta + y{t-1})) >
- // xi * beta**2 + xi-1 * beta + xi-2
- //
- // do q{i-t-1} -= 1;
-
- mp_int t1, t2;
- t1.grow_capacity(3);
- t2.grow_capacity(3);
-
- ++q[i - t - 1];
- do
- {
- --q[i - t - 1];
-
- // find left hand
- t1.zero();
- t1[0] = (t == 0) ? 0 : y[t - 1];
- t1[1] = y[t];
- t1.size_ = 2;
- t1.multiply_by_digit(q[i - t - 1]);
-
- // find right hand
- t2[0] = (i < 2) ? 0 : x[i - 2];
- t2[1] = (i == 0) ? 0 : x[i - 1];
- t2[2] = x[i];
- t2.size_ = 3;
- } while (t1.compare_magnitude(t2) == 1);
-
- // step 3.3 x = x - q{i-t-1} * y * beta**{i-t-1}
- t1 = y;
- t1.multiply_by_digit(q[i - t -1]);
- t1.shift_digits_left(i - t - 1);
- x -= t1;
-
- // if x < 0 then { x = x + y*beta**{i-t-1}; q{i-t-1} -= 1; }
- if (x.is_negative())
- {
- t1 = y;
- t1.shift_digits_left(i - t -1);
- x += t1;
-
- --q[i - t - 1] = q[i - t - 1];
- }
- }
-
- // now q is the quotient and x is the remainder [which we have to normalize]
-
- // get sign before writing to *this
- x.set_sign(x.is_zero() ? 1 : sign());
-
- q.clamp();
- swap(q);
- set_sign(neg);
+ clamp();
- if (remainder)
- {
- x >>= norm;
- remainder->swap(x);
- }
+ if (!*this)
+ set_sign(1);
}
-/*
template<class A, class T>
void divide(const mp_int<A,T>& x, const mp_int<A,T>& y, mp_int<A,T>& q, mp_int<A,T>& r)
{
- divide(x, y, &q, &r);
+ detail::classic_divide(x, y, q, &r);
}
-*/
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 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -23,7 +23,7 @@
if (p <= digit_type(0))
throw std::domain_error("jacobi: p must be greater than 0");
- if (a.is_zero())
+ if (!a)
return 0;
if (a == digit_type(1))
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 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -20,7 +20,8 @@
digits_[b / valid_bits] &= mask;
clamp();
- if (is_zero())
+
+ if (!*this)
set_sign(1);
}
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 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -1,4 +1,4 @@
-// Copyright Kevin Sopp 2008.
+// Copyright Kevin Sopp 2008 - 2009.
// 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)
@@ -20,6 +20,7 @@
#include <boost/random.hpp>
#include <boost/type_traits/is_integral.hpp>
#include <boost/utility/enable_if.hpp>
+#include <boost/mp_math/mp_int/detail/div.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>
@@ -29,8 +30,6 @@
namespace boost {
namespace mp_math {
-// digits are stored in least significant order
-
template<
class Allocator,
class Traits
@@ -39,40 +38,60 @@
:
Allocator::template rebind<typename Traits::digit_type>::other
{
- typedef Allocator allocator_type;
- typedef Traits traits_type;
- typedef typename Allocator::size_type size_type;
+private:
+
+ typedef typename Allocator::template
+ rebind<typename Traits::digit_type>::other base_allocator_type;
+
+public:
+
+ typedef Allocator allocator_type;
+ typedef Traits traits_type;
+ typedef typename base_allocator_type::size_type size_type;
mp_int();
+ explicit mp_int(const allocator_type& a);
+
template<typename IntegralT>
- mp_int(IntegralT, typename enable_if<is_integral<IntegralT> >::type* dummy = 0);
+ mp_int(IntegralT,
+ const allocator_type& a = allocator_type(),
+ typename enable_if<is_integral<IntegralT> >::type* dummy = 0);
template<typename charT>
- mp_int(const charT*);
+ mp_int(const charT*, const allocator_type& a = allocator_type());
template<typename charT>
- mp_int(const charT*, std::ios_base::fmtflags);
+ mp_int(const charT*,
+ std::ios_base::fmtflags,
+ const allocator_type& a = allocator_type());
template<typename charT, class traits, class Alloc>
- mp_int(const std::basic_string<charT,traits,Alloc>&);
-
+ mp_int(const std::basic_string<charT,traits,Alloc>&,
+ const allocator_type& a = allocator_type());
+
template<typename charT, class traits, class Alloc>
- mp_int(const std::basic_string<charT,traits,Alloc>&, std::ios_base::fmtflags);
-
+ mp_int(const std::basic_string<charT,traits,Alloc>&,
+ std::ios_base::fmtflags,
+ const allocator_type& a = allocator_type());
+
template<typename RandomAccessIterator>
- mp_int(RandomAccessIterator first, RandomAccessIterator last);
-
+ mp_int(RandomAccessIterator first,
+ RandomAccessIterator last,
+ const allocator_type& a = allocator_type());
+
template<typename RandomAccessIterator>
- mp_int(RandomAccessIterator first, RandomAccessIterator last,
- std::ios_base::fmtflags f);
+ mp_int(RandomAccessIterator first,
+ RandomAccessIterator last,
+ std::ios_base::fmtflags f,
+ const allocator_type& a = allocator_type());
mp_int(const mp_int& copy);
#ifdef BOOST_HAS_RVALUE_REFS
mp_int(mp_int&& copy);
#endif
-
+
~mp_int();
mp_int& operator = (const mp_int& rhs);
@@ -80,7 +99,7 @@
#ifdef BOOST_HAS_RVALUE_REFS
mp_int& operator = (mp_int&& rhs);
#endif
-
+
template<typename IntegralT>
mp_int& operator = (IntegralT rhs);
@@ -95,8 +114,8 @@
template<typename charT, class traits, class Alloc>
void assign(const std::basic_string<charT,traits,Alloc>&,
- std::ios_base::fmtflags);
-
+ std::ios_base::fmtflags);
+
template<typename RandomAccessIterator>
void assign(RandomAccessIterator first, RandomAccessIterator last,
std::ios_base::fmtflags);
@@ -169,7 +188,7 @@
operator unspecified_bool_type() const
{
- return is_zero() ? 0 : &mp_int::size_;
+ return !(size_ == 1 && digits_[0] == 0) ? &mp_int::size_ : 0;
}
bool is_even() const { return (digits_[0] & digit_type(1)) == 0; }
@@ -224,13 +243,14 @@
digit_type& operator[](size_type i) { return digits_[i]; }
const digit_type& operator[](size_type i) const { return digits_[i]; }
-
+
digit_type& at(size_type i)
{
if (i >= size_)
throw std::out_of_range("mp_int::at: array subscript out of range");
return digits_[i];
}
+
const digit_type& at(size_type i) const
{
if (i >= size_)
@@ -247,7 +267,7 @@
void print(bool all=false) const;
bool test_invariants() const;
- bool is_uninitialized() const { return size_ == 0U; }
+ bool is_uninitialized() const { return !size_; }
size_type size() const { return size_; }
size_type capacity() const
@@ -281,6 +301,7 @@
void grow_capacity(size_type n);
void clamp();
+ void clamp_high_digit();
int compare_magnitude(const mp_int& rhs) const;
int compare_to_digit(digit_type) const;
@@ -289,7 +310,6 @@
void add_magnitude(const mp_int& rhs);
void sub_smaller_magnitude(const mp_int& rhs);
- bool is_zero() const;
bool is_power_of_two() const;
void add_digit(digit_type);
void sub_digit(digit_type);
@@ -312,7 +332,6 @@
void comba_sqr();
digit_type divide_by_digit(digit_type); // returns remainder
- void divide(const mp_int& divisor, mp_int* remainder);
void divide_by_2();
digit_type divide_by_3();
void modulo_2_to_the_power_of(size_type);
@@ -326,13 +345,22 @@
{
digits_[0] |= digit_type(1);
}
-
+
void set_bit(size_type bit)
{
- digits_[bit / valid_bits] |=
- digit_type(1) << digit_type(bit % valid_bits);
+ digits_[bit / valid_bits] |= digit_type(1) << (bit % valid_bits);
}
+ void clear_bit(size_type bit)
+ {
+ digits_[bit / valid_bits] &= ~(digit_type(1) << (bit % valid_bits));
+ }
+
+ void set_bits(size_type beg, size_type end);
+ void clear_bits(size_type beg, size_type end);
+
+ void truncate(size_type prec);
+
template<class A, class T>
friend bool operator == (const mp_int<A,T>&, const mp_int<A,T>&);
@@ -343,7 +371,7 @@
void from_string(Iter first, Iter last, unsigned radix);
private:
-
+
digit_type* digits_;
size_type size_, capacity_;
};
@@ -364,7 +392,7 @@
cout << ",";
}
cout << "}";
-
+
if (all)
{
cout << capacity() - size_ << "{";
@@ -388,7 +416,7 @@
return false;
if (digits_[size_-1] == 0U)
return false;
- if (is_zero() && sign() != 1)
+ if (!*this && sign() != 1)
return false;
}
return true;
@@ -417,7 +445,8 @@
{
if (this != &rhs)
{
- this->deallocate(digits_, capacity());
+ if (digits_)
+ this->deallocate(digits_, capacity());
digits_ = 0;
size_ = 0;
capacity_ = 0;
@@ -592,15 +621,16 @@
template<class A, class T>
void mp_int<A,T>::clamp()
{
- // decrease size_ while the most significant digit is zero.
while (size_ > 1 && digits_[size_-1] == 0)
--size_;
}
+// For when we know that only one leading zero digit may exist.
template<class A, class T>
-inline bool mp_int<A,T>::is_zero() const
+inline void mp_int<A,T>::clamp_high_digit()
{
- return size_ == 1 && digits_[0] == 0;
+ if (size_ > 1 && digits_[size_-1] == 0)
+ --size_;
}
// disregards the sign of the numbers
@@ -613,7 +643,7 @@
// compare based on # of non-zero digits
if (size_ > rhs.size_)
return 1;
-
+
if (size_ < rhs.size_)
return -1;
@@ -660,7 +690,7 @@
else
return 1;
}
-
+
if (is_negative())
// if negative compare opposite direction
return rhs.compare_magnitude(*this);
@@ -705,7 +735,7 @@
// zero the top digits
std::memset(digits_ + size_ - b, 0, b * sizeof(digit_type));
-
+
// remove excess digits
size_ -= b;
}
@@ -716,7 +746,7 @@
{
// get number of digits and add that
size_type r = (size_ - 1) * valid_bits;
-
+
// take the last digit and count the bits in it
digit_type q = digits_[size_ - 1];
while (q > 0U)
@@ -736,7 +766,7 @@
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
};
- if (is_zero())
+ if (!*this)
return 0;
// scan lower digits until non-zero
@@ -767,6 +797,67 @@
return detail::integral_ops<IntegralT>::convert(*this);
}
+template<class A, class T>
+void mp_int<A,T>::set_bits(size_type beg, size_type end)
+{
+ const size_type beg_index = beg / digit_bits;
+ const size_type end_index = end / digit_bits;
+ const size_type first_bits = beg % digit_bits;
+ const size_type last_bits = end % digit_bits;
+
+ static const digit_type z = ~digit_type(0);
+
+ digit_type mask = z << first_bits;
+ if (beg_index == end_index && last_bits)
+ mask &= z >> (digit_bits - last_bits);
+
+ digits_[beg_index] |= mask;
+
+ for (size_type i = beg_index + ((beg % digit_bits) ? 1 : 0); i < end_index; ++i)
+ digits_[i] = digit_max;
+
+ if (beg_index != end_index && last_bits)
+ digits_[end_index] |= z >> (digit_bits - last_bits);
+}
+
+template<class A, class T>
+void mp_int<A,T>::clear_bits(size_type beg, size_type end)
+{
+ const size_type beg_index = beg / digit_bits;
+ const size_type end_index = end / digit_bits;
+ const size_type first_bits = beg % digit_bits;
+ const size_type last_bits = end % digit_bits;
+
+ static const digit_type z = ~digit_type(0);
+
+ digit_type mask;
+ if (first_bits)
+ mask = z >> (digit_bits - first_bits);
+ else
+ mask = 0;
+
+ if (beg_index == end_index)
+ mask |= z << last_bits;
+
+ digits_[beg_index] &= mask;
+
+ if (beg_index != end_index)
+ {
+ std::memset(digits_ + beg_index + 1, 0,
+ sizeof(digit_type) * (end_index - beg_index - 1));
+
+ digits_[end_index] &= z << last_bits;
+ }
+}
+
+template<class A, class T>
+void mp_int<A,T>::truncate(size_type prec)
+{
+ set_size((prec + digit_bits - 1)/digit_bits);
+ const size_type last_bits = prec % digit_bits;
+ digits_[size()] &= ~digit_type(0) << (digit_bits - last_bits);
+}
+
template<class A, class T>
inline void swap(mp_int<A,T>& lhs, mp_int<A,T>& rhs)
Modified: sandbox/mp_math/boost/mp_math/mp_int/mul.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/mul.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/mul.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -18,7 +18,8 @@
// make sure we can hold the result
grow_capacity(size_ + 1);
- const digit_type carry = ops_type::multiply_by_digit(digits_, digits_, size_, x);
+ const digit_type carry =
+ ops_type::multiply_by_digit(digits(), digits(), size(), x);
if (carry)
push(carry);
@@ -30,27 +31,11 @@
{
grow_capacity(size_ + 1);
- digit_type carry = 0;
- for (iterator d = begin(); d != end(); ++d)
- {
- // get what will be the *next* carry bit from the
- // MSB of the current digit
- const digit_type rr = *d >> (static_cast<digit_type>(valid_bits - 1));
-
- // now shift up this digit, add in the carry [from the previous]
- *d = (*d << digit_type(1)) | carry;
-
- // copy the carry that would be from the source
- // digit into the next iteration
- carry = rr;
- }
+ const digit_type carry =
+ ops_type::multiply_by_two(digits(), digits(), size());
- // new leading digit?
if (carry)
- {
- // add a MSB which is always 1 at this point
- push(1);
- }
+ push(carry);
}
@@ -310,8 +295,10 @@
}
tmp.clamp();
- if (tmp.is_zero())
+
+ if (!tmp)
tmp.set_sign(1);
+
swap(tmp);
}
@@ -350,8 +337,10 @@
}
tmp.clamp();
- if (tmp.is_zero())
+
+ if (!tmp)
tmp.set_sign(1);
+
swap(tmp);
}
Modified: sandbox/mp_math/boost/mp_math/mp_int/operators.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/operators.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/operators.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -283,7 +283,7 @@
}
if (carry)
- digits_[size_++] = carry;
+ push(carry);
}
clamp();
@@ -301,7 +301,7 @@
template<class A, class T>
inline mp_int<A,T>& mp_int<A,T>::operator - ()
{
- if (!is_zero())
+ if (*this)
set_sign(sign() == 1 ? -1 : 1);
return *this;
}
@@ -341,7 +341,8 @@
size_ = x->size_;
clamp();
- if (is_zero())
+
+ if (!*this)
set_sign(1);
}
return *this;
@@ -378,7 +379,8 @@
size_ = x->size_;
clamp();
- if (is_zero())
+
+ if (!*this)
set_sign(1);
}
return *this;
@@ -420,18 +422,19 @@
tmp.size_ = size_ + rhs.size_;
tmp.clamp();
- swap(tmp);
+ *this = tmp;
}
- set_sign(is_zero() ? 1 : neg);
+ set_sign(!*this ? 1 : neg);
return *this;
}
template<class A, class T>
-inline mp_int<A,T>& mp_int<A,T>::operator /= (const mp_int<A,T>& rhs)
+mp_int<A,T>& mp_int<A,T>::operator /= (const mp_int<A,T>& rhs)
{
- divide(rhs, 0);
+ const mp_int<A,T> tmp(*this);
+ detail::classic_divide(tmp, rhs, *this);
return *this;
}
@@ -439,9 +442,8 @@
template<class A, class T>
mp_int<A,T>& mp_int<A,T>::operator %= (const mp_int<A,T>& rhs)
{
- mp_int<A,T> remainder;
- divide(rhs, &remainder);
- swap(remainder);
+ mp_int<A,T> quotient;
+ detail::classic_divide(*this, rhs, quotient, this);
return *this;
}
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 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -1,4 +1,4 @@
-// Copyright Kevin Sopp 2008.
+// Copyright Kevin Sopp 2008 - 2009.
// 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)
@@ -20,7 +20,7 @@
mp_int<A,T> t1;
- if (x.is_zero())
+ if (!x)
{
t1.zero();
return t1;
@@ -49,35 +49,39 @@
}
-// 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].
+// Uses Newton-Raphson approximation.
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)
{
if ((n & 1) == 0 && x.is_negative())
throw std::domain_error("nth_root: argument must be positive if n is even");
+ if (n == 0U)
+ throw std::domain_error("nth_root: n must not be zero");
+ else if (n == 1U)
+ return x;
+
// 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);
mp_int<A,T> t1, t2, t3;
-
- t2 = typename mp_int<A,T>::digit_type(2);
+
+ typedef typename mp_int<A,T>::size_type size_type;
+
+ // initial approximation
+ const size_type result_precision = (x.precision() - 1) / n + 1;
+ t2.grow_capacity(1);
+ t2.set_size(1);
+ t2[0] = 0;
+ t2.set_bits(0, result_precision + 1);
do
{
t1 = t2;
// t2 = t1 - ((t1**n - x) / (n * t1**(n-1)))
-
+
// t3 = t1**(n-1)
t3 = pow(t1, n-1);
@@ -97,7 +101,7 @@
t2 = t1 - t3;
} while (t1 != t2);
-
+
// result can be off by a few so check
for (;;)
{
@@ -132,15 +136,26 @@
const_cast<mp_int<A,T>*>(&x)->set_sign(1);
mp_int<A,T> t1, t2, t3;
-
- t2 = typename mp_int<A,T>::digit_type(2);
+
+ typedef typename mp_int<A,T>::size_type size_type;
+
+ static const size_type digit_bits = mp_int<A,T>::digit_bits;
+
+ const size_type result_precision = (x.precision() - 1)
+ / n.template to_integral<size_type>() + 1;
+
+ t2.grow_capacity((result_precision + digit_bits - 1) / digit_bits);
+ t2.set_size ((result_precision + digit_bits - 1) / digit_bits);
+
+ t2[t2.size()-1] = 0;
+ t2.set_bits(0, result_precision + 1);
do
{
t1 = t2;
// t2 = t1 - ((t1**n - x) / (n * t1**(n-1)))
-
+
// t3 = t1**(n-1)
t3 = pow(t1, n-1);
@@ -160,7 +175,7 @@
t2 = t1 - t3;
} while (t1 != t2);
-
+
// result can be off by a few so check
for (;;)
{
Modified: sandbox/mp_math/boost/mp_math/mp_int/sqr.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/sqr.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/sqr.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -174,15 +174,35 @@
template<class A, class T>
void mp_int<A,T>::comba_sqr()
{
- mp_int tmp;
- tmp.grow_capacity(size_ + size_);
+ if (size() < 16U)
+ {
+ grow_capacity(size() + size());
+
+ digit_type Z[32];
+
+ ops_type::comba_sqr(Z, digits(), size());
+
+ std::memcpy(digits(), Z, (size() + size()) * sizeof(digit_type));
+
+ set_size(size() + size());
+
+ if (!digits()[size()-1])
+ pop();
+ }
+ else
+ {
+ mp_int tmp;
+ tmp.grow_capacity(size() + size());
+
+ ops_type::comba_sqr(tmp.digits(), digits(), size());
+
+ tmp.set_size(size() + size());
- ops_type::comba_sqr(tmp.digits(), digits(), size_);
+ if (!tmp[tmp.size()-1])
+ tmp.pop();
- tmp.size_ = size_ + size_;
-
- tmp.clamp();
- swap(tmp);
+ *this = tmp;
+ }
}
// computes *this = *this * *this
Modified: sandbox/mp_math/boost/mp_math/mp_int/string_conversion.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/string_conversion.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/string_conversion.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -64,17 +64,18 @@
if (offset >= valid_bits)
{
- digits_[size_++] = result;
+ push(result);
offset -= valid_bits;
result = static_cast<digit_type>(x >> (sc.radix_storage_bits - offset));
}
}
- if (result || !size_)
- digits_[size_++] = result;
+ if (result || is_uninitialized())
+ push(result);
clamp();
- if (is_zero())
+
+ if (!*this)
set_sign(1);
}
else // radix can only be 10 at this point
@@ -112,16 +113,16 @@
carry += ops_type::add_single_digit(digits_, digits_, size_, result);
if (carry)
- digits_[size_++] = carry;
+ push(carry);
}
else
- digits_[size_++] = result;
+ push(result);
}
// one last round for the remaining decimal digits
if (first != last)
{
- word_type radix_power = 1U;
+ digit_type radix_power = 1U;
digit_type result = 0U;
while (first != last)
@@ -141,10 +142,10 @@
carry += ops_type::add_single_digit(digits_, digits_, size_, result);
if (carry)
- digits_[size_++] = carry;
+ push(carry);
}
else
- digits_[size_++] = result;
+ push(result);
}
}
}
@@ -174,7 +175,7 @@
StringT s;
- if (!size_)
+ if (is_uninitialized())
return s;
digit_type radix;
@@ -212,7 +213,7 @@
const int prefix_offset = p - prefix;
- if (is_zero())
+ if (!*this)
{
s.reserve(prefix_offset + 1);
for (int i = 0; i < prefix_offset; ++i)
@@ -258,7 +259,7 @@
? uppercase_tab
: lowercase_tab;
- const digit_type mask = (1 << sc.radix_storage_bits) - 1;
+ const digit_type mask = (digit_type(1) << sc.radix_storage_bits) - 1;
int offset = total_bits % valid_bits;
if (!offset)
@@ -288,14 +289,18 @@
mp_int tmp = abs(*this);
- while (!tmp.is_zero())
+ while (tmp)
{
- digit_type remainder = tmp.divide_by_digit(sc.max_power_value);
-
+ digit_type remainder = ops_type::divide_by_digit(tmp.digits(),
+ tmp.digits(),
+ tmp.size(),
+ sc.max_power_value);
+ tmp.clamp_high_digit();
+
for (digit_type i = 0; i < m; ++i)
{
- if (remainder || !tmp.is_zero())
- *c++ = '0' + remainder % 10U;
+ if (remainder || tmp)
+ *c++ = static_cast<char_type>('0' + remainder % 10U);
remainder /= 10U;
}
}
Added: sandbox/mp_math/libs/mp_math/test/bitmanipulation.cpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/libs/mp_math/test/bitmanipulation.cpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -0,0 +1,78 @@
+// Copyright Kevin Sopp 2009.
+// 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(set_bits1, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0xff00000000ff");
+ x.set_bits(8, 40);
+ BOOST_CHECK_EQUAL(x, "0xffffffffffff");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(set_bits2, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0x8000");
+ x.set_bits(2, 7);
+ BOOST_CHECK_EQUAL(x, "0x807c");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(set_bits3, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0x80000000000000");
+ x.set_bits(12, 13);
+ BOOST_CHECK_EQUAL(x, "0x80000000001000");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(set_bits4, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0x8000000000000000");
+ x.set_bits(0, 18);
+ BOOST_CHECK_EQUAL(x, "0x800000000003FFFF");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(set_bits5, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0x80000000");
+ x.set_bits(8, 16);
+ BOOST_CHECK_EQUAL(x, "0x8000FF00");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(clear_bits1, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0xffffffffffff");
+ x.clear_bits(8, 40);
+ BOOST_CHECK_EQUAL(x, "0xff00000000ff");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(clear_bits2, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0x807c");
+ x.clear_bits(2, 7);
+ BOOST_CHECK_EQUAL(x, "0x8000");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(clear_bits3, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0x80000000001000");
+ x.clear_bits(12, 13);
+ BOOST_CHECK_EQUAL(x, "0x80000000000000");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(clear_bits4, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0x800000000003FFFF");
+ x.clear_bits(0, 18);
+ BOOST_CHECK_EQUAL(x, "0x8000000000000000");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(clear_bits5, mp_int_type, mp_int_types)
+{
+ mp_int_type x("0x8000FF00");
+ x.clear_bits(8, 16);
+ BOOST_CHECK_EQUAL(x, "0x80000000");
+}
+
Modified: sandbox/mp_math/libs/mp_math/test/cmp.cpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/cmp.cpp (original)
+++ sandbox/mp_math/libs/mp_math/test/cmp.cpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -130,6 +130,12 @@
BOOST_CHECK_LE(x, 32102);
}
+BOOST_AUTO_TEST_CASE_TEMPLATE(cmp_mp_int_le_integral_type3, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("0x100000000");
+ BOOST_CHECK_LE(1, x);
+}
+
BOOST_AUTO_TEST_CASE_TEMPLATE(cmp_mp_int_le_unsigned_integral_type, mp_int_type, mp_int_types)
{
Modified: sandbox/mp_math/libs/mp_math/test/ctors.cpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/ctors.cpp (original)
+++ sandbox/mp_math/libs/mp_math/test/ctors.cpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -37,7 +37,7 @@
BOOST_AUTO_TEST_CASE_TEMPLATE(ctor_from_decimal_string4, mp_int_type, mp_int_types)
{
const mp_int_type y("0");
- BOOST_CHECK_EQUAL(y.is_zero(), true);
+ BOOST_CHECK_EQUAL(!y, true);
}
BOOST_AUTO_TEST_CASE_TEMPLATE(ctor_from_positive_decimal_string1, mp_int_type, mp_int_types)
@@ -67,7 +67,7 @@
BOOST_AUTO_TEST_CASE_TEMPLATE(ctor_from_octal_string2, mp_int_type, mp_int_types)
{
const mp_int_type y("000000000000000000000000000000000");
- BOOST_CHECK_EQUAL(y.is_zero(), true);
+ BOOST_CHECK_EQUAL(!y, true);
}
BOOST_AUTO_TEST_CASE_TEMPLATE(ctor_from_octal_string3, mp_int_type, mp_int_types)
Modified: sandbox/mp_math/libs/mp_math/test/div.cpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/div.cpp (original)
+++ sandbox/mp_math/libs/mp_math/test/div.cpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -30,6 +30,22 @@
BOOST_CHECK_EQUAL(z, "2137");
}
+BOOST_AUTO_TEST_CASE_TEMPLATE(divide4, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("0x89ab89745cc653de58eecc8f8a874120065ea545f6f5f");
+ const mp_int_type y("0x1889ba8a789456adfc8005b1");
+ const mp_int_type z = x / y;
+ BOOST_CHECK_EQUAL(z, "0x59c48aa7a1446110b31f70");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(divide5, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("0x1889ba8a789456adfc8005b1");
+ const mp_int_type y("0x1889ba8a789456adfc8005b2");
+ const mp_int_type z = x / y;
+ BOOST_CHECK_EQUAL(z, "0");
+}
+
BOOST_AUTO_TEST_CASE_TEMPLATE(divide_by_2, mp_int_type, mp_int_types)
{
mp_int_type x("263825472");
Modified: sandbox/mp_math/libs/mp_math/test/integral_ops.cpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/integral_ops.cpp (original)
+++ sandbox/mp_math/libs/mp_math/test/integral_ops.cpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -9,7 +9,7 @@
BOOST_AUTO_TEST_CASE_TEMPLATE(construct_from_zero, mp_int_type, mp_int_types)
{
const mp_int_type x(0);
- BOOST_CHECK_EQUAL(x.is_zero(), true);
+ BOOST_CHECK_EQUAL(!x, true);
}
BOOST_AUTO_TEST_CASE_TEMPLATE(equal_signed_char_min, mp_int_type, mp_int_types)
@@ -190,6 +190,21 @@
BOOST_CHECK_EQUAL(z, "0");
}
+BOOST_AUTO_TEST_CASE_TEMPLATE(divide_by_unsigned_char1, mp_int_type, mp_int_types)
+{
+ const unsigned char y = 16;
+ const mp_int_type x("10000001");
+ const mp_int_type z = x / y;
+ BOOST_CHECK_EQUAL(z, "625000");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(divide_by_unsigned_char2, mp_int_type, mp_int_types)
+{
+ const unsigned char y = 128;
+ const mp_int_type x("14222200");
+ const mp_int_type z = x / y;
+ BOOST_CHECK_EQUAL(z, "111110");
+}
BOOST_AUTO_TEST_CASE_TEMPLATE(divide_by_signed_integral1, mp_int_type, mp_int_types)
{
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 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -1,4 +1,4 @@
-# Copyright Kevin Sopp 2008.
+# Copyright Kevin Sopp 2008 - 2009.
# 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)
@@ -7,18 +7,19 @@
: $(BOOST_ROOT)/libs/test/build//boost_unit_test_framework/<link>shared ;
alias boost_serialization
-: $(BOOST_ROOT)/libs/serialization/build//boost_serialization/<link>shared ;
+: $(BOOST_ROOT)/libs/serialization/build//boost_serialization/<link>shared ;
project
- : requirements
- <include>../../..
- <link>static
- <define>BOOST_TEST_DYN_LINK
- <define>BOOST_TEST_MAIN
- #<define>BOOST_MP_MATH_MP_INT_USE_ASM
- ;
+ : requirements
+ <include>../../..
+ <link>static
+ <define>BOOST_TEST_DYN_LINK
+ <define>BOOST_TEST_MAIN
+ <define>BOOST_MP_MATH_MP_INT_USE_ASM
+ ;
unit-test add : add.cpp boost_test ;
+unit-test bitmanipulation : bitmanipulation.cpp boost_test ;
unit-test bitwise_ops : bitwise_ops.cpp boost_test ;
unit-test bool_conversion : bool_conversion.cpp boost_test ;
unit-test cmp : cmp.cpp boost_test ;
Modified: sandbox/mp_math/libs/mp_math/test/root.cpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/root.cpp (original)
+++ sandbox/mp_math/libs/mp_math/test/root.cpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -27,3 +27,17 @@
BOOST_CHECK_EQUAL(y, "441");
}
+BOOST_AUTO_TEST_CASE_TEMPLATE(nth_root2, mp_int_type, mp_int_types)
+{
+ const mp_int_type x(
+ "0x2b93d251afa09c5481f4522279f7c19ca08124199621dfd18342a16c7303b31ccea8176b"
+ "d4a7a9bf991e30d8bde1e08356a728b9f5729c35d29884050101341228c5df3f98354d42b7"
+ "a0d7fdfbe8d5270b09ee89ba1eeab61be67eb4471d92fdffa88d1ca494ed3eec58a34ff958"
+ "b518a588584a2505c9c2b19ce1eb21cba36c7a5297cb6e532884e89451f4406b993582f3cd"
+ "b75cab98f8c4c6f3837977db2a594dfa16943062187ca95babc9da78bdd73ca7233eefc047"
+ "8d882e0d4f09a5083a31b801964343d47b6ce9e937df8c44a9a02bac5101da1823373e663c"
+ "1329ece1eb89fc178355660fe1c92c7d8ff11524702fad6e2255447946442356b00810101");
+ const mp_int_type y = nth_root(x, mp_int_type("257"));
+ BOOST_CHECK_EQUAL(y, "257");
+}
+
Modified: sandbox/mp_math/libs/mp_math/test/shift.cpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/shift.cpp (original)
+++ sandbox/mp_math/libs/mp_math/test/shift.cpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -40,3 +40,10 @@
x >>= 17;
BOOST_CHECK_EQUAL(x, "0");
}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(right_shift3, mp_int_type, mp_int_types)
+{
+ mp_int_type x("14222200");
+ x >>= 8;
+ BOOST_CHECK_EQUAL(x, "55555");
+}
Modified: sandbox/mp_math/libs/mp_math/test/string_ops.cpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/string_ops.cpp (original)
+++ sandbox/mp_math/libs/mp_math/test/string_ops.cpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -119,7 +119,7 @@
{
mp_int_type x;
x = "0";
- BOOST_CHECK_EQUAL(x.is_zero(), true);
+ BOOST_CHECK_EQUAL(!x, true);
}
BOOST_AUTO_TEST_CASE_TEMPLATE(op_assign5, mp_int_type, mp_int_types)
Modified: sandbox/mp_math/libs/mp_math/tools/benchmark/benchmark_gmp.hpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/tools/benchmark/benchmark_gmp.hpp (original)
+++ sandbox/mp_math/libs/mp_math/tools/benchmark/benchmark_gmp.hpp 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -84,7 +84,7 @@
explicit square_op(base& ba) : b(ba) {}
void operator()(unsigned int i) const
{
- mpz_pow_ui(b.dst[i].get_mpz_t(), b.src1[i].get_mpz_t(), 2);
+ b.dst[i] = b.src1[i] * b.src1[i];
}
};
Modified: sandbox/mp_math/project-root.jam
==============================================================================
--- sandbox/mp_math/project-root.jam (original)
+++ sandbox/mp_math/project-root.jam 2009-01-10 19:30:26 EST (Sat, 10 Jan 2009)
@@ -10,13 +10,11 @@
## IMPORTANT NOTE: This file MUST NOT be copied over a boost installation
##
-path-constant top : . ;
import modules ;
import path ;
local boost-root = [ modules.peek : BOOST_ROOT ] ;
-local math-header-include = $(top)/../.. ;
if ! $(boost-root)
{