From 5b2ff9570fa75d56e797184c694536d3b08f9977 Mon Sep 17 00:00:00 2001 From: Bruno Haible Date: Thu, 22 Mar 2007 11:52:34 +0000 Subject: [PATCH] Remove Paolo's variant of the algorithm, with his permission. --- ChangeLog | 4 ++ lib/frexp.c | 198 ++++++++++++++++++++++++------------------------------------ 2 files changed, 83 insertions(+), 119 deletions(-) diff --git a/ChangeLog b/ChangeLog index b84bd1c65..b280df2c3 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,7 @@ +2007-03-22 Bruno Haible + + * lib/frexp.c: Remove older implementation that uses divisions. + 2007-03-21 Bruno Haible * modules/frexp-tests: New file. diff --git a/lib/frexp.c b/lib/frexp.c index b4b9f8385..3074f2e15 100644 --- 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 , 2003, and + Bruno Haible , 2007. */ + #include #if !(defined USE_LONG_DOUBLE && !HAVE_LONG_DOUBLE) @@ -65,129 +68,86 @@ FUNC (DOUBLE x, int *exp) 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; - - 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++); - - for (; next >= exponents; bit >>= 1, next--) - if (x + x >= *next) + { + /* 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] */ + + /* 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) + 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] */ + + /* 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) { - x /= *next; - exponent |= bit; + exponent -= (1 << i); + x *= pow2_i; } - 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] */ - - /* 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] */ - - /* 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. */ - } - - /* 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. */ - } + else + break; + + 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. */ + } *exp = exponent; return (sign < 0 ? - x : x); -- 2.11.0