Chris@16: //  (C) Copyright John Maddock 2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SF_DIGAMMA_HPP
Chris@16: #define BOOST_MATH_SF_DIGAMMA_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@101: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/math/tools/rational.hpp>
Chris@101: #include <boost/math/tools/series.hpp>
Chris@16: #include <boost/math/tools/promotion.hpp>
Chris@16: #include <boost/math/policies/error_handling.hpp>
Chris@16: #include <boost/math/constants/constants.hpp>
Chris@16: #include <boost/mpl/comparison.hpp>
Chris@16: #include <boost/math/tools/big_constant.hpp>
Chris@16: 
Chris@16: namespace boost{
Chris@16: namespace math{
Chris@16: namespace detail{
Chris@16: //
Chris@16: // Begin by defining the smallest value for which it is safe to
Chris@16: // use the asymptotic expansion for digamma:
Chris@16: //
Chris@16: inline unsigned digamma_large_lim(const mpl::int_<0>*)
Chris@16: {  return 20;  }
Chris@101: inline unsigned digamma_large_lim(const mpl::int_<113>*)
Chris@101: {  return 20;  }
Chris@16: inline unsigned digamma_large_lim(const void*)
Chris@16: {  return 10;  }
Chris@16: //
Chris@16: // Implementations of the asymptotic expansion come next,
Chris@16: // the coefficients of the series have been evaluated
Chris@16: // in advance at high precision, and the series truncated
Chris@16: // at the first term that's too small to effect the result.
Chris@16: // Note that the series becomes divergent after a while
Chris@16: // so truncation is very important.
Chris@16: //
Chris@16: // This first one gives 34-digit precision for x >= 20:
Chris@16: //
Chris@16: template <class T>
Chris@101: inline T digamma_imp_large(T x, const mpl::int_<113>*)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std functions.
Chris@16:    static const T P[] = {
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.003968253968253968253968253968253968253968253968254),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0041666666666666666666666666666666666666666666666667),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0075757575757575757575757575757575757575757575757576),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.021092796092796092796092796092796092796092796092796),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.44325980392156862745098039215686274509803921568627),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 3.0539543302701197438039543302701197438039543302701),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -26.456212121212121212121212121212121212121212121212),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 281.4601449275362318840579710144927536231884057971),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -3607.510546398046398046398046398046398046398046398),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 54827.583333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -974936.82385057471264367816091954022988505747126437),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 20052695.796688078946143462272494530559046688078946),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -472384867.72162990196078431372549019607843137254902),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 12635724795.916666666666666666666666666666666666667)
Chris@16:    };
Chris@16:    x -= 1;
Chris@16:    T result = log(x);
Chris@16:    result += 1 / (2 * x);
Chris@16:    T z = 1 / (x*x);
Chris@16:    result -= z * tools::evaluate_polynomial(P, z);
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // 19-digit precision for x >= 10:
Chris@16: //
Chris@16: template <class T>
Chris@16: inline T digamma_imp_large(T x, const mpl::int_<64>*)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std functions.
Chris@16:    static const T P[] = {
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.003968253968253968253968253968253968253968253968254),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0041666666666666666666666666666666666666666666666667),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0075757575757575757575757575757575757575757575757576),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.021092796092796092796092796092796092796092796092796),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.44325980392156862745098039215686274509803921568627),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 3.0539543302701197438039543302701197438039543302701),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -26.456212121212121212121212121212121212121212121212),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 281.4601449275362318840579710144927536231884057971),
Chris@16:    };
Chris@16:    x -= 1;
Chris@16:    T result = log(x);
Chris@16:    result += 1 / (2 * x);
Chris@16:    T z = 1 / (x*x);
Chris@16:    result -= z * tools::evaluate_polynomial(P, z);
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // 17-digit precision for x >= 10:
Chris@16: //
Chris@16: template <class T>
Chris@16: inline T digamma_imp_large(T x, const mpl::int_<53>*)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std functions.
Chris@16:    static const T P[] = {
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 0.003968253968253968253968253968253968253968253968254),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0041666666666666666666666666666666666666666666666667),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 0.0075757575757575757575757575757575757575757575757576),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.021092796092796092796092796092796092796092796092796),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.44325980392156862745098039215686274509803921568627)
Chris@16:    };
Chris@16:    x -= 1;
Chris@16:    T result = log(x);
Chris@16:    result += 1 / (2 * x);
Chris@16:    T z = 1 / (x*x);
Chris@16:    result -= z * tools::evaluate_polynomial(P, z);
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // 9-digit precision for x >= 10:
Chris@16: //
Chris@16: template <class T>
Chris@16: inline T digamma_imp_large(T x, const mpl::int_<24>*)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std functions.
Chris@16:    static const T P[] = {
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 24, 0.083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 24, -0.0083333333333333333333333333333333333333333333333333),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 24, 0.003968253968253968253968253968253968253968253968254)
Chris@16:    };
Chris@16:    x -= 1;
Chris@16:    T result = log(x);
Chris@16:    result += 1 / (2 * x);
Chris@16:    T z = 1 / (x*x);
Chris@16:    result -= z * tools::evaluate_polynomial(P, z);
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@101: // Fully generic asymptotic expansion in terms of Bernoulli numbers, see:
Chris@101: // http://functions.wolfram.com/06.14.06.0012.01
Chris@101: //
Chris@101: template <class T>
Chris@101: struct digamma_series_func
Chris@101: {
Chris@101: private:
Chris@101:    int k;
Chris@101:    T xx;
Chris@101:    T term;
Chris@101: public:
Chris@101:    digamma_series_func(T x) : k(1), xx(x * x), term(1 / (x * x)) {}
Chris@101:    T operator()()
Chris@101:    {
Chris@101:       T result = term * boost::math::bernoulli_b2n<T>(k) / (2 * k);
Chris@101:       term /= xx;
Chris@101:       ++k;
Chris@101:       return result;
Chris@101:    }
Chris@101:    typedef T result_type;
Chris@101: };
Chris@101: 
Chris@101: template <class T, class Policy>
Chris@101: inline T digamma_imp_large(T x, const Policy& pol, const mpl::int_<0>*)
Chris@101: {
Chris@101:    BOOST_MATH_STD_USING
Chris@101:    digamma_series_func<T> s(x);
Chris@101:    T result = log(x) - 1 / (2 * x);
Chris@101:    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@101:    result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, -result);
Chris@101:    result = -result;
Chris@101:    policies::check_series_iterations<T>("boost::math::digamma<%1%>(%1%)", max_iter, pol);
Chris@101:    return result;
Chris@101: }
Chris@101: //
Chris@16: // Now follow rational approximations over the range [1,2].
Chris@16: //
Chris@16: // 35-digit precision:
Chris@16: //
Chris@16: template <class T>
Chris@101: T digamma_imp_1_2(T x, const mpl::int_<113>*)
Chris@16: {
Chris@16:    //
Chris@16:    // Now the approximation, we use the form:
Chris@16:    //
Chris@16:    // digamma(x) = (x - root) * (Y + R(x-1))
Chris@16:    //
Chris@16:    // Where root is the location of the positive root of digamma,
Chris@16:    // Y is a constant, and R is optimised for low absolute error
Chris@16:    // compared to Y.
Chris@16:    //
Chris@16:    // Max error found at 128-bit long double precision:  5.541e-35
Chris@16:    // Maximum Deviation Found (approximation error):     1.965e-35
Chris@16:    //
Chris@16:    static const float Y = 0.99558162689208984375F;
Chris@16: 
Chris@16:    static const T root1 = T(1569415565) / 1073741824uL;
Chris@16:    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
Chris@16:    static const T root3 = ((T(111616537) / 1073741824uL) / 1073741824uL) / 1073741824uL;
Chris@16:    static const T root4 = (((T(503992070) / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL;
Chris@16:    static const T root5 = BOOST_MATH_BIG_CONSTANT(T, 113, 0.52112228569249997894452490385577338504019838794544e-36);
Chris@16: 
Chris@16:    static const T P[] = {    
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.25479851061131551526977464225335883769),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.18684290534374944114622235683619897417),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.80360876047931768958995775910991929922),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.67227342794829064330498117008564270136),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.26569010991230617151285010695543858005),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.05775672694575986971640757748003553385),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0071432147823164975485922555833274240665),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00048740753910766168912364555706064993274),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.16454996865214115723416538844975174761e-4),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.20327832297631728077731148515093164955e-6)
Chris@16:    };
Chris@16:    static const T Q[] = {    
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 2.6210924610812025425088411043163287646),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 2.6850757078559596612621337395886392594),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 1.4320913706209965531250495490639289418),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.4410872083455009362557012239501953402),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.081385727399251729505165509278152487225),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0089478633066857163432104815183858149496),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00055861622855066424871506755481997374154),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.1760168552357342401304462967950178554e-4),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.20585454493572473724556649516040874384e-6),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, -0.90745971844439990284514121823069162795e-11),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 113, 0.48857673606545846774761343500033283272e-13),
Chris@16:    };
Chris@16:    T g = x - root1;
Chris@16:    g -= root2;
Chris@16:    g -= root3;
Chris@16:    g -= root4;
Chris@16:    g -= root5;
Chris@16:    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
Chris@16:    T result = g * Y + g * r;
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // 19-digit precision:
Chris@16: //
Chris@16: template <class T>
Chris@16: T digamma_imp_1_2(T x, const mpl::int_<64>*)
Chris@16: {
Chris@16:    //
Chris@16:    // Now the approximation, we use the form:
Chris@16:    //
Chris@16:    // digamma(x) = (x - root) * (Y + R(x-1))
Chris@16:    //
Chris@16:    // Where root is the location of the positive root of digamma,
Chris@16:    // Y is a constant, and R is optimised for low absolute error
Chris@16:    // compared to Y.
Chris@16:    //
Chris@16:    // Max error found at 80-bit long double precision:   5.016e-20
Chris@16:    // Maximum Deviation Found (approximation error):     3.575e-20
Chris@16:    //
Chris@16:    static const float Y = 0.99558162689208984375F;
Chris@16: 
Chris@16:    static const T root1 = T(1569415565) / 1073741824uL;
Chris@16:    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
Chris@16:    static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 64, 0.9016312093258695918615325266959189453125e-19);
Chris@16: 
Chris@16:    static const T P[] = {    
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.254798510611315515235),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.314628554532916496608),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.665836341559876230295),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.314767657147375752913),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0541156266153505273939),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00289268368333918761452)
Chris@16:    };
Chris@16:    static const T Q[] = {    
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 2.1195759927055347547),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 1.54350554664961128724),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.486986018231042975162),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0660481487173569812846),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00298999662592323990972),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, -0.165079794012604905639e-5),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 64, 0.317940243105952177571e-7)
Chris@16:    };
Chris@16:    T g = x - root1;
Chris@16:    g -= root2;
Chris@16:    g -= root3;
Chris@16:    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
Chris@16:    T result = g * Y + g * r;
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // 18-digit precision:
Chris@16: //
Chris@16: template <class T>
Chris@16: T digamma_imp_1_2(T x, const mpl::int_<53>*)
Chris@16: {
Chris@16:    //
Chris@16:    // Now the approximation, we use the form:
Chris@16:    //
Chris@16:    // digamma(x) = (x - root) * (Y + R(x-1))
Chris@16:    //
Chris@16:    // Where root is the location of the positive root of digamma,
Chris@16:    // Y is a constant, and R is optimised for low absolute error
Chris@16:    // compared to Y.
Chris@16:    //
Chris@16:    // Maximum Deviation Found:               1.466e-18
Chris@16:    // At double precision, max error found:  2.452e-17
Chris@16:    //
Chris@16:    static const float Y = 0.99558162689208984F;
Chris@16: 
Chris@16:    static const T root1 = T(1569415565) / 1073741824uL;
Chris@16:    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
Chris@16:    static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 53, 0.9016312093258695918615325266959189453125e-19);
Chris@16: 
Chris@16:    static const T P[] = {    
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 0.25479851061131551),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.32555031186804491),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.65031853770896507),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.28919126444774784),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.045251321448739056),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0020713321167745952)
Chris@16:    };
Chris@16:    static const T Q[] = {    
Chris@101:       BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 2.0767117023730469),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 1.4606242909763515),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 0.43593529692665969),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 0.054151797245674225),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, 0.0021284987017821144),
Chris@16:       BOOST_MATH_BIG_CONSTANT(T, 53, -0.55789841321675513e-6)
Chris@16:    };
Chris@16:    T g = x - root1;
Chris@16:    g -= root2;
Chris@16:    g -= root3;
Chris@16:    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
Chris@16:    T result = g * Y + g * r;
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // 9-digit precision:
Chris@16: //
Chris@16: template <class T>
Chris@16: inline T digamma_imp_1_2(T x, const mpl::int_<24>*)
Chris@16: {
Chris@16:    //
Chris@16:    // Now the approximation, we use the form:
Chris@16:    //
Chris@16:    // digamma(x) = (x - root) * (Y + R(x-1))
Chris@16:    //
Chris@16:    // Where root is the location of the positive root of digamma,
Chris@16:    // Y is a constant, and R is optimised for low absolute error
Chris@16:    // compared to Y.
Chris@16:    //
Chris@16:    // Maximum Deviation Found:              3.388e-010
Chris@16:    // At float precision, max error found:  2.008725e-008
Chris@16:    //
Chris@16:    static const float Y = 0.99558162689208984f;
Chris@16:    static const T root = 1532632.0f / 1048576;
Chris@16:    static const T root_minor = static_cast<T>(0.3700660185912626595423257213284682051735604e-6L);
Chris@16:    static const T P[] = {    
Chris@101:       0.25479851023250261e0f,
Chris@101:       -0.44981331915268368e0f,
Chris@101:       -0.43916936919946835e0f,
Chris@101:       -0.61041765350579073e-1f
Chris@16:    };
Chris@16:    static const T Q[] = {    
Chris@16:       0.1e1,
Chris@101:       0.15890202430554952e1f,
Chris@101:       0.65341249856146947e0f,
Chris@101:       0.63851690523355715e-1f
Chris@16:    };
Chris@16:    T g = x - root;
Chris@16:    g -= root_minor;
Chris@16:    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
Chris@16:    T result = g * Y + g * r;
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Tag, class Policy>
Chris@16: T digamma_imp(T x, const Tag* t, const Policy& pol)
Chris@16: {
Chris@16:    //
Chris@16:    // This handles reflection of negative arguments, and all our
Chris@16:    // error handling, then forwards to the T-specific approximation.
Chris@16:    //
Chris@16:    BOOST_MATH_STD_USING // ADL of std functions.
Chris@16: 
Chris@16:    T result = 0;
Chris@16:    //
Chris@16:    // Check for negative arguments and use reflection:
Chris@16:    //
Chris@101:    if(x <= -1)
Chris@16:    {
Chris@16:       // Reflect:
Chris@16:       x = 1 - x;
Chris@16:       // Argument reduction for tan:
Chris@16:       T remainder = x - floor(x);
Chris@16:       // Shift to negative if > 0.5:
Chris@16:       if(remainder > 0.5)
Chris@16:       {
Chris@16:          remainder -= 1;
Chris@16:       }
Chris@16:       //
Chris@16:       // check for evaluation at a negative pole:
Chris@16:       //
Chris@16:       if(remainder == 0)
Chris@16:       {
Chris@16:          return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol);
Chris@16:       }
Chris@16:       result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
Chris@16:    }
Chris@101:    if(x == 0)
Chris@101:       return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, x, pol);
Chris@16:    //
Chris@16:    // If we're above the lower-limit for the
Chris@16:    // asymptotic expansion then use it:
Chris@16:    //
Chris@16:    if(x >= digamma_large_lim(t))
Chris@16:    {
Chris@16:       result += digamma_imp_large(x, t);
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       //
Chris@16:       // If x > 2 reduce to the interval [1,2]:
Chris@16:       //
Chris@16:       while(x > 2)
Chris@16:       {
Chris@16:          x -= 1;
Chris@16:          result += 1/x;
Chris@16:       }
Chris@16:       //
Chris@16:       // If x < 1 use recurrance to shift to > 1:
Chris@16:       //
Chris@101:       while(x < 1)
Chris@16:       {
Chris@101:          result -= 1/x;
Chris@16:          x += 1;
Chris@16:       }
Chris@16:       result += digamma_imp_1_2(x, t);
Chris@16:    }
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@101: template <class T, class Policy>
Chris@101: T digamma_imp(T x, const mpl::int_<0>* t, const Policy& pol)
Chris@101: {
Chris@101:    //
Chris@101:    // This handles reflection of negative arguments, and all our
Chris@101:    // error handling, then forwards to the T-specific approximation.
Chris@101:    //
Chris@101:    BOOST_MATH_STD_USING // ADL of std functions.
Chris@101: 
Chris@101:    T result = 0;
Chris@101:    //
Chris@101:    // Check for negative arguments and use reflection:
Chris@101:    //
Chris@101:    if(x <= -1)
Chris@101:    {
Chris@101:       // Reflect:
Chris@101:       x = 1 - x;
Chris@101:       // Argument reduction for tan:
Chris@101:       T remainder = x - floor(x);
Chris@101:       // Shift to negative if > 0.5:
Chris@101:       if(remainder > 0.5)
Chris@101:       {
Chris@101:          remainder -= 1;
Chris@101:       }
Chris@101:       //
Chris@101:       // check for evaluation at a negative pole:
Chris@101:       //
Chris@101:       if(remainder == 0)
Chris@101:       {
Chris@101:          return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1 - x), pol);
Chris@101:       }
Chris@101:       result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
Chris@101:    }
Chris@101:    if(x == 0)
Chris@101:       return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, x, pol);
Chris@101:    //
Chris@101:    // If we're above the lower-limit for the
Chris@101:    // asymptotic expansion then use it, the
Chris@101:    // limit is a linear interpolation with
Chris@101:    // limit = 10 at 50 bit precision and
Chris@101:    // limit = 250 at 1000 bit precision.
Chris@101:    //
Chris@101:    T lim = 10 + (tools::digits<T>() - 50) * 240 / 950;
Chris@101:    T two_x = ldexp(x, 1);
Chris@101:    if(x >= lim)
Chris@101:    {
Chris@101:       result += digamma_imp_large(x, pol, t);
Chris@101:    }
Chris@101:    else if(floor(x) == x)
Chris@101:    {
Chris@101:       //
Chris@101:       // Special case for integer arguments, see
Chris@101:       // http://functions.wolfram.com/06.14.03.0001.01
Chris@101:       //
Chris@101:       result = -constants::euler<T, Policy>();
Chris@101:       T val = 1;
Chris@101:       while(val < x)
Chris@101:       {
Chris@101:          result += 1 / val;
Chris@101:          val += 1;
Chris@101:       }
Chris@101:    }
Chris@101:    else if(floor(two_x) == two_x)
Chris@101:    {
Chris@101:       //
Chris@101:       // Special case for half integer arguments, see:
Chris@101:       // http://functions.wolfram.com/06.14.03.0007.01
Chris@101:       //
Chris@101:       result = -2 * constants::ln_two<T, Policy>() - constants::euler<T, Policy>();
Chris@101:       int n = itrunc(x);
Chris@101:       if(n)
Chris@101:       {
Chris@101:          for(int k = 1; k < n; ++k)
Chris@101:             result += 1 / T(k);
Chris@101:          for(int k = n; k <= 2 * n - 1; ++k)
Chris@101:             result += 2 / T(k);
Chris@101:       }
Chris@101:    }
Chris@101:    else
Chris@101:    {
Chris@101:       //
Chris@101:       // Rescale so we can use the asymptotic expansion:
Chris@101:       //
Chris@101:       while(x < lim)
Chris@101:       {
Chris@101:          result -= 1 / x;
Chris@101:          x += 1;
Chris@101:       }
Chris@101:       result += digamma_imp_large(x, pol, t);
Chris@101:    }
Chris@101:    return result;
Chris@101: }
Chris@16: //
Chris@16: // Initializer: ensure all our constants are initialized prior to the first call of main:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: struct digamma_initializer
Chris@16: {
Chris@16:    struct init
Chris@16:    {
Chris@16:       init()
Chris@16:       {
Chris@101:          typedef typename policies::precision<T, Policy>::type precision_type;
Chris@101:          do_init(mpl::bool_<precision_type::value && (precision_type::value <= 113)>());
Chris@101:       }
Chris@101:       void do_init(const mpl::true_&)
Chris@101:       {
Chris@16:          boost::math::digamma(T(1.5), Policy());
Chris@16:          boost::math::digamma(T(500), Policy());
Chris@16:       }
Chris@101:       void do_init(const mpl::false_&){}
Chris@16:       void force_instantiate()const{}
Chris@16:    };
Chris@16:    static const init initializer;
Chris@16:    static void force_instantiate()
Chris@16:    {
Chris@16:       initializer.force_instantiate();
Chris@16:    }
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: const typename digamma_initializer<T, Policy>::init digamma_initializer<T, Policy>::initializer;
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@101:    digamma(T x, const Policy&)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::precision<T, Policy>::type precision_type;
Chris@16:    typedef typename mpl::if_<
Chris@16:       mpl::or_<
Chris@16:          mpl::less_equal<precision_type, mpl::int_<0> >,
Chris@101:          mpl::greater<precision_type, mpl::int_<114> >
Chris@16:       >,
Chris@16:       mpl::int_<0>,
Chris@16:       typename mpl::if_<
Chris@16:          mpl::less<precision_type, mpl::int_<25> >,
Chris@16:          mpl::int_<24>,
Chris@16:          typename mpl::if_<
Chris@16:             mpl::less<precision_type, mpl::int_<54> >,
Chris@16:             mpl::int_<53>,
Chris@101:             typename mpl::if_<
Chris@101:                mpl::less<precision_type, mpl::int_<65> >,
Chris@101:                mpl::int_<64>,
Chris@101:                mpl::int_<113>
Chris@101:             >::type
Chris@16:          >::type
Chris@16:       >::type
Chris@16:    >::type tag_type;
Chris@16: 
Chris@101:    typedef typename policies::normalise<
Chris@101:       Policy,
Chris@101:       policies::promote_float<false>,
Chris@101:       policies::promote_double<false>,
Chris@101:       policies::discrete_quantile<>,
Chris@101:       policies::assert_undefined<> >::type forwarding_policy;
Chris@101: 
Chris@16:    // Force initialization of constants:
Chris@101:    detail::digamma_initializer<value_type, forwarding_policy>::force_instantiate();
Chris@16: 
Chris@16:    return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp(
Chris@16:       static_cast<value_type>(x),
Chris@101:       static_cast<const tag_type*>(0), forwarding_policy()), "boost::math::digamma<%1%>(%1%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@16:    digamma(T x)
Chris@16: {
Chris@16:    return digamma(x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: #endif
Chris@16: