_hdtoa.c

Go to the documentation of this file.
00001 
00007 /*
00008  *    Copyright 2004-2006 Intel Corporation
00009  * 
00010  *    Licensed under the Apache License, Version 2.0 (the "License");
00011  *    you may not use this file except in compliance with the License.
00012  *    You may obtain a copy of the License at
00013  * 
00014  *        http://www.apache.org/licenses/LICENSE-2.0
00015  * 
00016  *    Unless required by applicable law or agreed to in writing, software
00017  *    distributed under the License is distributed on an "AS IS" BASIS,
00018  *    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00019  *    See the License for the specific language governing permissions and
00020  *    limitations under the License.
00021  */
00022 /*-
00023  * Copyright (c) 2004, 2005 David Schultz <das@FreeBSD.ORG>
00024  * All rights reserved.
00025  *
00026  * Redistribution and use in source and binary forms, with or without
00027  * modification, are permitted provided that the following conditions
00028  * are met:
00029  * 1. Redistributions of source code must retain the above copyright
00030  *    notice, this list of conditions and the following disclaimer.
00031  * 2. Redistributions in binary form must reproduce the above copyright
00032  *    notice, this list of conditions and the following disclaimer in the
00033  *    documentation and/or other materials provided with the distribution.
00034  *
00035  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
00036  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00037  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00038  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
00039  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
00040  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
00041  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00042  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00043  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
00044  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
00045  * SUCH DAMAGE.
00046  */
00047 
00048 // demmer: we need _GNU_SOURCE for linux systems to pull in fpclassify
00049 #define _GNU_SOURCE
00050 #include "compat/fpclassify.h"
00051 
00052 #include <sys/cdefs.h>
00053 #define __FBSDID(x) 
00054 __FBSDID("$FreeBSD: src/lib/libc/gdtoa/_hdtoa.c,v 1.3 2005/01/18 18:44:07 das Exp $");
00055 
00056 #include <float.h>
00057 #include <limits.h>
00058 #include <math.h>
00059 #include "fpmath.h"
00060 #include "gdtoaimp.h"
00061 
00062 /* Strings values used by dtoa() */
00063 #define INFSTR  "Infinity"
00064 #define NANSTR  "NaN"
00065 
00066 #define DBL_ADJ         (DBL_MAX_EXP - 2 + ((DBL_MANT_DIG - 1) % 4))
00067 #define LDBL_ADJ        (LDBL_MAX_EXP - 2 + ((LDBL_MANT_DIG - 1) % 4))
00068 
00069 /*
00070  * Round up the given digit string.  If the digit string is fff...f,
00071  * this procedure sets it to 100...0 and returns 1 to indicate that
00072  * the exponent needs to be bumped.  Otherwise, 0 is returned.
00073  */
00074 static int
00075 roundup(char *s0, int ndigits)
00076 {
00077         char *s;
00078 
00079         for (s = s0 + ndigits - 1; *s == 0xf; s--) {
00080                 if (s == s0) {
00081                         *s = 1;
00082                         return (1);
00083                 }
00084                 ++*s;
00085         }
00086         ++*s;
00087         return (0);
00088 }
00089 
00090 /*
00091  * Round the given digit string to ndigits digits according to the
00092  * current rounding mode.  Note that this could produce a string whose
00093  * value is not representable in the corresponding floating-point
00094  * type.  The exponent pointed to by decpt is adjusted if necessary.
00095  */
00096 static void
00097 dorounding(char *s0, int ndigits, int sign, int *decpt)
00098 {
00099         int adjust = 0; /* do we need to adjust the exponent? */
00100 
00101         switch (FLT_ROUNDS) {
00102         case 0:         /* toward zero */
00103         default:        /* implementation-defined */
00104                 break;
00105         case 1:         /* to nearest, halfway rounds to even */
00106                 if ((s0[ndigits] > 8) ||
00107                     (s0[ndigits] == 8 && s0[ndigits - 1] & 1))
00108                         adjust = roundup(s0, ndigits);
00109                 break;
00110         case 2:         /* toward +inf */
00111                 if (sign == 0)
00112                         adjust = roundup(s0, ndigits);
00113                 break;
00114         case 3:         /* toward -inf */
00115                 if (sign != 0)
00116                         adjust = roundup(s0, ndigits);
00117                 break;
00118         }
00119 
00120         if (adjust)
00121                 *decpt += 4;
00122 }
00123 
00124 /*
00125  * This procedure converts a double-precision number in IEEE format
00126  * into a string of hexadecimal digits and an exponent of 2.  Its
00127  * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
00128  * following exceptions:
00129  *
00130  * - An ndigits < 0 causes it to use as many digits as necessary to
00131  *   represent the number exactly.
00132  * - The additional xdigs argument should point to either the string
00133  *   "0123456789ABCDEF" or the string "0123456789abcdef", depending on
00134  *   which case is desired.
00135  * - This routine does not repeat dtoa's mistake of setting decpt
00136  *   to 9999 in the case of an infinity or NaN.  INT_MAX is used
00137  *   for this purpose instead.
00138  *
00139  * Note that the C99 standard does not specify what the leading digit
00140  * should be for non-zero numbers.  For instance, 0x1.3p3 is the same
00141  * as 0x2.6p2 is the same as 0x4.cp3.  This implementation chooses the
00142  * first digit so that subsequent digits are aligned on nibble
00143  * boundaries (before rounding).
00144  *
00145  * Inputs:      d, xdigs, ndigits
00146  * Outputs:     decpt, sign, rve
00147  */
00148 char *
00149 __hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
00150     char **rve)
00151 {
00152         static const int sigfigs = (DBL_MANT_DIG + 3) / 4;
00153         union IEEEd2bits u;
00154         char *s, *s0;
00155         int bufsize;
00156 
00157         u.d = d;
00158         *sign = u.bits.sign;
00159 
00160         switch (fpclassify(d)) {
00161         case FP_NORMAL:
00162                 *decpt = u.bits.exp - DBL_ADJ;
00163                 break;
00164         case FP_ZERO:
00165                 *decpt = 1;
00166                 return (nrv_alloc("0", rve, 1));
00167         case FP_SUBNORMAL:
00168                 u.d *= 0x1p514;
00169                 *decpt = u.bits.exp - (514 + DBL_ADJ);
00170                 break;
00171         case FP_INFINITE:
00172                 *decpt = INT_MAX;
00173                 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
00174         case FP_NAN:
00175                 *decpt = INT_MAX;
00176                 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
00177         default:
00178                 abort();
00179         }
00180 
00181         /* FP_NORMAL or FP_SUBNORMAL */
00182 
00183         if (ndigits == 0)               /* dtoa() compatibility */
00184                 ndigits = 1;
00185 
00186         /*
00187          * For simplicity, we generate all the digits even if the
00188          * caller has requested fewer.
00189          */
00190         bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
00191         s0 = rv_alloc(bufsize);
00192 
00193         /*
00194          * We work from right to left, first adding any requested zero
00195          * padding, then the least significant portion of the
00196          * mantissa, followed by the most significant.  The buffer is
00197          * filled with the byte values 0x0 through 0xf, which are
00198          * converted to xdigs[0x0] through xdigs[0xf] after the
00199          * rounding phase.
00200          */
00201         for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
00202                 *s = 0;
00203         for (; s > s0 + sigfigs - (DBL_MANL_SIZE / 4) - 1 && s > s0; s--) {
00204                 *s = u.bits.manl & 0xf;
00205                 u.bits.manl >>= 4;
00206         }
00207         for (; s > s0; s--) {
00208                 *s = u.bits.manh & 0xf;
00209                 u.bits.manh >>= 4;
00210         }
00211 
00212         /*
00213          * At this point, we have snarfed all the bits in the
00214          * mantissa, with the possible exception of the highest-order
00215          * (partial) nibble, which is dealt with by the next
00216          * statement.  We also tack on the implicit normalization bit.
00217          */
00218         *s = u.bits.manh | (1U << ((DBL_MANT_DIG - 1) % 4));
00219 
00220         /* If ndigits < 0, we are expected to auto-size the precision. */
00221         if (ndigits < 0) {
00222                 for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
00223                         ;
00224         }
00225 
00226         if (sigfigs > ndigits && s0[ndigits] != 0)
00227                 dorounding(s0, ndigits, u.bits.sign, decpt);
00228 
00229         s = s0 + ndigits;
00230         if (rve != NULL)
00231                 *rve = s;
00232         *s-- = '\0';
00233         for (; s >= s0; s--)
00234                 *s = xdigs[(unsigned int)*s];
00235 
00236         return (s0);
00237 }
00238 
00239 #if (LDBL_MANT_DIG > DBL_MANT_DIG)
00240 
00241 /*
00242  * This is the long double version of __hdtoa().
00243  */
00244 char *
00245 __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
00246     char **rve)
00247 {
00248         static const int sigfigs = (LDBL_MANT_DIG + 3) / 4;
00249         union IEEEl2bits u;
00250         char *s, *s0;
00251         int bufsize;
00252 
00253         u.e = e;
00254         *sign = u.bits.sign;
00255 
00256         switch (fpclassify(e)) {
00257         case FP_NORMAL:
00258                 *decpt = u.bits.exp - LDBL_ADJ;
00259                 break;
00260         case FP_ZERO:
00261                 *decpt = 1;
00262                 return (nrv_alloc("0", rve, 1));
00263         case FP_SUBNORMAL:
00264                 u.e *= 0x1p514L;
00265                 *decpt = u.bits.exp - (514 + LDBL_ADJ);
00266                 break;
00267         case FP_INFINITE:
00268                 *decpt = INT_MAX;
00269                 return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
00270         case FP_NAN:
00271                 *decpt = INT_MAX;
00272                 return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
00273         default:
00274                 abort();
00275         }
00276 
00277         /* FP_NORMAL or FP_SUBNORMAL */
00278 
00279         if (ndigits == 0)               /* dtoa() compatibility */
00280                 ndigits = 1;
00281 
00282         /*
00283          * For simplicity, we generate all the digits even if the
00284          * caller has requested fewer.
00285          */
00286         bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
00287         s0 = rv_alloc(bufsize);
00288 
00289         /*
00290          * We work from right to left, first adding any requested zero
00291          * padding, then the least significant portion of the
00292          * mantissa, followed by the most significant.  The buffer is
00293          * filled with the byte values 0x0 through 0xf, which are
00294          * converted to xdigs[0x0] through xdigs[0xf] after the
00295          * rounding phase.
00296          */
00297         for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
00298                 *s = 0;
00299         for (; s > s0 + sigfigs - (LDBL_MANL_SIZE / 4) - 1 && s > s0; s--) {
00300                 *s = u.bits.manl & 0xf;
00301                 u.bits.manl >>= 4;
00302         }
00303         for (; s > s0; s--) {
00304                 *s = u.bits.manh & 0xf;
00305                 u.bits.manh >>= 4;
00306         }
00307 
00308         /*
00309          * At this point, we have snarfed all the bits in the
00310          * mantissa, with the possible exception of the highest-order
00311          * (partial) nibble, which is dealt with by the next
00312          * statement.  We also tack on the implicit normalization bit.
00313          */
00314         *s = u.bits.manh | (1U << ((LDBL_MANT_DIG - 1) % 4));
00315 
00316         /* If ndigits < 0, we are expected to auto-size the precision. */
00317         if (ndigits < 0) {
00318                 for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
00319                         ;
00320         }
00321 
00322         if (sigfigs > ndigits && s0[ndigits] != 0)
00323                 dorounding(s0, ndigits, u.bits.sign, decpt);
00324 
00325         s = s0 + ndigits;
00326         if (rve != NULL)
00327                 *rve = s;
00328         *s-- = '\0';
00329         for (; s >= s0; s--)
00330                 *s = xdigs[(unsigned int)*s];
00331 
00332         return (s0);
00333 }
00334 
00335 #else   /* (LDBL_MANT_DIG == DBL_MANT_DIG) */
00336 
00337 char *
00338 __hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
00339     char **rve)
00340 {
00341 
00342         return (__hdtoa((double)e, xdigs, ndigits, decpt, sign, rve));
00343 }
00344 
00345 #endif  /* (LDBL_MANT_DIG == DBL_MANT_DIG) */

Generated on Thu Jun 7 16:56:48 2007 for DTN Reference Implementation by  doxygen 1.5.1