1 2 /////////////////////////////////////////////////////////////////////////////// 3 // Copyright 2013 Nikhar Agrawal 4 // Copyright 2013 Christopher Kormanyos 5 // Copyright 2014 John Maddock 6 // Copyright 2013 Paul Bristow 7 // Distributed under the Boost 8 // Software License, Version 1.0. (See accompanying file 9 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 10 11 #ifndef _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ 12 #define _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ 13 14 #include <cmath> 15 #include <limits> 16 #include <boost/cstdint.hpp> 17 #include <boost/math/policies/policy.hpp> 18 #include <boost/math/special_functions/bernoulli.hpp> 19 #include <boost/math/special_functions/trunc.hpp> 20 #include <boost/math/special_functions/zeta.hpp> 21 #include <boost/math/special_functions/digamma.hpp> 22 #include <boost/math/special_functions/sin_pi.hpp> 23 #include <boost/math/special_functions/cos_pi.hpp> 24 #include <boost/math/special_functions/pow.hpp> 25 #include <boost/mpl/if.hpp> 26 #include <boost/mpl/int.hpp> 27 #include <boost/static_assert.hpp> 28 #include <boost/type_traits/is_convertible.hpp> 29 30 #ifdef _MSC_VER 31 #pragma once 32 #pragma warning(push) 33 #pragma warning(disable:4702) // Unreachable code (release mode only warning) 34 #endif 35 36 namespace boost { namespace math { namespace detail{ 37 38 template<class T, class Policy> polygamma_atinfinityplus(const int n,const T & x,const Policy & pol,const char * function)39 T polygamma_atinfinityplus(const int n, const T& x, const Policy& pol, const char* function) // for large values of x such as for x> 400 40 { 41 // See http://functions.wolfram.com/GammaBetaErf/PolyGamma2/06/02/0001/ 42 BOOST_MATH_STD_USING 43 // 44 // sum == current value of accumulated sum. 45 // term == value of current term to be added to sum. 46 // part_term == value of current term excluding the Bernoulli number part 47 // 48 if(n + x == x) 49 { 50 // x is crazy large, just concentrate on the first part of the expression and use logs: 51 if(n == 1) return 1 / x; 52 T nlx = n * log(x); 53 if((nlx < tools::log_max_value<T>()) && (n < (int)max_factorial<T>::value)) 54 return ((n & 1) ? 1 : -1) * boost::math::factorial<T>(n - 1) * pow(x, -n); 55 else 56 return ((n & 1) ? 1 : -1) * exp(boost::math::lgamma(T(n), pol) - n * log(x)); 57 } 58 T term, sum, part_term; 59 T x_squared = x * x; 60 // 61 // Start by setting part_term to: 62 // 63 // (n-1)! / x^(n+1) 64 // 65 // which is common to both the first term of the series (with k = 1) 66 // and to the leading part. 67 // We can then get to the leading term by: 68 // 69 // part_term * (n + 2 * x) / 2 70 // 71 // and to the first term in the series 72 // (excluding the Bernoulli number) by: 73 // 74 // part_term n * (n + 1) / (2x) 75 // 76 // If either the factorial would overflow, 77 // or the power term underflows, this just gets set to 0 and then we 78 // know that we have to use logs for the initial terms: 79 // 80 part_term = ((n > (int)boost::math::max_factorial<T>::value) && (T(n) * n > tools::log_max_value<T>())) 81 ? T(0) : static_cast<T>(boost::math::factorial<T>(n - 1, pol) * pow(x, -n - 1)); 82 if(part_term == 0) 83 { 84 // Either n is very large, or the power term underflows, 85 // set the initial values of part_term, term and sum via logs: 86 part_term = static_cast<T>(boost::math::lgamma(n, pol) - (n + 1) * log(x)); 87 sum = exp(part_term + log(n + 2 * x) - boost::math::constants::ln_two<T>()); 88 part_term += log(T(n) * (n + 1)) - boost::math::constants::ln_two<T>() - log(x); 89 part_term = exp(part_term); 90 } 91 else 92 { 93 sum = part_term * (n + 2 * x) / 2; 94 part_term *= (T(n) * (n + 1)) / 2; 95 part_term /= x; 96 } 97 // 98 // If the leading term is 0, so is the result: 99 // 100 if(sum == 0) 101 return sum; 102 103 for(unsigned k = 1;;) 104 { 105 term = part_term * boost::math::bernoulli_b2n<T>(k, pol); 106 sum += term; 107 // 108 // Normal termination condition: 109 // 110 if(fabs(term / sum) < tools::epsilon<T>()) 111 break; 112 // 113 // Increment our counter, and move part_term on to the next value: 114 // 115 ++k; 116 part_term *= T(n + 2 * k - 2) * (n - 1 + 2 * k); 117 part_term /= (2 * k - 1) * 2 * k; 118 part_term /= x_squared; 119 // 120 // Emergency get out termination condition: 121 // 122 if(k > policies::get_max_series_iterations<Policy>()) 123 { 124 return policies::raise_evaluation_error(function, "Series did not converge, closest value was %1%", sum, pol); 125 } 126 } 127 128 if((n - 1) & 1) 129 sum = -sum; 130 131 return sum; 132 } 133 134 template<class T, class Policy> polygamma_attransitionplus(const int n,const T & x,const Policy & pol,const char * function)135 T polygamma_attransitionplus(const int n, const T& x, const Policy& pol, const char* function) 136 { 137 // See: http://functions.wolfram.com/GammaBetaErf/PolyGamma2/16/01/01/0017/ 138 139 // Use N = (0.4 * digits) + (4 * n) for target value for x: 140 BOOST_MATH_STD_USING 141 const int d4d = static_cast<int>(0.4F * policies::digits_base10<T, Policy>()); 142 const int N = d4d + (4 * n); 143 const int m = n; 144 const int iter = N - itrunc(x); 145 146 if(iter > (int)policies::get_max_series_iterations<Policy>()) 147 return policies::raise_evaluation_error<T>(function, ("Exceeded maximum series evaluations evaluating at n = " + boost::lexical_cast<std::string>(n) + " and x = %1%").c_str(), x, pol); 148 149 const int minus_m_minus_one = -m - 1; 150 151 T z(x); 152 T sum0(0); 153 T z_plus_k_pow_minus_m_minus_one(0); 154 155 // Forward recursion to larger x, need to check for overflow first though: 156 if(log(z + iter) * minus_m_minus_one > -tools::log_max_value<T>()) 157 { 158 for(int k = 1; k <= iter; ++k) 159 { 160 z_plus_k_pow_minus_m_minus_one = pow(z, minus_m_minus_one); 161 sum0 += z_plus_k_pow_minus_m_minus_one; 162 z += 1; 163 } 164 sum0 *= boost::math::factorial<T>(n); 165 } 166 else 167 { 168 for(int k = 1; k <= iter; ++k) 169 { 170 T log_term = log(z) * minus_m_minus_one + boost::math::lgamma(T(n + 1), pol); 171 sum0 += exp(log_term); 172 z += 1; 173 } 174 } 175 if((n - 1) & 1) 176 sum0 = -sum0; 177 178 return sum0 + polygamma_atinfinityplus(n, z, pol, function); 179 } 180 181 template <class T, class Policy> polygamma_nearzero(int n,T x,const Policy & pol,const char * function)182 T polygamma_nearzero(int n, T x, const Policy& pol, const char* function) 183 { 184 BOOST_MATH_STD_USING 185 // 186 // If we take this expansion for polygamma: http://functions.wolfram.com/06.15.06.0003.02 187 // and substitute in this expression for polygamma(n, 1): http://functions.wolfram.com/06.15.03.0009.01 188 // we get an alternating series for polygamma when x is small in terms of zeta functions of 189 // integer arguments (which are easy to evaluate, at least when the integer is even). 190 // 191 // In order to avoid spurious overflow, save the n! term for later, and rescale at the end: 192 // 193 T scale = boost::math::factorial<T>(n, pol); 194 // 195 // "factorial_part" contains everything except the zeta function 196 // evaluations in each term: 197 // 198 T factorial_part = 1; 199 // 200 // "prefix" is what we'll be adding the accumulated sum to, it will 201 // be n! / z^(n+1), but since we're scaling by n! it's just 202 // 1 / z^(n+1) for now: 203 // 204 T prefix = pow(x, n + 1); 205 if(prefix == 0) 206 return boost::math::policies::raise_overflow_error<T>(function, 0, pol); 207 prefix = 1 / prefix; 208 // 209 // First term in the series is necessarily < zeta(2) < 2, so 210 // ignore the sum if it will have no effect on the result anyway: 211 // 212 if(prefix > 2 / policies::get_epsilon<T, Policy>()) 213 return ((n & 1) ? 1 : -1) * 214 (tools::max_value<T>() / prefix < scale ? policies::raise_overflow_error<T>(function, 0, pol) : prefix * scale); 215 // 216 // As this is an alternating series we could accelerate it using 217 // "Convergence Acceleration of Alternating Series", 218 // Henri Cohen, Fernando Rodriguez Villegas, and Don Zagier, Experimental Mathematics, 1999. 219 // In practice however, it appears not to make any difference to the number of terms 220 // required except in some edge cases which are filtered out anyway before we get here. 221 // 222 T sum = prefix; 223 for(unsigned k = 0;;) 224 { 225 // Get the k'th term: 226 T term = factorial_part * boost::math::zeta(T(k + n + 1), pol); 227 sum += term; 228 // Termination condition: 229 if(fabs(term) < fabs(sum * boost::math::policies::get_epsilon<T, Policy>())) 230 break; 231 // 232 // Move on k and factorial_part: 233 // 234 ++k; 235 factorial_part *= (-x * (n + k)) / k; 236 // 237 // Last chance exit: 238 // 239 if(k > policies::get_max_series_iterations<Policy>()) 240 return policies::raise_evaluation_error<T>(function, "Series did not converge, best value is %1%", sum, pol); 241 } 242 // 243 // We need to multiply by the scale, at each stage checking for overflow: 244 // 245 if(boost::math::tools::max_value<T>() / scale < sum) 246 return boost::math::policies::raise_overflow_error<T>(function, 0, pol); 247 sum *= scale; 248 return n & 1 ? sum : T(-sum); 249 } 250 251 // 252 // Helper function which figures out which slot our coefficient is in 253 // given an angle multiplier for the cosine term of power: 254 // 255 template <class Table> dereference_table(Table & table,unsigned row,unsigned power)256 typename Table::value_type::reference dereference_table(Table& table, unsigned row, unsigned power) 257 { 258 return table[row][power / 2]; 259 } 260 261 262 263 template <class T, class Policy> poly_cot_pi(int n,T x,T xc,const Policy & pol,const char * function)264 T poly_cot_pi(int n, T x, T xc, const Policy& pol, const char* function) 265 { 266 BOOST_MATH_STD_USING 267 // Return n'th derivative of cot(pi*x) at x, these are simply 268 // tabulated for up to n = 9, beyond that it is possible to 269 // calculate coefficients as follows: 270 // 271 // The general form of each derivative is: 272 // 273 // pi^n * SUM{k=0, n} C[k,n] * cos^k(pi * x) * csc^(n+1)(pi * x) 274 // 275 // With constant C[0,1] = -1 and all other C[k,n] = 0; 276 // Then for each k < n+1: 277 // C[k-1, n+1] -= k * C[k, n]; 278 // C[k+1, n+1] += (k-n-1) * C[k, n]; 279 // 280 // Note that there are many different ways of representing this derivative thanks to 281 // the many trigonometric identies available. In particular, the sum of powers of 282 // cosines could be replaced by a sum of cosine multiple angles, and indeed if you 283 // plug the derivative into Mathematica this is the form it will give. The two 284 // forms are related via the Chebeshev polynomials of the first kind and 285 // T_n(cos(x)) = cos(n x). The polynomial form has the great advantage that 286 // all the cosine terms are zero at half integer arguments - right where this 287 // function has it's minimum - thus avoiding cancellation error in this region. 288 // 289 // And finally, since every other term in the polynomials is zero, we can save 290 // space by only storing the non-zero terms. This greatly complexifies 291 // subscripting the tables in the calculation, but halves the storage space 292 // (and complexity for that matter). 293 // 294 T s = fabs(x) < fabs(xc) ? boost::math::sin_pi(x, pol) : boost::math::sin_pi(xc, pol); 295 T c = boost::math::cos_pi(x, pol); 296 switch(n) 297 { 298 case 1: 299 return -constants::pi<T, Policy>() / (s * s); 300 case 2: 301 { 302 return 2 * constants::pi<T, Policy>() * constants::pi<T, Policy>() * c / boost::math::pow<3>(s, pol); 303 } 304 case 3: 305 { 306 int P[] = { -2, -4 }; 307 return boost::math::pow<3>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<4>(s, pol); 308 } 309 case 4: 310 { 311 int P[] = { 16, 8 }; 312 return boost::math::pow<4>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<5>(s, pol); 313 } 314 case 5: 315 { 316 int P[] = { -16, -88, -16 }; 317 return boost::math::pow<5>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<6>(s, pol); 318 } 319 case 6: 320 { 321 int P[] = { 272, 416, 32 }; 322 return boost::math::pow<6>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<7>(s, pol); 323 } 324 case 7: 325 { 326 int P[] = { -272, -2880, -1824, -64 }; 327 return boost::math::pow<7>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<8>(s, pol); 328 } 329 case 8: 330 { 331 int P[] = { 7936, 24576, 7680, 128 }; 332 return boost::math::pow<8>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<9>(s, pol); 333 } 334 case 9: 335 { 336 int P[] = { -7936, -137216, -185856, -31616, -256 }; 337 return boost::math::pow<9>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<10>(s, pol); 338 } 339 case 10: 340 { 341 int P[] = { 353792, 1841152, 1304832, 128512, 512 }; 342 return boost::math::pow<10>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<11>(s, pol); 343 } 344 case 11: 345 { 346 int P[] = { -353792, -9061376, -21253376, -8728576, -518656, -1024}; 347 return boost::math::pow<11>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<12>(s, pol); 348 } 349 case 12: 350 { 351 int P[] = { 22368256, 175627264, 222398464, 56520704, 2084864, 2048 }; 352 return boost::math::pow<12>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<13>(s, pol); 353 } 354 #ifndef BOOST_NO_LONG_LONG 355 case 13: 356 { 357 long long P[] = { -22368256LL, -795300864LL, -2868264960LL, -2174832640LL, -357888000LL, -8361984LL, -4096 }; 358 return boost::math::pow<13>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<14>(s, pol); 359 } 360 case 14: 361 { 362 long long P[] = { 1903757312LL, 21016670208LL, 41731645440LL, 20261765120LL, 2230947840LL, 33497088LL, 8192 }; 363 return boost::math::pow<14>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<15>(s, pol); 364 } 365 case 15: 366 { 367 long long P[] = { -1903757312LL, -89702612992LL, -460858269696LL, -559148810240LL, -182172651520LL, -13754155008LL, -134094848LL, -16384 }; 368 return boost::math::pow<15>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<16>(s, pol); 369 } 370 case 16: 371 { 372 long long P[] = { 209865342976LL, 3099269660672LL, 8885192097792LL, 7048869314560LL, 1594922762240LL, 84134068224LL, 536608768LL, 32768 }; 373 return boost::math::pow<16>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<17>(s, pol); 374 } 375 case 17: 376 { 377 long long P[] = { -209865342976LL, -12655654469632LL, -87815735738368LL, -155964390375424LL, -84842998005760LL, -13684856848384LL, -511780323328LL, -2146926592LL, -65536 }; 378 return boost::math::pow<17>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<18>(s, pol); 379 } 380 case 18: 381 { 382 long long P[] = { 29088885112832LL, 553753414467584LL, 2165206642589696LL, 2550316668551168LL, 985278548541440LL, 115620218667008LL, 3100738912256LL, 8588754944LL, 131072 }; 383 return boost::math::pow<18>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<19>(s, pol); 384 } 385 case 19: 386 { 387 long long P[] = { -29088885112832LL, -2184860175433728LL, -19686087844429824LL, -48165109676113920LL, -39471306959486976LL, -11124607890751488LL, -965271355195392LL, -18733264797696LL, -34357248000LL, -262144 }; 388 return boost::math::pow<19>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<20>(s, pol); 389 } 390 case 20: 391 { 392 long long P[] = { 4951498053124096LL, 118071834535526400LL, 603968063567560704LL, 990081991141490688LL, 584901762421358592LL, 122829335169859584LL, 7984436548730880LL, 112949304754176LL, 137433710592LL, 524288 }; 393 return boost::math::pow<20>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<21>(s, pol); 394 } 395 #endif 396 } 397 398 // 399 // We'll have to compute the coefficients up to n, 400 // complexity is O(n^2) which we don't worry about for now 401 // as the values are computed once and then cached. 402 // However, if the final evaluation would have too many 403 // terms just bail out right away: 404 // 405 if((unsigned)n / 2u > policies::get_max_series_iterations<Policy>()) 406 return policies::raise_evaluation_error<T>(function, "The value of n is so large that we're unable to compute the result in reasonable time, best guess is %1%", 0, pol); 407 #ifdef BOOST_HAS_THREADS 408 static boost::detail::lightweight_mutex m; 409 boost::detail::lightweight_mutex::scoped_lock l(m); 410 #endif 411 static int digits = tools::digits<T>(); 412 static std::vector<std::vector<T> > table(1, std::vector<T>(1, T(-1))); 413 414 int current_digits = tools::digits<T>(); 415 416 if(digits != current_digits) 417 { 418 // Oh my... our precision has changed! 419 table = std::vector<std::vector<T> >(1, std::vector<T>(1, T(-1))); 420 digits = current_digits; 421 } 422 423 int index = n - 1; 424 425 if(index >= (int)table.size()) 426 { 427 for(int i = (int)table.size() - 1; i < index; ++i) 428 { 429 int offset = i & 1; // 1 if the first cos power is 0, otherwise 0. 430 int sin_order = i + 2; // order of the sin term 431 int max_cos_order = sin_order - 1; // largest order of the polynomial of cos terms 432 int max_columns = (max_cos_order - offset) / 2; // How many entries there are in the current row. 433 int next_offset = offset ? 0 : 1; 434 int next_max_columns = (max_cos_order + 1 - next_offset) / 2; // How many entries there will be in the next row 435 table.push_back(std::vector<T>(next_max_columns + 1, T(0))); 436 437 for(int column = 0; column <= max_columns; ++column) 438 { 439 int cos_order = 2 * column + offset; // order of the cosine term in entry "column" 440 BOOST_ASSERT(column < (int)table[i].size()); 441 BOOST_ASSERT((cos_order + 1) / 2 < (int)table[i + 1].size()); 442 table[i + 1][(cos_order + 1) / 2] += ((cos_order - sin_order) * table[i][column]) / (sin_order - 1); 443 if(cos_order) 444 table[i + 1][(cos_order - 1) / 2] += (-cos_order * table[i][column]) / (sin_order - 1); 445 } 446 } 447 448 } 449 T sum = boost::math::tools::evaluate_even_polynomial(&table[index][0], c, table[index].size()); 450 if(index & 1) 451 sum *= c; // First coefficient is order 1, and really an odd polynomial. 452 if(sum == 0) 453 return sum; 454 // 455 // The remaining terms are computed using logs since the powers and factorials 456 // get real large real quick: 457 // 458 T power_terms = n * log(boost::math::constants::pi<T>()); 459 if(s == 0) 460 return sum * boost::math::policies::raise_overflow_error<T>(function, 0, pol); 461 power_terms -= log(fabs(s)) * (n + 1); 462 power_terms += boost::math::lgamma(T(n)); 463 power_terms += log(fabs(sum)); 464 465 if(power_terms > boost::math::tools::log_max_value<T>()) 466 return sum * boost::math::policies::raise_overflow_error<T>(function, 0, pol); 467 468 return exp(power_terms) * ((s < 0) && ((n + 1) & 1) ? -1 : 1) * boost::math::sign(sum); 469 } 470 471 template <class T, class Policy> 472 struct polygamma_initializer 473 { 474 struct init 475 { initboost::math::detail::polygamma_initializer::init476 init() 477 { 478 // Forces initialization of our table of coefficients and mutex: 479 boost::math::polygamma(30, T(-2.5f), Policy()); 480 } force_instantiateboost::math::detail::polygamma_initializer::init481 void force_instantiate()const{} 482 }; 483 static const init initializer; force_instantiateboost::math::detail::polygamma_initializer484 static void force_instantiate() 485 { 486 initializer.force_instantiate(); 487 } 488 }; 489 490 template <class T, class Policy> 491 const typename polygamma_initializer<T, Policy>::init polygamma_initializer<T, Policy>::initializer; 492 493 template<class T, class Policy> polygamma_imp(const int n,T x,const Policy & pol)494 inline T polygamma_imp(const int n, T x, const Policy &pol) 495 { 496 BOOST_MATH_STD_USING 497 static const char* function = "boost::math::polygamma<%1%>(int, %1%)"; 498 polygamma_initializer<T, Policy>::initializer.force_instantiate(); 499 if(n < 0) 500 return policies::raise_domain_error<T>(function, "Order must be >= 0, but got %1%", static_cast<T>(n), pol); 501 if(x < 0) 502 { 503 if(floor(x) == x) 504 { 505 // 506 // Result is infinity if x is odd, and a pole error if x is even. 507 // 508 if(lltrunc(x) & 1) 509 return policies::raise_overflow_error<T>(function, 0, pol); 510 else 511 return policies::raise_pole_error<T>(function, "Evaluation at negative integer %1%", x, pol); 512 } 513 T z = 1 - x; 514 T result = polygamma_imp(n, z, pol) + constants::pi<T, Policy>() * poly_cot_pi(n, z, x, pol, function); 515 return n & 1 ? T(-result) : result; 516 } 517 // 518 // Limit for use of small-x-series is chosen 519 // so that the series doesn't go too divergent 520 // in the first few terms. Ordinarily this 521 // would mean setting the limit to ~ 1 / n, 522 // but we can tolerate a small amount of divergence: 523 // 524 T small_x_limit = (std::min)(T(T(5) / n), T(0.25f)); 525 if(x < small_x_limit) 526 { 527 return polygamma_nearzero(n, x, pol, function); 528 } 529 else if(x > 0.4F * policies::digits_base10<T, Policy>() + 4.0f * n) 530 { 531 return polygamma_atinfinityplus(n, x, pol, function); 532 } 533 else if(x == 1) 534 { 535 return (n & 1 ? 1 : -1) * boost::math::factorial<T>(n, pol) * boost::math::zeta(T(n + 1), pol); 536 } 537 else if(x == 0.5f) 538 { 539 T result = (n & 1 ? 1 : -1) * boost::math::factorial<T>(n, pol) * boost::math::zeta(T(n + 1), pol); 540 if(fabs(result) >= ldexp(tools::max_value<T>(), -n - 1)) 541 return boost::math::sign(result) * policies::raise_overflow_error<T>(function, 0, pol); 542 result *= ldexp(T(1), n + 1) - 1; 543 return result; 544 } 545 else 546 { 547 return polygamma_attransitionplus(n, x, pol, function); 548 } 549 } 550 551 } } } // namespace boost::math::detail 552 553 #ifdef _MSC_VER 554 #pragma warning(pop) 555 #endif 556 557 #endif // _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ 558 559