gl/frexp.c

   1 /* Split a double into fraction and mantissa.
   2    Copyright (C) 2007-2012 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 3 of the License, or
   7    (at your option) 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
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
  16
  17 /* Written by Paolo Bonzini <bonzini@gnu.org>, 2003, and
  18    Bruno Haible <bruno@clisp.org>, 2007.  */
  19
  20 #if ! defined USE_LONG_DOUBLE
  21 # include <config.h>
  22 #endif
  23
  24 /* Specification.  */
  25 #include <math.h>
  26
  27 #include <float.h>
  28 #ifdef USE_LONG_DOUBLE
  29 # include "isnanl-nolibm.h"
  30 # include "fpucw.h"
  31 #else
  32 # include "isnand-nolibm.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 DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING
  44 # define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING ()
  45 # define END_ROUNDING() END_LONG_DOUBLE_ROUNDING ()
  46 # define L_(literal) literal##L
  47 #else
  48 # define FUNC frexp
  49 # define DOUBLE double
  50 # define ISNAN isnand
  51 # define DECL_ROUNDING
  52 # define BEGIN_ROUNDING()
  53 # define END_ROUNDING()
  54 # define L_(literal) literal
  55 #endif
  56
  57 DOUBLE
  58 FUNC (DOUBLE x, int *expptr)
  59 {
  60   int sign;
  61   int exponent;
  62   DECL_ROUNDING
  63
  64   /* Test for NaN, infinity, and zero.  */
  65   if (ISNAN (x) || x + x == x)
  66     {
  67       *expptr = 0;
  68       return x;
  69     }
  70
  71   sign = 0;
  72   if (x < 0)
  73     {
  74       x = - x;
  75       sign = -1;
  76     }
  77
  78   BEGIN_ROUNDING ();
  79
  80   {
  81     /* Since the exponent is an 'int', it fits in 64 bits.  Therefore the
  82        loops are executed no more than 64 times.  */
  83     DOUBLE pow2[64]; /* pow2[i] = 2^2^i */
  84     DOUBLE powh[64]; /* powh[i] = 2^-2^i */
  85     int i;
  86
  87     exponent = 0;
  88     if (x >= L_(1.0))
  89       {
  90         /* A positive exponent.  */
  91         DOUBLE pow2_i; /* = pow2[i] */
  92         DOUBLE powh_i; /* = powh[i] */
  93
  94         /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i,
  95            x * 2^exponent = argument, x >= 1.0.  */
  96         for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5);
  97              ;
  98              i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i)
  99           {
 100             if (x >= pow2_i)
 101               {
 102                 exponent += (1 << i);
 103                 x *= powh_i;
 104               }
 105             else
 106               break;
 107
 108             pow2[i] = pow2_i;
 109             powh[i] = powh_i;
 110           }
 111         /* Avoid making x too small, as it could become a denormalized
 112            number and thus lose precision.  */
 113         while (i > 0 && x < pow2[i - 1])
 114           {
 115             i--;
 116             powh_i = powh[i];
 117           }
 118         exponent += (1 << i);
 119         x *= powh_i;
 120         /* Here 2^-2^i <= x < 1.0.  */
 121       }
 122     else
 123       {
 124         /* A negative or zero exponent.  */
 125         DOUBLE pow2_i; /* = pow2[i] */
 126         DOUBLE powh_i; /* = powh[i] */
 127
 128         /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i,
 129            x * 2^exponent = argument, x < 1.0.  */
 130         for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5);
 131              ;
 132              i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i)
 133           {
 134             if (x < powh_i)
 135               {
 136                 exponent -= (1 << i);
 137                 x *= pow2_i;
 138               }
 139             else
 140               break;
 141
 142             pow2[i] = pow2_i;
 143             powh[i] = powh_i;
 144           }
 145         /* Here 2^-2^i <= x < 1.0.  */
 146       }
 147
 148     /* Invariants: x * 2^exponent = argument, and 2^-2^i <= x < 1.0.  */
 149     while (i > 0)
 150       {
 151         i--;
 152         if (x < powh[i])
 153           {
 154             exponent -= (1 << i);
 155             x *= pow2[i];
 156           }
 157       }
 158     /* Here 0.5 <= x < 1.0.  */
 159   }
 160
 161   if (sign < 0)
 162     x = - x;
 163
 164   END_ROUNDING ();
 165
 166   *expptr = exponent;
 167   return x;
 168 }