1 /////////////////////////////////////////////////////////////////////////////// 2 // Copyright 2018 John Maddock 3 // Distributed under the Boost 4 // Software License, Version 1.0. (See accompanying file 5 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 6 7 #ifndef BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_ 8 #define BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_ 9 10 #ifndef BOOST_MATH_PFQ_MAX_B_TERMS 11 # define BOOST_MATH_PFQ_MAX_B_TERMS 5 12 #endif 13 14 #include <array> 15 #include <cstdint> 16 #include <boost/math/special_functions/detail/hypergeometric_series.hpp> 17 18 namespace boost { namespace math { namespace detail { 19 20 template <class Seq, class Real> set_crossover_locations(const Seq & aj,const Seq & bj,const Real & z,unsigned int * crossover_locations)21 unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations) 22 { 23 BOOST_MATH_STD_USING 24 unsigned N_terms = 0; 25 26 if(aj.size() == 1 && bj.size() == 1) 27 { 28 // 29 // For 1F1 we can work out where the peaks in the series occur, 30 // which is to say when: 31 // 32 // (a + k)z / (k(b + k)) == +-1 33 // 34 // Then we are at either a maxima or a minima in the series, and the 35 // last such point must be a maxima since the series is globally convergent. 36 // Potentially then we are solving 2 quadratic equations and have up to 4 37 // solutions, any solutions which are complex or negative are discarded, 38 // leaving us with 4 options: 39 // 40 // 0 solutions: The series is directly convergent. 41 // 1 solution : The series diverges to a maxima before converging. 42 // 2 solutions: The series is initially convergent, followed by divergence to a maxima before final convergence. 43 // 3 solutions: The series diverges to a maxima, converges to a minima before diverging again to a second maxima before final convergence. 44 // 4 solutions: The series converges to a minima before diverging to a maxima, converging to a minima, diverging to a second maxima and then converging. 45 // 46 // The first 2 situations are adequately handled by direct series evaluation, while the 2,3 and 4 solutions are not. 47 // 48 Real a = *aj.begin(); 49 Real b = *bj.begin(); 50 Real sq = 4 * a * z + b * b - 2 * b * z + z * z; 51 if (sq >= 0) 52 { 53 Real t = (-sqrt(sq) - b + z) / 2; 54 if (t >= 0) 55 { 56 crossover_locations[N_terms] = itrunc(t); 57 ++N_terms; 58 } 59 t = (sqrt(sq) - b + z) / 2; 60 if (t >= 0) 61 { 62 crossover_locations[N_terms] = itrunc(t); 63 ++N_terms; 64 } 65 } 66 sq = -4 * a * z + b * b + 2 * b * z + z * z; 67 if (sq >= 0) 68 { 69 Real t = (-sqrt(sq) - b - z) / 2; 70 if (t >= 0) 71 { 72 crossover_locations[N_terms] = itrunc(t); 73 ++N_terms; 74 } 75 t = (sqrt(sq) - b - z) / 2; 76 if (t >= 0) 77 { 78 crossover_locations[N_terms] = itrunc(t); 79 ++N_terms; 80 } 81 } 82 std::sort(crossover_locations, crossover_locations + N_terms, std::less<Real>()); 83 // 84 // Now we need to discard every other terms, as these are the minima: 85 // 86 switch (N_terms) 87 { 88 case 0: 89 case 1: 90 break; 91 case 2: 92 crossover_locations[0] = crossover_locations[1]; 93 --N_terms; 94 break; 95 case 3: 96 crossover_locations[1] = crossover_locations[2]; 97 --N_terms; 98 break; 99 case 4: 100 crossover_locations[0] = crossover_locations[1]; 101 crossover_locations[1] = crossover_locations[3]; 102 N_terms -= 2; 103 break; 104 } 105 } 106 else 107 { 108 unsigned n = 0; 109 for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n) 110 { 111 crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1; 112 } 113 std::sort(crossover_locations, crossover_locations + bj.size(), std::less<Real>()); 114 N_terms = (unsigned)bj.size(); 115 } 116 return N_terms; 117 } 118 119 template <class Seq, class Real, class Policy, class Terminal> hypergeometric_pFq_checked_series_impl(const Seq & aj,const Seq & bj,const Real & z,const Policy & pol,const Terminal & termination,long long & log_scale)120 std::pair<Real, Real> hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, long long& log_scale) 121 { 122 BOOST_MATH_STD_USING 123 Real result = 1; 124 Real abs_result = 1; 125 Real term = 1; 126 Real term0 = 0; 127 Real tol = boost::math::policies::get_epsilon<Real, Policy>(); 128 std::uintmax_t k = 0; 129 Real upper_limit(sqrt(boost::math::tools::max_value<Real>())), diff; 130 Real lower_limit(1 / upper_limit); 131 long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>()) - 2; 132 Real scaling_factor = exp(Real(log_scaling_factor)); 133 Real term_m1; 134 long long local_scaling = 0; 135 136 if ((aj.size() == 1) && (bj.size() == 0)) 137 { 138 if (fabs(z) > 1) 139 { 140 if ((z > 0) && (floor(*aj.begin()) != *aj.begin())) 141 { 142 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == 1 and q == 0 and |z| > 1, result is imaginary", z, pol); 143 return std::make_pair(r, r); 144 } 145 std::pair<Real, Real> r = hypergeometric_pFq_checked_series_impl(aj, bj, Real(1 / z), pol, termination, log_scale); 146 Real mul = pow(-z, -*aj.begin()); 147 r.first *= mul; 148 r.second *= mul; 149 return r; 150 } 151 } 152 153 if (aj.size() > bj.size()) 154 { 155 if (aj.size() == bj.size() + 1) 156 { 157 if (fabs(z) > 1) 158 { 159 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| > 1, series does not converge", z, pol); 160 return std::make_pair(r, r); 161 } 162 if (fabs(z) == 1) 163 { 164 Real s = 0; 165 for (auto i = bj.begin(); i != bj.end(); ++i) 166 s += *i; 167 for (auto i = aj.begin(); i != aj.end(); ++i) 168 s -= *i; 169 if ((z == 1) && (s <= 0)) 170 { 171 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| == 1, in a situation where the series does not converge", z, pol); 172 return std::make_pair(r, r); 173 } 174 if ((z == -1) && (s <= -1)) 175 { 176 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p == q+1 and |z| == 1, in a situation where the series does not converge", z, pol); 177 return std::make_pair(r, r); 178 } 179 } 180 } 181 else 182 { 183 Real r = policies::raise_domain_error("boost::math::hypergeometric_pFq", "Got p > q+1, series does not converge", z, pol); 184 return std::make_pair(r, r); 185 } 186 } 187 188 while (!termination(k)) 189 { 190 for (auto ai = aj.begin(); ai != aj.end(); ++ai) 191 { 192 term *= *ai + k; 193 } 194 if (term == 0) 195 { 196 // There is a negative integer in the aj's: 197 return std::make_pair(result, abs_result); 198 } 199 for (auto bi = bj.begin(); bi != bj.end(); ++bi) 200 { 201 if (*bi + k == 0) 202 { 203 // The series is undefined: 204 result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol); 205 return std::make_pair(result, result); 206 } 207 term /= *bi + k; 208 } 209 term *= z; 210 ++k; 211 term /= k; 212 //std::cout << k << " " << *bj.begin() + k << " " << result << " " << term << /*" " << term_at_k(*aj.begin(), *bj.begin(), z, k, pol) <<*/ std::endl; 213 result += term; 214 abs_result += abs(term); 215 //std::cout << "k = " << k << " term = " << term * exp(log_scale) << " result = " << result * exp(log_scale) << " abs_result = " << abs_result * exp(log_scale) << std::endl; 216 217 // 218 // Rescaling: 219 // 220 if (fabs(abs_result) >= upper_limit) 221 { 222 abs_result /= scaling_factor; 223 result /= scaling_factor; 224 term /= scaling_factor; 225 log_scale += log_scaling_factor; 226 local_scaling += log_scaling_factor; 227 } 228 if (fabs(abs_result) < lower_limit) 229 { 230 abs_result *= scaling_factor; 231 result *= scaling_factor; 232 term *= scaling_factor; 233 log_scale -= log_scaling_factor; 234 local_scaling -= log_scaling_factor; 235 } 236 237 if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term))) 238 break; 239 if (abs_result * tol > abs(result)) 240 { 241 // We have no correct bits in the result... just give up! 242 result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the reuslt are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol); 243 return std::make_pair(result, result); 244 } 245 term0 = term; 246 } 247 //std::cout << "result = " << result << std::endl; 248 //std::cout << "local_scaling = " << local_scaling << std::endl; 249 //std::cout << "Norm result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(result) * exp(boost::multiprecision::mpfr_float_50(local_scaling)) << std::endl; 250 // 251 // We have to be careful when one of the b's crosses the origin: 252 // 253 if(bj.size() > BOOST_MATH_PFQ_MAX_B_TERMS) 254 policies::raise_domain_error<Real>("boost::math::hypergeometric_pFq<%1%>(Seq, Seq, %1%)", 255 "The number of b terms must be less than the value of BOOST_MATH_PFQ_MAX_B_TERMS (" BOOST_STRINGIZE(BOOST_MATH_PFQ_MAX_B_TERMS) "), but got %1%.", 256 Real(bj.size()), pol); 257 258 unsigned crossover_locations[BOOST_MATH_PFQ_MAX_B_TERMS]; 259 260 unsigned N_crossovers = set_crossover_locations(aj, bj, z, crossover_locations); 261 262 bool terminate = false; // Set to true if one of the a's passes through the origin and terminates the series. 263 264 for (unsigned n = 0; n < N_crossovers; ++n) 265 { 266 if (k < crossover_locations[n]) 267 { 268 for (auto ai = aj.begin(); ai != aj.end(); ++ai) 269 { 270 if ((*ai < 0) && (floor(*ai) == *ai) && (*ai > crossover_locations[n])) 271 return std::make_pair(result, abs_result); // b's will never cross the origin! 272 } 273 // 274 // local results: 275 // 276 Real loop_result = 0; 277 Real loop_abs_result = 0; 278 long long loop_scale = 0; 279 // 280 // loop_error_scale will be used to increase the size of the error 281 // estimate (absolute sum), based on the errors inherent in calculating 282 // the pochhammer symbols. 283 // 284 Real loop_error_scale = 0; 285 //boost::multiprecision::mpfi_float err_est = 0; 286 // 287 // b hasn't crossed the origin yet and the series may spring back into life at that point 288 // so we need to jump forward to that term and then evaluate forwards and backwards from there: 289 // 290 unsigned s = crossover_locations[n]; 291 std::uintmax_t backstop = k; 292 long long s1(1), s2(1); 293 term = 0; 294 for (auto ai = aj.begin(); ai != aj.end(); ++ai) 295 { 296 if ((floor(*ai) == *ai) && (*ai < 0) && (-*ai <= s)) 297 { 298 // One of the a terms has passed through zero and terminated the series: 299 terminate = true; 300 break; 301 } 302 else 303 { 304 int ls = 1; 305 Real p = log_pochhammer(*ai, s, pol, &ls); 306 s1 *= ls; 307 term += p; 308 loop_error_scale = (std::max)(p, loop_error_scale); 309 //err_est += boost::multiprecision::mpfi_float(p); 310 } 311 } 312 //std::cout << "term = " << term << std::endl; 313 if (terminate) 314 break; 315 for (auto bi = bj.begin(); bi != bj.end(); ++bi) 316 { 317 int ls = 1; 318 Real p = log_pochhammer(*bi, s, pol, &ls); 319 s2 *= ls; 320 term -= p; 321 loop_error_scale = (std::max)(p, loop_error_scale); 322 //err_est -= boost::multiprecision::mpfi_float(p); 323 } 324 //std::cout << "term = " << term << std::endl; 325 Real p = lgamma(Real(s + 1), pol); 326 term -= p; 327 loop_error_scale = (std::max)(p, loop_error_scale); 328 //err_est -= boost::multiprecision::mpfi_float(p); 329 p = s * log(fabs(z)); 330 term += p; 331 loop_error_scale = (std::max)(p, loop_error_scale); 332 //err_est += boost::multiprecision::mpfi_float(p); 333 //err_est = exp(err_est); 334 //std::cout << err_est << std::endl; 335 // 336 // Convert loop_error scale to the absolute error 337 // in term after exp is applied: 338 // 339 loop_error_scale *= tools::epsilon<Real>(); 340 // 341 // Convert to relative error after exp: 342 // 343 loop_error_scale = fabs(expm1(loop_error_scale, pol)); 344 // 345 // Convert to multiplier for the error term: 346 // 347 loop_error_scale /= tools::epsilon<Real>(); 348 349 if (z < 0) 350 s1 *= (s & 1 ? -1 : 1); 351 352 if (term <= tools::log_min_value<Real>()) 353 { 354 // rescale if we can: 355 long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2); 356 term -= scale; 357 loop_scale += scale; 358 } 359 if (term > 10) 360 { 361 int scale = itrunc(floor(term)); 362 term -= scale; 363 loop_scale += scale; 364 } 365 //std::cout << "term = " << term << std::endl; 366 term = s1 * s2 * exp(term); 367 //std::cout << "term = " << term << std::endl; 368 //std::cout << "loop_scale = " << loop_scale << std::endl; 369 k = s; 370 term0 = term; 371 long long saved_loop_scale = loop_scale; 372 bool terms_are_growing = true; 373 bool trivial_small_series_check = false; 374 do 375 { 376 loop_result += term; 377 loop_abs_result += fabs(term); 378 //std::cout << "k = " << k << " term = " << term * exp(loop_scale) << " result = " << loop_result * exp(loop_scale) << " abs_result = " << loop_abs_result * exp(loop_scale) << std::endl; 379 if (fabs(loop_result) >= upper_limit) 380 { 381 loop_result /= scaling_factor; 382 loop_abs_result /= scaling_factor; 383 term /= scaling_factor; 384 loop_scale += log_scaling_factor; 385 } 386 if (fabs(loop_result) < lower_limit) 387 { 388 loop_result *= scaling_factor; 389 loop_abs_result *= scaling_factor; 390 term *= scaling_factor; 391 loop_scale -= log_scaling_factor; 392 } 393 term_m1 = term; 394 for (auto ai = aj.begin(); ai != aj.end(); ++ai) 395 { 396 term *= *ai + k; 397 } 398 if (term == 0) 399 { 400 // There is a negative integer in the aj's: 401 return std::make_pair(result, abs_result); 402 } 403 for (auto bi = bj.begin(); bi != bj.end(); ++bi) 404 { 405 if (*bi + k == 0) 406 { 407 // The series is undefined: 408 result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol); 409 return std::make_pair(result, result); 410 } 411 term /= *bi + k; 412 } 413 term *= z / (k + 1); 414 415 ++k; 416 diff = fabs(term / loop_result); 417 terms_are_growing = fabs(term) > fabs(term_m1); 418 if (!trivial_small_series_check && !terms_are_growing) 419 { 420 // 421 // Now that we have started to converge, check to see if the value of 422 // this local sum is trivially small compared to the result. If so 423 // abort this part of the series. 424 // 425 trivial_small_series_check = true; 426 Real d; 427 if (loop_scale > local_scaling) 428 { 429 long long rescale = local_scaling - loop_scale; 430 if (rescale < tools::log_min_value<Real>()) 431 d = 1; // arbitrary value, we want to keep going 432 else 433 d = fabs(term / (result * exp(Real(rescale)))); 434 } 435 else 436 { 437 long long rescale = loop_scale - local_scaling; 438 if (rescale < tools::log_min_value<Real>()) 439 d = 0; // terminate this loop 440 else 441 d = fabs(term * exp(Real(rescale)) / result); 442 } 443 if (d < boost::math::policies::get_epsilon<Real, Policy>()) 444 break; 445 } 446 } while (!termination(k - s) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || terms_are_growing)); 447 448 //std::cout << "Norm loop result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(loop_result)* exp(boost::multiprecision::mpfr_float_50(loop_scale)) << std::endl; 449 // 450 // We now need to combine the results of the first series summation with whatever 451 // local results we have now. First though, rescale abs_result by loop_error_scale 452 // to factor in the error in the pochhammer terms at the start of this block: 453 // 454 std::uintmax_t next_backstop = k; 455 loop_abs_result += loop_error_scale * fabs(loop_result); 456 if (loop_scale > local_scaling) 457 { 458 // 459 // Need to shrink previous result: 460 // 461 long long rescale = local_scaling - loop_scale; 462 local_scaling = loop_scale; 463 log_scale -= rescale; 464 Real ex = exp(Real(rescale)); 465 result *= ex; 466 abs_result *= ex; 467 result += loop_result; 468 abs_result += loop_abs_result; 469 } 470 else if (local_scaling > loop_scale) 471 { 472 // 473 // Need to shrink local result: 474 // 475 long long rescale = loop_scale - local_scaling; 476 Real ex = exp(Real(rescale)); 477 loop_result *= ex; 478 loop_abs_result *= ex; 479 result += loop_result; 480 abs_result += loop_abs_result; 481 } 482 else 483 { 484 result += loop_result; 485 abs_result += loop_abs_result; 486 } 487 // 488 // Now go backwards as well: 489 // 490 k = s; 491 term = term0; 492 loop_result = 0; 493 loop_abs_result = 0; 494 loop_scale = saved_loop_scale; 495 trivial_small_series_check = false; 496 do 497 { 498 --k; 499 if (k == backstop) 500 break; 501 term_m1 = term; 502 for (auto ai = aj.begin(); ai != aj.end(); ++ai) 503 { 504 term /= *ai + k; 505 } 506 for (auto bi = bj.begin(); bi != bj.end(); ++bi) 507 { 508 if (*bi + k == 0) 509 { 510 // The series is undefined: 511 result = boost::math::policies::raise_domain_error("boost::math::hypergeometric_pFq<%1%>", "One of the b values was the negative integer %1%", *bi, pol); 512 return std::make_pair(result, result); 513 } 514 term *= *bi + k; 515 } 516 term *= (k + 1) / z; 517 loop_result += term; 518 loop_abs_result += fabs(term); 519 520 if (!trivial_small_series_check && (fabs(term) < fabs(term_m1))) 521 { 522 // 523 // Now that we have started to converge, check to see if the value of 524 // this local sum is trivially small compared to the result. If so 525 // abort this part of the series. 526 // 527 trivial_small_series_check = true; 528 Real d; 529 if (loop_scale > local_scaling) 530 { 531 long long rescale = local_scaling - loop_scale; 532 if (rescale < tools::log_min_value<Real>()) 533 d = 1; // keep going 534 else 535 d = fabs(term / (result * exp(Real(rescale)))); 536 } 537 else 538 { 539 long long rescale = loop_scale - local_scaling; 540 if (rescale < tools::log_min_value<Real>()) 541 d = 0; // stop, underflow 542 else 543 d = fabs(term * exp(Real(rescale)) / result); 544 } 545 if (d < boost::math::policies::get_epsilon<Real, Policy>()) 546 break; 547 } 548 549 //std::cout << "k = " << k << " result = " << result << " abs_result = " << abs_result << std::endl; 550 if (fabs(loop_result) >= upper_limit) 551 { 552 loop_result /= scaling_factor; 553 loop_abs_result /= scaling_factor; 554 term /= scaling_factor; 555 loop_scale += log_scaling_factor; 556 } 557 if (fabs(loop_result) < lower_limit) 558 { 559 loop_result *= scaling_factor; 560 loop_abs_result *= scaling_factor; 561 term *= scaling_factor; 562 loop_scale -= log_scaling_factor; 563 } 564 diff = fabs(term / loop_result); 565 } while (!termination(s - k) && ((diff > boost::math::policies::get_epsilon<Real, Policy>()) || (fabs(term) > fabs(term_m1)))); 566 567 //std::cout << "Norm loop result = " << std::setprecision(35) << boost::multiprecision::mpfr_float_50(loop_result)* exp(boost::multiprecision::mpfr_float_50(loop_scale)) << std::endl; 568 // 569 // We now need to combine the results of the first series summation with whatever 570 // local results we have now. First though, rescale abs_result by loop_error_scale 571 // to factor in the error in the pochhammer terms at the start of this block: 572 // 573 loop_abs_result += loop_error_scale * fabs(loop_result); 574 // 575 if (loop_scale > local_scaling) 576 { 577 // 578 // Need to shrink previous result: 579 // 580 long long rescale = local_scaling - loop_scale; 581 local_scaling = loop_scale; 582 log_scale -= rescale; 583 Real ex = exp(Real(rescale)); 584 result *= ex; 585 abs_result *= ex; 586 result += loop_result; 587 abs_result += loop_abs_result; 588 } 589 else if (local_scaling > loop_scale) 590 { 591 // 592 // Need to shrink local result: 593 // 594 long long rescale = loop_scale - local_scaling; 595 Real ex = exp(Real(rescale)); 596 loop_result *= ex; 597 loop_abs_result *= ex; 598 result += loop_result; 599 abs_result += loop_abs_result; 600 } 601 else 602 { 603 result += loop_result; 604 abs_result += loop_abs_result; 605 } 606 // 607 // Reset k to the largest k we reached 608 // 609 k = next_backstop; 610 } 611 } 612 613 return std::make_pair(result, abs_result); 614 } 615 616 struct iteration_terminator 617 { iteration_terminatorboost::math::detail::iteration_terminator618 iteration_terminator(std::uintmax_t i) : m(i) {} 619 operator ()boost::math::detail::iteration_terminator620 bool operator()(std::uintmax_t v) const { return v >= m; } 621 622 std::uintmax_t m; 623 }; 624 625 template <class Seq, class Real, class Policy> 626 Real hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, long long& log_scale) 627 { 628 BOOST_MATH_STD_USING 629 iteration_terminator term(boost::math::policies::get_max_series_iterations<Policy>()); 630 std::pair<Real, Real> result = hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, term, log_scale); 631 // 632 // Check to see how many digits we've lost, if it's more than half, raise an evaluation error - 633 // this is an entirely arbitrary cut off, but not unreasonable. 634 // 635 if (result.second * sqrt(boost::math::policies::get_epsilon<Real, Policy>()) > abs(result.first)) 636 { 637 return boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that fewer than half the bits in the result are correct, last result was %1%", Real(result.first * exp(Real(log_scale))), pol); 638 } 639 return result.first; 640 } 641 642 template <class Real, class Policy> hypergeometric_1F1_checked_series_impl(const Real & a,const Real & b,const Real & z,const Policy & pol,long long & log_scale)643 inline Real hypergeometric_1F1_checked_series_impl(const Real& a, const Real& b, const Real& z, const Policy& pol, long long& log_scale) 644 { 645 std::array<Real, 1> aj = { a }; 646 std::array<Real, 1> bj = { b }; 647 return hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, log_scale); 648 } 649 650 } } } // namespaces 651 652 #endif // BOOST_HYPERGEOMETRIC_PFQ_SERIES_HPP_ 653