source: trunk/libs/newlib/src/newlib/libm/math/sf_erf.c @ 543

Last change on this file since 543 was 444, checked in by satin@…, 6 years ago

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