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 # else
  32 #  include "isnan.h"
  33 # endif
  34
  35 /* This file assumes FLT_RADIX = 2.  If FLT_RADIX is a power of 2 greater
  36    than 2, or not even a power of 2, some rounding errors can occur, so that
  37    then the returned mantissa is only guaranteed to be <= 1.0, not < 1.0.  */
  38
  39 # ifdef USE_LONG_DOUBLE
  40 #  define FUNC frexpl
  41 #  define DOUBLE long double
  42 #  define ISNAN isnanl
  43 #  define L_(literal) literal##L
  44 # else
  45 #  define FUNC frexp
  46 #  define DOUBLE double
  47 #  define ISNAN isnan
  48 #  define L_(literal) literal
  49 # endif
  50
  51 DOUBLE
  52 FUNC (DOUBLE x, int *exp)
  53 {
  54   int sign;
  55   int exponent;
  56
  57   /* Test for NaN, infinity, and zero.  */
  58   if (ISNAN (x) || x + x == x)
  59     {
  60       *exp = 0;
  61       return x;
  62     }
  63
  64   sign = 0;
  65   if (x < 0)
  66     {
  67       x = - x;
  68       sign = -1;
  69     }
  70
  71   {
  72     /* Since the exponent is an 'int', it fits in 64 bits.  Therefore the
  73        loops are executed no more than 64 times.  */
  74     DOUBLE pow2[64]; /* pow2[i] = 2^2^i */
  75     DOUBLE powh[64]; /* powh[i] = 2^-2^i */
  76     int i;
  77
  78     exponent = 0;
  79     if (x >= L_(1.0))
  80       {
  81         /* A positive exponent.  */
  82         DOUBLE pow2_i; /* = pow2[i] */
  83         DOUBLE powh_i; /* = powh[i] */
  84
  85         /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i,
  86            x * 2^exponent = argument, x >= 1.0.  */
  87         for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5);
  88              ;
  89              i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i)
  90           {
  91             if (x >= pow2_i)
  92               {
  93                 exponent += (1 << i);
  94                 x *= powh_i;
  95               }
  96             else
  97               break;
  98
  99             pow2[i] = pow2_i;
 100             powh[i] = powh_i;
 101           }
 102         /* Avoid making x too small, as it could become a denormalized
 103            number and thus lose precision.  */
 104         while (i > 0 && x < pow2[i - 1])
 105           {
 106             i--;
 107             powh_i = powh[i];
 108           }
 109         exponent += (1 << i);
 110         x *= powh_i;
 111         /* Here 2^-2^i <= x < 1.0.  */
 112       }
 113     else
 114       {
 115         /* A negative or zero exponent.  */
 116         DOUBLE pow2_i; /* = pow2[i] */
 117         DOUBLE powh_i; /* = powh[i] */
 118
 119         /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i,
 120            x * 2^exponent = argument, x < 1.0.  */
 121         for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5);
 122              ;
 123              i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i)
 124           {
 125             if (x < powh_i)
 126               {
 127                 exponent -= (1 << i);
 128                 x *= pow2_i;
 129               }
 130             else
 131               break;
 132
 133             pow2[i] = pow2_i;
 134             powh[i] = powh_i;
 135           }
 136         /* Here 2^-2^i <= x < 1.0.  */
 137       }
 138
 139     /* Invariants: x * 2^exponent = argument, and 2^-2^i <= x < 1.0.  */
 140     while (i > 0)
 141       {
 142         i--;
 143         if (x < powh[i])
 144           {
 145             exponent -= (1 << i);
 146             x *= pow2[i];
 147           }
 148       }
 149     /* Here 0.5 <= x < 1.0.  */
 150   }
 151
 152   *exp = exponent;
 153   return (sign < 0 ? - x : x);
 154 }
 155
 156 #else
 157
 158 /* This declaration is solely to ensure that after preprocessing
 159    this file is never empty.  */
 160 typedef int dummy;
 161
 162 #endif