[lmi] use boost::math:expm1() and log1p()

[Vaclav Slavik]

Hi,

as promised, here's a patch to use expm1() and log1p() instead of C99
expm1l() and log1pl() (removal of expm1.{c,h} is not included in the
patch for simplicity). 

Vaclav

=== modified file 'Makefile.am'
--- Makefile.am 2009-03-26 16:47:29 +0000
+++ Makefile.am 2009-04-02 18:27:39 +0000
@@ -245,7 +245,6 @@
     dbnames.cpp \
     death_benefits.cpp \
     emit_ledger.cpp \
-    expm1.c \
     facets.cpp \
     fenv_guard.cpp \
     fenv_lmi.cpp \
@@ -641,8 +640,7 @@
 test_loads_SOURCES = \
   $(common_test_objects) \
   loads.cpp \
-  loads_test.cpp \
-  expm1.c
+  loads_test.cpp
 test_loads_CXXFLAGS = $(AM_CXXFLAGS)
 
 test_map_lookup_SOURCES = \

=== modified file 'config_ming323.hpp'
--- config_ming323.hpp  2008-12-27 02:57:00 +0000
+++ config_ming323.hpp  2009-04-02 18:25:52 +0000
@@ -43,18 +43,10 @@
 // Version numbers are in 'include/_mingw.h' here:
 //   http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/?cvsroot=src
 
-#if 308 <= LMI_MINGW_VERSION
-#   define LMI_COMPILER_PROVIDES_EXPM1L
-#endif // 308 <= LMI_MINGW_VERSION
-
 #if 200 <= LMI_MINGW_VERSION
 #   define LMI_COMPILER_PROVIDES_ISNAN
 #endif // 200 <= LMI_MINGW_VERSION
 
-#if 202 <= LMI_MINGW_VERSION
-#   define LMI_COMPILER_PROVIDES_LOG1PL
-#endif // 202 <= LMI_MINGW_VERSION
-
 #if 200 <= LMI_MINGW_VERSION
 #   define LMI_COMPILER_PROVIDES_RINT
 #endif // 200 <= LMI_MINGW_VERSION
@@ -67,5 +59,20 @@
 #   define LMI_COMPILER_PROVIDES_STRTOLD
 #endif // 204 <= LMI_MINGW_VERSION
 
+// Configure Boost library for this compiler in cases where it lacks
+// features autodetection.
+
+#if defined BOOST_MATH_EXPM1_INCLUDED || defined BOOST_MATH_LOG1P_INCLUDED
+#   error config.hpp must be included before boost::math headers.
+#endif
+
+#if 308 <= LMI_MINGW_VERSION
+#   define BOOST_HAS_EXP1M
+#endif // 308 <= LMI_MINGW_VERSION
+
+#if 202 <= LMI_MINGW_VERSION
+#   define BOOST_HAS_LOG1P
+#endif // 202 <= LMI_MINGW_VERSION
+
 #endif // config_ming323_hpp
 

=== modified file 'configure.ac'
--- configure.ac        2008-12-27 02:57:00 +0000
+++ configure.ac        2009-04-02 18:42:40 +0000
@@ -255,11 +255,9 @@
 AM_CONDITIONAL(LMI_WITH_CGI, [test "x$lmi_cgicc_option" == "xyes"])
 
 dnl existing code already uses these macros so continue to use them instead of
-dnl the usual HAVE_EXPM1/HAVE_LOG1P
+dnl the usual HAVE_STRTOF etc.
 AC_CHECK_FUNC(strtof, AC_DEFINE(LMI_COMPILER_PROVIDES_STRTOF, [1], [Define 
this if you have strtof() function]))
 AC_CHECK_FUNC(strtold, AC_DEFINE(LMI_COMPILER_PROVIDES_STRTOLD, [1], [Define 
this if you have strtold() function]))
-AC_CHECK_FUNC(expm1l, AC_DEFINE(LMI_COMPILER_PROVIDES_EXPM1L, [1], [Define 
this if you have expm1l() function]))
-AC_CHECK_FUNC(log1pl, AC_DEFINE(LMI_COMPILER_PROVIDES_LOG1PL, [1], [Define 
this if you have log1pl() function]))
 AC_CHECK_FUNC(rint, AC_DEFINE(LMI_COMPILER_PROVIDES_RINT, [1], [Define this if 
you have rint() function]))
 
 dnl === Library checks ===
@@ -342,6 +340,11 @@
     [AC_MSG_ERROR([Boost headers not found, $errmsg])]
 )
 
+AC_CHECK_HEADER([boost/math/special_functions/expm1.hpp],
+    [],
+    [AC_MSG_ERROR([Boost version with boost::math::expm1() required, $errmsg])]
+)
+
 dnl By default no particular toolkit suffix should be needed when linking
 dnl against boost libraries. Except for mingw where it seems to be
 dnl required. Because under mingw default installation does _not_ create

=== modified file 'interest_rates.cpp'
--- interest_rates.cpp  2009-03-08 20:14:44 +0000
+++ interest_rates.cpp  2009-04-02 18:36:15 +0000
@@ -104,7 +104,7 @@
 // For the annual-effective method, transformation from annual to
 // daily and back again by naive methods would lose considerable
 // precision even when the spread and fee are zero, because i is
-// small relative to (1 + i). That is why expm1l() and log1pl() are
+// small relative to (1 + i). That is why expm1() and log1p() are
 // used instead of pow().
 //
 // If both spread and fee are zero, then the net rate should exactly

=== modified file 'math_functors.hpp'
--- math_functors.hpp   2008-12-27 02:57:00 +0000
+++ math_functors.hpp   2009-04-02 17:20:48 +0000
@@ -33,19 +33,13 @@
 #   define BOOST_STATIC_ASSERT(deliberately_ignored) class IgNoRe
 #endif // Defined __BORLANDC__ .
 
+#include <boost/math/special_functions/expm1.hpp>
+#include <boost/math/special_functions/log1p.hpp>
+
 #include <algorithm>
-#include <cmath>
 #include <functional>
 #include <stdexcept>
 
-#if !defined LMI_COMPILER_PROVIDES_EXPM1L
-extern "C" long double expm1l(long double);
-#endif // !defined LMI_COMPILER_PROVIDES_EXPM1L
-
-#if !defined LMI_COMPILER_PROVIDES_LOG1PL
-extern "C" long double log1pl(long double);
-#endif // !defined LMI_COMPILER_PROVIDES_LOG1PL
-
 // TODO ?? Write functors here for other refactorable uses of
 // std::pow() found throughout the program.
 
@@ -113,7 +107,7 @@
 //
 // Implementation note: greater accuracy and speed are obtained by
 // applying the transformation
-//   (1+i)^n - 1 <-> expm1l(log1pl(i) * n)
+//   (1+i)^n - 1 <-> expm1(log1p(i) * n)
 // to naive power-based formulas.
 
 template<typename T, int n>
@@ -124,7 +118,7 @@
     T operator()(T const& i) const
         {
         static long double const reciprocal_n = 1.0L / n;
-        long double z = expm1l(log1pl(i) * reciprocal_n);
+        long double z = boost::math::expm1(boost::math::log1p(i) * 
reciprocal_n);
         return static_cast<T>(z);
         }
 };
@@ -148,7 +142,7 @@
     BOOST_STATIC_ASSERT(0 < n);
     T operator()(T const& i) const
         {
-        long double z = expm1l(log1pl(i) * n);
+        long double z = boost::math::expm1(boost::math::log1p(i) * n);
         return static_cast<T>(z);
         }
 };
@@ -173,7 +167,10 @@
     T operator()(T const& i) const
         {
         static long double const reciprocal_n = 1.0L / n;
-        long double z = -n * expm1l(log1pl(i) * -reciprocal_n);
+        long double z = -n * boost::math::expm1
+                             (  boost::math::log1p(i)
+                             *  -reciprocal_n
+                             );
         return static_cast<T>(z);
         }
 };
@@ -203,11 +200,11 @@
     T operator()(T const& i, T const& spread, T const& fee) const
         {
         static long double const reciprocal_n = 1.0L / n;
-        long double z = expm1l
+        long double z = boost::math::expm1
             (
-            n * log1pl
-                (   expm1l(reciprocal_n * log1pl(i))
-                -   expm1l(reciprocal_n * log1pl(spread))
+            n * boost::math::log1p
+                (   boost::math::expm1(reciprocal_n * boost::math::log1p(i))
+                -   boost::math::expm1(reciprocal_n * 
boost::math::log1p(spread))
                 -          reciprocal_n * fee
                 )
             );
@@ -264,7 +261,10 @@
         else
             {
             static long double const reciprocal_12 = 1.0L / 12;
-            long double monthly_q = -expm1l(log1pl(-q) * reciprocal_12);
+            long double monthly_q = -boost::math::expm1
+                                     (   boost::math::log1p(-q)
+                                     *   reciprocal_12
+                                     );
             if(1.0L == monthly_q)
                 {
                 throw std::logic_error("Monthly q equals unity.");

=== modified file 'math_functors_test.cpp'
--- math_functors_test.cpp      2008-12-27 02:57:00 +0000
+++ math_functors_test.cpp      2009-04-02 18:38:40 +0000
@@ -153,7 +153,7 @@
     std::cout
         << "\n  -0.004 upper 365 by various methods\n"
         << std::setprecision(20)
-        << "    long double precision, expm1l and log1pl\n      "
+        << "    long double precision, expm1 and log1p\n      "
         << net_i_from_gross<double,365>()(0.0, 0.004, 0.0) << '\n'
         << "    long double precision, pow\n      "
         << net_i_from_gross_naive<double,365>()(0.0, 0.004, 0.0) << '\n'
@@ -163,7 +163,7 @@
     fenv_precision(fe_dblprec);
     std::cout
         << std::setprecision(20)
-        << "    double precision, expm1l and log1pl\n      "
+        << "    double precision, expm1 and log1p\n      "
         << net_i_from_gross<double,365>      ()(0.0, 0.004, 0.0) << '\n'
         << "    double precision, pow\n      "
         << net_i_from_gross_naive<double,365>()(0.0, 0.004, 0.0) << '\n'
@@ -176,7 +176,7 @@
 // different implementations.
 
 // This implementation naively uses std::pow(); it is both slower and
-// less inaccurate than an alternative using expm1l() and log1pl().
+// less inaccurate than an alternative using expm1() and log1p().
 void mete0()
 {
     volatile double x;
@@ -199,7 +199,7 @@
 void assay_speed()
 {
     std::cout << "  Speed test: pow   \n    " << TimeAnAliquot(mete0) << '\n';
-    std::cout << "  Speed test: expm1l\n    " << TimeAnAliquot(mete1) << '\n';
+    std::cout << "  Speed test: expm1\n     " << TimeAnAliquot(mete1) << '\n';
 }
 
 int test_main(int, char*[])

=== modified file 'platform_dependent.hpp'
--- platform_dependent.hpp      2008-12-27 02:57:00 +0000
+++ platform_dependent.hpp      2009-04-02 18:37:57 +0000
@@ -48,9 +48,8 @@
 // that are useful but absent from the standard language: for example,
 // it is difficult to implement cgi-bin without putenv(). Some others:
 //   _wcsdup(), fileno(), strcasecmp(), strdup()
-// should be avoided in general, but are required by wx. Still others,
-// like expm1l(), are in C99 but not C++98; the way their prototypes
-// are provided for gcc varies by platform.
+// should be avoided in general, but are required by wx.; the way their
+// prototypes are provided for gcc varies by platform.
 
 #if defined __GNUC__ && defined __STRICT_ANSI__
 #   define LMI_GNUC_STRICT_ANSI
@@ -97,16 +96,6 @@
 // and therefore a prototype for getch() is instead provided by other
 // means, locally, wherever it's needed.
 
-// GNU/Linux (but not MinGW) requires including certain headers while
-// __STRICT_ANSI__ is not defined in order to get prototypes for
-// certain functions, for C++ with '-std=c++98':
-//   math.h:  expm1l() and log1pl()
-// Use the C instead of the C++ system header so that the present file
-// can be included in C as well as C++ translation units, which is
-// temporarily useful until 'expm1.c' can be removed.
-
-#include <math.h>
-
 // Although gcc may once have defined __STRICT_ANSI__ differently:
 //   http://gcc.gnu.org/bugzilla/show_bug.cgi?id=3199
 // its intended value now is 1: