lib/libc/gdtoa/hdtoa.c

   1 /*      $NetBSD: hdtoa.c,v 1.9 2011/07/04 11:46:41 mrg Exp $    */
   2
   3 /*-
   4  * Copyright (c) 2004, 2005 David Schultz <das@FreeBSD.ORG>
   5  * All rights reserved.
   6  *
   7  * Redistribution and use in source and binary forms, with or without
   8  * modification, are permitted provided that the following conditions
   9  * are met:
  10  * 1. Redistributions of source code must retain the above copyright
  11  *    notice, this list of conditions and the following disclaimer.
  12  * 2. Redistributions in binary form must reproduce the above copyright
  13  *    notice, this list of conditions and the following disclaimer in the
  14  *    documentation and/or other materials provided with the distribution.
  15  *
  16  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  18  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  19  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  20  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  21  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  22  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  23  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  24  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  25  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  26  * SUCH DAMAGE.
  27  */
  28
  29 #include <sys/cdefs.h>
  30 #if 0
  31 __FBSDID("$FreeBSD: src/lib/libc/gdtoa/_hdtoa.c,v 1.4 2007/01/03 04:57:58 das Exp $");
  32 #else
  33 __RCSID("$NetBSD: hdtoa.c,v 1.9 2011/07/04 11:46:41 mrg Exp $");
  34 #endif
  35
  36 #include <float.h>
  37 #include <limits.h>
  38 #include <math.h>
  39 #ifndef __vax__
  40 #include <machine/ieee.h>
  41 #else
  42 #include <machine/vaxfp.h>
  43 #define ieee_double_u vax_dfloating_u
  44 #define dblu_d dfltu_d
  45 #define dblu_dbl dfltu_dflt
  46 #define dbl_sign dflt_sign
  47 #define dbl_exp dflt_exp
  48 #define dbl_frach dflt_frach
  49 #define dbl_fracm dflt_fracm
  50 #define dbl_fracl dflt_fracl
  51 #define DBL_FRACHBITS   DFLT_FRACHBITS
  52 #define DBL_FRACMBITS   DFLT_FRACMBITS
  53 #define DBL_FRACLBITS   DFLT_FRACLBITS
  54 #define DBL_EXPBITS     DFLT_EXPBITS
  55 #endif
  56 #include "gdtoaimp.h"
  57
  58 /* Strings values used by dtoa() */
  59 #define INFSTR  "Infinity"
  60 #define NANSTR  "NaN"
  61
  62 #define DBL_ADJ         (DBL_MAX_EXP - 2 + ((DBL_MANT_DIG - 1) % 4))
  63 #define LDBL_ADJ        (LDBL_MAX_EXP - 2 + ((LDBL_MANT_DIG - 1) % 4))
  64
  65 /*
  66  * Round up the given digit string.  If the digit string is fff...f,
  67  * this procedure sets it to 100...0 and returns 1 to indicate that
  68  * the exponent needs to be bumped.  Otherwise, 0 is returned.
  69  */
  70 static int
  71 roundup(char *s0, int ndigits)
  72 {
  73         char *s;
  74
  75         for (s = s0 + ndigits - 1; *s == 0xf; s--) {
  76                 if (s == s0) {
  77                         *s = 1;
  78                         return (1);
  79                 }
  80                 *s = 0;
  81         }
  82         ++*s;
  83         return (0);
  84 }
  85
  86 /*
  87  * Round the given digit string to ndigits digits according to the
  88  * current rounding mode.  Note that this could produce a string whose
  89  * value is not representable in the corresponding floating-point
  90  * type.  The exponent pointed to by decpt is adjusted if necessary.
  91  */
  92 static void
  93 dorounding(char *s0, int ndigits, int sign, int *decpt)
  94 {
  95         int adjust = 0; /* do we need to adjust the exponent? */
  96
  97         switch (FLT_ROUNDS) {
  98         case 0:         /* toward zero */
  99         default:        /* implementation-defined */
 100                 break;
 101         case 1:         /* to nearest, halfway rounds to even */
 102                 if ((s0[ndigits] > 8) ||
 103                     (s0[ndigits] == 8 && s0[ndigits - 1] & 1))
 104                         adjust = roundup(s0, ndigits);
 105                 break;
 106         case 2:         /* toward +inf */
 107                 if (sign == 0)
 108                         adjust = roundup(s0, ndigits);
 109                 break;
 110         case 3:         /* toward -inf */
 111                 if (sign != 0)
 112                         adjust = roundup(s0, ndigits);
 113                 break;
 114         }
 115
 116         if (adjust)
 117                 *decpt += 4;
 118 }
 119
 120 /*
 121  * This procedure converts a double-precision number in IEEE format
 122  * into a string of hexadecimal digits and an exponent of 2.  Its
 123  * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
 124  * following exceptions:
 125  *
 126  * - An ndigits < 0 causes it to use as many digits as necessary to
 127  *   represent the number exactly.
 128  * - The additional xdigs argument should point to either the string
 129  *   "0123456789ABCDEF" or the string "0123456789abcdef", depending on
 130  *   which case is desired.
 131  * - This routine does not repeat dtoa's mistake of setting decpt
 132  *   to 9999 in the case of an infinity or NaN.  INT_MAX is used
 133  *   for this purpose instead.
 134  *
 135  * Note that the C99 standard does not specify what the leading digit
 136  * should be for non-zero numbers.  For instance, 0x1.3p3 is the same
 137  * as 0x2.6p2 is the same as 0x4.cp3.  This implementation chooses the
 138  * first digit so that subsequent digits are aligned on nibble
 139  * boundaries (before rounding).
 140  *
 141  * Inputs:      d, xdigs, ndigits
 142  * Outputs:     decpt, sign, rve
 143  */
 144 char *
 145 hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
 146     char **rve)
 147 {
 148         static const int sigfigs = (DBL_MANT_DIG + 3) / 4;
 149         union ieee_double_u u;
 150         char *s, *s0;
 151         size_t bufsize;
 152
 153         u.dblu_d = d;
 154         *sign = u.dblu_dbl.dbl_sign;
 155
 156         switch (fpclassify(d)) {
 157         case FP_NORMAL:
 158                 *decpt = u.dblu_dbl.dbl_exp - DBL_ADJ;
 159                 break;
 160         case FP_ZERO:
 161                 *decpt = 1;
 162                 return (nrv_alloc("0", rve, 1));
 163         case FP_SUBNORMAL:
 164 #ifdef __vax__
 165                 /* (DBL_MAX_EXP=127 / 2) + 2 = 65? */
 166                 u.dblu_d *= 0x1p65;
 167                 *decpt = u.dblu_dbl.dbl_exp - (65 + DBL_ADJ);
 168 #else
 169                 /* (DBL_MAX_EXP=1024 / 2) + 2 = 514? */
 170                 u.dblu_d *= 0x1p514;
 171                 *decpt = u.dblu_dbl.dbl_exp - (514 + DBL_ADJ);
 172 #endif
 173                 break;
 174         case FP_INFINITE:
 175                 *decpt = INT_MAX;
 176                 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
 177         case FP_NAN:
 178                 *decpt = INT_MAX;
 179                 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
 180         default:
 181                 abort();
 182         }
 183
 184         /* FP_NORMAL or FP_SUBNORMAL */
 185
 186         if (ndigits == 0)               /* dtoa() compatibility */
 187                 ndigits = 1;
 188
 189         /*
 190          * For simplicity, we generate all the digits even if the
 191          * caller has requested fewer.
 192          */
 193         bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
 194         s0 = rv_alloc(bufsize);
 195         if (s0 == NULL)
 196                 return NULL;
 197
 198         /*
 199          * We work from right to left, first adding any requested zero
 200          * padding, then the least significant portion of the
 201          * mantissa, followed by the most significant.  The buffer is
 202          * filled with the byte values 0x0 through 0xf, which are
 203          * converted to xdigs[0x0] through xdigs[0xf] after the
 204          * rounding phase.
 205          */
 206         for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
 207                 *s = 0;
 208         for (; s > s0 + sigfigs - (DBL_FRACLBITS / 4) - 1 && s > s0; s--) {
 209                 *s = u.dblu_dbl.dbl_fracl & 0xf;
 210                 u.dblu_dbl.dbl_fracl >>= 4;
 211         }
 212 #ifdef DBL_FRACMBITS
 213         for (; s > s0; s--) {
 214                 *s = u.dblu_dbl.dbl_fracm & 0xf;
 215                 u.dblu_dbl.dbl_fracm >>= 4;
 216         }
 217 #endif
 218         for (; s > s0; s--) {
 219                 *s = u.dblu_dbl.dbl_frach & 0xf;
 220                 u.dblu_dbl.dbl_frach >>= 4;
 221         }
 222
 223         /*
 224          * At this point, we have snarfed all the bits in the
 225          * mantissa, with the possible exception of the highest-order
 226          * (partial) nibble, which is dealt with by the next
 227          * statement.  We also tack on the implicit normalization bit.
 228          */
 229         *s = u.dblu_dbl.dbl_frach | (1U << ((DBL_MANT_DIG - 1) % 4));
 230
 231         /* If ndigits < 0, we are expected to auto-size the precision. */
 232         if (ndigits < 0) {
 233                 for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
 234                         continue;
 235         }
 236
 237         if (sigfigs > ndigits && s0[ndigits] != 0)
 238                 dorounding(s0, ndigits, u.dblu_dbl.dbl_sign, decpt);
 239
 240         s = s0 + ndigits;
 241         if (rve != NULL)
 242                 *rve = s;
 243         *s-- = '\0';
 244         for (; s >= s0; s--)
 245                 *s = xdigs[(unsigned int)*s];
 246
 247         return (s0);
 248 }
 249
 250 #if (LDBL_MANT_DIG > DBL_MANT_DIG)
 251
 252 /*
 253  * This is the long double version of hdtoa().
 254  */
 255 char *
 256 hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
 257     char **rve)
 258 {
 259         static const int sigfigs = (LDBL_MANT_DIG + 3) / 4;
 260         union ieee_ext_u u;
 261         char *s, *s0;
 262         size_t bufsize;
 263
 264         memset(&u, 0, sizeof u);
 265         u.extu_ld = e;
 266         *sign = u.extu_ext.ext_sign;
 267
 268         switch (fpclassify(e)) {
 269         case FP_NORMAL:
 270                 *decpt = u.extu_ext.ext_exp - LDBL_ADJ;
 271                 break;
 272         case FP_ZERO:
 273                 *decpt = 1;
 274                 return (nrv_alloc("0", rve, 1));
 275         case FP_SUBNORMAL:
 276                 u.extu_ld *= 0x1p514L;
 277                 *decpt = u.extu_ext.ext_exp - (514 + LDBL_ADJ);
 278                 break;
 279         case FP_INFINITE:
 280                 *decpt = INT_MAX;
 281                 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
 282         case FP_NAN:
 283                 *decpt = INT_MAX;
 284                 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
 285         default:
 286                 abort();
 287         }
 288
 289         /* FP_NORMAL or FP_SUBNORMAL */
 290
 291         if (ndigits == 0)               /* dtoa() compatibility */
 292                 ndigits = 1;
 293
 294         /*
 295          * For simplicity, we generate all the digits even if the
 296          * caller has requested fewer.
 297          */
 298         bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
 299         s0 = rv_alloc(bufsize);
 300         if (s0 == NULL)
 301                 return NULL;
 302
 303         /*
 304          * We work from right to left, first adding any requested zero
 305          * padding, then the least significant portion of the
 306          * mantissa, followed by the most significant.  The buffer is
 307          * filled with the byte values 0x0 through 0xf, which are
 308          * converted to xdigs[0x0] through xdigs[0xf] after the
 309          * rounding phase.
 310          */
 311         for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
 312                 *s = 0;
 313         for (; s > s0 + sigfigs - (EXT_FRACLBITS / 4) - 1 && s > s0; s--) {
 314                 *s = u.extu_ext.ext_fracl & 0xf;
 315                 u.extu_ext.ext_fracl >>= 4;
 316         }
 317 #ifdef EXT_FRACHMBITS
 318         for (; s > s0; s--) {
 319                 *s = u.extu_ext.ext_frachm & 0xf;
 320                 u.extu_ext.ext_frachm >>= 4;
 321         }
 322 #endif
 323 #ifdef EXT_FRACLMBITS
 324         for (; s > s0; s--) {
 325                 *s = u.extu_ext.ext_fraclm & 0xf;
 326                 u.extu_ext.ext_fraclm >>= 4;
 327         }
 328 #endif
 329         for (; s > s0; s--) {
 330                 *s = u.extu_ext.ext_frach & 0xf;
 331                 u.extu_ext.ext_frach >>= 4;
 332         }
 333
 334         /*
 335          * At this point, we have snarfed all the bits in the
 336          * mantissa, with the possible exception of the highest-order
 337          * (partial) nibble, which is dealt with by the next
 338          * statement.  We also tack on the implicit normalization bit.
 339          */
 340         *s = u.extu_ext.ext_frach | (1U << ((LDBL_MANT_DIG - 1) % 4));
 341
 342         /* If ndigits < 0, we are expected to auto-size the precision. */
 343         if (ndigits < 0) {
 344                 for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
 345                         continue;
 346         }
 347
 348         if (sigfigs > ndigits && s0[ndigits] != 0)
 349                 dorounding(s0, ndigits, u.extu_ext.ext_sign, decpt);
 350
 351         s = s0 + ndigits;
 352         if (rve != NULL)
 353                 *rve = s;
 354         *s-- = '\0';
 355         for (; s >= s0; s--)
 356                 *s = xdigs[(unsigned int)*s];
 357
 358         return (s0);
 359 }
 360
 361 #else   /* (LDBL_MANT_DIG == DBL_MANT_DIG) */
 362
 363 char *
 364 hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
 365     char **rve)
 366 {
 367
 368         return (hdtoa((double)e, xdigs, ndigits, decpt, sign, rve));
 369 }
 370
 371 #endif  /* (LDBL_MANT_DIG == DBL_MANT_DIG) */