Boost logo

Boost-Commit :

From: john_at_[hidden]
Date: 2007-10-02 05:30:55


Author: johnmaddock
Date: 2007-10-02 05:30:41 EDT (Tue, 02 Oct 2007)
New Revision: 39650
URL: http://svn.boost.org/trac/boost/changeset/39650

Log:
Added workaround for apparently broken std::fmod(long double,long double) on Darwin.
Added more tracing macros to try and track down remaining Darwin issues.
Text files modified:
   sandbox/math_toolkit/boost/math/concepts/real_concept.hpp | 2 +-
   sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp | 2 +-
   sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp | 23 +++++++++++++++++++++--
   sandbox/math_toolkit/boost/math/special_functions/ellint_1.hpp | 4 ++--
   sandbox/math_toolkit/boost/math/special_functions/ellint_2.hpp | 4 ++--
   sandbox/math_toolkit/boost/math/special_functions/ellint_3.hpp | 4 ++--
   sandbox/math_toolkit/boost/math/special_functions/spherical_harmonic.hpp | 6 +++---
   sandbox/math_toolkit/boost/math/tools/config.hpp | 23 ++++++++++++++++++++++-
   8 files changed, 54 insertions(+), 14 deletions(-)

Modified: sandbox/math_toolkit/boost/math/concepts/real_concept.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/concepts/real_concept.hpp (original)
+++ sandbox/math_toolkit/boost/math/concepts/real_concept.hpp 2007-10-02 05:30:41 EDT (Tue, 02 Oct 2007)
@@ -196,7 +196,7 @@
 inline real_concept ceil(real_concept a)
 { return std::ceil(a.value()); }
 inline real_concept fmod(real_concept a, real_concept b)
-{ return std::fmod(a.value(), b.value()); }
+{ return boost::math::tools::fmod_workaround(a.value(), b.value()); }
 inline real_concept cosh(real_concept a)
 { return std::cosh(a.value()); }
 inline real_concept exp(real_concept a)

Modified: sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp (original)
+++ sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp 2007-10-02 05:30:41 EDT (Tue, 02 Oct 2007)
@@ -195,7 +195,7 @@
 inline boost::math::concepts::std_real_concept ceil(boost::math::concepts::std_real_concept a)
 { return std::ceil(a.value()); }
 inline boost::math::concepts::std_real_concept fmod(boost::math::concepts::std_real_concept a, boost::math::concepts::std_real_concept b)
-{ return std::fmod(a.value(), b.value()); }
+{ return boost::math::tools::fmod_workaround(a.value(), b.value()); }
 inline boost::math::concepts::std_real_concept cosh(boost::math::concepts::std_real_concept a)
 { return std::cosh(a.value()); }
 inline boost::math::concepts::std_real_concept exp(boost::math::concepts::std_real_concept a)

Modified: sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp 2007-10-02 05:30:41 EDT (Tue, 02 Oct 2007)
@@ -100,7 +100,9 @@
     // modified Lentz's method, see
     // Lentz, Applied Optics, vol 15, 668 (1976)
     tolerance = 2 * tools::epsilon<T>();
+ BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
     tiny = sqrt(tools::min_value<T>());
+ BOOST_MATH_INSTRUMENT_VARIABLE(tiny);
     C = f = tiny; // b0 = 0, replace with tiny
     D = 0;
     for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
@@ -114,11 +116,13 @@
         D = 1 / D;
         delta = C * D;
         f *= delta;
- if (abs(delta - 1) < tolerance)
+ BOOST_MATH_INSTRUMENT_VARIABLE(delta-1);
+ if (abs(delta - 1) <= tolerance)
         {
            break;
         }
     }
+ BOOST_MATH_INSTRUMENT_VARIABLE(k);
     policies::check_series_iterations("boost::math::bessel_ik<%1%>(%1%,%1%)", k, pol);
 
     *fv = f;
@@ -154,6 +158,11 @@
     current = 1; // q1
     Q = C = -a; // Q1 = C1 because q1 = 1
     S = 1 + Q * delta; // S1
+ BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
+ BOOST_MATH_INSTRUMENT_VARIABLE(a);
+ BOOST_MATH_INSTRUMENT_VARIABLE(b);
+ BOOST_MATH_INSTRUMENT_VARIABLE(D);
+ BOOST_MATH_INSTRUMENT_VARIABLE(f);
     for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) // starting from 2
     {
         // continued fraction f = z1 / z0
@@ -172,6 +181,8 @@
         S += Q * delta;
 
         // S converges slower than f
+ BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta);
+ BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance);
         if (abs(Q * delta) < abs(S) * tolerance)
         {
            break;
@@ -181,6 +192,8 @@
 
     *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S;
     *Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;
+ BOOST_MATH_INSTRUMENT_VARIABLE(*Kv);
+ BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1);
 
     return 0;
 }
@@ -201,6 +214,9 @@
     T W, current, prev, next;
     bool reflect = false;
     unsigned n, k;
+ BOOST_MATH_INSTRUMENT_VARIABLE(v);
+ BOOST_MATH_INSTRUMENT_VARIABLE(x);
+ BOOST_MATH_INSTRUMENT_VARIABLE(kind);
 
     BOOST_MATH_STD_USING
     using namespace boost::math::tools;
@@ -216,6 +232,8 @@
     }
     n = tools::real_cast<unsigned>(v + 0.5f);
     u = v - n; // -1/2 <= u < 1/2
+ BOOST_MATH_INSTRUMENT_VARIABLE(n);
+ BOOST_MATH_INSTRUMENT_VARIABLE(u);
 
     if (x < 0)
     {
@@ -303,7 +321,8 @@
         *I = Iv;
         *K = Kv;
     }
-
+ BOOST_MATH_INSTRUMENT_VARIABLE(*I);
+ BOOST_MATH_INSTRUMENT_VARIABLE(*K);
     return 0;
 }
 

Modified: sandbox/math_toolkit/boost/math/special_functions/ellint_1.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/ellint_1.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/ellint_1.hpp 2007-10-02 05:30:41 EDT (Tue, 02 Oct 2007)
@@ -80,12 +80,12 @@
        // so rewritten to use fmod instead:
        //
        BOOST_MATH_INSTRUMENT_CODE("pi/2 = " << constants::pi<T>() / 2);
- T rphi = fmod(phi, constants::pi<T>() / 2);
+ T rphi = boost::math::tools::fmod_workaround(phi, constants::pi<T>() / 2);
        BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
        T m = 2 * (phi - rphi) / constants::pi<T>();
        BOOST_MATH_INSTRUMENT_VARIABLE(m);
        int s = 1;
- if(fmod(m, T(2)) > 0.5)
+ if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
        {
           m += 1;
           s = -1;

Modified: sandbox/math_toolkit/boost/math/special_functions/ellint_2.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/ellint_2.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/ellint_2.hpp 2007-10-02 05:30:41 EDT (Tue, 02 Oct 2007)
@@ -69,10 +69,10 @@
        // but that fails if T has more digits than a long long,
        // so rewritten to use fmod instead:
        //
- T rphi = fmod(phi, constants::pi<T>() / 2);
+ T rphi = boost::math::tools::fmod_workaround(phi, constants::pi<T>() / 2);
        T m = 2 * (phi - rphi) / constants::pi<T>();
        int s = 1;
- if(fmod(m, T(2)) > 0.5)
+ if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
        {
           m += 1;
           s = -1;

Modified: sandbox/math_toolkit/boost/math/special_functions/ellint_3.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/ellint_3.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/ellint_3.hpp 2007-10-02 05:30:41 EDT (Tue, 02 Oct 2007)
@@ -180,10 +180,10 @@
     }
     else
     {
- T rphi = fmod(fabs(phi), constants::pi<T>() / 2);
+ T rphi = boost::math::tools::fmod_workaround(fabs(phi), constants::pi<T>() / 2);
        T m = 2 * (fabs(phi) - rphi) / constants::pi<T>();
        int sign = 1;
- if(fmod(m, T(2)) > 0.5)
+ if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
        {
           m += 1;
           sign = -1;

Modified: sandbox/math_toolkit/boost/math/special_functions/spherical_harmonic.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/spherical_harmonic.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/spherical_harmonic.hpp 2007-10-02 05:30:41 EDT (Tue, 02 Oct 2007)
@@ -56,7 +56,7 @@
    if(m&1)
    {
       // Check phase if theta is outside [0, PI]:
- T mod = fmod(theta, 2 * constants::pi<T>());
+ T mod = boost::math::tools::fmod_workaround(theta, 2 * constants::pi<T>());
       if(mod < 0)
          mod += 2 * constants::pi<T>();
       if(mod > constants::pi<T>())
@@ -83,7 +83,7 @@
    if(m&1)
    {
       // Check phase if theta is outside [0, PI]:
- T mod = fmod(theta, 2 * constants::pi<T>());
+ T mod = boost::math::tools::fmod_workaround(theta, 2 * constants::pi<T>());
       if(mod < 0)
          mod += 2 * constants::pi<T>();
       if(mod > constants::pi<T>())
@@ -114,7 +114,7 @@
    if(m&1)
    {
       // Check phase if theta is outside [0, PI]:
- U mod = fmod(theta, 2 * constants::pi<U>());
+ U mod = boost::math::tools::fmod_workaround(theta, 2 * constants::pi<U>());
       if(mod < 0)
          mod += 2 * constants::pi<U>();
       if(mod > constants::pi<U>())

Modified: sandbox/math_toolkit/boost/math/tools/config.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/tools/config.hpp (original)
+++ sandbox/math_toolkit/boost/math/tools/config.hpp 2007-10-02 05:30:41 EDT (Tue, 02 Oct 2007)
@@ -6,6 +6,10 @@
 #include <boost/detail/workaround.hpp>
 #include <algorithm> // for min and max
 #include <cmath>
+#include <climits>
+#if (defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__))
+# include <math.h>
+#endif
 
 #include <boost/math/tools/user.hpp>
 
@@ -101,7 +105,24 @@
 {
    return (std::max)((std::max)(a, b), (std::max)(c, d));
 }
-
+//
+// We call this short forwarding function so that we can work around a bug
+// on Darwin that causes std::fmod to return a NaN. The test case is:
+// std::fmod(1185.0L, 1.5L);
+//
+template <class T>
+inline T fmod_workaround(T a, T b)
+{
+ BOOST_MATH_STD_USING
+ return fmod(a, b);
+}
+#if (defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__)) && ((LDBL_MANT_DIG == 106) || (__LDBL_MANT_DIG__ == 106))
+template <>
+inline long double fmod_workaround(long double a, long double b)
+{
+ return ::fmodl(a, b);
+}
+#endif
 } // namespace tools
 }} // namespace boost namespace math
 


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