new module printf-frexp

[Bruno Haible]

This module is also for the 'a' and 'A' conversion of printf: Splitting a
'double' into mantissa and exponent, without requiring libm.

This module assumes that the radix of floating-point numbers is 2. A safe
assumption, I hope, 18 years after IEEE 754.

2007-02-25  Bruno Haible  <address@hidden>

        * modules/printf-frexp: New file.
        * lib/printf-frexp.h: New file.
        * lib/printf-frexp.c: New file.
        * m4/printf-frexp.m4: New file.

============================== modules/printf-frexp ==========================
Description:
printf_frexp() function: split a double into fraction and mantissa, for
hexadecimal printf, without requiring libm.

Files:
lib/printf-frexp.h
lib/printf-frexp.c
m4/printf-frexp.m4

Depends-on:
verify

configure.ac:
gl_FUNC_PRINTF_FREXP

Makefile.am:
lib_SOURCES += printf-frexp.c

Include:
#include "printf-frexp.h"

License:
LGPL

Maintainer:
Bruno Haible

=============================== lib/printf-frexp.h ===========================
/* Split a double into fraction and mantissa, for hexadecimal printf.
   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.  */

/* Write a finite, positive number x as
     x = mantissa * 2^exp
   where exp >= DBL_MIN_EXP - 1,
         mantissa < 2.0,
         if x is not a denormalized number then mantissa >= 1.0.
   Store exp and return mantissa.  */
extern double printf_frexp (double x, int *exp);
=============================== lib/printf-frexp.c ===========================
/* Split a double into fraction and mantissa, for hexadecimal printf.
   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.  */
# ifdef USE_LONG_DOUBLE
#  include "printf-frexpl.h"
# else
#  include "printf-frexp.h"
# endif

# include <float.h>
# include <math.h>
# include "verify.h"

/* This file assumes FLT_RADIX is 2.  */
verify (FLT_RADIX == 2);

# ifdef USE_LONG_DOUBLE
#  define FUNC printf_frexpl
#  define DOUBLE long double
#  define MIN_EXP LDBL_MIN_EXP
#  if HAVE_FREXPL_IN_LIBC && HAVE_LDEXPL_IN_LIBC
#   define USE_FREXP_LDEXP
#   define FREXP frexpl
#   define LDEXP ldexpl
#  endif
#  define L_(literal) literal##L
# else
#  define FUNC printf_frexp
#  define DOUBLE double
#  define MIN_EXP DBL_MIN_EXP
#  if HAVE_FREXP_IN_LIBC && HAVE_LDEXP_IN_LIBC
#   define USE_FREXP_LDEXP
#   define FREXP frexp
#   define LDEXP ldexp
#  endif
#  define L_(literal) literal
# endif

DOUBLE
FUNC (DOUBLE x, int *exp)
{
  int exponent;

# ifdef USE_FREXP_LDEXP
  /* frexp and ldexp are usually faster than the loop below.  */
  x = FREXP (x, &exponent);

  x = x + x;
  exponent -= 1;

  if (exponent < MIN_EXP - 1)
    {
      x = LDEXP (x, exponent - (MIN_EXP - 1));
      exponent = MIN_EXP - 1;
    }
# else
  /* 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 nonnegative 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;
          }
      }
      /* Here 1.0 <= x < 2^2^i.  */
    }
  else
    {
      /* A negative 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, exponent >= MIN_EXP - 1.  */
        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 (exponent - (1 << i) < MIN_EXP - 1)
              break;

            exponent -= (1 << i);
            x *= pow2_i;
            if (x >= L_(1.0))
              break;

            pow2[i] = pow2_i;
            powh[i] = powh_i;
          }
      }
      /* Here either x < 1.0 and exponent - 2^i < MIN_EXP - 1 <= exponent,
         or 1.0 <= x < 2^2^i and exponent >= MIN_EXP - 1.  */

      if (x < L_(1.0))
        /* Invariants: x * 2^exponent = argument, x < 1.0 and
           exponent - 2^i < MIN_EXP - 1 <= exponent.  */
        while (i > 0)
          {
            i--;
            if (exponent - (1 << i) >= MIN_EXP - 1)
              {
                exponent -= (1 << i);
                x *= pow2[i];
                if (x >= L_(1.0))
                  break;
              }
          }

      /* Here either x < 1.0 and exponent = MIN_EXP - 1,
         or 1.0 <= x < 2^2^i and exponent >= MIN_EXP - 1.  */
    }

  /* Invariants: x * 2^exponent = argument, and
     either x < 1.0 and exponent = MIN_EXP - 1,
     or 1.0 <= x < 2^2^i and exponent >= MIN_EXP - 1.  */
  while (i > 0)
    {
      i--;
      if (x >= pow2[i])
        {
          exponent += (1 << i);
          x *= powh[i];
        }
    }
  /* Here either x < 1.0 and exponent = MIN_EXP - 1,
     or 1.0 <= x < 2.0 and exponent >= MIN_EXP - 1.  */
# endif

  *exp = exponent;
  return x;
}

#else

/* This declaration is solely to ensure that after preprocessing
   this file is never empty.  */
typedef int dummy;

#endif
=============================== m4/printf-frexp.m4 ===========================
# printf-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.

dnl Check how to define printf_frexp() without linking with libm.

AC_DEFUN([gl_FUNC_PRINTF_FREXP],
[
  AC_CACHE_CHECK([whether frexp can be used without linking with libm],
    [gl_cv_func_frexp_no_libm],
    [
      AC_TRY_LINK([#include <math.h>
                   double x;
                   int y;],
                  [return frexp (x, &y) < 1;],
        [gl_cv_func_frexp_no_libm=yes],
        [gl_cv_func_frexp_no_libm=no])
    ])
  if test $gl_cv_func_frexp_no_libm = yes; then
    AC_DEFINE([HAVE_FREXP_IN_LIBC], 1,
      [Define if the frexp function is available in libc.])
  fi

  AC_CACHE_CHECK([whether ldexp can be used without linking with libm],
    [gl_cv_func_ldexp_no_libm],
    [
      AC_TRY_LINK([#include <math.h>
                   double x;
                   int y;],
                  [return ldexp (x, y) < 1;],
        [gl_cv_func_ldexp_no_libm=yes],
        [gl_cv_func_ldexp_no_libm=no])
    ])
  if test $gl_cv_func_ldexp_no_libm = yes; then
    AC_DEFINE([HAVE_LDEXP_IN_LIBC], 1,
      [Define if the ldexp function is available in libc.])
  fi
])