Boost math library on PPC64 has thousands of errors

[Tony Reix <tony.reix@…>]

Hi,

I'm trying to port Boost 1.55 on PowerPC64 BE and LE. Thousands of errors appear with math library.

Investigation one test (bessel_zeros_example_1) which loops, I see by means of gdb the following:

boost::math::cyl_bessel_j_zero< double, ..... boost::math::detail::cyl_bessel_j_zero_imp<float128, .....

This "float128" seems erroneous. However, looking at boost/math/special_functions/bessel.hpp code, I see that this type is chosen at compilation time, by means of:

typedef typename policies::evaluation<result_type, Policy>::type value_type;

and it is a nightmare to understand why the compiler did choose float128. Unless gdb is showing nuts ?

Exact test is: bin.v2/libs/math/example/bessel_zeros_example_1.test/gcc-4.8.2/debug/bessel_zeros_example_1

Comparing execution path of this test code in PPC64/BE to Intel64, I see some floor() value computed differently on the 2 systems, and leading the test program to enter a wrong path and stay stucked in an infinite loop.

There is some fundamental error, hidden by compiler which compiles templates and gives me no hint.

I need help to determine what is wrong with Boost math on PPC64.

Tony

[Tony Reix <tony.reix@…>]

The issue may be due to the value of LDBL_EPSILON on Linux on PPC64.

Lets use: gcc -E -dM file.cpp which prints the constants. Then:

gcc -E -dM file.cpp | grep DBL_EPSILON

gives the value of constants ending like this.

AMD64: #define LDBL_EPSILON 1.08420217248550443401e-19L #define DBL_EPSILON double(2.22044604925031308085e-16L)

AIX 61S PPC64: #define LDBL_EPSILON 2.2204460492503131e-16L #define DBL_EPSILON double(2.2204460492503131e-16L)

Ubuntu/PPC64-LE: #define LDBL_EPSILON 4.94065645841246544176568792868221e-324L #define DBL_EPSILON double(2.22044604925031308084726333618164e-16L)

Fedora/PPC64-BE: #define LDBL_EPSILON 4.94065645841246544176568792868221e-324L #define DBL_EPSILON double(2.22044604925031308084726333618164e-16L)

IBM said: LDBL_EPSILON is the smallest nonnegative number such that 1 + LDBL_EPSILON is representable. Given that PowerPC (on Linux) uses the "double + double" representation, where the value of the long double is determined as the sum of two doubles, LDBL_EPSILON ends up very small, namely the smallest representable "double". (Because you can then represent 1 + LDBL_EPSILON by setting the first double of the pair to 1, and the second to LDBL_EPSILON.)

The unusually small epsilon is just one aspect of the fact that the PowerPC long double format simply is not IEEE compliant. (It may well be that the fact that long double is not IEEE compliant gives trouble to Boost ... Maybe there is a way to tell Boost to avoid using long double?)

On Intel, long double is 80-bit IEEE, which reflects the epsilon value you see there. On AIX, they're apparently using "long double" == "double".

I think that something must be changed in Boost. However, I'm unable to find what and where.

[John Maddock]

The simplest solution is to define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS which will disable all long double support. That define would go in boost/math/tools/config.hpp near the top.

The problem with a "double double" type is that not only is it not IEEE conforming, but it's basically impossible to reason about it's behavior in any kind of sensible way. As none of us have any access to a system that uses such a type we've never tried to support it: if you wanted to try and patch Boost to support that type I have some ideas where to start, but it wouldn't be a small job I suspect :-(

[Tony Reix <tony.reix@…>]

Replying to johnmaddock:

Thanks for your comments !
I'm now recompiling/testing Boost with such a change:

# g++ -E -dM boost/math/tools/config.hpp | grep LONG_DOUBLE
#define _GLIBCXX_LONG_DOUBLE_COMPAT 1
#define SIZEOF_LONG_DOUBLE 16
#define BOOST_MATH_NO_NATIVE_LONG_DOUBLE_FP_CLASSIFY
#define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
#define LONG_DOUBLE_128 1
#define LONG_DOUBLE_MATH_OPTIONAL 1

I should have the answer in some hours.

About your suggestion to provide me with some ideas on how to support "double double", I'm OK to work with you. I may spend days and maybe some weeks on this. I have access to both PPC64 BIG/LITTLE Endian machines, installed either with Ubuntu or Fedora (for now).
I'm new to Boost.
Just let me know your suggestions, plus some explanations and pages to read, and I'll experiment them.

:)

[John Maddock]

Actually I see there's at least a partial patch for "double double" support in there already, it's in boost/math/tools/precision.hpp:

#if (defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__)) && ((LDBL_MANT_DIG == 106) || (__LDBL_MANT_DIG__ == 106))
template <>
inline long double epsilon<long double>(const mpl::true_& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE(long double))
{
   // numeric_limits on Darwin tells lies here.
   // This static assert fails for some unknown reason, so
   // disabled for now...
   // BOOST_STATIC_ASSERT(std::numeric_limits<long double>::digits == 106);
   return 2.4651903288156618919116517665087e-32L;
}
#endif

So that preprocessor branch need to be enabled for your platform as well. Then it's a case of seeing what still breaks.

Are you familiar with Git? If so the best thing would be for you to clone Boost.Math develop branch (​https://github.com/boostorg/math/tree/develop) and submit pull requests.

Cheers, John.

[Tony Reix <tony.reix@…>]

Replying to johnmaddock:

I've the result for the previous change (BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS): I now have MUCH less errors. However, I simply did a quick rebuild/retest from existing state and not a complete rebuilt/retest after a full clean . Anyway, I now see that this long double issue is the root cause for MANY of the errors I've seen.

About your second suggestion, I've removed the previous change (NO "long double") and added a defined(PPC) at beginning of the #if . However, I have yet no idea about how long is the mantissa... I've asked for some paper summarizing the choices made for Linux on PPC64 (because they are different from the ones made for AIX).

About github, I'm no more familiar, but I have a colleague who knows it very well and should be able to refresh/update my memory.

Last: are there some Boost documentation explaining the main files used by Boost for its basic configuration, like boost/math/tools/config.hpp ?

Thanks/Merci

Tony

[Tony Reix <tony.reix@…>]

#define LDBL_MAX 1.79769313486231580793728971405301e+308L
#define LDBL_MAX_EXP 1024
#define LDBL_HAS_INFINITY 1
#define LDBL_MIN 2.00416836000897277799610805135016e-292L
#define LDBL_HAS_QUIET_NAN 1
#define LDBL_HAS_DENORM 1
#define LDBL_EPSILON 4.94065645841246544176568792868221e-324L
#define LDBL_MANT_DIG 106
#define LDBL_MIN_EXP (-968)
#define LDBL_DENORM_MIN 4.94065645841246544176568792868221e-324L
#define LDBL_MIN_10_EXP (-291)
#define LDBL_DIG 31

[Tony Reix <tony.reix@…>]

Replying to johnmaddock:

Are you familiar with Git? If so the best thing would be for you to clone Boost.Math develop branch (​https://github.com/boostorg/math/tree/develop) and submit pull requests.

I have now a clone of Boost.Math develop branch on my PC.
What can I do with it ? I even don't know how to compile it.

Thx Tony

[Tony Reix <tony.reix@…>]

Replying to johnmaddock:

About your seconde proposal (I've removed the previous change (NO "long double") and I've added a defined(PPC) at beginning of the #if in boost/math/tools/precision.hpp), I still have errors, but much less than at first. And, more important, I no more have the tests looping infinitely and no core in math library.

4300 tests (whole Boost)
83 failures

(inc. 48 failure in boost/math library)

Only 7046 failures in the sub-tests (instead of 50,000)

That seems to be a good start !

Here are the failing math tests:

161 test_tr1_c_long_double.test/gcc-4.8.2/debug/link-static/test_tr1_c_long_double.output
201 complex_test.test/gcc-4.8.2/debug/link-static/complex_test.output
201 octonion_test.test/gcc-4.8.2/debug/link-static/octonion_test.output
201 quaternion_test.test/gcc-4.8.2/debug/link-static/quaternion_test.output
201 test_airy.test/gcc-4.8.2/debug/link-static/test_airy.output
201 test_bessel_airy_zeros.test/gcc-4.8.2/debug/link-static/test_bessel_airy_zeros.output
201 test_bessel_j.test/gcc-4.8.2/debug/link-static/test_bessel_j.output
201 test_bessel_k.test/gcc-4.8.2/debug/link-static/test_bessel_k.output
201 test_bessel_y.test/gcc-4.8.2/debug/link-static/test_bessel_y.output
201 test_beta.test/gcc-4.8.2/debug/link-static/test_beta.output
201 test_binomial_coeff.test/gcc-4.8.2/debug/link-static/test_binomial_coeff.output
201 test_binomial_real_concept1.test/gcc-4.8.2/debug/link-static/test_binomial_real_concept1.output
201 test_binomial_real_concept2.test/gcc-4.8.2/debug/link-static/test_binomial_real_concept2.output
201 test_binomial_real_concept4.test/gcc-4.8.2/debug/link-static/test_binomial_real_concept4.output
201 test_binomial_real_concept6.test/gcc-4.8.2/debug/link-static/test_binomial_real_concept6.output
201 test_constants.test/gcc-4.8.2/debug/link-static/test_constants.output
201 test_expint.test/gcc-4.8.2/debug/link-static/test_expint.output
201 test_factorials.test/gcc-4.8.2/debug/link-static/test_factorials.output
201 test_hankel.test/gcc-4.8.2/debug/link-static/test_hankel.output
201 test_hypergeometric_dist0.test/gcc-4.8.2/debug/link-static/test_hypergeometric_dist0.output
201 test_hypergeometric_dist1.test/gcc-4.8.2/debug/link-static/test_hypergeometric_dist1.output
201 test_hypergeometric_dist2.test/gcc-4.8.2/debug/link-static/test_hypergeometric_dist2.output
201 test_hypergeometric_dist3.test/gcc-4.8.2/debug/link-static/test_hypergeometric_dist3.output
201 test_hypergeometric_dist4.test/gcc-4.8.2/debug/link-static/test_hypergeometric_dist4.output
201 test_hypergeometric_dist5.test/gcc-4.8.2/debug/link-static/test_hypergeometric_dist5.output
201 test_ibeta_double.test/gcc-4.8.2/debug/link-static/test_ibeta_double.output
201 test_ibeta_inv_ab_double.test/gcc-4.8.2/debug/link-static/test_ibeta_inv_ab_double.output
201 test_ibeta_inv_double.test/gcc-4.8.2/debug/link-static/test_ibeta_inv_double.output
201 test_ibeta_inv_long_double.test/gcc-4.8.2/debug/link-static/test_ibeta_inv_long_double.output
201 test_ibeta_real_concept3.test/gcc-4.8.2/debug/link-static/test_ibeta_real_concept3.output
201 test_igamma_inva_real_concept.test/gcc-4.8.2/debug/link-static/test_igamma_inva_real_concept.output
201 test_igamma_inv_real_concept.test/gcc-4.8.2/debug/link-static/test_igamma_inv_real_concept.output
201 test_jacobi.test/gcc-4.8.2/debug/link-static/test_jacobi.output
201 test_legendre.test/gcc-4.8.2/debug/link-static/test_legendre.output
201 test_nc_beta_double.test/gcc-4.8.2/debug/link-static/test_nc_beta_double.output
201 test_nc_beta_long_double.test/gcc-4.8.2/debug/link-static/test_nc_beta_long_double.output
201 test_nc_chi_squared_long_double.test/gcc-4.8.2/debug/link-static/test_nc_chi_squared_long_double.output
201 test_nc_chi_squared_real_concept.test/gcc-4.8.2/debug/link-static/test_nc_chi_squared_real_concept.output
201 test_nc_t_long_double.test/gcc-4.8.2/debug/link-static/test_nc_t_long_double.output
201 test_next.test/gcc-4.8.2/debug/link-static/test_next.output
201 test_owens_t.test/gcc-4.8.2/debug/link-static/test_owens_t.output
201 test_rationals.test/gcc-4.8.2/debug/link-static/test_rationals.output
201 test_round.test/gcc-4.8.2/debug/link-static/test_round.output
201 test_skew_normal.test/gcc-4.8.2/debug/link-static/test_skew_normal.output
201 test_spherical_harmonic.test/gcc-4.8.2/debug/link-static/test_spherical_harmonic.output
201 test_students_t.test/gcc-4.8.2/debug/link-static/test_students_t.output
201 test_tr1_long_double.test/gcc-4.8.2/debug/link-static/test_tr1_long_double.output
201 test_zeta.test/gcc-4.8.2/debug/link-static/test_zeta.output

Let me know if you have a recommendation about which one I should look at at first.

Thx ! Tony

[Tony Reix <tony.reix@…>]

Discussion is continuing by exchanging emails.

IBM implements "long double" as 2 doubles (128 bits). It seems that Boost does not (or not entirely) manage this case.

John and I are working together to determine how to have "long double" correctly managed on Linux On Power.

[Tony Reix <tony.reix@…>]

Hi John Are you back ? so that we continue experimenting and fixing this issue ? Thx Tony