diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/e_pow.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_pow.c | 172 |
1 files changed, 24 insertions, 148 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index f6e5fcdbce..542d03a7e3 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -20,13 +20,9 @@ /* MODULE_NAME: upow.c */ /* */ /* FUNCTIONS: upow */ -/* power1 */ -/* my_log2 */ /* log1 */ /* checkint */ /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h */ -/* halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c */ -/* uexp.c upow.c */ /* root.tbl uexp.tbl upow.tbl */ /* An ultimate power routine. Given two IEEE double machine numbers y,x */ /* it computes the correctly rounded (to nearest) value of x^y. */ @@ -50,11 +46,8 @@ static const double huge = 1.0e300, tiny = 1.0e-300; -double __exp1 (double x, double xx, double error); -static double log1 (double x, double *delta, double *error); -static double my_log2 (double x, double *delta, double *error); -double __slowpow (double x, double y, double z); -static double power1 (double x, double y); +double __exp1 (double x, double xx); +static double log1 (double x, double *delta); static int checkint (double x); /* An ultimate power routine. Given two IEEE double machine numbers y, x it @@ -63,7 +56,7 @@ double SECTION __ieee754_pow (double x, double y) { - double z, a, aa, error, t, a1, a2, y1, y2; + double z, a, aa, t, a1, a2, y1, y2; mynumber u, v; int k; int4 qx, qy; @@ -100,7 +93,7 @@ __ieee754_pow (double x, double y) not matter if |y| <= 2**-64. */ if (fabs (y) < 0x1p-64) y = y < 0 ? -0x1p-64 : 0x1p-64; - z = log1 (x, &aa, &error); /* x^y =e^(y log (X)) */ + z = log1 (x, &aa); /* x^y =e^(y log (X)) */ t = y * CN; y1 = t - (t - y); y2 = y - y1; @@ -111,9 +104,16 @@ __ieee754_pow (double x, double y) aa = y2 * a1 + y * a2; a1 = a + aa; a2 = (a - a1) + aa; - error = error * fabs (y); - t = __exp1 (a1, a2, 1.9e16 * error); /* return -10 or 0 if wasn't computed exactly */ - retval = (t > 0) ? t : power1 (x, y); + + /* Maximum relative error RElog of log1 is 1.0e-21 (69.7 bits). + Maximum relative error REexp of __exp1 is 8.8e-22 (69.9 bits). + We actually compute exp ((1 + RElog) * log (x) * y) * (1 + REexp). + Since RElog/REexp are tiny and log (x) * y is at most log (DBL_MAX), + this is equivalent to pow (x, y) * (1 + 710 * RElog + REexp). + So the relative error is 710 * 1.0e-21 + 8.8e-22 = 7.1e-19 + (60.2 bits). The worst-case ULP error is 0.5064. */ + + retval = __exp1 (a1, a2); } if (isinf (retval)) @@ -218,33 +218,11 @@ __ieee754_pow (double x, double y) strong_alias (__ieee754_pow, __pow_finite) #endif -/* Compute x^y using more accurate but more slow log routine. */ -static double -SECTION -power1 (double x, double y) -{ - double z, a, aa, error, t, a1, a2, y1, y2; - z = my_log2 (x, &aa, &error); - t = y * CN; - y1 = t - (t - y); - y2 = y - y1; - t = z * CN; - a1 = t - (t - z); - a2 = z - a1; - a = y * z; - aa = ((y1 * a1 - a) + y1 * a2 + y2 * a1) + y2 * a2 + aa * y; - a1 = a + aa; - a2 = (a - a1) + aa; - error = error * fabs (y); - t = __exp1 (a1, a2, 1.9e16 * error); - return (t >= 0) ? t : __slowpow (x, y, z); -} - /* Compute log(x) (x is left argument). The result is the returned double + the - parameter DELTA. The result is bounded by ERROR. */ + parameter DELTA. */ static double SECTION -log1 (double x, double *delta, double *error) +log1 (double x, double *delta) { unsigned int i, j; int m; @@ -260,9 +238,7 @@ log1 (double x, double *delta, double *error) u.x = x; m = u.i[HIGH_HALF]; - *error = 0; - *delta = 0; - if (m < 0x00100000) /* 1<x<2^-1007 */ + if (m < 0x00100000) /* Handle denormal x. */ { x = x * t52.x; add = -52.0; @@ -284,7 +260,7 @@ log1 (double x, double *delta, double *error) v.x = u.x + bigu.x; uu = v.x - bigu.x; i = (v.i[LOW_HALF] & 0x000003ff) << 2; - if (two52.i[LOW_HALF] == 1023) /* nx = 0 */ + if (two52.i[LOW_HALF] == 1023) /* Exponent of x is 0. */ { if (i > 1192 && i < 1208) /* |x-1| < 1.5*2**-10 */ { @@ -296,8 +272,8 @@ log1 (double x, double *delta, double *error) * (r7 + t * r8))))) - 0.5 * t2 * (t + t1)); res = e1 + e2; - *error = 1.0e-21 * fabs (t); *delta = (e1 - res) + e2; + /* Max relative error is 1.464844e-24, so accurate to 79.1 bits. */ return res; } /* |x-1| < 1.5*2**-10 */ else @@ -316,12 +292,12 @@ log1 (double x, double *delta, double *error) t2 = ((((t - t1) + e) + (ui.x[i + 3] + vj.x[j + 2])) + e2 + e * e * (p2 + e * (p3 + e * p4))); res = t1 + t2; - *error = 1.0e-24; *delta = (t1 - res) + t2; + /* Max relative error is 1.0e-24, so accurate to 79.7 bits. */ return res; } - } /* nx = 0 */ - else /* nx != 0 */ + } + else /* Exponent of x != 0. */ { eps = u.x - uu; nx = (two52.x - two52e.x) + add; @@ -334,113 +310,13 @@ log1 (double x, double *delta, double *error) t2 = ((((t - t1) + e) + nx * ln2b.x + ui.x[i + 3] + e2) + e * e * (q2 + e * (q3 + e * (q4 + e * (q5 + e * q6))))); res = t1 + t2; - *error = 1.0e-21; - *delta = (t1 - res) + t2; - return res; - } /* nx != 0 */ -} - -/* Slower but more accurate routine of log. The returned result is double + - DELTA. The result is bounded by ERROR. */ -static double -SECTION -my_log2 (double x, double *delta, double *error) -{ - unsigned int i, j; - int m; - double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0; - double ou1, ou2, lu1, lu2, ov, lv1, lv2, a, a1, a2; - double y, yy, z, zz, j1, j2, j7, j8; -#ifndef DLA_FMS - double j3, j4, j5, j6; -#endif - mynumber u, v; -#ifdef BIG_ENDI - mynumber /**/ two52 = {{0x43300000, 0x00000000}}; /* 2**52 */ -#else -# ifdef LITTLE_ENDI - mynumber /**/ two52 = {{0x00000000, 0x43300000}}; /* 2**52 */ -# endif -#endif - - u.x = x; - m = u.i[HIGH_HALF]; - *error = 0; - *delta = 0; - add = 0; - if (m < 0x00100000) - { /* x < 2^-1022 */ - x = x * t52.x; - add = -52.0; - u.x = x; - m = u.i[HIGH_HALF]; - } - - if ((m & 0x000fffff) < 0x0006a09e) - { - u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3ff00000; - two52.i[LOW_HALF] = (m >> 20); - } - else - { - u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3fe00000; - two52.i[LOW_HALF] = (m >> 20) + 1; - } - - v.x = u.x + bigu.x; - uu = v.x - bigu.x; - i = (v.i[LOW_HALF] & 0x000003ff) << 2; - /*------------------------------------- |x-1| < 2**-11------------------------------- */ - if ((two52.i[LOW_HALF] == 1023) && (i == 1200)) - { - t = x - 1.0; - EMULV (t, s3, y, yy, j1, j2, j3, j4, j5); - ADD2 (-0.5, 0, y, yy, z, zz, j1, j2); - MUL2 (t, 0, z, zz, y, yy, j1, j2, j3, j4, j5, j6, j7, j8); - MUL2 (t, 0, y, yy, z, zz, j1, j2, j3, j4, j5, j6, j7, j8); - - e1 = t + z; - e2 = ((((t - e1) + z) + zz) + t * t * t - * (ss3 + t * (s4 + t * (s5 + t * (s6 + t * (s7 + t * s8)))))); - res = e1 + e2; - *error = 1.0e-25 * fabs (t); - *delta = (e1 - res) + e2; - return res; - } - /*----------------------------- |x-1| > 2**-11 -------------------------- */ - else - { /*Computing log(x) according to log table */ - nx = (two52.x - two52e.x) + add; - ou1 = ui.x[i]; - ou2 = ui.x[i + 1]; - lu1 = ui.x[i + 2]; - lu2 = ui.x[i + 3]; - v.x = u.x * (ou1 + ou2) + bigv.x; - vv = v.x - bigv.x; - j = v.i[LOW_HALF] & 0x0007ffff; - j = j + j + j; - eps = u.x - uu * vv; - ov = vj.x[j]; - lv1 = vj.x[j + 1]; - lv2 = vj.x[j + 2]; - a = (ou1 + ou2) * (1.0 + ov); - a1 = (a + 1.0e10) - 1.0e10; - a2 = a * (1.0 - a1 * uu * vv); - e1 = eps * a1; - e2 = eps * a2; - e = e1 + e2; - e2 = (e1 - e) + e2; - t = nx * ln2a.x + lu1 + lv1; - t1 = t + e; - t2 = ((((t - t1) + e) + (lu2 + lv2 + nx * ln2b.x + e2)) + e * e - * (p2 + e * (p3 + e * p4))); - res = t1 + t2; - *error = 1.0e-27; *delta = (t1 - res) + t2; + /* Max relative error is 1.0e-21, so accurate to 69.7 bits. */ return res; } } + /* This function receives a double x and checks if it is an integer. If not, it returns 0, else it returns 1 if even or -1 if odd. */ static int |