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