Mercurial > hg > octave-lojdl > gnulib-hg
changeset 8491:af98e0b44387
New module 'frexp'.
author | Bruno Haible <bruno@clisp.org> |
---|---|
date | Thu, 22 Mar 2007 03:59:38 +0000 |
parents | abafc0fdc7c6 |
children | 1dcdd2ec146b |
files | ChangeLog lib/frexp.c lib/math_.h m4/frexp.m4 m4/math_h.m4 modules/frexp modules/math |
diffstat | 7 files changed, 356 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,13 @@ +2007-03-21 Bruno Haible <bruno@clisp.org> + + * modules/frexp: New file. + * lib/frexp.c: New file. + * m4/frexp.m4: New file. + * lib/math_.h (frexp): New declaration. + * m4/math_h.m4 (gl_MATH_H_DEFAULTS): Also initialize GNULIB_FREXP and + REPLACE_FREXP. + * modules/math (math.h): Also substitute GNULIB_FREXP, REPLACE_FREXP. + 2007-03-21 Bruno Haible <bruno@clisp.org> * modules/isnanl-tests: New file.
new file mode 100644 --- /dev/null +++ b/lib/frexp.c @@ -0,0 +1,202 @@ +/* Split a double into fraction and mantissa. + Copyright (C) 2007 Free Software Foundation, Inc. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, + Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ + +#include <config.h> + +#if !(defined USE_LONG_DOUBLE && !HAVE_LONG_DOUBLE) + +/* Specification. */ +# include <math.h> + +# include <float.h> +# ifdef USE_LONG_DOUBLE +# include "isnanl-nolibm.h" +# else +# include "isnan.h" +# endif + +/* This file assumes FLT_RADIX = 2. If FLT_RADIX is a power of 2 greater + than 2, or not even a power of 2, some rounding errors can occur, so that + then the returned mantissa is only guaranteed to be <= 1.0, not < 1.0. */ + +# ifdef USE_LONG_DOUBLE +# define FUNC frexpl +# define DOUBLE long double +# define ISNAN isnanl +# define L_(literal) literal##L +# else +# define FUNC frexp +# define DOUBLE double +# define ISNAN isnan +# define L_(literal) literal +# endif + +DOUBLE +FUNC (DOUBLE x, int *exp) +{ + int sign; + int exponent; + + /* Test for NaN, infinity, and zero. */ + if (ISNAN (x) || x + x == x) + { + *exp = 0; + return x; + } + + sign = 0; + if (x < 0) + { + x = - x; + 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) + { + x /= *next; + exponent |= bit; + } + } + 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] */ + + /* 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. */ + } + + *exp = exponent; + return (sign < 0 ? - x : x); +} + +#else + +/* This declaration is solely to ensure that after preprocessing + this file is never empty. */ +typedef int dummy; + +#endif
--- a/lib/math_.h +++ b/lib/math_.h @@ -30,6 +30,27 @@ #endif +/* Write x as + x = mantissa * 2^exp + where + If x finite and nonzero: 0.5 <= |mantissa| < 1.0. + If x is zero: mantissa = x, exp = 0. + If x is infinite or NaN: mantissa = x, exp unspecified. + Store exp and return mantissa. */ +#if @GNULIB_FREXP@ +# if @REPLACE_FREXP@ +# define frexp rpl_frexp +extern double frexp (double x, int *exp); +# endif +#elif defined GNULIB_POSIXCHECK +# undef frexp +# define frexp(x,e) \ + (GL_LINK_WARNING ("frexp is unportable - " \ + "use gnulib module frexp for portability"), \ + frexp (x, e)) +#endif + + #if @GNULIB_MATHL@ || !@HAVE_DECL_ACOSL@ extern long double acosl (long double x); #endif
new file mode 100644 --- /dev/null +++ b/m4/frexp.m4 @@ -0,0 +1,92 @@ +# frexp.m4 serial 1 +dnl Copyright (C) 2007 Free Software Foundation, Inc. +dnl This file is free software; the Free Software Foundation +dnl gives unlimited permission to copy and/or distribute it, +dnl with or without modifications, as long as this notice is preserved. + +AC_DEFUN([gl_FUNC_FREXP], +[ + AC_REQUIRE([gl_MATH_H_DEFAULTS]) + FREXP_LIBM= + AC_CACHE_CHECK([whether frexp() can be used without linking with libm], + [gl_cv_func_frexp_no_libm], + [ + AC_TRY_LINK([#include <math.h> + long double x;], + [int e; return frexp (x, &e) > 0;], + [gl_cv_func_frexp_no_libm=yes], + [gl_cv_func_frexp_no_libm=no]) + ]) + if test $gl_cv_func_frexp_no_libm = no; then + AC_CACHE_CHECK([whether frexp() can be used with libm], + [gl_cv_func_frexp_in_libm], + [ + save_LIBS="$LIBS" + LIBS="$LIBS -lm" + AC_TRY_LINK([#include <math.h> + long double x;], + [int e; return frexp (x, &e) > 0;], + [gl_cv_func_frexp_in_libm=yes], + [gl_cv_func_frexp_in_libm=no]) + LIBS="$save_LIBS" + ]) + if test $gl_cv_func_frexp_in_libm = yes; then + FREXP_LIBM=-lm + fi + fi + if test $gl_cv_func_frexp_no_libm = yes \ + || test $gl_cv_func_frexp_in_libm = yes; then + save_LIBS="$LIBS" + LIBS="$LIBS $FREXP_LIBM" + gl_FUNC_FREXP_WORKS + LIBS="$save_LIBS" + case "$gl_cv_func_frexp_works" in + *yes) gl_func_frexp=yes ;; + *) gl_func_frexp=no; REPLACE_FREXP=1; FREXP_LIBM= ;; + esac + else + gl_func_frexp=no + fi + if test $gl_func_frexp = yes; then + AC_DEFINE([HAVE_FREXP], 1, + [Define if the frexp() function is available and works.]) + else + AC_LIBOBJ([frexp]) + fi +]) + +dnl Test whether frexp() works also on denormalized numbers. +dnl This test fails e.g. on NetBSD. +AC_DEFUN([gl_FUNC_FREXP_WORKS], +[ + AC_REQUIRE([AC_PROG_CC]) + AC_REQUIRE([AC_CANONICAL_HOST]) dnl for cross-compiles + AC_CACHE_CHECK([whether frexp works], [gl_cv_func_frexp_works], + [ + AC_TRY_RUN([ +#include <float.h> +#include <math.h> +int main() +{ + int i; + volatile double x; + for (i = 1, x = 1.0; i >= DBL_MIN_EXP; i--, x *= 0.5) + ; + if (x > 0.0) + { + int exp; + double y = frexp (x, &exp); + /* On machines with IEEE754 arithmetic: x = 1.11254e-308, exp = -1022. + On NetBSD: y = 0.75. Correct: y = 0.5. */ + if (y != 0.5) + return 1; + } + return 0; +}], [gl_cv_func_frexp_works=yes], [gl_cv_func_frexp_works=no], + [case "$host_os" in + netbsd*) gl_cv_func_frexp_works="guessing no";; + *) gl_cv_func_frexp_works="guessing yes";; + esac + ]) + ]) +])
--- a/m4/math_h.m4 +++ b/m4/math_h.m4 @@ -1,4 +1,4 @@ -# math_h.m4 serial 2 +# math_h.m4 serial 3 dnl Copyright (C) 2007 Free Software Foundation, Inc. dnl This file is free software; the Free Software Foundation dnl gives unlimited permission to copy and/or distribute it, @@ -21,6 +21,7 @@ AC_DEFUN([gl_MATH_H_DEFAULTS], [ + GNULIB_FREXP=0; AC_SUBST([GNULIB_FREXP]) GNULIB_MATHL=0; AC_SUBST([GNULIB_MATHL]) dnl Assume proper GNU behavior unless another module says otherwise. HAVE_DECL_ACOSL=1; AC_SUBST([HAVE_DECL_ACOSL]) @@ -36,4 +37,5 @@ HAVE_DECL_SINL=1; AC_SUBST([HAVE_DECL_SINL]) HAVE_DECL_SQRTL=1; AC_SUBST([HAVE_DECL_SQRTL]) HAVE_DECL_TANL=1; AC_SUBST([HAVE_DECL_TANL]) + REPLACE_FREXP=0; AC_SUBST([REPLACE_FREXP]) ])
new file mode 100644 --- /dev/null +++ b/modules/frexp @@ -0,0 +1,26 @@ +Description: +frexp() function: split a double into its constituents. + +Files: +lib/frexp.c +m4/frexp.m4 + +Depends-on: +math +isnan-nolibm + +configure.ac: +gl_FUNC_FREXP +gl_MATH_MODULE_INDICATOR([frexp]) + +Makefile.am: + +Include: +<math.h> + +License: +LGPL + +Maintainer: +Bruno Haible, Paolo Bonzini +
--- a/modules/math +++ b/modules/math @@ -21,6 +21,7 @@ rm -f $@-t $@ { echo '/* DO NOT EDIT! GENERATED AUTOMATICALLY! */' && \ sed -e 's|@''ABSOLUTE_MATH_H''@|$(ABSOLUTE_MATH_H)|g' \ + -e 's|@''GNULIB_FREXP''@|$(GNULIB_FREXP)|g' \ -e 's|@''GNULIB_MATHL''@|$(GNULIB_MATHL)|g' \ -e 's|@''HAVE_DECL_ACOSL''@|$(HAVE_DECL_ACOSL)|g' \ -e 's|@''HAVE_DECL_ASINL''@|$(HAVE_DECL_ASINL)|g' \ @@ -35,6 +36,7 @@ -e 's|@''HAVE_DECL_SINL''@|$(HAVE_DECL_SINL)|g' \ -e 's|@''HAVE_DECL_SQRTL''@|$(HAVE_DECL_SQRTL)|g' \ -e 's|@''HAVE_DECL_TANL''@|$(HAVE_DECL_TANL)|g' \ + -e 's|@''REPLACE_FREXP''@|$(REPLACE_FREXP)|g' \ -e '/definition of GL_LINK_WARNING/r $(LINK_WARNING_H)' \ < $(srcdir)/math_.h; \ } > $@-t