New module 'frexp'.
authorBruno Haible <bruno@clisp.org>
Thu, 22 Mar 2007 03:59:38 +0000 (03:59 +0000)
committerBruno Haible <bruno@clisp.org>
Thu, 22 Mar 2007 03:59:38 +0000 (03:59 +0000)
ChangeLog
lib/frexp.c [new file with mode: 0644]
lib/math_.h
m4/frexp.m4 [new file with mode: 0644]
m4/math_h.m4
modules/frexp [new file with mode: 0644]
modules/math

index ac28398..ec09c26 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,15 @@
 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.
        * tests/test-isnanl.c: New file.
 
diff --git a/lib/frexp.c b/lib/frexp.c
new file mode 100644 (file)
index 0000000..b4b9f83
--- /dev/null
@@ -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
index be56101..18f3807 100644 (file)
@@ -30,6 +30,27 @@ extern "C" {
 #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
diff --git a/m4/frexp.m4 b/m4/frexp.m4
new file mode 100644 (file)
index 0000000..b0f0d73
--- /dev/null
@@ -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
+      ])
+    ])
+])
index 38cda1c..5b393f4 100644 (file)
@@ -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_MODULE_INDICATOR],
 
 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 @@ AC_DEFUN([gl_MATH_H_DEFAULTS],
   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])
 ])
diff --git a/modules/frexp b/modules/frexp
new file mode 100644 (file)
index 0000000..03407fd
--- /dev/null
@@ -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
+
index 6e175cb..7d265bd 100644 (file)
@@ -21,6 +21,7 @@ math.h: math_.h
        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 @@ math.h: math_.h
              -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