[Boost-commit] svn:boost r85416 - in sandbox/multiprecision.cpp_bin_float: boost/multiprecision libs/multiprecision/test

Subject: [Boost-commit] svn:boost r85416 - in sandbox/multiprecision.cpp_bin_float: boost/multiprecision libs/multiprecision/test
From: john_at_[hidden]
Date: 2013-08-21 12:56:33

Author: johnmaddock
Date: 2013-08-21 12:56:33 EDT (Wed, 21 Aug 2013)
New Revision: 85416
URL: http://svn.boost.org/trac/boost/changeset/85416

Log:
Initial commit of correctly rounded binary-decimal and decimal-binary conversion routines.

Added:
   sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/test_cpp_bin_float_io.cpp
      - copied, changed from r85253, trunk/libs/multiprecision/test/test_float_io.cpp
Text files modified: 
   sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float.hpp             |   293 ++++++++++++++++++++++++++++++++++++++- 
   sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/test_cpp_bin_float_io.cpp |   161 ++-------------------                   
   2 files changed, 301 insertions(+), 153 deletions(-)

Modified: sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float.hpp
==============================================================================
--- sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float.hpp	Wed Aug 21 09:34:02 2013	(r85415)
+++ sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float.hpp	2013-08-21 12:56:33 EDT (Wed, 21 Aug 2013)	(r85416)
@@ -120,7 +120,7 @@
    typename boost::enable_if_c<
       (number_category<Float>::value == number_kind_floating_point) 
          && !is_floating_point<Float>::value
-         && (std::numeric_limits<number<Float> >::radix == 2), 
+         /*&& (std::numeric_limits<number<Float> >::radix == 2)*/, 
       cpp_bin_float&>::type assign_float(Float f)
    {
       BOOST_MATH_STD_USING
@@ -219,7 +219,165 @@
 
    cpp_bin_float& operator=(const char *s)
    {
-      boost::multiprecision::detail::convert_from_string(*this, s);
+      const char* org_s = s;
+      cpp_int n;
+      int decimal_exp = 0;
+      bool ss = false;
+      //
+      // Extract the sign:
+      //
+      if(*s == '-')
+      {
+         ss = true;
+         ++s;
+      }
+      else if(*s == '+')
+         ++s;
+      //
+      // Special cases first:
+      //
+      if((std::strcmp(s, "nan") == 0) || (std::strcmp(s, "NaN") == 0) || (std::strcmp(s, "NAN") == 0))
+      {
+         return *this = std::numeric_limits<number<cpp_bin_float<Bits> > >::quiet_NaN().backend();
+      }
+      if((std::strcmp(s, "inf") == 0) || (std::strcmp(s, "Inf") == 0) || (std::strcmp(s, "INF") == 0) || (std::strcmp(s, "infinity") == 0) || (std::strcmp(s, "Infinity") == 0) || (std::strcmp(s, "INFINITY") == 0))
+      {
+         *this = std::numeric_limits<number<cpp_bin_float<Bits> > >::infinity().backend();
+         if(ss)
+            negate();
+         return *this;
+      }
+      //
+      // Digits before the point:
+      //
+      while(*s && (*s >= '0') && (*s <= '9'))
+      {
+         n *= 10u;
+         n += *s - '0';
+         ++s;
+      }
+      // The decimal point (we really should localise this!!)
+      if(*s && (*s == '.'))
+         ++s;
+      //
+      // Digits after the point:
+      //
+      while(*s && (*s >= '0') && (*s <= '9'))
+      {
+         n *= 10u;
+         n += *s - '0';
+         --decimal_exp;
+         ++s;
+      }
+      //
+      // See if there's an exponent:
+      //
+      if(*s && ((*s == 'e') || (*s == 'E')))
+      {
+         ++s;
+         int e = 0;
+         int es = false;
+         if(*s && (*s == '-'))
+         {
+            es = true;
+            ++s;
+         }
+         else if(*s && (*s == '+'))
+            ++s;
+         while(*s && (*s >= '0') && (*s <= '9'))
+         {
+            e *= 10u;
+            e += *s - '0';
+            ++s;
+         }
+         if(es)
+            e = -e;
+         decimal_exp += e;
+      }
+      if(*s)
+      {
+         //
+         // Oops unexpected input at the end of the number:
+         //
+         BOOST_THROW_EXCEPTION(std::runtime_error("Unable to parse string as a valid floating point number."));
+      }
+      if(n == 0)
+      {
+         // Result is necessarily zero:
+         *this = 0;
+         return *this;
+      }
+      if(decimal_exp > 300)
+      {
+         //
+         // TODO, FIXME, temporary hack!!
+         boost::multiprecision::detail::convert_from_string(*this, org_s);
+      }
+      else if(decimal_exp >= 0)
+      {
+         // Nice and simple, the result is an integer...
+         n *= pow(cpp_int(10), decimal_exp);
+         exponent() = (int)Bits - 1;
+         copy_and_round(*this, n.backend());
+         if(ss != sign())
+            negate();
+      }
+      else if(decimal_exp > -300)
+      {
+         // Result is the ratio of two integers: we need to organise the
+         // division so as to produce at least an N-bit result which we can
+         // round according to the remainder.
+         cpp_int d = pow(cpp_int(10), -decimal_exp);
+         int shift = (int)Bits - msb(n) + msb(d);
+         exponent() = Bits - 1;
+         if(shift > 0)
+         {
+            n <<= shift;
+            exponent() -= shift;
+         }
+         cpp_int q, r;
+         divide_qr(n, d, q, r);
+         int gb = msb(q);
+         BOOST_ASSERT(gb >= Bits - 1);
+         //
+         // Check for rounding conditions we have to
+         // handle ourselves:
+         //
+         if(gb == Bits - 1)
+         {
+            // Exactly the right number of bits, use the remainder to round:
+            r *= 2;
+            int c = r.compare(d);
+            if(c == 0)
+            {
+               // Tie:
+               if(q.backend().limbs()[0] & 1)
+                  ++q;
+            }
+            else if(c > 0)
+               ++q;
+         }
+         else if(bit_test(q, gb - (int)Bits) && ((int)lsb(q) == (gb - (int)Bits)))
+         {
+            // Too many bits in q and the bits in q indicate a tie, but we can break that using r:
+            q >>= gb - (int)Bits + 1;
+            BOOST_ASSERT(msb(q) >= Bits - 1);
+            if(r)
+               ++q;
+            else if(q.backend().limbs()[0] & 1)
+               ++q;
+            exponent() += gb - (int)Bits + 1;
+         }
+         copy_and_round(*this, q.backend());
+         if(ss != sign())
+            negate();
+      }
+      else
+      {
+         // TODO, FIXME, temporary hack!!!
+         boost::multiprecision::detail::convert_from_string(*this, org_s);
+      }
+
       return *this;
    }
 
@@ -235,21 +393,138 @@
       if(dig == 0)
          dig = std::numeric_limits<number<cpp_bin_float<Bits> > >::max_digits10;
 
+      bool scientific = (f & std::ios_base::scientific) == std::ios_base::scientific;
+      bool fixed = !scientific && (f & std::ios_base::fixed);
+
       if(exponent() <= cpp_bin_float<Bits>::max_exponent)
       {
+         // How far to left-shift in order to demormalise the mantissa:
          int shift = (int)Bits - exponent() - 1;
-         bool fractional = ((int)eval_lsb(bits()) < shift) || (shift < 0);
-         if(!fractional)
+         if(std::abs(exponent()) < 1000)
          {
-            rep_type r(bits());
-            eval_right_shift(r, shift);
-            std::string s = r.str(0, std::ios_base::fmtflags(0));
-            boost::multiprecision::detail::format_float_string(s, s.size() - 1, dig, f, false);
+            //
+            // With a smallish exponent we can use exact integer arithmetic
+            // to figure out what to print, basically we create an N digit
+            // integer plus a remainder:
+            //
+            int digits_wanted = static_cast<int>(dig);
+            int base10_exp = exponent() >= 0 ? static_cast<int>(std::floor(0.30103 * exponent())) : static_cast<int>(std::ceil(0.30103 * exponent()));
+            //
+            // For fixed formatting we want /dig/ digits after the decimal point,
+            // so if the exponent is zero, allowing for the one digit before the
+            // decimal point, we want 1 + dig digits etc.
+            //
+            if(fixed)
+               digits_wanted += 1 + base10_exp;
+            if(scientific)
+               digits_wanted += 1;
+            if(digits_wanted < -1)
+            {
+               // Fixed precision, no significant digits, and nothing to round!
+               std::string s("0");
+               if(sign())
+                  s.insert(0, 1, '-');
+               boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, true);
+               return s;
+            }
+            //
+            // power10 is the base10 exponent we need to multiply/divide by in order
+            // to convert our denormalised number to an integer with the right number of digits:
+            //
+            int power10 = digits_wanted - base10_exp - 1;
+            cpp_int i;
+            std::string s;
+            int roundup = 0; // 0=no rounding, 1=tie, 2=up
+            do
+            {
+               //
+               // Our integer is: bits() * 2^-shift * 10^power10
+               //
+               i = bits();
+               if(shift < 0)
+               {
+                  i <<= -shift;
+                  if(power10 > 0)
+                     i *= pow(cpp_int(10), power10);
+                  else if(power10 < 0)
+                  {
+                     cpp_int r;
+                     cpp_int d = pow(cpp_int(10), -power10);
+                     divide_qr(i, d, i, r);
+                     r <<= 1;
+                     int c = r.compare(d);
+                     roundup = c < 0 ? 0 : c == 0 ? 1 : 2;
+                  }
+               }
+               else
+               {
+                  //
+                  // Our integer is bits() * 2^-shift * 10^power10
+                  //
+                  if(power10 >= 0)
+                  {
+                     if(power10)
+                        i *= pow(cpp_int(10), power10);
+                     if(shift && bit_test(i, shift - 1))
+                     {
+                        if((int)lsb(i) == shift - 1)
+                           roundup = 1;
+                        else
+                           roundup = 2;
+                     }
+                     i >>= shift;
+                  }
+                  else
+                  {
+                     cpp_int r;
+                     cpp_int d = pow(cpp_int(10), -power10);
+                     d <<= shift;
+                     divide_qr(i, d, i, r);
+                     r <<= 1;
+                     int c = r.compare(d);
+                     roundup = c < 0 ? 0 : c == 0 ? 1 : 2;
+                  }
+               }
+               s = i.str(0, std::ios_base::fmtflags(0));
+               //
+               // Check if we got the right number of digits, this
+               // is really a test of whether we calculated the
+               // decimal exponent correctly:
+               //
+               int digits_got = i ? s.size() : 0;
+               if(digits_got != digits_wanted)
+               {
+                  base10_exp += digits_got - digits_wanted;
+                  if(fixed)
+                     digits_wanted = digits_got;
+                  power10 = digits_wanted - base10_exp - 1;
+                  if(fixed)
+                     break;
+                  roundup = 0;
+               }
+               else
+                  break;
+            }
+            while(true);
+            //
+            // Check whether we need to round up: note that we could equally round up
+            // the integer /i/ above, but since we need to perform the rounding *after*
+            // the conversion to a string and the digit count check, we might as well
+            // do it hear:
+            //
+            if((roundup == 2) || ((roundup == 1) && ((s[s.size() - 1] - '0') & 1)))
+            {
+               boost::multiprecision::detail::round_string_up_at(s, s.size() - 1, base10_exp);
+            }
+
             if(sign())
-               s.insert(s.begin(), '-');
+               s.insert(0, 1, '-');
+
+            boost::multiprecision::detail::format_float_string(s, base10_exp, dig, f, false);
             return s;
          }
       }
+      // TODO, FIXME, temporary hack!!!
       return boost::multiprecision::detail::convert_to_string(*this, dig, f);
    }
 

Copied and modified: sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/test_cpp_bin_float_io.cpp (from r85253, trunk/libs/multiprecision/test/test_float_io.cpp)
==============================================================================
--- trunk/libs/multiprecision/test/test_float_io.cpp	Fri Aug  9 08:27:11 2013	(r85253, copy source)
+++ sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/test_cpp_bin_float_io.cpp	2013-08-21 12:56:33 EDT (Wed, 21 Aug 2013)	(r85416)
@@ -1,4 +1,4 @@
-// Copyright John Maddock 2011.
+// Copyright John Maddock 2013.
 
 // Use, modification and distribution are subject to the
 // Boost Software License, Version 1.0.
@@ -9,38 +9,7 @@
 #  define _SCL_SECURE_NO_WARNINGS
 #endif
 
-#if !defined(TEST_MPF_50) && !defined(TEST_CPP_DEC_FLOAT) && !defined(TEST_MPFR_50) \
-      && !defined(TEST_MPFI_50) && !defined(TEST_FLOAT128)
-#  define TEST_MPF_50
-#  define TEST_CPP_DEC_FLOAT
-#  define TEST_MPFR_50
-#  define TEST_MPFI_50
-#  define TEST_FLOAT128
-
-#ifdef _MSC_VER
-#pragma message("CAUTION!!: No backend type specified so testing everything.... this will take some time!!")
-#endif
-#ifdef __GNUC__
-#pragma warning "CAUTION!!: No backend type specified so testing everything.... this will take some time!!"
-#endif
-
-#endif
-
-#if defined(TEST_MPF_50)
-#include <boost/multiprecision/gmp.hpp>
-#endif
-#if defined(TEST_MPFR_50)
-#include <boost/multiprecision/mpfr.hpp>
-#endif
-#if defined(TEST_MPFI_50)
-#include <boost/multiprecision/mpfi.hpp>
-#endif
-#ifdef TEST_CPP_DEC_FLOAT
-#include <boost/multiprecision/cpp_dec_float.hpp>
-#endif
-#ifdef TEST_FLOAT128
-#include <boost/multiprecision/float128.hpp>
-#endif
+#include <boost/multiprecision/cpp_bin_float.hpp>
 
 #include <boost/random/mersenne_twister.hpp>
 #include <boost/random/uniform_int.hpp>
@@ -53,36 +22,6 @@
 #pragma warning(disable:4127)
 #endif
 
-#if defined(TEST_MPF_50)
-template <unsigned N, boost::multiprecision::expression_template_option ET>
-bool has_bad_bankers_rounding(const boost::multiprecision::number<boost::multiprecision::gmp_float<N>, ET>&)
-{  return true;  }
-#endif
-#if defined(TEST_FLOAT128) && defined(BOOST_INTEL)
-bool has_bad_bankers_rounding(const boost::multiprecision::float128&)
-{  return true;  }
-#endif
-template <class T>
-bool has_bad_bankers_rounding(const T&) { return false; }
-
-bool is_bankers_rounding_error(const std::string& s, const char* expect)
-{
-   // This check isn't foolproof: that would require *much* more sophisticated code!!!
-   std::string::size_type l = std::strlen(expect);
-   if(l != s.size())
-      return false;
-   std::string::size_type len = s.find('e');
-   if(len == std::string::npos)
-      len = l - 1;
-   else
-      --len;
-   if(s.compare(0, len, expect, len))
-      return false;
-   if(s[len] != expect[len] + 1)
-      return false;
-   return true;
-}
-
 void print_flags(std::ios_base::fmtflags f)
 {
    std::cout << "Formatting flags were: ";
@@ -130,20 +69,13 @@
             const char* expect = string_data[j][col];
             if(ss.str() != expect)
             {
-               if(has_bad_bankers_rounding(mp_t()) && is_bankers_rounding_error(ss.str(), expect))
-               {
-                  std::cout << "Ignoring bankers-rounding error with GMP mp_f.\n";
-               }
-               else
-               {
-                  std::cout << std::setprecision(20) << "Testing value " << val << std::endl;
-                  print_flags(f[i]);
-                  std::cout << "Precision is: " << prec << std::endl;
-                  std::cout << "Got:      " << ss.str() << std::endl;
-                  std::cout << "Expected: " << expect << std::endl;
-                  ++boost::detail::test_errors();
-                  mp_t(val).str(prec, f[i]); // for debugging
-               }
+               std::cout << std::setprecision(20) << "Testing value " << val << std::endl;
+               print_flags(f[i]);
+               std::cout << "Precision is: " << prec << std::endl;
+               std::cout << "Got:      " << ss.str() << std::endl;
+               std::cout << "Expected: " << expect << std::endl;
+               ++boost::detail::test_errors();
+               mp_t(val).str(prec, f[i]); // for debugging
             }
          }
       }
@@ -221,38 +153,11 @@
    e_type e;
    val = frexp(val, &e);
 
-   static boost::random::uniform_int_distribution<e_type> ui(0, std::numeric_limits<T>::max_exponent - 10);
+   static boost::random::uniform_int_distribution<e_type> ui(0, 250);
    return ldexp(val, ui(gen));
 }
 
 template <class T>
-struct max_digits10_proxy
-{
-   static const unsigned value = std::numeric_limits<T>::digits10 + 5;
-};
-#ifdef TEST_CPP_DEC_FLOAT
-template <unsigned D, boost::multiprecision::expression_template_option ET>
-struct max_digits10_proxy<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<D>, ET> >
-{
-   static const unsigned value = std::numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<D>, ET> >::max_digits10;
-};
-#endif
-#ifdef TEST_MPF_50
-template <unsigned D, boost::multiprecision::expression_template_option ET>
-struct max_digits10_proxy<boost::multiprecision::number<boost::multiprecision::gmp_float<D>, ET> >
-{
-   static const unsigned value = std::numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<D>, ET> >::max_digits10;
-};
-#endif
-#ifdef TEST_MPFR_50
-template <unsigned D, boost::multiprecision::expression_template_option ET>
-struct max_digits10_proxy<boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<D>, ET> >
-{
-   static const unsigned value = std::numeric_limits<boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<D>, ET> >::max_digits10;
-};
-#endif
-
-template <class T>
 void do_round_trip(const T& val, std::ios_base::fmtflags f)
 {
    std::stringstream ss;
@@ -281,7 +186,10 @@
 template <class T>
 void test_round_trip()
 {
-   for(unsigned i = 0; i < 1000; ++i)
+   std::cout << "digits = " << std::numeric_limits<T>::digits << std::endl;
+   std::cout << "digits10 = " << std::numeric_limits<T>::digits10 << std::endl;
+   std::cout << "max_digits10 = " << std::numeric_limits<T>::max_digits10 << std::endl;
+   for(unsigned i = 0; i < 10000; ++i)
    {
       T val = generate_random<T>();
       do_round_trip(val);
@@ -293,44 +201,9 @@
 
 int main()
 {
-#ifdef TEST_MPFR_50
-   test<boost::multiprecision::mpfr_float_50>();
-   test<boost::multiprecision::mpfr_float_100>();
-
-   test_round_trip<boost::multiprecision::mpfr_float_50>();
-   test_round_trip<boost::multiprecision::mpfr_float_100>();
-#endif
-#ifdef TEST_MPFI_50
-   test<boost::multiprecision::mpfr_float_50>();
-   test<boost::multiprecision::mpfr_float_100>();
-
-   test_round_trip<boost::multiprecision::mpfr_float_50>();
-   test_round_trip<boost::multiprecision::mpfr_float_100>();
-#endif
-#ifdef TEST_CPP_DEC_FLOAT
-   test<boost::multiprecision::cpp_dec_float_50>();
-   test<boost::multiprecision::cpp_dec_float_100>();
-   
-   // cpp_dec_float has extra guard digits that messes this up:
-   test_round_trip<boost::multiprecision::cpp_dec_float_50>();
-   test_round_trip<boost::multiprecision::cpp_dec_float_100>();
-#endif
-#ifdef TEST_MPF_50
-   test<boost::multiprecision::mpf_float_50>();
-   test<boost::multiprecision::mpf_float_100>();
-   /*
-   // I can't get this to work with mpf_t - mpf_str appears
-   // not to actually print enough decimal digits:
-   test_round_trip<boost::multiprecision::mpf_float_50>();
-   test_round_trip<boost::multiprecision::mpf_float_100>();
-   */
-#endif
-#ifdef TEST_FLOAT128
-   test<boost::multiprecision::float128>();
-#ifndef BOOST_INTEL
-   test_round_trip<boost::multiprecision::float128>();
-#endif
-#endif
+   using namespace boost::multiprecision;
+   test<number<cpp_bin_float<113> > >();
+   test_round_trip<number<cpp_bin_float<113> > >();
    return boost::report_errors();
 }