1 /* sf_erf.c -- float version of s_erf.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3 */
4
5 /*
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 *
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
12 * is preserved.
13 * ====================================================
14 */
15
16 #include "fdlibm.h"
17 #include "math_config.h"
18
19 #ifdef __v810__
20 #define const
21 #endif
22
23 static const float tiny = 1e-30, half = 5.0000000000e-01, /* 0x3F000000 */
24 one = 1.0000000000e+00, /* 0x3F800000 */
25 two = 2.0000000000e+00, /* 0x40000000 */
26 /* c = (subfloat)0.84506291151 */
27 erx = 8.4506291151e-01, /* 0x3f58560b */
28 /*
29 * Coefficients for approximation to erf on [0,0.84375]
30 */
31 efx = 1.2837916613e-01, /* 0x3e0375d4 */
32 efx8 = 1.0270333290e+00, /* 0x3f8375d4 */
33 pp0 = 1.2837916613e-01, /* 0x3e0375d4 */
34 pp1 = -3.2504209876e-01, /* 0xbea66beb */
35 pp2 = -2.8481749818e-02, /* 0xbce9528f */
36 pp3 = -5.7702702470e-03, /* 0xbbbd1489 */
37 pp4 = -2.3763017452e-05, /* 0xb7c756b1 */
38 qq1 = 3.9791721106e-01, /* 0x3ecbbbce */
39 qq2 = 6.5022252500e-02, /* 0x3d852a63 */
40 qq3 = 5.0813062117e-03, /* 0x3ba68116 */
41 qq4 = 1.3249473704e-04, /* 0x390aee49 */
42 qq5 = -3.9602282413e-06, /* 0xb684e21a */
43 /*
44 * Coefficients for approximation to erf in [0.84375,1.25]
45 */
46 pa0 = -2.3621185683e-03, /* 0xbb1acdc6 */
47 pa1 = 4.1485610604e-01, /* 0x3ed46805 */
48 pa2 = -3.7220788002e-01, /* 0xbebe9208 */
49 pa3 = 3.1834661961e-01, /* 0x3ea2fe54 */
50 pa4 = -1.1089469492e-01, /* 0xbde31cc2 */
51 pa5 = 3.5478305072e-02, /* 0x3d1151b3 */
52 pa6 = -2.1663755178e-03, /* 0xbb0df9c0 */
53 qa1 = 1.0642088205e-01, /* 0x3dd9f331 */
54 qa2 = 5.4039794207e-01, /* 0x3f0a5785 */
55 qa3 = 7.1828655899e-02, /* 0x3d931ae7 */
56 qa4 = 1.2617121637e-01, /* 0x3e013307 */
57 qa5 = 1.3637083583e-02, /* 0x3c5f6e13 */
58 qa6 = 1.1984500103e-02, /* 0x3c445aa3 */
59 /*
60 * Coefficients for approximation to erfc in [1.25,1/0.35]
61 */
62 ra0 = -9.8649440333e-03, /* 0xbc21a093 */
63 ra1 = -6.9385856390e-01, /* 0xbf31a0b7 */
64 ra2 = -1.0558626175e+01, /* 0xc128f022 */
65 ra3 = -6.2375331879e+01, /* 0xc2798057 */
66 ra4 = -1.6239666748e+02, /* 0xc322658c */
67 ra5 = -1.8460508728e+02, /* 0xc3389ae7 */
68 ra6 = -8.1287437439e+01, /* 0xc2a2932b */
69 ra7 = -9.8143291473e+00, /* 0xc11d077e */
70 sa1 = 1.9651271820e+01, /* 0x419d35ce */
71 sa2 = 1.3765776062e+02, /* 0x4309a863 */
72 sa3 = 4.3456588745e+02, /* 0x43d9486f */
73 sa4 = 6.4538726807e+02, /* 0x442158c9 */
74 sa5 = 4.2900814819e+02, /* 0x43d6810b */
75 sa6 = 1.0863500214e+02, /* 0x42d9451f */
76 sa7 = 6.5702495575e+00, /* 0x40d23f7c */
77 sa8 = -6.0424413532e-02, /* 0xbd777f97 */
78 /*
79 * Coefficients for approximation to erfc in [1/.35,28]
80 */
81 rb0 = -9.8649431020e-03, /* 0xbc21a092 */
82 rb1 = -7.9928326607e-01, /* 0xbf4c9dd4 */
83 rb2 = -1.7757955551e+01, /* 0xc18e104b */
84 rb3 = -1.6063638306e+02, /* 0xc320a2ea */
85 rb4 = -6.3756646729e+02, /* 0xc41f6441 */
86 rb5 = -1.0250950928e+03, /* 0xc480230b */
87 rb6 = -4.8351919556e+02, /* 0xc3f1c275 */
88 sb1 = 3.0338060379e+01, /* 0x41f2b459 */
89 sb2 = 3.2579251099e+02, /* 0x43a2e571 */
90 sb3 = 1.5367296143e+03, /* 0x44c01759 */
91 sb4 = 3.1998581543e+03, /* 0x4547fdbb */
92 sb5 = 2.5530502930e+03, /* 0x451f90ce */
93 sb6 = 4.7452853394e+02, /* 0x43ed43a7 */
94 sb7 = -2.2440952301e+01; /* 0xc1b38712 */
95
96 float
erff(float x)97 erff(float x)
98 {
99 __int32_t hx, ix, i;
100 float R, S, P, Q, s, y, z, r;
101 GET_FLOAT_WORD(hx, x);
102 ix = hx & 0x7fffffff;
103 if (!FLT_UWORD_IS_FINITE(ix)) { /* erf(nan)=nan */
104 i = ((__uint32_t)hx >> 31) << 1;
105 return (float)(1 - i) + one / x; /* erf(+-inf)=+-1 */
106 }
107
108 if (ix < 0x3f580000) { /* |x|<0.84375 */
109 if (ix < 0x31800000) { /* |x|<2**-28 */
110 if (ix < 0x04000000)
111 /*avoid underflow */
112 return (float)0.125 * ((float)8.0 * x + efx8 * x);
113 return x + efx * x;
114 }
115 z = x * x;
116 r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
117 s = one + z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z * qq5))));
118 y = r / s;
119 return x + x * y;
120 }
121 if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
122 s = fabsf(x) - one;
123 P = pa0 +
124 s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 + s * pa6)))));
125 Q = one +
126 s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 + s * qa6)))));
127 if (hx >= 0)
128 return erx + P / Q;
129 else
130 return -erx - P / Q;
131 }
132 if (ix >= 0x40c00000) { /* inf>|x|>=6 */
133 if (hx >= 0)
134 return one - tiny;
135 else
136 return tiny - one;
137 }
138 x = fabsf(x);
139 s = one / (x * x);
140 if (ix < 0x4036DB6E) { /* |x| < 1/0.35 */
141 R = ra0 +
142 s * (ra1 +
143 s * (ra2 +
144 s * (ra3 + s * (ra4 + s * (ra5 + s * (ra6 + s * ra7))))));
145 S = one +
146 s * (sa1 +
147 s * (sa2 +
148 s * (sa3 +
149 s * (sa4 +
150 s * (sa5 + s * (sa6 + s * (sa7 + s * sa8)))))));
151 } else { /* |x| >= 1/0.35 */
152 R = rb0 +
153 s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 + s * (rb5 + s * rb6)))));
154 S = one +
155 s * (sb1 +
156 s * (sb2 +
157 s * (sb3 + s * (sb4 + s * (sb5 + s * (sb6 + s * sb7))))));
158 }
159 GET_FLOAT_WORD(ix, x);
160 SET_FLOAT_WORD(z, ix & 0xfffff000);
161 r = expf(-z * z - (float)0.5625) * expf((z - x) * (z + x) + R / S);
162 if (hx >= 0)
163 return one - r / x;
164 else
165 return r / x - one;
166 }
167
168 float
erfcf(float x)169 erfcf(float x)
170 {
171 __int32_t hx, ix;
172 float R, S, P, Q, s, y, z, r;
173 GET_FLOAT_WORD(hx, x);
174 ix = hx & 0x7fffffff;
175 if (!FLT_UWORD_IS_FINITE(ix)) { /* erfc(nan)=nan */
176 /* erfc(+-inf)=0,2 */
177 return (float)(((__uint32_t)hx >> 31) << 1) + one / x;
178 }
179
180 if (ix < 0x3f580000) { /* |x|<0.84375 */
181 if (ix < 0x23800000) /* |x|<2**-56 */
182 return one - x;
183 z = x * x;
184 r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
185 s = one + z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z * qq5))));
186 y = r / s;
187 if (hx < 0x3e800000) { /* x<1/4 */
188 return one - (x + x * y);
189 } else {
190 r = x * y;
191 r += (x - half);
192 return half - r;
193 }
194 }
195 if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
196 s = fabsf(x) - one;
197 P = pa0 +
198 s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 + s * pa6)))));
199 Q = one +
200 s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 + s * qa6)))));
201 if (hx >= 0) {
202 z = one - erx;
203 return z - P / Q;
204 } else {
205 z = erx + P / Q;
206 return one + z;
207 }
208 }
209 if (ix < 0x41e00000) { /* |x|<28 */
210 x = fabsf(x);
211 s = one / (x * x);
212 if (ix < 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
213 R = ra0 +
214 s * (ra1 +
215 s * (ra2 +
216 s * (ra3 +
217 s * (ra4 + s * (ra5 + s * (ra6 + s * ra7))))));
218 S = one +
219 s * (sa1 +
220 s * (sa2 +
221 s * (sa3 +
222 s * (sa4 +
223 s * (sa5 +
224 s * (sa6 + s * (sa7 + s * sa8)))))));
225 } else { /* |x| >= 1/.35 ~ 2.857143 */
226 if (hx < 0 && ix >= 0x40c00000)
227 return two - tiny; /* x < -6 */
228 R = rb0 +
229 s * (rb1 +
230 s * (rb2 + s * (rb3 + s * (rb4 + s * (rb5 + s * rb6)))));
231 S = one +
232 s * (sb1 +
233 s * (sb2 +
234 s * (sb3 +
235 s * (sb4 + s * (sb5 + s * (sb6 + s * sb7))))));
236 }
237 GET_FLOAT_WORD(ix, x);
238 SET_FLOAT_WORD(z, ix & 0xfffff000);
239 r = expf(-z * z - (float)0.5625) * expf((z - x) * (z + x) + R / S);
240 if (hx > 0)
241 return r / x;
242 else
243 return two - r / x;
244 } else {
245 if (hx > 0)
246 return __math_uflowf(0);
247 else
248 return two - tiny;
249 }
250 }
251
252 _MATH_ALIAS_f_f(erf)
253
254 _MATH_ALIAS_f_f(erfc)
255
