bug-gnulib
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

new module 'frexp'


From: Bruno Haible
Subject: new module 'frexp'
Date: Thu, 22 Mar 2007 05:01:05 +0100
User-agent: KMail/1.5.4

NetBSD does not only have a broken isnanl() function. It also has a broken
frexp(): it does not treat denormalized numbers correctly.

Here is a replacement module for frexp(). Paolo, I tried using your
algorithm from frexpl.c, but
  - I cannot enable it since it is GPL and I want an LGPL module (because
    vasnprintf uses it, and vasnprintf is LGPL),
  - it produces wrong results for x between 2^(DBL_MAX_EXP-1) and 2^DBL_MAX_EXP.
    You can debug it by changing the 'if (0)' to 'if (1)' and running the
    unit test.

2007-03-21  Bruno Haible  <address@hidden>

        * 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.

============================= modules/frexp ===============================
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

============================== lib/frexp.c ================================
/* 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
============================== m4/frexp.m4 ================================
# 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
      ])
    ])
])
===========================================================================
*** lib/math_.h 7 Mar 2007 01:12:01 -0000       1.3
--- lib/math_.h 22 Mar 2007 03:56:01 -0000
***************
*** 30,35 ****
--- 30,56 ----
  #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@ || address@hidden@
  extern long double acosl (long double x);
  #endif
*** m4/math_h.m4        7 Mar 2007 01:12:01 -0000       1.2
--- m4/math_h.m4        22 Mar 2007 03:56:01 -0000
***************
*** 1,4 ****
! # math_h.m4 serial 2
  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,
--- 1,4 ----
! # 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,26 ****
--- 21,27 ----
  
  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,39 ****
--- 37,41 ----
    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])
  ])
*** modules/math        7 Mar 2007 01:12:01 -0000       1.2
--- modules/math        22 Mar 2007 03:56:01 -0000
***************
*** 21,26 ****
--- 21,27 ----
        rm -f address@hidden $@
        { 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,40 ****
--- 36,42 ----
              -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; \
        } > address@hidden





reply via email to

[Prev in Thread] Current Thread [Next in Thread]