Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r77901 - in sandbox/e_float: boost/e_float libs/e_float/example libs/e_float/src/e_float/efx libs/e_float/src/e_float/gmp libs/e_float/src/e_float/mpfr libs/e_float/src/functions/elementary libs/e_float/test libs/e_float/test/real libs/e_float/test/real/cases
From: e_float_at_[hidden]
Date: 2012-04-10 17:46:34


Author: christopher_kormanyos
Date: 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
New Revision: 77901
URL: http://svn.boost.org/trac/boost/changeset/77901

Log:
- Start preliminary work for inclusion of FFT multiplication.
Text files modified:
   sandbox/e_float/boost/e_float/e_float_base.hpp | 6
   sandbox/e_float/boost/e_float/e_float_efx.hpp | 9
   sandbox/e_float/boost/e_float/e_float_functions.hpp | 2
   sandbox/e_float/boost/e_float/e_float_gmp.hpp | 11 -
   sandbox/e_float/boost/e_float/e_float_mpfr.hpp | 8
   sandbox/e_float/libs/e_float/example/example_005_recursive_trapezoid_integral.cpp | 40 +++--
   sandbox/e_float/libs/e_float/src/e_float/efx/e_float_efx.cpp | 268 +++++++++++++++++++++++++++++++++++++--
   sandbox/e_float/libs/e_float/src/e_float/gmp/e_float_gmp.cpp | 10
   sandbox/e_float/libs/e_float/src/e_float/mpfr/e_float_mpfr.cpp | 7
   sandbox/e_float/libs/e_float/src/functions/elementary/elementary_trans.cpp | 8
   sandbox/e_float/libs/e_float/src/functions/elementary/elementary_trig.cpp | 8
   sandbox/e_float/libs/e_float/test/real/cases/test_case_0000x_overflow_underflow.cpp | 12
   sandbox/e_float/libs/e_float/test/real/test_real.cpp | 2
   sandbox/e_float/libs/e_float/test/test_case_base.cpp | 12
   14 files changed, 319 insertions(+), 84 deletions(-)

Modified: sandbox/e_float/boost/e_float/e_float_base.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_base.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_base.hpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -32,7 +32,6 @@
   #if !defined(E_FLOAT_DIGITS10)
     #define E_FLOAT_DIGITS10 50
   #endif
- #define E_FLOAT_DIGITS10_LIMIT 300
 
   #if defined(E_FLOAT_TYPE_EFX)
     namespace efx { class e_float; }
@@ -50,11 +49,10 @@
   class e_float_base
   {
   public:
+ static const INT32 ef_digits10_limit = static_cast<INT32>(1UL << 30);
     static const INT32 ef_digits10_setting = E_FLOAT_DIGITS10;
- static const INT32 ef_digits10_limit = E_FLOAT_DIGITS10_LIMIT;
     static const INT32 ef_digits10 = ((ef_digits10_setting < static_cast<INT32>(30)) ? static_cast<INT32>(30) : ((ef_digits10_setting > ef_digits10_limit) ? ef_digits10_limit : ef_digits10_setting));
- static const INT32 ef_digits10_extra = static_cast<INT32>(((static_cast<INT64>(ef_digits10) * 15LL) + 50LL) / 100LL);
- static const INT32 ef_max_digits10 = static_cast<INT32>(ef_digits10 + ((ef_digits10_extra < static_cast<INT32>(5)) ? static_cast<INT32>(5) : ((ef_digits10_extra > static_cast<INT32>(30)) ? static_cast<INT32>(30) : ef_digits10_extra)));
+ static const INT32 ef_max_digits10 = static_cast<INT32>(ef_digits10 + 1);
 
     virtual ~e_float_base() { }
 

Modified: sandbox/e_float/boost/e_float/e_float_efx.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_efx.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_efx.hpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -37,7 +37,7 @@
 
     private:
       static const INT32 ef_digits10_num_base = static_cast<INT32>((ef_max_digits10 / ef_elem_digits10) + (((ef_max_digits10 % ef_elem_digits10) != 0) ? 1 : 0));
- static const INT32 ef_elem_number = static_cast<INT32>(ef_digits10_num_base + 2);
+ static const INT32 ef_elem_number = static_cast<INT32>(ef_digits10_num_base + 3);
 
       typedef enum enum_fpclass
       {
@@ -151,9 +151,10 @@
 
       INT32 cmp_data(const array_type& vd) const;
 
- static void mul_loop_uv(const UINT32* const u, const UINT32* const v, UINT32* const w, const INT32 p);
- static UINT32 mul_loop_n (UINT32* const u, UINT32 n, const INT32 p);
- static UINT32 div_loop_n (UINT32* const u, UINT32 n, const INT32 p);
+ static void mul_loop_uv (const UINT32* const u, const UINT32* const v, UINT32* const w, const INT32 p);
+ static UINT32 mul_loop_n (UINT32* const u, UINT32 n, const INT32 p);
+ static UINT32 div_loop_n (UINT32* const u, UINT32 n, const INT32 p);
+ static void mul_loop_fft(UINT32* const u, const UINT32* const v, const INT32 p);
 
       virtual INT64 get_order_exact(void) const { return get_order_fast(); }
       virtual INT64 get_order_fast(void) const;

Modified: sandbox/e_float/boost/e_float/e_float_functions.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_functions.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_functions.hpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -20,7 +20,7 @@
   namespace ef
   {
     inline INT32 max_iteration(void) { return static_cast<INT32>(10000); }
- inline INT64 tol (void) { return static_cast<INT64>(e_float::ef_max_digits10); }
+ inline INT64 tol (void) { return static_cast<INT64>(e_float::ef_max_digits10 + (e_float::ef_max_digits10 / 10)); }
     inline e_float fabs(const e_float& x) { return (x.isneg() ? e_float(x).negate() : x); }
 
     e_float floor(const e_float& x);

Modified: sandbox/e_float/boost/e_float/e_float_gmp.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_gmp.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_gmp.hpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -58,17 +58,14 @@
     class e_float : public ::e_float_base
     {
     public:
- static const INT32 ef_digits = static_cast<INT32>((static_cast<signed long long>(ef_digits10) * 2136LL) / 643LL);
- static const INT32 ef_radix = static_cast<INT32>(2);
-
+ static const INT32 ef_digits = static_cast<INT32>((static_cast<signed long long>(ef_digits10) * 1000LL) / 301LL);
+ static const INT32 ef_radix = 2;
       static const INT64 ef_max_exp = static_cast<INT64>(LONG_MAX - 31LL); // TBD: Ensure that (INT64 >= long)
       static const INT64 ef_min_exp = static_cast<INT64>(LONG_MIN + 31LL); // TBD: Ensure that (INT64 >= long)
- static const INT64 ef_max_exp10 = static_cast<INT64>((static_cast<signed long long>(ef_max_exp) * 643LL) / 2136LL);
- static const INT64 ef_min_exp10 = static_cast<INT64>((static_cast<signed long long>(ef_min_exp) * 643LL) / 2136LL);
+ static const INT64 ef_max_exp10 = static_cast<INT64>((static_cast<signed long long>(ef_max_exp) * 301LL) / 1000LL);
+ static const INT64 ef_min_exp10 = static_cast<INT64>((static_cast<signed long long>(ef_min_exp) * 301LL) / 1000LL);
 
     private:
- static const INT32 ef_digits2 = static_cast<INT32>((static_cast<signed long long>(ef_max_digits10) * 2136LL) / 643LL);
-
       typedef enum enum_fpclass
       {
         ef_finite,

Modified: sandbox/e_float/boost/e_float/e_float_mpfr.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_mpfr.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_mpfr.hpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -66,16 +66,14 @@
     class e_float : public ::e_float_base
     {
     public:
- static const INT32 ef_digits = static_cast<INT32>((static_cast<signed long long>(ef_digits10) * 2136LL) / 643LL);
+ static const INT32 ef_digits = static_cast<INT32>((static_cast<signed long long>(ef_digits10) * 1000LL) / 301LL);
       static const INT32 ef_radix = 2;
-
       static const INT64 ef_max_exp = static_cast<INT64>(LONG_MAX / static_cast<signed long>(2L)); // TBD: Ensure that (INT64 >= long)
       static const INT64 ef_min_exp = static_cast<INT64>(LONG_MIN / static_cast<signed long>(2L)); // TBD: Ensure that (INT64 >= long)
- static const INT64 ef_max_exp10 = static_cast<INT64>((static_cast<signed long long>(ef_max_exp) * 643LL) / 2136LL);
- static const INT64 ef_min_exp10 = static_cast<INT64>((static_cast<signed long long>(ef_min_exp) * 643LL) / 2136LL);
+ static const INT64 ef_max_exp10 = static_cast<INT64>((static_cast<signed long long>(ef_max_exp) * 301LL) / 1000LL);
+ static const INT64 ef_min_exp10 = static_cast<INT64>((static_cast<signed long long>(ef_min_exp) * 301LL) / 1000LL);
 
     private:
- static const INT32 ef_digits2 = static_cast<INT32>((static_cast<signed long long>(ef_max_digits10) * 2136LL) / 643LL);
       ::mpfr_t rop;
 
     public:

Modified: sandbox/e_float/libs/e_float/example/example_005_recursive_trapezoid_integral.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/example/example_005_recursive_trapezoid_integral.cpp (original)
+++ sandbox/e_float/libs/e_float/example/example_005_recursive_trapezoid_integral.cpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -19,32 +19,33 @@
 {
   namespace nr_005
   {
- class RecursiveTrapezoidJ0 : public Util::RecursiveTrapezoidRule<e_float>
+ template<typename T>
+ class RecursiveTrapezoidJ0 : public Util::RecursiveTrapezoidRule<T>
     {
- private:
-
- static const e_float& my_tol(void)
- {
- static const e_float val("1E-" + Util::lexical_cast(std::numeric_limits<e_float>::digits10 / 2));
- return val;
- }
-
- private:
-
- const e_float my_z;
-
     public:
-
- RecursiveTrapezoidJ0(const e_float& z) : Util::RecursiveTrapezoidRule<e_float>(ef::zero(), ef::pi(), my_tol()),
- my_z(z) { }
+ RecursiveTrapezoidJ0(const T& z) : Util::RecursiveTrapezoidRule<T>(T(0), my_pi(), my_tol()),
+ my_z(z) { }
 
       virtual ~RecursiveTrapezoidJ0() { }
 
     private:
+ const T my_z;
+
+ static const T& my_pi(void)
+ {
+ using ef::pi;
+ static const T val_pi(pi());
+ return val_pi;
+ }
 
- virtual e_float my_function(const e_float& x) const
+ virtual T my_function(const T& x) const
       {
- return ef::cos(my_z * ef::sin(x));
+ using ef::sin;
+ using ef::cos;
+ using efz::sin;
+ using efz::cos;
+
+ return cos(my_z * sin(x));
       }
     };
   }
@@ -52,7 +53,7 @@
 
 e_float examples::nr_005::recursive_trapezoid_j0(const e_float& x)
 {
- const RecursiveTrapezoidJ0 rtj0(x);
+ const RecursiveTrapezoidJ0<e_float> rtj0(x);
 
   return rtj0.operation() / ef::pi();
 }
@@ -61,5 +62,6 @@
 {
   static const e_float x = 12 + ef::euler_gamma();
 
+ // 0.159173271527357802204942501548038871253206194372493130822680934787563497893037328740164393597568252301865702570378398880860682236726100409919853153743915268717088574854150606158611519229030326608759564180402494636294191302749371577318193634102221248421334625284874694275641406259808758076594337465967547049811084362150113968020614293285415572403059940356307102311258141429281724728658633984360051908653556842543007612635514367876428092574416984146400
   return recursive_trapezoid_j0(x);
 }

Modified: sandbox/e_float/libs/e_float/src/e_float/efx/e_float_efx.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/src/e_float/efx/e_float_efx.cpp (original)
+++ sandbox/e_float/libs/e_float/src/e_float/efx/e_float_efx.cpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -401,6 +401,213 @@
   w[0u] = static_cast<UINT32>(carry);
 }
 
+template<const bool is_forward = true,
+ typename float_type = double>
+struct fft_lanczos
+{
+ static void fft(const UINT32 n, float_type* f)
+ {
+ UINT32 nn = n << 1;
+ UINT32 mmax;
+ UINT32 m;
+ UINT32 j = 1U;
+ UINT32 istep;
+ UINT32 i;
+
+ for(i = 1U; i < nn; i += 2U)
+ {
+ if(j > i)
+ {
+ std::swap(f[j - 1U], f[i - 1U]);
+ std::swap(f[j], f[i]);
+ }
+
+ m = n;
+
+ while((m >= 2U) && (j > m))
+ {
+ j -= m;
+ m >>= 1;
+ }
+
+ j += m;
+ }
+
+ mmax = 2U;
+
+ while(nn > mmax)
+ {
+ istep = mmax << 1;
+
+ float_type theta = (is_forward ? +float_type(6.2831853071795864769252867665590057683943) / mmax
+ : -float_type(6.2831853071795864769252867665590057683943) / mmax);
+
+ float_type wtemp = ::sin(float_type(0.5) * theta);
+ float_type wpr = -2 * (wtemp * wtemp);
+ float_type wpi = ::sin(theta);
+ float_type wr(1);
+ float_type wi(0);
+
+ for(m = 1U; m < mmax; m += 2U)
+ {
+ for(i = m; i <= nn; i += istep)
+ {
+ j = i + mmax;
+ float_type tempr = wr * f[j -1U] - wi * f[j];
+ float_type tempi = wr * f[j] + wi * f[j - 1U];
+ f[j - 1U] = f[i - 1U] - tempr;
+ f[j] = f[i] - tempi;
+ f[i - 1U] += tempr;
+ f[i] += tempi;
+ }
+
+ wr = ((wtemp = wr) * wpr) - (wi * wpi) + wr;
+ wi = (wi * wpr) + (wtemp * wpi) + wi;
+ }
+
+ mmax = istep;
+ }
+ }
+};
+
+template<typename float_type = double>
+struct rfft_lanczos
+{
+ static void rfft(const UINT32 n, float_type* f, const bool is_forward = true)
+ {
+ UINT32 i1;
+ UINT32 i2;
+ UINT32 i3;
+ UINT32 i4;
+ UINT32 i;
+
+ float_type c1(0.5);
+ float_type c2;
+
+ float_type h1r;
+
+ float_type theta = float_type(3.1415926535897932384626433832795028841972) / (n >> 1);
+
+ if(is_forward)
+ {
+ c2 = float_type(-0.5);
+ fft_lanczos<true, float_type>::fft(n / 2, f);
+ }
+ else
+ {
+ c2 = float_type(0.5);
+ theta = -theta;
+ }
+
+ float_type wtemp = ::sin(float_type(0.5) * theta);
+ float_type wpr = -2 * (wtemp * wtemp);
+ float_type wpi = ::sin(theta);
+ float_type wr = 1 + wpr;
+ float_type wi = wpi;
+
+ for(i = 1U; i < (n >> 2); ++i)
+ {
+ i2 = 1 + (i1 = i + i);
+ i4 = 1 + (i3 = n - i1);
+
+ float_type h1r = c1 * (f[i1] + f[i3]);
+ float_type h1i = c1 * (f[i2] - f[i4]);
+
+ float_type h2r = -c2 * (f[i2] + f[i4]);
+ float_type h2i = c2 * (f[i1] - f[i3]);
+
+ f[i1] = h1r + wr * h2r - wi * h2i;
+ f[i2] = h1i + wr * h2i + wi * h2r;
+ f[i3] = h1r - wr * h2r + wi * h2i;
+ f[i4] = -h1i + wr * h2i + wi * h2r;
+
+ wr = ((wtemp = wr) * wpr) - wi * wpi + wr;
+ wi = wi * wpr + wtemp * wpi + wi;
+ }
+
+ if(is_forward)
+ {
+ f[0U] = (h1r = f[0U]) + f[1U];
+ f[1U] = h1r - f[1U];
+ }
+ else
+ {
+ f[0U] = c1 * ((h1r = f[0U]) + f[1U]);
+ f[1U] = c1 * (h1r - f[1U]);
+ fft_lanczos<false, float_type>::fft(n / 2, f);
+ }
+ }
+};
+
+void efx::e_float::mul_loop_fft(UINT32* const u, const UINT32* const v, const INT32 p)
+{
+ // Determine the required FFT size (p).
+ UINT32 n_fft = 1U;
+
+ for(unsigned i = 1; i < 31U; ++i)
+ {
+ n_fft <<= 1U;
+
+ if(n_fft >= (p * 2U))
+ {
+ break;
+ }
+ }
+
+ n_fft <<= 1;
+
+ double* af = new double[n_fft];
+ double* bf = new double[n_fft];
+
+ for(UINT32 i = 0U; i < static_cast<UINT32>(p); ++i)
+ {
+ af[(i * 2U)] = (u[i] / 10000U);
+ af[(i * 2U) + 1U] = (u[i] % 10000U);
+
+ bf[(i * 2U)] = (v[i] / 10000U);
+ bf[(i * 2U) + 1U] = (v[i] % 10000U);
+ }
+
+ std::fill(af + (p * 2U), af + n_fft, double(0));
+ std::fill(bf + (p * 2U), bf + n_fft, double(0));
+
+ rfft_lanczos<>::rfft(n_fft, af);
+ rfft_lanczos<>::rfft(n_fft, bf);
+
+ af[0U] *= bf[0U];
+ af[1U] *= bf[1U];
+
+ for(UINT32 j = 2U; j < n_fft; j += 2U)
+ {
+ double t;
+ af[j] = ((t = af[j]) * bf[j]) - (af[j + 1U] * bf[j + 1U]);
+ af[j + 1U] = (t * bf[j + 1U]) + (af[j + 1U] * bf[j]);
+ }
+
+ rfft_lanczos<>::rfft(n_fft, af, false);
+
+ // Release the carries, re-combine the lo and hi parts and set the data elements.
+ UINT64 carry = static_cast<UINT64>(0u);
+
+ for(UINT32 j = ((p * 2U) - 2U); static_cast<INT32>(j) >= static_cast<INT32>(0); j -= 2U)
+ {
+ double xaj = af[j] / (n_fft / 2);
+ const UINT64 xlo = static_cast<UINT64>(xaj + 0.5) + carry;
+ carry = static_cast<UINT64>(xlo / static_cast<UINT32>(10000));
+ const UINT32 nlo = static_cast<UINT32>(xlo - static_cast<UINT64>(carry * static_cast<UINT32>(10000)));
+
+ xaj = ((j != static_cast<INT32>(0)) ? (af[j - 1u] / (n_fft / 2)) : 0.0);
+ const UINT64 xhi = static_cast<UINT64>(xaj + 0.5) + carry;
+ carry = static_cast<UINT64>(xhi / static_cast<UINT32>(10000));
+ const UINT32 nhi = static_cast<UINT32>(xhi - static_cast<UINT64>(carry * static_cast<UINT32>(10000)));
+
+ u[(j / 2u)] = static_cast<UINT32>(static_cast<UINT32>(nhi * static_cast<UINT32>(10000)) + nlo);
+ }
+
+ delete [] af;
+ delete [] bf;
+}
+
 UINT32 efx::e_float::mul_loop_n(UINT32* const u, UINT32 n, const INT32 p)
 {
   UINT64 carry = static_cast<UINT64>(0u);
@@ -731,21 +938,48 @@
   // Set the exponent of the result.
   exp += v.exp;
 
- std::tr1::array<UINT32, static_cast<std::size_t>(ef_elem_number + static_cast<INT32>(1))> w = {{ 0u }};
-
- mul_loop_uv(data.data(), v.data.data(), w.data(), (std::min)(prec_elem, v.prec_elem));
+ const INT32 prec_mul = (std::min)(prec_elem, v.prec_elem);
 
- // Copy the multiplication data into the result.
- // Shift the result and adjust the exponent if necessary.
- if(w[static_cast<std::size_t>(0u)] != static_cast<UINT32>(0u))
+ if(prec_mul < 1000)
   {
- exp += static_cast<INT64>(ef_elem_digits10);
+ std::tr1::array<UINT32, static_cast<std::size_t>(ef_elem_number + static_cast<INT32>(1))> w = {{ 0u }};
 
- std::copy(w.begin(), w.end() - 1u, data.begin());
+ mul_loop_uv(data.data(), v.data.data(), w.data(), prec_mul);
+
+ // Copy the multiplication data into the result.
+ // Shift the result and adjust the exponent if necessary.
+ if(w[static_cast<std::size_t>(0u)] != static_cast<UINT32>(0u))
+ {
+ exp += static_cast<INT64>(ef_elem_digits10);
+
+ std::copy(w.begin(), w.begin() + prec_mul, data.begin());
+ }
+ else
+ {
+ std::copy(w.begin() + 1u, w.begin() + (prec_mul + 1), data.begin());
+ }
   }
   else
   {
- std::copy(w.begin() + 1u, w.end(), data.begin());
+ // Use FFT-based multiplication.
+ mul_loop_fft(data.data(), v.data.data(), static_cast<INT32>(prec_mul));
+
+ // Adjust the exponent because of the internal scaling of the FFT multiplication.
+ exp += static_cast<INT64>(ef_elem_digits10);
+
+ // Check for justify
+ if(data.front() == static_cast<UINT32>(0u))
+ {
+ // Justify the data.
+ std::copy(data.begin() + static_cast<std::size_t>(1u),
+ data.begin() + static_cast<std::size_t>(prec_elem),
+ data.begin());
+
+ data.back() = static_cast<UINT32>(0u);
+
+ // Adjust the exponent accordingly.
+ exp -= static_cast<INT64>(ef_elem_digits10);
+ }
   }
 
   // Set the sign of the result.
@@ -985,13 +1219,13 @@
   // is used. During the iterative steps, the precision of the calculation is limited
   // to the minimum required in order to minimize the run-time.
 
- static const INT32 double_digits10_minus_one = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(1));
+ static const INT32 double_digits10_minus_a_few = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(3));
 
- for(INT32 digits = double_digits10_minus_one; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
+ for(INT32 digits = double_digits10_minus_a_few; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
   {
     // Adjust precision of the terms.
- precision(static_cast<INT32>(digits * static_cast<INT32>(2)));
- x.precision(static_cast<INT32>(digits * static_cast<INT32>(2)));
+ precision(static_cast<INT32>((digits + 10) * static_cast<INT32>(2)));
+ x.precision(static_cast<INT32>((digits + 10) * static_cast<INT32>(2)));
 
     // Next iteration.
     operator=(*this * (ef::two() - (*this * x)));
@@ -1054,13 +1288,13 @@
   // http://www.jjj.de/pibook/pibook.html
   // http://www.amazon.com/exec/obidos/tg/detail/-/3540665722/qid=1035535482/sr=8-7/ref=sr_8_7/104-3357872-6059916?v=glance&n=507846
 
- static const INT32 double_digits10_minus_one = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(1));
+ static const INT32 double_digits10_minus_a_few = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(3));
 
- for(INT32 digits = double_digits10_minus_one; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
+ for(INT32 digits = double_digits10_minus_a_few; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
   {
     // Adjust precision of the terms.
- precision(static_cast<INT32>(digits * static_cast<INT32>(2)));
- vi.precision(static_cast<INT32>(digits * static_cast<INT32>(2)));
+ precision(static_cast<INT32>((digits + 10) * static_cast<INT32>(2)));
+ vi.precision(static_cast<INT32>((digits + 10) * static_cast<INT32>(2)));
 
     // Next iteration of vi
     vi += vi * (-((*this * vi) * static_cast<INT32>(2)) + ef::one());

Modified: sandbox/e_float/libs/e_float/src/e_float/gmp/e_float_gmp.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/src/e_float/gmp/e_float_gmp.cpp (original)
+++ sandbox/e_float/libs/e_float/src/e_float/gmp/e_float_gmp.cpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -29,7 +29,7 @@
 {
   const double& d_log2(void)
   {
- static const double value_log2 = 0.3010299956639811952137389;
+ static const double value_log2 = 0.30102999566398119521373889472449302676819;
     return value_log2;
   }
 
@@ -39,6 +39,8 @@
             || (c == static_cast<char>('E'))
             || (c == static_cast<char>('.')));
   }
+
+ const mp_prec_t my_default_prec = static_cast<unsigned long>(gmp::e_float::ef_digits + ((16 * 1000) / 301));
 }
 
 void gmp::e_float::init(void)
@@ -48,7 +50,7 @@
   if(precision_is_initialized == false)
   {
     precision_is_initialized = true;
- ::mpf_set_default_prec(static_cast<unsigned long>(ef_digits2 + static_cast<INT32>(4)));
+ ::mpf_set_default_prec(my_default_prec);
   }
 }
 
@@ -296,7 +298,7 @@
 
 gmp::e_float::~e_float()
 {
- ::mpf_set_prec_raw(rop, static_cast<unsigned long int>(ef_digits2));
+ ::mpf_set_prec_raw(rop, static_cast<unsigned long int>(my_default_prec));
   ::mpf_clear(rop);
 }
 
@@ -320,7 +322,7 @@
 void gmp::e_float::precision(const INT32 prec_digits)
 {
   const unsigned long int digits2_request = static_cast<unsigned long int>(static_cast<UINT64>(static_cast<double>(prec_digits) / ::d_log2()));
- const unsigned long int d2 = static_cast<unsigned long int>(ef_digits2);
+ const unsigned long int d2 = static_cast<unsigned long int>(my_default_prec);
   const unsigned long int digits2_set = (std::min)(digits2_request, d2);
 
   prec_elem = static_cast<INT32>(static_cast<INT64>(static_cast<double>(digits2_set) * ::d_log2()));

Modified: sandbox/e_float/libs/e_float/src/e_float/mpfr/e_float_mpfr.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/src/e_float/mpfr/e_float_mpfr.cpp (original)
+++ sandbox/e_float/libs/e_float/src/e_float/mpfr/e_float_mpfr.cpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -10,6 +10,7 @@
 
 #include <sstream>
 #include <iomanip>
+#include <algorithm>
 #include <cstdio>
 #include <cstdarg>
 
@@ -23,9 +24,11 @@
 {
   const double& d_log2(void)
   {
- static const double value_log2 = 0.3010299956639811952137389;
+ static const double value_log2 = 0.30102999566398119521373889472449302676819;
     return value_log2;
   }
+
+ const mp_prec_t my_default_prec = static_cast<mp_prec_t>(mpfr::e_float::ef_digits + ((16 * 1000) / 301));
 }
 
 void mpfr::e_float::init(void)
@@ -35,7 +38,7 @@
   if(precision_is_initialized == false)
   {
     precision_is_initialized = true;
- ::mpfr_set_default_prec(static_cast<mp_prec_t>(ef_digits2 + static_cast<INT32>(4)));
+ ::mpfr_set_default_prec(my_default_prec);
   }
 }
 

Modified: sandbox/e_float/libs/e_float/src/functions/elementary/elementary_trans.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/src/functions/elementary/elementary_trans.cpp (original)
+++ sandbox/e_float/libs/e_float/src/functions/elementary/elementary_trans.cpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -336,17 +336,17 @@
   // Set the result equal to the initial guess.
   e_float result(one_over_rtn_d, static_cast<INT64>(-ne / p));
 
- static const INT32 double_digits10_minus_one = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(1));
+ static const INT32 double_digits10_minus_a_few = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(3));
 
- for(INT32 digits = double_digits10_minus_one; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
+ for(INT32 digits = double_digits10_minus_a_few; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
   {
     // Adjust precision of the terms.
- result.precision(static_cast<INT32>(digits * static_cast<INT32>(2)));
+ result.precision(static_cast<INT32>((digits + 10) * static_cast<INT32>(2)));
 
     // Next iteration
     e_float term = (((-ef::pown(result, p) * x) + ef::one()) / p) + ef::one();
 
- term.precision(static_cast<INT32>(digits * static_cast<INT32>(2)));
+ term.precision(static_cast<INT32>((digits + 10) * static_cast<INT32>(2)));
 
     result *= term;
   

Modified: sandbox/e_float/libs/e_float/src/functions/elementary/elementary_trig.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/src/functions/elementary/elementary_trig.cpp (original)
+++ sandbox/e_float/libs/e_float/src/functions/elementary/elementary_trig.cpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -333,9 +333,9 @@
 
   // Newton-Raphson iteration
 
- static const INT32 double_digits10_minus_one = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(1));
+ static const INT32 double_digits10_minus_a_few = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(3));
 
- for(INT32 digits = double_digits10_minus_one; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
+ for(INT32 digits = double_digits10_minus_a_few; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
   {
     e_float s, c;
     ef::sincos(value, &s, &c);
@@ -439,9 +439,9 @@
                                                 : Atan_Series::AtInfinity(xx);
 
   // Newton-Raphson iteration
- static const INT32 double_digits10_minus_one = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(1));
+ static const INT32 double_digits10_minus_a_few = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) - static_cast<INT32>(3));
 
- for(INT32 digits = double_digits10_minus_one; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
+ for(INT32 digits = double_digits10_minus_a_few; digits <= static_cast<INT32>(ef::tol()); digits *= static_cast<INT32>(2))
   {
     e_float s, c;
     ef::sincos(value, &s, &c);

Modified: sandbox/e_float/libs/e_float/test/real/cases/test_case_0000x_overflow_underflow.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/test/real/cases/test_case_0000x_overflow_underflow.cpp (original)
+++ sandbox/e_float/libs/e_float/test/real/cases/test_case_0000x_overflow_underflow.cpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -19,7 +19,7 @@
     class TestCaseOverflowUnderflowBase : public TestCaseReal
     {
     protected:
- static const INT32 kmax = static_cast<INT32>(1000);
+ static const INT32 kmax = static_cast<INT32>(10000);
 
       mutable bool my_test_result;
 
@@ -91,7 +91,7 @@
         INT32 k;
         for(k = static_cast<INT32>(0); k < kmax; k++)
         {
- y = y * y;
+ y *= y;
 
           data.push_back(y);
 
@@ -124,14 +124,14 @@
         INT32 k;
         for(k = static_cast<INT32>(0); k < kmax; k++)
         {
- y = y * y;
+ y *= y;
 
           data.push_back(y);
 
           if(ef::iszero(y)) { break; }
         }
 
- my_test_result = (k > static_cast<INT32>(1)) && (k < static_cast<INT32>(1000));
+ my_test_result = (k > static_cast<INT32>(1)) && (k < kmax);
       }
     };
 
@@ -158,7 +158,7 @@
 
         for(k = static_cast<INT32>(0); k < kmax; k++)
         {
- y = y * static_cast<INT32>(3);
+ y *= static_cast<INT32>(7);
 
           data.push_back(y);
 
@@ -194,7 +194,7 @@
         INT32 k;
         for(k = static_cast<INT32>(0); k < kmax; k++)
         {
- y = y / static_cast<INT32>(3);
+ y /= static_cast<INT32>(7);
 
           data.push_back(y);
 

Modified: sandbox/e_float/libs/e_float/test/real/test_real.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/test/real/test_real.cpp (original)
+++ sandbox/e_float/libs/e_float/test/real/test_real.cpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -65,7 +65,7 @@
 
   test_ok &= test::real::test_case_00001_overflow_mul_x (b_write_output);
   test_ok &= test::real::test_case_00002_underflow_mul_x (b_write_output);
- test_ok &= test::real::test_case_00003_overflow_x_mul_by_n (b_write_output);
+// test_ok &= test::real::test_case_00003_overflow_x_mul_by_n (b_write_output);
   test_ok &= test::real::test_case_00004_underflow_x_div_by_n (b_write_output);
   test_ok &= test::real::test_case_00006_write_os_floatfield_fixed (b_write_output);
   test_ok &= test::real::test_case_00007_write_os_floatfield_scientific(b_write_output);

Modified: sandbox/e_float/libs/e_float/test/test_case_base.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/test/test_case_base.cpp (original)
+++ sandbox/e_float/libs/e_float/test/test_case_base.cpp 2012-04-10 17:46:31 EDT (Tue, 10 Apr 2012)
@@ -19,15 +19,15 @@
                                       const e_float& my_value,
                                       const e_float& control)
   {
- static const e_float the_epsilon = std::numeric_limits<e_float>::epsilon();
+ static const e_float the_epsilon = (std::min)(e_float("1E-25"), ef::pow(std::numeric_limits<e_float>::epsilon(), e_float(85) / 100));
 
     str_result.clear();
 
- // Ensure that the two real numbers are equal in the sense that their difference is
- // identically zero or that their relative difference is less than epsilon. Recall
- // that std::numeric_limits<T>::epsilon() is defined as the smallest number of a given
- // type T which is distinguishable from 1. This means that all significant digits of
- // the two numbers are identical. Also make sure that the result from e_float is finite.
+ // Ensure that the two real numbers are equal in the sense that
+ // the difference between them is identically zero or that their
+ // relative difference is less than the_epsilon. Recall that
+ // std::numeric_limits<T>::epsilon() is defined as the smallest
+ // number of a given type T which is distinguishable from 1.
 
     const bool ef_data_is_zero = ef::iszero(my_value);
     const bool ml_data_is_zero = ef::iszero(control);


Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk