[math toolkit] Quick and Late Review
Executive Summary: Accept this library :-) I didn't have much time for this review (and sorry it's late), but I'll add my little bit to the other reviews. First, I'm not an expert in numerics, but I've been around enough to appreciate the importance and difficulty of doing this well -- which I believe this library clearly does. My last serious brush writing this sort of code was in college in a computer-math class -- and to date myself, we mostly used Fortran 4 at the time. Since I'm probably the least qualified to evaluate the math, my review will mostly focus on 'boost/library-related' mechanics. My general impression is that this library is, quite simply, an incredible piece of work. In general, the documentation is wonderfully written, illustrated beautifully, and stoked with references to all the details one could want. The effort applied to the test cases in this library is clearly immense -- and I'm sure the effort in making them all work was immense as well. I had a few failures when I ran the tests on a couple different machines (see below), but I'm certain those can be cleaned up. My one real issue with the library construction relates to error handling. I think it could be improved as I outline below. Here's what I did for my review. Downloaded, compiled and ran tests on 2 machines, read thru parts of the documentation, browsed over the tests/implementation. 1) Error handling This is one part of the implementation that I have issues with. Given that there the library supports user replaceable error handling, they can do pretty much anything they want....so, in the end, the user can do what they want. Anyway, this is a multi-part thing for me, so lets go into details: a) Macros If I understand correctly, to get 'fully signaling' functions I have to write the following code: #define BOOST_MATH_THROW_ON_DOMAIN_ERROR #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR #define BOOST_MATH_THROW_ON_UNDERFLOW_ERROR #define BOOST_MATH_THROW_ON_DENORM_ERROR #define BOOST_MATH_THROW_ON_LOGIC_ERROR #include <boost/math/...whatever..> This isn't 'pretty', to say the least. Shouldn't there at least be a single macro to convert from 'NAN-based' exception handling to exceptions? Maybe: #define BOOST_MATH_THROW_ON_ALL_ERROR b) I suggest that exceptions on all be the default with no macro required. Of course, numerics experts might have a different feeling here, but many c++ programmers now expect exceptions by default when something goes wrong. I believe it simplifies library usage and those that really want to use C style programming with all kinds of 'errno' checks are free to set the macro (that's one macro) back the other way. c) Use of 'errno'. I see that this actually is how errors are transmitted in many cases for the non-signaling version of the error functions. It's never mentioned in the docs -- it should be. And what the errors are set to needs to be documented. And probably an example. d) Handling 'Not Implemented' as 'domain error' Handling of Functions that are Not Implemented Functions that are not implemented for any reason, usually because they are not defined, (or we were unable to find a definition), are handled as domain errors. I guess I'm wondering why the not-implemented cases aren't a compilation errors instead of runtime errors instead? c) Uses of logic_error I think of logic_error as primarily for this sort of code: switch (something) { //... default: //can't ever possibly happen, really... throw std::logic_error(...). } But in the implementation I see code like: // detail/gamma_inva.hpp if(max_iter >= 200) tools::logic_error<T>(BOOST_CURRENT_FUNCTION, "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first); return (r.first + r.second) / 2; This seems like a unique numeric error for some types of algorithms. So, I'd suggest this isn't so much a logic error as 'convergence_error' or something. BTW, this is one of those cases where you can't check the result in the non-signaling case -- you need to check errno to see that there's an error which is not documented. d) documentation -- figuring out how to set it up Let's start with the docs on the configuration macro. There is says: "BOOST_MATH_THROW_ON_DOMAIN_ERROR If the macro BOOST_MATH_THROW_ON_DOMAIN_ERROR is define when building the library, " Of course, there's no library built, so I'm already a bit confused. Then continuing, later the paragraph says: defining this macro is generally recommended to get helpful error messages, which leaves me wondering what the default is. It is clarified later in the Error Handling Example where it says: But, by default, none of these exceptions will be raised, and the most appropriate value, usually a NaN, or zero, will be returned. Also in the macro config section it says: • BOOST_MATH_THROW_ON_POLE_ERROR By default, pole error is the same as domain_error. Meaning? So, my first issue is that I need to look at too many different places to really figure out how the error handling works. One section up front that deals with it and how to configure it would be good. Also, mention of the fact that errno needs to be checked in some cases for non-signaling error setup. 2) NTL patch (docs p249) In order to do so you will need to apply the following patch to NTL: libs/math/tools/ntl.diff. This patch adds trivial converting constructors to NTL::RR and NTL::quad_float, and forces conversions to RR to proceed via long double rather than double. The latter change Sounds kinda inconvenient. Isn't there a way this could be done without actually changing the NTL library? 3) 64bit Linux, gcc 4.0 -- 2 Test failures (see below) a) special functions b) remez 4) 32bit Linux, gcc 4.1 -- 1 Test failure (see below) a) test_beta_dist 5) Documentation nits (pdf) p 41 -- ON_OVERFLOW is listed twice BOOST_MATH_THROW_ON_DOMAIN_ERROR BOOST_MATH_THROW_ON_OVERFLOW_ERROR BOOST_MATH_THROW_ON_OVERFLOW_ERROR BOOST_MATH_THROW_ON_UNDERFLOW_ERROR BOOST_MATH_THROW_ON_DENORM_ERROR BOOST_MATH_THROW_ON_LOGIC_ERROR p 261 Formula after "With individual coefficients defined in closed form by:" doesn't display well... p 5 sentence ends with an apparent blown reference: 'See' BOOST_MATH_THROW_ON_DOMAIN_ERROR.....defining this macro is generally recommended to get helpful error messages, but using a try and catch blocks are also recommended. See Ok, that's it -- nice job to all! Jeff **************************************************************** * Test failure output * **************************************************************** MkDir1 ../bin.v2 MkDir1 ../bin.v2/test MkDir1 ../bin.v2/test/special_functions_test.test MkDir1 ../bin.v2/test/special_functions_test.test/gcc-4.0 MkDir1 ../bin.v2/test/special_functions_test.test/gcc-4.0/debug gcc.compile.c++ ../bin.v2/test/special_functions_test.test/gcc-4.0/debug/special_functions_test.o gcc.link ../bin.v2/test/special_functions_test.test/gcc-4.0/debug/special_functions_test testing.capture-output ../bin.v2/test/special_functions_test.test/gcc-4.0/debug/special_functions_test.run ====== BEGIN OUTPUT ====== Results of special functions test. (C) Copyright Hubert Holin 2003-2005. 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_0.txt) Running 21 test cases... Testing atanh in the real domain for float. Testing atanh in the real domain for double. Testing atanh in the real domain for long double. Testing asinh in the real domain for float. asinh_test.hpp(56): error in "asinh_test<f>": check ::std::less_equal<T>()( asinh_error_evaluator(x), static_cast<T>(4) ) failed for ( 44.7235527, 4 ) asinh_test.hpp(56): error in "asinh_test<f>": check ::std::less_equal<T>()( asinh_error_evaluator(x), static_cast<T>(4) ) failed for ( 44.7235527, 4 ) Testing asinh in the real domain for double. Testing asinh in the real domain for long double. Testing acosh in the real domain for float. Testing acosh in the real domain for double. Testing acosh in the real domain for long double. Testing sinc_pi in the real domain for float. Testing sinc_pi in the real domain for double. Testing sinc_pi in the real domain for long double. Testing sinhc_pi in the real domain for float. Testing sinhc_pi in the real domain for double. Testing sinhc_pi in the real domain for long double. Testing sinc_pi in the complex domain for float. Testing sinc_pi in the complex domain for double. Testing sinc_pi in the complex domain for long double. Testing sinhc_pi in the complex domain for float. Testing sinhc_pi in the complex domain for double. Testing sinhc_pi in the complex domain for long double. *** 2 failures detected in test suite "Master Test Suite" EXIT STATUS: 201 ********************************************************************* b) remez ../bin.v2/test/test_rationals.test/gcc-4.0/debug/test_rationals.test MkDir1 ../bin.v2/test/test_remez.test MkDir1 ../bin.v2/test/test_remez.test/gcc-4.0 MkDir1 ../bin.v2/test/test_remez.test/gcc-4.0/debug ...on 300th target... gcc.compile.c++ ../bin.v2/test/test_remez.test/gcc-4.0/debug/test_remez.o gcc.link ../bin.v2/test/test_remez.test/gcc-4.0/debug/test_remez testing.capture-output ../bin.v2/test/test_remez.test/gcc-4.0/debug/test_remez.run ====== BEGIN OUTPUT ====== Running 1 test case... Testing expm1 approximation, pinned to origin, abolute error, 6 term polynomial Interpolation Error: 7.21585e-06 4.67237e-06 1.10153e-05 0.655374 3.94117e-06 5.87427e-06 0.334035 3.63169e-06 4.84225e-06 0.0741202 3.58908e-06 3.73314e-06 0.018902 3.58717e-06 3.58806e-06 0.00154359 3.58716e-06 3.58716e-06 1.01598e-05 3.58716e-06 3.58716e-06 2.31932e-06 ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing expm1 approximation, pinned to origin, relative error, 6 term polynomial Interpolation Error: 3.58716e-06 5.05911e-05 0.000105139 7.56861e+293 8.99588e-06 1.66451e-05 1 7.52236e-06 1.05985e-05 0.976718 5.96499e-06 1.07151e-05 0.988659 5.14076e-06 8.48977e-06 0.995751 4.74167e-06 5.85447e-06 0.998221 4.56107e-06 4.75215e-06 0.999149 ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing exp approximation, not pinned to origin, abolute error, 6 term polynomial Interpolation Error: 3.58716e-06 -3.20174e-06 3.21961e-06 0.0694688 -3.21087e-06 3.21088e-06 0.000783811 -3.21087e-06 3.21087e-06 4.70208e-06 -3.21087e-06 3.21087e-06 2.39339e-06 -3.21087e-06 3.21087e-06 2.29707e-06 -3.21087e-06 3.21087e-06 8.69113e-07 -3.21087e-06 3.21087e-06 2.03257e-06 ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing exp approximation, not pinned to origin, relative error, 6 term polynomial Interpolation Error: 3.58716e-06 -2.64244e-06 3.45878e-06 0.849443 -3.00933e-06 3.0251e-06 0.114149 -3.0165e-06 3.0165e-06 0.00224181 -3.0165e-06 3.0165e-06 1.27282e-06 -3.0165e-06 3.0165e-06 6.6682e-06 -3.0165e-06 3.0165e-06 6.36906e-06 -3.0165e-06 3.0165e-06 3.28809e-06 ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing cos approximation, not pinned to origin, abolute error, 5 term polynomial Interpolation Error: 3.58716e-06 -4.18771e-05 4.18777e-05 0.0557121 -4.18775e-05 4.18775e-05 0.00135094 -4.18775e-05 4.18775e-05 6.44868e-06 -4.18775e-05 4.18775e-05 2.29391e-06 -4.18775e-05 4.18775e-05 7.84055e-07 -4.18775e-05 4.18775e-05 6.24391e-07 -4.18775e-05 4.18775e-05 6.24397e-07 ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing cos approximation, not pinned to origin, relative error, 5 term polynomial -5.53899e-05 5.6258e-05 0.669216 -5.59351e-05 5.59354e-05 0.265433 -5.59352e-05 5.59352e-05 0.241548 -5.59352e-05 5.59352e-05 1.71512e-06 -5.59352e-05 5.59352e-05 2.22885e-06 -5.59352e-05 5.59352e-05 1.83167e-06 -5.59352e-05 5.59352e-05 2.00456e-06 ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing sin approximation, pinned to origin, abolute error, 4 term polynomial 7.02899e-06 2.21518e-05 0.721841 1.57522e-05 1.6717e-05 0.0987357 1.60752e-05 1.60789e-05 0.00438332 1.60764e-05 1.60764e-05 3.04903e-05 1.60764e-05 1.60764e-05 8.6383e-07 1.60764e-05 1.60764e-05 1.0611e-06 1.60764e-05 1.60764e-05 1.0611e-06 ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing sin approximation, pinned to origin, relative error, 4 term polynomial 6.26553e-05 6.26949e-05 0.00949821 6.26764e-05 6.26764e-05 0.000194848 6.26764e-05 6.26764e-05 1.53971e-06 6.26764e-05 6.26764e-05 1.02629e-06 6.26764e-05 6.26764e-05 5.10901e-07 6.26764e-05 6.26764e-05 1.52291e-06 6.26764e-05 6.26764e-05 6.96439e-07 ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing expm1 approximation, pinned to origin, abolute error, 3+3 term rational Interpolation Error: 8.61648e-07 -9.83816e-07 1.86722e-06 0.617146 -5.72238e-07 6.18375e-07 0.204662 -5.33136e-07 5.34022e-07 0.0162229 -5.32772e-07 5.32912e-07 8.61675e-05 -5.32584e-07 5.32656e-07 4.1111e-05 -5.32486e-07 5.32524e-07 2.06807e-05 -5.32436e-07 5.32455e-07 1.04893e-05 ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing expm1 approximation, pinned to origin, relative error, 3+3 term rational Interpolation Error: 5.32455e-07 5.47611e-06 1.13109e-05 285666 Check failed in file ../../../boost/numeric/ublas/lu.hpp at line 272: detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2) unknown location(0): fatal error in "test_main_caller( argc, argv )": std::logic_error: internal logic *** 1 failure detected in test suite "Test Program" EXIT STATUS: 201 ***************************************************************** **passed** ../bin.v2/test/test_beta.test/gcc-4.1/debug/test_beta.test MkDir1 ../bin.v2/test/test_beta_dist.test MkDir1 ../bin.v2/test/test_beta_dist.test/gcc-4.1 MkDir1 ../bin.v2/test/test_beta_dist.test/gcc-4.1/debug gcc.compile.c++ ../bin.v2/test/test_beta_dist.test/gcc-4.1/debug/test_beta_dist.o gcc.link ../bin.v2/test/test_beta_dist.test/gcc-4.1/debug/test_beta_dist testing.capture-output ../bin.v2/test/test_beta_dist.test/gcc-4.1/debug/test_beta_dist.run ====== BEGIN OUTPUT ====== Running 1 test case... BOOST_MATH_THROW_ON_DOMAIN_ERROR is defined to throw on domain error. test_beta_dist.cpp(538): error in "test_main_caller( argc, argv )": check beta_distribution<double>::estimate_alpha(mean(mybeta22), variance(mybeta22)) == mybeta22.alpha() failed [1.9999999999999998 != 2] test_beta_dist.cpp(539): error in "test_main_caller( argc, argv )": check beta_distribution<double>::estimate_beta(mean(mybeta22), variance(mybeta22)) == mybeta22.beta() failed [1.9999999999999998 != 2] numeric_limits<real_concept>::is_specialized 0 numeric_limits<real_concept>::digits 0 numeric_limits<real_concept>::digits10 0 numeric_limits<real_concept>::epsilon 0 Boost::math::tools::epsilon = 1.19209e-07 std::numeric_limits::epsilon = 1.19209e-07 epsilon = 1.19209e-07, Tolerance = 0.0119209%. Boost::math::tools::epsilon = 2.22045e-16 std::numeric_limits::epsilon = 2.22045e-16 epsilon = 2.22045e-16, Tolerance = 2.22045e-11%. Boost::math::tools::epsilon = 1.0842e-19 std::numeric_limits::epsilon = 1.0842e-19 epsilon = 2.22045e-16, Tolerance = 2.22045e-11%. Boost::math::tools::epsilon = 1.0842e-19 std::numeric_limits::epsilon = 0 epsilon = 2.22045e-16, Tolerance = 2.22045e-11%. *** 2 failures detected in test suite "Test Program"
Paris (U.E.), le 29/04/2007 Bonsoir (Sorry, I did not have time to write the second part of my review, and now it is too late... at least the most important part went in). In article <4633B09B.50903@crystalclearsoftware.com>, Jeff Garland <jeff@crystalclearsoftware.com> wrote: [SNIP]
3) 64bit Linux, gcc 4.0 -- 2 Test failures (see below) a) special functions
[SNIP]
Jeff
**************************************************************** * Test failure output * **************************************************************** MkDir1 ../bin.v2 MkDir1 ../bin.v2/test MkDir1 ../bin.v2/test/special functions test.test MkDir1 ../bin.v2/test/special functions test.test/gcc-4.0 MkDir1 ../bin.v2/test/special functions test.test/gcc-4.0/debug gcc.compile.c++ ../bin.v2/test/special functions test.test/gcc-4.0/debug/special functions test.o gcc.link ../bin.v2/test/special functions test.test/gcc-4.0/debug/special functions test testing.capture-output ../bin.v2/test/special functions test.test/gcc-4.0/debug/special functions test.run ====== BEGIN OUTPUT ====== Results of special functions test.
(C) Copyright Hubert Holin 2003-2005. 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 0.txt)
Running 21 test cases... Testing atanh in the real domain for float. Testing atanh in the real domain for double. Testing atanh in the real domain for long double. Testing asinh in the real domain for float. asinh test.hpp(56): error in "asinh test<f>": check ::std::less equal<T>()( asinh error evaluator(x), static cast<T>(4) ) failed for ( 44.7235527, 4 ) asinh test.hpp(56): error in "asinh test<f>": check ::std::less equal<T>()( asinh error evaluator(x), static cast<T>(4) ) failed for ( 44.7235527, 4 ) Testing asinh in the real domain for double. Testing asinh in the real domain for long double. Testing acosh in the real domain for float. Testing acosh in the real domain for double. Testing acosh in the real domain for long double. Testing sinc pi in the real domain for float. Testing sinc pi in the real domain for double. Testing sinc pi in the real domain for long double. Testing sinhc pi in the real domain for float. Testing sinhc pi in the real domain for double. Testing sinhc pi in the real domain for long double. Testing sinc pi in the complex domain for float. Testing sinc pi in the complex domain for double. Testing sinc pi in the complex domain for long double. Testing sinhc pi in the complex domain for float. Testing sinhc pi in the complex domain for double. Testing sinhc pi in the complex domain for long double.
*** 2 failures detected in test suite "Master Test Suite"
[SNIP] OK, these errors I know and hate. I am the culprit here. The problem, as far as I could tell, boils down to a QOI issue on 64 bits platforms, using Gcc: the implementation of transcendental functions for 64 bits behaves just like if it were a 32 bits value being passed. As far as I know, this is legal (if horrendous) behavior. I still also find the idea of not being able to use 64 bits floating point variables on platforms that purport to support them (my platform of choice included) abhorrent. Perhaps the only clean way out would be to re-implement these basic functions ourselves, but that, of course, would be the source of other problems (overriding part of the user's standard library). Merci Hubert Holin
Hubert Holin wrote:
*** 2 failures detected in test suite "Master Test Suite"
[SNIP]
OK, these errors I know and hate. I am the culprit here. The problem, as far as I could tell, boils down to a QOI issue on 64 bits platforms, using Gcc: the implementation of transcendental functions for 64 bits behaves just like if it were a 32 bits value being passed. As far as I know, this is legal (if horrendous) behavior. I still also find the idea of not being able to use 64 bits floating point variables on platforms that purport to support them (my platform of choice included) abhorrent. Perhaps the only clean way out would be to re-implement these basic functions ourselves, but that, of course, would be the source of other problems (overriding part of the user's standard library).
Have you submitted a gcc bug report? If you re-implement them, the thing to do would be to contribute them as a patch to the gcc c++ std lib. Jeff
Jeff Garland wrote:
1) Error handling
This is one part of the implementation that I have issues with. Given that there the library supports user replaceable error handling, they can do pretty much anything they want....so, in the end, the user can do what they want.
Anyway, this is a multi-part thing for me, so lets go into details:
a) Macros
If I understand correctly, to get 'fully signaling' functions I have to write the following code:
#define BOOST_MATH_THROW_ON_DOMAIN_ERROR #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR #define BOOST_MATH_THROW_ON_UNDERFLOW_ERROR #define BOOST_MATH_THROW_ON_DENORM_ERROR #define BOOST_MATH_THROW_ON_LOGIC_ERROR #include <boost/math/...whatever..>
This isn't 'pretty', to say the least. Shouldn't there at least be a single macro to convert from 'NAN-based' exception handling to exceptions? Maybe:
#define BOOST_MATH_THROW_ON_ALL_ERROR
b) I suggest that exceptions on all be the default with no macro required. Of course, numerics experts might have a different feeling here, but many c++ programmers now expect exceptions by default when something goes wrong. I believe it simplifies library usage and those that really want to use C style programming with all kinds of 'errno' checks are free to set the macro (that's one macro) back the other way.
Hmmm, I'm not sure about this, underflow and denormalised results aren't necessarily errors that need an exception. There's another issue here: of all these error conditions domain errors are probably the only ones that we can guarentee to detect and signal with 100% certainty. I've tried really hard to ensure the others get detected as well, but it's almost certainly impossible to be 100% sure that every possible case is handled correctly. I must document this better as well BTW.
c) Use of 'errno'. I see that this actually is how errors are transmitted in many cases for the non-signaling version of the error functions. It's never mentioned in the docs -- it should be. And what the errors are set to needs to be documented. And probably an example.
Will do.
d) Handling 'Not Implemented' as 'domain error'
Handling of Functions that are Not Implemented
Functions that are not implemented for any reason, usually because they are not defined, (or we were unable to find a definition), are handled as domain errors.
I guess I'm wondering why the not-implemented cases aren't a compilation errors instead of runtime errors instead?
Hmmm, probably an unfortunate choice of words: there are some statistical properties for some distributions that aren't defined: for example the Cauchy distribution doesn't have a mean. You can still compile code that asks for the mean of the Cauchy, but you will get a domain error if it's called. This leaway is actually useful: I used it to create an "any_distribution" class that virtualises away the actual type of the distribution so that the distribution type can be set at runtime. Without the "every property compiles for every distribution" rule this would be very hard to do. In any case, it's a legitimate question to ask what the mean of any distribution is, you just can't always expect a finite answer :-)
c) Uses of logic_error
I think of logic_error as primarily for this sort of code: switch (something) { //... default: //can't ever possibly happen, really... throw std::logic_error(...). }
But in the implementation I see code like: // detail/gamma_inva.hpp
if(max_iter >= 200) tools::logic_error<T>(BOOST_CURRENT_FUNCTION, "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first); return (r.first + r.second) / 2;
This seems like a unique numeric error for some types of algorithms. So, I'd suggest this isn't so much a logic error as 'convergence_error' or something. BTW, this is one of those cases where you can't check the result in the non-signaling case -- you need to check errno to see that there's an error which is not documented.
I'll look at this again, but I believe that logic errors are always signalled as exceptions in the current code, that might be a bug :-) If the answer to the question is "we can't compute this" then what kind of error is it if not a logic_error? runtime_error I guess, but that seems too general? All suggestions welcomed.
d) documentation -- figuring out how to set it up
Let's start with the docs on the configuration macro. There is says:
"BOOST_MATH_THROW_ON_DOMAIN_ERROR If the macro BOOST_MATH_THROW_ON_DOMAIN_ERROR is define when building the library, "
Of course, there's no library built, so I'm already a bit confused. Then continuing, later the paragraph says:
defining this macro is generally recommended to get helpful error messages,
which leaves me wondering what the default is. It is clarified later in the Error Handling Example where it says:
But, by default, none of these exceptions will be raised, and the most appropriate value, usually a NaN, or zero, will be returned.
Also in the macro config section it says:
• BOOST_MATH_THROW_ON_POLE_ERROR By default, pole error is the same as domain_error.
Meaning?
So, my first issue is that I need to look at too many different places to really figure out how the error handling works. One section up front that deals with it and how to configure it would be good.
Also, mention of the fact that errno needs to be checked in some cases for non-signaling error setup.
Will do.
2) NTL patch (docs p249)
In order to do so you will need to apply the following patch to NTL: libs/math/tools/ntl.diff. This patch adds trivial converting constructors to NTL::RR and NTL::quad_float, and forces conversions to RR to proceed via long double rather than double. The latter change
Sounds kinda inconvenient. Isn't there a way this could be done without actually changing the NTL library?
I'd love too. Probably by writing our own wrapper around the library I guess, but it was quicker/easier to patch it. Unfortunately the library seems not to be actively maintained at present :-(
Running 21 test cases... Testing atanh in the real domain for float. Testing atanh in the real domain for double. Testing atanh in the real domain for long double. Testing asinh in the real domain for float. asinh_test.hpp(56): error in "asinh_test<f>": check ::std::less_equal<T>()( asinh_error_evaluator(x), static_cast<T>(4) ) failed for ( 44.7235527, 4 ) asinh_test.hpp(56): error in "asinh_test<f>": check ::std::less_equal<T>()( asinh_error_evaluator(x), static_cast<T>(4) ) failed for ( 44.7235527, 4 )
I'll leave these to Hubert to sort out :-)
b) remez ~~~~~~~~~~~~~~~~~~~~~~~~~ Testing expm1 approximation, pinned to origin, relative error, 3+3 term rational Interpolation Error: 5.32455e-07 5.47611e-06 1.13109e-05 285666 Check failed in file ../../../boost/numeric/ublas/lu.hpp at line 272: detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2) unknown location(0): fatal error in "test_main_caller( argc, argv )": std::logic_error: internal logic
*** 1 failure detected in test suite "Test Program"
It's a failure in ublas that occurs on some platforms and not others (to be fair the problem matrix being solved is very "stiff", maybe I'll just remove that part of the test for now).
====== BEGIN OUTPUT ====== Running 1 test case... BOOST_MATH_THROW_ON_DOMAIN_ERROR is defined to throw on domain error. test_beta_dist.cpp(538): error in "test_main_caller( argc, argv )": check beta_distribution<double>::estimate_alpha(mean(mybeta22), variance(mybeta22)) == mybeta22.alpha() failed [1.9999999999999998 != 2] test_beta_dist.cpp(539): error in "test_main_caller( argc, argv )": check beta_distribution<double>::estimate_beta(mean(mybeta22), variance(mybeta22)) == mybeta22.beta() failed [1.9999999999999998 != 2]
Looks like the tolerances are set too low, Paul? Many thanks for the positive review! John.
John Maddock wrote:
Jeff Garland wrote:
1) Error handling
This is one part of the implementation that I have issues with. Given that there the library supports user replaceable error handling, they can do pretty much anything they want....so, in the end, the user can do what they want.
Anyway, this is a multi-part thing for me, so lets go into details:
a) Macros
If I understand correctly, to get 'fully signaling' functions I have to write the following code:
#define BOOST_MATH_THROW_ON_DOMAIN_ERROR #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR #define BOOST_MATH_THROW_ON_UNDERFLOW_ERROR #define BOOST_MATH_THROW_ON_DENORM_ERROR #define BOOST_MATH_THROW_ON_LOGIC_ERROR #include <boost/math/...whatever..>
This isn't 'pretty', to say the least. Shouldn't there at least be a single macro to convert from 'NAN-based' exception handling to exceptions? Maybe:
#define BOOST_MATH_THROW_ON_ALL_ERROR
b) I suggest that exceptions on all be the default with no macro required. Of course, numerics experts might have a different feeling here, but many c++ programmers now expect exceptions by default when something goes wrong. I believe it simplifies library usage and those that really want to use C style programming with all kinds of 'errno' checks are free to set the macro (that's one macro) back the other way.
Hmmm, I'm not sure about this, underflow and denormalised results aren't necessarily errors that need an exception.
Then are they really 'errors'? With underflow I would imagine that it would only matter if the result was impacted. As for denormalised, I'm not sure why I'd need to know -- likely this is my ignorance of numeric processing. But still, maybe there's some standard set of policies here that could be grouped logically into 'signalling' and 'non-signalling' error policies?
There's another issue here: of all these error conditions domain errors are probably the only ones that we can guarentee to detect and signal with 100% certainty. I've tried really hard to ensure the others get detected as well, but it's almost certainly impossible to be 100% sure that every possible case is handled correctly. I must document this better as well BTW.
That's a QOI issue to be sure.
d) Handling 'Not Implemented' as 'domain error'
Handling of Functions that are Not Implemented
Functions that are not implemented for any reason, usually because they are not defined, (or we were unable to find a definition), are handled as domain errors.
I guess I'm wondering why the not-implemented cases aren't a compilation errors instead of runtime errors instead?
Hmmm, probably an unfortunate choice of words: there are some statistical properties for some distributions that aren't defined: for example the Cauchy distribution doesn't have a mean. You can still compile code that asks for the mean of the Cauchy, but you will get a domain error if it's called. This leaway is actually useful: I used it to create an "any_distribution" class that virtualises away the actual type of the distribution so that the distribution type can be set at runtime. Without the "every property compiles for every distribution" rule this would be very hard to do. In any case, it's a legitimate question to ask what the mean of any distribution is, you just can't always expect a finite answer :-)
I guess I'm still uncomfortable here. Cauchy doesn't have a mean and that isn't going to change at runtime -- it's fixed for all time. If I call the function I always get an error -- really I should never call the function. So, at a minimum, this should be a unique and different error from 'domain error' which tells me I've called a valid function with parameters that are 'out of range'. However, If a distribution doesn't define a particular function why not exclude it using a trait, enable-if, or other compile-time logic to avoid making the call in the first place? I think you can still achieve your wrapper without ever calling the function at runtime and then handling an error. I'm guessing what's really going on is there are multiple distribution concepts that need to be defined with different groups of valid operations.
2) NTL patch (docs p249)
In order to do so you will need to apply the following patch to NTL: libs/math/tools/ntl.diff. This patch adds trivial converting constructors to NTL::RR and NTL::quad_float, and forces conversions to RR to proceed via long double rather than double. The latter change
Sounds kinda inconvenient. Isn't there a way this could be done without actually changing the NTL library?
I'd love too. Probably by writing our own wrapper around the library I guess, but it was quicker/easier to patch it. Unfortunately the library seems not to be actively maintained at present :-(
Isn't NTL mostly a wrapper around gmp? Anyway, Arseny is building a BigInt this summer. One implementation strategy he's using it to provide one implementation as a wrapper of GMP. Might be fairly easy to do something similar with a floating point type given that you've defined the concepts well -- SoC 2008 ;-)
Many thanks for the positive review!
Sure, but the real thx go to you guys for the thousands of hours of work on this :-) Jeff
On 4/28/07, Jeff Garland <jeff@crystalclearsoftware.com> wrote:
Executive Summary: Accept this library :-)
1) Error handling
c) Use of 'errno'. I see that this actually is how errors are transmitted in many cases for the non-signaling version of the error functions. It's never mentioned in the docs -- it should be. And what the errors are set to needs to be documented. And probably an example.
If it is using errno, can I assume it is not thread-safe? Any other threading issues?
Gottlob Frege wrote:
If it is using errno, can I assume it is not thread-safe? Any other threading issues?
Darn, "thread safety" is another page of docs I've been meaning to write ! As Stefan has said already, errno should be thread safe already (otherwise none of the std lib math functions would be thread safe?). If you're still not happy, turn on reporting of errors via exceptions and errno gets left unused anyway. Other than that, the code is intended to be thread safe *for built in real-number types* : so float, double and long double are all thread safe. For non-built-in types - NTL::RR for example - initialisation of the various constants used in the implementation is potentially *not* thread safe. I don't like this at all, but it would be a signficant challenge to fix it. I believe that some compilers are now starting to offer the option of having static-constants initialised in a thread safe manner (Commeau, and maybe others?), if that's the case then the problem is solved. This is a topic of hot debate for the next C++ std revision BTW, so hopefully all compilers will be required to do the right thing here at some point. HTH, John.
participants (5)
-
Gottlob Frege -
Hubert Holin -
Jeff Garland -
John Maddock -
Stefan Seefeld