diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/mplog.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mplog.c | 44 |
1 files changed, 23 insertions, 21 deletions
diff --git a/sysdeps/ieee754/dbl-64/mplog.c b/sysdeps/ieee754/dbl-64/mplog.c index e3d10846e2..f8d5c1095f 100644 --- a/sysdeps/ieee754/dbl-64/mplog.c +++ b/sysdeps/ieee754/dbl-64/mplog.c @@ -1,4 +1,3 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. @@ -37,27 +36,30 @@ #include "endian.h" #include "mpa.h" -void __mpexp(mp_no *, mp_no *, int); - -void __mplog(mp_no *x, mp_no *y, int p) { - int i,m; - static const int mp[33] = {0,0,0,0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3, - 4,4,4,4,4,4,4,4,4,4,4,4,4,4}; - mp_no mpt1,mpt2; +void +__mplog (mp_no *x, mp_no *y, int p) +{ + int i, m; + static const int mp[33] = + { + 0, 0, 0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 + }; + mp_no mpt1, mpt2; - /* Choose m */ + /* Choose m. */ m = mp[p]; - /* Perform m newton iterations to solve for y: exp(y)-x=0. */ - /* The iterations formula is: y(n+1)=y(n)+(x*exp(-y(n))-1). */ - __cpy(y,&mpt1,p); - for (i=0; i<m; i++) { - mpt1.d[0]=-mpt1.d[0]; - __mpexp(&mpt1,&mpt2,p); - __mul(x,&mpt2,&mpt1,p); - __sub(&mpt1,&mpone,&mpt2,p); - __add(y,&mpt2,&mpt1,p); - __cpy(&mpt1,y,p); - } - return; + /* Perform m newton iterations to solve for y: exp(y) - x = 0. The + iterations formula is: y(n + 1) = y(n) + (x * exp(-y(n)) - 1). */ + __cpy (y, &mpt1, p); + for (i = 0; i < m; i++) + { + mpt1.d[0] = -mpt1.d[0]; + __mpexp (&mpt1, &mpt2, p); + __mul (x, &mpt2, &mpt1, p); + __sub (&mpt1, &mpone, &mpt2, p); + __add (y, &mpt2, &mpt1, p); + __cpy (&mpt1, y, p); + } } |