Mercurial > hg > octave-nkf > gnulib-hg
changeset 8589:3334000e84a6
Set the FPU control word as needed for 'long double' computations.
| author | Bruno Haible <bruno@clisp.org> |
|---|---|
| date | Fri, 30 Mar 2007 00:19:34 +0000 |
| parents | 84d2dbfe5855 |
| children | ebe4e154c5db |
| files | ChangeLog lib/ldexpl.c modules/ldexpl |
| diffstat | 3 files changed, 34 insertions(+), 11 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,10 @@ +2007-03-29 Bruno Haible <bruno@clisp.org> + + * lib/ldexpl.c: Include fpucw.h. + (ldexpl): Use BEGIN/END_LONG_DOUBLE_ROUNDING. Skip the last unneeded + multiplication. + * modules/ldexpl (Depends-on): Add fpucw. + 2007-03-29 Bruno Haible <bruno@clisp.org> * modules/ldexpl: New file.
--- a/lib/ldexpl.c +++ b/lib/ldexpl.c @@ -25,6 +25,7 @@ #include <math.h> #include <float.h> +#include "fpucw.h" #include "isnanl.h" long double @@ -32,22 +33,36 @@ { long double factor; int bit; + DECL_LONG_DOUBLE_ROUNDING + + BEGIN_LONG_DOUBLE_ROUNDING (); /* Check for zero, nan and infinity. */ - if (isnanl (x) || x + x == x) - return x; - - if (exp < 0) + if (!(isnanl (x) || x + x == x)) { - exp = -exp; - factor = 0.5L; + if (exp < 0) + { + exp = -exp; + factor = 0.5L; + } + else + factor = 2.0L; + + if (exp > 0) + for (bit = 1;;) + { + /* Invariant: Here bit = 2^i, factor = 2^-2^i or = 2^2^i, + and bit <= exp. */ + if (exp & bit) + x *= factor; + bit <<= 1; + if (bit > exp) + break; + factor = factor * factor; + } } - else - factor = 2.0L; - for (bit = 1; bit <= exp; bit <<= 1, factor *= factor) - if (exp & bit) - x *= factor; + END_LONG_DOUBLE_ROUNDING (); return x; }