gonzui

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