gdtoa/FreeBSD/_hdtoa.c

   1 /*-
   2  * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice, this list of conditions and the following disclaimer.
  10  * 2. Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  *
  14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  15  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  16  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  17  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  18  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  19  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  20  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  21  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  22  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  23  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  24  * SUCH DAMAGE.
  25  */
  26
  27 #include <sys/cdefs.h>
  28 __FBSDID("$FreeBSD: src/lib/libc/gdtoa/_hdtoa.c,v 1.2 2004/01/21 04:51:50 grehan Exp $");
  29
  30 #include <float.h>
  31 #include <inttypes.h>
  32 #include <limits.h>
  33 #include <math.h>
  34 #include <stdlib.h>
  35 #include "fpmath.h"
  36 #include "gdtoaimp.h"
  37
  38 /* Strings values used by dtoa() */
  39 #define INFSTR  "Infinity"
  40 #define NANSTR  "NaN"
  41
  42 #define DBL_BIAS        (DBL_MAX_EXP - 1)
  43 #define LDBL_BIAS       (LDBL_MAX_EXP - 1)
  44
  45 #ifdef  LDBL_IMPLICIT_NBIT
  46 #define LDBL_NBIT_ADJ   0
  47 #else
  48 #define LDBL_NBIT_ADJ   1
  49 #endif
  50
  51 /*
  52  * Efficiently compute the log2 of an integer.  Uses a combination of
  53  * arcane tricks found in fortune and arcane tricks not (yet) in
  54  * fortune.  This routine behaves similarly to fls(9).
  55  */
  56 static int
  57 log2_32(uint32_t n)
  58 {
  59
  60         n |= (n >> 1);
  61         n |= (n >> 2);
  62         n |= (n >> 4);
  63         n |= (n >> 8);
  64         n |= (n >> 16);
  65
  66         n = (n & 0x55555555) + ((n & 0xaaaaaaaa) >> 1);
  67         n = (n & 0x33333333) + ((n & 0xcccccccc) >> 2);
  68         n = (n & 0x0f0f0f0f) + ((n & 0xf0f0f0f0) >> 4);
  69         n = (n & 0x00ff00ff) + ((n & 0xff00ff00) >> 8);
  70         n = (n & 0x0000ffff) + ((n & 0xffff0000) >> 16);
  71         return (n - 1);
  72 }
  73
  74 #if (LDBL_MANH_SIZE > 32 || LDBL_MANL_SIZE > 32)
  75
  76 static int
  77 log2_64(uint64_t n)
  78 {
  79
  80         if (n >> 32 != 0)
  81                 return (log2_32((uint32_t)(n >> 32)) + 32);
  82         else
  83                 return (log2_32((uint32_t)n));
  84 }
  85
  86 #endif  /* (LDBL_MANH_SIZE > 32 || LDBL_MANL_SIZE > 32) */
  87
  88 /*
  89  * Round up the given digit string.  If the digit string is fff...f,
  90  * this procedure sets it to 100...0 and returns 1 to indicate that
  91  * the exponent needs to be bumped.  Otherwise, 0 is returned.
  92  */
  93 static int
  94 roundup(char *s0, int ndigits)
  95 {
  96         char *s;
  97
  98         for (s = s0 + ndigits - 1; *s == 0xf; s--) {
  99                 if (s == s0) {
 100                         *s = 1;
 101                         return (1);
 102                 }
 103                 ++*s;
 104         }
 105         ++*s;
 106         return (0);
 107 }
 108
 109 /*
 110  * Round the given digit string to ndigits digits according to the
 111  * current rounding mode.  Note that this could produce a string whose
 112  * value is not representable in the corresponding floating-point
 113  * type.  The exponent pointed to by decpt is adjusted if necessary.
 114  */
 115 static void
 116 dorounding(char *s0, int ndigits, int sign, int *decpt)
 117 {
 118         int adjust = 0; /* do we need to adjust the exponent? */
 119
 120         switch (FLT_ROUNDS) {
 121         case 0:         /* toward zero */
 122         default:        /* implementation-defined */
 123                 break;
 124         case 1:         /* to nearest, halfway rounds to even */
 125                 if ((s0[ndigits] > 8) ||
 126                     (s0[ndigits] == 8 && s0[ndigits - 1] & 1))
 127                         adjust = roundup(s0, ndigits);
 128                 break;
 129         case 2:         /* toward +inf */
 130                 if (sign == 0)
 131                         adjust = roundup(s0, ndigits);
 132                 break;
 133         case 3:         /* toward -inf */
 134                 if (sign != 0)
 135                         adjust = roundup(s0, ndigits);
 136                 break;
 137         }
 138
 139         if (adjust)
 140                 *decpt += 4;
 141 }
 142
 143 /*
 144  * This procedure converts a double-precision number in IEEE format
 145  * into a string of hexadecimal digits and an exponent of 2.  Its
 146  * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
 147  * following exceptions:
 148  *
 149  * - An ndigits < 0 causes it to use as many digits as necessary to
 150  *   represent the number exactly.
 151  * - The additional xdigs argument should point to either the string
 152  *   "0123456789ABCDEF" or the string "0123456789abcdef", depending on
 153  *   which case is desired.
 154  * - This routine does not repeat dtoa's mistake of setting decpt
 155  *   to 9999 in the case of an infinity or NaN.  INT_MAX is used
 156  *   for this purpose instead.
 157  *
 158  * Note that the C99 standard does not specify what the leading digit
 159  * should be for non-zero numbers.  For instance, 0x1.3p3 is the same
 160  * as 0x2.6p2 is the same as 0x4.cp3.  This implementation chooses the
 161  * first digit so that subsequent digits are aligned on nibble
 162  * boundaries (before rounding).
 163  *
 164  * Inputs:      d, xdigs, ndigits
 165  * Outputs:     decpt, sign, rve
 166  */
 167 char *
 168 __hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
 169     char **rve)
 170 {
 171         union IEEEd2bits u;
 172         char *s, *s0;
 173         int bufsize;
 174         int impnbit;    /* implicit normalization bit */
 175         int pos;
 176         int shift;      /* for subnormals, # of shifts required to normalize */
 177         int sigfigs;    /* number of significant hex figures in result */
 178
 179         u.d = d;
 180         *sign = u.bits.sign;
 181
 182         switch (fpclassify(d)) {
 183         case FP_NORMAL:
 184                 sigfigs = (DBL_MANT_DIG + 3) / 4;
 185                 impnbit = 1 << ((DBL_MANT_DIG - 1) % 4);
 186                 *decpt = u.bits.exp - DBL_BIAS + 1 -
 187                     ((DBL_MANT_DIG - 1) % 4);
 188                 break;
 189         case FP_ZERO:
 190                 *decpt = 1;
 191                 return (nrv_alloc("0", rve, 1));
 192         case FP_SUBNORMAL:
 193                 /*
 194                  * The position of the highest-order bit tells us by
 195                  * how much to adjust the exponent (decpt).  The
 196                  * adjustment is raised to the next nibble boundary
 197                  * since we will later choose the leftmost hexadecimal
 198                  * digit so that all subsequent digits align on nibble
 199                  * boundaries.
 200                  */
 201                 if (u.bits.manh != 0) {
 202                         pos = log2_32(u.bits.manh);
 203                         shift = DBL_MANH_SIZE - pos;
 204                 } else {
 205                         pos = log2_32(u.bits.manl);
 206                         shift = DBL_MANH_SIZE + DBL_MANL_SIZE - pos;
 207                 }
 208                 sigfigs = (3 + DBL_MANT_DIG - shift) / 4;
 209                 impnbit = 0;
 210                 *decpt = DBL_MIN_EXP - ((shift + 3) & ~(4 - 1));
 211                 break;
 212         case FP_INFINITE:
 213                 *decpt = INT_MAX;
 214                 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
 215         case FP_NAN:
 216                 *decpt = INT_MAX;
 217                 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
 218         default:
 219                 abort();
 220         }
 221
 222         /* FP_NORMAL or FP_SUBNORMAL */
 223
 224         if (ndigits == 0)               /* dtoa() compatibility */
 225                 ndigits = 1;
 226
 227         /*
 228          * For simplicity, we generate all the digits even if the
 229          * caller has requested fewer.
 230          */
 231         bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
 232         s0 = rv_alloc(bufsize);
 233
 234         /*
 235          * We work from right to left, first adding any requested zero
 236          * padding, then the least significant portion of the
 237          * mantissa, followed by the most significant.  The buffer is
 238          * filled with the byte values 0x0 through 0xf, which are
 239          * converted to xdigs[0x0] through xdigs[0xf] after the
 240          * rounding phase.
 241          */
 242         for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
 243                 *s = 0;
 244         for (; s > s0 + sigfigs - (DBL_MANL_SIZE / 4) - 1 && s > s0; s--) {
 245                 *s = u.bits.manl & 0xf;
 246                 u.bits.manl >>= 4;
 247         }
 248         for (; s > s0; s--) {
 249                 *s = u.bits.manh & 0xf;
 250                 u.bits.manh >>= 4;
 251         }
 252
 253         /*
 254          * At this point, we have snarfed all the bits in the
 255          * mantissa, with the possible exception of the highest-order
 256          * (partial) nibble, which is dealt with by the next
 257          * statement.  That nibble is usually in manh, but it could be
 258          * in manl instead for small subnormals.  We also tack on the
 259          * implicit normalization bit if appropriate.
 260          */
 261         *s = u.bits.manh | u.bits.manl | impnbit;
 262
 263         /* If ndigits < 0, we are expected to auto-size the precision. */
 264         if (ndigits < 0) {
 265                 for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
 266                         ;
 267         }
 268
 269         if (sigfigs > ndigits && s0[ndigits] != 0)
 270                 dorounding(s0, ndigits, u.bits.sign, decpt);
 271
 272         s = s0 + ndigits;
 273         if (rve != NULL)
 274                 *rve = s;
 275         *s-- = '\0';
 276         for (; s >= s0; s--)
 277                 *s = xdigs[(unsigned int)*s];
 278
 279         return (s0);
 280 }
 281
 282 #if (LDBL_MANT_DIG > DBL_MANT_DIG)
 283
 284 /*
 285  * This is the long double version of __hdtoa().
 286  *
 287  * On architectures that have an explicit integer bit, unnormals and
 288  * pseudo-denormals cause problems in the conversion routine, so they
 289  * are ``fixed'' by effectively toggling the integer bit.  Although
 290  * this is not correct behavior, the hardware will not produce these
 291  * formats externally.
 292  */
 293 char *
 294 __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
 295     char **rve)
 296 {
 297         union IEEEl2bits u;
 298         char *s, *s0;
 299         int bufsize;
 300         int impnbit;    /* implicit normalization bit */
 301         int pos;
 302         int shift;      /* for subnormals, # of shifts required to normalize */
 303         int sigfigs;    /* number of significant hex figures in result */
 304
 305         u.e = e;
 306         *sign = u.bits.sign;
 307
 308         switch (fpclassify(e)) {
 309         case FP_NORMAL:
 310                 sigfigs = (LDBL_MANT_DIG + 3) / 4;
 311                 impnbit = 1 << ((LDBL_MANT_DIG - 1) % 4);
 312                 *decpt = u.bits.exp - LDBL_BIAS + 1 -
 313                     ((LDBL_MANT_DIG - 1) % 4);
 314                 break;
 315         case FP_ZERO:
 316                 *decpt = 1;
 317                 return (nrv_alloc("0", rve, 1));
 318         case FP_SUBNORMAL:
 319                 /*
 320                  * The position of the highest-order bit tells us by
 321                  * how much to adjust the exponent (decpt).  The
 322                  * adjustment is raised to the next nibble boundary
 323                  * since we will later choose the leftmost hexadecimal
 324                  * digit so that all subsequent digits align on nibble
 325                  * boundaries.
 326                  */
 327 #ifdef  LDBL_IMPLICIT_NBIT
 328                 /* Don't trust the normalization bit to be off. */
 329                 u.bits.manh &= ~(~0ULL << (LDBL_MANH_SIZE - 1));
 330 #endif
 331                 if (u.bits.manh != 0) {
 332 #if LDBL_MANH_SIZE > 32
 333                         pos = log2_64(u.bits.manh);
 334 #else
 335                         pos = log2_32(u.bits.manh);
 336 #endif
 337                         shift = LDBL_MANH_SIZE - LDBL_NBIT_ADJ - pos;
 338                 } else {
 339 #if LDBL_MANL_SIZE > 32
 340                         pos = log2_64(u.bits.manl);
 341 #else
 342                         pos = log2_32(u.bits.manl);
 343 #endif
 344                         shift = LDBL_MANH_SIZE + LDBL_MANL_SIZE -
 345                             LDBL_NBIT_ADJ - pos;
 346                 }
 347                 sigfigs = (3 + LDBL_MANT_DIG - LDBL_NBIT_ADJ - shift) / 4;
 348                 *decpt = LDBL_MIN_EXP + LDBL_NBIT_ADJ -
 349                     ((shift + 3) & ~(4 - 1));
 350                 impnbit = 0;
 351                 break;
 352         case FP_INFINITE:
 353                 *decpt = INT_MAX;
 354                 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
 355         case FP_NAN:
 356                 *decpt = INT_MAX;
 357                 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
 358         default:
 359                 abort();
 360         }
 361
 362         /* FP_NORMAL or FP_SUBNORMAL */
 363
 364         if (ndigits == 0)               /* dtoa() compatibility */
 365                 ndigits = 1;
 366
 367         /*
 368          * For simplicity, we generate all the digits even if the
 369          * caller has requested fewer.
 370          */
 371         bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
 372         s0 = rv_alloc(bufsize);
 373
 374         /*
 375          * We work from right to left, first adding any requested zero
 376          * padding, then the least significant portion of the
 377          * mantissa, followed by the most significant.  The buffer is
 378          * filled with the byte values 0x0 through 0xf, which are
 379          * converted to xdigs[0x0] through xdigs[0xf] after the
 380          * rounding phase.
 381          */
 382         for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
 383                 *s = 0;
 384         for (; s > s0 + sigfigs - (LDBL_MANL_SIZE / 4) - 1 && s > s0; s--) {
 385                 *s = u.bits.manl & 0xf;
 386                 u.bits.manl >>= 4;
 387         }
 388         for (; s > s0; s--) {
 389                 *s = u.bits.manh & 0xf;
 390                 u.bits.manh >>= 4;
 391         }
 392
 393         /*
 394          * At this point, we have snarfed all the bits in the
 395          * mantissa, with the possible exception of the highest-order
 396          * (partial) nibble, which is dealt with by the next
 397          * statement.  That nibble is usually in manh, but it could be
 398          * in manl instead for small subnormals.  We also tack on the
 399          * implicit normalization bit if appropriate.
 400          */
 401         *s = u.bits.manh | u.bits.manl | impnbit;
 402
 403         /* If ndigits < 0, we are expected to auto-size the precision. */
 404         if (ndigits < 0) {
 405                 for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
 406                         ;
 407         }
 408
 409         if (sigfigs > ndigits && s0[ndigits] != 0)
 410                 dorounding(s0, ndigits, u.bits.sign, decpt);
 411
 412         s = s0 + ndigits;
 413         if (rve != NULL)
 414                 *rve = s;
 415         *s-- = '\0';
 416         for (; s >= s0; s--)
 417                 *s = xdigs[(unsigned int)*s];
 418
 419         return (s0);
 420 }
 421
 422 #else   /* (LDBL_MANT_DIG == DBL_MANT_DIG) */
 423
 424 char *
 425 __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
 426     char **rve)
 427 {
 428
 429         return (__hdtoa((double)e, xdigs, ndigits, decpt, sign, rve));
 430 }
 431
 432 #endif  /* (LDBL_MANT_DIG == DBL_MANT_DIG) */