1 /*
2 * Copyright 2018-2020 Raffaello Giulietti
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 * THE SOFTWARE.
21 */
22
23 package jdk.internal.math;
24
25 import java.io.IOException;
26
27 import static java.lang.Double.*;
28 import static java.lang.Long.*;
29 import static java.lang.Math.multiplyHigh;
30 import static jdk.internal.math.MathUtils.*;
31
32 /**
33 * This class exposes a method to render a {@code double} as a string.
34 *
35 * @author Raffaello Giulietti
36 */
37 final public class DoubleToDecimal {
38 /*
39 For full details about this code see the following references:
40
41 [1] Giulietti, "The Schubfach way to render doubles",
42 https://drive.google.com/open?id=1luHhyQF9zKlM8yJ1nebU0OgVYhfC6CBN
43
44 [2] IEEE Computer Society, "IEEE Standard for Floating-Point Arithmetic"
45
46 [3] Bouvier & Zimmermann, "Division-Free Binary-to-Decimal Conversion"
47
48 Divisions are avoided altogether for the benefit of those architectures
49 that do not provide specific machine instructions or where they are slow.
50 This is discussed in section 10 of [1].
51 */
52
53 // The precision in bits.
54 static final int P = 53;
55
56 // Exponent width in bits.
57 private static final int W = (Double.SIZE - 1) - (P - 1);
58
59 // Minimum value of the exponent: -(2^(W-1)) - P + 3.
60 static final int Q_MIN = (-1 << W - 1) - P + 3;
61
62 // Maximum value of the exponent: 2^(W-1) - P.
63 static final int Q_MAX = (1 << W - 1) - P;
64
65 // 10^(E_MIN - 1) <= MIN_VALUE < 10^E_MIN
66 static final int E_MIN = -323;
67
68 // 10^(E_MAX - 1) <= MAX_VALUE < 10^E_MAX
69 static final int E_MAX = 309;
70
71 // Threshold to detect tiny values, as in section 8.1.1 of [1]
72 static final long C_TINY = 3;
73
74 // The minimum and maximum k, as in section 8 of [1]
75 static final int K_MIN = -324;
76 static final int K_MAX = 292;
77
78 // H is as in section 8 of [1].
79 static final int H = 17;
80
81 // Minimum value of the significand of a normal value: 2^(P-1).
82 private static final long C_MIN = 1L << P - 1;
83
84 // Mask to extract the biased exponent.
85 private static final int BQ_MASK = (1 << W) - 1;
86
87 // Mask to extract the fraction bits.
88 private static final long T_MASK = (1L << P - 1) - 1;
89
90 // Used in rop().
91 private static final long MASK_63 = (1L << 63) - 1;
92
93 // Used for left-to-tight digit extraction.
94 private static final int MASK_28 = (1 << 28) - 1;
95
96 private static final int NON_SPECIAL = 0;
97 private static final int PLUS_ZERO = 1;
98 private static final int MINUS_ZERO = 2;
99 private static final int PLUS_INF = 3;
100 private static final int MINUS_INF = 4;
101 private static final int NAN = 5;
102
103 // For thread-safety, each thread gets its own instance of this class.
104 private static final ThreadLocal<DoubleToDecimal> threadLocal =
105 ThreadLocal.withInitial(DoubleToDecimal::new);
106
107 /*
108 Room for the longer of the forms
109 -ddddd.dddddddddddd H + 2 characters
110 -0.00ddddddddddddddddd H + 5 characters
111 -d.ddddddddddddddddE-eee H + 7 characters
112 where there are H digits d
113 */
114 public final int MAX_CHARS = H + 7;
115
116 // Numerical results are created here...
117 private final byte[] bytes = new byte[MAX_CHARS];
118
119 // ... and copied here in appendTo()
120 private final char[] chars = new char[MAX_CHARS];
121
122 // Index into bytes of rightmost valid character.
123 private int index;
124
125 private DoubleToDecimal() {
126 }
127
128 /**
129 * Returns a string rendering of the {@code double} argument.
130 *
131 * <p>The characters of the result are all drawn from the ASCII set.
132 * <ul>
133 * <li> Any NaN, whether quiet or signaling, is rendered as
134 * {@code "NaN"}, regardless of the sign bit.
135 * <li> The infinities +∞ and -∞ are rendered as
136 * {@code "Infinity"} and {@code "-Infinity"}, respectively.
137 * <li> The positive and negative zeroes are rendered as
138 * {@code "0.0"} and {@code "-0.0"}, respectively.
139 * <li> A finite negative {@code v} is rendered as the sign
140 * '{@code -}' followed by the rendering of the magnitude -{@code v}.
141 * <li> A finite positive {@code v} is rendered in two stages:
142 * <ul>
143 * <li> <em>Selection of a decimal</em>: A well-defined
144 * decimal <i>d</i><sub><code>v</code></sub> is selected
145 * to represent {@code v}.
146 * <li> <em>Formatting as a string</em>: The decimal
147 * <i>d</i><sub><code>v</code></sub> is formatted as a string,
148 * either in plain or in computerized scientific notation,
149 * depending on its value.
150 * </ul>
151 * </ul>
152 *
153 * <p>A <em>decimal</em> is a number of the form
154 * <i>d</i>×10<sup><i>i</i></sup>
155 * for some (unique) integers <i>d</i> > 0 and <i>i</i> such that
156 * <i>d</i> is not a multiple of 10.
157 * These integers are the <em>significand</em> and
158 * the <em>exponent</em>, respectively, of the decimal.
159 * The <em>length</em> of the decimal is the (unique)
160 * integer <i>n</i> meeting
161 * 10<sup><i>n</i>-1</sup> ≤ <i>d</i> < 10<sup><i>n</i></sup>.
162 *
163 * <p>The decimal <i>d</i><sub><code>v</code></sub>
164 * for a finite positive {@code v} is defined as follows:
165 * <ul>
166 * <li>Let <i>R</i> be the set of all decimals that round to {@code v}
167 * according to the usual round-to-closest rule of
168 * IEEE 754 floating-point arithmetic.
169 * <li>Let <i>m</i> be the minimal length over all decimals in <i>R</i>.
170 * <li>When <i>m</i> ≥ 2, let <i>T</i> be the set of all decimals
171 * in <i>R</i> with length <i>m</i>.
172 * Otherwise, let <i>T</i> be the set of all decimals
173 * in <i>R</i> with length 1 or 2.
174 * <li>Define <i>d</i><sub><code>v</code></sub> as
175 * the decimal in <i>T</i> that is closest to {@code v}.
176 * Or if there are two such decimals in <i>T</i>,
177 * select the one with the even significand (there is exactly one).
178 * </ul>
179 *
180 * <p>The (uniquely) selected decimal <i>d</i><sub><code>v</code></sub>
181 * is then formatted.
182 *
183 * <p>Let <i>d</i>, <i>i</i> and <i>n</i> be the significand, exponent and
184 * length of <i>d</i><sub><code>v</code></sub>, respectively.
185 * Further, let <i>e</i> = <i>n</i> + <i>i</i> - 1 and let
186 * <i>d</i><sub>1</sub>…<i>d</i><sub><i>n</i></sub>
187 * be the usual decimal expansion of the significand.
188 * Note that <i>d</i><sub>1</sub> ≠ 0 ≠ <i>d</i><sub><i>n</i></sub>.
189 * <ul>
190 * <li>Case -3 ≤ <i>e</i> < 0:
191 * <i>d</i><sub><code>v</code></sub> is formatted as
192 * <code>0.0</code>…<code>0</code><!--
193 * --><i>d</i><sub>1</sub>…<i>d</i><sub><i>n</i></sub>,
194 * where there are exactly -(<i>n</i> + <i>i</i>) zeroes between
195 * the decimal point and <i>d</i><sub>1</sub>.
196 * For example, 123 × 10<sup>-4</sup> is formatted as
197 * {@code 0.0123}.
198 * <li>Case 0 ≤ <i>e</i> < 7:
199 * <ul>
200 * <li>Subcase <i>i</i> ≥ 0:
201 * <i>d</i><sub><code>v</code></sub> is formatted as
202 * <i>d</i><sub>1</sub>…<i>d</i><sub><i>n</i></sub><!--
203 * --><code>0</code>…<code>0.0</code>,
204 * where there are exactly <i>i</i> zeroes
205 * between <i>d</i><sub><i>n</i></sub> and the decimal point.
206 * For example, 123 × 10<sup>2</sup> is formatted as
207 * {@code 12300.0}.
208 * <li>Subcase <i>i</i> < 0:
209 * <i>d</i><sub><code>v</code></sub> is formatted as
210 * <i>d</i><sub>1</sub>…<!--
211 * --><i>d</i><sub><i>n</i>+<i>i</i></sub>.<!--
212 * --><i>d</i><sub><i>n</i>+<i>i</i>+1</sub>…<!--
213 * --><i>d</i><sub><i>n</i></sub>.
214 * There are exactly -<i>i</i> digits to the right of
215 * the decimal point.
216 * For example, 123 × 10<sup>-1</sup> is formatted as
217 * {@code 12.3}.
218 * </ul>
219 * <li>Case <i>e</i> < -3 or <i>e</i> ≥ 7:
220 * computerized scientific notation is used to format
221 * <i>d</i><sub><code>v</code></sub>.
222 * Here <i>e</i> is formatted as by {@link Integer#toString(int)}.
223 * <ul>
224 * <li>Subcase <i>n</i> = 1:
225 * <i>d</i><sub><code>v</code></sub> is formatted as
226 * <i>d</i><sub>1</sub><code>.0E</code><i>e</i>.
227 * For example, 1 × 10<sup>23</sup> is formatted as
228 * {@code 1.0E23}.
229 * <li>Subcase <i>n</i> > 1:
230 * <i>d</i><sub><code>v</code></sub> is formatted as
231 * <i>d</i><sub>1</sub><code>.</code><i>d</i><sub>2</sub><!--
232 * -->…<i>d</i><sub><i>n</i></sub><code>E</code><i>e</i>.
233 * For example, 123 × 10<sup>-21</sup> is formatted as
234 * {@code 1.23E-19}.
235 * </ul>
236 * </ul>
237 *
238 * @param v the {@code double} to be rendered.
239 * @return a string rendering of the argument.
240 */
241 public static String toString(double v) {
242 return threadLocalInstance().toDecimalString(v);
243 }
244
245 /**
246 * Appends the rendering of the {@code v} to {@code app}.
247 *
248 * <p>The outcome is the same as if {@code v} were first
249 * {@link #toString(double) rendered} and the resulting string were then
250 * {@link Appendable#append(CharSequence) appended} to {@code app}.
251 *
252 * @param v the {@code double} whose rendering is appended.
253 * @param app the {@link Appendable} to append to.
254 * @throws IOException If an I/O error occurs
255 */
256 public static Appendable appendTo(double v, Appendable app)
257 throws IOException {
258 return threadLocalInstance().appendDecimalTo(v, app);
259 }
260
261 private static DoubleToDecimal threadLocalInstance() {
262 return threadLocal.get();
263 }
264
265 private String toDecimalString(double v) {
266 switch (toDecimal(v)) {
267 case NON_SPECIAL: return charsToString();
268 case PLUS_ZERO: return "0.0";
269 case MINUS_ZERO: return "-0.0";
270 case PLUS_INF: return "Infinity";
271 case MINUS_INF: return "-Infinity";
272 default: return "NaN";
273 }
274 }
275
276 private Appendable appendDecimalTo(double v, Appendable app)
277 throws IOException {
278 switch (toDecimal(v)) {
279 case NON_SPECIAL:
280 for (int i = 0; i <= index; ++i) {
281 chars[i] = (char) bytes[i];
282 }
283 if (app instanceof StringBuilder) {
284 return ((StringBuilder) app).append(chars, 0, index + 1);
285 }
286 if (app instanceof StringBuffer) {
287 return ((StringBuffer) app).append(chars, 0, index + 1);
288 }
289 for (int i = 0; i <= index; ++i) {
290 app.append(chars[i]);
291 }
292 return app;
293 case PLUS_ZERO: return app.append("0.0");
294 case MINUS_ZERO: return app.append("-0.0");
295 case PLUS_INF: return app.append("Infinity");
296 case MINUS_INF: return app.append("-Infinity");
297 default: return app.append("NaN");
298 }
299 }
300
301 /*
302 Returns
303 PLUS_ZERO iff v is 0.0
304 MINUS_ZERO iff v is -0.0
305 PLUS_INF iff v is POSITIVE_INFINITY
306 MINUS_INF iff v is NEGATIVE_INFINITY
307 NAN iff v is NaN
308 */
309 private int toDecimal(double v) {
310 /*
311 For full details see references [2] and [1].
312
313 For finite v != 0, determine integers c and q such that
314 |v| = c 2^q and
315 Q_MIN <= q <= Q_MAX and
316 either 2^(P-1) <= c < 2^P (normal)
317 or 0 < c < 2^(P-1) and q = Q_MIN (subnormal)
318 */
319 long bits = doubleToRawLongBits(v);
320 long t = bits & T_MASK;
321 int bq = (int) (bits >>> P - 1) & BQ_MASK;
322 if (bq < BQ_MASK) {
323 index = -1;
324 if (bits < 0) {
325 append('-');
326 }
327 if (bq != 0) {
328 // normal value. Here mq = -q
329 int mq = -Q_MIN + 1 - bq;
330 long c = C_MIN | t;
331 // The fast path discussed in section 8.2 of [1].
332 if (0 < mq & mq < P) {
333 long f = c >> mq;
334 if (f << mq == c) {
335 return toChars(f, 0);
336 }
337 }
338 return toDecimal(-mq, c, 0);
339 }
340 if (t != 0) {
341 // subnormal value
342 return t < C_TINY
343 ? toDecimal(Q_MIN, 10 * t, -1)
344 : toDecimal(Q_MIN, t, 0);
345 }
346 return bits == 0 ? PLUS_ZERO : MINUS_ZERO;
347 }
348 if (t != 0) {
349 return NAN;
350 }
351 return bits > 0 ? PLUS_INF : MINUS_INF;
352 }
353
354 private int toDecimal(int q, long c, int dk) {
355 /*
356 The skeleton corresponds to figure 4 of [1].
357 The efficient computations are those summarized in figure 7.
358
359 Here's a correspondence between Java names and names in [1],
360 expressed as approximate LaTeX source code and informally.
361 Other names are identical.
362 cb: \bar{c} "c-bar"
363 cbr: \bar{c}_r "c-bar-r"
364 cbl: \bar{c}_l "c-bar-l"
365
366 vb: \bar{v} "v-bar"
367 vbr: \bar{v}_r "v-bar-r"
368 vbl: \bar{v}_l "v-bar-l"
369
370 rop: r_o' "r-o-prime"
371 */
372 int out = (int) c & 0x1;
373 long cb = c << 2;
374 long cbr = cb + 2;
375 long cbl;
376 int k;
377 /*
378 flog10pow2(e) = floor(log_10(2^e))
379 flog10threeQuartersPow2(e) = floor(log_10(3/4 2^e))
380 flog2pow10(e) = floor(log_2(10^e))
381 */
382 if (c != C_MIN | q == Q_MIN) {
383 // regular spacing
384 cbl = cb - 2;
385 k = flog10pow2(q);
386 } else {
387 // irregular spacing
388 cbl = cb - 1;
389 k = flog10threeQuartersPow2(q);
390 }
391 int h = q + flog2pow10(-k) + 2;
392
393 // g1 and g0 are as in section 9.9.3 of [1], so g = g1 2^63 + g0
394 long g1 = g1(k);
395 long g0 = g0(k);
396
397 long vb = rop(g1, g0, cb << h);
398 long vbl = rop(g1, g0, cbl << h);
399 long vbr = rop(g1, g0, cbr << h);
400
401 long s = vb >> 2;
402 if (s >= 100) {
403 /*
404 For n = 17, m = 1 the table in section 10 of [1] shows
405 s' = floor(s / 10) = floor(s 115_292_150_460_684_698 / 2^60)
406 = floor(s 115_292_150_460_684_698 2^4 / 2^64)
407
408 sp10 = 10 s'
409 tp10 = 10 t'
410 upin iff u' = sp10 10^k in Rv
411 wpin iff w' = tp10 10^k in Rv
412 See section 9.4 of [1].
413 */
414 long sp10 = 10 * multiplyHigh(s, 115_292_150_460_684_698L << 4);
415 long tp10 = sp10 + 10;
416 boolean upin = vbl + out <= sp10 << 2;
417 boolean wpin = (tp10 << 2) + out <= vbr;
418 if (upin != wpin) {
419 return toChars(upin ? sp10 : tp10, k);
420 }
421 }
422
423 /*
424 10 <= s < 100 or s >= 100 and u', w' not in Rv
425 uin iff u = s 10^k in Rv
426 win iff w = t 10^k in Rv
427 See section 9.4 of [1].
428 */
429 long t = s + 1;
430 boolean uin = vbl + out <= s << 2;
431 boolean win = (t << 2) + out <= vbr;
432 if (uin != win) {
433 // Exactly one of u or w lies in Rv.
434 return toChars(uin ? s : t, k + dk);
435 }
436 /*
437 Both u and w lie in Rv: determine the one closest to v.
438 See section 9.4 of [1].
439 */
440 long cmp = vb - (s + t << 1);
441 return toChars(cmp < 0 || cmp == 0 && (s & 0x1) == 0 ? s : t, k + dk);
442 }
443
444 /*
445 Computes rop(cp g 2^(-127)), where g = g1 2^63 + g0
446 See section 9.10 and figure 5 of [1].
447 */
448 private static long rop(long g1, long g0, long cp) {
449 long x1 = multiplyHigh(g0, cp);
450 long y0 = g1 * cp;
451 long y1 = multiplyHigh(g1, cp);
452 long z = (y0 >>> 1) + x1;
453 long vbp = y1 + (z >>> 63);
454 return vbp | (z & MASK_63) + MASK_63 >>> 63;
455 }
456
457 /*
458 Formats the decimal f 10^e.
459 */
460 private int toChars(long f, int e) {
461 /*
462 For details not discussed here see section 10 of [1].
463
464 Determine len such that
465 10^(len-1) <= f < 10^len
466 */
467 int len = flog10pow2(Long.SIZE - numberOfLeadingZeros(f));
468 if (f >= pow10(len)) {
469 len += 1;
470 }
471
472 /*
473 Let fp and ep be the original f and e, respectively.
474 Transform f and e to ensure
475 10^(H-1) <= f < 10^H
476 fp 10^ep = f 10^(e-H) = 0.f 10^e
477 */
478 f *= pow10(H - len);
479 e += len;
480
481 /*
482 The toChars?() methods perform left-to-right digits extraction
483 using ints, provided that the arguments are limited to 8 digits.
484 Therefore, split the H = 17 digits of f into:
485 h = the most significant digit of f
486 m = the next 8 most significant digits of f
487 l = the last 8, least significant digits of f
488
489 For n = 17, m = 8 the table in section 10 of [1] shows
490 floor(f / 10^8) = floor(193_428_131_138_340_668 f / 2^84) =
491 floor(floor(193_428_131_138_340_668 f / 2^64) / 2^20)
492 and for n = 9, m = 8
493 floor(hm / 10^8) = floor(1_441_151_881 hm / 2^57)
494 */
495 long hm = multiplyHigh(f, 193_428_131_138_340_668L) >>> 20;
496 int l = (int) (f - 100_000_000L * hm);
497 int h = (int) (hm * 1_441_151_881L >>> 57);
498 int m = (int) (hm - 100_000_000 * h);
499
500 if (0 < e && e <= 7) {
501 return toChars1(h, m, l, e);
502 }
503 if (-3 < e && e <= 0) {
504 return toChars2(h, m, l, e);
505 }
506 return toChars3(h, m, l, e);
507 }
508
509 private int toChars1(int h, int m, int l, int e) {
510 /*
511 0 < e <= 7: plain format without leading zeroes.
512 Left-to-right digits extraction:
513 algorithm 1 in [3], with b = 10, k = 8, n = 28.
514 */
515 appendDigit(h);
516 int y = y(m);
517 int t;
518 int i = 1;
519 for (; i < e; ++i) {
520 t = 10 * y;
521 appendDigit(t >>> 28);
522 y = t & MASK_28;
523 }
524 append('.');
525 for (; i <= 8; ++i) {
526 t = 10 * y;
527 appendDigit(t >>> 28);
528 y = t & MASK_28;
529 }
530 lowDigits(l);
531 return NON_SPECIAL;
532 }
533
534 private int toChars2(int h, int m, int l, int e) {
535 // -3 < e <= 0: plain format with leading zeroes.
536 appendDigit(0);
537 append('.');
538 for (; e < 0; ++e) {
539 appendDigit(0);
540 }
541 appendDigit(h);
542 append8Digits(m);
543 lowDigits(l);
544 return NON_SPECIAL;
545 }
546
547 private int toChars3(int h, int m, int l, int e) {
548 // -3 >= e | e > 7: computerized scientific notation
549 appendDigit(h);
550 append('.');
551 append8Digits(m);
552 lowDigits(l);
553 exponent(e - 1);
554 return NON_SPECIAL;
555 }
556
557 private void lowDigits(int l) {
558 if (l != 0) {
559 append8Digits(l);
560 }
561 removeTrailingZeroes();
562 }
563
564 private void append8Digits(int m) {
565 /*
566 Left-to-right digits extraction:
567 algorithm 1 in [3], with b = 10, k = 8, n = 28.
568 */
569 int y = y(m);
570 for (int i = 0; i < 8; ++i) {
571 int t = 10 * y;
572 appendDigit(t >>> 28);
573 y = t & MASK_28;
574 }
575 }
576
577 private void removeTrailingZeroes() {
578 while (bytes[index] == '0') {
579 --index;
580 }
581 // ... but do not remove the one directly to the right of '.'
582 if (bytes[index] == '.') {
583 ++index;
584 }
585 }
586
587 private int y(int a) {
588 /*
589 Algorithm 1 in [3] needs computation of
590 floor((a + 1) 2^n / b^k) - 1
591 with a < 10^8, b = 10, k = 8, n = 28.
592 Noting that
593 (a + 1) 2^n <= 10^8 2^28 < 10^17
594 For n = 17, m = 8 the table in section 10 of [1] leads to:
595 */
596 return (int) (multiplyHigh(
597 (long) (a + 1) << 28,
598 193_428_131_138_340_668L) >>> 20) - 1;
599 }
600
601 private void exponent(int e) {
602 append('E');
603 if (e < 0) {
604 append('-');
605 e = -e;
606 }
607 if (e < 10) {
608 appendDigit(e);
609 return;
610 }
611 int d;
612 if (e >= 100) {
613 /*
614 For n = 3, m = 2 the table in section 10 of [1] shows
615 floor(e / 100) = floor(1_311 e / 2^17)
616 */
617 d = e * 1_311 >>> 17;
618 appendDigit(d);
619 e -= 100 * d;
620 }
621 /*
622 For n = 2, m = 1 the table in section 10 of [1] shows
623 floor(e / 10) = floor(103 e / 2^10)
624 */
625 d = e * 103 >>> 10;
626 appendDigit(d);
627 appendDigit(e - 10 * d);
628 }
629
630 private void append(int c) {
631 bytes[++index] = (byte) c;
632 }
633
634 private void appendDigit(int d) {
635 bytes[++index] = (byte) ('0' + d);
636 }
637
638 // Using the deprecated constructor enhances performance.
639 @SuppressWarnings("deprecation")
640 private String charsToString() {
641 return new String(bytes, 0, 0, index + 1);
642 }
643
644 }