]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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 | ||
1e59de90 TL |
14 | #include <array> |
15 | #include <cstdint> | |
92f5a8d4 TL |
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> | |
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. | |
f67539c2 | 41 | // 1 solution : The series diverges to a maxima before converging. |
92f5a8d4 TL |
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 | // | |
f67539c2 | 46 | // The first 2 situations are adequately handled by direct series evaluation, while the 2,3 and 4 solutions are not. |
92f5a8d4 TL |
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> | |
1e59de90 | 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) |
92f5a8d4 TL |
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>(); | |
1e59de90 | 128 | std::uintmax_t k = 0; |
92f5a8d4 TL |
129 | Real upper_limit(sqrt(boost::math::tools::max_value<Real>())), diff; |
130 | Real lower_limit(1 / upper_limit); | |
1e59de90 | 131 | long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value<Real>()) - 2; |
92f5a8d4 TL |
132 | Real scaling_factor = exp(Real(log_scaling_factor)); |
133 | Real term_m1; | |
1e59de90 | 134 | long long local_scaling = 0; |
92f5a8d4 TL |
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; | |
1e59de90 | 278 | long long loop_scale = 0; |
92f5a8d4 TL |
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]; | |
1e59de90 TL |
291 | std::uintmax_t backstop = k; |
292 | long long s1(1), s2(1); | |
92f5a8d4 TL |
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 | // | |
20effc67 | 343 | loop_error_scale = fabs(expm1(loop_error_scale, pol)); |
92f5a8d4 TL |
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: | |
1e59de90 | 355 | long long scale = lltrunc(floor(term - tools::log_min_value<Real>()) - 2); |
92f5a8d4 TL |
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; | |
1e59de90 | 371 | long long saved_loop_scale = loop_scale; |
92f5a8d4 TL |
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 | { | |
1e59de90 | 429 | long long rescale = local_scaling - loop_scale; |
92f5a8d4 | 430 | if (rescale < tools::log_min_value<Real>()) |
f67539c2 | 431 | d = 1; // arbitrary value, we want to keep going |
92f5a8d4 TL |
432 | else |
433 | d = fabs(term / (result * exp(Real(rescale)))); | |
434 | } | |
435 | else | |
436 | { | |
1e59de90 | 437 | long long rescale = loop_scale - local_scaling; |
92f5a8d4 TL |
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 | // | |
1e59de90 | 454 | std::uintmax_t next_backstop = k; |
92f5a8d4 TL |
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 | // | |
1e59de90 | 461 | long long rescale = local_scaling - loop_scale; |
92f5a8d4 TL |
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 | // | |
1e59de90 | 475 | long long rescale = loop_scale - local_scaling; |
92f5a8d4 TL |
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 | { | |
1e59de90 | 531 | long long rescale = local_scaling - loop_scale; |
92f5a8d4 TL |
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 | { | |
1e59de90 | 539 | long long rescale = loop_scale - local_scaling; |
92f5a8d4 TL |
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 | // | |
1e59de90 | 580 | long long rescale = local_scaling - loop_scale; |
92f5a8d4 TL |
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 | // | |
1e59de90 | 594 | long long rescale = loop_scale - local_scaling; |
92f5a8d4 TL |
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 | { | |
1e59de90 | 618 | iteration_terminator(std::uintmax_t i) : m(i) {} |
92f5a8d4 | 619 | |
1e59de90 | 620 | bool operator()(std::uintmax_t v) const { return v >= m; } |
92f5a8d4 | 621 | |
1e59de90 | 622 | std::uintmax_t m; |
92f5a8d4 TL |
623 | }; |
624 | ||
625 | template <class Seq, class Real, class Policy> | |
1e59de90 | 626 | Real hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, long long& log_scale) |
92f5a8d4 TL |
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> | |
1e59de90 | 643 | inline Real hypergeometric_1F1_checked_series_impl(const Real& a, const Real& b, const Real& z, const Policy& pol, long long& log_scale) |
92f5a8d4 | 644 | { |
1e59de90 TL |
645 | std::array<Real, 1> aj = { a }; |
646 | std::array<Real, 1> bj = { b }; | |
92f5a8d4 TL |
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_ |