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