10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #pragma weak erf = __erf
31 #pragma weak erfc = __erfc
32
33 /* INDENT OFF */
34 /*
35 * double erf(double x)
36 * double erfc(double x)
37 * x
38 * 2 |\
39 * erf(x) = --------- | exp(-t*t)dt
40 * sqrt(pi) \|
41 * 0
42 *
43 * erfc(x) = 1-erf(x)
44 * Note that
45 * erf(-x) = -erf(x)
46 * erfc(-x) = 2 - erfc(x)
47 *
48 * Method:
49 * 1. For |x| in [0, 0.84375]
50 * erf(x) = x + x*R(x^2)
51 * erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
110 * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
111 * x*sqrt(pi)
112 * We use rational approximation to approximate
113 * g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
114 * Here is the error bound for R1/S1 and R2/S2
115 * |R1/S1 - f(x)| < 2**(-62.57)
116 * |R2/S2 - f(x)| < 2**(-61.52)
117 *
118 * 5. For inf > x >= 28
119 * erf(x) = sign(x) *(1 - tiny) (raise inexact)
120 * erfc(x) = tiny*tiny (raise underflow) if x > 0
121 * = 2 - tiny if x<0
122 *
123 * 7. Special case:
124 * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
125 * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
126 * erfc/erf(NaN) is NaN
127 */
128 /* INDENT ON */
129
130 #include "libm_synonyms.h" /* __erf, __erfc, __exp */
131 #include "libm_macros.h"
132 #include <math.h>
133
134 static const double xxx[] = {
135 /* tiny */ 1e-300,
136 /* half */ 5.00000000000000000000e-01, /* 3FE00000, 00000000 */
137 /* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */
138 /* two */ 2.00000000000000000000e+00, /* 40000000, 00000000 */
139 /* erx */ 8.45062911510467529297e-01, /* 3FEB0AC1, 60000000 */
140 /*
141 * Coefficients for approximation to erf on [0,0.84375]
142 */
143 /* efx */ 1.28379167095512586316e-01, /* 3FC06EBA, 8214DB69 */
144 /* efx8 */ 1.02703333676410069053e+00, /* 3FF06EBA, 8214DB69 */
145 /* pp0 */ 1.28379167095512558561e-01, /* 3FC06EBA, 8214DB68 */
146 /* pp1 */ -3.25042107247001499370e-01, /* BFD4CD7D, 691CB913 */
147 /* pp2 */ -2.84817495755985104766e-02, /* BF9D2A51, DBD7194F */
148 /* pp3 */ -5.77027029648944159157e-03, /* BF77A291, 236668E4 */
149 /* pp4 */ -2.37630166566501626084e-05, /* BEF8EAD6, 120016AC */
150 /* qq1 */ 3.97917223959155352819e-01, /* 3FD97779, CDDADC09 */
|
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #pragma weak __erf = erf
31 #pragma weak __erfc = erfc
32
33 /* INDENT OFF */
34 /*
35 * double erf(double x)
36 * double erfc(double x)
37 * x
38 * 2 |\
39 * erf(x) = --------- | exp(-t*t)dt
40 * sqrt(pi) \|
41 * 0
42 *
43 * erfc(x) = 1-erf(x)
44 * Note that
45 * erf(-x) = -erf(x)
46 * erfc(-x) = 2 - erfc(x)
47 *
48 * Method:
49 * 1. For |x| in [0, 0.84375]
50 * erf(x) = x + x*R(x^2)
51 * erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
110 * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
111 * x*sqrt(pi)
112 * We use rational approximation to approximate
113 * g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
114 * Here is the error bound for R1/S1 and R2/S2
115 * |R1/S1 - f(x)| < 2**(-62.57)
116 * |R2/S2 - f(x)| < 2**(-61.52)
117 *
118 * 5. For inf > x >= 28
119 * erf(x) = sign(x) *(1 - tiny) (raise inexact)
120 * erfc(x) = tiny*tiny (raise underflow) if x > 0
121 * = 2 - tiny if x<0
122 *
123 * 7. Special case:
124 * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
125 * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
126 * erfc/erf(NaN) is NaN
127 */
128 /* INDENT ON */
129
130 #include "libm_macros.h"
131 #include <math.h>
132
133 static const double xxx[] = {
134 /* tiny */ 1e-300,
135 /* half */ 5.00000000000000000000e-01, /* 3FE00000, 00000000 */
136 /* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */
137 /* two */ 2.00000000000000000000e+00, /* 40000000, 00000000 */
138 /* erx */ 8.45062911510467529297e-01, /* 3FEB0AC1, 60000000 */
139 /*
140 * Coefficients for approximation to erf on [0,0.84375]
141 */
142 /* efx */ 1.28379167095512586316e-01, /* 3FC06EBA, 8214DB69 */
143 /* efx8 */ 1.02703333676410069053e+00, /* 3FF06EBA, 8214DB69 */
144 /* pp0 */ 1.28379167095512558561e-01, /* 3FC06EBA, 8214DB68 */
145 /* pp1 */ -3.25042107247001499370e-01, /* BFD4CD7D, 691CB913 */
146 /* pp2 */ -2.84817495755985104766e-02, /* BF9D2A51, DBD7194F */
147 /* pp3 */ -5.77027029648944159157e-03, /* BF77A291, 236668E4 */
148 /* pp4 */ -2.37630166566501626084e-05, /* BEF8EAD6, 120016AC */
149 /* qq1 */ 3.97917223959155352819e-01, /* 3FD97779, CDDADC09 */
|