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