← Back to branch summary

~ubuntu-branches/ubuntu/lucid/bogofilter/lucid-updates

« back to all changes in this revision

Viewing changes to gsl/specfunc/.svn/text-base/erfc.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/erfc.c.svn-base
1
 
/* specfunc/erfc.c
2
 
 * 
3
 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003 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:  J. Theiler (modifications by G. Jungman) */
21
 
 
22
 
/*
23
 
 * See Hart et al, Computer Approximations, John Wiley and Sons, New York (1968)
24
 
 * (This applies only to the erfc8 stuff, which is the part
25
 
 *  of the original code that survives. I have replaced much of
26
 
 *  the other stuff with Chebyshev fits. These are simpler and
27
 
 *  more precise than the original approximations. [GJ])
28
 
 */
29
 
#include <config.h>
30
 
#include <gsl/gsl_math.h>
31
 
#include <gsl/gsl_errno.h>
32
 
#include <gsl/gsl_sf_exp.h>
33
 
#include <gsl/gsl_sf_erf.h>
34
 
 
35
 
#include "check.h"
36
 
 
37
 
#include "chebyshev.h"
38
 
#include "cheb_eval.c"
39
 
 
40
 
#define LogRootPi_  0.57236494292470008706
41
 
 
42
 
 
43
 
static double erfc8_sum(double x)
44
 
{
45
 
  /* estimates erfc(x) valid for 8 < x < 100 */
46
 
  /* This is based on index 5725 in Hart et al */
47
 
 
48
 
  static double P[] = {
49
 
      2.97886562639399288862,
50
 
      7.409740605964741794425,
51
 
      6.1602098531096305440906,
52
 
      5.019049726784267463450058,
53
 
      1.275366644729965952479585264,
54
 
      0.5641895835477550741253201704
55
 
  };
56
 
  static double Q[] = {
57
 
      3.3690752069827527677,
58
 
      9.608965327192787870698,
59
 
      17.08144074746600431571095,
60
 
      12.0489519278551290360340491,
61
 
      9.396034016235054150430579648,
62
 
      2.260528520767326969591866945,
63
 
      1.0
64
 
  };
65
 
  double num=0.0, den=0.0;
66
 
  int i;
67
 
 
68
 
  num = P[5];
69
 
  for (i=4; i>=0; --i) {
70
 
      num = x*num + P[i];
71
 
  }
72
 
  den = Q[6];
73
 
  for (i=5; i>=0; --i) {
74
 
      den = x*den + Q[i];
75
 
  }
76
 
 
77
 
  return num/den;
78
 
}
79
 
 
80
 
inline
81
 
static double erfc8(double x)
82
 
{
83
 
  double e;
84
 
  e = erfc8_sum(x);
85
 
  e *= exp(-x*x);
86
 
  return e;
87
 
}
88
 
 
89
 
inline
90
 
static double log_erfc8(double x)
91
 
{
92
 
  double e;
93
 
  e = erfc8_sum(x);
94
 
  e = log(e) - x*x;
95
 
  return e;
96
 
}
97
 
 
98
 
#if 0
99
 
/* Abramowitz+Stegun, 7.2.14 */
100
 
static double erfcasympsum(double x)
101
 
{
102
 
  int i;
103
 
  double e = 1.;
104
 
  double coef = 1.;
105
 
  for (i=1; i<5; ++i) {
106
 
    /* coef *= -(2*i-1)/(2*x*x); ??? [GJ] */
107
 
    coef *= -(2*i+1)/(i*(4*x*x*x*x));
108
 
    e += coef;
109
 
    /*
110
 
    if (fabs(coef) < 1.0e-15) break;
111
 
    if (fabs(coef) > 1.0e10) break;
112
 
    
113
 
    [GJ]: These tests are not useful. This function is only
114
 
    used below. Took them out; they gum up the pipeline.
115
 
    */
116
 
  }
117
 
  return e;
118
 
}
119
 
#endif /* 0 */
120
 
 
121
 
 
122
 
/* Abramowitz+Stegun, 7.1.5 */
123
 
static int erfseries(double x, gsl_sf_result * result)
124
 
{
125
 
  double coef = x;
126
 
  double e    = coef;
127
 
  double del;
128
 
  int k;
129
 
  for (k=1; k<30; ++k) {
130
 
    coef *= -x*x/k;
131
 
    del   = coef/(2.0*k+1.0);
132
 
    e += del;
133
 
  }
134
 
  result->val = 2.0 / M_SQRTPI * e;
135
 
  result->err = 2.0 / M_SQRTPI * (fabs(del) + GSL_DBL_EPSILON);
136
 
  return GSL_SUCCESS;
137
 
}
138
 
 
139
 
 
140
 
/* Chebyshev fit for erfc((t+1)/2), -1 < t < 1
141
 
 */
142
 
static double erfc_xlt1_data[20] = {
143
 
  1.06073416421769980345174155056,
144
 
 -0.42582445804381043569204735291,
145
 
  0.04955262679620434040357683080,
146
 
  0.00449293488768382749558001242,
147
 
 -0.00129194104658496953494224761,
148
 
 -0.00001836389292149396270416979,
149
 
  0.00002211114704099526291538556,
150
 
 -5.23337485234257134673693179020e-7,
151
 
 -2.78184788833537885382530989578e-7,
152
 
  1.41158092748813114560316684249e-8,
153
 
  2.72571296330561699984539141865e-9,
154
 
 -2.06343904872070629406401492476e-10,
155
 
 -2.14273991996785367924201401812e-11,
156
 
  2.22990255539358204580285098119e-12,
157
 
  1.36250074650698280575807934155e-13,
158
 
 -1.95144010922293091898995913038e-14,
159
 
 -6.85627169231704599442806370690e-16,
160
 
  1.44506492869699938239521607493e-16,
161
 
  2.45935306460536488037576200030e-18,
162
 
 -9.29599561220523396007359328540e-19
163
 
};
164
 
static cheb_series erfc_xlt1_cs = {
165
 
  erfc_xlt1_data,
166
 
  19,
167
 
  -1, 1,
168
 
  12
169
 
};
170
 
 
171
 
/* Chebyshev fit for erfc(x) exp(x^2), 1 < x < 5, x = 2t + 3, -1 < t < 1
172
 
 */
173
 
static double erfc_x15_data[25] = {
174
 
  0.44045832024338111077637466616,
175
 
 -0.143958836762168335790826895326,
176
 
  0.044786499817939267247056666937,
177
 
 -0.013343124200271211203618353102,
178
 
  0.003824682739750469767692372556,
179
 
 -0.001058699227195126547306482530,
180
 
  0.000283859419210073742736310108,
181
 
 -0.000073906170662206760483959432,
182
 
  0.000018725312521489179015872934,
183
 
 -4.62530981164919445131297264430e-6,
184
 
  1.11558657244432857487884006422e-6,
185
 
 -2.63098662650834130067808832725e-7,
186
 
  6.07462122724551777372119408710e-8,
187
 
 -1.37460865539865444777251011793e-8,
188
 
  3.05157051905475145520096717210e-9,
189
 
 -6.65174789720310713757307724790e-10,
190
 
  1.42483346273207784489792999706e-10,
191
 
 -3.00141127395323902092018744545e-11,
192
 
  6.22171792645348091472914001250e-12,
193
 
 -1.26994639225668496876152836555e-12,
194
 
  2.55385883033257575402681845385e-13,
195
 
 -5.06258237507038698392265499770e-14,
196
 
  9.89705409478327321641264227110e-15,
197
 
 -1.90685978789192181051961024995e-15,
198
 
  3.50826648032737849245113757340e-16
199
 
};
200
 
static cheb_series erfc_x15_cs = {
201
 
  erfc_x15_data,
202
 
  24,
203
 
  -1, 1,
204
 
  16
205
 
};
206
 
 
207
 
/* Chebyshev fit for erfc(x) x exp(x^2), 5 < x < 10, x = (5t + 15)/2, -1 < t < 1
208
 
 */
209
 
static double erfc_x510_data[20] = {
210
 
  1.11684990123545698684297865808,
211
 
  0.003736240359381998520654927536,
212
 
 -0.000916623948045470238763619870,
213
 
  0.000199094325044940833965078819,
214
 
 -0.000040276384918650072591781859,
215
 
  7.76515264697061049477127605790e-6,
216
 
 -1.44464794206689070402099225301e-6,
217
 
  2.61311930343463958393485241947e-7,
218
 
 -4.61833026634844152345304095560e-8,
219
 
  8.00253111512943601598732144340e-9,
220
 
 -1.36291114862793031395712122089e-9,
221
 
  2.28570483090160869607683087722e-10,
222
 
 -3.78022521563251805044056974560e-11,
223
 
  6.17253683874528285729910462130e-12,
224
 
 -9.96019290955316888445830597430e-13,
225
 
  1.58953143706980770269506726000e-13,
226
 
 -2.51045971047162509999527428316e-14,
227
 
  3.92607828989125810013581287560e-15,
228
 
 -6.07970619384160374392535453420e-16,
229
 
  9.12600607264794717315507477670e-17
230
 
};
231
 
static cheb_series erfc_x510_cs = {
232
 
  erfc_x510_data,
233
 
  19,
234
 
  -1, 1,
235
 
  12
236
 
};
237
 
 
238
 
#if 0
239
 
inline
240
 
static double
241
 
erfc_asymptotic(double x)
242
 
{
243
 
  return exp(-x*x)/x * erfcasympsum(x) / M_SQRTPI;
244
 
}
245
 
inline
246
 
static double
247
 
log_erfc_asymptotic(double x)
248
 
{
249
 
  return log(erfcasympsum(x)/x) - x*x - LogRootPi_;
250
 
}
251
 
#endif /* 0 */
252
 
 
253
 
 
254
 
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
255
 
 
256
 
int gsl_sf_erfc_e(double x, gsl_sf_result * result)
257
 
{
258
 
  const double ax = fabs(x);
259
 
  double e_val, e_err;
260
 
 
261
 
  /* CHECK_POINTER(result) */
262
 
 
263
 
  if(ax <= 1.0) {
264
 
    double t = 2.0*ax - 1.0;
265
 
    gsl_sf_result c;
266
 
    cheb_eval_e(&erfc_xlt1_cs, t, &c);
267
 
    e_val = c.val;
268
 
    e_err = c.err;
269
 
  }
270
 
  else if(ax <= 5.0) {
271
 
    double ex2 = exp(-x*x);
272
 
    double t = 0.5*(ax-3.0);
273
 
    gsl_sf_result c;
274
 
    cheb_eval_e(&erfc_x15_cs, t, &c);
275
 
    e_val = ex2 * c.val;
276
 
    e_err = ex2 * (c.err + 2.0*fabs(x)*GSL_DBL_EPSILON);
277
 
  }
278
 
  else if(ax < 10.0) {
279
 
    double exterm = exp(-x*x) / ax;
280
 
    double t = (2.0*ax - 15.0)/5.0;
281
 
    gsl_sf_result c;
282
 
    cheb_eval_e(&erfc_x510_cs, t, &c);
283
 
    e_val = exterm * c.val;
284
 
    e_err = exterm * (c.err + 2.0*fabs(x)*GSL_DBL_EPSILON + GSL_DBL_EPSILON);
285
 
  }
286
 
  else {
287
 
    e_val = erfc8(ax);
288
 
    e_err = (x*x + 1.0) * GSL_DBL_EPSILON * fabs(e_val);
289
 
  }
290
 
 
291
 
  if(x < 0.0) {
292
 
    result->val  = 2.0 - e_val;
293
 
    result->err  = e_err;
294
 
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
295
 
  }
296
 
  else {
297
 
    result->val  = e_val;
298
 
    result->err  = e_err;
299
 
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
300
 
  }
301
 
 
302
 
  return GSL_SUCCESS;
303
 
}
304
 
 
305
 
 
306
 
int gsl_sf_log_erfc_e(double x, gsl_sf_result * result)
307
 
{
308
 
  /* CHECK_POINTER(result) */
309
 
 
310
 
  if(x*x < 10.0*GSL_ROOT6_DBL_EPSILON) {
311
 
    const double y = x / M_SQRTPI;
312
 
    /* series for -1/2 Log[Erfc[Sqrt[Pi] y]] */
313
 
    const double c3 = (4.0 - M_PI)/3.0;
314
 
    const double c4 = 2.0*(1.0 - M_PI/3.0);
315
 
    const double c5 = -0.001829764677455021;  /* (96.0 - 40.0*M_PI + 3.0*M_PI*M_PI)/30.0  */
316
 
    const double c6 =  0.02629651521057465;   /* 2.0*(120.0 - 60.0*M_PI + 7.0*M_PI*M_PI)/45.0 */
317
 
    const double c7 = -0.01621575378835404;
318
 
    const double c8 =  0.00125993961762116;
319
 
    const double c9 =  0.00556964649138;
320
 
    const double c10 = -0.0045563339802;
321
 
    const double c11 =  0.0009461589032;
322
 
    const double c12 =  0.0013200243174;
323
 
    const double c13 = -0.00142906;
324
 
    const double c14 =  0.00048204;
325
 
    double series = c8 + y*(c9 + y*(c10 + y*(c11 + y*(c12 + y*(c13 + c14*y)))));
326
 
    series = y*(1.0 + y*(1.0 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*(c7 + y*series)))))));
327
 
    result->val = -2.0 * series;
328
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
329
 
    return GSL_SUCCESS;
330
 
  }
331
 
  /*
332
 
  don't like use of log1p(); added above series stuff for small x instead, should be ok [GJ]
333
 
  else if (fabs(x) < 1.0) {
334
 
    gsl_sf_result result_erf;
335
 
    gsl_sf_erf_e(x, &result_erf);
336
 
    result->val  = log1p(-result_erf.val);
337
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
338
 
    return GSL_SUCCESS;
339
 
  }
340
 
  */
341
 
  else if(x > 8.0) {
342
 
    result->val = log_erfc8(x);
343
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
344
 
    return GSL_SUCCESS;
345
 
  }
346
 
  else {
347
 
    gsl_sf_result result_erfc;
348
 
    gsl_sf_erfc_e(x, &result_erfc);
349
 
    result->val  = log(result_erfc.val);
350
 
    result->err  = fabs(result_erfc.err / result_erfc.val);
351
 
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
352
 
    return GSL_SUCCESS;
353
 
  }
354
 
}
355
 
 
356
 
 
357
 
int gsl_sf_erf_e(double x, gsl_sf_result * result)
358
 
{
359
 
  /* CHECK_POINTER(result) */
360
 
 
361
 
  if(fabs(x) < 1.0) {
362
 
    return erfseries(x, result);
363
 
  }
364
 
  else {
365
 
    gsl_sf_result result_erfc;
366
 
    gsl_sf_erfc_e(x, &result_erfc);
367
 
    result->val  = 1.0 - result_erfc.val;
368
 
    result->err  = result_erfc.err;
369
 
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
370
 
    return GSL_SUCCESS;
371
 
  }
372
 
}
373
 
 
374
 
 
375
 
int gsl_sf_erf_Z_e(double x, gsl_sf_result * result)
376
 
{
377
 
  /* CHECK_POINTER(result) */
378
 
 
379
 
  {
380
 
    const double ex2 = exp(-x*x/2.0);
381
 
    result->val  = ex2 / (M_SQRT2 * M_SQRTPI);
382
 
    result->err  = fabs(x * result->val) * GSL_DBL_EPSILON;
383
 
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
384
 
    CHECK_UNDERFLOW(result);
385
 
    return GSL_SUCCESS;
386
 
  }
387
 
}
388
 
 
389
 
 
390
 
int gsl_sf_erf_Q_e(double x, gsl_sf_result * result)
391
 
{
392
 
  /* CHECK_POINTER(result) */
393
 
 
394
 
  {
395
 
    gsl_sf_result result_erfc;
396
 
    int stat = gsl_sf_erfc_e(x/M_SQRT2, &result_erfc);
397
 
    result->val  = 0.5 * result_erfc.val;
398
 
    result->err  = 0.5 * result_erfc.err;
399
 
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
400
 
    return stat;
401
 
  }
402
 
}
403
 
 
404
 
 
405
 
int gsl_sf_hazard_e(double x, gsl_sf_result * result)
406
 
{
407
 
  if(x < 25.0)
408
 
  {
409
 
    gsl_sf_result result_ln_erfc;
410
 
    const int stat_l = gsl_sf_log_erfc_e(x/M_SQRT2, &result_ln_erfc);
411
 
    const double lnc = -0.22579135264472743236; /* ln(sqrt(2/pi)) */
412
 
    const double arg = lnc - 0.5*x*x - result_ln_erfc.val;
413
 
    const int stat_e = gsl_sf_exp_e(arg, result);
414
 
    result->err += 3.0 * (1.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
415
 
    result->err += fabs(result_ln_erfc.err * result->val);
416
 
    return GSL_ERROR_SELECT_2(stat_l, stat_e);
417
 
  }
418
 
  else
419
 
  {
420
 
    const double ix2 = 1.0/(x*x);
421
 
    const double corrB = 1.0 - 9.0*ix2 * (1.0 - 11.0*ix2);
422
 
    const double corrM = 1.0 - 5.0*ix2 * (1.0 - 7.0*ix2 * corrB);
423
 
    const double corrT = 1.0 - ix2 * (1.0 - 3.0*ix2*corrM);
424
 
    result->val = x / corrT;
425
 
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
426
 
    return GSL_SUCCESS;
427
 
  }
428
 
}
429
 
 
430
 
 
431
 
 
432
 
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
433
 
 
434
 
#include "eval.h"
435
 
 
436
 
double gsl_sf_erfc(double x)
437
 
{
438
 
  EVAL_RESULT(gsl_sf_erfc_e(x, &result));
439
 
}
440
 
 
441
 
double gsl_sf_log_erfc(double x)
442
 
{
443
 
  EVAL_RESULT(gsl_sf_log_erfc_e(x, &result));
444
 
}
445
 
 
446
 
double gsl_sf_erf(double x)
447
 
{
448
 
  EVAL_RESULT(gsl_sf_erf_e(x, &result));
449
 
}
450
 
 
451
 
double gsl_sf_erf_Z(double x)
452
 
{
453
 
  EVAL_RESULT(gsl_sf_erf_Z_e(x, &result));
454
 
}
455
 
 
456
 
double gsl_sf_erf_Q(double x)
457
 
{
458
 
  EVAL_RESULT(gsl_sf_erf_Q_e(x, &result));
459
 
}
460
 
 
461
 
double gsl_sf_hazard(double x)
462
 
{
463
 
  EVAL_RESULT(gsl_sf_hazard_e(x, &result));
464
 
}
465