[Boost-commit] svn:boost r86550 - in trunk: boost/math/special_functions libs/math/test

Subject: [Boost-commit] svn:boost r86550 - in trunk: boost/math/special_functions libs/math/test
From: john_at_[hidden]
Date: 2013-11-03 12:55:22

Author: johnmaddock
Date: 2013-11-03 12:55:22 EST (Sun, 03 Nov 2013)
New Revision: 86550
URL: http://svn.boost.org/trac/boost/changeset/86550

Log:
Fixes for digamma and zeta based on bug reports by Rocco Romeo.

Text files modified: 
   trunk/boost/math/special_functions/digamma.hpp |     8 +++++---                                
   trunk/boost/math/special_functions/zeta.hpp    |    26 +++++++++++++++++++++-----              
   trunk/libs/math/test/Jamfile.v2                |     2 +-                                      
   trunk/libs/math/test/test_digamma.hpp          |    13 +++++++++++++                           
   trunk/libs/math/test/test_zeta.cpp             |     2 +-                                      
   trunk/libs/math/test/test_zeta.hpp             |    25 ++++++++++++++++++++++++-               
   6 files changed, 65 insertions(+), 11 deletions(-)

Modified: trunk/boost/math/special_functions/digamma.hpp
==============================================================================
--- trunk/boost/math/special_functions/digamma.hpp	Sun Nov  3 12:33:00 2013	(r86549)
+++ trunk/boost/math/special_functions/digamma.hpp	2013-11-03 12:55:22 EST (Sun, 03 Nov 2013)	(r86550)
@@ -356,7 +356,7 @@
    //
    // Check for negative arguments and use reflection:
    //
-   if(x < 0)
+   if(x <= -1)
    {
       // Reflect:
       x = 1 - x;
@@ -376,6 +376,8 @@
       }
       result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
    }
+   if(x == 0)
+      return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, x, pol);
    //
    // If we're above the lower-limit for the
    // asymptotic expansion then use it:
@@ -397,9 +399,9 @@
       //
       // If x < 1 use recurrance to shift to > 1:
       //
-      if(x < 1)
+      while(x < 1)
       {
-         result = -1/x;
+         result -= 1/x;
          x += 1;
       }
       result += digamma_imp_1_2(x, t);

Modified: trunk/boost/math/special_functions/zeta.hpp
==============================================================================
--- trunk/boost/math/special_functions/zeta.hpp	Sun Nov  3 12:33:00 2013	(r86549)
+++ trunk/boost/math/special_functions/zeta.hpp	2013-11-03 12:55:22 EST (Sun, 03 Nov 2013)	(r86550)
@@ -866,9 +866,10 @@
 T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag)
 {
    BOOST_MATH_STD_USING
+   static const char* function = "boost::math::zeta<%1%>";
    if(sc == 0)
       return policies::raise_pole_error<T>(
-         "boost::math::zeta<%1%>", 
+         function, 
          "Evaluation of zeta function at pole %1%", 
          s, pol);
    T result;
@@ -883,10 +884,25 @@
          result = 0;
       else
       {
-         result = boost::math::sin_pi(0.5f * sc, pol)
-            * 2 * pow(2 * constants::pi<T>(), -s) 
-            * boost::math::tgamma(s, pol) 
-            * zeta_imp(s, sc, pol, tag);
+         if(s > max_factorial<T>::value)
+         {
+            T mult = boost::math::sin_pi(0.5f * sc, pol) * 2 * zeta_imp(s, sc, pol, tag);
+            result = boost::math::lgamma(s, pol);
+            result -= s * log(2 * constants::pi<T>());
+            if(result > tools::log_max_value<T>())
+               return sign(mult) * policies::raise_overflow_error<T>(function, 0, pol);
+            result = exp(result);
+            if(tools::max_value<T>() / fabs(mult) < result)
+               return boost::math::sign(mult) * policies::raise_overflow_error<T>(function, 0, pol);
+            result *= mult;
+         }
+         else
+         {
+            result = boost::math::sin_pi(0.5f * sc, pol)
+               * 2 * pow(2 * constants::pi<T>(), -s) 
+               * boost::math::tgamma(s, pol) 
+               * zeta_imp(s, sc, pol, tag);
+         }
       }
    }
    else

Modified: trunk/libs/math/test/Jamfile.v2
==============================================================================
--- trunk/libs/math/test/Jamfile.v2	Sun Nov  3 12:33:00 2013	(r86549)
+++ trunk/libs/math/test/Jamfile.v2	2013-11-03 12:55:22 EST (Sun, 03 Nov 2013)	(r86550)
@@ -648,7 +648,7 @@
 run test_triangular.cpp pch ../../test/build//boost_unit_test_framework ;
 run test_uniform.cpp pch ../../test/build//boost_unit_test_framework ;
 run test_weibull.cpp ../../test/build//boost_unit_test_framework ;
-run test_zeta.cpp test_instances//test_instances ../../test/build//boost_unit_test_framework pch_light ;
+run test_zeta.cpp ../../test/build//boost_unit_test_framework pch ;
 
 run test_policy.cpp ../../test/build//boost_unit_test_framework ;
 run test_policy_2.cpp ../../test/build//boost_unit_test_framework ;

Modified: trunk/libs/math/test/test_digamma.hpp
==============================================================================
--- trunk/libs/math/test/test_digamma.hpp	Sun Nov  3 12:33:00 2013	(r86549)
+++ trunk/libs/math/test/test_digamma.hpp	2013-11-03 12:55:22 EST (Sun, 03 Nov 2013)	(r86550)
@@ -72,5 +72,18 @@
 
    do_test_digamma<T>(digamma_neg_data, name, "Digamma Function: Negative Values");
 
+    static const boost::array<boost::array<T, 2>, 5> digamma_bugs = {{
+       // Test cases from Rocco Romeo:
+        {{ static_cast<T>(std::ldexp(1.0, -100)), SC_(-1.26765060022822940149670320537657721566490153286060651209008e30) }},
+        {{ static_cast<T>(-std::ldexp(1.0, -100)), SC_(1.26765060022822940149670320537542278433509846713939348790992e30) }},
+        {{ static_cast<T>(1), SC_(-0.577215664901532860606512090082402431042159335939923598805767) }},
+        {{ SC_(-1) + static_cast<T>(std::ldexp(1.0, -20)), SC_(-1.04857557721314249602848739817764518743062133735858753112190e6) }},
+        {{ SC_(-1) - static_cast<T>(std::ldexp(1.0, -20)), SC_(1.04857642278181269259522681939281063878220298942888100442172e6) }},
+    }};
+   do_test_digamma<T>(digamma_bugs, name, "Digamma Function: Values near 0");
+
+   BOOST_CHECK_THROW(boost::math::digamma(T(0)), std::domain_error);
+   BOOST_CHECK_THROW(boost::math::digamma(T(-1)), std::domain_error);
+   BOOST_CHECK_THROW(boost::math::digamma(T(-2)), std::domain_error);
 }
 

Modified: trunk/libs/math/test/test_zeta.cpp
==============================================================================
--- trunk/libs/math/test/test_zeta.cpp	Sun Nov  3 12:33:00 2013	(r86549)
+++ trunk/libs/math/test/test_zeta.cpp	2013-11-03 12:55:22 EST (Sun, 03 Nov 2013)	(r86550)
@@ -3,7 +3,7 @@
 //  Boost Software License, Version 1.0. (See accompanying file
 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 
-#include "pch_light.hpp"
+#include "pch.hpp"
 #include "test_zeta.hpp"
 
 //

Modified: trunk/libs/math/test/test_zeta.hpp
==============================================================================
--- trunk/libs/math/test/test_zeta.hpp	Sun Nov  3 12:33:00 2013	(r86549)
+++ trunk/libs/math/test/test_zeta.hpp	2013-11-03 12:55:22 EST (Sun, 03 Nov 2013)	(r86550)
@@ -5,7 +5,7 @@
 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 
 #include <boost/math/concepts/real_concept.hpp>
-#include <boost/math/special_functions/math_fwd.hpp>
+#include <boost/math/special_functions/zeta.hpp>
 #define BOOST_TEST_MAIN
 #include <boost/test/unit_test.hpp>
 #include <boost/test/floating_point_comparison.hpp>
@@ -150,5 +150,28 @@
    BOOST_CHECK_CLOSE(::boost::math::zeta(-ldexp(static_cast<T>(1), -32)), static_cast<T>(-0.499999999786042949889995597926798240562852438685508646794693L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::zeta(-ldexp(static_cast<T>(1), -33)), static_cast<T>(-0.499999999893021474931402198791408471637626205588681812641711L), tolerance);
    BOOST_CHECK_CLOSE(::boost::math::zeta(-ldexp(static_cast<T>(1), -34)), static_cast<T>(-0.499999999946510737462302199352114463422268928922372277519378L), tolerance);
+#ifdef BOOST_MSVC
+#pragma warning(push)
+#pragma warning(disable:4127 4756)
+#endif
+   //
+   // Very large negative values need special handling in our code, test them here, due to a bug report by Rocco Romeo:
+   //
+   BOOST_CHECK_EQUAL(::boost::math::zeta(static_cast<T>(-200)), static_cast<T>(0));
+   if(std::numeric_limits<T>::max_exponent >= 1024)
+   {
+      BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(-171)), static_cast<T>(1.28194898634822427378088228065956967928127061276520385040051e172L), tolerance * 10);
+      BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(-171.5)), static_cast<T>(4.73930233055054501360661283732419615206017226423071857829425e172L), tolerance * 1000);
+      BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(-172.5)), static_cast<T>(-1.30113885243175165293156588942160456456090687128236657847674e174L), tolerance * 100);
+      BOOST_CHECK_CLOSE(::boost::math::zeta(static_cast<T>(-173)), static_cast<T>(-9.66241211085609184243169684777934860657838245104636064505158e174L), tolerance * 100);
+   }
+   if(std::numeric_limits<T>::has_infinity)
+   {
+      BOOST_CHECK_EQUAL(::boost::math::zeta(-10007, boost::math::policies::make_policy(boost::math::policies::overflow_error<boost::math::policies::ignore_error>())), std::numeric_limits<T>::infinity());
+      BOOST_CHECK_EQUAL(::boost::math::zeta(-10009, boost::math::policies::make_policy(boost::math::policies::overflow_error<boost::math::policies::ignore_error>())), -std::numeric_limits<T>::infinity());
+   }
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
 }