dns response handling: don't treat too many addresses as an error
[musl] / src / math / erfl.c
index 2fd3c44..e267c23 100644 (file)
@@ -253,8 +253,8 @@ static long double erfc1(long double x)
 
 static long double erfc2(uint32_t ix, long double x)
 {
+       union ldshape u;
        long double s,z,R,S;
-       uint32_t i0,i1;
 
        if (ix < 0x3fffa000)  /* 0.84375 <= |x| < 1.25 */
                return erfc1(x);
@@ -266,50 +266,33 @@ static long double erfc2(uint32_t ix, long double x)
                     s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
                S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
                     s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
-       } else {  /* 2.857 <= |x| */
+       } else if (ix < 0x4001d555) {  /* 2.857 <= |x| < 6.6666259765625 */
                R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
                     s * (rb[5] + s * (rb[6] + s * rb[7]))))));
                S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
                     s * (sb[5] + s * (sb[6] + s))))));
-       }
-       if (ix < 0x4000b6db) {  /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
-               R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
-                    s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
-               S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
-                    s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
-       } else if (ix < 0x4001d555) {  /* 6.6666259765625 > |x| >= 1/.35 ~ 2.857143 */
-               R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
-                    s * (rb[5] + s * (rb[6] + s * rb[7]))))));
-               S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
-                    s * (sb[5] + s * (sb[6] + s))))));
-       } else { /* 107 > |x| >= 6.666 */
+       } else { /* 6.666 <= |x| < 107 (erfc only) */
                R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
                     s * (rc[4] + s * rc[5]))));
                S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
                     s * (sc[4] + s))));
        }
-       z = x;
-       GET_LDOUBLE_WORDS(ix, i0, i1, z);
-       i1 = 0;
-       i0 &= 0xffffff00;
-       SET_LDOUBLE_WORDS(z, ix, i0, i1);
+       u.f = x;
+       u.i.m &= -1ULL << 40;
+       z = u.f;
        return expl(-z*z - 0.5625) * expl((z - x) * (z + x) + R / S) / x;
 }
 
 long double erfl(long double x)
 {
        long double r, s, z, y;
-       uint32_t i0, i1, ix;
-       int sign;
+       union ldshape u = {x};
+       uint32_t ix = (u.i.se & 0x7fffU)<<16 | u.i.m>>48;
+       int sign = u.i.se >> 15;
 
-       GET_LDOUBLE_WORDS(ix, i0, i1, x);
-       sign = ix >> 15;
-       ix &= 0x7fff;
-       if (ix == 0x7fff) {
+       if (ix >= 0x7fff0000)
                /* erf(nan)=nan, erf(+-inf)=+-1 */
                return 1 - 2*sign + 1/x;
-       }
-       ix = (ix << 16) | (i0 >> 16);
        if (ix < 0x3ffed800) {  /* |x| < 0.84375 */
                if (ix < 0x3fde8000) {  /* |x| < 2**-33 */
                        return 0.125 * (8 * x + efx8 * x);  /* avoid underflow */
@@ -332,17 +315,13 @@ long double erfl(long double x)
 long double erfcl(long double x)
 {
        long double r, s, z, y;
-       uint32_t i0, i1, ix;
-       int sign;
+       union ldshape u = {x};
+       uint32_t ix = (u.i.se & 0x7fffU)<<16 | u.i.m>>48;
+       int sign = u.i.se >> 15;
 
-       GET_LDOUBLE_WORDS(ix, i0, i1, x);
-       sign = ix>>15;
-       ix &= 0x7fff;
-       if (ix == 0x7fff)
+       if (ix >= 0x7fff0000)
                /* erfc(nan) = nan, erfc(+-inf) = 0,2 */
                return 2*sign + 1/x;
-
-       ix = (ix << 16) | (i0 >> 16);
        if (ix < 0x3ffed800) {  /* |x| < 0.84375 */
                if (ix < 0x3fbe0000)  /* |x| < 2**-65 */
                        return 1.0 - x;
@@ -358,6 +337,17 @@ long double erfcl(long double x)
        }
        if (ix < 0x4005d600)  /* |x| < 107 */
                return sign ? 2 - erfc2(ix,x) : erfc2(ix,x);
-       return sign ? 2 - 0x1p-16382L : 0x1p-16382L*0x1p-16382L;
+       y = 0x1p-16382L;
+       return sign ? 2 - y : y*y;
+}
+#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
+// TODO: broken implementation to make things compile
+long double erfl(long double x)
+{
+       return erf(x);
+}
+long double erfcl(long double x)
+{
+       return erfc(x);
 }
 #endif