]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | |
2 | /////////////////////////////////////////////////////////////////////////////// | |
3 | // Copyright Christopher Kormanyos 2013 - 2014. | |
4 | // Copyright John Maddock 2013. | |
5 | // Distributed under the Boost Software License, | |
6 | // Version 1.0. (See accompanying file LICENSE_1_0.txt | |
7 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
8 | // | |
9 | // This work is based on an earlier work: | |
10 | // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations", | |
11 | // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469 | |
12 | // | |
13 | ||
14 | #include <algorithm> | |
15 | #include <cstdint> | |
16 | #include <deque> | |
17 | #include <functional> | |
18 | #include <iostream> | |
19 | #include <limits> | |
20 | #include <numeric> | |
21 | #include <vector> | |
22 | #include <boost/math/constants/constants.hpp> | |
23 | #include <boost/noncopyable.hpp> | |
24 | ||
25 | //#define USE_CPP_BIN_FLOAT | |
26 | #define USE_CPP_DEC_FLOAT | |
27 | //#define USE_MPFR | |
28 | ||
29 | #if !defined(DIGIT_COUNT) | |
30 | #define DIGIT_COUNT 100 | |
31 | #endif | |
32 | ||
33 | #if !defined(BOOST_NO_CXX11_HDR_CHRONO) | |
34 | #include <chrono> | |
35 | #define STD_CHRONO std::chrono | |
36 | #else | |
37 | #include <boost/chrono.hpp> | |
38 | #define STD_CHRONO boost::chrono | |
39 | #endif | |
40 | ||
41 | #if defined(USE_CPP_BIN_FLOAT) | |
42 | #include <boost/multiprecision/cpp_bin_float.hpp> | |
43 | typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<DIGIT_COUNT + 10> > mp_type; | |
44 | #elif defined(USE_CPP_DEC_FLOAT) | |
45 | #include <boost/multiprecision/cpp_dec_float.hpp> | |
46 | typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<DIGIT_COUNT + 10> > mp_type; | |
47 | #elif defined(USE_MPFR) | |
48 | #include <boost/multiprecision/mpfr.hpp> | |
49 | typedef boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<DIGIT_COUNT + 10> > mp_type; | |
50 | #else | |
51 | #error no multiprecision floating type is defined | |
52 | #endif | |
53 | ||
54 | template <class clock_type> | |
55 | struct stopwatch | |
56 | { | |
57 | public: | |
58 | typedef typename clock_type::duration duration_type; | |
59 | ||
60 | stopwatch() : m_start(clock_type::now()) { } | |
61 | ||
62 | stopwatch(const stopwatch& other) : m_start(other.m_start) { } | |
63 | ||
64 | stopwatch& operator=(const stopwatch& other) | |
65 | { | |
66 | m_start = other.m_start; | |
67 | return *this; | |
68 | } | |
69 | ||
70 | ~stopwatch() { } | |
71 | ||
72 | duration_type elapsed() const | |
73 | { | |
74 | return (clock_type::now() - m_start); | |
75 | } | |
76 | ||
77 | void reset() | |
78 | { | |
79 | m_start = clock_type::now(); | |
80 | } | |
81 | ||
82 | private: | |
83 | typename clock_type::time_point m_start; | |
84 | }; | |
85 | ||
86 | namespace my_math | |
87 | { | |
88 | template<class T> T chebyshev_t(const std::int32_t n, const T& x); | |
89 | ||
90 | template<class T> T chebyshev_t(const std::uint32_t n, const T& x, std::vector<T>* vp); | |
91 | ||
92 | template<class T> bool isneg(const T& x) { return (x < T(0)); } | |
93 | ||
94 | template<class T> const T& zero() { static const T value_zero(0); return value_zero; } | |
95 | template<class T> const T& one () { static const T value_one (1); return value_one; } | |
96 | template<class T> const T& two () { static const T value_two (2); return value_two; } | |
97 | } | |
98 | ||
99 | namespace orthogonal_polynomial_series | |
100 | { | |
101 | template<typename T> static inline T orthogonal_polynomial_template(const T& x, const std::uint32_t n, std::vector<T>* const vp = static_cast<std::vector<T>*>(0u)) | |
102 | { | |
103 | // Compute the value of an orthogonal chebyshev polinomial. | |
104 | // Use stable upward recursion. | |
105 | ||
106 | if(vp != nullptr) | |
107 | { | |
108 | vp->clear(); | |
109 | vp->reserve(static_cast<std::size_t>(n + 1u)); | |
110 | } | |
111 | ||
112 | T y0 = my_math::one<T>(); | |
113 | ||
114 | if(vp != nullptr) { vp->push_back(y0); } | |
115 | ||
116 | if(n == static_cast<std::uint32_t>(0u)) | |
117 | { | |
118 | return y0; | |
119 | } | |
120 | ||
121 | T y1 = x; | |
122 | ||
123 | if(vp != nullptr) { vp->push_back(y1); } | |
124 | ||
125 | if(n == static_cast<std::uint32_t>(1u)) | |
126 | { | |
127 | return y1; | |
128 | } | |
129 | ||
130 | T a = my_math::two <T>(); | |
131 | T b = my_math::zero<T>(); | |
132 | T c = my_math::one <T>(); | |
133 | ||
134 | T yk; | |
135 | ||
136 | // Calculate higher orders using the recurrence relation. | |
137 | // The direction of stability is upward recursion. | |
138 | for(std::int32_t k = static_cast<std::int32_t>(2); k <= static_cast<std::int32_t>(n); ++k) | |
139 | { | |
140 | yk = (((a * x) + b) * y1) - (c * y0); | |
141 | ||
142 | y0 = y1; | |
143 | y1 = yk; | |
144 | ||
145 | if(vp != nullptr) { vp->push_back(yk); } | |
146 | } | |
147 | ||
148 | return yk; | |
149 | } | |
150 | } | |
151 | ||
152 | template<class T> T my_math::chebyshev_t(const std::int32_t n, const T& x) | |
153 | { | |
154 | if(my_math::isneg(x)) | |
155 | { | |
156 | const bool b_negate = ((n % static_cast<std::int32_t>(2)) != static_cast<std::int32_t>(0)); | |
157 | ||
158 | const T y = chebyshev_t(n, -x); | |
159 | ||
160 | return (!b_negate ? y : -y); | |
161 | } | |
162 | ||
163 | if(n < static_cast<std::int32_t>(0)) | |
164 | { | |
165 | const std::int32_t nn = static_cast<std::int32_t>(-n); | |
166 | ||
167 | return chebyshev_t(nn, x); | |
168 | } | |
169 | else | |
170 | { | |
171 | return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::uint32_t>(n)); | |
172 | } | |
173 | } | |
174 | ||
175 | template<class T> T my_math::chebyshev_t(const std::uint32_t n, const T& x, std::vector<T>* const vp) { return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::int32_t>(n), vp); } | |
176 | ||
177 | namespace util | |
178 | { | |
179 | template <class T> float digit_scale() | |
180 | { | |
181 | const int d = ((std::max)(std::numeric_limits<T>::digits10, 15)); | |
182 | return static_cast<float>(d) / 300.0F; | |
183 | } | |
184 | } | |
185 | ||
186 | namespace examples | |
187 | { | |
188 | namespace nr_006 | |
189 | { | |
190 | template<typename T> class hypergeometric_pfq_base : private boost::noncopyable | |
191 | { | |
192 | public: | |
193 | virtual ~hypergeometric_pfq_base() { } | |
194 | ||
195 | virtual void ccoef() const = 0; | |
196 | ||
197 | virtual T series() const | |
198 | { | |
199 | using my_math::chebyshev_t; | |
200 | ||
201 | // Compute the Chebyshev coefficients. | |
202 | // Get the values of the shifted Chebyshev polynomials. | |
203 | std::vector<T> chebyshev_t_shifted_values; | |
204 | const T z_shifted = ((Z / W) * static_cast<std::int32_t>(2)) - static_cast<std::int32_t>(1); | |
205 | ||
206 | chebyshev_t(static_cast<std::uint32_t>(C.size()), | |
207 | z_shifted, | |
208 | &chebyshev_t_shifted_values); | |
209 | ||
210 | // Luke: C ---------- COMPUTE SCALE FACTOR ---------- | |
211 | // Luke: C | |
212 | // Luke: C ---------- SCALE THE COEFFICIENTS ---------- | |
213 | // Luke: C | |
214 | ||
215 | // The coefficient scaling is preformed after the Chebyshev summation, | |
216 | // and it is carried out with a single division operation. | |
217 | bool b_neg = false; | |
218 | ||
219 | const T scale = std::accumulate(C.begin(), | |
220 | C.end(), | |
221 | T(0), | |
222 | [&b_neg](T& scale_sum, const T& ck) -> T | |
223 | { | |
224 | ((!b_neg) ? (scale_sum += ck) : (scale_sum -= ck)); | |
225 | b_neg = (!b_neg); | |
226 | return scale_sum; | |
227 | }); | |
228 | ||
229 | // Compute the result of the series expansion using unscaled coefficients. | |
230 | const T sum = std::inner_product(C.begin(), | |
231 | C.end(), | |
232 | chebyshev_t_shifted_values.begin(), | |
233 | T(0)); | |
234 | ||
235 | // Return the properly scaled result. | |
236 | return sum / scale; | |
237 | } | |
238 | ||
239 | protected: | |
240 | const T Z; | |
241 | const T W; | |
242 | mutable std::deque<T> C; | |
243 | ||
244 | hypergeometric_pfq_base(const T& z, | |
245 | const T& w) : Z(z), | |
246 | W(w), | |
247 | C(0u) { } | |
248 | ||
249 | virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale<T>() * 500.0F); } | |
250 | }; | |
251 | ||
252 | template<typename T> class ccoef4_hypergeometric_0f1 : public hypergeometric_pfq_base<T> | |
253 | { | |
254 | public: | |
255 | ccoef4_hypergeometric_0f1(const T& c, | |
256 | const T& z, | |
257 | const T& w) : hypergeometric_pfq_base<T>(z, w), | |
258 | CP(c) { } | |
259 | ||
260 | virtual ~ccoef4_hypergeometric_0f1() { } | |
261 | ||
262 | virtual void ccoef() const | |
263 | { | |
264 | // See Luke 1977 page 80. | |
265 | const std::int32_t N1 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1)); | |
266 | const std::int32_t N2 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(2)); | |
267 | ||
268 | // Luke: C ---------- START COMPUTING COEFFICIENTS USING ---------- | |
269 | // Luke: C ---------- BACKWARD RECURRENCE SCHEME ---------- | |
270 | // Luke: C | |
271 | T A3(0); | |
272 | T A2(0); | |
273 | T A1(boost::math::tools::root_epsilon<T>()); | |
274 | ||
275 | hypergeometric_pfq_base<T>::C.resize(1u, A1); | |
276 | ||
277 | std::int32_t X1 = N2; | |
278 | ||
279 | T C1 = T(1) - CP; | |
280 | ||
281 | const T Z1 = T(4) / hypergeometric_pfq_base<T>::W; | |
282 | ||
283 | for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k) | |
284 | { | |
285 | const T DIVFAC = T(1) / X1; | |
286 | ||
287 | --X1; | |
288 | ||
289 | // The terms have been slightly re-arranged resulting in lower complexity. | |
290 | // Parentheses have been added to avoid reliance on operator precedence. | |
291 | const T term = (A2 - ((A3 * DIVFAC) * X1)) | |
292 | + ((A2 * X1) * ((1 + (C1 + X1)) * Z1)) | |
293 | + ((A1 * X1) * ((DIVFAC - (C1 * Z1)) + (X1 * Z1))); | |
294 | ||
295 | hypergeometric_pfq_base<T>::C.push_front(term); | |
296 | ||
297 | A3 = A2; | |
298 | A2 = A1; | |
299 | A1 = hypergeometric_pfq_base<T>::C.front(); | |
300 | } | |
301 | ||
302 | hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2); | |
303 | } | |
304 | ||
305 | private: | |
306 | const T CP; | |
307 | }; | |
308 | ||
309 | template<typename T> class ccoef1_hypergeometric_1f0 : public hypergeometric_pfq_base<T> | |
310 | { | |
311 | public: | |
312 | ccoef1_hypergeometric_1f0(const T& a, | |
313 | const T& z, | |
314 | const T& w) : hypergeometric_pfq_base<T>(z, w), | |
315 | AP(a) { } | |
316 | ||
317 | virtual ~ccoef1_hypergeometric_1f0() { } | |
318 | ||
319 | virtual void ccoef() const | |
320 | { | |
321 | // See Luke 1977 page 67. | |
322 | const std::int32_t N1 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1)); | |
323 | const std::int32_t N2 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2)); | |
324 | ||
325 | // Luke: C ---------- START COMPUTING COEFFICIENTS USING ---------- | |
326 | // Luke: C ---------- BACKWARD RECURRENCE SCHEME ---------- | |
327 | // Luke: C | |
328 | T A2(0); | |
329 | T A1(boost::math::tools::root_epsilon<T>()); | |
330 | ||
331 | hypergeometric_pfq_base<T>::C.resize(1u, A1); | |
332 | ||
333 | std::int32_t X1 = N2; | |
334 | ||
335 | T V1 = T(1) - AP; | |
336 | ||
337 | // Here, we have corrected what appears to be an error in Luke's code. | |
338 | ||
339 | // Luke's original code listing has: | |
340 | // AFAC = 2 + FOUR/W | |
341 | // But it appears as though the correct form is: | |
342 | // AFAC = 2 - FOUR/W. | |
343 | ||
344 | const T AFAC = 2 - (T(4) / hypergeometric_pfq_base<T>::W); | |
345 | ||
346 | for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k) | |
347 | { | |
348 | --X1; | |
349 | ||
350 | // The terms have been slightly re-arranged resulting in lower complexity. | |
351 | // Parentheses have been added to avoid reliance on operator precedence. | |
352 | const T term = -(((X1 * AFAC) * A1) + ((X1 + V1) * A2)) / (X1 - V1); | |
353 | ||
354 | hypergeometric_pfq_base<T>::C.push_front(term); | |
355 | ||
356 | A2 = A1; | |
357 | A1 = hypergeometric_pfq_base<T>::C.front(); | |
358 | } | |
359 | ||
360 | hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2); | |
361 | } | |
362 | ||
363 | private: | |
364 | const T AP; | |
365 | ||
366 | virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale<T>() * 1600.0F); } | |
367 | }; | |
368 | ||
369 | template<typename T> class ccoef3_hypergeometric_1f1 : public hypergeometric_pfq_base<T> | |
370 | { | |
371 | public: | |
372 | ccoef3_hypergeometric_1f1(const T& a, | |
373 | const T& c, | |
374 | const T& z, | |
375 | const T& w) : hypergeometric_pfq_base<T>(z, w), | |
376 | AP(a), | |
377 | CP(c) { } | |
378 | ||
379 | virtual ~ccoef3_hypergeometric_1f1() { } | |
380 | ||
381 | virtual void ccoef() const | |
382 | { | |
383 | // See Luke 1977 page 74. | |
384 | const std::int32_t N1 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1)); | |
385 | const std::int32_t N2 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(2)); | |
386 | ||
387 | // Luke: C ---------- START COMPUTING COEFFICIENTS USING ---------- | |
388 | // Luke: C ---------- BACKWARD RECURRENCE SCHEME ---------- | |
389 | // Luke: C | |
390 | T A3(0); | |
391 | T A2(0); | |
392 | T A1(boost::math::tools::root_epsilon<T>()); | |
393 | ||
394 | hypergeometric_pfq_base<T>::C.resize(1u, A1); | |
395 | ||
396 | std::int32_t X = N1; | |
397 | std::int32_t X1 = N2; | |
398 | ||
399 | T XA = X + AP; | |
400 | T X3A = (X + 3) - AP; | |
401 | ||
402 | const T Z1 = T(4) / hypergeometric_pfq_base<T>::W; | |
403 | ||
404 | for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k) | |
405 | { | |
406 | --X; | |
407 | --X1; | |
408 | --XA; | |
409 | --X3A; | |
410 | ||
411 | const T X3A_over_X2 = X3A / static_cast<std::int32_t>(X + 2); | |
412 | ||
413 | // The terms have been slightly re-arranged resulting in lower complexity. | |
414 | // Parentheses have been added to avoid reliance on operator precedence. | |
415 | const T PART1 = A1 * (((X + CP) * Z1) - X3A_over_X2); | |
416 | const T PART2 = A2 * (Z1 * ((X + 3) - CP) + (XA / X1)); | |
417 | const T PART3 = A3 * X3A_over_X2; | |
418 | ||
419 | const T term = (((PART1 + PART2) + PART3) * X1) / XA; | |
420 | ||
421 | hypergeometric_pfq_base<T>::C.push_front(term); | |
422 | ||
423 | A3 = A2; | |
424 | A2 = A1; | |
425 | A1 = hypergeometric_pfq_base<T>::C.front(); | |
426 | } | |
427 | ||
428 | hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2); | |
429 | } | |
430 | ||
431 | private: | |
432 | const T AP; | |
433 | const T CP; | |
434 | }; | |
435 | ||
436 | template<typename T> class ccoef6_hypergeometric_1f2 : public hypergeometric_pfq_base<T> | |
437 | { | |
438 | public: | |
439 | ccoef6_hypergeometric_1f2(const T& a, | |
440 | const T& b, | |
441 | const T& c, | |
442 | const T& z, | |
443 | const T& w) : hypergeometric_pfq_base<T>(z, w), | |
444 | AP(a), | |
445 | BP(b), | |
446 | CP(c) { } | |
447 | ||
448 | virtual ~ccoef6_hypergeometric_1f2() { } | |
449 | ||
450 | virtual void ccoef() const | |
451 | { | |
452 | // See Luke 1977 page 85. | |
453 | const std::int32_t N1 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1)); | |
454 | ||
455 | // Luke: C ---------- START COMPUTING COEFFICIENTS USING ---------- | |
456 | // Luke: C ---------- BACKWARD RECURRENCE SCHEME ---------- | |
457 | // Luke: C | |
458 | T A4(0); | |
459 | T A3(0); | |
460 | T A2(0); | |
461 | T A1(boost::math::tools::root_epsilon<T>()); | |
462 | ||
463 | hypergeometric_pfq_base<T>::C.resize(1u, A1); | |
464 | ||
465 | std::int32_t X = N1; | |
466 | T PP = X + AP; | |
467 | ||
468 | const T Z1 = T(4) / hypergeometric_pfq_base<T>::W; | |
469 | ||
470 | for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k) | |
471 | { | |
472 | --X; | |
473 | --PP; | |
474 | ||
475 | const std::int32_t TWO_X = static_cast<std::int32_t>(X * 2); | |
476 | const std::int32_t X_PLUS_1 = static_cast<std::int32_t>(X + 1); | |
477 | const std::int32_t X_PLUS_3 = static_cast<std::int32_t>(X + 3); | |
478 | const std::int32_t X_PLUS_4 = static_cast<std::int32_t>(X + 4); | |
479 | ||
480 | const T QQ = T(TWO_X + 3) / static_cast<std::int32_t>(TWO_X + static_cast<std::int32_t>(5)); | |
481 | const T SS = (X + BP) * (X + CP); | |
482 | ||
483 | // The terms have been slightly re-arranged resulting in lower complexity. | |
484 | // Parentheses have been added to avoid reliance on operator precedence. | |
485 | const T PART1 = A1 * (((PP - (QQ * (PP + 1))) * 2) + (SS * Z1)); | |
486 | const T PART2 = (A2 * (X + 2)) * ((((TWO_X + 1) * PP) / X_PLUS_1) - ((QQ * 4) * (PP + 1)) + (((TWO_X + 3) * (PP + 2)) / X_PLUS_3) + ((Z1 * 2) * (SS - (QQ * (X_PLUS_1 + BP)) * (X_PLUS_1 + CP)))); | |
487 | const T PART3 = A3 * ((((X_PLUS_3 - AP) - (QQ * (X_PLUS_4 - AP))) * 2) + (((QQ * Z1) * (X_PLUS_4 - BP)) * (X_PLUS_4 - CP))); | |
488 | const T PART4 = ((A4 * QQ) * (X_PLUS_4 - AP)) / X_PLUS_3; | |
489 | ||
490 | const T term = (((PART1 - PART2) + (PART3 - PART4)) * X_PLUS_1) / PP; | |
491 | ||
492 | hypergeometric_pfq_base<T>::C.push_front(term); | |
493 | ||
494 | A4 = A3; | |
495 | A3 = A2; | |
496 | A2 = A1; | |
497 | A1 = hypergeometric_pfq_base<T>::C.front(); | |
498 | } | |
499 | ||
500 | hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2); | |
501 | } | |
502 | ||
503 | private: | |
504 | const T AP; | |
505 | const T BP; | |
506 | const T CP; | |
507 | }; | |
508 | ||
509 | template<typename T> class ccoef2_hypergeometric_2f1 : public hypergeometric_pfq_base<T> | |
510 | { | |
511 | public: | |
512 | ccoef2_hypergeometric_2f1(const T& a, | |
513 | const T& b, | |
514 | const T& c, | |
515 | const T& z, | |
516 | const T& w) : hypergeometric_pfq_base<T>(z, w), | |
517 | AP(a), | |
518 | BP(b), | |
519 | CP(c) { } | |
520 | ||
521 | virtual ~ccoef2_hypergeometric_2f1() { } | |
522 | ||
523 | virtual void ccoef() const | |
524 | { | |
525 | // See Luke 1977 page 59. | |
526 | const std::int32_t N1 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1)); | |
527 | const std::int32_t N2 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2)); | |
528 | ||
529 | // Luke: C ---------- START COMPUTING COEFFICIENTS USING ---------- | |
530 | // Luke: C ---------- BACKWARD RECURRENCE SCHEME ---------- | |
531 | // Luke: C | |
532 | T A3(0); | |
533 | T A2(0); | |
534 | T A1(boost::math::tools::root_epsilon<T>()); | |
535 | ||
536 | hypergeometric_pfq_base<T>::C.resize(1u, A1); | |
537 | ||
538 | std::int32_t X = N1; | |
539 | std::int32_t X1 = N2; | |
540 | std::int32_t X3 = static_cast<std::int32_t>((X * 2) + 3); | |
541 | ||
542 | T X3A = (X + 3) - AP; | |
543 | T X3B = (X + 3) - BP; | |
544 | ||
545 | const T Z1 = T(4) / hypergeometric_pfq_base<T>::W; | |
546 | ||
547 | for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k) | |
548 | { | |
549 | --X; | |
550 | --X1; | |
551 | --X3A; | |
552 | --X3B; | |
553 | X3 -= 2; | |
554 | ||
555 | const std::int32_t X_PLUS_2 = static_cast<std::int32_t>(X + 2); | |
556 | ||
557 | const T XAB = T(1) / ((X + AP) * (X + BP)); | |
558 | ||
559 | // The terms have been slightly re-arranged resulting in lower complexity. | |
560 | // Parentheses have been added to avoid reliance on operator precedence. | |
561 | const T PART1 = (A1 * X1) * (2 - (((AP + X1) * (BP + X1)) * ((T(X3) / X_PLUS_2) * XAB)) + ((CP + X) * (XAB * Z1))); | |
562 | const T PART2 = (A2 * XAB) * ((X3A * X3B) - (X3 * ((X3A + X3B) - 1)) + (((3 - CP) + X) * (X1 * Z1))); | |
563 | const T PART3 = (A3 * X1) * (X3A / X_PLUS_2) * (X3B * XAB); | |
564 | ||
565 | const T term = (PART1 + PART2) - PART3; | |
566 | ||
567 | hypergeometric_pfq_base<T>::C.push_front(term); | |
568 | ||
569 | A3 = A2; | |
570 | A2 = A1; | |
571 | A1 = hypergeometric_pfq_base<T>::C.front(); | |
572 | } | |
573 | ||
574 | hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2); | |
575 | } | |
576 | ||
577 | private: | |
578 | const T AP; | |
579 | const T BP; | |
580 | const T CP; | |
581 | ||
582 | virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale<T>() * 1600.0F); } | |
583 | }; | |
584 | ||
585 | template<class T> T luke_ccoef4_hypergeometric_0f1(const T& a, const T& x); | |
586 | template<class T> T luke_ccoef1_hypergeometric_1f0(const T& a, const T& x); | |
587 | template<class T> T luke_ccoef3_hypergeometric_1f1(const T& a, const T& b, const T& x); | |
588 | template<class T> T luke_ccoef6_hypergeometric_1f2(const T& a, const T& b, const T& c, const T& x); | |
589 | template<class T> T luke_ccoef2_hypergeometric_2f1(const T& a, const T& b, const T& c, const T& x); | |
590 | } | |
591 | } | |
592 | ||
593 | template<class T> | |
594 | T examples::nr_006::luke_ccoef4_hypergeometric_0f1(const T& a, const T& x) | |
595 | { | |
596 | const ccoef4_hypergeometric_0f1<T> hypergeometric_0f1_object(a, x, T(-20)); | |
597 | ||
598 | hypergeometric_0f1_object.ccoef(); | |
599 | ||
600 | return hypergeometric_0f1_object.series(); | |
601 | } | |
602 | ||
603 | template<class T> | |
604 | T examples::nr_006::luke_ccoef1_hypergeometric_1f0(const T& a, const T& x) | |
605 | { | |
606 | const ccoef1_hypergeometric_1f0<T> hypergeometric_1f0_object(a, x, T(-20)); | |
607 | ||
608 | hypergeometric_1f0_object.ccoef(); | |
609 | ||
610 | return hypergeometric_1f0_object.series(); | |
611 | } | |
612 | ||
613 | template<class T> | |
614 | T examples::nr_006::luke_ccoef3_hypergeometric_1f1(const T& a, const T& b, const T& x) | |
615 | { | |
616 | const ccoef3_hypergeometric_1f1<T> hypergeometric_1f1_object(a, b, x, T(-20)); | |
617 | ||
618 | hypergeometric_1f1_object.ccoef(); | |
619 | ||
620 | return hypergeometric_1f1_object.series(); | |
621 | } | |
622 | ||
623 | template<class T> | |
624 | T examples::nr_006::luke_ccoef6_hypergeometric_1f2(const T& a, const T& b, const T& c, const T& x) | |
625 | { | |
626 | const ccoef6_hypergeometric_1f2<T> hypergeometric_1f2_object(a, b, c, x, T(-20)); | |
627 | ||
628 | hypergeometric_1f2_object.ccoef(); | |
629 | ||
630 | return hypergeometric_1f2_object.series(); | |
631 | } | |
632 | ||
633 | template<class T> | |
634 | T examples::nr_006::luke_ccoef2_hypergeometric_2f1(const T& a, const T& b, const T& c, const T& x) | |
635 | { | |
636 | const ccoef2_hypergeometric_2f1<T> hypergeometric_2f1_object(a, b, c, x, T(-20)); | |
637 | ||
638 | hypergeometric_2f1_object.ccoef(); | |
639 | ||
640 | return hypergeometric_2f1_object.series(); | |
641 | } | |
642 | ||
643 | int main() | |
644 | { | |
645 | stopwatch<STD_CHRONO::high_resolution_clock> my_stopwatch; | |
646 | float total_time = 0.0F; | |
647 | ||
648 | std::vector<mp_type> hypergeometric_0f1_results(20U); | |
649 | std::vector<mp_type> hypergeometric_1f0_results(20U); | |
650 | std::vector<mp_type> hypergeometric_1f1_results(20U); | |
651 | std::vector<mp_type> hypergeometric_2f1_results(20U); | |
652 | std::vector<mp_type> hypergeometric_1f2_results(20U); | |
653 | ||
654 | const mp_type a(mp_type(3) / 7); | |
655 | const mp_type b(mp_type(2) / 3); | |
656 | const mp_type c(mp_type(1) / 4); | |
657 | ||
658 | std::int_least16_t i; | |
659 | ||
660 | std::cout << "test hypergeometric_0f1." << std::endl; | |
661 | i = 1U; | |
662 | my_stopwatch.reset(); | |
663 | ||
664 | // Generate a table of values of Hypergeometric0F1. | |
665 | // Compare with the Mathematica command: | |
666 | // Table[N[HypergeometricPFQ[{}, {3/7}, -(i*EulerGamma)], 100], {i, 1, 20, 1}] | |
667 | std::for_each(hypergeometric_0f1_results.begin(), | |
668 | hypergeometric_0f1_results.end(), | |
669 | [&i, &a](mp_type& new_value) | |
670 | { | |
671 | const mp_type x(-(boost::math::constants::euler<mp_type>() * i)); | |
672 | ||
673 | new_value = examples::nr_006::luke_ccoef4_hypergeometric_0f1(a, x); | |
674 | ||
675 | ++i; | |
676 | }); | |
677 | ||
678 | total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count(); | |
679 | ||
680 | // Print the values of Hypergeometric0F1. | |
681 | std::for_each(hypergeometric_0f1_results.begin(), | |
682 | hypergeometric_0f1_results.end(), | |
683 | [](const mp_type& h) | |
684 | { | |
685 | std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl; | |
686 | }); | |
687 | ||
688 | std::cout << "test hypergeometric_1f0." << std::endl; | |
689 | i = 1U; | |
690 | my_stopwatch.reset(); | |
691 | ||
692 | // Generate a table of values of Hypergeometric1F0. | |
693 | // Compare with the Mathematica command: | |
694 | // Table[N[HypergeometricPFQ[{3/7}, {}, -(i*EulerGamma)], 100], {i, 1, 20, 1}] | |
695 | std::for_each(hypergeometric_1f0_results.begin(), | |
696 | hypergeometric_1f0_results.end(), | |
697 | [&i, &a](mp_type& new_value) | |
698 | { | |
699 | const mp_type x(-(boost::math::constants::euler<mp_type>() * i)); | |
700 | ||
701 | new_value = examples::nr_006::luke_ccoef1_hypergeometric_1f0(a, x); | |
702 | ||
703 | ++i; | |
704 | }); | |
705 | ||
706 | total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count(); | |
707 | ||
708 | // Print the values of Hypergeometric1F0. | |
709 | std::for_each(hypergeometric_1f0_results.begin(), | |
710 | hypergeometric_1f0_results.end(), | |
711 | [](const mp_type& h) | |
712 | { | |
713 | std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl; | |
714 | }); | |
715 | ||
716 | std::cout << "test hypergeometric_1f1." << std::endl; | |
717 | i = 1U; | |
718 | my_stopwatch.reset(); | |
719 | ||
720 | // Generate a table of values of Hypergeometric1F1. | |
721 | // Compare with the Mathematica command: | |
722 | // Table[N[HypergeometricPFQ[{3/7}, {2/3}, -(i*EulerGamma)], 100], {i, 1, 20, 1}] | |
723 | std::for_each(hypergeometric_1f1_results.begin(), | |
724 | hypergeometric_1f1_results.end(), | |
725 | [&i, &a, &b](mp_type& new_value) | |
726 | { | |
727 | const mp_type x(-(boost::math::constants::euler<mp_type>() * i)); | |
728 | ||
729 | new_value = examples::nr_006::luke_ccoef3_hypergeometric_1f1(a, b, x); | |
730 | ||
731 | ++i; | |
732 | }); | |
733 | ||
734 | total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count(); | |
735 | ||
736 | // Print the values of Hypergeometric1F1. | |
737 | std::for_each(hypergeometric_1f1_results.begin(), | |
738 | hypergeometric_1f1_results.end(), | |
739 | [](const mp_type& h) | |
740 | { | |
741 | std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl; | |
742 | }); | |
743 | ||
744 | std::cout << "test hypergeometric_1f2." << std::endl; | |
745 | i = 1U; | |
746 | my_stopwatch.reset(); | |
747 | ||
748 | // Generate a table of values of Hypergeometric1F2. | |
749 | // Compare with the Mathematica command: | |
750 | // Table[N[HypergeometricPFQ[{3/7}, {2/3, 1/4}, -(i*EulerGamma)], 100], {i, 1, 20, 1}] | |
751 | std::for_each(hypergeometric_1f2_results.begin(), | |
752 | hypergeometric_1f2_results.end(), | |
753 | [&i, &a, &b, &c](mp_type& new_value) | |
754 | { | |
755 | const mp_type x(-(boost::math::constants::euler<mp_type>() * i)); | |
756 | ||
757 | new_value = examples::nr_006::luke_ccoef6_hypergeometric_1f2(a, b, c, x); | |
758 | ||
759 | ++i; | |
760 | }); | |
761 | ||
762 | total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count(); | |
763 | ||
764 | // Print the values of Hypergeometric1F2. | |
765 | std::for_each(hypergeometric_1f2_results.begin(), | |
766 | hypergeometric_1f2_results.end(), | |
767 | [](const mp_type& h) | |
768 | { | |
769 | std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl; | |
770 | }); | |
771 | ||
772 | std::cout << "test hypergeometric_2f1." << std::endl; | |
773 | i = 1U; | |
774 | my_stopwatch.reset(); | |
775 | ||
776 | // Generate a table of values of Hypergeometric2F1. | |
777 | // Compare with the Mathematica command: | |
778 | // Table[N[HypergeometricPFQ[{3/7, 2/3}, {1/4}, -(i * EulerGamma)], 100], {i, 1, 20, 1}] | |
779 | std::for_each(hypergeometric_2f1_results.begin(), | |
780 | hypergeometric_2f1_results.end(), | |
781 | [&i, &a, &b, &c](mp_type& new_value) | |
782 | { | |
783 | const mp_type x(-(boost::math::constants::euler<mp_type>() * i)); | |
784 | ||
785 | new_value = examples::nr_006::luke_ccoef2_hypergeometric_2f1(a, b, c, x); | |
786 | ||
787 | ++i; | |
788 | }); | |
789 | ||
790 | total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count(); | |
791 | ||
792 | // Print the values of Hypergeometric2F1. | |
793 | std::for_each(hypergeometric_2f1_results.begin(), | |
794 | hypergeometric_2f1_results.end(), | |
795 | [](const mp_type& h) | |
796 | { | |
797 | std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl; | |
798 | }); | |
799 | ||
800 | std::cout << "Total execution time = " << std::setprecision(3) << total_time << "s" << std::endl; | |
801 | } |