]>
Commit | Line | Data |
---|---|---|
b34f60ac | 1 | // Special functions -*- C++ -*- |
2 | ||
f1717362 | 3 | // Copyright (C) 2006-2016 Free Software Foundation, Inc. |
b34f60ac | 4 | // |
5 | // This file is part of the GNU ISO C++ Library. This library is free | |
6 | // software; you can redistribute it and/or modify it under the | |
7 | // terms of the GNU General Public License as published by the | |
6bc9506f | 8 | // Free Software Foundation; either version 3, or (at your option) |
b34f60ac | 9 | // any later version. |
10 | // | |
11 | // This library is distributed in the hope that it will be useful, | |
12 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
14 | // GNU General Public License for more details. | |
15 | // | |
6bc9506f | 16 | // Under Section 7 of GPL version 3, you are granted additional |
17 | // permissions described in the GCC Runtime Library Exception, version | |
18 | // 3.1, as published by the Free Software Foundation. | |
19 | ||
20 | // You should have received a copy of the GNU General Public License and | |
21 | // a copy of the GCC Runtime Library Exception along with this program; | |
22 | // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see | |
23 | // <http://www.gnu.org/licenses/>. | |
b34f60ac | 24 | |
25 | /** @file tr1/exp_integral.tcc | |
26 | * This is an internal header file, included by other library headers. | |
5846aeac | 27 | * Do not attempt to use it directly. @headername{tr1/cmath} |
b34f60ac | 28 | */ |
29 | ||
30 | // | |
31 | // ISO C++ 14882 TR1: 5.2 Special functions | |
32 | // | |
33 | ||
34 | // Written by Edward Smith-Rowland based on: | |
35 | // | |
36 | // (1) Handbook of Mathematical Functions, | |
37 | // Ed. by Milton Abramowitz and Irene A. Stegun, | |
38 | // Dover Publications, New-York, Section 5, pp. 228-251. | |
39 | // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl | |
40 | // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, | |
41 | // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), | |
42 | // 2nd ed, pp. 222-225. | |
43 | // | |
44 | ||
c17b0a1c | 45 | #ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC |
46 | #define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1 | |
b34f60ac | 47 | |
48 | #include "special_function_util.h" | |
49 | ||
2948dd21 | 50 | namespace std _GLIBCXX_VISIBILITY(default) |
b34f60ac | 51 | { |
c17b0a1c | 52 | namespace tr1 |
53 | { | |
b34f60ac | 54 | // [5.2] Special functions |
55 | ||
b34f60ac | 56 | // Implementation-space details. |
b34f60ac | 57 | namespace __detail |
58 | { | |
2948dd21 | 59 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
b34f60ac | 60 | |
cd7f5f45 | 61 | template<typename _Tp> _Tp __expint_E1(_Tp); |
8411500a | 62 | |
b34f60ac | 63 | /** |
64 | * @brief Return the exponential integral @f$ E_1(x) @f$ | |
65 | * by series summation. This should be good | |
66 | * for @f$ x < 1 @f$. | |
67 | * | |
68 | * The exponential integral is given by | |
69 | * \f[ | |
70 | * E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt | |
71 | * \f] | |
72 | * | |
73 | * @param __x The argument of the exponential integral function. | |
74 | * @return The exponential integral. | |
75 | */ | |
76 | template<typename _Tp> | |
77 | _Tp | |
cd7f5f45 | 78 | __expint_E1_series(_Tp __x) |
b34f60ac | 79 | { |
80 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
81 | _Tp __term = _Tp(1); | |
82 | _Tp __esum = _Tp(0); | |
83 | _Tp __osum = _Tp(0); | |
84 | const unsigned int __max_iter = 100; | |
85 | for (unsigned int __i = 1; __i < __max_iter; ++__i) | |
86 | { | |
87 | __term *= - __x / __i; | |
88 | if (std::abs(__term) < __eps) | |
89 | break; | |
90 | if (__term >= _Tp(0)) | |
91 | __esum += __term / __i; | |
92 | else | |
93 | __osum += __term / __i; | |
94 | } | |
95 | ||
96 | return - __esum - __osum | |
97 | - __numeric_constants<_Tp>::__gamma_e() - std::log(__x); | |
98 | } | |
99 | ||
100 | ||
101 | /** | |
102 | * @brief Return the exponential integral @f$ E_1(x) @f$ | |
103 | * by asymptotic expansion. | |
104 | * | |
105 | * The exponential integral is given by | |
106 | * \f[ | |
107 | * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt | |
108 | * \f] | |
109 | * | |
110 | * @param __x The argument of the exponential integral function. | |
111 | * @return The exponential integral. | |
112 | */ | |
113 | template<typename _Tp> | |
114 | _Tp | |
cd7f5f45 | 115 | __expint_E1_asymp(_Tp __x) |
b34f60ac | 116 | { |
117 | _Tp __term = _Tp(1); | |
118 | _Tp __esum = _Tp(1); | |
119 | _Tp __osum = _Tp(0); | |
120 | const unsigned int __max_iter = 1000; | |
121 | for (unsigned int __i = 1; __i < __max_iter; ++__i) | |
122 | { | |
123 | _Tp __prev = __term; | |
124 | __term *= - __i / __x; | |
125 | if (std::abs(__term) > std::abs(__prev)) | |
126 | break; | |
127 | if (__term >= _Tp(0)) | |
128 | __esum += __term; | |
129 | else | |
130 | __osum += __term; | |
131 | } | |
132 | ||
133 | return std::exp(- __x) * (__esum + __osum) / __x; | |
134 | } | |
135 | ||
136 | ||
137 | /** | |
138 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
139 | * by series summation. | |
140 | * | |
141 | * The exponential integral is given by | |
142 | * \f[ | |
143 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
144 | * \f] | |
145 | * | |
146 | * @param __n The order of the exponential integral function. | |
147 | * @param __x The argument of the exponential integral function. | |
148 | * @return The exponential integral. | |
149 | */ | |
150 | template<typename _Tp> | |
151 | _Tp | |
cd7f5f45 | 152 | __expint_En_series(unsigned int __n, _Tp __x) |
b34f60ac | 153 | { |
154 | const unsigned int __max_iter = 100; | |
155 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
156 | const int __nm1 = __n - 1; | |
157 | _Tp __ans = (__nm1 != 0 | |
158 | ? _Tp(1) / __nm1 : -std::log(__x) | |
159 | - __numeric_constants<_Tp>::__gamma_e()); | |
160 | _Tp __fact = _Tp(1); | |
161 | for (int __i = 1; __i <= __max_iter; ++__i) | |
162 | { | |
163 | __fact *= -__x / _Tp(__i); | |
164 | _Tp __del; | |
165 | if ( __i != __nm1 ) | |
166 | __del = -__fact / _Tp(__i - __nm1); | |
167 | else | |
168 | { | |
cbbbc10c | 169 | _Tp __psi = -__numeric_constants<_Tp>::gamma_e(); |
b34f60ac | 170 | for (int __ii = 1; __ii <= __nm1; ++__ii) |
171 | __psi += _Tp(1) / _Tp(__ii); | |
172 | __del = __fact * (__psi - std::log(__x)); | |
173 | } | |
174 | __ans += __del; | |
175 | if (std::abs(__del) < __eps * std::abs(__ans)) | |
176 | return __ans; | |
177 | } | |
178 | std::__throw_runtime_error(__N("Series summation failed " | |
179 | "in __expint_En_series.")); | |
180 | } | |
181 | ||
182 | ||
183 | /** | |
184 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
185 | * by continued fractions. | |
186 | * | |
187 | * The exponential integral is given by | |
188 | * \f[ | |
189 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
190 | * \f] | |
191 | * | |
192 | * @param __n The order of the exponential integral function. | |
193 | * @param __x The argument of the exponential integral function. | |
194 | * @return The exponential integral. | |
195 | */ | |
196 | template<typename _Tp> | |
197 | _Tp | |
cd7f5f45 | 198 | __expint_En_cont_frac(unsigned int __n, _Tp __x) |
b34f60ac | 199 | { |
200 | const unsigned int __max_iter = 100; | |
201 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
202 | const _Tp __fp_min = std::numeric_limits<_Tp>::min(); | |
203 | const int __nm1 = __n - 1; | |
204 | _Tp __b = __x + _Tp(__n); | |
205 | _Tp __c = _Tp(1) / __fp_min; | |
206 | _Tp __d = _Tp(1) / __b; | |
207 | _Tp __h = __d; | |
208 | for ( unsigned int __i = 1; __i <= __max_iter; ++__i ) | |
209 | { | |
210 | _Tp __a = -_Tp(__i * (__nm1 + __i)); | |
211 | __b += _Tp(2); | |
212 | __d = _Tp(1) / (__a * __d + __b); | |
213 | __c = __b + __a / __c; | |
214 | const _Tp __del = __c * __d; | |
215 | __h *= __del; | |
216 | if (std::abs(__del - _Tp(1)) < __eps) | |
217 | { | |
218 | const _Tp __ans = __h * std::exp(-__x); | |
219 | return __ans; | |
220 | } | |
221 | } | |
222 | std::__throw_runtime_error(__N("Continued fraction failed " | |
223 | "in __expint_En_cont_frac.")); | |
224 | } | |
225 | ||
226 | ||
227 | /** | |
228 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
229 | * by recursion. Use upward recursion for @f$ x < n @f$ | |
230 | * and downward recursion (Miller's algorithm) otherwise. | |
231 | * | |
232 | * The exponential integral is given by | |
233 | * \f[ | |
234 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
235 | * \f] | |
236 | * | |
237 | * @param __n The order of the exponential integral function. | |
238 | * @param __x The argument of the exponential integral function. | |
239 | * @return The exponential integral. | |
240 | */ | |
241 | template<typename _Tp> | |
242 | _Tp | |
cd7f5f45 | 243 | __expint_En_recursion(unsigned int __n, _Tp __x) |
b34f60ac | 244 | { |
245 | _Tp __En; | |
246 | _Tp __E1 = __expint_E1(__x); | |
247 | if (__x < _Tp(__n)) | |
248 | { | |
249 | // Forward recursion is stable only for n < x. | |
250 | __En = __E1; | |
251 | for (unsigned int __j = 2; __j < __n; ++__j) | |
252 | __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1); | |
253 | } | |
254 | else | |
255 | { | |
256 | // Backward recursion is stable only for n >= x. | |
257 | __En = _Tp(1); | |
258 | const int __N = __n + 20; // TODO: Check this starting number. | |
259 | _Tp __save = _Tp(0); | |
260 | for (int __j = __N; __j > 0; --__j) | |
261 | { | |
262 | __En = (std::exp(-__x) - __j * __En) / __x; | |
263 | if (__j == __n) | |
264 | __save = __En; | |
265 | } | |
266 | _Tp __norm = __En / __E1; | |
267 | __En /= __norm; | |
268 | } | |
269 | ||
270 | return __En; | |
271 | } | |
272 | ||
273 | /** | |
274 | * @brief Return the exponential integral @f$ Ei(x) @f$ | |
275 | * by series summation. | |
276 | * | |
277 | * The exponential integral is given by | |
278 | * \f[ | |
279 | * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt | |
280 | * \f] | |
281 | * | |
282 | * @param __x The argument of the exponential integral function. | |
283 | * @return The exponential integral. | |
284 | */ | |
285 | template<typename _Tp> | |
286 | _Tp | |
cd7f5f45 | 287 | __expint_Ei_series(_Tp __x) |
b34f60ac | 288 | { |
289 | _Tp __term = _Tp(1); | |
290 | _Tp __sum = _Tp(0); | |
291 | const unsigned int __max_iter = 1000; | |
292 | for (unsigned int __i = 1; __i < __max_iter; ++__i) | |
293 | { | |
294 | __term *= __x / __i; | |
295 | __sum += __term / __i; | |
296 | if (__term < std::numeric_limits<_Tp>::epsilon() * __sum) | |
297 | break; | |
298 | } | |
299 | ||
300 | return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x); | |
301 | } | |
302 | ||
303 | ||
304 | /** | |
305 | * @brief Return the exponential integral @f$ Ei(x) @f$ | |
306 | * by asymptotic expansion. | |
307 | * | |
308 | * The exponential integral is given by | |
309 | * \f[ | |
310 | * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt | |
311 | * \f] | |
312 | * | |
313 | * @param __x The argument of the exponential integral function. | |
314 | * @return The exponential integral. | |
315 | */ | |
316 | template<typename _Tp> | |
317 | _Tp | |
cd7f5f45 | 318 | __expint_Ei_asymp(_Tp __x) |
b34f60ac | 319 | { |
320 | _Tp __term = _Tp(1); | |
321 | _Tp __sum = _Tp(1); | |
322 | const unsigned int __max_iter = 1000; | |
323 | for (unsigned int __i = 1; __i < __max_iter; ++__i) | |
324 | { | |
325 | _Tp __prev = __term; | |
326 | __term *= __i / __x; | |
327 | if (__term < std::numeric_limits<_Tp>::epsilon()) | |
328 | break; | |
329 | if (__term >= __prev) | |
330 | break; | |
331 | __sum += __term; | |
332 | } | |
333 | ||
334 | return std::exp(__x) * __sum / __x; | |
335 | } | |
336 | ||
337 | ||
338 | /** | |
339 | * @brief Return the exponential integral @f$ Ei(x) @f$. | |
340 | * | |
341 | * The exponential integral is given by | |
342 | * \f[ | |
343 | * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt | |
344 | * \f] | |
345 | * | |
346 | * @param __x The argument of the exponential integral function. | |
347 | * @return The exponential integral. | |
348 | */ | |
349 | template<typename _Tp> | |
350 | _Tp | |
cd7f5f45 | 351 | __expint_Ei(_Tp __x) |
b34f60ac | 352 | { |
353 | if (__x < _Tp(0)) | |
354 | return -__expint_E1(-__x); | |
355 | else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon())) | |
356 | return __expint_Ei_series(__x); | |
357 | else | |
358 | return __expint_Ei_asymp(__x); | |
359 | } | |
360 | ||
361 | ||
362 | /** | |
363 | * @brief Return the exponential integral @f$ E_1(x) @f$. | |
364 | * | |
365 | * The exponential integral is given by | |
366 | * \f[ | |
367 | * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt | |
368 | * \f] | |
369 | * | |
370 | * @param __x The argument of the exponential integral function. | |
371 | * @return The exponential integral. | |
372 | */ | |
373 | template<typename _Tp> | |
374 | _Tp | |
cd7f5f45 | 375 | __expint_E1(_Tp __x) |
b34f60ac | 376 | { |
377 | if (__x < _Tp(0)) | |
378 | return -__expint_Ei(-__x); | |
379 | else if (__x < _Tp(1)) | |
380 | return __expint_E1_series(__x); | |
381 | else if (__x < _Tp(100)) // TODO: Find a good asymptotic switch point. | |
382 | return __expint_En_cont_frac(1, __x); | |
383 | else | |
384 | return __expint_E1_asymp(__x); | |
385 | } | |
386 | ||
387 | ||
388 | /** | |
389 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
390 | * for large argument. | |
391 | * | |
392 | * The exponential integral is given by | |
393 | * \f[ | |
394 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
395 | * \f] | |
396 | * | |
397 | * This is something of an extension. | |
398 | * | |
399 | * @param __n The order of the exponential integral function. | |
400 | * @param __x The argument of the exponential integral function. | |
401 | * @return The exponential integral. | |
402 | */ | |
403 | template<typename _Tp> | |
404 | _Tp | |
cd7f5f45 | 405 | __expint_asymp(unsigned int __n, _Tp __x) |
b34f60ac | 406 | { |
407 | _Tp __term = _Tp(1); | |
408 | _Tp __sum = _Tp(1); | |
409 | for (unsigned int __i = 1; __i <= __n; ++__i) | |
410 | { | |
411 | _Tp __prev = __term; | |
412 | __term *= -(__n - __i + 1) / __x; | |
413 | if (std::abs(__term) > std::abs(__prev)) | |
414 | break; | |
415 | __sum += __term; | |
416 | } | |
417 | ||
418 | return std::exp(-__x) * __sum / __x; | |
419 | } | |
420 | ||
421 | ||
422 | /** | |
423 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
424 | * for large order. | |
425 | * | |
426 | * The exponential integral is given by | |
427 | * \f[ | |
428 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
429 | * \f] | |
430 | * | |
431 | * This is something of an extension. | |
432 | * | |
433 | * @param __n The order of the exponential integral function. | |
434 | * @param __x The argument of the exponential integral function. | |
435 | * @return The exponential integral. | |
436 | */ | |
437 | template<typename _Tp> | |
438 | _Tp | |
cd7f5f45 | 439 | __expint_large_n(unsigned int __n, _Tp __x) |
b34f60ac | 440 | { |
441 | const _Tp __xpn = __x + __n; | |
442 | const _Tp __xpn2 = __xpn * __xpn; | |
443 | _Tp __term = _Tp(1); | |
444 | _Tp __sum = _Tp(1); | |
445 | for (unsigned int __i = 1; __i <= __n; ++__i) | |
446 | { | |
447 | _Tp __prev = __term; | |
448 | __term *= (__n - 2 * (__i - 1) * __x) / __xpn2; | |
449 | if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon()) | |
450 | break; | |
451 | __sum += __term; | |
452 | } | |
453 | ||
454 | return std::exp(-__x) * __sum / __xpn; | |
455 | } | |
456 | ||
457 | ||
458 | /** | |
459 | * @brief Return the exponential integral @f$ E_n(x) @f$. | |
460 | * | |
461 | * The exponential integral is given by | |
462 | * \f[ | |
463 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
464 | * \f] | |
465 | * This is something of an extension. | |
466 | * | |
467 | * @param __n The order of the exponential integral function. | |
468 | * @param __x The argument of the exponential integral function. | |
469 | * @return The exponential integral. | |
470 | */ | |
471 | template<typename _Tp> | |
472 | _Tp | |
cd7f5f45 | 473 | __expint(unsigned int __n, _Tp __x) |
b34f60ac | 474 | { |
475 | // Return NaN on NaN input. | |
476 | if (__isnan(__x)) | |
477 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
478 | else if (__n <= 1 && __x == _Tp(0)) | |
479 | return std::numeric_limits<_Tp>::infinity(); | |
480 | else | |
481 | { | |
482 | _Tp __E0 = std::exp(__x) / __x; | |
483 | if (__n == 0) | |
484 | return __E0; | |
485 | ||
486 | _Tp __E1 = __expint_E1(__x); | |
487 | if (__n == 1) | |
488 | return __E1; | |
489 | ||
490 | if (__x == _Tp(0)) | |
491 | return _Tp(1) / static_cast<_Tp>(__n - 1); | |
492 | ||
493 | _Tp __En = __expint_En_recursion(__n, __x); | |
494 | ||
495 | return __En; | |
496 | } | |
497 | } | |
498 | ||
499 | ||
500 | /** | |
048ff85f | 501 | * @brief Return the exponential integral @f$ Ei(x) @f$. |
b34f60ac | 502 | * |
503 | * The exponential integral is given by | |
504 | * \f[ | |
505 | * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt | |
506 | * \f] | |
507 | * | |
508 | * @param __x The argument of the exponential integral function. | |
509 | * @return The exponential integral. | |
510 | */ | |
511 | template<typename _Tp> | |
512 | inline _Tp | |
cd7f5f45 | 513 | __expint(_Tp __x) |
b34f60ac | 514 | { |
515 | if (__isnan(__x)) | |
516 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
517 | else | |
518 | return __expint_Ei(__x); | |
519 | } | |
520 | ||
2948dd21 | 521 | _GLIBCXX_END_NAMESPACE_VERSION |
b34f60ac | 522 | } // namespace std::tr1::__detail |
c17b0a1c | 523 | } |
b34f60ac | 524 | } |
525 | ||
c17b0a1c | 526 | #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC |