diff options
| -rw-r--r-- | ChangeLog | 5 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpa.h | 30 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpexp.c | 33 |
3 files changed, 39 insertions, 29 deletions
@@ -1,5 +1,10 @@ 2013-01-18 Siddhesh Poyarekar <siddhesh@redhat.com> + * sysdeps/ieee754/dbl-64/mpa.h (__pow_mp): New function to get an + mp_no from a power of two. + * sysdeps/ieee754/dbl-64/mpexp.c (__mpexp): Remove + __mpexp_twomm1. Use __pow_mp. + * sysdeps/ieee754/dbl-64/mpexp.c (__mpexp): Remove unnecessary multiplication. diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h index debb3b2..06343d4 100644 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ b/sysdeps/ieee754/dbl-64/mpa.h @@ -123,3 +123,33 @@ extern void __mpsqrt (mp_no *, mp_no *, int); extern void __mpexp (mp_no *, mp_no *, int); extern void __c32 (mp_no *, mp_no *, mp_no *, int); extern int __mpranred (double, mp_no *, int); + +/* Given a power POW, build a multiprecision number 2^POW. */ +static inline void +__pow_mp (int pow, mp_no *y, int p) +{ + int i, rem; + + /* The exponent is E such that E is a factor of 2^24. The remainder (of the + form 2^x) goes entirely into the first digit of the mantissa as it is + always less than 2^24. */ + EY = pow / 24; + rem = pow - EY * 24; + EY++; + + /* If the remainder is negative, it means that POW was negative since + |EY * 24| <= |pow|. Adjust so that REM is positive and still less than + 24 because of which, the mantissa digit is less than 2^24. */ + if (rem < 0) + { + EY--; + rem += 24; + } + /* The sign of any 2^x is always positive. */ + Y[0] = ONE; + Y[1] = 1 << rem; + + /* Everything else is ZERO. */ + for (i = 2; i <= p; i++) + Y[i] = ZERO; +} diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c index 811fb46..8d288ff 100644 --- a/sysdeps/ieee754/dbl-64/mpexp.c +++ b/sysdeps/ieee754/dbl-64/mpexp.c @@ -43,7 +43,7 @@ SECTION __mpexp (mp_no *x, mp_no *y, int p) { int i, j, k, m, m1, m2, n; - double a, b; + double b; static const int np[33] = { 0, 0, 0, 0, 3, 3, 4, 4, 5, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6, @@ -61,19 +61,6 @@ __mpexp (mp_no *x, mp_no *y, int p) 70, 73, 76, 78, 81 }; - /* Stored values for 2^-m, where values of m are defined in M1P above. */ - static const double __mpexp_twomm1[33] = - { - 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, - 0x1.0p-17, 0x1.0p-23, 0x1.0p-23, 0x1.0p-28, - 0x1.0p-27, 0x1.0p-38, 0x1.0p-42, 0x1.0p-39, - 0x1.0p-43, 0x1.0p-47, 0x1.0p-43, 0x1.0p-47, - 0x1.0p-50, 0x1.0p-54, 0x1.0p-57, 0x1.0p-60, - 0x1.0p-64, 0x1.0p-67, 0x1.0p-71, 0x1.0p-74, - 0x1.0p-68, 0x1.0p-71, 0x1.0p-74, 0x1.0p-77, - 0x1.0p-70, 0x1.0p-73, 0x1.0p-76, 0x1.0p-78, - 0x1.0p-81 - }; static const int m1np[7][18] = { {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, @@ -98,18 +85,10 @@ __mpexp (mp_no *x, mp_no *y, int p) /* Choose m,n and compute a=2**(-m). */ n = np[p]; m1 = m1p[p]; - a = __mpexp_twomm1[p]; - for (i = 0; i < EX; i++) - a *= RADIXI; - for (; i > EX; i--) - a *= RADIX; b = X[1]; m2 = 24 * EX; for (; b < HALFRAD; m2--) - { - a *= TWO; - b *= TWO; - } + b *= TWO; if (b == HALFRAD) { for (i = 2; i <= p; i++) @@ -118,10 +97,7 @@ __mpexp (mp_no *x, mp_no *y, int p) break; } if (i == p + 1) - { - m2--; - a *= TWO; - } + m2--; } m = m1 + m2; @@ -134,14 +110,13 @@ __mpexp (mp_no *x, mp_no *y, int p) than 2^-55. */ assert (p < 18); m = 0; - a = ONE; for (i = n - 1; i > 0; i--, n--) if (m1np[i][p] + m2 > 0) break; } /* Compute s=x*2**(-m). Put result in mps. */ - __dbl_mp (a, &mpt1, p); + __pow_mp (-m, &mpt1, p); __mul (x, &mpt1, &mps, p); /* Evaluate the polynomial. Put result in mpt2. */ |