1*86d7f5d3SJohn Marino /*-
2*86d7f5d3SJohn Marino * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG>
3*86d7f5d3SJohn Marino * All rights reserved.
4*86d7f5d3SJohn Marino *
5*86d7f5d3SJohn Marino * Redistribution and use in source and binary forms, with or without
6*86d7f5d3SJohn Marino * modification, are permitted provided that the following conditions
7*86d7f5d3SJohn Marino * are met:
8*86d7f5d3SJohn Marino * 1. Redistributions of source code must retain the above copyright
9*86d7f5d3SJohn Marino * notice, this list of conditions and the following disclaimer.
10*86d7f5d3SJohn Marino * 2. Redistributions in binary form must reproduce the above copyright
11*86d7f5d3SJohn Marino * notice, this list of conditions and the following disclaimer in the
12*86d7f5d3SJohn Marino * documentation and/or other materials provided with the distribution.
13*86d7f5d3SJohn Marino *
14*86d7f5d3SJohn Marino * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15*86d7f5d3SJohn Marino * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16*86d7f5d3SJohn Marino * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17*86d7f5d3SJohn Marino * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18*86d7f5d3SJohn Marino * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19*86d7f5d3SJohn Marino * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20*86d7f5d3SJohn Marino * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21*86d7f5d3SJohn Marino * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22*86d7f5d3SJohn Marino * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23*86d7f5d3SJohn Marino * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24*86d7f5d3SJohn Marino * SUCH DAMAGE.
25*86d7f5d3SJohn Marino *
26*86d7f5d3SJohn Marino * $FreeBSD: src/lib/libc/gdtoa/_hdtoa.c,v 1.7 2008/04/12 14:53:52 das Exp $
27*86d7f5d3SJohn Marino */
28*86d7f5d3SJohn Marino
29*86d7f5d3SJohn Marino #include <float.h>
30*86d7f5d3SJohn Marino #include <limits.h>
31*86d7f5d3SJohn Marino #include <math.h>
32*86d7f5d3SJohn Marino
33*86d7f5d3SJohn Marino #include "../stdio/floatio.h"
34*86d7f5d3SJohn Marino #include "fpmath.h"
35*86d7f5d3SJohn Marino #include "gdtoaimp.h"
36*86d7f5d3SJohn Marino
37*86d7f5d3SJohn Marino /* Strings values used by dtoa() */
38*86d7f5d3SJohn Marino #define INFSTR "Infinity"
39*86d7f5d3SJohn Marino #define NANSTR "NaN"
40*86d7f5d3SJohn Marino
41*86d7f5d3SJohn Marino #define DBL_ADJ (DBL_MAX_EXP - 2)
42*86d7f5d3SJohn Marino #define SIGFIGS ((DBL_MANT_DIG + 3) / 4 + 1)
43*86d7f5d3SJohn Marino
44*86d7f5d3SJohn Marino static const float one[] = { 1.0f, -1.0f };
45*86d7f5d3SJohn Marino
46*86d7f5d3SJohn Marino /*
47*86d7f5d3SJohn Marino * This procedure converts a double-precision number in IEEE format
48*86d7f5d3SJohn Marino * into a string of hexadecimal digits and an exponent of 2. Its
49*86d7f5d3SJohn Marino * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
50*86d7f5d3SJohn Marino * following exceptions:
51*86d7f5d3SJohn Marino *
52*86d7f5d3SJohn Marino * - An ndigits < 0 causes it to use as many digits as necessary to
53*86d7f5d3SJohn Marino * represent the number exactly.
54*86d7f5d3SJohn Marino * - The additional xdigs argument should point to either the string
55*86d7f5d3SJohn Marino * "0123456789ABCDEF" or the string "0123456789abcdef", depending on
56*86d7f5d3SJohn Marino * which case is desired.
57*86d7f5d3SJohn Marino * - This routine does not repeat dtoa's mistake of setting decpt
58*86d7f5d3SJohn Marino * to 9999 in the case of an infinity or NaN. INT_MAX is used
59*86d7f5d3SJohn Marino * for this purpose instead.
60*86d7f5d3SJohn Marino *
61*86d7f5d3SJohn Marino * Note that the C99 standard does not specify what the leading digit
62*86d7f5d3SJohn Marino * should be for non-zero numbers. For instance, 0x1.3p3 is the same
63*86d7f5d3SJohn Marino * as 0x2.6p2 is the same as 0x4.cp3. This implementation always makes
64*86d7f5d3SJohn Marino * the leading digit a 1. This ensures that the exponent printed is the
65*86d7f5d3SJohn Marino * actual base-2 exponent, i.e., ilogb(d).
66*86d7f5d3SJohn Marino *
67*86d7f5d3SJohn Marino * Inputs: d, xdigs, ndigits
68*86d7f5d3SJohn Marino * Outputs: decpt, sign, rve
69*86d7f5d3SJohn Marino */
70*86d7f5d3SJohn Marino char *
__hdtoa(double d,const char * xdigs,int ndigits,int * decpt,int * sign,char ** rve)71*86d7f5d3SJohn Marino __hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
72*86d7f5d3SJohn Marino char **rve)
73*86d7f5d3SJohn Marino {
74*86d7f5d3SJohn Marino union IEEEd2bits u;
75*86d7f5d3SJohn Marino char *s, *s0;
76*86d7f5d3SJohn Marino int bufsize;
77*86d7f5d3SJohn Marino uint32_t manh, manl;
78*86d7f5d3SJohn Marino
79*86d7f5d3SJohn Marino u.d = d;
80*86d7f5d3SJohn Marino *sign = u.bits.sign;
81*86d7f5d3SJohn Marino
82*86d7f5d3SJohn Marino switch (fpclassify(d)) {
83*86d7f5d3SJohn Marino case FP_NORMAL:
84*86d7f5d3SJohn Marino *decpt = u.bits.exp - DBL_ADJ;
85*86d7f5d3SJohn Marino break;
86*86d7f5d3SJohn Marino case FP_ZERO:
87*86d7f5d3SJohn Marino *decpt = 1;
88*86d7f5d3SJohn Marino return (nrv_alloc("0", rve, 1));
89*86d7f5d3SJohn Marino case FP_SUBNORMAL:
90*86d7f5d3SJohn Marino u.d *= 0x1p514;
91*86d7f5d3SJohn Marino *decpt = u.bits.exp - (514 + DBL_ADJ);
92*86d7f5d3SJohn Marino break;
93*86d7f5d3SJohn Marino case FP_INFINITE:
94*86d7f5d3SJohn Marino *decpt = INT_MAX;
95*86d7f5d3SJohn Marino return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
96*86d7f5d3SJohn Marino default: /* FP_NAN or unrecognized */
97*86d7f5d3SJohn Marino *decpt = INT_MAX;
98*86d7f5d3SJohn Marino return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
99*86d7f5d3SJohn Marino }
100*86d7f5d3SJohn Marino
101*86d7f5d3SJohn Marino /* FP_NORMAL or FP_SUBNORMAL */
102*86d7f5d3SJohn Marino
103*86d7f5d3SJohn Marino if (ndigits == 0) /* dtoa() compatibility */
104*86d7f5d3SJohn Marino ndigits = 1;
105*86d7f5d3SJohn Marino
106*86d7f5d3SJohn Marino /*
107*86d7f5d3SJohn Marino * If ndigits < 0, we are expected to auto-size, so we allocate
108*86d7f5d3SJohn Marino * enough space for all the digits.
109*86d7f5d3SJohn Marino */
110*86d7f5d3SJohn Marino bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
111*86d7f5d3SJohn Marino s0 = rv_alloc(bufsize);
112*86d7f5d3SJohn Marino
113*86d7f5d3SJohn Marino /* Round to the desired number of digits. */
114*86d7f5d3SJohn Marino if (SIGFIGS > ndigits && ndigits > 0) {
115*86d7f5d3SJohn Marino float redux = one[u.bits.sign];
116*86d7f5d3SJohn Marino int offset = 4 * ndigits + DBL_MAX_EXP - 4 - DBL_MANT_DIG;
117*86d7f5d3SJohn Marino u.bits.exp = offset;
118*86d7f5d3SJohn Marino u.d += redux;
119*86d7f5d3SJohn Marino u.d -= redux;
120*86d7f5d3SJohn Marino *decpt += u.bits.exp - offset;
121*86d7f5d3SJohn Marino }
122*86d7f5d3SJohn Marino
123*86d7f5d3SJohn Marino manh = u.bits.manh;
124*86d7f5d3SJohn Marino manl = u.bits.manl;
125*86d7f5d3SJohn Marino *s0 = '1';
126*86d7f5d3SJohn Marino for (s = s0 + 1; s < s0 + bufsize; s++) {
127*86d7f5d3SJohn Marino *s = xdigs[(manh >> (DBL_MANH_SIZE - 4)) & 0xf];
128*86d7f5d3SJohn Marino manh = (manh << 4) | (manl >> (DBL_MANL_SIZE - 4));
129*86d7f5d3SJohn Marino manl <<= 4;
130*86d7f5d3SJohn Marino }
131*86d7f5d3SJohn Marino
132*86d7f5d3SJohn Marino /* If ndigits < 0, we are expected to auto-size the precision. */
133*86d7f5d3SJohn Marino if (ndigits < 0) {
134*86d7f5d3SJohn Marino for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--)
135*86d7f5d3SJohn Marino ;
136*86d7f5d3SJohn Marino }
137*86d7f5d3SJohn Marino
138*86d7f5d3SJohn Marino s = s0 + ndigits;
139*86d7f5d3SJohn Marino *s = '\0';
140*86d7f5d3SJohn Marino if (rve != NULL)
141*86d7f5d3SJohn Marino *rve = s;
142*86d7f5d3SJohn Marino return (s0);
143*86d7f5d3SJohn Marino }
144