/src/ruby/ext/json/vendor/fpconv.c
Line | Count | Source |
1 | | // Boost Software License - Version 1.0 - August 17th, 2003 |
2 | | // |
3 | | // Permission is hereby granted, free of charge, to any person or organization |
4 | | // obtaining a copy of the software and accompanying documentation covered by |
5 | | // this license (the "Software") to use, reproduce, display, distribute, |
6 | | // execute, and transmit the Software, and to prepare derivative works of the |
7 | | // Software, and to permit third-parties to whom the Software is furnished to |
8 | | // do so, all subject to the following: |
9 | | // |
10 | | // The copyright notices in the Software and this entire statement, including |
11 | | // the above license grant, this restriction and the following disclaimer, |
12 | | // must be included in all copies of the Software, in whole or in part, and |
13 | | // all derivative works of the Software, unless such copies or derivative |
14 | | // works are solely in the form of machine-executable object code generated by |
15 | | // a source language processor. |
16 | | // |
17 | | // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
18 | | // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
19 | | // FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT |
20 | | // SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE |
21 | | // FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, |
22 | | // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER |
23 | | // DEALINGS IN THE SOFTWARE. |
24 | | |
25 | | // The contents of this file is extracted from https://github.com/night-shift/fpconv |
26 | | // It was slightly modified to append ".0" to plain floats, for use with the https://github.com/ruby/json package. |
27 | | |
28 | | #include <stdbool.h> |
29 | | #include <string.h> |
30 | | #include <stdint.h> |
31 | | |
32 | | #if JSON_DEBUG |
33 | | #include <assert.h> |
34 | | #endif |
35 | | |
36 | | #define npowers 87 |
37 | | #define steppowers 8 |
38 | | #define firstpower -348 /* 10 ^ -348 */ |
39 | | |
40 | | #define expmax -32 |
41 | | #define expmin -60 |
42 | | |
43 | | typedef struct Fp { |
44 | | uint64_t frac; |
45 | | int exp; |
46 | | } Fp; |
47 | | |
48 | | static const Fp powers_ten[] = { |
49 | | { 18054884314459144840U, -1220 }, { 13451937075301367670U, -1193 }, |
50 | | { 10022474136428063862U, -1166 }, { 14934650266808366570U, -1140 }, |
51 | | { 11127181549972568877U, -1113 }, { 16580792590934885855U, -1087 }, |
52 | | { 12353653155963782858U, -1060 }, { 18408377700990114895U, -1034 }, |
53 | | { 13715310171984221708U, -1007 }, { 10218702384817765436U, -980 }, |
54 | | { 15227053142812498563U, -954 }, { 11345038669416679861U, -927 }, |
55 | | { 16905424996341287883U, -901 }, { 12595523146049147757U, -874 }, |
56 | | { 9384396036005875287U, -847 }, { 13983839803942852151U, -821 }, |
57 | | { 10418772551374772303U, -794 }, { 15525180923007089351U, -768 }, |
58 | | { 11567161174868858868U, -741 }, { 17236413322193710309U, -715 }, |
59 | | { 12842128665889583758U, -688 }, { 9568131466127621947U, -661 }, |
60 | | { 14257626930069360058U, -635 }, { 10622759856335341974U, -608 }, |
61 | | { 15829145694278690180U, -582 }, { 11793632577567316726U, -555 }, |
62 | | { 17573882009934360870U, -529 }, { 13093562431584567480U, -502 }, |
63 | | { 9755464219737475723U, -475 }, { 14536774485912137811U, -449 }, |
64 | | { 10830740992659433045U, -422 }, { 16139061738043178685U, -396 }, |
65 | | { 12024538023802026127U, -369 }, { 17917957937422433684U, -343 }, |
66 | | { 13349918974505688015U, -316 }, { 9946464728195732843U, -289 }, |
67 | | { 14821387422376473014U, -263 }, { 11042794154864902060U, -236 }, |
68 | | { 16455045573212060422U, -210 }, { 12259964326927110867U, -183 }, |
69 | | { 18268770466636286478U, -157 }, { 13611294676837538539U, -130 }, |
70 | | { 10141204801825835212U, -103 }, { 15111572745182864684U, -77 }, |
71 | | { 11258999068426240000U, -50 }, { 16777216000000000000U, -24 }, |
72 | | { 12500000000000000000U, 3 }, { 9313225746154785156U, 30 }, |
73 | | { 13877787807814456755U, 56 }, { 10339757656912845936U, 83 }, |
74 | | { 15407439555097886824U, 109 }, { 11479437019748901445U, 136 }, |
75 | | { 17105694144590052135U, 162 }, { 12744735289059618216U, 189 }, |
76 | | { 9495567745759798747U, 216 }, { 14149498560666738074U, 242 }, |
77 | | { 10542197943230523224U, 269 }, { 15709099088952724970U, 295 }, |
78 | | { 11704190886730495818U, 322 }, { 17440603504673385349U, 348 }, |
79 | | { 12994262207056124023U, 375 }, { 9681479787123295682U, 402 }, |
80 | | { 14426529090290212157U, 428 }, { 10748601772107342003U, 455 }, |
81 | | { 16016664761464807395U, 481 }, { 11933345169920330789U, 508 }, |
82 | | { 17782069995880619868U, 534 }, { 13248674568444952270U, 561 }, |
83 | | { 9871031767461413346U, 588 }, { 14708983551653345445U, 614 }, |
84 | | { 10959046745042015199U, 641 }, { 16330252207878254650U, 667 }, |
85 | | { 12166986024289022870U, 694 }, { 18130221999122236476U, 720 }, |
86 | | { 13508068024458167312U, 747 }, { 10064294952495520794U, 774 }, |
87 | | { 14996968138956309548U, 800 }, { 11173611982879273257U, 827 }, |
88 | | { 16649979327439178909U, 853 }, { 12405201291620119593U, 880 }, |
89 | | { 9242595204427927429U, 907 }, { 13772540099066387757U, 933 }, |
90 | | { 10261342003245940623U, 960 }, { 15290591125556738113U, 986 }, |
91 | | { 11392378155556871081U, 1013 }, { 16975966327722178521U, 1039 }, |
92 | | { 12648080533535911531U, 1066 } |
93 | | }; |
94 | | |
95 | | static Fp find_cachedpow10(int exp, int* k) |
96 | 0 | { |
97 | 0 | const double one_log_ten = 0.30102999566398114; |
98 | 0 |
|
99 | 0 | int approx = (int)(-(exp + npowers) * one_log_ten); |
100 | 0 | int idx = (approx - firstpower) / steppowers; |
101 | 0 |
|
102 | 0 | while(1) { |
103 | 0 | int current = exp + powers_ten[idx].exp + 64; |
104 | 0 |
|
105 | 0 | if(current < expmin) { |
106 | 0 | idx++; |
107 | 0 | continue; |
108 | 0 | } |
109 | 0 |
|
110 | 0 | if(current > expmax) { |
111 | 0 | idx--; |
112 | 0 | continue; |
113 | 0 | } |
114 | 0 |
|
115 | 0 | *k = (firstpower + idx * steppowers); |
116 | 0 |
|
117 | 0 | return powers_ten[idx]; |
118 | 0 | } |
119 | 0 | } |
120 | | |
121 | | #define fracmask 0x000FFFFFFFFFFFFFU |
122 | | #define expmask 0x7FF0000000000000U |
123 | | #define hiddenbit 0x0010000000000000U |
124 | | #define signmask 0x8000000000000000U |
125 | | #define expbias (1023 + 52) |
126 | | |
127 | | #define absv(n) ((n) < 0 ? -(n) : (n)) |
128 | | #define minv(a, b) ((a) < (b) ? (a) : (b)) |
129 | | |
130 | | static const uint64_t tens[] = { |
131 | | 10000000000000000000U, 1000000000000000000U, 100000000000000000U, |
132 | | 10000000000000000U, 1000000000000000U, 100000000000000U, |
133 | | 10000000000000U, 1000000000000U, 100000000000U, |
134 | | 10000000000U, 1000000000U, 100000000U, |
135 | | 10000000U, 1000000U, 100000U, |
136 | | 10000U, 1000U, 100U, |
137 | | 10U, 1U |
138 | | }; |
139 | | |
140 | | static inline uint64_t get_dbits(double d) |
141 | 0 | { |
142 | 0 | union { |
143 | 0 | double dbl; |
144 | 0 | uint64_t i; |
145 | 0 | } dbl_bits = { d }; |
146 | 0 |
|
147 | 0 | return dbl_bits.i; |
148 | 0 | } |
149 | | |
150 | | static Fp build_fp(double d) |
151 | 0 | { |
152 | 0 | uint64_t bits = get_dbits(d); |
153 | 0 |
|
154 | 0 | Fp fp; |
155 | 0 | fp.frac = bits & fracmask; |
156 | 0 | fp.exp = (bits & expmask) >> 52; |
157 | 0 |
|
158 | 0 | if(fp.exp) { |
159 | 0 | fp.frac += hiddenbit; |
160 | 0 | fp.exp -= expbias; |
161 | 0 |
|
162 | 0 | } else { |
163 | 0 | fp.exp = -expbias + 1; |
164 | 0 | } |
165 | 0 |
|
166 | 0 | return fp; |
167 | 0 | } |
168 | | |
169 | | static void normalize(Fp* fp) |
170 | 0 | { |
171 | 0 | while ((fp->frac & hiddenbit) == 0) { |
172 | 0 | fp->frac <<= 1; |
173 | 0 | fp->exp--; |
174 | 0 | } |
175 | 0 |
|
176 | 0 | int shift = 64 - 52 - 1; |
177 | 0 | fp->frac <<= shift; |
178 | 0 | fp->exp -= shift; |
179 | 0 | } |
180 | | |
181 | | static void get_normalized_boundaries(Fp* fp, Fp* lower, Fp* upper) |
182 | 0 | { |
183 | 0 | upper->frac = (fp->frac << 1) + 1; |
184 | 0 | upper->exp = fp->exp - 1; |
185 | 0 |
|
186 | 0 | while ((upper->frac & (hiddenbit << 1)) == 0) { |
187 | 0 | upper->frac <<= 1; |
188 | 0 | upper->exp--; |
189 | 0 | } |
190 | 0 |
|
191 | 0 | int u_shift = 64 - 52 - 2; |
192 | 0 |
|
193 | 0 | upper->frac <<= u_shift; |
194 | 0 | upper->exp = upper->exp - u_shift; |
195 | 0 |
|
196 | 0 |
|
197 | 0 | int l_shift = fp->frac == hiddenbit ? 2 : 1; |
198 | 0 |
|
199 | 0 | lower->frac = (fp->frac << l_shift) - 1; |
200 | 0 | lower->exp = fp->exp - l_shift; |
201 | 0 |
|
202 | 0 |
|
203 | 0 | lower->frac <<= lower->exp - upper->exp; |
204 | 0 | lower->exp = upper->exp; |
205 | 0 | } |
206 | | |
207 | | static Fp multiply(Fp* a, Fp* b) |
208 | 0 | { |
209 | 0 | const uint64_t lomask = 0x00000000FFFFFFFF; |
210 | 0 |
|
211 | 0 | uint64_t ah_bl = (a->frac >> 32) * (b->frac & lomask); |
212 | 0 | uint64_t al_bh = (a->frac & lomask) * (b->frac >> 32); |
213 | 0 | uint64_t al_bl = (a->frac & lomask) * (b->frac & lomask); |
214 | 0 | uint64_t ah_bh = (a->frac >> 32) * (b->frac >> 32); |
215 | 0 |
|
216 | 0 | uint64_t tmp = (ah_bl & lomask) + (al_bh & lomask) + (al_bl >> 32); |
217 | 0 | /* round up */ |
218 | 0 | tmp += 1U << 31; |
219 | 0 |
|
220 | 0 | Fp fp = { |
221 | 0 | ah_bh + (ah_bl >> 32) + (al_bh >> 32) + (tmp >> 32), |
222 | 0 | a->exp + b->exp + 64 |
223 | 0 | }; |
224 | 0 |
|
225 | 0 | return fp; |
226 | 0 | } |
227 | | |
228 | | static void round_digit(char* digits, int ndigits, uint64_t delta, uint64_t rem, uint64_t kappa, uint64_t frac) |
229 | 0 | { |
230 | 0 | while (rem < frac && delta - rem >= kappa && |
231 | 0 | (rem + kappa < frac || frac - rem > rem + kappa - frac)) { |
232 | 0 |
|
233 | 0 | digits[ndigits - 1]--; |
234 | 0 | rem += kappa; |
235 | 0 | } |
236 | 0 | } |
237 | | |
238 | | static int generate_digits(Fp* fp, Fp* upper, Fp* lower, char* digits, int* K) |
239 | 0 | { |
240 | 0 | uint64_t wfrac = upper->frac - fp->frac; |
241 | 0 | uint64_t delta = upper->frac - lower->frac; |
242 | 0 |
|
243 | 0 | Fp one; |
244 | 0 | one.frac = 1ULL << -upper->exp; |
245 | 0 | one.exp = upper->exp; |
246 | 0 |
|
247 | 0 | uint64_t part1 = upper->frac >> -one.exp; |
248 | 0 | uint64_t part2 = upper->frac & (one.frac - 1); |
249 | 0 |
|
250 | 0 | int idx = 0, kappa = 10; |
251 | 0 | const uint64_t* divp; |
252 | 0 | /* 1000000000 */ |
253 | 0 | for(divp = tens + 10; kappa > 0; divp++) { |
254 | 0 |
|
255 | 0 | uint64_t div = *divp; |
256 | 0 | unsigned digit = (unsigned) (part1 / div); |
257 | 0 |
|
258 | 0 | if (digit || idx) { |
259 | 0 | digits[idx++] = digit + '0'; |
260 | 0 | } |
261 | 0 |
|
262 | 0 | part1 -= digit * div; |
263 | 0 | kappa--; |
264 | 0 |
|
265 | 0 | uint64_t tmp = (part1 <<-one.exp) + part2; |
266 | 0 | if (tmp <= delta) { |
267 | 0 | *K += kappa; |
268 | 0 | round_digit(digits, idx, delta, tmp, div << -one.exp, wfrac); |
269 | 0 |
|
270 | 0 | return idx; |
271 | 0 | } |
272 | 0 | } |
273 | 0 |
|
274 | 0 | /* 10 */ |
275 | 0 | const uint64_t* unit = tens + 18; |
276 | 0 |
|
277 | 0 | while(true) { |
278 | 0 | part2 *= 10; |
279 | 0 | delta *= 10; |
280 | 0 | kappa--; |
281 | 0 |
|
282 | 0 | unsigned digit = (unsigned) (part2 >> -one.exp); |
283 | 0 | if (digit || idx) { |
284 | 0 | digits[idx++] = digit + '0'; |
285 | 0 | } |
286 | 0 |
|
287 | 0 | part2 &= one.frac - 1; |
288 | 0 | if (part2 < delta) { |
289 | 0 | *K += kappa; |
290 | 0 | round_digit(digits, idx, delta, part2, one.frac, wfrac * *unit); |
291 | 0 |
|
292 | 0 | return idx; |
293 | 0 | } |
294 | 0 |
|
295 | 0 | unit--; |
296 | 0 | } |
297 | 0 | } |
298 | | |
299 | | static int grisu2(double d, char* digits, int* K) |
300 | 0 | { |
301 | 0 | Fp w = build_fp(d); |
302 | 0 |
|
303 | 0 | Fp lower, upper; |
304 | 0 | get_normalized_boundaries(&w, &lower, &upper); |
305 | 0 |
|
306 | 0 | normalize(&w); |
307 | 0 |
|
308 | 0 | int k; |
309 | 0 | Fp cp = find_cachedpow10(upper.exp, &k); |
310 | 0 |
|
311 | 0 | w = multiply(&w, &cp); |
312 | 0 | upper = multiply(&upper, &cp); |
313 | 0 | lower = multiply(&lower, &cp); |
314 | 0 |
|
315 | 0 | lower.frac++; |
316 | 0 | upper.frac--; |
317 | 0 |
|
318 | 0 | *K = -k; |
319 | 0 |
|
320 | 0 | return generate_digits(&w, &upper, &lower, digits, K); |
321 | 0 | } |
322 | | |
323 | | static int emit_digits(char* digits, int ndigits, char* dest, int K, bool neg) |
324 | 0 | { |
325 | 0 | int exp = absv(K + ndigits - 1); |
326 | 0 |
|
327 | 0 | if(K >= 0 && exp < 15) { |
328 | 0 | memcpy(dest, digits, ndigits); |
329 | 0 | memset(dest + ndigits, '0', K); |
330 | 0 |
|
331 | 0 | /* add a .0 to mark this as a float. */ |
332 | 0 | dest[ndigits + K] = '.'; |
333 | 0 | dest[ndigits + K + 1] = '0'; |
334 | 0 |
|
335 | 0 | return ndigits + K + 2; |
336 | 0 | } |
337 | 0 |
|
338 | 0 | /* write decimal w/o scientific notation */ |
339 | 0 | if(K < 0 && (K > -7 || exp < 10)) { |
340 | 0 | int offset = ndigits - absv(K); |
341 | 0 | /* fp < 1.0 -> write leading zero */ |
342 | 0 | if(offset <= 0) { |
343 | 0 | offset = -offset; |
344 | 0 | dest[0] = '0'; |
345 | 0 | dest[1] = '.'; |
346 | 0 | memset(dest + 2, '0', offset); |
347 | 0 | memcpy(dest + offset + 2, digits, ndigits); |
348 | 0 |
|
349 | 0 | return ndigits + 2 + offset; |
350 | 0 |
|
351 | 0 | /* fp > 1.0 */ |
352 | 0 | } else { |
353 | 0 | memcpy(dest, digits, offset); |
354 | 0 | dest[offset] = '.'; |
355 | 0 | memcpy(dest + offset + 1, digits + offset, ndigits - offset); |
356 | 0 |
|
357 | 0 | return ndigits + 1; |
358 | 0 | } |
359 | 0 | } |
360 | 0 |
|
361 | 0 | /* write decimal w/ scientific notation */ |
362 | 0 | ndigits = minv(ndigits, 18 - neg); |
363 | 0 |
|
364 | 0 | int idx = 0; |
365 | 0 | dest[idx++] = digits[0]; |
366 | 0 |
|
367 | 0 | if(ndigits > 1) { |
368 | 0 | dest[idx++] = '.'; |
369 | 0 | memcpy(dest + idx, digits + 1, ndigits - 1); |
370 | 0 | idx += ndigits - 1; |
371 | 0 | } |
372 | 0 |
|
373 | 0 | dest[idx++] = 'e'; |
374 | 0 |
|
375 | 0 | char sign = K + ndigits - 1 < 0 ? '-' : '+'; |
376 | 0 | dest[idx++] = sign; |
377 | 0 |
|
378 | 0 | int cent = 0; |
379 | 0 |
|
380 | 0 | if(exp > 99) { |
381 | 0 | cent = exp / 100; |
382 | 0 | dest[idx++] = cent + '0'; |
383 | 0 | exp -= cent * 100; |
384 | 0 | } |
385 | 0 | if(exp > 9) { |
386 | 0 | int dec = exp / 10; |
387 | 0 | dest[idx++] = dec + '0'; |
388 | 0 | exp -= dec * 10; |
389 | 0 |
|
390 | 0 | } else if(cent) { |
391 | 0 | dest[idx++] = '0'; |
392 | 0 | } |
393 | 0 |
|
394 | 0 | dest[idx++] = exp % 10 + '0'; |
395 | 0 |
|
396 | 0 | return idx; |
397 | 0 | } |
398 | | |
399 | | static int filter_special(double fp, char* dest) |
400 | 0 | { |
401 | 0 | if(fp == 0.0) { |
402 | 0 | dest[0] = '0'; |
403 | 0 | dest[1] = '.'; |
404 | 0 | dest[2] = '0'; |
405 | 0 | return 3; |
406 | 0 | } |
407 | 0 |
|
408 | 0 | uint64_t bits = get_dbits(fp); |
409 | 0 |
|
410 | 0 | bool nan = (bits & expmask) == expmask; |
411 | 0 |
|
412 | 0 | if(!nan) { |
413 | 0 | return 0; |
414 | 0 | } |
415 | 0 |
|
416 | 0 | if(bits & fracmask) { |
417 | 0 | dest[0] = 'n'; dest[1] = 'a'; dest[2] = 'n'; |
418 | 0 |
|
419 | 0 | } else { |
420 | 0 | dest[0] = 'i'; dest[1] = 'n'; dest[2] = 'f'; |
421 | 0 | } |
422 | 0 |
|
423 | 0 | return 3; |
424 | 0 | } |
425 | | |
426 | | /* Fast and accurate double to string conversion based on Florian Loitsch's |
427 | | * Grisu-algorithm[1]. |
428 | | * |
429 | | * Input: |
430 | | * fp -> the double to convert, dest -> destination buffer. |
431 | | * The generated string will never be longer than 32 characters. |
432 | | * Make sure to pass a pointer to at least 32 bytes of memory. |
433 | | * The emitted string will not be null terminated. |
434 | | * |
435 | | * |
436 | | * |
437 | | * Output: |
438 | | * The number of written characters. |
439 | | * |
440 | | * Exemplary usage: |
441 | | * |
442 | | * void print(double d) |
443 | | * { |
444 | | * char buf[28 + 1] // plus null terminator |
445 | | * int str_len = fpconv_dtoa(d, buf); |
446 | | * |
447 | | * buf[str_len] = '\0'; |
448 | | * printf("%s", buf); |
449 | | * } |
450 | | * |
451 | | */ |
452 | | static int fpconv_dtoa(double d, char dest[32]) |
453 | 0 | { |
454 | 0 | char digits[18]; |
455 | 0 |
|
456 | 0 | int str_len = 0; |
457 | 0 | bool neg = false; |
458 | 0 |
|
459 | 0 | if(get_dbits(d) & signmask) { |
460 | 0 | dest[0] = '-'; |
461 | 0 | str_len++; |
462 | 0 | neg = true; |
463 | 0 | } |
464 | 0 |
|
465 | 0 | int spec = filter_special(d, dest + str_len); |
466 | 0 |
|
467 | 0 | if(spec) { |
468 | 0 | return str_len + spec; |
469 | 0 | } |
470 | 0 |
|
471 | 0 | int K = 0; |
472 | 0 | int ndigits = grisu2(d, digits, &K); |
473 | 0 |
|
474 | 0 | str_len += emit_digits(digits, ndigits, dest + str_len, K, neg); |
475 | 0 | #if JSON_DEBUG |
476 | 0 | assert(str_len <= 32); |
477 | 0 | #endif |
478 | 0 |
|
479 | 0 | return str_len; |
480 | 0 | } |