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