Subject: [Boost-commit] svn:boost r49448 - in sandbox/mp_math: boost/mp_math/mp_int libs/mp_math/test
From: baraclese_at_[hidden]
Date: 2008-10-24 20:23:09
Author: baraclese
Date: 2008-10-24 20:23:08 EDT (Fri, 24 Oct 2008)
New Revision: 49448
URL: http://svn.boost.org/trac/boost/changeset/49448
Log:
- fixed lcm, it was completely broken
- added tests for lcm
- added variadic gcd and lcm including tests
- minor optimization to gcd
- modinv
- renamed mp_int::slow_modinv to even_modinv
- renamed mp_int::fast_modinv to odd_modinv
- minor tweaks to both
- factored modinv functionality out into modinv.hpp
- added non-member modinv function, added tests
- some tweaks to mp_int::modulo_2_to_the_power_of
Added:
sandbox/mp_math/boost/mp_math/mp_int/modinv.hpp (contents, props changed)
sandbox/mp_math/libs/mp_math/test/lcm.cpp (contents, props changed)
sandbox/mp_math/libs/mp_math/test/modinv.cpp (contents, props changed)
Text files modified:
sandbox/mp_math/boost/mp_math/mp_int/gcd.hpp | 52 ++++-----
sandbox/mp_math/boost/mp_math/mp_int/lcm.hpp | 30 +++-
sandbox/mp_math/boost/mp_math/mp_int/mod.hpp | 209 ++-------------------------------------
sandbox/mp_math/boost/mp_math/mp_int/mp_int.hpp | 9
sandbox/mp_math/libs/mp_math/test/gcd.cpp | 30 +++++
sandbox/mp_math/libs/mp_math/test/jamfile.v2 | 3
6 files changed, 94 insertions(+), 239 deletions(-)
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-24 20:23:08 EDT (Fri, 24 Oct 2008)
@@ -3,58 +3,52 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
-/* Greatest Common Divisor using the binary method */
+// 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)
{
- typedef typename mp_int<A,T>::size_type size_type;
-
- /* either zero then gcd is the largest */
- if (a.is_zero())
+ // either zero then gcd is the largest
+ if (!a)
return abs(b);
- if (b.is_zero())
+ if (!b)
return abs(a);
- /* get copies of a and b we can modify */
+ // get copies of a and b we can modify
mp_int<A,T> u = abs(a);
mp_int<A,T> v = abs(b);
- /* B1. Find the common power of two for u and v */
+ typedef typename mp_int<A,T>::size_type size_type;
+
+ // Find the common power of two for u and v
const size_type u_lsb = u.count_lsb();
const size_type v_lsb = v.count_lsb();
const size_type k = std::min(u_lsb, v_lsb);
- if (k > 0)
- {
- /* divide the power of two out */
- u.shift_right(k,0);
- v.shift_right(k,0);
- }
-
- /* divide any remaining factors of two out */
- if (u_lsb != k)
- u.shift_right(u_lsb - k, 0);
-
- if (v_lsb != k)
- v.shift_right(v_lsb - k, 0);
+ // divide out powers of two
+ u >>= u_lsb;
+ v >>= v_lsb;
- while (!v.is_zero())
+ while (v)
{
- /* make sure v is the largest */
- if (u.compare_magnitude(v) == 1)
- /* swap u and v to make sure v is >= u */
+ if (u > v)
u.swap(v);
- /* subtract smallest from largest */
v.sub_smaller_magnitude(u);
- /* Divide out all factors of two */
- v.shift_right(v.count_lsb(), 0);
+ // Divide out all factors of two
+ v >>= v.count_lsb();
}
- /* multiply by 2**k which we divided out at the beginning */
+ // multiply by 2**k which we divided out at the beginning
u <<= k;
return u;
}
+#ifdef BOOST_HAS_VARIADIC_TMPL
+template<class A, class T, class... MpInts>
+mp_int<A,T> gcd(const mp_int<A,T>& a, const mp_int<A,T>& b, const MpInts&... args)
+{
+ return gcd(gcd(a, b), args...);
+}
+#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-24 20:23:08 EDT (Fri, 24 Oct 2008)
@@ -3,20 +3,30 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
-/* computes least common multiple as |a*b|/(a, b) */
+// 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)
{
- /* t1 = get the GCD of the two inputs */
- const mp_int<A,T> t1 = gcd(a,b);
-
- /* divide the smallest by the GCD */
- const mp_int<A,T>* smallest = a.compare_magnitude(b) == -1 ? &a : &b;
-
- mp_int<A,T> t2 = *smallest / t1;
+ mp_int<A,T> result;
+
+ if (!a || !b)
+ {
+ result.zero();
+ return result;
+ }
- t2.sign_ = 1;
+ result = a / gcd(a, b) * b;
+
+ result.set_sign(1);
- return t2;
+ return result;
+}
+
+#ifdef BOOST_HAS_VARIADIC_TMPL
+template<class A, class T, class... MpInts>
+mp_int<A,T> lcm(const mp_int<A,T>& a, const mp_int<A,T>& b, const MpInts&... args)
+{
+ return lcm(lcm(a, b), args...);
}
+#endif
Modified: sandbox/mp_math/boost/mp_math/mp_int/mod.hpp
==============================================================================
--- sandbox/mp_math/boost/mp_math/mp_int/mod.hpp (original)
+++ sandbox/mp_math/boost/mp_math/mp_int/mod.hpp 2008-10-24 20:23:08 EDT (Fri, 24 Oct 2008)
@@ -3,207 +3,24 @@
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
-/* calc a value mod 2**b */
+// *this % 2**b
template<class A, class T>
-void mp_int<A,T>::modulo_2_to_the_power_of(int b)
+void mp_int<A,T>::modulo_2_to_the_power_of(size_type b)
{
- /* if b is <= 0 then zero the int */
- if (b <= 0)
- {
- zero();
+ // if modulus >= *this then return
+ if (b >= used_ * valid_bits)
return;
- }
-
- /* if the modulus is larger than the value then return */
- if (b >= static_cast<int>(used_ * valid_bits))
- return;
-
- /* zero digits above the last digit of the modulus */
- for (size_type x = (b / valid_bits) + ((b % valid_bits) == 0 ? 0 : 1); x < used_; ++x)
- digits_[x] = 0;
-
- /* clear the digit that is not completely outside/inside the modulus */
- digits_[b / valid_bits] &= static_cast<digit_type>(
- ((digit_type(1)) << ((static_cast<digit_type>(b)) % valid_bits)) - (digit_type(1)));
+ // zero digits above the last digit of the modulus
+ const size_type offset = (b / valid_bits) + ((b % valid_bits) == 0 ? 0 : 1);
+ std::memset(digits_ + offset, 0, sizeof(digit_type) * (used_ - offset));
+
+ // clear remaining high bits
+ const digit_type mask = (1 << (static_cast<digit_type>(b % valid_bits))) - 1;
+ digits_[b / valid_bits] &= mask;
+
clamp();
if (is_zero())
sign_ = 1;
-}
-
-// hac 14.61, pp608
-// *this = *this**-1 (mod b)
-template<class A, class T>
-void mp_int<A,T>::modinv(const mp_int& b)
-{
- /* b cannot be negative */
- if (b.is_negative() || b.is_zero())
- throw std::domain_error("modinv: b is negative or zero");
-
- /* if the modulus is odd we can use a faster routine instead */
- if (b.is_odd())
- fast_modinv(b);
- else
- slow_modinv(b);
-}
-
-/* hac 14.61, pp608 */
-template<class Al, class T>
-void mp_int<Al,T>::slow_modinv(const mp_int& b)
-{
- /* b cannot be negative */
- if (b.is_negative() || b.is_zero())
- throw std::domain_error("mp_int::slow_modinv: b is negative or zero"); // XXX: this is already tested in modinv
-
- const mp_int x = *this % b;
- const mp_int y(b); // TODO no need to copy b here since b or y is never changed
-
- /* [modified] if x,y are both even then return an error! */
- if (x.is_even() && y.is_even())
- throw std::domain_error("mp_int::slow_modinv: no inverse exists");// TODO: different text?
-
- 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.is_zero())
- 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("mp_int::slow_modinv: no inverse exists"); // TODO return false if no inverse exists?
-
- /* if its too low */
- while (C.compare_to_digit(0) == -1)
- C += b;
-
- /* too big */
- while (C.compare_magnitude(b) != -1)
- C -= b;
-
- /* C is now the inverse */
- swap(C);
-}
-
-/* computes the modular inverse via binary extended euclidean algorithm,
- * that is *this = 1 / *this mod b
- *
- * Based on slow modinv except this is optimized for the case where b is
- * odd as per HAC Note 14.64 on pp. 610
- */
-template<class Al, class T>
-void mp_int<Al,T>::fast_modinv(const mp_int& b)
-{
- if (b.is_even())
- throw std::domain_error("mp_int::fast_modinv: b must be odd");
-
- /* x == modulus, y == value to invert */
- mp_int x = b;
-
- /* we need y = |a| */
- mp_int y = *this % b;
-
- /* 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.is_zero())
- 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("mp_int::fast_modinv: no inverse exists");
-
- /* D is now the inverse */
- while (D.sign_ == -1)
- D += b;
-
- swap(D);
-}
+ }
Added: sandbox/mp_math/boost/mp_math/mp_int/modinv.hpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/boost/mp_math/mp_int/modinv.hpp 2008-10-24 20:23:08 EDT (Fri, 24 Oct 2008)
@@ -0,0 +1,187 @@
+// 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)
+
+
+// hac 14.61, pp608
+// *this = *this**-1 (mod b)
+template<class A, class T>
+void mp_int<A,T>::modinv(const mp_int& b)
+{
+ if (b.is_negative() || !b)
+ throw std::domain_error("modinv: modulus is negative or zero");
+
+ // if the modulus is odd we can use a faster routine
+ if (b.is_odd())
+ odd_modinv(b);
+ else
+ even_modinv(b);
+}
+
+/* hac 14.61, pp608 */
+template<class A1, class T>
+void mp_int<A1,T>::even_modinv(const mp_int& y)
+{
+ assert(y.is_positive() && y);
+
+ 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;
+
+ swap(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;
+
+ swap(D);
+}
+
+
+// returns the modular multiplicative inverse x of a (mod m) such that
+// x*a = 1 (mod m) =>
+// a^-1 = x (mod m)
+// The inverse exists only if a and m are coprime (i.e. gcd(a,m) = 1).
+// If no inverse exists this function will throw std::domain_error.
+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;
+}
+
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-24 20:23:08 EDT (Fri, 24 Oct 2008)
@@ -312,7 +312,7 @@
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(int);
+ void modulo_2_to_the_power_of(size_type);
size_type precision() const;
size_type count_lsb() const;
void shift_right(size_type b, mp_int* remainder);
@@ -344,9 +344,9 @@
void barret_modpow(const mp_int& exp, const mp_int& m, int reduction_mode);
void fast_modpow(const mp_int& exp, const mp_int& m, int reduction_mode);
- void modinv(const mp_int& b);
- void slow_modinv(const mp_int& b);
- void fast_modinv(const mp_int& b);
+ 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()
{
@@ -873,6 +873,7 @@
#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/mul.hpp>
#include <boost/mp_math/mp_int/operators.hpp>
Modified: sandbox/mp_math/libs/mp_math/test/gcd.cpp
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/gcd.cpp (original)
+++ sandbox/mp_math/libs/mp_math/test/gcd.cpp 2008-10-24 20:23:08 EDT (Fri, 24 Oct 2008)
@@ -30,3 +30,33 @@
const mp_int_type z = boost::mp_math::gcd(x,y);
BOOST_CHECK_EQUAL(z, "16384");
}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(gcd4, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("0");
+ const mp_int_type y("0");
+ const mp_int_type z = boost::mp_math::gcd(x,y);
+ BOOST_CHECK_EQUAL(z, "0");
+}
+
+#ifdef BOOST_HAS_VARIADIC_TMPL
+BOOST_AUTO_TEST_CASE_TEMPLATE(variadic_gcd1, mp_int_type, mp_int_types)
+{
+ const mp_int_type a("42");
+ const mp_int_type b("56");
+ const mp_int_type c("140");
+ const mp_int_type z = boost::mp_math::gcd(a,b,c);
+ BOOST_CHECK_EQUAL(z, "14");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(variadic_gcd2, mp_int_type, mp_int_types)
+{
+ const mp_int_type a("1200000000");
+ const mp_int_type b("2400000000");
+ const mp_int_type c("3600000000");
+ const mp_int_type d("600000000000000");
+ const mp_int_type z = boost::mp_math::gcd(a,b,c,d);
+ BOOST_CHECK_EQUAL(z, "1200000000");
+}
+#endif
+
Modified: sandbox/mp_math/libs/mp_math/test/jamfile.v2
==============================================================================
--- sandbox/mp_math/libs/mp_math/test/jamfile.v2 (original)
+++ sandbox/mp_math/libs/mp_math/test/jamfile.v2 2008-10-24 20:23:08 EDT (Fri, 24 Oct 2008)
@@ -13,6 +13,7 @@
: requirements
<include>../../..
<link>static
+ <toolset>gcc:<cxxflags>-std=c++0x
<define>BOOST_TEST_DYN_LINK
<define>BOOST_TEST_MAIN
#<define>BOOST_MP_MATH_MP_INT_USE_ASM
@@ -27,6 +28,8 @@
unit-test gcd : gcd.cpp boost_test ;
unit-test integral_ops : integral_ops.cpp boost_test ;
unit-test jacobi : jacobi.cpp boost_test ;
+unit-test lcm : lcm.cpp boost_test ;
+unit-test modinv : modinv.cpp boost_test ;
unit-test mul : mul.cpp boost_test ;
unit-test pow : pow.cpp boost_test ;
unit-test prime : prime.cpp boost_test ;
Added: sandbox/mp_math/libs/mp_math/test/lcm.cpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/libs/mp_math/test/lcm.cpp 2008-10-24 20:23:08 EDT (Fri, 24 Oct 2008)
@@ -0,0 +1,54 @@
+// Copyright Kevin Sopp 2008.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/test/unit_test.hpp>
+#include "prerequisite.hpp"
+
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(lcm1, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("0");
+ const mp_int_type y("0");
+ const mp_int_type z = boost::mp_math::lcm(x,y);
+ BOOST_CHECK_EQUAL(z, "0");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(lcm2, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("51111");
+ const mp_int_type y("0");
+ const mp_int_type z = boost::mp_math::lcm(x,y);
+ BOOST_CHECK_EQUAL(z, "0");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(lcm3, mp_int_type, mp_int_types)
+{
+ const mp_int_type x("4");
+ const mp_int_type y("6");
+ const mp_int_type z = boost::mp_math::lcm(x,y);
+ BOOST_CHECK_EQUAL(z, "12");
+}
+
+#ifdef BOOST_HAS_VARIADIC_TMPL
+BOOST_AUTO_TEST_CASE_TEMPLATE(variadic_lcm1, mp_int_type, mp_int_types)
+{
+ const mp_int_type a("120");
+ const mp_int_type b("204");
+ const mp_int_type c("136");
+ const mp_int_type z = boost::mp_math::lcm(a,b,c);
+ BOOST_CHECK_EQUAL(z, "2040");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(variadic_lcm2, mp_int_type, mp_int_types)
+{
+ const mp_int_type a("12010");
+ const mp_int_type b("3299");
+ const mp_int_type c("84780");
+ const mp_int_type d("15");
+ const mp_int_type z = boost::mp_math::lcm(a,b,c,d);
+ BOOST_CHECK_EQUAL(z, "335906753220");
+}
+#endif
+
Added: sandbox/mp_math/libs/mp_math/test/modinv.cpp
==============================================================================
--- (empty file)
+++ sandbox/mp_math/libs/mp_math/test/modinv.cpp 2008-10-24 20:23:08 EDT (Fri, 24 Oct 2008)
@@ -0,0 +1,25 @@
+// Copyright Kevin Sopp 2008.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/test/unit_test.hpp>
+#include "prerequisite.hpp"
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modinv1, mp_int_type, mp_int_types)
+{
+ mp_int_type a("35");
+ mp_int_type m("33");
+ mp_int_type i = boost::mp_math::modinv(a, m);
+ BOOST_CHECK_EQUAL(i, "17");
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(modinv2, mp_int_type, mp_int_types)
+{
+ mp_int_type a("17");
+ mp_int_type m("26");
+ mp_int_type i = boost::mp_math::modinv(a, m);
+ BOOST_CHECK_EQUAL(i, "23");
+}
+
+