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