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