[Boost-commit] svn:boost r76398 - in sandbox/big_number: boost/multiprecision libs/multiprecision/test

Subject: [Boost-commit] svn:boost r76398 - in sandbox/big_number: boost/multiprecision libs/multiprecision/test
From: john_at_[hidden]
Date: 2012-01-10 08:20:31

Author: johnmaddock
Date: 2012-01-10 08:20:29 EST (Tue, 10 Jan 2012)
New Revision: 76398
URL: http://svn.boost.org/trac/boost/changeset/76398

Log:
Switch to faster division code and add better test case for packed int's.
Added:
   sandbox/big_number/libs/multiprecision/test/packed_int_test.cpp   (contents, props changed)
Text files modified: 
   sandbox/big_number/boost/multiprecision/packed_cpp_int.hpp |   166 +++++++++++++++++++++++++++++++++++++++ 
   sandbox/big_number/libs/multiprecision/test/Jamfile.v2     |     8 +                                       
   2 files changed, 173 insertions(+), 1 deletions(-)

Modified: sandbox/big_number/boost/multiprecision/packed_cpp_int.hpp
==============================================================================
--- sandbox/big_number/boost/multiprecision/packed_cpp_int.hpp	(original)
+++ sandbox/big_number/boost/multiprecision/packed_cpp_int.hpp	2012-01-10 08:20:29 EST (Tue, 10 Jan 2012)
@@ -504,7 +504,7 @@
    }
    return count;
 }
-
+/*
 template <unsigned Bits, bool Signed>
 void divide_unsigned_helper(packed_cpp_int<Bits, Signed>& result, const packed_cpp_int<Bits, Signed>& x, const packed_cpp_int<Bits, Signed>& y, packed_cpp_int<Bits, Signed>& r)
 {
@@ -588,6 +588,170 @@
    }
    BOOST_ASSERT(r.compare(y) < 0); // remainder must be less than the divisor or our code has failed
 }
+*/
+template <unsigned Bits, bool Signed>
+void divide_unsigned_helper(packed_cpp_int<Bits, Signed>& result, const packed_cpp_int<Bits, Signed>& x, const packed_cpp_int<Bits, Signed>& y, packed_cpp_int<Bits, Signed>& r)
+{
+   if((&result == &x) || (&r == &x))
+   {
+      packed_cpp_int<Bits, Signed> t(x);
+      divide_unsigned_helper(result, t, y, r);
+      return;
+   }
+   if((&result == &y) || (&r == &y))
+   {
+      packed_cpp_int<Bits, Signed> t(y);
+      divide_unsigned_helper(result, x, t, r);
+      return;
+   }
+
+
+   using default_ops::subtract;
+
+
+   if (is_zero(y))
+   {
+      volatile int i = 0;
+      i /= i;
+      return;
+   }
+
+   if (is_zero(x))
+   {
+      r = y;
+      result = x;
+      return;
+   }
+
+   if(&result == &r)
+   {
+      packed_cpp_int<Bits, Signed> rem;
+      divide_unsigned_helper(result, x, y, rem);
+      r = rem;
+      return;
+   }
+
+   r = x;
+   result = static_cast<boost::uint32_t>(0u);
+   if(x.compare(y) < 0)
+   {
+      return; // We already have the answer: zero.
+   }
+
+   //
+   // Find the most significant words of numerator and denominator.
+   // Note that this code can't run past the end of the array because
+   // we know already that neither are all zero:
+   //
+   boost::uint32_t r_order = 0;
+   while(r.data()[r_order] == 0)
+      ++r_order;
+   boost::uint32_t y_order = 0;
+   while(y.data()[y_order] == 0)
+      ++y_order;
+
+   packed_cpp_int<Bits, Signed> t;
+   bool r_neg = false;
+
+   //
+   // See if we can short-circuit long division, and use basic arithmetic instead:
+   //
+   if(r_order == packed_cpp_int<Bits, Signed>::limb_count - 1)
+   {
+      result = r.data()[packed_cpp_int<Bits, Signed>::limb_count - 1] / y.data()[packed_cpp_int<Bits, Signed>::limb_count - 1];
+      r = x.data()[packed_cpp_int<Bits, Signed>::limb_count - 1] % y.data()[packed_cpp_int<Bits, Signed>::limb_count - 1];
+      return;
+   }
+   else if(r_order == packed_cpp_int<Bits, Signed>::limb_count - 2)
+   {
+      unsigned long long a, b;
+      a = (static_cast<unsigned long long>(r.data()[r_order]) << packed_cpp_int<Bits, Signed>::limb_bits) | r.data()[r_order + 1];
+      b = y_order < packed_cpp_int<Bits, Signed>::limb_count - 1 ? 
+         (static_cast<unsigned long long>(y.data()[y_order]) << packed_cpp_int<Bits, Signed>::limb_bits) | y.data()[y_order + 1] 
+         : y.data()[y_order];
+      result = a / b;
+      r = a % b;
+      return;
+   }
+
+   do
+   {
+      //
+      // Update r_order, this can't run past the end as r must be non-zero at this point:
+      //
+      while(r.data()[r_order] == 0)
+         ++r_order;
+      //
+      // Calculate our best guess for how many times y divides into r:
+      //
+      boost::uint32_t guess;
+      if((r.data()[r_order] <= y.data()[y_order]) && (r_order < packed_cpp_int<Bits, Signed>::limb_count - 1))
+      {
+         unsigned long long a, b, v;
+         a = (static_cast<unsigned long long>(r.data()[r_order]) << packed_cpp_int<Bits, Signed>::limb_bits) | r.data()[r_order + 1];
+         b = y.data()[y_order];
+         v = a / b;
+         if(v > packed_cpp_int<Bits, Signed>::max_limb_value)
+            guess = 1;
+         else
+         {
+            guess = static_cast<boost::uint32_t>(v);
+            ++r_order;
+         }
+      }
+      else
+      {
+         guess = r.data()[r_order] / y.data()[y_order];
+      }
+      //
+      // Update result:
+      //
+      boost::uint32_t shift = y_order - r_order;
+      t = boost::uint32_t(0);
+      t.data()[packed_cpp_int<Bits, Signed>::limb_count - 1 - shift] = guess;
+      if(r_neg)
+         subtract(result, t);
+      else
+         add(result, t);
+      //
+      // Calculate guess * y, we use a fused mutiply-shift O(N) for this
+      // rather than a full O(N^2) multiply:
+      //
+      boost::uintmax_t carry = 0;
+      for(unsigned i = packed_cpp_int<Bits, Signed>::limb_count - 1; i > packed_cpp_int<Bits, Signed>::limb_count - shift - 1; --i)
+         t.data()[i] = 0;
+      for(int i = packed_cpp_int<Bits, Signed>::limb_count - 1; i >= static_cast<int>(shift); --i)
+      {
+         boost::uintmax_t v = static_cast<boost::uintmax_t>(y.data()[i]) * static_cast<boost::uintmax_t>(guess) + carry;
+         t.data()[i - shift] = static_cast<packed_cpp_int<Bits, Signed>::limb_type>(v);
+         carry = v >> sizeof(packed_cpp_int<Bits, Signed>::limb_type) * CHAR_BIT;
+      }
+      t.data()[0] &= packed_cpp_int<Bits, Signed>::upper_limb_mask;
+       //
+      // Update r:
+      //
+      subtract(r, t);
+      if(r.data()[0] & packed_cpp_int<Bits, Signed>::sign_bit_mask)
+      {
+         r.negate();
+         r_neg = !r_neg;
+      }
+   }
+   while(r.compare(y) > 0);
+
+   //
+   // We now just have to normalise the result:
+   //
+   if(r_neg)
+   {
+      // We have one too many in the result:
+      decrement(result);
+      r.negate();
+      add(r, y);
+   }
+
+   BOOST_ASSERT(r.compare(y) < 0); // remainder must be less than the divisor or our code has failed
+}
 
 template <unsigned Bits, bool Signed>
 inline void divide(packed_cpp_int<Bits, Signed>& result, const packed_cpp_int<Bits, Signed>& a, const packed_cpp_int<Bits, Signed>& b)

Modified: sandbox/big_number/libs/multiprecision/test/Jamfile.v2
==============================================================================
--- sandbox/big_number/libs/multiprecision/test/Jamfile.v2	(original)
+++ sandbox/big_number/libs/multiprecision/test/Jamfile.v2	2012-01-10 08:20:29 EST (Tue, 10 Jan 2012)
@@ -595,6 +595,14 @@
               <define>TEST_PACKED_INT2
         : test_int_io_packed_int2 ;
 
+run packed_int_test.cpp gmp
+        : # command line
+        : # input files
+        : # requirements
+         [ check-target-builds ../config//has_gmp : : <build>no ]
+         release  # otherwise runtime is too slow!!
+        : packed_int_test.cpp ;
+
 run test_rational_io.cpp $(TOMMATH)
         : # command line
         : # input files

Added: sandbox/big_number/libs/multiprecision/test/packed_int_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/big_number/libs/multiprecision/test/packed_int_test.cpp	2012-01-10 08:20:29 EST (Tue, 10 Jan 2012)
@@ -0,0 +1,151 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2012 John Maddock. 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_
+
+//
+// Compare arithmetic results using packed_cpp_int to GMP results.
+//
+
+#ifdef _MSC_VER
+#  define _SCL_SECURE_NO_WARNINGS
+#endif
+
+#define BOOST_CHRONO_HEADER_ONLY
+
+#include <boost/multiprecision/gmp.hpp>
+#include <boost/multiprecision/packed_cpp_int.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/uniform_int.hpp>
+#include <boost/chrono.hpp>
+#include "test.hpp"
+
+template <class Clock>
+struct stopwatch
+{
+   typedef typename Clock::duration duration;
+   stopwatch()
+   {
+      m_start = Clock::now();
+   }
+   duration elapsed()
+   {
+      return Clock::now() - m_start;
+   }
+   void reset()
+   {
+      m_start = Clock::now();
+   }
+
+private:
+   typename Clock::time_point m_start;
+};
+
+template <class T>
+T generate_random(unsigned bits_wanted)
+{
+   static boost::random::mt19937 gen;
+   typedef boost::random::mt19937::result_type random_type;
+
+   T max_val;
+   unsigned digits;
+   if(std::numeric_limits<T>::is_bounded && (bits_wanted == std::numeric_limits<T>::digits))
+   {
+      max_val = (std::numeric_limits<T>::max)();
+      digits = std::numeric_limits<T>::digits;
+   }
+   else
+   {
+      max_val = T(1) << bits_wanted;
+      digits = bits_wanted;
+   }
+
+   unsigned bits_per_r_val = std::numeric_limits<random_type>::digits - 1;
+   while((random_type(1) << bits_per_r_val) > (gen.max)()) --bits_per_r_val;
+
+   unsigned terms_needed = digits / bits_per_r_val + 1;
+
+   T val = 0;
+   for(unsigned i = 0; i < terms_needed; ++i)
+   {
+      val *= (gen.max)();
+      val += gen();
+   }
+   val %= max_val;
+   return val;
+}
+
+int main()
+{
+   using namespace boost::multiprecision;
+   typedef mp_number<packed_cpp_int<1024, true> > packed_type;
+   stopwatch<boost::chrono::high_resolution_clock> w;
+   unsigned last_error_count = 0;
+   for(int i = 0; i < 1000; ++i)
+   {
+      mpz_int a = generate_random<mpz_int>(1000);
+      mpz_int b = generate_random<mpz_int>(512);
+      mpz_int c = generate_random<mpz_int>(256);
+      mpz_int d = generate_random<mpz_int>(32);
+
+      packed_type a1 = a.str();
+      packed_type b1 = b.str();
+      packed_type c1 = c.str();
+      packed_type d1 = d.str();
+
+      BOOST_CHECK_EQUAL(a.str(), a1.str());
+      BOOST_CHECK_EQUAL(b.str(), b1.str());
+      BOOST_CHECK_EQUAL(c.str(), c1.str());
+      BOOST_CHECK_EQUAL(d.str(), d1.str());
+      BOOST_CHECK_EQUAL(mpz_int(a+b).str(), packed_type(a1 + b1).str());
+      BOOST_CHECK_EQUAL(mpz_int(a-b).str(), packed_type(a1 - b1).str());
+      BOOST_CHECK_EQUAL(mpz_int(mpz_int(-a)+b).str(), packed_type(packed_type(-a1) + b1).str());
+      BOOST_CHECK_EQUAL(mpz_int(mpz_int(-a)-b).str(), packed_type(packed_type(-a1) - b1).str());
+      BOOST_CHECK_EQUAL(mpz_int(c * d).str(), packed_type(c1 * d1).str());
+      BOOST_CHECK_EQUAL(mpz_int(b * c).str(), packed_type(b1 * c1).str());
+      BOOST_CHECK_EQUAL(mpz_int(a / b).str(), packed_type(a1 / b1).str());
+      BOOST_CHECK_EQUAL(mpz_int(a / d).str(), packed_type(a1 / d1).str());
+      BOOST_CHECK_EQUAL(mpz_int(a % b).str(), packed_type(a1 % b1).str());
+      BOOST_CHECK_EQUAL(mpz_int(a % d).str(), packed_type(a1 % d1).str());
+
+      if(last_error_count != boost::detail::test_errors())
+      {
+         last_error_count = boost::detail::test_errors();
+         std::cout << std::hex << std::showbase;
+
+         std::cout << "a    = " << a << std::endl;
+         std::cout << "a1   = " << a1 << std::endl;
+         std::cout << "b    = " << b << std::endl;
+         std::cout << "b1   = " << b1 << std::endl;
+         std::cout << "c    = " << c << std::endl;
+         std::cout << "c1   = " << c1 << std::endl;
+         std::cout << "d    = " << d << std::endl;
+         std::cout << "d1   = " << d1 << std::endl;
+         std::cout << "a + b   = " << a+b << std::endl;
+         std::cout << "a1 + b1 = " << a1+b1 << std::endl;
+         std::cout << "a - b   = " << a-b << std::endl;
+         std::cout << "a1 - b1 = " << a1-b1 << std::endl;
+         std::cout << "-a + b   = " << mpz_int(-a)+b << std::endl;
+         std::cout << "-a1 + b1 = " << packed_type(-a1)+b1 << std::endl;
+         std::cout << "-a - b   = " << mpz_int(-a)-b << std::endl;
+         std::cout << "-a1 - b1 = " << packed_type(-a1)-b1 << std::endl;
+         std::cout << "c*d    = " << c*d << std::endl;
+         std::cout << "c1*d1  = " << c1*d1 << std::endl;
+         std::cout << "b*c    = " << b*c << std::endl;
+         std::cout << "b1*c1  = " << b1*c1 << std::endl;
+         std::cout << "a/b    = " << a/b << std::endl;
+         std::cout << "a1/b1  = " << a1/b1 << std::endl;
+         std::cout << "a/d    = " << a/d << std::endl;
+         std::cout << "a1/d1  = " << a1/d1 << std::endl;
+         std::cout << "a%b    = " << a%b << std::endl;
+         std::cout << "a1%b1  = " << a1%b1 << std::endl;
+         std::cout << "a%d    = " << a%d << std::endl;
+         std::cout << "a1%d1  = " << a1%d1 << std::endl;
+      }
+   }
+   std::cout << w.elapsed() << std::endl;
+   return boost::report_errors();
+}
+
+
+