]>
Commit | Line | Data |
---|---|---|
e778b260 RN |
1 | #include <stdbool.h> |
2 | #include <string.h> | |
3 | ||
4 | #include "fpconv.h" | |
5 | #include "powers.h" | |
6 | ||
7 | #define fracmask 0x000FFFFFFFFFFFFFU | |
8 | #define expmask 0x7FF0000000000000U | |
9 | #define hiddenbit 0x0010000000000000U | |
10 | #define signmask 0x8000000000000000U | |
11 | #define expbias (1023 + 52) | |
12 | ||
13 | #define absv(n) ((n) < 0 ? -(n) : (n)) | |
14 | #define minv(a, b) ((a) < (b) ? (a) : (b)) | |
15 | ||
16 | static uint64_t tens[] = { | |
17 | 10000000000000000000U, 1000000000000000000U, 100000000000000000U, | |
18 | 10000000000000000U, 1000000000000000U, 100000000000000U, | |
19 | 10000000000000U, 1000000000000U, 100000000000U, | |
20 | 10000000000U, 1000000000U, 100000000U, | |
21 | 10000000U, 1000000U, 100000U, | |
22 | 10000U, 1000U, 100U, | |
23 | 10U, 1U | |
24 | }; | |
25 | ||
26 | static inline uint64_t get_dbits(double d) | |
27 | { | |
28 | union { | |
29 | double dbl; | |
30 | uint64_t i; | |
31 | } dbl_bits = { d }; | |
32 | ||
33 | return dbl_bits.i; | |
34 | } | |
35 | ||
36 | static Fp build_fp(double d) | |
37 | { | |
38 | uint64_t bits = get_dbits(d); | |
39 | ||
40 | Fp fp; | |
41 | fp.frac = bits & fracmask; | |
42 | fp.exp = (bits & expmask) >> 52; | |
43 | ||
44 | if(fp.exp) { | |
45 | fp.frac += hiddenbit; | |
46 | fp.exp -= expbias; | |
47 | ||
48 | } else { | |
49 | fp.exp = -expbias + 1; | |
50 | } | |
51 | ||
52 | return fp; | |
53 | } | |
54 | ||
55 | static void normalize(Fp* fp) | |
56 | { | |
57 | while ((fp->frac & hiddenbit) == 0) { | |
58 | fp->frac <<= 1; | |
59 | fp->exp--; | |
60 | } | |
61 | ||
62 | int shift = 64 - 52 - 1; | |
63 | fp->frac <<= shift; | |
64 | fp->exp -= shift; | |
65 | } | |
66 | ||
67 | static void get_normalized_boundaries(Fp* fp, Fp* lower, Fp* upper) | |
68 | { | |
69 | upper->frac = (fp->frac << 1) + 1; | |
70 | upper->exp = fp->exp - 1; | |
71 | ||
72 | while ((upper->frac & (hiddenbit << 1)) == 0) { | |
73 | upper->frac <<= 1; | |
74 | upper->exp--; | |
75 | } | |
76 | ||
77 | int u_shift = 64 - 52 - 2; | |
78 | ||
79 | upper->frac <<= u_shift; | |
80 | upper->exp = upper->exp - u_shift; | |
81 | ||
82 | ||
83 | int l_shift = fp->frac == hiddenbit ? 2 : 1; | |
84 | ||
85 | lower->frac = (fp->frac << l_shift) - 1; | |
86 | lower->exp = fp->exp - l_shift; | |
87 | ||
88 | ||
89 | lower->frac <<= lower->exp - upper->exp; | |
90 | lower->exp = upper->exp; | |
91 | } | |
92 | ||
93 | static Fp multiply(Fp* a, Fp* b) | |
94 | { | |
95 | const uint64_t lomask = 0x00000000FFFFFFFF; | |
96 | ||
97 | uint64_t ah_bl = (a->frac >> 32) * (b->frac & lomask); | |
98 | uint64_t al_bh = (a->frac & lomask) * (b->frac >> 32); | |
99 | uint64_t al_bl = (a->frac & lomask) * (b->frac & lomask); | |
100 | uint64_t ah_bh = (a->frac >> 32) * (b->frac >> 32); | |
101 | ||
102 | uint64_t tmp = (ah_bl & lomask) + (al_bh & lomask) + (al_bl >> 32); | |
103 | /* round up */ | |
104 | tmp += 1U << 31; | |
105 | ||
106 | Fp fp = { | |
107 | ah_bh + (ah_bl >> 32) + (al_bh >> 32) + (tmp >> 32), | |
108 | a->exp + b->exp + 64 | |
109 | }; | |
110 | ||
111 | return fp; | |
112 | } | |
113 | ||
114 | static void round_digit(char* digits, int ndigits, uint64_t delta, uint64_t rem, uint64_t kappa, uint64_t frac) | |
115 | { | |
116 | while (rem < frac && delta - rem >= kappa && | |
117 | (rem + kappa < frac || frac - rem > rem + kappa - frac)) { | |
118 | ||
119 | digits[ndigits - 1]--; | |
120 | rem += kappa; | |
121 | } | |
122 | } | |
123 | ||
124 | static int generate_digits(Fp* fp, Fp* upper, Fp* lower, char* digits, int* K) | |
125 | { | |
126 | uint64_t wfrac = upper->frac - fp->frac; | |
127 | uint64_t delta = upper->frac - lower->frac; | |
128 | ||
129 | Fp one; | |
130 | one.frac = 1ULL << -upper->exp; | |
131 | one.exp = upper->exp; | |
132 | ||
133 | uint64_t part1 = upper->frac >> -one.exp; | |
134 | uint64_t part2 = upper->frac & (one.frac - 1); | |
135 | ||
136 | int idx = 0, kappa = 10; | |
137 | uint64_t* divp; | |
138 | /* 1000000000 */ | |
139 | for(divp = tens + 10; kappa > 0; divp++) { | |
140 | ||
141 | uint64_t div = *divp; | |
142 | unsigned digit = part1 / div; | |
143 | ||
144 | if (digit || idx) { | |
145 | digits[idx++] = digit + '0'; | |
146 | } | |
147 | ||
148 | part1 -= digit * div; | |
149 | kappa--; | |
150 | ||
151 | uint64_t tmp = (part1 <<-one.exp) + part2; | |
152 | if (tmp <= delta) { | |
153 | *K += kappa; | |
154 | round_digit(digits, idx, delta, tmp, div << -one.exp, wfrac); | |
155 | ||
156 | return idx; | |
157 | } | |
158 | } | |
159 | ||
160 | /* 10 */ | |
161 | uint64_t* unit = tens + 18; | |
162 | ||
163 | while(true) { | |
164 | part2 *= 10; | |
165 | delta *= 10; | |
166 | kappa--; | |
167 | ||
168 | unsigned digit = part2 >> -one.exp; | |
169 | if (digit || idx) { | |
170 | digits[idx++] = digit + '0'; | |
171 | } | |
172 | ||
173 | part2 &= one.frac - 1; | |
174 | if (part2 < delta) { | |
175 | *K += kappa; | |
176 | round_digit(digits, idx, delta, part2, one.frac, wfrac * *unit); | |
177 | ||
178 | return idx; | |
179 | } | |
180 | ||
181 | unit--; | |
182 | } | |
183 | } | |
184 | ||
185 | static int grisu2(double d, char* digits, int* K) | |
186 | { | |
187 | Fp w = build_fp(d); | |
188 | ||
189 | Fp lower, upper; | |
190 | get_normalized_boundaries(&w, &lower, &upper); | |
191 | ||
192 | normalize(&w); | |
193 | ||
194 | int k; | |
195 | Fp cp = find_cachedpow10(upper.exp, &k); | |
196 | ||
197 | w = multiply(&w, &cp); | |
198 | upper = multiply(&upper, &cp); | |
199 | lower = multiply(&lower, &cp); | |
200 | ||
201 | lower.frac++; | |
202 | upper.frac--; | |
203 | ||
204 | *K = -k; | |
205 | ||
206 | return generate_digits(&w, &upper, &lower, digits, K); | |
207 | } | |
208 | ||
209 | static int emit_digits(char* digits, int ndigits, char* dest, int K, bool neg) | |
210 | { | |
211 | int exp = absv(K + ndigits - 1); | |
212 | ||
213 | /* write plain integer */ | |
214 | if(K >= 0 && (exp < (ndigits + 7))) { | |
215 | memcpy(dest, digits, ndigits); | |
216 | memset(dest + ndigits, '0', K); | |
217 | ||
218 | return ndigits + K; | |
219 | } | |
220 | ||
221 | /* write decimal w/o scientific notation */ | |
222 | if(K < 0 && (K > -7 || exp < 4)) { | |
223 | int offset = ndigits - absv(K); | |
224 | /* fp < 1.0 -> write leading zero */ | |
225 | if(offset <= 0) { | |
226 | offset = -offset; | |
227 | dest[0] = '0'; | |
228 | dest[1] = '.'; | |
229 | memset(dest + 2, '0', offset); | |
230 | memcpy(dest + offset + 2, digits, ndigits); | |
231 | ||
232 | return ndigits + 2 + offset; | |
233 | ||
234 | /* fp > 1.0 */ | |
235 | } else { | |
236 | memcpy(dest, digits, offset); | |
237 | dest[offset] = '.'; | |
238 | memcpy(dest + offset + 1, digits + offset, ndigits - offset); | |
239 | ||
240 | return ndigits + 1; | |
241 | } | |
242 | } | |
243 | ||
244 | /* write decimal w/ scientific notation */ | |
245 | ndigits = minv(ndigits, 18 - neg); | |
246 | ||
247 | int idx = 0; | |
248 | dest[idx++] = digits[0]; | |
249 | ||
250 | if(ndigits > 1) { | |
251 | dest[idx++] = '.'; | |
252 | memcpy(dest + idx, digits + 1, ndigits - 1); | |
253 | idx += ndigits - 1; | |
254 | } | |
255 | ||
256 | dest[idx++] = 'e'; | |
257 | ||
258 | char sign = K + ndigits - 1 < 0 ? '-' : '+'; | |
259 | dest[idx++] = sign; | |
260 | ||
261 | int cent = 0; | |
262 | ||
263 | if(exp > 99) { | |
264 | cent = exp / 100; | |
265 | dest[idx++] = cent + '0'; | |
266 | exp -= cent * 100; | |
267 | } | |
268 | if(exp > 9) { | |
269 | int dec = exp / 10; | |
270 | dest[idx++] = dec + '0'; | |
271 | exp -= dec * 10; | |
272 | ||
273 | } else if(cent) { | |
274 | dest[idx++] = '0'; | |
275 | } | |
276 | ||
277 | dest[idx++] = exp % 10 + '0'; | |
278 | ||
279 | return idx; | |
280 | } | |
281 | ||
282 | static int filter_special(double fp, char* dest) | |
283 | { | |
284 | if(fp == 0.0) { | |
285 | dest[0] = '0'; | |
286 | return 1; | |
287 | } | |
288 | ||
289 | uint64_t bits = get_dbits(fp); | |
290 | ||
291 | bool nan = (bits & expmask) == expmask; | |
292 | ||
293 | if(!nan) { | |
294 | return 0; | |
295 | } | |
296 | ||
297 | if(bits & fracmask) { | |
298 | dest[0] = 'n'; dest[1] = 'a'; dest[2] = 'n'; | |
299 | ||
300 | } else { | |
301 | dest[0] = 'i'; dest[1] = 'n'; dest[2] = 'f'; | |
302 | } | |
303 | ||
304 | return 3; | |
305 | } | |
306 | ||
307 | int fpconv_dtoa(double d, char dest[24]) | |
308 | { | |
309 | char digits[18]; | |
310 | ||
311 | int str_len = 0; | |
312 | bool neg = false; | |
313 | ||
314 | if(get_dbits(d) & signmask) { | |
315 | dest[0] = '-'; | |
316 | str_len++; | |
317 | neg = true; | |
318 | } | |
319 | ||
320 | int spec = filter_special(d, dest + str_len); | |
321 | ||
322 | if(spec) { | |
323 | return str_len + spec; | |
324 | } | |
325 | ||
326 | int K = 0; | |
327 | int ndigits = grisu2(d, digits, &K); | |
328 | ||
329 | str_len += emit_digits(digits, ndigits, dest + str_len, K, neg); | |
330 | ||
331 | return str_len; | |
332 | } |