diff options
| -rw-r--r-- | ChangeLog | 2 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpexp.c | 149 |
2 files changed, 97 insertions, 54 deletions
diff --git a/ChangeLog b/ChangeLog index 5cd7df0103..1e8a85c1d1 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,7 @@ 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 doubles __mpexp_twomm1. Adjust usage. * sysdeps/ieee754/dbl-64/mpexp.h (__mpexp_twomm1): diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c index ac41ecf66c..feca23c9b6 100644 --- a/sysdeps/ieee754/dbl-64/mpexp.c +++ b/sysdeps/ieee754/dbl-64/mpexp.c @@ -37,17 +37,20 @@ # define SECTION #endif -/* Multi-Precision exponential function subroutine (for p >= 4, */ -/* 2**(-55) <= abs(x) <= 1024). */ +/* Multi-Precision exponential function subroutine (for p >= 4, + 2**(-55) <= abs(x) <= 1024). */ void SECTION -__mpexp(mp_no *x, mp_no *y, int p) { - - int i,j,k,m,m1,m2,n; - 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, - 6,6,6,6,7,7,7,7,8,8,8,8,8}; - static const int m1p[33]= +__mpexp (mp_no *x, mp_no *y, int p) +{ + int i, j, k, m, m1, m2, n; + 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, + 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8 + }; + static const int m1p[33] = { 0, 0, 0, 0, 17, 23, 23, 28, @@ -72,29 +75,55 @@ __mpexp(mp_no *x, mp_no *y, int p) { 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}, - { 0, 0, 0, 0,36,48,60,72, 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,17,23,29,35,41,47,53,59,65, 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, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63}, - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}}; - 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,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; + 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, 24, 32, 40, 48, 56, 64, 72, 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, 0, 0, 23, 28, 33, 38, 42, 47, 52, 57, 62, 66, 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, 0, 0, 0, 0, 0, 0, 0, 0, 43, 47, 50, 54} + }; + 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, 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) */ - 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]*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; } - } + /* 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] * 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; if (__glibc_unlikely (m <= 0)) @@ -112,30 +141,42 @@ __mpexp(mp_no *x, mp_no *y, int p) { break; } - /* Compute s=x*2**(-m). Put result in mps */ - __dbl_mp(a,&mpt1,p); - __mul(x,&mpt1,&mps,p); + /* Compute s=x*2**(-m). Put result in mps. */ + __dbl_mp (a, &mpt1, p); + __mul (x, &mpt1, &mps, p); - /* Evaluate the polynomial. Put result in mpt2 */ - mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=n; - __dvd(&mps,&mpk,&mpt1,p); - __add(&mpone,&mpt1,&mpak,p); - for (k=n-1; k>1; k--) { - __mul(&mps,&mpak,&mpt1,p); - mpk.d[1] = k; - __dvd(&mpt1,&mpk,&mpt2,p); - __add(&mpone,&mpt2,&mpak,p); - } - __mul(&mps,&mpak,&mpt1,p); - __add(&mpone,&mpt1,&mpt2,p); + /* Evaluate the polynomial. Put result in mpt2. */ + mpk.e = 1; + mpk.d[0] = ONE; + mpk.d[1] = n; + __dvd (&mps, &mpk, &mpt1, p); + __add (&mpone, &mpt1, &mpak, p); + for (k = n - 1; k > 1; k--) + { + __mul (&mps, &mpak, &mpt1, p); + mpk.d[1] = k; + __dvd (&mpt1, &mpk, &mpt2, p); + __add (&mpone, &mpt2, &mpak, p); + } + __mul (&mps, &mpak, &mpt1, p); + __add (&mpone, &mpt1, &mpt2, p); - /* Raise polynomial value to the power of 2**m. Put result in y */ - for (k=0,j=0; k<m; ) { - __mul(&mpt2,&mpt2,&mpt1,p); k++; - if (k==m) { j=1; break; } - __mul(&mpt1,&mpt1,&mpt2,p); k++; - } - if (j) __cpy(&mpt1,y,p); - else __cpy(&mpt2,y,p); + /* Raise polynomial value to the power of 2**m. Put result in y. */ + for (k = 0, j = 0; k < m;) + { + __mul (&mpt2, &mpt2, &mpt1, p); + k++; + if (k == m) + { + j = 1; + break; + } + __mul (&mpt1, &mpt1, &mpt2, p); + k++; + } + if (j) + __cpy (&mpt1, y, p); + else + __cpy (&mpt2, y, p); return; } |