- GET_HIGH_WORD(hx, x);
- ix = hx & 0x7fffffff;
- if (ix >= 0x7ff00000) {
- /* erf(nan)=nan, erf(+-inf)=+-1 */
- i = ((uint32_t)hx>>31)<<1;
- return (double)(1-i) + 1.0/x;
- }
- if (ix < 0x3feb0000) { /* |x|<0.84375 */
- if (ix < 0x3e300000) { /* |x|<2**-28 */
- if (ix < 0x00800000)
- /* avoid underflow */
- return 0.125*(8.0*x + efx8*x);
- return x + efx*x;
- }
- z = x*x;
- r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
- s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
- y = r/s;
- return x + x*y;
- }
- if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
- s = fabs(x)-1.0;
- P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
- Q = 1.0+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
- if (hx >= 0)
- return erx + P/Q;
- return -erx - P/Q;
- }
- if (ix >= 0x40180000) { /* inf > |x| >= 6 */
- if (hx >= 0)
- return 1.0 - tiny;
- return tiny - 1.0;
- }