Home | History | Annotate | Download | only in lib
      1 /* Split a double into fraction and mantissa, for hexadecimal printf.
      2    Copyright (C) 2007, 2009-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 #if ! defined USE_LONG_DOUBLE
     18 # include <config.h>
     19 #endif
     20 
     21 /* Specification.  */
     22 #ifdef USE_LONG_DOUBLE
     23 # include "printf-frexpl.h"
     24 #else
     25 # include "printf-frexp.h"
     26 #endif
     27 
     28 #include <float.h>
     29 #include <math.h>
     30 #ifdef USE_LONG_DOUBLE
     31 # include "fpucw.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 <= 2.0, not < 2.0.  */
     37 
     38 #ifdef USE_LONG_DOUBLE
     39 # define FUNC printf_frexpl
     40 # define DOUBLE long double
     41 # define MIN_EXP LDBL_MIN_EXP
     42 # if HAVE_FREXPL_IN_LIBC && HAVE_LDEXPL_IN_LIBC
     43 #  define USE_FREXP_LDEXP
     44 #  define FREXP frexpl
     45 #  define LDEXP ldexpl
     46 # endif
     47 # define DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING
     48 # define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING ()
     49 # define END_ROUNDING() END_LONG_DOUBLE_ROUNDING ()
     50 # define L_(literal) literal##L
     51 #else
     52 # define FUNC printf_frexp
     53 # define DOUBLE double
     54 # define MIN_EXP DBL_MIN_EXP
     55 # if HAVE_FREXP_IN_LIBC && HAVE_LDEXP_IN_LIBC
     56 #  define USE_FREXP_LDEXP
     57 #  define FREXP frexp
     58 #  define LDEXP ldexp
     59 # endif
     60 # define DECL_ROUNDING
     61 # define BEGIN_ROUNDING()
     62 # define END_ROUNDING()
     63 # define L_(literal) literal
     64 #endif
     65 
     66 DOUBLE
     67 FUNC (DOUBLE x, int *expptr)
     68 {
     69   int exponent;
     70   DECL_ROUNDING
     71 
     72   BEGIN_ROUNDING ();
     73 
     74 #ifdef USE_FREXP_LDEXP
     75   /* frexp and ldexp are usually faster than the loop below.  */
     76   x = FREXP (x, &exponent);
     77 
     78   x = x + x;
     79   exponent -= 1;
     80 
     81   if (exponent < MIN_EXP - 1)
     82     {
     83       x = LDEXP (x, exponent - (MIN_EXP - 1));
     84       exponent = MIN_EXP - 1;
     85     }
     86 #else
     87   {
     88     /* Since the exponent is an 'int', it fits in 64 bits.  Therefore the
     89        loops are executed no more than 64 times.  */
     90     DOUBLE pow2[64]; /* pow2[i] = 2^2^i */
     91     DOUBLE powh[64]; /* powh[i] = 2^-2^i */
     92     int i;
     93 
     94     exponent = 0;
     95     if (x >= L_(1.0))
     96       {
     97         /* A nonnegative exponent.  */
     98         {
     99           DOUBLE pow2_i; /* = pow2[i] */
    100           DOUBLE powh_i; /* = powh[i] */
    101 
    102           /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i,
    103              x * 2^exponent = argument, x >= 1.0.  */
    104           for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5);
    105                ;
    106                i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i)
    107             {
    108               if (x >= pow2_i)
    109                 {
    110                   exponent += (1 << i);
    111                   x *= powh_i;
    112                 }
    113               else
    114                 break;
    115 
    116               pow2[i] = pow2_i;
    117               powh[i] = powh_i;
    118             }
    119         }
    120         /* Here 1.0 <= x < 2^2^i.  */
    121       }
    122     else
    123       {
    124         /* A negative exponent.  */
    125         {
    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, exponent >= MIN_EXP - 1.  */
    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 (exponent - (1 << i) < MIN_EXP - 1)
    136                 break;
    137 
    138               exponent -= (1 << i);
    139               x *= pow2_i;
    140               if (x >= L_(1.0))
    141                 break;
    142 
    143               pow2[i] = pow2_i;
    144               powh[i] = powh_i;
    145             }
    146         }
    147         /* Here either x < 1.0 and exponent - 2^i < MIN_EXP - 1 <= exponent,
    148            or 1.0 <= x < 2^2^i and exponent >= MIN_EXP - 1.  */
    149 
    150         if (x < L_(1.0))
    151           /* Invariants: x * 2^exponent = argument, x < 1.0 and
    152              exponent - 2^i < MIN_EXP - 1 <= exponent.  */
    153           while (i > 0)
    154             {
    155               i--;
    156               if (exponent - (1 << i) >= MIN_EXP - 1)
    157                 {
    158                   exponent -= (1 << i);
    159                   x *= pow2[i];
    160                   if (x >= L_(1.0))
    161                     break;
    162                 }
    163             }
    164 
    165         /* Here either x < 1.0 and exponent = MIN_EXP - 1,
    166            or 1.0 <= x < 2^2^i and exponent >= MIN_EXP - 1.  */
    167       }
    168 
    169     /* Invariants: x * 2^exponent = argument, and
    170        either x < 1.0 and exponent = MIN_EXP - 1,
    171        or 1.0 <= x < 2^2^i and exponent >= MIN_EXP - 1.  */
    172     while (i > 0)
    173       {
    174         i--;
    175         if (x >= pow2[i])
    176           {
    177             exponent += (1 << i);
    178             x *= powh[i];
    179           }
    180       }
    181     /* Here either x < 1.0 and exponent = MIN_EXP - 1,
    182        or 1.0 <= x < 2.0 and exponent >= MIN_EXP - 1.  */
    183   }
    184 #endif
    185 
    186   END_ROUNDING ();
    187 
    188   *expptr = exponent;
    189   return x;
    190 }
    191