lib/frexp.c

   1 /* Split a double into fraction and mantissa.
   2    Copyright (C) 2007 Free Software Foundation, Inc.
   3
   4    This program is free software; you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 2, or (at your option)
   7    any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License along
  15    with this program; if not, write to the Free Software Foundation,
  16    Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
  17
  18 /* Written by Paolo Bonzini <bonzini@gnu.org>, 2003, and
  19    Bruno Haible <bruno@clisp.org>, 2007.  */
  20
  21 #include <config.h>
  22
  23 #if !(defined USE_LONG_DOUBLE && !HAVE_LONG_DOUBLE)
  24
  25 /* Specification.  */
  26 # include <math.h>
  27
  28 # include <float.h>
  29 # ifdef USE_LONG_DOUBLE
  30 #  include "isnanl-nolibm.h"
  31 #  include "fpucw.h"
  32 # else
  33 #  include "isnan.h"
  34 # endif
  35
  36 /* This file assumes FLT_RADIX = 2.  If FLT_RADIX is a power of 2 greater
  37    than 2, or not even a power of 2, some rounding errors can occur, so that
  38    then the returned mantissa is only guaranteed to be <= 1.0, not < 1.0.  */
  39
  40 # ifdef USE_LONG_DOUBLE
  41 #  define FUNC frexpl
  42 #  define DOUBLE long double
  43 #  define ISNAN isnanl
  44 #  define DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING
  45 #  define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING ()
  46 #  define END_ROUNDING() END_LONG_DOUBLE_ROUNDING ()
  47 #  define L_(literal) literal##L
  48 # else
  49 #  define FUNC frexp
  50 #  define DOUBLE double
  51 #  define ISNAN isnan
  52 #  define DECL_ROUNDING
  53 #  define BEGIN_ROUNDING()
  54 #  define END_ROUNDING()
  55 #  define L_(literal) literal
  56 # endif
  57
  58 DOUBLE
  59 FUNC (DOUBLE x, int *exp)
  60 {
  61   int sign;
  62   int exponent;
  63   DECL_ROUNDING
  64
  65   /* Test for NaN, infinity, and zero.  */
  66   if (ISNAN (x) || x + x == x)
  67     {
  68       *exp = 0;
  69       return x;
  70     }
  71
  72   sign = 0;
  73   if (x < 0)
  74     {
  75       x = - x;
  76       sign = -1;
  77     }
  78
  79   BEGIN_ROUNDING ();
  80
  81   {
  82     /* Since the exponent is an 'int', it fits in 64 bits.  Therefore the
  83        loops are executed no more than 64 times.  */
  84     DOUBLE pow2[64]; /* pow2[i] = 2^2^i */
  85     DOUBLE powh[64]; /* powh[i] = 2^-2^i */
  86     int i;
  87
  88     exponent = 0;
  89     if (x >= L_(1.0))
  90       {
  91         /* A positive exponent.  */
  92         DOUBLE pow2_i; /* = pow2[i] */
  93         DOUBLE powh_i; /* = powh[i] */
  94
  95         /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i,
  96            x * 2^exponent = argument, x >= 1.0.  */
  97         for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5);
  98              ;
  99              i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i)
 100           {
 101             if (x >= pow2_i)
 102               {
 103                 exponent += (1 << i);
 104                 x *= powh_i;
 105               }
 106             else
 107               break;
 108
 109             pow2[i] = pow2_i;
 110             powh[i] = powh_i;
 111           }
 112         /* Avoid making x too small, as it could become a denormalized
 113            number and thus lose precision.  */
 114         while (i > 0 && x < pow2[i - 1])
 115           {
 116             i--;
 117             powh_i = powh[i];
 118           }
 119         exponent += (1 << i);
 120         x *= powh_i;
 121         /* Here 2^-2^i <= x < 1.0.  */
 122       }
 123     else
 124       {
 125         /* A negative or zero exponent.  */
 126         DOUBLE pow2_i; /* = pow2[i] */
 127         DOUBLE powh_i; /* = powh[i] */
 128
 129         /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i,
 130            x * 2^exponent = argument, x < 1.0.  */
 131         for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5);
 132              ;
 133              i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i)
 134           {
 135             if (x < powh_i)
 136               {
 137                 exponent -= (1 << i);
 138                 x *= pow2_i;
 139               }
 140             else
 141               break;
 142
 143             pow2[i] = pow2_i;
 144             powh[i] = powh_i;
 145           }
 146         /* Here 2^-2^i <= x < 1.0.  */
 147       }
 148
 149     /* Invariants: x * 2^exponent = argument, and 2^-2^i <= x < 1.0.  */
 150     while (i > 0)
 151       {
 152         i--;
 153         if (x < powh[i])
 154           {
 155             exponent -= (1 << i);
 156             x *= pow2[i];
 157           }
 158       }
 159     /* Here 0.5 <= x < 1.0.  */
 160   }
 161
 162   if (sign < 0)
 163     x = - x;
 164
 165   END_ROUNDING ();
 166
 167   *exp = exponent;
 168   return x;
 169 }
 170
 171 #else
 172
 173 /* This declaration is solely to ensure that after preprocessing
 174    this file is never empty.  */
 175 typedef int dummy;
 176
 177 #endif