OpenE2K/glibc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Simplify calculation of 2^-m in __mpexp

Browse Source
This commit is contained in:
Siddhesh Poyarekar 2013-01-18 11:18:13 +05:30
parent d3b9ea6148
commit caa99d06e7
3 changed files with 39 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
ChangeLog
View File

@ -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.

30
sysdeps/ieee754/dbl-64/mpa.h
View File

@ -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;
}

33
sysdeps/ieee754/dbl-64/mpexp.c
View File

@ -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. */