Mercurial > hg > octave-lojdl > gnulib-hg
changeset 8494:b556de32ce5a
Remove Paolo's variant of the algorithm, with his permission.
author | Bruno Haible <bruno@clisp.org> |
---|---|
date | Thu, 22 Mar 2007 11:52:34 +0000 |
parents | c62e99bf6912 |
children | eff68e94da66 |
files | ChangeLog lib/frexp.c |
diffstat | 2 files changed, 79 insertions(+), 115 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,7 @@ +2007-03-22 Bruno Haible <bruno@clisp.org> + + * lib/frexp.c: Remove older implementation that uses divisions. + 2007-03-21 Bruno Haible <bruno@clisp.org> * modules/frexp-tests: New file.
--- a/lib/frexp.c +++ b/lib/frexp.c @@ -15,6 +15,9 @@ with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ +/* Written by Paolo Bonzini <bonzini@gnu.org>, 2003, and + Bruno Haible <bruno@clisp.org>, 2007. */ + #include <config.h> #if !(defined USE_LONG_DOUBLE && !HAVE_LONG_DOUBLE) @@ -65,129 +68,86 @@ sign = -1; } - if (0) - { - /* Implementation contributed by Paolo Bonzini. - Disabled because it's under GPL and doesn't pass the tests. */ - - /* Since the exponent is an 'int', it fits in 64 bits. Therefore the - loops are executed no more than 64 times. */ - DOUBLE exponents[64]; - DOUBLE *next; - int bit; + { + /* Since the exponent is an 'int', it fits in 64 bits. Therefore the + loops are executed no more than 64 times. */ + DOUBLE pow2[64]; /* pow2[i] = 2^2^i */ + DOUBLE powh[64]; /* powh[i] = 2^-2^i */ + int i; - exponent = 0; - if (x >= L_(1.0)) - { - for (next = exponents, exponents[0] = L_(2.0), bit = 1; - *next <= x + x; - bit <<= 1, next[1] = next[0] * next[0], next++); + exponent = 0; + if (x >= L_(1.0)) + { + /* A positive exponent. */ + DOUBLE pow2_i; /* = pow2[i] */ + DOUBLE powh_i; /* = powh[i] */ - for (; next >= exponents; bit >>= 1, next--) - if (x + x >= *next) + /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, + x * 2^exponent = argument, x >= 1.0. */ + for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); + ; + i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) + { + if (x >= pow2_i) { - x /= *next; - exponent |= bit; + exponent += (1 << i); + x *= powh_i; } - } - else if (x < L_(0.5)) - { - for (next = exponents, exponents[0] = L_(0.5), bit = 1; - *next > x; - bit <<= 1, next[1] = next[0] * next[0], next++); - - for (; next >= exponents; bit >>= 1, next--) - if (x < *next) - { - x /= *next; - exponent |= bit; - } - exponent = - exponent; - } - } - else - { - /* Implementation contributed by Bruno Haible. */ - - /* Since the exponent is an 'int', it fits in 64 bits. Therefore the - loops are executed no more than 64 times. */ - DOUBLE pow2[64]; /* pow2[i] = 2^2^i */ - DOUBLE powh[64]; /* powh[i] = 2^-2^i */ - int i; - - exponent = 0; - if (x >= L_(1.0)) - { - /* A positive exponent. */ - DOUBLE pow2_i; /* = pow2[i] */ - DOUBLE powh_i; /* = powh[i] */ + else + break; - /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, - x * 2^exponent = argument, x >= 1.0. */ - for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); - ; - i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) - { - if (x >= pow2_i) - { - exponent += (1 << i); - x *= powh_i; - } - else - break; + pow2[i] = pow2_i; + powh[i] = powh_i; + } + /* Avoid making x too small, as it could become a denormalized + number and thus lose precision. */ + while (i > 0 && x < pow2[i - 1]) + { + i--; + powh_i = powh[i]; + } + exponent += (1 << i); + x *= powh_i; + /* Here 2^-2^i <= x < 1.0. */ + } + else + { + /* A negative or zero exponent. */ + DOUBLE pow2_i; /* = pow2[i] */ + DOUBLE powh_i; /* = powh[i] */ - pow2[i] = pow2_i; - powh[i] = powh_i; - } - /* Avoid making x too small, as it could become a denormalized - number and thus lose precision. */ - while (i > 0 && x < pow2[i - 1]) - { - i--; - powh_i = powh[i]; - } - exponent += (1 << i); - x *= powh_i; - /* Here 2^-2^i <= x < 1.0. */ - } - else - { - /* A negative or zero exponent. */ - DOUBLE pow2_i; /* = pow2[i] */ - DOUBLE powh_i; /* = powh[i] */ + /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, + x * 2^exponent = argument, x < 1.0. */ + for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); + ; + i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) + { + if (x < powh_i) + { + exponent -= (1 << i); + x *= pow2_i; + } + else + break; - /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, - x * 2^exponent = argument, x < 1.0. */ - for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); - ; - i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) - { - if (x < powh_i) - { - exponent -= (1 << i); - x *= pow2_i; - } - else - break; + pow2[i] = pow2_i; + powh[i] = powh_i; + } + /* Here 2^-2^i <= x < 1.0. */ + } - pow2[i] = pow2_i; - powh[i] = powh_i; - } - /* Here 2^-2^i <= x < 1.0. */ - } - - /* Invariants: x * 2^exponent = argument, and 2^-2^i <= x < 1.0. */ - while (i > 0) - { - i--; - if (x < powh[i]) - { - exponent -= (1 << i); - x *= pow2[i]; - } - } - /* Here 0.5 <= x < 1.0. */ - } + /* Invariants: x * 2^exponent = argument, and 2^-2^i <= x < 1.0. */ + while (i > 0) + { + i--; + if (x < powh[i]) + { + exponent -= (1 << i); + x *= pow2[i]; + } + } + /* Here 0.5 <= x < 1.0. */ + } *exp = exponent; return (sign < 0 ? - x : x);