GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: lib/libm/src/ld80/s_erfl.c Lines: 0 107 0.0 %
Date: 2017-11-13 Branches: 0 44 0.0 %

Line Branch Exec Source
1
/*
2
 * ====================================================
3
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4
 *
5
 * Developed at SunPro, a Sun Microsystems, Inc. business.
6
 * Permission to use, copy, modify, and distribute this
7
 * software is freely granted, provided that this notice
8
 * is preserved.
9
 * ====================================================
10
 */
11
12
/*
13
 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
14
 *
15
 * Permission to use, copy, modify, and distribute this software for any
16
 * purpose with or without fee is hereby granted, provided that the above
17
 * copyright notice and this permission notice appear in all copies.
18
 *
19
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
20
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
21
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
22
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
23
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
24
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
25
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
26
 */
27
28
/* double erf(double x)
29
 * double erfc(double x)
30
 *			     x
31
 *		      2      |\
32
 *     erf(x)  =  ---------  | exp(-t*t)dt
33
 *		   sqrt(pi) \|
34
 *			     0
35
 *
36
 *     erfc(x) =  1-erf(x)
37
 *  Note that
38
 *		erf(-x) = -erf(x)
39
 *		erfc(-x) = 2 - erfc(x)
40
 *
41
 * Method:
42
 *	1. For |x| in [0, 0.84375]
43
 *	    erf(x)  = x + x*R(x^2)
44
 *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
45
 *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
46
 *	   Remark. The formula is derived by noting
47
 *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
48
 *	   and that
49
 *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
50
 *	   is close to one. The interval is chosen because the fix
51
 *	   point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
52
 *	   near 0.6174), and by some experiment, 0.84375 is chosen to
53
 *	   guarantee the error is less than one ulp for erf.
54
 *
55
 *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
56
 *         c = 0.84506291151 rounded to single (24 bits)
57
 *	erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
58
 *	erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
59
 *			  1+(c+P1(s)/Q1(s))    if x < 0
60
 *	   Remark: here we use the taylor series expansion at x=1.
61
 *		erf(1+s) = erf(1) + s*Poly(s)
62
 *			 = 0.845.. + P1(s)/Q1(s)
63
 *	   Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
64
 *
65
 *      3. For x in [1.25,1/0.35(~2.857143)],
66
 *	erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
67
 *              z=1/x^2
68
 *	erf(x)  = 1 - erfc(x)
69
 *
70
 *      4. For x in [1/0.35,107]
71
 *	erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
72
 *			= 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z))
73
 *                             if -6.666<x<0
74
 *			= 2.0 - tiny		(if x <= -6.666)
75
 *              z=1/x^2
76
 *	erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
77
 *	erf(x)  = sign(x)*(1.0 - tiny)
78
 *      Note1:
79
 *	   To compute exp(-x*x-0.5625+R/S), let s be a single
80
 *	   precision number and s := x; then
81
 *		-x*x = -s*s + (s-x)*(s+x)
82
 *	        exp(-x*x-0.5626+R/S) =
83
 *			exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
84
 *      Note2:
85
 *	   Here 4 and 5 make use of the asymptotic series
86
 *			  exp(-x*x)
87
 *		erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
88
 *			  x*sqrt(pi)
89
 *
90
 *      5. For inf > x >= 107
91
 *	erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
92
 *	erfc(x) = tiny*tiny (raise underflow) if x > 0
93
 *			= 2 - tiny if x<0
94
 *
95
 *      7. Special case:
96
 *	erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
97
 *	erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
98
 *		erfc/erf(NaN) is NaN
99
 */
100
101
102
#include <math.h>
103
104
#include "math_private.h"
105
106
static const long double
107
tiny = 1e-4931L,
108
  half = 0.5L,
109
  one = 1.0L,
110
  two = 2.0L,
111
	/* c = (float)0.84506291151 */
112
  erx = 0.845062911510467529296875L,
113
/*
114
 * Coefficients for approximation to  erf on [0,0.84375]
115
 */
116
  /* 2/sqrt(pi) - 1 */
117
  efx = 1.2837916709551257389615890312154517168810E-1L,
118
  /* 8 * (2/sqrt(pi) - 1) */
119
  efx8 = 1.0270333367641005911692712249723613735048E0L,
120
121
  pp[6] = {
122
    1.122751350964552113068262337278335028553E6L,
123
    -2.808533301997696164408397079650699163276E6L,
124
    -3.314325479115357458197119660818768924100E5L,
125
    -6.848684465326256109712135497895525446398E4L,
126
    -2.657817695110739185591505062971929859314E3L,
127
    -1.655310302737837556654146291646499062882E2L,
128
  },
129
130
  qq[6] = {
131
    8.745588372054466262548908189000448124232E6L,
132
    3.746038264792471129367533128637019611485E6L,
133
    7.066358783162407559861156173539693900031E5L,
134
    7.448928604824620999413120955705448117056E4L,
135
    4.511583986730994111992253980546131408924E3L,
136
    1.368902937933296323345610240009071254014E2L,
137
    /* 1.000000000000000000000000000000000000000E0 */
138
  },
139
140
/*
141
 * Coefficients for approximation to  erf  in [0.84375,1.25]
142
 */
143
/* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
144
   -0.15625 <= x <= +.25
145
   Peak relative error 8.5e-22  */
146
147
  pa[8] = {
148
    -1.076952146179812072156734957705102256059E0L,
149
     1.884814957770385593365179835059971587220E2L,
150
    -5.339153975012804282890066622962070115606E1L,
151
     4.435910679869176625928504532109635632618E1L,
152
     1.683219516032328828278557309642929135179E1L,
153
    -2.360236618396952560064259585299045804293E0L,
154
     1.852230047861891953244413872297940938041E0L,
155
     9.394994446747752308256773044667843200719E-2L,
156
  },
157
158
  qa[7] =  {
159
    4.559263722294508998149925774781887811255E2L,
160
    3.289248982200800575749795055149780689738E2L,
161
    2.846070965875643009598627918383314457912E2L,
162
    1.398715859064535039433275722017479994465E2L,
163
    6.060190733759793706299079050985358190726E1L,
164
    2.078695677795422351040502569964299664233E1L,
165
    4.641271134150895940966798357442234498546E0L,
166
    /* 1.000000000000000000000000000000000000000E0 */
167
  },
168
169
/*
170
 * Coefficients for approximation to  erfc in [1.25,1/0.35]
171
 */
172
/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
173
   1/2.85711669921875 < 1/x < 1/1.25
174
   Peak relative error 3.1e-21  */
175
176
    ra[] = {
177
      1.363566591833846324191000679620738857234E-1L,
178
      1.018203167219873573808450274314658434507E1L,
179
      1.862359362334248675526472871224778045594E2L,
180
      1.411622588180721285284945138667933330348E3L,
181
      5.088538459741511988784440103218342840478E3L,
182
      8.928251553922176506858267311750789273656E3L,
183
      7.264436000148052545243018622742770549982E3L,
184
      2.387492459664548651671894725748959751119E3L,
185
      2.220916652813908085449221282808458466556E2L,
186
    },
187
188
    sa[] = {
189
      -1.382234625202480685182526402169222331847E1L,
190
      -3.315638835627950255832519203687435946482E2L,
191
      -2.949124863912936259747237164260785326692E3L,
192
      -1.246622099070875940506391433635999693661E4L,
193
      -2.673079795851665428695842853070996219632E4L,
194
      -2.880269786660559337358397106518918220991E4L,
195
      -1.450600228493968044773354186390390823713E4L,
196
      -2.874539731125893533960680525192064277816E3L,
197
      -1.402241261419067750237395034116942296027E2L,
198
      /* 1.000000000000000000000000000000000000000E0 */
199
    },
200
/*
201
 * Coefficients for approximation to  erfc in [1/.35,107]
202
 */
203
/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
204
   1/6.6666259765625 < 1/x < 1/2.85711669921875
205
   Peak relative error 4.2e-22  */
206
    rb[] = {
207
      -4.869587348270494309550558460786501252369E-5L,
208
      -4.030199390527997378549161722412466959403E-3L,
209
      -9.434425866377037610206443566288917589122E-2L,
210
      -9.319032754357658601200655161585539404155E-1L,
211
      -4.273788174307459947350256581445442062291E0L,
212
      -8.842289940696150508373541814064198259278E0L,
213
      -7.069215249419887403187988144752613025255E0L,
214
      -1.401228723639514787920274427443330704764E0L,
215
    },
216
217
    sb[] = {
218
      4.936254964107175160157544545879293019085E-3L,
219
      1.583457624037795744377163924895349412015E-1L,
220
      1.850647991850328356622940552450636420484E0L,
221
      9.927611557279019463768050710008450625415E0L,
222
      2.531667257649436709617165336779212114570E1L,
223
      2.869752886406743386458304052862814690045E1L,
224
      1.182059497870819562441683560749192539345E1L,
225
      /* 1.000000000000000000000000000000000000000E0 */
226
    },
227
/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
228
   1/107 <= 1/x <= 1/6.6666259765625
229
   Peak relative error 1.1e-21  */
230
    rc[] = {
231
      -8.299617545269701963973537248996670806850E-5L,
232
      -6.243845685115818513578933902532056244108E-3L,
233
      -1.141667210620380223113693474478394397230E-1L,
234
      -7.521343797212024245375240432734425789409E-1L,
235
      -1.765321928311155824664963633786967602934E0L,
236
      -1.029403473103215800456761180695263439188E0L,
237
    },
238
239
    sc[] = {
240
      8.413244363014929493035952542677768808601E-3L,
241
      2.065114333816877479753334599639158060979E-1L,
242
      1.639064941530797583766364412782135680148E0L,
243
      4.936788463787115555582319302981666347450E0L,
244
      5.005177727208955487404729933261347679090E0L,
245
      /* 1.000000000000000000000000000000000000000E0 */
246
    };
247
248
long double
249
erfl(long double x)
250
{
251
  long double R, S, P, Q, s, y, z, r;
252
  int32_t ix, i;
253
  u_int32_t se, i0, i1;
254
255
  GET_LDOUBLE_WORDS (se, i0, i1, x);
256
  ix = se & 0x7fff;
257
258
  if (ix >= 0x7fff)
259
    {				/* erf(nan)=nan */
260
      i = ((se & 0xffff) >> 15) << 1;
261
      return (long double) (1 - i) + one / x;	/* erf(+-inf)=+-1 */
262
    }
263
264
  ix = (ix << 16) | (i0 >> 16);
265
  if (ix < 0x3ffed800) /* |x|<0.84375 */
266
    {
267
      if (ix < 0x3fde8000) /* |x|<2**-33 */
268
	{
269
	  if (ix < 0x00080000)
270
	    return 0.125 * (8.0 * x + efx8 * x);	/*avoid underflow */
271
	  return x + efx * x;
272
	}
273
      z = x * x;
274
      r = pp[0] + z * (pp[1]
275
	+ z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
276
      s = qq[0] + z * (qq[1]
277
	+ z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
278
      y = r / s;
279
      return x + x * y;
280
    }
281
  if (ix < 0x3fffa000) /* 1.25 */
282
    {				/* 0.84375 <= |x| < 1.25 */
283
      s = fabsl (x) - one;
284
      P = pa[0] + s * (pa[1] + s * (pa[2]
285
	+ s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
286
      Q = qa[0] + s * (qa[1] + s * (qa[2]
287
	+ s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
288
      if ((se & 0x8000) == 0)
289
	return erx + P / Q;
290
      else
291
	return -erx - P / Q;
292
    }
293
  if (ix >= 0x4001d555) /* 6.6666259765625 */
294
    {				/* inf>|x|>=6.666 */
295
      if ((se & 0x8000) == 0)
296
	return one - tiny;
297
      else
298
	return tiny - one;
299
    }
300
  x = fabsl (x);
301
  s = one / (x * x);
302
  if (ix < 0x4000b6db) /* 2.85711669921875 */
303
    {
304
      R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
305
	s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
306
      S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
307
	s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
308
    }
309
  else
310
    {				/* |x| >= 1/0.35 */
311
      R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
312
	s * (rb[5] + s * (rb[6] + s * rb[7]))))));
313
      S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
314
	s * (sb[5] + s * (sb[6] + s))))));
315
    }
316
  z = x;
317
  GET_LDOUBLE_WORDS (i, i0, i1, z);
318
  i1 = 0;
319
  SET_LDOUBLE_WORDS (z, i, i0, i1);
320
  r =
321
    expl (-z * z - 0.5625) * expl ((z - x) * (z + x) + R / S);
322
  if ((se & 0x8000) == 0)
323
    return one - r / x;
324
  else
325
    return r / x - one;
326
}
327
DEF_STD(erfl);
328
329
long double
330
erfcl(long double x)
331
{
332
  int32_t hx, ix;
333
  long double R, S, P, Q, s, y, z, r;
334
  u_int32_t se, i0, i1;
335
336
  GET_LDOUBLE_WORDS (se, i0, i1, x);
337
  ix = se & 0x7fff;
338
  if (ix >= 0x7fff)
339
    {				/* erfc(nan)=nan */
340
      /* erfc(+-inf)=0,2 */
341
      return (long double) (((se & 0xffff) >> 15) << 1) + one / x;
342
    }
343
344
  ix = (ix << 16) | (i0 >> 16);
345
  if (ix < 0x3ffed800) /* |x|<0.84375 */
346
    {
347
      if (ix < 0x3fbe0000) /* |x|<2**-65 */
348
	return one - x;
349
      z = x * x;
350
      r = pp[0] + z * (pp[1]
351
	+ z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
352
      s = qq[0] + z * (qq[1]
353
	+ z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
354
      y = r / s;
355
      if (ix < 0x3ffd8000) /* x<1/4 */
356
	{
357
	  return one - (x + x * y);
358
	}
359
      else
360
	{
361
	  r = x * y;
362
	  r += (x - half);
363
	  return half - r;
364
	}
365
    }
366
  if (ix < 0x3fffa000) /* 1.25 */
367
    {				/* 0.84375 <= |x| < 1.25 */
368
      s = fabsl (x) - one;
369
      P = pa[0] + s * (pa[1] + s * (pa[2]
370
	+ s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
371
      Q = qa[0] + s * (qa[1] + s * (qa[2]
372
	+ s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
373
      if ((se & 0x8000) == 0)
374
	{
375
	  z = one - erx;
376
	  return z - P / Q;
377
	}
378
      else
379
	{
380
	  z = erx + P / Q;
381
	  return one + z;
382
	}
383
    }
384
  if (ix < 0x4005d600) /* 107 */
385
    {				/* |x|<107 */
386
      x = fabsl (x);
387
      s = one / (x * x);
388
      if (ix < 0x4000b6db) /* 2.85711669921875 */
389
	{			/* |x| < 1/.35 ~ 2.857143 */
390
	  R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
391
	    s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
392
	  S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
393
	    s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
394
	}
395
      else if (ix < 0x4001d555) /* 6.6666259765625 */
396
	{			/* 6.666 > |x| >= 1/.35 ~ 2.857143 */
397
	  R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
398
	    s * (rb[5] + s * (rb[6] + s * rb[7]))))));
399
	  S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
400
	    s * (sb[5] + s * (sb[6] + s))))));
401
	}
402
      else
403
	{			/* |x| >= 6.666 */
404
	  if (se & 0x8000)
405
	    return two - tiny;	/* x < -6.666 */
406
407
	  R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
408
						    s * (rc[4] + s * rc[5]))));
409
	  S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
410
						    s * (sc[4] + s))));
411
	}
412
      z = x;
413
      GET_LDOUBLE_WORDS (hx, i0, i1, z);
414
      i1 = 0;
415
      i0 &= 0xffffff00;
416
      SET_LDOUBLE_WORDS (z, hx, i0, i1);
417
      r = expl (-z * z - 0.5625) *
418
	expl ((z - x) * (z + x) + R / S);
419
      if ((se & 0x8000) == 0)
420
	return r / x;
421
      else
422
	return two - r / x;
423
    }
424
  else
425
    {
426
      if ((se & 0x8000) == 0)
427
	return tiny * tiny;
428
      else
429
	return two - tiny;
430
    }
431
}
432
DEF_STD(erfcl);