1 /* specfunc/exp.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* Author: G. Jungman */
21
22 #include "gsl__config.h"
23 #include "gsl_math.h"
24 #include "gsl_errno.h"
25 #include "gsl_sf_gamma.h"
26 #include "gsl_sf_exp.h"
27
28 #include "gsl_specfunc__error.h"
29
30 /* Evaluate the continued fraction for exprel.
31 * [Abramowitz+Stegun, 4.2.41]
32 */
33 static
34 int
exprel_n_CF(const int N,const double x,gsl_sf_result * result)35 exprel_n_CF(const int N, const double x, gsl_sf_result * result)
36 {
37 const double RECUR_BIG = GSL_SQRT_DBL_MAX;
38 const int maxiter = 5000;
39 int n = 1;
40 double Anm2 = 1.0;
41 double Bnm2 = 0.0;
42 double Anm1 = 0.0;
43 double Bnm1 = 1.0;
44 double a1 = 1.0;
45 double b1 = 1.0;
46 double a2 = -x;
47 double b2 = N+1;
48 double an, bn;
49
50 double fn;
51
52 double An = b1*Anm1 + a1*Anm2; /* A1 */
53 double Bn = b1*Bnm1 + a1*Bnm2; /* B1 */
54
55 /* One explicit step, before we get to the main pattern. */
56 n++;
57 Anm2 = Anm1;
58 Bnm2 = Bnm1;
59 Anm1 = An;
60 Bnm1 = Bn;
61 An = b2*Anm1 + a2*Anm2; /* A2 */
62 Bn = b2*Bnm1 + a2*Bnm2; /* B2 */
63
64 fn = An/Bn;
65
66 while(n < maxiter) {
67 double old_fn;
68 double del;
69 n++;
70 Anm2 = Anm1;
71 Bnm2 = Bnm1;
72 Anm1 = An;
73 Bnm1 = Bn;
74 an = ( GSL_IS_ODD(n) ? ((n-1)/2)*x : -(N+(n/2)-1)*x );
75 bn = N + n - 1;
76 An = bn*Anm1 + an*Anm2;
77 Bn = bn*Bnm1 + an*Bnm2;
78
79 if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
80 An /= RECUR_BIG;
81 Bn /= RECUR_BIG;
82 Anm1 /= RECUR_BIG;
83 Bnm1 /= RECUR_BIG;
84 Anm2 /= RECUR_BIG;
85 Bnm2 /= RECUR_BIG;
86 }
87
88 old_fn = fn;
89 fn = An/Bn;
90 del = old_fn/fn;
91
92 if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;
93 }
94
95 result->val = fn;
96 result->err = 2.0*(n+1.0)*GSL_DBL_EPSILON*fabs(fn);
97
98 if(n == maxiter)
99 GSL_ERROR ("error", GSL_EMAXITER);
100 else
101 return GSL_SUCCESS;
102 }
103
104
105 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
106
gsl_sf_exp_e(const double x,gsl_sf_result * result)107 int gsl_sf_exp_e(const double x, gsl_sf_result * result)
108 {
109 if(x > GSL_LOG_DBL_MAX) {
110 OVERFLOW_ERROR(result);
111 }
112 else if(x < GSL_LOG_DBL_MIN) {
113 UNDERFLOW_ERROR(result);
114 }
115 else {
116 result->val = exp(x);
117 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
118 return GSL_SUCCESS;
119 }
120 }
121
gsl_sf_exp_e10_e(const double x,gsl_sf_result_e10 * result)122 int gsl_sf_exp_e10_e(const double x, gsl_sf_result_e10 * result)
123 {
124 if(x > INT_MAX-1) {
125 OVERFLOW_ERROR_E10(result);
126 }
127 else if(x < INT_MIN+1) {
128 UNDERFLOW_ERROR_E10(result);
129 }
130 else {
131 const int N = (int) floor(x/M_LN10);
132 result->val = exp(x-N*M_LN10);
133 result->err = 2.0 * (fabs(x)+1.0) * GSL_DBL_EPSILON * fabs(result->val);
134 result->e10 = N;
135 return GSL_SUCCESS;
136 }
137 }
138
139
gsl_sf_exp_mult_e(const double x,const double y,gsl_sf_result * result)140 int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result)
141 {
142 const double ay = fabs(y);
143
144 if(y == 0.0) {
145 result->val = 0.0;
146 result->err = 0.0;
147 return GSL_SUCCESS;
148 }
149 else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN)
150 && (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN)
151 ) {
152 const double ex = exp(x);
153 result->val = y * ex;
154 result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
155 return GSL_SUCCESS;
156 }
157 else {
158 const double ly = log(ay);
159 const double lnr = x + ly;
160
161 if(lnr > GSL_LOG_DBL_MAX - 0.01) {
162 OVERFLOW_ERROR(result);
163 }
164 else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
165 UNDERFLOW_ERROR(result);
166 }
167 else {
168 const double sy = GSL_SIGN(y);
169 const double M = floor(x);
170 const double N = floor(ly);
171 const double a = x - M;
172 const double b = ly - N;
173 const double berr = 2.0 * GSL_DBL_EPSILON * (fabs(ly) + fabs(N));
174 result->val = sy * exp(M+N) * exp(a+b);
175 result->err = berr * fabs(result->val);
176 result->err += 2.0 * GSL_DBL_EPSILON * (M + N + 1.0) * fabs(result->val);
177 return GSL_SUCCESS;
178 }
179 }
180 }
181
182
gsl_sf_exp_mult_e10_e(const double x,const double y,gsl_sf_result_e10 * result)183 int gsl_sf_exp_mult_e10_e(const double x, const double y, gsl_sf_result_e10 * result)
184 {
185 const double ay = fabs(y);
186
187 if(y == 0.0) {
188 result->val = 0.0;
189 result->err = 0.0;
190 result->e10 = 0;
191 return GSL_SUCCESS;
192 }
193 else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN)
194 && (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN)
195 ) {
196 const double ex = exp(x);
197 result->val = y * ex;
198 result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
199 result->e10 = 0;
200 return GSL_SUCCESS;
201 }
202 else {
203 const double ly = log(ay);
204 const double l10_val = (x + ly)/M_LN10;
205
206 if(l10_val > INT_MAX-1) {
207 OVERFLOW_ERROR_E10(result);
208 }
209 else if(l10_val < INT_MIN+1) {
210 UNDERFLOW_ERROR_E10(result);
211 }
212 else {
213 const double sy = GSL_SIGN(y);
214 const int N = (int) floor(l10_val);
215 const double arg_val = (l10_val - N) * M_LN10;
216 const double arg_err = 2.0 * GSL_DBL_EPSILON * fabs(ly);
217
218 result->val = sy * exp(arg_val);
219 result->err = arg_err * fabs(result->val);
220 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
221 result->e10 = N;
222
223 return GSL_SUCCESS;
224 }
225 }
226 }
227
228
gsl_sf_exp_mult_err_e(const double x,const double dx,const double y,const double dy,gsl_sf_result * result)229 int gsl_sf_exp_mult_err_e(const double x, const double dx,
230 const double y, const double dy,
231 gsl_sf_result * result)
232 {
233 const double ay = fabs(y);
234
235 if(y == 0.0) {
236 result->val = 0.0;
237 result->err = fabs(dy * exp(x));
238 return GSL_SUCCESS;
239 }
240 else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN)
241 && (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN)
242 ) {
243 double ex = exp(x);
244 result->val = y * ex;
245 result->err = ex * (fabs(dy) + fabs(y*dx));
246 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
247 return GSL_SUCCESS;
248 }
249 else {
250 const double ly = log(ay);
251 const double lnr = x + ly;
252
253 if(lnr > GSL_LOG_DBL_MAX - 0.01) {
254 OVERFLOW_ERROR(result);
255 }
256 else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
257 UNDERFLOW_ERROR(result);
258 }
259 else {
260 const double sy = GSL_SIGN(y);
261 const double M = floor(x);
262 const double N = floor(ly);
263 const double a = x - M;
264 const double b = ly - N;
265 const double eMN = exp(M+N);
266 const double eab = exp(a+b);
267 result->val = sy * eMN * eab;
268 result->err = eMN * eab * 2.0*GSL_DBL_EPSILON;
269 result->err += eMN * eab * fabs(dy/y);
270 result->err += eMN * eab * fabs(dx);
271 return GSL_SUCCESS;
272 }
273 }
274 }
275
276
gsl_sf_exp_mult_err_e10_e(const double x,const double dx,const double y,const double dy,gsl_sf_result_e10 * result)277 int gsl_sf_exp_mult_err_e10_e(const double x, const double dx,
278 const double y, const double dy,
279 gsl_sf_result_e10 * result)
280 {
281 const double ay = fabs(y);
282
283 if(y == 0.0) {
284 result->val = 0.0;
285 result->err = fabs(dy * exp(x));
286 result->e10 = 0;
287 return GSL_SUCCESS;
288 }
289 else if( ( x < 0.5*GSL_LOG_DBL_MAX && x > 0.5*GSL_LOG_DBL_MIN)
290 && (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN)
291 ) {
292 const double ex = exp(x);
293 result->val = y * ex;
294 result->err = ex * (fabs(dy) + fabs(y*dx));
295 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
296 result->e10 = 0;
297 return GSL_SUCCESS;
298 }
299 else {
300 const double ly = log(ay);
301 const double l10_val = (x + ly)/M_LN10;
302
303 if(l10_val > INT_MAX-1) {
304 OVERFLOW_ERROR_E10(result);
305 }
306 else if(l10_val < INT_MIN+1) {
307 UNDERFLOW_ERROR_E10(result);
308 }
309 else {
310 const double sy = GSL_SIGN(y);
311 const int N = (int) floor(l10_val);
312 const double arg_val = (l10_val - N) * M_LN10;
313 const double arg_err = dy/fabs(y) + dx + 2.0*GSL_DBL_EPSILON*fabs(arg_val);
314
315 result->val = sy * exp(arg_val);
316 result->err = arg_err * fabs(result->val);
317 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
318 result->e10 = N;
319
320 return GSL_SUCCESS;
321 }
322 }
323 }
324
325
gsl_sf_expm1_e(const double x,gsl_sf_result * result)326 int gsl_sf_expm1_e(const double x, gsl_sf_result * result)
327 {
328 const double cut = 0.002;
329
330 if(x < GSL_LOG_DBL_MIN) {
331 result->val = -1.0;
332 result->err = GSL_DBL_EPSILON;
333 return GSL_SUCCESS;
334 }
335 else if(x < -cut) {
336 result->val = exp(x) - 1.0;
337 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
338 return GSL_SUCCESS;
339 }
340 else if(x < cut) {
341 result->val = x * (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
342 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
343 return GSL_SUCCESS;
344 }
345 else if(x < GSL_LOG_DBL_MAX) {
346 result->val = exp(x) - 1.0;
347 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
348 return GSL_SUCCESS;
349 }
350 else {
351 OVERFLOW_ERROR(result);
352 }
353 }
354
355
gsl_sf_exprel_e(const double x,gsl_sf_result * result)356 int gsl_sf_exprel_e(const double x, gsl_sf_result * result)
357 {
358 const double cut = 0.002;
359
360 if(x < GSL_LOG_DBL_MIN) {
361 result->val = -1.0/x;
362 result->err = GSL_DBL_EPSILON * fabs(result->val);
363 return GSL_SUCCESS;
364 }
365 else if(x < -cut) {
366 result->val = (exp(x) - 1.0)/x;
367 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
368 return GSL_SUCCESS;
369 }
370 else if(x < cut) {
371 result->val = (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
372 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
373 return GSL_SUCCESS;
374 }
375 else if(x < GSL_LOG_DBL_MAX) {
376 result->val = (exp(x) - 1.0)/x;
377 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
378 return GSL_SUCCESS;
379 }
380 else {
381 OVERFLOW_ERROR(result);
382 }
383 }
384
385
gsl_sf_exprel_2_e(double x,gsl_sf_result * result)386 int gsl_sf_exprel_2_e(double x, gsl_sf_result * result)
387 {
388 const double cut = 0.002;
389
390 if(x < GSL_LOG_DBL_MIN) {
391 result->val = -2.0/x*(1.0 + 1.0/x);
392 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
393 return GSL_SUCCESS;
394 }
395 else if(x < -cut) {
396 result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
397 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
398 return GSL_SUCCESS;
399 }
400 else if(x < cut) {
401 result->val = (1.0 + 1.0/3.0*x*(1.0 + 0.25*x*(1.0 + 0.2*x*(1.0 + 1.0/6.0*x))));
402 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
403 return GSL_SUCCESS;
404 }
405 else if(x < GSL_LOG_DBL_MAX) {
406 result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
407 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
408 return GSL_SUCCESS;
409 }
410 else {
411 OVERFLOW_ERROR(result);
412 }
413 }
414
415
416 int
gsl_sf_exprel_n_e(const int N,const double x,gsl_sf_result * result)417 gsl_sf_exprel_n_e(const int N, const double x, gsl_sf_result * result)
418 {
419 if(N < 0) {
420 DOMAIN_ERROR(result);
421 }
422 else if(x == 0.0) {
423 result->val = 1.0;
424 result->err = 0.0;
425 return GSL_SUCCESS;
426 }
427 else if(fabs(x) < GSL_ROOT3_DBL_EPSILON * N) {
428 result->val = 1.0 + x/(N+1) * (1.0 + x/(N+2));
429 result->err = 2.0 * GSL_DBL_EPSILON;
430 return GSL_SUCCESS;
431 }
432 else if(N == 0) {
433 return gsl_sf_exp_e(x, result);
434 }
435 else if(N == 1) {
436 return gsl_sf_exprel_e(x, result);
437 }
438 else if(N == 2) {
439 return gsl_sf_exprel_2_e(x, result);
440 }
441 else {
442 if(x > N && (-x + N*(1.0 + log(x/N)) < GSL_LOG_DBL_EPSILON)) {
443 /* x is much larger than n.
444 * Ignore polynomial part, so
445 * exprel_N(x) ~= e^x N!/x^N
446 */
447 gsl_sf_result lnf_N;
448 double lnr_val;
449 double lnr_err;
450 double lnterm;
451 gsl_sf_lnfact_e(N, &lnf_N);
452 lnterm = N*log(x);
453 lnr_val = x + lnf_N.val - lnterm;
454 lnr_err = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(lnterm));
455 lnr_err += lnf_N.err;
456 return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
457 }
458 else if(x > N) {
459 /* Write the identity
460 * exprel_n(x) = e^x n! / x^n (1 - Gamma[n,x]/Gamma[n])
461 * then use the asymptotic expansion
462 * Gamma[n,x] ~ x^(n-1) e^(-x) (1 + (n-1)/x + (n-1)(n-2)/x^2 + ...)
463 */
464 double ln_x = log(x);
465 gsl_sf_result lnf_N;
466 double lg_N;
467 double lnpre_val;
468 double lnpre_err;
469 gsl_sf_lnfact_e(N, &lnf_N); /* log(N!) */
470 lg_N = lnf_N.val - log(N); /* log(Gamma(N)) */
471 lnpre_val = x + lnf_N.val - N*ln_x;
472 lnpre_err = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(N*ln_x));
473 lnpre_err += lnf_N.err;
474 if(lnpre_val < GSL_LOG_DBL_MAX - 5.0) {
475 int stat_eG;
476 gsl_sf_result bigG_ratio;
477 gsl_sf_result pre;
478 int stat_ex = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &pre);
479 double ln_bigG_ratio_pre = -x + (N-1)*ln_x - lg_N;
480 double bigGsum = 1.0;
481 double term = 1.0;
482 int k;
483 for(k=1; k<N; k++) {
484 term *= (N-k)/x;
485 bigGsum += term;
486 }
487 stat_eG = gsl_sf_exp_mult_e(ln_bigG_ratio_pre, bigGsum, &bigG_ratio);
488 if(stat_eG == GSL_SUCCESS) {
489 result->val = pre.val * (1.0 - bigG_ratio.val);
490 result->err = pre.val * (2.0*GSL_DBL_EPSILON + bigG_ratio.err);
491 result->err += pre.err * fabs(1.0 - bigG_ratio.val);
492 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
493 return stat_ex;
494 }
495 else {
496 result->val = 0.0;
497 result->err = 0.0;
498 return stat_eG;
499 }
500 }
501 else {
502 OVERFLOW_ERROR(result);
503 }
504 }
505 else if(x > -10.0*N) {
506 return exprel_n_CF(N, x, result);
507 }
508 else {
509 /* x -> -Inf asymptotic:
510 * exprel_n(x) ~ e^x n!/x^n - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
511 * ~ - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
512 */
513 double sum = 1.0;
514 double term = 1.0;
515 int k;
516 for(k=1; k<N; k++) {
517 term *= (N-k)/x;
518 sum += term;
519 }
520 result->val = -N/x * sum;
521 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
522 return GSL_SUCCESS;
523 }
524 }
525 }
526
527
528 int
gsl_sf_exp_err_e(const double x,const double dx,gsl_sf_result * result)529 gsl_sf_exp_err_e(const double x, const double dx, gsl_sf_result * result)
530 {
531 const double adx = fabs(dx);
532
533 /* CHECK_POINTER(result) */
534
535 if(x + adx > GSL_LOG_DBL_MAX) {
536 OVERFLOW_ERROR(result);
537 }
538 else if(x - adx < GSL_LOG_DBL_MIN) {
539 UNDERFLOW_ERROR(result);
540 }
541 else {
542 const double ex = exp(x);
543 const double edx = exp(adx);
544 result->val = ex;
545 result->err = ex * GSL_MAX_DBL(GSL_DBL_EPSILON, edx - 1.0/edx);
546 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
547 return GSL_SUCCESS;
548 }
549 }
550
551
552 int
gsl_sf_exp_err_e10_e(const double x,const double dx,gsl_sf_result_e10 * result)553 gsl_sf_exp_err_e10_e(const double x, const double dx, gsl_sf_result_e10 * result)
554 {
555 const double adx = fabs(dx);
556
557 /* CHECK_POINTER(result) */
558
559 if(x + adx > INT_MAX - 1) {
560 OVERFLOW_ERROR_E10(result);
561 }
562 else if(x - adx < INT_MIN + 1) {
563 UNDERFLOW_ERROR_E10(result);
564 }
565 else {
566 const int N = (int)floor(x/M_LN10);
567 const double ex = exp(x-N*M_LN10);
568 result->val = ex;
569 result->err = ex * (2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) + adx);
570 result->e10 = N;
571 return GSL_SUCCESS;
572 }
573 }
574
575
576 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
577
578 #include "gsl_specfunc__eval.h"
579
gsl_sf_exp(const double x)580 double gsl_sf_exp(const double x)
581 {
582 EVAL_RESULT(gsl_sf_exp_e(x, &result));
583 }
584
gsl_sf_exp_mult(const double x,const double y)585 double gsl_sf_exp_mult(const double x, const double y)
586 {
587 EVAL_RESULT(gsl_sf_exp_mult_e(x, y, &result));
588 }
589
gsl_sf_expm1(const double x)590 double gsl_sf_expm1(const double x)
591 {
592 EVAL_RESULT(gsl_sf_expm1_e(x, &result));
593 }
594
gsl_sf_exprel(const double x)595 double gsl_sf_exprel(const double x)
596 {
597 EVAL_RESULT(gsl_sf_exprel_e(x, &result));
598 }
599
gsl_sf_exprel_2(const double x)600 double gsl_sf_exprel_2(const double x)
601 {
602 EVAL_RESULT(gsl_sf_exprel_2_e(x, &result));
603 }
604
gsl_sf_exprel_n(const int n,const double x)605 double gsl_sf_exprel_n(const int n, const double x)
606 {
607 EVAL_RESULT(gsl_sf_exprel_n_e(n, x, &result));
608 }
609