diff options
Diffstat (limited to 'libc/sysdeps/ieee754/dbl-64/mpexp.c')
| -rw-r--r-- | libc/sysdeps/ieee754/dbl-64/mpexp.c | 31 |
1 files changed, 20 insertions, 11 deletions
diff --git a/libc/sysdeps/ieee754/dbl-64/mpexp.c b/libc/sysdeps/ieee754/dbl-64/mpexp.c index 6b1fcf255..c4048207e 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpexp.c +++ b/libc/sysdeps/ieee754/dbl-64/mpexp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2011 Free Software Foundation + * Copyright (C) 2001-2013 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -31,6 +31,7 @@ #include "endian.h" #include "mpa.h" #include "mpexp.h" +#include <assert.h> #ifndef SECTION # define SECTION @@ -56,9 +57,6 @@ __mpexp(mp_no *x, mp_no *y, int p) { { 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 mpone = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 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 = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; @@ -74,23 +72,34 @@ __mpexp(mp_no *x, mp_no *y, int p) { for (i=2; i<=p; i++) { if (X[i]!=ZERO) break; } if (i==p+1) { m2--; a *= TWO; } } - if ((m=m1+m2) <= 0) { - m=0; a=ONE; - for (i=n-1; i>0; i--,n--) { if (m1np[i][p]+m2>0) break; } - } + + m = m1 + m2; + if (__glibc_unlikely (m <= 0)) + { + /* The m1np array which is used to determine if we can reduce the + polynomial expansion iterations, has only 18 elements. Besides, + numbers smaller than those required by p >= 18 should not come here + at all since the fast phase of exp returns 1.0 for anything less + 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); __mul(x,&mpt1,&mps,p); /* Evaluate the polynomial. Put result in mpt2 */ - mpone.e=1; mpone.d[0]=ONE; mpone.d[1]=ONE; - mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=__mpexp_nn[n].d; + 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]=__mpexp_nn[k].d; + mpk.d[1] = k; __dvd(&mpt1,&mpk,&mpt2,p); __add(&mpone,&mpt2,&mpak,p); } |