]>
Commit | Line | Data |
---|---|---|
7c62b943 BK |
1 | // Special functions -*- C++ -*- |
2 | ||
748086b7 | 3 | // Copyright (C) 2006, 2007, 2008, 2009 |
7c62b943 BK |
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 | |
748086b7 | 9 | // Free Software Foundation; either version 3, or (at your option) |
7c62b943 BK |
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 | // | |
748086b7 JJ |
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/>. | |
7c62b943 BK |
25 | |
26 | /** @file tr1/exp_integral.tcc | |
27 | * This is an internal header file, included by other library headers. | |
28 | * You should not attempt to use it directly. | |
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 | ||
e133ace8 PC |
46 | #ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC |
47 | #define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1 | |
7c62b943 BK |
48 | |
49 | #include "special_function_util.h" | |
50 | ||
51 | namespace std | |
52 | { | |
e133ace8 PC |
53 | namespace tr1 |
54 | { | |
7c62b943 BK |
55 | |
56 | // [5.2] Special functions | |
57 | ||
7c62b943 | 58 | // Implementation-space details. |
7c62b943 BK |
59 | namespace __detail |
60 | { | |
61 | ||
62 | /** | |
63 | * @brief Return the exponential integral @f$ E_1(x) @f$ | |
64 | * by series summation. This should be good | |
65 | * for @f$ x < 1 @f$. | |
66 | * | |
67 | * The exponential integral is given by | |
68 | * \f[ | |
69 | * E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt | |
70 | * \f] | |
71 | * | |
72 | * @param __x The argument of the exponential integral function. | |
73 | * @return The exponential integral. | |
74 | */ | |
75 | template<typename _Tp> | |
76 | _Tp | |
77 | __expint_E1_series(const _Tp __x) | |
78 | { | |
79 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
80 | _Tp __term = _Tp(1); | |
81 | _Tp __esum = _Tp(0); | |
82 | _Tp __osum = _Tp(0); | |
83 | const unsigned int __max_iter = 100; | |
84 | for (unsigned int __i = 1; __i < __max_iter; ++__i) | |
85 | { | |
86 | __term *= - __x / __i; | |
87 | if (std::abs(__term) < __eps) | |
88 | break; | |
89 | if (__term >= _Tp(0)) | |
90 | __esum += __term / __i; | |
91 | else | |
92 | __osum += __term / __i; | |
93 | } | |
94 | ||
95 | return - __esum - __osum | |
96 | - __numeric_constants<_Tp>::__gamma_e() - std::log(__x); | |
97 | } | |
98 | ||
99 | ||
100 | /** | |
101 | * @brief Return the exponential integral @f$ E_1(x) @f$ | |
102 | * by asymptotic expansion. | |
103 | * | |
104 | * The exponential integral is given by | |
105 | * \f[ | |
106 | * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt | |
107 | * \f] | |
108 | * | |
109 | * @param __x The argument of the exponential integral function. | |
110 | * @return The exponential integral. | |
111 | */ | |
112 | template<typename _Tp> | |
113 | _Tp | |
114 | __expint_E1_asymp(const _Tp __x) | |
115 | { | |
116 | _Tp __term = _Tp(1); | |
117 | _Tp __esum = _Tp(1); | |
118 | _Tp __osum = _Tp(0); | |
119 | const unsigned int __max_iter = 1000; | |
120 | for (unsigned int __i = 1; __i < __max_iter; ++__i) | |
121 | { | |
122 | _Tp __prev = __term; | |
123 | __term *= - __i / __x; | |
124 | if (std::abs(__term) > std::abs(__prev)) | |
125 | break; | |
126 | if (__term >= _Tp(0)) | |
127 | __esum += __term; | |
128 | else | |
129 | __osum += __term; | |
130 | } | |
131 | ||
132 | return std::exp(- __x) * (__esum + __osum) / __x; | |
133 | } | |
134 | ||
135 | ||
136 | /** | |
137 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
138 | * by series summation. | |
139 | * | |
140 | * The exponential integral is given by | |
141 | * \f[ | |
142 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
143 | * \f] | |
144 | * | |
145 | * @param __n The order of the exponential integral function. | |
146 | * @param __x The argument of the exponential integral function. | |
147 | * @return The exponential integral. | |
148 | */ | |
149 | template<typename _Tp> | |
150 | _Tp | |
151 | __expint_En_series(const unsigned int __n, const _Tp __x) | |
152 | { | |
153 | const unsigned int __max_iter = 100; | |
154 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
155 | const int __nm1 = __n - 1; | |
156 | _Tp __ans = (__nm1 != 0 | |
157 | ? _Tp(1) / __nm1 : -std::log(__x) | |
158 | - __numeric_constants<_Tp>::__gamma_e()); | |
159 | _Tp __fact = _Tp(1); | |
160 | for (int __i = 1; __i <= __max_iter; ++__i) | |
161 | { | |
162 | __fact *= -__x / _Tp(__i); | |
163 | _Tp __del; | |
164 | if ( __i != __nm1 ) | |
165 | __del = -__fact / _Tp(__i - __nm1); | |
166 | else | |
167 | { | |
168 | _Tp __psi = -_TR1_GAMMA_TCC; | |
169 | for (int __ii = 1; __ii <= __nm1; ++__ii) | |
170 | __psi += _Tp(1) / _Tp(__ii); | |
171 | __del = __fact * (__psi - std::log(__x)); | |
172 | } | |
173 | __ans += __del; | |
174 | if (std::abs(__del) < __eps * std::abs(__ans)) | |
175 | return __ans; | |
176 | } | |
177 | std::__throw_runtime_error(__N("Series summation failed " | |
178 | "in __expint_En_series.")); | |
179 | } | |
180 | ||
181 | ||
182 | /** | |
183 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
184 | * by continued fractions. | |
185 | * | |
186 | * The exponential integral is given by | |
187 | * \f[ | |
188 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
189 | * \f] | |
190 | * | |
191 | * @param __n The order of the exponential integral function. | |
192 | * @param __x The argument of the exponential integral function. | |
193 | * @return The exponential integral. | |
194 | */ | |
195 | template<typename _Tp> | |
196 | _Tp | |
197 | __expint_En_cont_frac(const unsigned int __n, const _Tp __x) | |
198 | { | |
199 | const unsigned int __max_iter = 100; | |
200 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
201 | const _Tp __fp_min = std::numeric_limits<_Tp>::min(); | |
202 | const int __nm1 = __n - 1; | |
203 | _Tp __b = __x + _Tp(__n); | |
204 | _Tp __c = _Tp(1) / __fp_min; | |
205 | _Tp __d = _Tp(1) / __b; | |
206 | _Tp __h = __d; | |
207 | for ( unsigned int __i = 1; __i <= __max_iter; ++__i ) | |
208 | { | |
209 | _Tp __a = -_Tp(__i * (__nm1 + __i)); | |
210 | __b += _Tp(2); | |
211 | __d = _Tp(1) / (__a * __d + __b); | |
212 | __c = __b + __a / __c; | |
213 | const _Tp __del = __c * __d; | |
214 | __h *= __del; | |
215 | if (std::abs(__del - _Tp(1)) < __eps) | |
216 | { | |
217 | const _Tp __ans = __h * std::exp(-__x); | |
218 | return __ans; | |
219 | } | |
220 | } | |
221 | std::__throw_runtime_error(__N("Continued fraction failed " | |
222 | "in __expint_En_cont_frac.")); | |
223 | } | |
224 | ||
225 | ||
226 | /** | |
227 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
228 | * by recursion. Use upward recursion for @f$ x < n @f$ | |
229 | * and downward recursion (Miller's algorithm) otherwise. | |
230 | * | |
231 | * The exponential integral is given by | |
232 | * \f[ | |
233 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
234 | * \f] | |
235 | * | |
236 | * @param __n The order of the exponential integral function. | |
237 | * @param __x The argument of the exponential integral function. | |
238 | * @return The exponential integral. | |
239 | */ | |
240 | template<typename _Tp> | |
241 | _Tp | |
242 | __expint_En_recursion(const unsigned int __n, const _Tp __x) | |
243 | { | |
244 | _Tp __En; | |
245 | _Tp __E1 = __expint_E1(__x); | |
246 | if (__x < _Tp(__n)) | |
247 | { | |
248 | // Forward recursion is stable only for n < x. | |
249 | __En = __E1; | |
250 | for (unsigned int __j = 2; __j < __n; ++__j) | |
251 | __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1); | |
252 | } | |
253 | else | |
254 | { | |
255 | // Backward recursion is stable only for n >= x. | |
256 | __En = _Tp(1); | |
257 | const int __N = __n + 20; // TODO: Check this starting number. | |
258 | _Tp __save = _Tp(0); | |
259 | for (int __j = __N; __j > 0; --__j) | |
260 | { | |
261 | __En = (std::exp(-__x) - __j * __En) / __x; | |
262 | if (__j == __n) | |
263 | __save = __En; | |
264 | } | |
265 | _Tp __norm = __En / __E1; | |
266 | __En /= __norm; | |
267 | } | |
268 | ||
269 | return __En; | |
270 | } | |
271 | ||
272 | /** | |
273 | * @brief Return the exponential integral @f$ Ei(x) @f$ | |
274 | * by series summation. | |
275 | * | |
276 | * The exponential integral is given by | |
277 | * \f[ | |
278 | * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt | |
279 | * \f] | |
280 | * | |
281 | * @param __x The argument of the exponential integral function. | |
282 | * @return The exponential integral. | |
283 | */ | |
284 | template<typename _Tp> | |
285 | _Tp | |
286 | __expint_Ei_series(const _Tp __x) | |
287 | { | |
288 | _Tp __term = _Tp(1); | |
289 | _Tp __sum = _Tp(0); | |
290 | const unsigned int __max_iter = 1000; | |
291 | for (unsigned int __i = 1; __i < __max_iter; ++__i) | |
292 | { | |
293 | __term *= __x / __i; | |
294 | __sum += __term / __i; | |
295 | if (__term < std::numeric_limits<_Tp>::epsilon() * __sum) | |
296 | break; | |
297 | } | |
298 | ||
299 | return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x); | |
300 | } | |
301 | ||
302 | ||
303 | /** | |
304 | * @brief Return the exponential integral @f$ Ei(x) @f$ | |
305 | * by asymptotic expansion. | |
306 | * | |
307 | * The exponential integral is given by | |
308 | * \f[ | |
309 | * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt | |
310 | * \f] | |
311 | * | |
312 | * @param __x The argument of the exponential integral function. | |
313 | * @return The exponential integral. | |
314 | */ | |
315 | template<typename _Tp> | |
316 | _Tp | |
317 | __expint_Ei_asymp(const _Tp __x) | |
318 | { | |
319 | _Tp __term = _Tp(1); | |
320 | _Tp __sum = _Tp(1); | |
321 | const unsigned int __max_iter = 1000; | |
322 | for (unsigned int __i = 1; __i < __max_iter; ++__i) | |
323 | { | |
324 | _Tp __prev = __term; | |
325 | __term *= __i / __x; | |
326 | if (__term < std::numeric_limits<_Tp>::epsilon()) | |
327 | break; | |
328 | if (__term >= __prev) | |
329 | break; | |
330 | __sum += __term; | |
331 | } | |
332 | ||
333 | return std::exp(__x) * __sum / __x; | |
334 | } | |
335 | ||
336 | ||
337 | /** | |
338 | * @brief Return the exponential integral @f$ Ei(x) @f$. | |
339 | * | |
340 | * The exponential integral is given by | |
341 | * \f[ | |
342 | * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt | |
343 | * \f] | |
344 | * | |
345 | * @param __x The argument of the exponential integral function. | |
346 | * @return The exponential integral. | |
347 | */ | |
348 | template<typename _Tp> | |
349 | _Tp | |
350 | __expint_Ei(const _Tp __x) | |
351 | { | |
352 | if (__x < _Tp(0)) | |
353 | return -__expint_E1(-__x); | |
354 | else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon())) | |
355 | return __expint_Ei_series(__x); | |
356 | else | |
357 | return __expint_Ei_asymp(__x); | |
358 | } | |
359 | ||
360 | ||
361 | /** | |
362 | * @brief Return the exponential integral @f$ E_1(x) @f$. | |
363 | * | |
364 | * The exponential integral is given by | |
365 | * \f[ | |
366 | * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt | |
367 | * \f] | |
368 | * | |
369 | * @param __x The argument of the exponential integral function. | |
370 | * @return The exponential integral. | |
371 | */ | |
372 | template<typename _Tp> | |
373 | _Tp | |
374 | __expint_E1(const _Tp __x) | |
375 | { | |
376 | if (__x < _Tp(0)) | |
377 | return -__expint_Ei(-__x); | |
378 | else if (__x < _Tp(1)) | |
379 | return __expint_E1_series(__x); | |
380 | else if (__x < _Tp(100)) // TODO: Find a good asymptotic switch point. | |
381 | return __expint_En_cont_frac(1, __x); | |
382 | else | |
383 | return __expint_E1_asymp(__x); | |
384 | } | |
385 | ||
386 | ||
387 | /** | |
388 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
389 | * for large argument. | |
390 | * | |
391 | * The exponential integral is given by | |
392 | * \f[ | |
393 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
394 | * \f] | |
395 | * | |
396 | * This is something of an extension. | |
397 | * | |
398 | * @param __n The order of the exponential integral function. | |
399 | * @param __x The argument of the exponential integral function. | |
400 | * @return The exponential integral. | |
401 | */ | |
402 | template<typename _Tp> | |
403 | _Tp | |
404 | __expint_asymp(const unsigned int __n, const _Tp __x) | |
405 | { | |
406 | _Tp __term = _Tp(1); | |
407 | _Tp __sum = _Tp(1); | |
408 | for (unsigned int __i = 1; __i <= __n; ++__i) | |
409 | { | |
410 | _Tp __prev = __term; | |
411 | __term *= -(__n - __i + 1) / __x; | |
412 | if (std::abs(__term) > std::abs(__prev)) | |
413 | break; | |
414 | __sum += __term; | |
415 | } | |
416 | ||
417 | return std::exp(-__x) * __sum / __x; | |
418 | } | |
419 | ||
420 | ||
421 | /** | |
422 | * @brief Return the exponential integral @f$ E_n(x) @f$ | |
423 | * for large order. | |
424 | * | |
425 | * The exponential integral is given by | |
426 | * \f[ | |
427 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
428 | * \f] | |
429 | * | |
430 | * This is something of an extension. | |
431 | * | |
432 | * @param __n The order of the exponential integral function. | |
433 | * @param __x The argument of the exponential integral function. | |
434 | * @return The exponential integral. | |
435 | */ | |
436 | template<typename _Tp> | |
437 | _Tp | |
438 | __expint_large_n(const unsigned int __n, const _Tp __x) | |
439 | { | |
440 | const _Tp __xpn = __x + __n; | |
441 | const _Tp __xpn2 = __xpn * __xpn; | |
442 | _Tp __term = _Tp(1); | |
443 | _Tp __sum = _Tp(1); | |
444 | for (unsigned int __i = 1; __i <= __n; ++__i) | |
445 | { | |
446 | _Tp __prev = __term; | |
447 | __term *= (__n - 2 * (__i - 1) * __x) / __xpn2; | |
448 | if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon()) | |
449 | break; | |
450 | __sum += __term; | |
451 | } | |
452 | ||
453 | return std::exp(-__x) * __sum / __xpn; | |
454 | } | |
455 | ||
456 | ||
457 | /** | |
458 | * @brief Return the exponential integral @f$ E_n(x) @f$. | |
459 | * | |
460 | * The exponential integral is given by | |
461 | * \f[ | |
462 | * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt | |
463 | * \f] | |
464 | * This is something of an extension. | |
465 | * | |
466 | * @param __n The order of the exponential integral function. | |
467 | * @param __x The argument of the exponential integral function. | |
468 | * @return The exponential integral. | |
469 | */ | |
470 | template<typename _Tp> | |
471 | _Tp | |
472 | __expint(const unsigned int __n, const _Tp __x) | |
473 | { | |
474 | // Return NaN on NaN input. | |
475 | if (__isnan(__x)) | |
476 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
477 | else if (__n <= 1 && __x == _Tp(0)) | |
478 | return std::numeric_limits<_Tp>::infinity(); | |
479 | else | |
480 | { | |
481 | _Tp __E0 = std::exp(__x) / __x; | |
482 | if (__n == 0) | |
483 | return __E0; | |
484 | ||
485 | _Tp __E1 = __expint_E1(__x); | |
486 | if (__n == 1) | |
487 | return __E1; | |
488 | ||
489 | if (__x == _Tp(0)) | |
490 | return _Tp(1) / static_cast<_Tp>(__n - 1); | |
491 | ||
492 | _Tp __En = __expint_En_recursion(__n, __x); | |
493 | ||
494 | return __En; | |
495 | } | |
496 | } | |
497 | ||
498 | ||
499 | /** | |
6165bbdd | 500 | * @brief Return the exponential integral @f$ Ei(x) @f$. |
7c62b943 BK |
501 | * |
502 | * The exponential integral is given by | |
503 | * \f[ | |
504 | * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt | |
505 | * \f] | |
506 | * | |
507 | * @param __x The argument of the exponential integral function. | |
508 | * @return The exponential integral. | |
509 | */ | |
510 | template<typename _Tp> | |
511 | inline _Tp | |
512 | __expint(const _Tp __x) | |
513 | { | |
514 | if (__isnan(__x)) | |
515 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
516 | else | |
517 | return __expint_Ei(__x); | |
518 | } | |
519 | ||
520 | } // namespace std::tr1::__detail | |
e133ace8 | 521 | } |
7c62b943 BK |
522 | } |
523 | ||
e133ace8 | 524 | #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC |