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