[Boost-commit] svn:boost r86473 - in sandbox/multiprecision.cpp_bin_float: boost/multiprecision boost/multiprecision/cpp_bin_float libs/multiprecision/doc libs/multiprecision/test/math libs/multiprecision/test/math/instances

Subject: [Boost-commit] svn:boost r86473 - in sandbox/multiprecision.cpp_bin_float: boost/multiprecision boost/multiprecision/cpp_bin_float libs/multiprecision/doc libs/multiprecision/test/math libs/multiprecision/test/math/instances
From: john_at_[hidden]
Date: 2013-10-27 05:27:01

Author: johnmaddock
Date: 2013-10-27 05:27:01 EDT (Sun, 27 Oct 2013)
New Revision: 86473
URL: http://svn.boost.org/trac/boost/changeset/86473

Log:
Lots of small fixes for better overflow/underflow handling.
Add some optimised sign manipulation routines.
Add a new version of exp optimised for binary floating point types.
Tweak test types.

Added:
   sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float/transcendental.hpp   (contents, props changed)
Text files modified: 
   sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float.hpp                |    54 ++++++++++++++++---                     
   sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float/io.hpp             |    70 +++++++++++++++++++++++---              
   sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float/transcendental.hpp |   106 ++++++++++++++++++++++++++++++++++++++++
   sandbox/multiprecision.cpp_bin_float/libs/multiprecision/doc/multiprecision.qbk            |     2                                         
   sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/math/instances/instances.cpp |    14 ++++                                    
   sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/math/setup.hpp               |     4                                         
   6 files changed, 227 insertions(+), 23 deletions(-)

Modified: sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float.hpp
==============================================================================
--- sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float.hpp	Sun Oct 27 05:23:16 2013	(r86472)
+++ sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float.hpp	2013-10-27 05:27:01 EDT (Sun, 27 Oct 2013)	(r86473)
@@ -163,7 +163,7 @@
          m_exponent += bits;
          eval_add(*this, ipart);
       }
-      m_exponent += e;
+      m_exponent += static_cast<Exponent>(e);
       return *this;
    }
 
@@ -263,12 +263,12 @@
          unsigned shift = msb(fi);
          if(shift >= bit_count)
          {
-            m_exponent = shift;
+            m_exponent = static_cast<Exponent>(shift);
             m_data = static_cast<limb_type>(fi >> (shift + 1 - bit_count));
          }
          else
          {
-            m_exponent = shift;
+            m_exponent = static_cast<Exponent>(shift);
             eval_left_shift(m_data, bit_count - shift - 1);
          }
          BOOST_ASSERT(eval_bit_test(m_data, bit_count-1));
@@ -370,12 +370,12 @@
       res.bits() = static_cast<limb_type>(0u);
       return;
    }
-   unsigned msb = eval_msb(arg);
+   int msb = eval_msb(arg);
    if(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count > msb + 1)
    {
       // Must have had cancellation in subtraction, shift left and copy:
       eval_left_shift(arg, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb - 1);
-      res.exponent() -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb - 1;
+      res.exponent() -= static_cast<Exponent>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb - 1);
    }
    else if(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count < msb + 1)
    {
@@ -391,7 +391,7 @@
       }
       // Shift off the cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count we don't need:
       eval_right_shift(arg, msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1);
-      res.exponent() += msb - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
+      res.exponent() += static_cast<Exponent>(msb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1);
       if(roundup)
       {
          eval_increment(arg);
@@ -912,7 +912,7 @@
    // We can set the exponent and sign of the result up front:
    //
    int gb = msb(v);
-   res.exponent() = u.exponent() - gb - 1;
+   res.exponent() = u.exponent() - static_cast<Exponent>(gb) - static_cast<Exponent>(1);
    res.sign() = u.sign();
    //
    // Now get the quotient and remainder:
@@ -982,6 +982,27 @@
 }
 
 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
+inline int eval_get_sign(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
+{
+   return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero ? 0 : arg.sign() ? -1 : 1;
+}
+
+template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
+inline bool eval_is_zero(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
+{
+   return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
+}
+
+template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
+inline bool eval_eq(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &a, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &b)
+{
+   return (a.exponent() == b.exponent())
+      && (a.sign() == b.sign())
+      && (a.bits().compare(b.bits()) == 0)
+      && (a.exponent() != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan);
+}
+
+template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
 inline void eval_convert_to(long long *res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
 {
    switch(arg.exponent())
@@ -1120,8 +1141,22 @@
       res = arg;
       return;
    }
-   res = arg;
-   res.exponent() += e;
+   if((e > 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent - e < arg.exponent()))
+   {
+      // Overflow:
+      res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
+      res.sign() = arg.sign();
+   }
+   else if((e < 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - e > arg.exponent()))
+   {
+      // Underflow:
+      res = limb_type(0);
+   }
+   else
+   {
+      res = arg;
+      res.exponent() += e;
+   }
 }
 
 template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
@@ -1339,6 +1374,7 @@
 }} // namespaces
 
 #include <boost/multiprecision/cpp_bin_float/io.hpp>
+#include <boost/multiprecision/cpp_bin_float/transcendental.hpp>
 
 namespace std{
 

Modified: sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float/io.hpp
==============================================================================
--- sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float/io.hpp	Sun Oct 27 05:23:16 2013	(r86472)
+++ sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float/io.hpp	2013-10-27 05:27:01 EDT (Sun, 27 Oct 2013)	(r86473)
@@ -286,6 +286,7 @@
 #endif
    boost::int64_t error = 0;
    boost::intmax_t calc_exp = 0;
+   boost::intmax_t final_exponent = 0;
 
    if(decimal_exp >= 0)
    {
@@ -300,13 +301,11 @@
          }
          else
             t = n;
-         exponent() = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1;
-         exponent() += static_cast<exponent_type>(decimal_exp);
-         exponent() += static_cast<exponent_type>(calc_exp);
+         final_exponent = (boost::int64_t)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp + calc_exp;
          int rshift = msb(t) - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
          if(rshift > 0)
          {
-            exponent() += rshift;
+            final_exponent += rshift;
             int roundup = boost::multiprecision::cpp_bf_io_detail::get_round_mode(t, rshift - 1, error);
             t >>= rshift;
             if((roundup == 2) || ((roundup == 1) && t.backend().limbs()[0] & 1))
@@ -326,6 +325,22 @@
          {
             BOOST_ASSERT(!error);
          }
+         if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
+         {
+            exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
+            final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
+         }
+         else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
+         {
+            // Underflow:
+            exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
+            final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
+         }
+         else
+         {
+            exponent() = static_cast<Exponent>(final_exponent);
+            final_exponent = 0;
+         }
          copy_and_round(*this, t.backend());
          break;
       }
@@ -345,11 +360,11 @@
          cpp_int d;
          calc_exp = boost::multiprecision::cpp_bf_io_detail::restricted_pow(d, cpp_int(5), -decimal_exp, max_bits, error);
          int shift = (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - msb(n) + msb(d);
-         exponent() = static_cast<exponent_type>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp - calc_exp);
+         final_exponent = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 + decimal_exp - calc_exp;
          if(shift > 0)
          {
             n <<= shift;
-            exponent() -= shift;
+            final_exponent -= static_cast<Exponent>(shift);
          }
          cpp_int q, r;
          divide_qr(n, d, q, r);
@@ -371,7 +386,7 @@
             // note that the radius of error in r is error/2 * q:
             int shift = gb - (int)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
             q >>= shift;
-            exponent() += shift;
+            final_exponent += static_cast<Exponent>(shift);
             BOOST_ASSERT((msb(q) >= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
             if(error && (r < (error / 2) * q))
                roundup = -1;
@@ -396,12 +411,29 @@
             if(shift > 0)
             {
                n >>= shift;
-               exponent() += shift;
+               final_exponent += static_cast<Exponent>(shift);
             }
             continue;
          }
          else if((roundup == 2) || ((roundup == 1) && q.backend().limbs()[0] & 1))
             ++q;
+         if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
+         {
+            // Overflow:
+            exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
+            final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
+         }
+         else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
+         {
+            // Underflow:
+            exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
+            final_exponent -= cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
+         }
+         else
+         {
+            exponent() = static_cast<Exponent>(final_exponent);
+            final_exponent = 0;
+         }
          copy_and_round(*this, q.backend());
          if(ss != sign())
             negate();
@@ -409,7 +441,27 @@
       }
       while(true);
    }
-
+   //
+   // Check for scaling and/or over/under-flow:
+   //
+   final_exponent += exponent();
+   if(final_exponent > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
+   {
+      // Overflow:
+      exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
+      bits() = limb_type(0);
+   }
+   else if(final_exponent < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
+   {
+      // Underflow:
+      exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
+      bits() = limb_type(0);
+      sign() = 0;
+   }
+   else
+   {
+      exponent() = static_cast<Exponent>(final_exponent);
+   }
    return *this;
 }
 

Added: sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float/transcendental.hpp
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ sandbox/multiprecision.cpp_bin_float/boost/multiprecision/cpp_bin_float/transcendental.hpp	2013-10-27 05:27:01 EDT (Sun, 27 Oct 2013)	(r86473)
@@ -0,0 +1,106 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2013 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_
+
+#ifndef BOOST_MULTIPRECISION_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
+#define BOOST_MULTIPRECISION_CPP_BIN_FLOAT_TRANSCENDENTAL_HPP
+
+namespace boost{ namespace multiprecision{ namespace backends{
+
+template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
+void eval_exp_taylor(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
+{
+   static const int bits = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count;
+   //
+   // Taylor series for small argument:
+   //
+   res = limb_type(0);
+   cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> num(arg), denom, t;
+   denom = limb_type(1);
+   eval_add(res, num);
+
+   for(unsigned k = 2; ; ++k)
+   {
+      eval_multiply(denom, k);
+      eval_multiply(num, arg);
+      eval_divide(t, num, denom);
+      eval_add(res, t);
+      if(eval_is_zero(t) || (res.exponent() - bits > t.exponent()))
+         break;
+   }
+}
+
+template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
+void eval_exp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> &arg)
+{
+   using default_ops::eval_multiply;
+   using default_ops::eval_subtract;
+   using default_ops::eval_add;
+   using default_ops::eval_convert_to;
+
+   int type = eval_fpclassify(arg);
+   bool isneg = eval_get_sign(arg) < 0;
+   if(type == FP_NAN)
+   {
+      res = arg;
+      return;
+   }
+   else if(type == FP_INFINITE)
+   {
+      res = arg;
+      if(isneg)
+         res = limb_type(0u);
+      else 
+         res = arg;
+      return;
+   }
+   else if(type == FP_ZERO)
+   {
+      res = limb_type(1);
+      return;
+   }
+   cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> t, n;
+   if(isneg)
+   {
+      t = arg;
+      t.negate();
+      eval_exp(res, t);
+      t.swap(res);
+      res = limb_type(1);
+      eval_divide(res, t);
+      return;
+   }
+
+   eval_divide(n, arg, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >());
+   eval_floor(n, n);
+   eval_multiply(t, n, default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >());
+   eval_subtract(t, arg);
+   t.negate();
+   BOOST_ASSERT(t.compare(limb_type(0)) >= 0);
+   BOOST_ASSERT(t.compare(default_ops::get_constant_ln2<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> >()) < 0);
+
+   Exponent k, nn;
+   eval_convert_to(&nn, n);
+   k = nn ? Exponent(1) << (msb(nn) / 2) : 0;
+   eval_ldexp(t, t, -k);
+
+   eval_exp_taylor(res, t);
+   //
+   // Square 1 + res k times:
+   //
+   for(int s = 0; s < k; ++s)
+   {
+      t.swap(res);
+      eval_multiply(res, t, t);
+      eval_ldexp(t, t, 1);
+      eval_add(res, t);
+   }
+   eval_add(res, limb_type(1));
+   eval_ldexp(res, res, nn);
+}
+
+}}} // namespaces
+
+#endif
+

Modified: sandbox/multiprecision.cpp_bin_float/libs/multiprecision/doc/multiprecision.qbk
==============================================================================
--- sandbox/multiprecision.cpp_bin_float/libs/multiprecision/doc/multiprecision.qbk	Sun Oct 27 05:23:16 2013	(r86472)
+++ sandbox/multiprecision.cpp_bin_float/libs/multiprecision/doc/multiprecision.qbk	2013-10-27 05:27:01 EDT (Sun, 27 Oct 2013)	(r86473)
@@ -694,7 +694,7 @@
 Normally `cpp_bin_float` allocates no memory: all of the space required for its digits are allocated
 directly within the class.  As a result care should be taken not to use the class with too high a digit count
 as stack space requirements can grow out of control.  If that represents a problem then providing an allocator
-as the final template parameter causes `cpp_bin_float` to dynamically allocate the memory it needs: this
+as a template parameter causes `cpp_bin_float` to dynamically allocate the memory it needs: this
 significantly reduces the size of `cpp_bin_float` and increases the viable upper limit on the number of digits
 at the expense of performance.  However, please bear in mind that arithmetic operations rapidly become ['very] expensive
 as the digit count grows: the current implementation really isn't optimized or designed for large digit counts.

Modified: sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/math/instances/instances.cpp
==============================================================================
--- sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/math/instances/instances.cpp	Sun Oct 27 05:23:16 2013	(r86472)
+++ sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/math/instances/instances.cpp	2013-10-27 05:27:01 EDT (Sun, 27 Oct 2013)	(r86473)
@@ -6,9 +6,19 @@
 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
 #include <boost/math/special_functions.hpp>
 
+#define BOOST_MATH_TEST_TYPE boost::multiprecision::float128_t
+
 #include "../setup.hpp"
 
-//#define BOOST_MATH_TEST_TYPE double
-//#define TEST_GROUP_1
+#define TEST_GROUP_1
+#define TEST_GROUP_2
+#define TEST_GROUP_3
+#define TEST_GROUP_4
+#define TEST_GROUP_5
+#define TEST_GROUP_6
+#define TEST_GROUP_7
+#define TEST_GROUP_8
+#define TEST_GROUP_9
+#define TEST_GROUP_10
 #include "libs/math/test/test_instances/test_instances.hpp"
 

Modified: sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/math/setup.hpp
==============================================================================
--- sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/math/setup.hpp	Sun Oct 27 05:23:16 2013	(r86472)
+++ sandbox/multiprecision.cpp_bin_float/libs/multiprecision/test/math/setup.hpp	2013-10-27 05:27:01 EDT (Sun, 27 Oct 2013)	(r86473)
@@ -95,9 +95,9 @@
 #ifdef TEST_CPP_BIN_FLOAT
 #include <boost/multiprecision/cpp_bin_float.hpp>
 
-#define CPP_BIN_FLOAT_TESTS test(cpp_bin_float_50>(), "cpp_bin_float_50");
+#define CPP_BIN_FLOAT_TESTS test(float128_t(), "float128_t");
 
-typedef boost::multiprecision::cpp_bin_float_50 test_type_1;
+typedef boost::multiprecision::float128_t test_type_1;
 
 #else