Rev 1906 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed
| Rev 1906 | Rev 3362 | ||
|---|---|---|---|
| Line 11... | Line 11... | ||
| 11 | * software is freely granted, provided that this notice |
11 | * software is freely granted, provided that this notice |
| 12 | * is preserved. |
12 | * is preserved. |
| 13 | * ==================================================== |
13 | * ==================================================== |
| 14 | */ |
14 | */ |
| Line 15... | Line -... | ||
| 15 | - | ||
| 16 | /* |
15 | |
| 17 | #include "fdlibm.h" |
- | |
| 18 | */ |
- | |
| 19 | #include |
- | |
| 20 | #include |
- | |
| 21 | #include |
- | |
| 22 | - | ||
| 23 | #define __ieee754_expf expf |
- | |
| 24 | - | ||
| 25 | - | ||
| 26 | - | ||
| 27 | typedef union |
- | |
| 28 | { |
- | |
| 29 | float value; |
- | |
| 30 | uint32_t word; |
- | |
| 31 | } ieee_float_shape_type; |
- | |
| 32 | - | ||
| 33 | /* Get a 32 bit int from a float. */ |
- | |
| 34 | - | ||
| 35 | static inline int |
- | |
| 36 | __get_float_word(float d) |
- | |
| 37 | { |
- | |
| 38 | ieee_float_shape_type u; |
- | |
| 39 | u.value = d; |
- | |
| 40 | return u.word; |
- | |
| 41 | } |
- | |
| 42 | - | ||
| 43 | /* Set a float from a 32 bit int. */ |
- | |
| 44 | - | ||
| 45 | #define SET_FLOAT_WORD(d,i) \ |
- | |
| 46 | do { \ |
- | |
| 47 | ieee_float_shape_type sf_u; \ |
- | |
| 48 | sf_u.word = (i); \ |
- | |
| 49 | (d) = sf_u.value; \ |
- | |
| 50 | } while (0) |
- | |
| 51 | - | ||
| 52 | static inline void __trunc_float_word(float * x) |
- | |
| 53 | { |
- | |
| 54 | ieee_float_shape_type u; |
- | |
| 55 | u.value = * x; |
- | |
| 56 | u.word &= 0xfffff000; |
- | |
| Line 57... | Line 16... | ||
| 57 | } |
16 | #include "fdlibm.h" |
| 58 | 17 | ||
| 59 | #ifdef __v810__ |
18 | #ifdef __v810__ |
| Line 144... | Line 103... | ||
| 144 | #else |
103 | #else |
| 145 | float erff(x) |
104 | float erff(x) |
| 146 | float x; |
105 | float x; |
| 147 | #endif |
106 | #endif |
| 148 | { |
107 | { |
| 149 | int32_t hx,ix,i; |
108 | __int32_t hx,ix,i; |
| 150 | float R,S,P,Q,s,y,z,r; |
109 | float R,S,P,Q,s,y,z,r; |
| 151 | hx = __get_float_word(x); |
110 | GET_FLOAT_WORD(hx,x); |
| 152 | ix = hx&0x7fffffff; |
111 | ix = hx&0x7fffffff; |
| 153 | if(!(ix<0x7f800000L)) { /* erf(nan)=nan */ |
112 | if(!FLT_UWORD_IS_FINITE(ix)) { /* erf(nan)=nan */ |
| 154 | i = ((uint32_t)hx>>31)<<1; |
113 | i = ((__uint32_t)hx>>31)<<1; |
| 155 | return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */ |
114 | return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */ |
| 156 | } |
115 | } |
| Line 157... | Line 116... | ||
| 157 | 116 | ||
| 158 | if(ix < 0x3f580000) { /* |x|<0.84375 */ |
117 | if(ix < 0x3f580000) { /* |x|<0.84375 */ |
| Line 188... | Line 147... | ||
| 188 | R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( |
147 | R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( |
| 189 | rb5+s*rb6))))); |
148 | rb5+s*rb6))))); |
| 190 | S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( |
149 | S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( |
| 191 | sb5+s*(sb6+s*sb7)))))); |
150 | sb5+s*(sb6+s*sb7)))))); |
| 192 | } |
151 | } |
| 193 | - | ||
| 194 | z = x; |
152 | GET_FLOAT_WORD(ix,x); |
| 195 | __trunc_float_word (&z); |
153 | SET_FLOAT_WORD(z,ix&0xfffff000); |
| 196 | r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S); |
154 | r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S); |
| 197 | if(hx>=0) return one-r/x; else return r/x-one; |
155 | if(hx>=0) return one-r/x; else return r/x-one; |
| 198 | } |
156 | } |
| Line 199... | Line 157... | ||
| 199 | 157 | ||
| Line 202... | Line 160... | ||
| 202 | #else |
160 | #else |
| 203 | float erfcf(x) |
161 | float erfcf(x) |
| 204 | float x; |
162 | float x; |
| 205 | #endif |
163 | #endif |
| 206 | { |
164 | { |
| 207 | int32_t hx,ix; |
165 | __int32_t hx,ix; |
| 208 | float R,S,P,Q,s,y,z,r; |
166 | float R,S,P,Q,s,y,z,r; |
| 209 | hx = __get_float_word(x); |
167 | GET_FLOAT_WORD(hx,x); |
| 210 | ix = hx&0x7fffffff; |
168 | ix = hx&0x7fffffff; |
| 211 | if(!(ix<0x7f800000L)) { /* erfc(nan)=nan */ |
169 | if(!FLT_UWORD_IS_FINITE(ix)) { /* erfc(nan)=nan */ |
| 212 | /* erfc(+-inf)=0,2 */ |
170 | /* erfc(+-inf)=0,2 */ |
| 213 | return (float)(((uint32_t)hx>>31)<<1)+one/x; |
171 | return (float)(((__uint32_t)hx>>31)<<1)+one/x; |
| 214 | } |
172 | } |
| Line 215... | Line 173... | ||
| 215 | 173 | ||
| 216 | if(ix < 0x3f580000) { /* |x|<0.84375 */ |
174 | if(ix < 0x3f580000) { /* |x|<0.84375 */ |
| 217 | if(ix < 0x23800000) /* |x|<2**-56 */ |
175 | if(ix < 0x23800000) /* |x|<2**-56 */ |
| Line 236... | Line 194... | ||
| 236 | z = one-erx; return z - P/Q; |
194 | z = one-erx; return z - P/Q; |
| 237 | } else { |
195 | } else { |
| 238 | z = erx+P/Q; return one+z; |
196 | z = erx+P/Q; return one+z; |
| 239 | } |
197 | } |
| 240 | } |
198 | } |
| 241 | - | ||
| 242 | if (ix < 0x41e00000) { /* |x|<28 */ |
199 | if (ix < 0x41e00000) { /* |x|<28 */ |
| 243 | x = fabsf(x); |
200 | x = fabsf(x); |
| 244 | s = one/(x*x); |
201 | s = one/(x*x); |
| 245 | if(ix< 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/ |
202 | if(ix< 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/ |
| 246 | R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( |
203 | R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( |
| Line 252... | Line 209... | ||
| 252 | R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( |
209 | R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( |
| 253 | rb5+s*rb6))))); |
210 | rb5+s*rb6))))); |
| 254 | S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( |
211 | S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( |
| 255 | sb5+s*(sb6+s*sb7)))))); |
212 | sb5+s*(sb6+s*sb7)))))); |
| 256 | } |
213 | } |
| 257 | z = x; |
214 | GET_FLOAT_WORD(ix,x); |
| 258 | __trunc_float_word (&z); |
215 | SET_FLOAT_WORD(z,ix&0xfffff000); |
| 259 | r = __ieee754_expf(-z*z-(float)0.5625)* |
216 | r = __ieee754_expf(-z*z-(float)0.5625)* |
| 260 | __ieee754_expf((z-x)*(z+x)+R/S); |
217 | __ieee754_expf((z-x)*(z+x)+R/S); |
| 261 | if(hx>0) return r/x; else return two-r/x; |
218 | if(hx>0) return r/x; else return two-r/x; |
| 262 | } else { |
219 | } else { |
| 263 | /* set range error */ |
- | |
| 264 | errno = ERANGE; |
- | |
| 265 | if(hx>0) return tiny*tiny; else return two-tiny; |
220 | if(hx>0) return tiny*tiny; else return two-tiny; |
| 266 | } |
221 | } |
| 267 | } |
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) */ |
|