1 /* Copyright (c) 2005, Dmitry Xmelkov
2 All rights reserved.
3 Rewritten in C by Soren Kuula
4 Redistribution and use in source and binary forms, with or without
5 modification, are permitted provided that the following conditions are met:
6 * Redistributions of source code must retain the above copyright
7 notice, this list of conditions and the following disclaimer.
8 * Redistributions in binary form must reproduce the above copyright
9 notice, this list of conditions and the following disclaimer in
10 the documentation and/or other materials provided with the
11 distribution.
12 * Neither the name of the copyright holders nor the names of
13 contributors may be used to endorse or promote products derived
14 from this software without specific prior written permission.
15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
16 AND 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 COPYRIGHT OWNER OR CONTRIBUTORS BE
19 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
20 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
21 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
22 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
23 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
24 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
25 POSSIBILITY OF SUCH DAMAGE. */
26
27 #define _NEED_IO_FLOAT32
28 #include "dtoa.h"
29
30 /*
31 * 2^b ~= f * r * 10^e
32 * where
33 * i = b div 8
34 * r = 2^(b mod 8)
35 * f = factorTable[i]
36 * e = exponentTable[i]
37 */
38 static const int8_t exponentTable[32] = {
39 -36, -33, -31, -29, -26, -24, -21, -19,
40 -17, -14, -12, -9, -7, -4, -2, 0,
41 3, 5, 8, 10, 12, 15, 17, 20,
42 22, 24, 27, 29, 32, 34, 36, 39
43 };
44
45 static const uint32_t factorTable[32] = {
46 2295887404UL,
47 587747175UL,
48 1504632769UL,
49 3851859889UL,
50 986076132UL,
51 2524354897UL,
52 646234854UL,
53 1654361225UL,
54 4235164736UL,
55 1084202172UL,
56 2775557562UL,
57 710542736UL,
58 1818989404UL,
59 465661287UL,
60 1192092896UL,
61 3051757813UL,
62 781250000UL,
63 2000000000UL,
64 512000000UL,
65 1310720000UL,
66 3355443200UL,
67 858993459UL,
68 2199023256UL,
69 562949953UL,
70 1441151881UL,
71 3689348815UL,
72 944473297UL,
73 2417851639UL,
74 618970020UL,
75 1584563250UL,
76 4056481921UL,
77 1038459372UL
78 };
79
80 #define max(a, b) ({\
81 __typeof(a) _a = a;\
82 __typeof(b) _b = b;\
83 _a > _b ? _a : _b; })
84
85 #define min(a, b) ({\
86 __typeof(a) _a = a;\
87 __typeof(b) _b = b;\
88 _a < _b ? _a : _b; })
89
__ftoa_engine(float val,struct dtoa * ftoa,int maxDigits,bool fmode,int maxDecimals)90 int __ftoa_engine(float val, struct dtoa *ftoa, int maxDigits, bool fmode, int maxDecimals)
91 {
92 uint8_t flags = 0;
93
94 union {
95 float v;
96 uint32_t u;
97 } x;
98
99 x.v = val;
100
101 uint32_t frac = x.u & 0x007fffffUL;
102
103 if (maxDigits > FTOA_MAX_DIG)
104 maxDigits = FTOA_MAX_DIG;
105
106 /* Read the sign, shift the exponent in place and delete it from frac.
107 */
108 if (x.u & ((uint32_t)1 << 31))
109 flags = DTOA_MINUS;
110
111 uint8_t exp = x.u >> 23;
112
113 ftoa->exp = 0;
114
115 /*
116 * Test for easy cases, zero and NaN
117 */
118 if(exp==0 && frac==0)
119 {
120 flags |= DTOA_ZERO;
121 ftoa->digits[0] = '0';
122 maxDigits = 1;
123 } else if(exp == 0xff) {
124 if(frac == 0)
125 flags |= DTOA_INF;
126 else
127 flags |= DTOA_NAN;
128 } else {
129
130 /* The implicit leading 1 is made explicit, except if value
131 * subnormal, in which case exp needs to be adjusted by one
132 */
133 if (exp == 0)
134 exp = 1;
135 else
136 frac |= (1UL<<23);
137
138 uint8_t idx = exp>>3;
139 int8_t exp10 = exponentTable[idx];
140
141 /*
142 * We COULD try making the multiplication in situ, where we make
143 * frac and a 64 bit int overlap in memory and select/weigh the
144 * upper 32 bits that way. For starters, this is less risky:
145 */
146 int64_t prod = (int64_t)frac * (int64_t)factorTable[idx];
147
148 /*
149 * The expConvFactorTable are factor are correct iff the lower 3 exponent
150 * bits are 1 (=7). Else we need to compensate by divding frac.
151 * If the lower 3 bits are 7 we are right.
152 * If the lower 3 bits are 6 we right-shift once
153 * ..
154 * If the lower 3 bits are 0 we right-shift 7x
155 */
156 prod >>= (15-(exp & 7));
157
158 /*
159 * Now convert to decimal.
160 */
161
162 uint8_t hadNonzeroDigit = 0; /* have seen a non-zero digit flag */
163 uint8_t outputIdx = 0;
164 int64_t decimal = 100000000000000ull;
165 uint8_t saveMaxDigits = maxDigits;
166
167 do {
168 /* Compute next digit */
169 char digit = prod / decimal + '0';
170
171 if(!hadNonzeroDigit)
172 {
173 /* Don't return results with a leading zero! Instead
174 * skip those and decrement exp10 accordingly.
175 */
176 if (digit == '0') {
177 exp10--;
178 prod = prod % decimal;
179 decimal /= 10;
180 continue;
181 }
182
183 /* Found the first non-zero digit */
184 hadNonzeroDigit = 1;
185
186 /* If limiting decimals... */
187 if(fmode)
188 {
189 maxDigits = min(maxDigits, max(1, maxDecimals + exp10 + 1));
190 if (maxDigits == 0)
191 break;
192 }
193 }
194 prod = prod % decimal;
195 decimal /= 10;
196
197 /* Now we have a digit. */
198 if(digit < '0' + 10) {
199 /* normal case. */
200 ftoa->digits[outputIdx] = digit;
201 } else {
202 /*
203 * Uh, oh. Something went wrong with our conversion. Write
204 * 9s and use the round-up code to report back to the caller
205 * that we've carried
206 */
207 for(outputIdx = 0; outputIdx < maxDigits; outputIdx++)
208 ftoa->digits[outputIdx] = '9';
209 goto round_up;
210 }
211 outputIdx++;
212 } while (outputIdx<maxDigits);
213
214 /* Rounding: */
215 decimal *= 10;
216 if (prod - (decimal >> 1) >= 0)
217 {
218 uint8_t rounded;
219 round_up:
220
221 rounded = 0;
222 while(outputIdx != 0) {
223
224 /* Increment digit, check if we're done */
225 if(++ftoa->digits[outputIdx-1] < '0' + 10) {
226 rounded = 1;
227 break;
228 }
229
230 ftoa->digits[--outputIdx] = '0';
231 /* and the loop continues, carrying to next digit. */
232 }
233
234 if (!rounded) {
235
236 /* Rounded past the first digit;
237 * reset the leading digit to 1,
238 * bump exp and tell the caller we've carried.
239 * The remaining digits will already be '0',
240 * so we don't need to mess with them
241 */
242 ftoa->digits[0] = '1';
243 exp10++;
244 if (fmode)
245 maxDigits = min(saveMaxDigits, max(1, maxDecimals + exp10 + 1));
246 flags |= DTOA_CARRY;
247 }
248 }
249 ftoa->exp = exp10;
250 }
251 ftoa->flags = flags;
252 return maxDigits;
253 }
254
