[boost] [math toolkit] Quick and Late Review

From: Jeff Garland (jeff_at_[hidden])
Date: 2007-04-28 16:37:47

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"