From 7490eb81ae6b3d9e2c9784fd845e0fc2d5484292 Mon Sep 17 00:00:00 2001 From: Siddhesh Poyarekar Date: Thu, 10 Jan 2013 14:56:04 +0530 Subject: Fix formatting in mpexp.c --- sysdeps/ieee754/dbl-64/mpexp.c | 149 ++++++++++++++++++++++++++--------------- 1 file changed, 95 insertions(+), 54 deletions(-) (limited to 'sysdeps/ieee754') 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; iEX; i--) a *= RADIX; - b = X[1]*RADIXI; m2 = 24*EX; - for (; b 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