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 <[EMAIL PROTECTED]>
* 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@ || [EMAIL PROTECTED]@
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 [EMAIL PROTECTED] $@
{ 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; \
} > [EMAIL PROTECTED]
