← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/bogofilter/oneiric-updates

« back to all changes in this revision

Viewing changes to gsl/specfunc/.svn/text-base/exp.c.svn-base

  • Committer: Bazaar Package Importer
  • Author(s): Loïc Minier
  • Date: 2010-04-10 11:08:53 UTC
  • mfrom: (1.1.12 upstream)
  • Revision ID: james.westby@ubuntu.com-20100410110853-5w0vep9aa9qqr9lx
Tags: 1.2.1-0ubuntu1
https://launchpad.net/bugs/557468
https://launchpad.net/bugs/320829
* New upstream bugfix release; LP: #557468.
  - Fixes parsing of the first line of the body in MIME messages;
    LP: #320829.

Show diffs side-by-side

added added

removed removed

Lines of Context:
gsl/specfunc/.svn/text-base/exp.c.svn-base
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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
18
 
 */
19
 
 
20
 
/* Author:  G. Jungman */
21
 
 
22
 
#include <config.h>
23
 
#include <gsl/gsl_math.h>
24
 
#include <gsl/gsl_errno.h>
25
 
#include <gsl/gsl_sf_gamma.h>
26
 
#include <gsl/gsl_sf_exp.h>
27
 
 
28
 
#include "error.h"
29
 
 
30
 
/* Evaluate the continued fraction for exprel.
31
 
 * [Abramowitz+Stegun, 4.2.41]
32
 
 */
33
 
static
34
 
int
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
 
 
107
 
#ifndef HIDE_INLINE_STATIC
108
 
int gsl_sf_exp_e(const double x, gsl_sf_result * result)
109
 
{
110
 
  if(x > GSL_LOG_DBL_MAX) {
111
 
    OVERFLOW_ERROR(result);
112
 
  }
113
 
  else if(x < GSL_LOG_DBL_MIN) {
114
 
    UNDERFLOW_ERROR(result);
115
 
  }
116
 
  else {
117
 
    result->val = exp(x);
118
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
119
 
    return GSL_SUCCESS;
120
 
  }
121
 
}
122
 
#endif
123
 
 
124
 
int gsl_sf_exp_e10_e(const double x, gsl_sf_result_e10 * result)
125
 
{
126
 
  if(x > INT_MAX-1) {
127
 
    OVERFLOW_ERROR_E10(result);
128
 
  }
129
 
  else if(x < INT_MIN+1) {
130
 
    UNDERFLOW_ERROR_E10(result);
131
 
  }
132
 
  else {
133
 
    const int N = (int) floor(x/M_LN10);
134
 
    result->val = exp(x-N*M_LN10);
135
 
    result->err = 2.0 * (fabs(x)+1.0) * GSL_DBL_EPSILON * fabs(result->val);
136
 
    result->e10 = N;
137
 
    return GSL_SUCCESS;
138
 
  }
139
 
}
140
 
 
141
 
 
142
 
int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result)
143
 
{
144
 
  const double ay  = fabs(y);
145
 
 
146
 
  if(y == 0.0) {
147
 
    result->val = 0.0;
148
 
    result->err = 0.0;
149
 
    return GSL_SUCCESS;
150
 
  }
151
 
  else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
152
 
          && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
153
 
    ) {
154
 
    const double ex = exp(x);
155
 
    result->val = y * ex;
156
 
    result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
157
 
    return GSL_SUCCESS;
158
 
  }
159
 
  else {
160
 
    const double ly  = log(ay);
161
 
    const double lnr = x + ly;
162
 
 
163
 
    if(lnr > GSL_LOG_DBL_MAX - 0.01) {
164
 
      OVERFLOW_ERROR(result);
165
 
    }
166
 
    else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
167
 
      UNDERFLOW_ERROR(result);
168
 
    }
169
 
    else {
170
 
      const double sy   = GSL_SIGN(y);
171
 
      const double M    = floor(x);
172
 
      const double N    = floor(ly);
173
 
      const double a    = x  - M;
174
 
      const double b    = ly - N;
175
 
      const double berr = 2.0 * GSL_DBL_EPSILON * (fabs(ly) + fabs(N));
176
 
      result->val  = sy * exp(M+N) * exp(a+b);
177
 
      result->err  = berr * fabs(result->val);
178
 
      result->err += 2.0 * GSL_DBL_EPSILON * (M + N + 1.0) * fabs(result->val);
179
 
      return GSL_SUCCESS;
180
 
    }
181
 
  }
182
 
}
183
 
 
184
 
 
185
 
int gsl_sf_exp_mult_e10_e(const double x, const double y, gsl_sf_result_e10 * result)
186
 
{
187
 
  const double ay  = fabs(y);
188
 
 
189
 
  if(y == 0.0) {
190
 
    result->val = 0.0;
191
 
    result->err = 0.0;
192
 
    result->e10 = 0;
193
 
    return GSL_SUCCESS;
194
 
  }
195
 
  else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
196
 
          && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
197
 
    ) {
198
 
    const double ex = exp(x);
199
 
    result->val = y * ex;
200
 
    result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
201
 
    result->e10 = 0;
202
 
    return GSL_SUCCESS;
203
 
  }
204
 
  else {
205
 
    const double ly  = log(ay);
206
 
    const double l10_val = (x + ly)/M_LN10;
207
 
 
208
 
    if(l10_val > INT_MAX-1) {
209
 
      OVERFLOW_ERROR_E10(result);
210
 
    }
211
 
    else if(l10_val < INT_MIN+1) {
212
 
      UNDERFLOW_ERROR_E10(result);
213
 
    }
214
 
    else {
215
 
      const double sy  = GSL_SIGN(y);
216
 
      const int    N   = (int) floor(l10_val);
217
 
      const double arg_val = (l10_val - N) * M_LN10;
218
 
      const double arg_err = 2.0 * GSL_DBL_EPSILON * fabs(ly);
219
 
 
220
 
      result->val  = sy * exp(arg_val);
221
 
      result->err  = arg_err * fabs(result->val);
222
 
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
223
 
      result->e10 = N;
224
 
 
225
 
      return GSL_SUCCESS;
226
 
    }
227
 
  }
228
 
}
229
 
 
230
 
 
231
 
int gsl_sf_exp_mult_err_e(const double x, const double dx,
232
 
                             const double y, const double dy,
233
 
                             gsl_sf_result * result)
234
 
{
235
 
  const double ay  = fabs(y);
236
 
 
237
 
  if(y == 0.0) {
238
 
    result->val = 0.0;
239
 
    result->err = fabs(dy * exp(x));
240
 
    return GSL_SUCCESS;
241
 
  }
242
 
  else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
243
 
          && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
244
 
    ) {
245
 
    double ex = exp(x);
246
 
    result->val  = y * ex;
247
 
    result->err  = ex * (fabs(dy) + fabs(y*dx));
248
 
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
249
 
    return GSL_SUCCESS;
250
 
  }
251
 
  else {
252
 
    const double ly  = log(ay);
253
 
    const double lnr = x + ly;
254
 
 
255
 
    if(lnr > GSL_LOG_DBL_MAX - 0.01) {
256
 
      OVERFLOW_ERROR(result);
257
 
    }
258
 
    else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
259
 
      UNDERFLOW_ERROR(result);
260
 
    }
261
 
    else {
262
 
      const double sy  = GSL_SIGN(y);
263
 
      const double M   = floor(x);
264
 
      const double N   = floor(ly);
265
 
      const double a   = x  - M;
266
 
      const double b   = ly - N;
267
 
      const double eMN = exp(M+N);
268
 
      const double eab = exp(a+b);
269
 
      result->val  = sy * eMN * eab;
270
 
      result->err  = eMN * eab * 2.0*GSL_DBL_EPSILON;
271
 
      result->err += eMN * eab * fabs(dy/y);
272
 
      result->err += eMN * eab * fabs(dx);
273
 
      return GSL_SUCCESS;
274
 
    }
275
 
  }
276
 
}
277
 
 
278
 
 
279
 
int gsl_sf_exp_mult_err_e10_e(const double x, const double dx,
280
 
                             const double y, const double dy,
281
 
                             gsl_sf_result_e10 * result)
282
 
{
283
 
  const double ay  = fabs(y);
284
 
 
285
 
  if(y == 0.0) {
286
 
    result->val = 0.0;
287
 
    result->err = fabs(dy * exp(x));
288
 
    result->e10 = 0;
289
 
    return GSL_SUCCESS;
290
 
  }
291
 
  else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
292
 
          && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
293
 
    ) {
294
 
    const double ex = exp(x);
295
 
    result->val  = y * ex;
296
 
    result->err  = ex * (fabs(dy) + fabs(y*dx));
297
 
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
298
 
    result->e10 = 0;
299
 
    return GSL_SUCCESS;
300
 
  }
301
 
  else {
302
 
    const double ly  = log(ay);
303
 
    const double l10_val = (x + ly)/M_LN10;
304
 
 
305
 
    if(l10_val > INT_MAX-1) {
306
 
      OVERFLOW_ERROR_E10(result);
307
 
    }
308
 
    else if(l10_val < INT_MIN+1) {
309
 
      UNDERFLOW_ERROR_E10(result);
310
 
    }
311
 
    else {
312
 
      const double sy  = GSL_SIGN(y);
313
 
      const int    N   = (int) floor(l10_val);
314
 
      const double arg_val = (l10_val - N) * M_LN10;
315
 
      const double arg_err = dy/fabs(y) + dx + 2.0*GSL_DBL_EPSILON*fabs(arg_val);
316
 
 
317
 
      result->val  = sy * exp(arg_val);
318
 
      result->err  = arg_err * fabs(result->val);
319
 
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
320
 
      result->e10 = N;
321
 
 
322
 
      return GSL_SUCCESS;
323
 
    }
324
 
  }
325
 
}
326
 
 
327
 
 
328
 
int gsl_sf_expm1_e(const double x, gsl_sf_result * result)
329
 
{
330
 
  const double cut = 0.002;
331
 
 
332
 
  if(x < GSL_LOG_DBL_MIN) {
333
 
    result->val = -1.0;
334
 
    result->err = GSL_DBL_EPSILON;
335
 
    return GSL_SUCCESS;
336
 
  }
337
 
  else if(x < -cut) {
338
 
    result->val = exp(x) - 1.0;
339
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
340
 
    return GSL_SUCCESS;
341
 
  }
342
 
  else if(x < cut) {
343
 
    result->val = x * (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
344
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
345
 
    return GSL_SUCCESS;
346
 
  } 
347
 
  else if(x < GSL_LOG_DBL_MAX) {
348
 
    result->val = exp(x) - 1.0;
349
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
350
 
    return GSL_SUCCESS;
351
 
  }
352
 
  else {
353
 
    OVERFLOW_ERROR(result);
354
 
  }
355
 
}
356
 
 
357
 
 
358
 
int gsl_sf_exprel_e(const double x, gsl_sf_result * result)
359
 
{
360
 
  const double cut = 0.002;
361
 
 
362
 
  if(x < GSL_LOG_DBL_MIN) {
363
 
    result->val = -1.0/x;
364
 
    result->err = GSL_DBL_EPSILON * fabs(result->val);
365
 
    return GSL_SUCCESS;
366
 
  }
367
 
  else if(x < -cut) {
368
 
    result->val = (exp(x) - 1.0)/x;
369
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
370
 
    return GSL_SUCCESS;
371
 
  }
372
 
  else if(x < cut) {
373
 
    result->val = (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
374
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
375
 
    return GSL_SUCCESS;
376
 
  } 
377
 
  else if(x < GSL_LOG_DBL_MAX) {
378
 
    result->val = (exp(x) - 1.0)/x;
379
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
380
 
    return GSL_SUCCESS;
381
 
  }
382
 
  else {
383
 
    OVERFLOW_ERROR(result);
384
 
  }
385
 
}
386
 
 
387
 
 
388
 
int gsl_sf_exprel_2_e(double x, gsl_sf_result * result)
389
 
{
390
 
  const double cut = 0.002;
391
 
 
392
 
  if(x < GSL_LOG_DBL_MIN) {
393
 
    result->val = -2.0/x*(1.0 + 1.0/x);
394
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
395
 
    return GSL_SUCCESS;
396
 
  }
397
 
  else if(x < -cut) {
398
 
    result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
399
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
400
 
    return GSL_SUCCESS;
401
 
  }
402
 
  else if(x < cut) {
403
 
    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))));
404
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
405
 
    return GSL_SUCCESS;
406
 
  } 
407
 
  else if(x < GSL_LOG_DBL_MAX) {
408
 
    result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
409
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
410
 
    return GSL_SUCCESS;
411
 
  }
412
 
  else {
413
 
    OVERFLOW_ERROR(result);
414
 
  }
415
 
}
416
 
 
417
 
 
418
 
int
419
 
gsl_sf_exprel_n_e(const int N, const double x, gsl_sf_result * result)
420
 
{
421
 
  if(N < 0) {
422
 
    DOMAIN_ERROR(result);
423
 
  }
424
 
  else if(x == 0.0) {
425
 
    result->val = 1.0;
426
 
    result->err = 0.0;
427
 
    return GSL_SUCCESS;
428
 
  }
429
 
  else if(fabs(x) < GSL_ROOT3_DBL_EPSILON * N) {
430
 
    result->val = 1.0 + x/(N+1) * (1.0 + x/(N+2));
431
 
    result->err = 2.0 * GSL_DBL_EPSILON;
432
 
    return GSL_SUCCESS;
433
 
  }
434
 
  else if(N == 0) {
435
 
    return gsl_sf_exp_e(x, result);
436
 
  }
437
 
  else if(N == 1) {
438
 
    return gsl_sf_exprel_e(x, result);
439
 
  }
440
 
  else if(N == 2) {
441
 
    return gsl_sf_exprel_2_e(x, result);
442
 
  }
443
 
  else {
444
 
    if(x > N && (-x + N*(1.0 + log(x/N)) < GSL_LOG_DBL_EPSILON)) {
445
 
      /* x is much larger than n.
446
 
       * Ignore polynomial part, so
447
 
       * exprel_N(x) ~= e^x N!/x^N
448
 
       */
449
 
      gsl_sf_result lnf_N;
450
 
      double lnr_val;
451
 
      double lnr_err;
452
 
      double lnterm;
453
 
      gsl_sf_lnfact_e(N, &lnf_N);
454
 
      lnterm = N*log(x);
455
 
      lnr_val  = x + lnf_N.val - lnterm;
456
 
      lnr_err  = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(lnterm));
457
 
      lnr_err += lnf_N.err;
458
 
      return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
459
 
    }
460
 
    else if(x > N) {
461
 
      /* Write the identity
462
 
       *   exprel_n(x) = e^x n! / x^n (1 - Gamma[n,x]/Gamma[n])
463
 
       * then use the asymptotic expansion
464
 
       * Gamma[n,x] ~ x^(n-1) e^(-x) (1 + (n-1)/x + (n-1)(n-2)/x^2 + ...)
465
 
       */
466
 
      double ln_x = log(x);
467
 
      gsl_sf_result lnf_N;
468
 
      double lg_N;
469
 
      double lnpre_val;
470
 
      double lnpre_err;
471
 
      gsl_sf_lnfact_e(N, &lnf_N);    /* log(N!)       */
472
 
      lg_N  = lnf_N.val - log(N);       /* log(Gamma(N)) */
473
 
      lnpre_val  = x + lnf_N.val - N*ln_x;
474
 
      lnpre_err  = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(N*ln_x));
475
 
      lnpre_err += lnf_N.err;
476
 
      if(lnpre_val < GSL_LOG_DBL_MAX - 5.0) {
477
 
        int stat_eG;
478
 
        gsl_sf_result bigG_ratio;
479
 
        gsl_sf_result pre;
480
 
        int stat_ex = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &pre);
481
 
        double ln_bigG_ratio_pre = -x + (N-1)*ln_x - lg_N;
482
 
        double bigGsum = 1.0;
483
 
        double term = 1.0;
484
 
        int k;
485
 
        for(k=1; k<N; k++) {
486
 
          term *= (N-k)/x;
487
 
          bigGsum += term;
488
 
        }
489
 
        stat_eG = gsl_sf_exp_mult_e(ln_bigG_ratio_pre, bigGsum, &bigG_ratio);
490
 
        if(stat_eG == GSL_SUCCESS) {
491
 
          result->val  = pre.val * (1.0 - bigG_ratio.val);
492
 
          result->err  = pre.val * (2.0*GSL_DBL_EPSILON + bigG_ratio.err);
493
 
          result->err += pre.err * fabs(1.0 - bigG_ratio.val);
494
 
          result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
495
 
          return stat_ex;
496
 
        }
497
 
        else {
498
 
          result->val = 0.0;
499
 
          result->err = 0.0;
500
 
          return stat_eG;
501
 
        }
502
 
      }
503
 
      else {
504
 
        OVERFLOW_ERROR(result);
505
 
      }
506
 
    }
507
 
    else if(x > -10.0*N) {
508
 
      return exprel_n_CF(N, x, result);
509
 
    }
510
 
    else {
511
 
      /* x -> -Inf asymptotic:
512
 
       * exprel_n(x) ~ e^x n!/x^n - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
513
 
       *             ~ - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
514
 
       */
515
 
      double sum  = 1.0;
516
 
      double term = 1.0;
517
 
      int k;
518
 
      for(k=1; k<N; k++) {
519
 
        term *= (N-k)/x;
520
 
        sum  += term;
521
 
      }
522
 
      result->val = -N/x * sum;
523
 
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
524
 
      return GSL_SUCCESS;
525
 
    }
526
 
  }
527
 
}
528
 
 
529
 
 
530
 
int
531
 
gsl_sf_exp_err_e(const double x, const double dx, gsl_sf_result * result)
532
 
{
533
 
  const double adx = fabs(dx);
534
 
 
535
 
  /* CHECK_POINTER(result) */
536
 
 
537
 
  if(x + adx > GSL_LOG_DBL_MAX) {
538
 
    OVERFLOW_ERROR(result);
539
 
  }
540
 
  else if(x - adx < GSL_LOG_DBL_MIN) {
541
 
    UNDERFLOW_ERROR(result);
542
 
  }
543
 
  else {
544
 
    const double ex  = exp(x);
545
 
    const double edx = exp(adx);
546
 
    result->val  = ex;
547
 
    result->err  = ex * GSL_MAX_DBL(GSL_DBL_EPSILON, edx - 1.0/edx);
548
 
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
549
 
    return GSL_SUCCESS;
550
 
  }
551
 
}
552
 
 
553
 
 
554
 
int
555
 
gsl_sf_exp_err_e10_e(const double x, const double dx, gsl_sf_result_e10 * result)
556
 
{
557
 
  const double adx = fabs(dx);
558
 
 
559
 
  /* CHECK_POINTER(result) */
560
 
 
561
 
  if(x + adx > INT_MAX - 1) {
562
 
    OVERFLOW_ERROR_E10(result);
563
 
  }
564
 
  else if(x - adx < INT_MIN + 1) {
565
 
    UNDERFLOW_ERROR_E10(result);
566
 
  }
567
 
  else {
568
 
    const int    N  = (int)floor(x/M_LN10);
569
 
    const double ex = exp(x-N*M_LN10);
570
 
    result->val = ex;
571
 
    result->err = ex * (2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) + adx);
572
 
    result->e10 = N;
573
 
    return GSL_SUCCESS;
574
 
  }
575
 
}
576
 
 
577
 
 
578
 
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
579
 
 
580
 
#include "eval.h"
581
 
 
582
 
double gsl_sf_exp(const double x)
583
 
{
584
 
  EVAL_RESULT(gsl_sf_exp_e(x, &result));
585
 
}
586
 
 
587
 
double gsl_sf_exp_mult(const double x, const double y)
588
 
{
589
 
  EVAL_RESULT(gsl_sf_exp_mult_e(x, y, &result));
590
 
}
591
 
 
592
 
double gsl_sf_expm1(const double x)
593
 
{
594
 
  EVAL_RESULT(gsl_sf_expm1_e(x, &result));
595
 
}
596
 
 
597
 
double gsl_sf_exprel(const double x)
598
 
{
599
 
  EVAL_RESULT(gsl_sf_exprel_e(x, &result));
600
 
}
601
 
 
602
 
double gsl_sf_exprel_2(const double x)
603
 
{
604
 
  EVAL_RESULT(gsl_sf_exprel_2_e(x, &result));
605
 
}
606
 
 
607
 
double gsl_sf_exprel_n(const int n, const double x)
608
 
{
609
 
  EVAL_RESULT(gsl_sf_exprel_n_e(n, x, &result));
610
 
}