[Polygamma]

Fix real_concept compilation and runtime.
Add digits_base10 support function to policies.

include/boost/math/policies/policy.hpp

@@ -850,6 +850,11 @@ inline int digits(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))
    typedef mpl::bool_< std::numeric_limits<T>::is_specialized > tag_type;
    return detail::digits_imp<T, Policy>(tag_type());
 }
+template <class T, class Policy>
+inline int digits_base10(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))
+{
+   return boost::math::policies::digits<T, Policy>() * 1000L / 301;
+}
 
 template <class Policy>
 inline unsigned long get_max_series_iterations()

include/boost/math/special_functions/detail/polygamma.hpp

@@ -23,7 +23,9 @@
   #include <boost/static_assert.hpp>
   #include <boost/type_traits/is_convertible.hpp>
 
-  namespace boost { namespace math {
+  namespace boost { namespace math { namespace detail{
+
+#if 0
 
   template<class T, class Policy>
   T digamma_atinfinityplus(const int, const T &x, const Policy&)
@@ -66,7 +68,7 @@
       sum_base_10_exp = T(exponent_value) * 0.303F;
 
       long int order_check =  boost::math::ltrunc(term_base_10_exp) - boost::math::ltrunc(sum_base_10_exp);
-      long int tol         =  std::numeric_limits<T>::digits10;
+      long int tol         =  std::numeric_limits<T>::digits_base10;
 
 
       if((two_k > 24) && (order_check < -tol))
@@ -77,6 +79,7 @@
 
     return (log_z - one_over_2z) - sum;
   }
+#endif
 
   template<class T, class Policy>
   T polygamma_atinfinityplus(const int n, const T& x, const Policy& pol) // for large values of x such as for x> 400
@@ -154,9 +157,7 @@
         //
         if(k > policies::get_max_series_iterations<Policy>())
         {
-           policies::raise_evaluation_error(
-              "polygamma<%1%>(int, %1%)",
-              "Series did not converge, closest value was %1%", sum, pol);
+           policies::raise_evaluation_error("polygamma<%1%>(int, %1%)", "Series did not converge, closest value was %1%", sum, pol);
            break;
         }
      }
@@ -174,7 +175,7 @@
 
     // Use N = (0.4 * digits) + (4 * n) for target value for x:
     BOOST_MATH_STD_USING
-    const int d4d  = static_cast<boost::int32_t>(0.4F * std::numeric_limits<T>::digits10);
+    const int d4d  = static_cast<boost::int32_t>(0.4F * policies::digits_base10<T, Policy>());
     const int N4dn = static_cast<boost::int32_t>(d4d + (4 * n));
     const int N    = static_cast<boost::int32_t>((std::min)(N4dn, (std::numeric_limits<int>::max)()));
     const int m    = n;
@@ -201,7 +202,7 @@
   }
 
   template<class T, class Policy>
-  T polygamma_nearzero(const int n, const T& x, const Policy&)
+  T polygamma_nearzero(const int n, const T& x, const Policy& pol)
   {
     BOOST_MATH_STD_USING
     // not defined for digamma
@@ -224,9 +225,10 @@
           bool    b_neg_term     =  ((n % 2) == 0);
           T sum                  =  ((!b_neg_term) ? pg_kn : -pg_kn);
 
-    for(int k = 1; k < max_iteration<T>::value; k++)
+    for(unsigned k = 1;; k++)
     {
-      k_plus_n_fact   *= k_plus_n_plus_one++;
+      k_plus_n_fact   *= k_plus_n_plus_one;
+      k_plus_n_plus_one += 1;
       one_over_k_fact /= k;
       z_pow_k         *= x;
 
@@ -234,21 +236,7 @@
 
       const T term = (pg * z_pow_k) * one_over_k_fact;
 
-      T term_base_10_exp = ((term < 0) ? -term : term);
-      T sum_base_10_exp  = ((sum  < 0) ? -sum  : sum);
-
-      int exponent_value;
-
-      static_cast<void>(frexp(term_base_10_exp, &exponent_value));
-      term_base_10_exp = T(exponent_value) * 0.303F;
-
-      static_cast<void>(frexp(sum_base_10_exp, &exponent_value));
-      sum_base_10_exp = T(exponent_value) * 0.303F;
-
-      long int order_check =  boost::math::ltrunc(term_base_10_exp) - boost::math::ltrunc(sum_base_10_exp);
-      long int tol         =  std::numeric_limits<T>::digits10;
-
-      if((k > 12) && (order_check < -tol))
+      if(fabs(term / sum) < tools::epsilon<T>())
       {
         break;
       }
@@ -256,6 +244,12 @@
       b_neg_term = !b_neg_term;
 
       ((!b_neg_term) ? sum += term : sum -= term);
+
+      if(k > policies::get_max_series_iterations<Policy>())
+      {
+         policies::raise_evaluation_error("polygamma<%1%>(int, %1%)", "Series did not converge, closest value was %1%", sum, pol);
+         break;
+      }
     }
 
     return term0 + sum;
@@ -269,7 +263,7 @@
     {
       return polygamma_nearzero(n, x, pol);
     }
-    else if(x > 0.4F * std::numeric_limits<T>::digits10 + 4 * n)
+    else if(x > 0.4F * policies::digits_base10<T, Policy>() + 4 * n)
     {
       return polygamma_atinfinityplus(n, x, pol);
     }

include/boost/math/special_functions/polygamma.hpp

@@ -18,33 +18,57 @@
 
   namespace boost { namespace math {
 
-  template<class T>	  
-  struct promoteftod
-  {
-	  typedef T type;
-  };  
-
-  template<>
-  struct promoteftod<float>
-  {
-	  typedef double type;
-  };
-
+  
   template<class T, class Policy>
-  inline T polygamma(const int n, T x, const Policy &pol)
+  inline typename tools::promote_args<T>::type polygamma(const int n, T x, const Policy &pol)
   {
-	typedef typename promoteftod<T>::type result_type;
-//	std::cout<<"~:"<<typeid(T).name()<<std::endl;
-//	std::cout<<"~:"<<typeid(result_type).name()<<std::endl;
-	result_type xx=result_type(x);
-        result_type result= boost::math::detail::polygamma_imp(n,xx,pol);
-	return T(result);
+     //
+     // We've found some standard library functions to misbehave if any FPU exception flags
+     // are set prior to their call, this code will clear those flags, then reset them
+     // on exit:
+     //
+     BOOST_FPU_EXCEPTION_GUARD
+     //
+     // The type of the result - the common type of T and U after
+     // any integer types have been promoted to double:
+     //
+     typedef typename tools::promote_args<T>::type result_type;
+     //
+     // The type used for the calculation.  This may be a wider type than
+     // the result in order to ensure full precision:
+     //
+     typedef typename policies::evaluation<result_type, Policy>::type value_type;
+     //
+     // The type of the policy to forward to the actual implementation.
+     // We disable promotion of float and double as that's [possibly]
+     // happened already in the line above.  Also reset to the default
+     // any policies we don't use (reduces code bloat if we're called
+     // multiple times with differing policies we don't actually use).
+     // Also normalise the type, again to reduce code bloat in case we're
+     // called multiple times with functionally identical policies that happen
+     // to be different types.
+     //
+     typedef typename policies::normalise<
+        Policy,
+        policies::promote_float<false>,
+        policies::promote_double<false>,
+        policies::discrete_quantile<>,
+        policies::assert_undefined<> >::type forwarding_policy;
+     //
+     // Whew.  Now we can make the actual call to the implementation.
+     // Arguments are explicitly cast to the evaluation type, and the result
+     // passed through checked_narrowing_cast which handles things like overflow
+     // according to the policy passed:
+     //
+     return policies::checked_narrowing_cast<result_type, forwarding_policy>(
+        detail::polygamma_imp(n, static_cast<value_type>(x), forwarding_policy()),
+        "boost::math::polygamma<%1%>(int, %1%)");
   }
 
   template<class T>
-  inline T polygamma(const int n, T x)
+  inline typename tools::promote_args<T>::type polygamma(const int n, T x)
   {
-      return boost::math::polygamma(n,x,policies::policy<>());
+      return boost::math::polygamma(n, x, policies::policy<>());
   }
 
   template<class T, class Policy>

test/Jamfile.v2

@@ -603,6 +603,7 @@ run test_nc_t.cpp  pch ../../test/build//boost_unit_test_framework
 run test_normal.cpp pch ../../test/build//boost_unit_test_framework ;
 run test_owens_t.cpp ../../test/build//boost_unit_test_framework ;
 run test_pareto.cpp ../../test/build//boost_unit_test_framework ;
+run test_polygamma.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ;
 run test_poisson.cpp ../../test/build//boost_unit_test_framework
         : # command line
         : # input files

test/test_polygamma.cpp

@@ -5,7 +5,7 @@
 
 #include <pch_light.hpp>
 #include "test_polygamma.hpp"
-
+//#include <boost/cstdfloat.hpp>
 
 void expected_results()
 {
@@ -17,16 +17,30 @@ void expected_results()
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
    {
-      largest_type = "(long\\s+)?double";
+      largest_type = "(long\\s+)?double|real_concept";
    }
    else
    {
-      largest_type = "long double";
+      largest_type = "long double|real_concept";
    }
 #else
-   largest_type = "(long\\s+)?double";
+   largest_type = "(long\\s+)?double|real_concept";
 #endif
 
+   add_expected_result(
+      ".*",                          // compiler
+      ".*",                          // stdlib
+      ".*",                          // platform
+      largest_type,                  // test type(s)
+      ".*large arguments",           // test data group
+      ".*", 400, 200);               // test function
+   add_expected_result(
+      ".*",                          // compiler
+      ".*",                          // stdlib
+      ".*",                          // platform
+      largest_type,                  // test type(s)
+      ".*",                          // test data group
+      ".*", 20, 10);                 // test function
    //
    // Finish off by printing out the compiler/stdlib/platform names,
    // we do this to make it easier to mark up expected error rates.
@@ -40,11 +54,14 @@ BOOST_AUTO_TEST_CASE( test_main )
    expected_results();
    BOOST_MATH_CONTROL_FP;
 
-   //test_polygamma(0.0F, "float");
+   test_polygamma(0.0F, "float");
    test_polygamma(0.0, "double");
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
    test_polygamma(0.0L, "long double");
-   //test_polygamma(boost::math::concepts::real_concept(0.1), "real_concept");
+   test_polygamma(boost::math::concepts::real_concept(0.1), "real_concept");
+#endif
+#ifdef BOOST_FLOAT128_C
+   //test_polygamma(BOOST_FLOAT128_C(0.0), "float128_t");
 #endif
 }