mirror/glibc
2
0
Fork 0
mirror of git://sourceware.org/git/glibc.git synced 2025-01-18 12:16:13 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix formatting in mpexp.c

Browse Source
This commit is contained in:
Siddhesh Poyarekar 2013-01-10 14:56:04 +05:30
parent 751b85f795
commit 7490eb81ae
2 changed files with 97 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ChangeLog
View File

@ -1,5 +1,7 @@
2013-01-10 Siddhesh Poyarekar <siddhesh@redhat.com> 2013-01-10 Siddhesh Poyarekar <siddhesh@redhat.com>
* sysdeps/ieee754/dbl-64/mpexp.c: Fix formatting.
* sysdeps/ieee754/dbl-64/mpexp.c (__mpexp): New array of * sysdeps/ieee754/dbl-64/mpexp.c (__mpexp): New array of
doubles __mpexp_twomm1. Adjust usage. doubles __mpexp_twomm1. Adjust usage.
* sysdeps/ieee754/dbl-64/mpexp.h (__mpexp_twomm1): * sysdeps/ieee754/dbl-64/mpexp.h (__mpexp_twomm1):

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

@ -37,17 +37,20 @@
# define SECTION # define SECTION
#endif #endif
/* Multi-Precision exponential function subroutine (for p >= 4, */ /* Multi-Precision exponential function subroutine (for p >= 4,
/* 2**(-55) <= abs(x) <= 1024). */ 2**(-55) <= abs(x) <= 1024). */
void void
SECTION SECTION
__mpexp(mp_no *x, mp_no *y, int p) { __mpexp (mp_no *x, mp_no *y, int p)
{
int i,j,k,m,m1,m2,n; int i, j, k, m, m1, m2, n;
double a,b; double a, 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, static const int np[33] =
6,6,6,6,7,7,7,7,8,8,8,8,8}; {
static const int m1p[33]= 0, 0, 0, 0, 3, 3, 4, 4, 5, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8
};
static const int m1p[33] =
{ {
0, 0, 0, 0, 0, 0, 0, 0,
17, 23, 23, 28, 17, 23, 23, 28,
@ -72,28 +75,54 @@ __mpexp(mp_no *x, mp_no *y, int p) {
0x1.0p-70, 0x1.0p-73, 0x1.0p-76, 0x1.0p-78, 0x1.0p-70, 0x1.0p-73, 0x1.0p-76, 0x1.0p-78,
0x1.0p-81 0x1.0p-81
}; };
static const int m1np[7][18] = { static const int m1np[7][18] =
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {
{ 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 36, 48, 60, 72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 24, 32, 40, 48, 56, 64, 72, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0}, {0, 0, 0, 0, 17, 23, 29, 35, 41, 47, 53, 59, 65, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63}, {0, 0, 0, 0, 0, 0, 23, 28, 33, 38, 42, 47, 52, 57, 62, 66, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}}; {0, 0, 0, 0, 0, 0, 0, 0, 27, 0, 0, 39, 43, 47, 51, 55, 59, 63},
mp_no mpk = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 43, 47, 50, 54}
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, };
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; mp_no mpk =
mp_no mps,mpak,mpt1,mpt2; {
0,
{
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
}
};
mp_no mps, mpak, mpt1, mpt2;
/* Choose m,n and compute a=2**(-m) */ /* Choose m,n and compute a=2**(-m). */
n = np[p]; m1 = m1p[p]; a = __mpexp_twomm1[p]; n = np[p];
for (i=0; i<EX; i++) a *= RADIXI; m1 = m1p[p];
for ( ; i>EX; i--) a *= RADIX; a = __mpexp_twomm1[p];
b = X[1]*RADIXI; m2 = 24*EX; for (i = 0; i < EX; i++)
for (; b<HALF; m2--) { a *= TWO; b *= TWO; } a *= RADIXI;
if (b == HALF) { for (; i > EX; i--)
for (i=2; i<=p; i++) { if (X[i]!=ZERO) break; } a *= RADIX;
if (i==p+1) { m2--; a *= TWO; } b = X[1] * RADIXI;
m2 = 24 * EX;
for (; b < HALF; m2--)
{
a *= TWO;
b *= TWO;
}
if (b == HALF)
{
for (i = 2; i <= p; i++)
{
if (X[i] != ZERO)
break;
}
if (i == p + 1)
{
m2--;
a *= TWO;
}
} }
m = m1 + m2; m = m1 + m2;
@ -112,30 +141,42 @@ __mpexp(mp_no *x, mp_no *y, int p) {
break; break;
} }
/* Compute s=x*2**(-m). Put result in mps */ /* Compute s=x*2**(-m). Put result in mps. */
__dbl_mp(a,&mpt1,p); __dbl_mp (a, &mpt1, p);
__mul(x,&mpt1,&mps,p); __mul (x, &mpt1, &mps, p);
/* Evaluate the polynomial. Put result in mpt2 */ /* Evaluate the polynomial. Put result in mpt2. */
mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=n; mpk.e = 1;
__dvd(&mps,&mpk,&mpt1,p); mpk.d[0] = ONE;
__add(&mpone,&mpt1,&mpak,p); mpk.d[1] = n;
for (k=n-1; k>1; k--) { __dvd (&mps, &mpk, &mpt1, p);
__mul(&mps,&mpak,&mpt1,p); __add (&mpone, &mpt1, &mpak, p);
for (k = n - 1; k > 1; k--)
{
__mul (&mps, &mpak, &mpt1, p);
mpk.d[1] = k; mpk.d[1] = k;
__dvd(&mpt1,&mpk,&mpt2,p); __dvd (&mpt1, &mpk, &mpt2, p);
__add(&mpone,&mpt2,&mpak,p); __add (&mpone, &mpt2, &mpak, p);
} }
__mul(&mps,&mpak,&mpt1,p); __mul (&mps, &mpak, &mpt1, p);
__add(&mpone,&mpt1,&mpt2,p); __add (&mpone, &mpt1, &mpt2, p);
/* Raise polynomial value to the power of 2**m. Put result in y */ /* Raise polynomial value to the power of 2**m. Put result in y. */
for (k=0,j=0; k<m; ) { for (k = 0, j = 0; k < m;)
__mul(&mpt2,&mpt2,&mpt1,p); k++; {
if (k==m) { j=1; break; } __mul (&mpt2, &mpt2, &mpt1, p);
__mul(&mpt1,&mpt1,&mpt2,p); k++; k++;
if (k == m)
{
j = 1;
break;
} }
if (j) __cpy(&mpt1,y,p); __mul (&mpt1, &mpt1, &mpt2, p);
else __cpy(&mpt2,y,p); k++;
}
if (j)
__cpy (&mpt1, y, p);
else
__cpy (&mpt2, y, p);
return; return;
} }