1 /* 2 * Copyright (c) 2003, 2011, Oracle and/or its affiliates. All rights reserved. 3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. 4 * 5 * This code is free software; you can redistribute it and/or modify it 6 * under the terms of the GNU General Public License version 2 only, as 7 * published by the Free Software Foundation. Oracle designates this 8 * particular file as subject to the "Classpath" exception as provided 9 * by Oracle in the LICENSE file that accompanied this code. 10 * 11 * This code is distributed in the hope that it will be useful, but WITHOUT 12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 14 * version 2 for more details (a copy is included in the LICENSE file that 15 * accompanied this code). 16 * 17 * You should have received a copy of the GNU General Public License version 18 * 2 along with this work; if not, write to the Free Software Foundation, 19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 20 * 21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA 22 * or visit www.oracle.com if you need additional information or have any 23 * questions. 24 */ 25 26 package sun.misc; 27 28 import sun.misc.FpUtils; 29 import sun.misc.DoubleConsts; 30 import sun.misc.FloatConsts; 31 import java.util.regex.*; 32 33 public class FormattedFloatingDecimal{ 34 boolean isExceptional; 35 boolean isNegative; 36 int decExponent; // value set at construction, then immutable 37 int decExponentRounded; 38 char digits[]; 39 int nDigits; 40 int bigIntExp; 41 int bigIntNBits; 42 boolean mustSetRoundDir = false; 43 boolean fromHex = false; 44 int roundDir = 0; // set by doubleValue 45 int precision; // number of digits to the right of decimal 46 47 public enum Form { SCIENTIFIC, COMPATIBLE, DECIMAL_FLOAT, GENERAL }; 48 49 private Form form; 50 FormattedFloatingDecimal( boolean negSign, int decExponent, char []digits, int n, boolean e, int precision, Form form )51 private FormattedFloatingDecimal( boolean negSign, int decExponent, char []digits, int n, boolean e, int precision, Form form ) 52 { 53 isNegative = negSign; 54 isExceptional = e; 55 this.decExponent = decExponent; 56 this.digits = digits; 57 this.nDigits = n; 58 this.precision = precision; 59 this.form = form; 60 } 61 62 /* 63 * Constants of the implementation 64 * Most are IEEE-754 related. 65 * (There are more really boring constants at the end.) 66 */ 67 static final long signMask = 0x8000000000000000L; 68 static final long expMask = 0x7ff0000000000000L; 69 static final long fractMask= ~(signMask|expMask); 70 static final int expShift = 52; 71 static final int expBias = 1023; 72 static final long fractHOB = ( 1L<<expShift ); // assumed High-Order bit 73 static final long expOne = ((long)expBias)<<expShift; // exponent of 1.0 74 static final int maxSmallBinExp = 62; 75 static final int minSmallBinExp = -( 63 / 3 ); 76 static final int maxDecimalDigits = 15; 77 static final int maxDecimalExponent = 308; 78 static final int minDecimalExponent = -324; 79 static final int bigDecimalExponent = 324; // i.e. abs(minDecimalExponent) 80 81 static final long highbyte = 0xff00000000000000L; 82 static final long highbit = 0x8000000000000000L; 83 static final long lowbytes = ~highbyte; 84 85 static final int singleSignMask = 0x80000000; 86 static final int singleExpMask = 0x7f800000; 87 static final int singleFractMask = ~(singleSignMask|singleExpMask); 88 static final int singleExpShift = 23; 89 static final int singleFractHOB = 1<<singleExpShift; 90 static final int singleExpBias = 127; 91 static final int singleMaxDecimalDigits = 7; 92 static final int singleMaxDecimalExponent = 38; 93 static final int singleMinDecimalExponent = -45; 94 95 static final int intDecimalDigits = 9; 96 97 98 /* 99 * count number of bits from high-order 1 bit to low-order 1 bit, 100 * inclusive. 101 */ 102 private static int countBits( long v )103 countBits( long v ){ 104 // 105 // the strategy is to shift until we get a non-zero sign bit 106 // then shift until we have no bits left, counting the difference. 107 // we do byte shifting as a hack. Hope it helps. 108 // 109 if ( v == 0L ) return 0; 110 111 while ( ( v & highbyte ) == 0L ){ 112 v <<= 8; 113 } 114 while ( v > 0L ) { // i.e. while ((v&highbit) == 0L ) 115 v <<= 1; 116 } 117 118 int n = 0; 119 while (( v & lowbytes ) != 0L ){ 120 v <<= 8; 121 n += 8; 122 } 123 while ( v != 0L ){ 124 v <<= 1; 125 n += 1; 126 } 127 return n; 128 } 129 130 /* 131 * Keep big powers of 5 handy for future reference. 132 */ 133 private static FDBigInt b5p[]; 134 135 private static synchronized FDBigInt big5pow( int p )136 big5pow( int p ){ 137 assert p >= 0 : p; // negative power of 5 138 if ( b5p == null ){ 139 b5p = new FDBigInt[ p+1 ]; 140 }else if (b5p.length <= p ){ 141 FDBigInt t[] = new FDBigInt[ p+1 ]; 142 System.arraycopy( b5p, 0, t, 0, b5p.length ); 143 b5p = t; 144 } 145 if ( b5p[p] != null ) 146 return b5p[p]; 147 else if ( p < small5pow.length ) 148 return b5p[p] = new FDBigInt( small5pow[p] ); 149 else if ( p < long5pow.length ) 150 return b5p[p] = new FDBigInt( long5pow[p] ); 151 else { 152 // construct the value. 153 // recursively. 154 int q, r; 155 // in order to compute 5^p, 156 // compute its square root, 5^(p/2) and square. 157 // or, let q = p / 2, r = p -q, then 158 // 5^p = 5^(q+r) = 5^q * 5^r 159 q = p >> 1; 160 r = p - q; 161 FDBigInt bigq = b5p[q]; 162 if ( bigq == null ) 163 bigq = big5pow ( q ); 164 if ( r < small5pow.length ){ 165 return (b5p[p] = bigq.mult( small5pow[r] ) ); 166 }else{ 167 FDBigInt bigr = b5p[ r ]; 168 if ( bigr == null ) 169 bigr = big5pow( r ); 170 return (b5p[p] = bigq.mult( bigr ) ); 171 } 172 } 173 } 174 175 // 176 // a common operation 177 // 178 private static FDBigInt multPow52( FDBigInt v, int p5, int p2 )179 multPow52( FDBigInt v, int p5, int p2 ){ 180 if ( p5 != 0 ){ 181 if ( p5 < small5pow.length ){ 182 v = v.mult( small5pow[p5] ); 183 } else { 184 v = v.mult( big5pow( p5 ) ); 185 } 186 } 187 if ( p2 != 0 ){ 188 v.lshiftMe( p2 ); 189 } 190 return v; 191 } 192 193 // 194 // another common operation 195 // 196 private static FDBigInt constructPow52( int p5, int p2 )197 constructPow52( int p5, int p2 ){ 198 FDBigInt v = new FDBigInt( big5pow( p5 ) ); 199 if ( p2 != 0 ){ 200 v.lshiftMe( p2 ); 201 } 202 return v; 203 } 204 205 /* 206 * Make a floating double into a FDBigInt. 207 * This could also be structured as a FDBigInt 208 * constructor, but we'd have to build a lot of knowledge 209 * about floating-point representation into it, and we don't want to. 210 * 211 * AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES 212 * bigIntExp and bigIntNBits 213 * 214 */ 215 private FDBigInt doubleToBigInt( double dval )216 doubleToBigInt( double dval ){ 217 long lbits = Double.doubleToLongBits( dval ) & ~signMask; 218 int binexp = (int)(lbits >>> expShift); 219 lbits &= fractMask; 220 if ( binexp > 0 ){ 221 lbits |= fractHOB; 222 } else { 223 assert lbits != 0L : lbits; // doubleToBigInt(0.0) 224 binexp +=1; 225 while ( (lbits & fractHOB ) == 0L){ 226 lbits <<= 1; 227 binexp -= 1; 228 } 229 } 230 binexp -= expBias; 231 int nbits = countBits( lbits ); 232 /* 233 * We now know where the high-order 1 bit is, 234 * and we know how many there are. 235 */ 236 int lowOrderZeros = expShift+1-nbits; 237 lbits >>>= lowOrderZeros; 238 239 bigIntExp = binexp+1-nbits; 240 bigIntNBits = nbits; 241 return new FDBigInt( lbits ); 242 } 243 244 /* 245 * Compute a number that is the ULP of the given value, 246 * for purposes of addition/subtraction. Generally easy. 247 * More difficult if subtracting and the argument 248 * is a normalized a power of 2, as the ULP changes at these points. 249 */ ulp( double dval, boolean subtracting )250 private static double ulp( double dval, boolean subtracting ){ 251 long lbits = Double.doubleToLongBits( dval ) & ~signMask; 252 int binexp = (int)(lbits >>> expShift); 253 double ulpval; 254 if ( subtracting && ( binexp >= expShift ) && ((lbits&fractMask) == 0L) ){ 255 // for subtraction from normalized, powers of 2, 256 // use next-smaller exponent 257 binexp -= 1; 258 } 259 if ( binexp > expShift ){ 260 ulpval = Double.longBitsToDouble( ((long)(binexp-expShift))<<expShift ); 261 } else if ( binexp == 0 ){ 262 ulpval = Double.MIN_VALUE; 263 } else { 264 ulpval = Double.longBitsToDouble( 1L<<(binexp-1) ); 265 } 266 if ( subtracting ) ulpval = - ulpval; 267 268 return ulpval; 269 } 270 271 /* 272 * Round a double to a float. 273 * In addition to the fraction bits of the double, 274 * look at the class instance variable roundDir, 275 * which should help us avoid double-rounding error. 276 * roundDir was set in hardValueOf if the estimate was 277 * close enough, but not exact. It tells us which direction 278 * of rounding is preferred. 279 */ 280 float stickyRound( double dval )281 stickyRound( double dval ){ 282 long lbits = Double.doubleToLongBits( dval ); 283 long binexp = lbits & expMask; 284 if ( binexp == 0L || binexp == expMask ){ 285 // what we have here is special. 286 // don't worry, the right thing will happen. 287 return (float) dval; 288 } 289 lbits += (long)roundDir; // hack-o-matic. 290 return (float)Double.longBitsToDouble( lbits ); 291 } 292 293 294 /* 295 * This is the easy subcase -- 296 * all the significant bits, after scaling, are held in lvalue. 297 * negSign and decExponent tell us what processing and scaling 298 * has already been done. Exceptional cases have already been 299 * stripped out. 300 * In particular: 301 * lvalue is a finite number (not Inf, nor NaN) 302 * lvalue > 0L (not zero, nor negative). 303 * 304 * The only reason that we develop the digits here, rather than 305 * calling on Long.toString() is that we can do it a little faster, 306 * and besides want to treat trailing 0s specially. If Long.toString 307 * changes, we should re-evaluate this strategy! 308 */ 309 private void developLongDigits( int decExponent, long lvalue, long insignificant )310 developLongDigits( int decExponent, long lvalue, long insignificant ){ 311 char digits[]; 312 int ndigits; 313 int digitno; 314 int c; 315 // 316 // Discard non-significant low-order bits, while rounding, 317 // up to insignificant value. 318 int i; 319 for ( i = 0; insignificant >= 10L; i++ ) 320 insignificant /= 10L; 321 if ( i != 0 ){ 322 long pow10 = long5pow[i] << i; // 10^i == 5^i * 2^i; 323 long residue = lvalue % pow10; 324 lvalue /= pow10; 325 decExponent += i; 326 if ( residue >= (pow10>>1) ){ 327 // round up based on the low-order bits we're discarding 328 lvalue++; 329 } 330 } 331 if ( lvalue <= Integer.MAX_VALUE ){ 332 assert lvalue > 0L : lvalue; // lvalue <= 0 333 // even easier subcase! 334 // can do int arithmetic rather than long! 335 int ivalue = (int)lvalue; 336 ndigits = 10; 337 digits = (char[])(perThreadBuffer.get()); 338 digitno = ndigits-1; 339 c = ivalue%10; 340 ivalue /= 10; 341 while ( c == 0 ){ 342 decExponent++; 343 c = ivalue%10; 344 ivalue /= 10; 345 } 346 while ( ivalue != 0){ 347 digits[digitno--] = (char)(c+'0'); 348 decExponent++; 349 c = ivalue%10; 350 ivalue /= 10; 351 } 352 digits[digitno] = (char)(c+'0'); 353 } else { 354 // same algorithm as above (same bugs, too ) 355 // but using long arithmetic. 356 ndigits = 20; 357 digits = (char[])(perThreadBuffer.get()); 358 digitno = ndigits-1; 359 c = (int)(lvalue%10L); 360 lvalue /= 10L; 361 while ( c == 0 ){ 362 decExponent++; 363 c = (int)(lvalue%10L); 364 lvalue /= 10L; 365 } 366 while ( lvalue != 0L ){ 367 digits[digitno--] = (char)(c+'0'); 368 decExponent++; 369 c = (int)(lvalue%10L); 370 lvalue /= 10; 371 } 372 digits[digitno] = (char)(c+'0'); 373 } 374 char result []; 375 ndigits -= digitno; 376 result = new char[ ndigits ]; 377 System.arraycopy( digits, digitno, result, 0, ndigits ); 378 this.digits = result; 379 this.decExponent = decExponent+1; 380 this.nDigits = ndigits; 381 } 382 383 // 384 // add one to the least significant digit. 385 // in the unlikely event there is a carry out, 386 // deal with it. 387 // assert that this will only happen where there 388 // is only one digit, e.g. (float)1e-44 seems to do it. 389 // 390 private void roundup()391 roundup(){ 392 int i; 393 int q = digits[ i = (nDigits-1)]; 394 if ( q == '9' ){ 395 while ( q == '9' && i > 0 ){ 396 digits[i] = '0'; 397 q = digits[--i]; 398 } 399 if ( q == '9' ){ 400 // carryout! High-order 1, rest 0s, larger exp. 401 decExponent += 1; 402 digits[0] = '1'; 403 return; 404 } 405 // else fall through. 406 } 407 digits[i] = (char)(q+1); 408 } 409 410 // Given the desired number of digits predict the result's exponent. checkExponent(int length)411 private int checkExponent(int length) { 412 if (length >= nDigits || length < 0) 413 return decExponent; 414 415 for (int i = 0; i < length; i++) 416 if (digits[i] != '9') 417 // a '9' anywhere in digits will absorb the round 418 return decExponent; 419 return decExponent + (digits[length] >= '5' ? 1 : 0); 420 } 421 422 // Unlike roundup(), this method does not modify digits. It also 423 // rounds at a particular precision. applyPrecision(int length)424 private char [] applyPrecision(int length) { 425 char [] result = new char[nDigits]; 426 for (int i = 0; i < result.length; i++) result[i] = '0'; 427 428 if (length >= nDigits || length < 0) { 429 // no rounding necessary 430 System.arraycopy(digits, 0, result, 0, nDigits); 431 return result; 432 } 433 if (length == 0) { 434 // only one digit (0 or 1) is returned because the precision 435 // excludes all significant digits 436 if (digits[0] >= '5') { 437 result[0] = '1'; 438 } 439 return result; 440 } 441 442 int i = length; 443 int q = digits[i]; 444 if (q >= '5' && i > 0) { 445 q = digits[--i]; 446 if ( q == '9' ) { 447 while ( q == '9' && i > 0 ){ 448 q = digits[--i]; 449 } 450 if ( q == '9' ){ 451 // carryout! High-order 1, rest 0s, larger exp. 452 result[0] = '1'; 453 return result; 454 } 455 } 456 result[i] = (char)(q + 1); 457 } 458 while (--i >= 0) { 459 result[i] = digits[i]; 460 } 461 return result; 462 } 463 464 /* 465 * FIRST IMPORTANT CONSTRUCTOR: DOUBLE 466 */ FormattedFloatingDecimal( double d )467 public FormattedFloatingDecimal( double d ) 468 { 469 this(d, Integer.MAX_VALUE, Form.COMPATIBLE); 470 } 471 FormattedFloatingDecimal( double d, int precision, Form form )472 public FormattedFloatingDecimal( double d, int precision, Form form ) 473 { 474 long dBits = Double.doubleToLongBits( d ); 475 long fractBits; 476 int binExp; 477 int nSignificantBits; 478 479 this.precision = precision; 480 this.form = form; 481 482 // discover and delete sign 483 if ( (dBits&signMask) != 0 ){ 484 isNegative = true; 485 dBits ^= signMask; 486 } else { 487 isNegative = false; 488 } 489 // Begin to unpack 490 // Discover obvious special cases of NaN and Infinity. 491 binExp = (int)( (dBits&expMask) >> expShift ); 492 fractBits = dBits&fractMask; 493 if ( binExp == (int)(expMask>>expShift) ) { 494 isExceptional = true; 495 if ( fractBits == 0L ){ 496 digits = infinity; 497 } else { 498 digits = notANumber; 499 isNegative = false; // NaN has no sign! 500 } 501 nDigits = digits.length; 502 return; 503 } 504 isExceptional = false; 505 // Finish unpacking 506 // Normalize denormalized numbers. 507 // Insert assumed high-order bit for normalized numbers. 508 // Subtract exponent bias. 509 if ( binExp == 0 ){ 510 if ( fractBits == 0L ){ 511 // not a denorm, just a 0! 512 decExponent = 0; 513 digits = zero; 514 nDigits = 1; 515 return; 516 } 517 while ( (fractBits&fractHOB) == 0L ){ 518 fractBits <<= 1; 519 binExp -= 1; 520 } 521 nSignificantBits = expShift + binExp +1; // recall binExp is - shift count. 522 binExp += 1; 523 } else { 524 fractBits |= fractHOB; 525 nSignificantBits = expShift+1; 526 } 527 binExp -= expBias; 528 // call the routine that actually does all the hard work. 529 dtoa( binExp, fractBits, nSignificantBits ); 530 } 531 532 /* 533 * SECOND IMPORTANT CONSTRUCTOR: SINGLE 534 */ FormattedFloatingDecimal( float f )535 public FormattedFloatingDecimal( float f ) 536 { 537 this(f, Integer.MAX_VALUE, Form.COMPATIBLE); 538 } FormattedFloatingDecimal( float f, int precision, Form form )539 public FormattedFloatingDecimal( float f, int precision, Form form ) 540 { 541 int fBits = Float.floatToIntBits( f ); 542 int fractBits; 543 int binExp; 544 int nSignificantBits; 545 546 this.precision = precision; 547 this.form = form; 548 549 // discover and delete sign 550 if ( (fBits&singleSignMask) != 0 ){ 551 isNegative = true; 552 fBits ^= singleSignMask; 553 } else { 554 isNegative = false; 555 } 556 // Begin to unpack 557 // Discover obvious special cases of NaN and Infinity. 558 binExp = (int)( (fBits&singleExpMask) >> singleExpShift ); 559 fractBits = fBits&singleFractMask; 560 if ( binExp == (int)(singleExpMask>>singleExpShift) ) { 561 isExceptional = true; 562 if ( fractBits == 0L ){ 563 digits = infinity; 564 } else { 565 digits = notANumber; 566 isNegative = false; // NaN has no sign! 567 } 568 nDigits = digits.length; 569 return; 570 } 571 isExceptional = false; 572 // Finish unpacking 573 // Normalize denormalized numbers. 574 // Insert assumed high-order bit for normalized numbers. 575 // Subtract exponent bias. 576 if ( binExp == 0 ){ 577 if ( fractBits == 0 ){ 578 // not a denorm, just a 0! 579 decExponent = 0; 580 digits = zero; 581 nDigits = 1; 582 return; 583 } 584 while ( (fractBits&singleFractHOB) == 0 ){ 585 fractBits <<= 1; 586 binExp -= 1; 587 } 588 nSignificantBits = singleExpShift + binExp +1; // recall binExp is - shift count. 589 binExp += 1; 590 } else { 591 fractBits |= singleFractHOB; 592 nSignificantBits = singleExpShift+1; 593 } 594 binExp -= singleExpBias; 595 // call the routine that actually does all the hard work. 596 dtoa( binExp, ((long)fractBits)<<(expShift-singleExpShift), nSignificantBits ); 597 } 598 599 private void dtoa( int binExp, long fractBits, int nSignificantBits )600 dtoa( int binExp, long fractBits, int nSignificantBits ) 601 { 602 int nFractBits; // number of significant bits of fractBits; 603 int nTinyBits; // number of these to the right of the point. 604 int decExp; 605 606 // Examine number. Determine if it is an easy case, 607 // which we can do pretty trivially using float/long conversion, 608 // or whether we must do real work. 609 nFractBits = countBits( fractBits ); 610 nTinyBits = Math.max( 0, nFractBits - binExp - 1 ); 611 if ( binExp <= maxSmallBinExp && binExp >= minSmallBinExp ){ 612 // Look more closely at the number to decide if, 613 // with scaling by 10^nTinyBits, the result will fit in 614 // a long. 615 if ( (nTinyBits < long5pow.length) && ((nFractBits + n5bits[nTinyBits]) < 64 ) ){ 616 /* 617 * We can do this: 618 * take the fraction bits, which are normalized. 619 * (a) nTinyBits == 0: Shift left or right appropriately 620 * to align the binary point at the extreme right, i.e. 621 * where a long int point is expected to be. The integer 622 * result is easily converted to a string. 623 * (b) nTinyBits > 0: Shift right by expShift-nFractBits, 624 * which effectively converts to long and scales by 625 * 2^nTinyBits. Then multiply by 5^nTinyBits to 626 * complete the scaling. We know this won't overflow 627 * because we just counted the number of bits necessary 628 * in the result. The integer you get from this can 629 * then be converted to a string pretty easily. 630 */ 631 long halfULP; 632 if ( nTinyBits == 0 ) { 633 if ( binExp > nSignificantBits ){ 634 halfULP = 1L << ( binExp-nSignificantBits-1); 635 } else { 636 halfULP = 0L; 637 } 638 if ( binExp >= expShift ){ 639 fractBits <<= (binExp-expShift); 640 } else { 641 fractBits >>>= (expShift-binExp) ; 642 } 643 developLongDigits( 0, fractBits, halfULP ); 644 return; 645 } 646 /* 647 * The following causes excess digits to be printed 648 * out in the single-float case. Our manipulation of 649 * halfULP here is apparently not correct. If we 650 * better understand how this works, perhaps we can 651 * use this special case again. But for the time being, 652 * we do not. 653 * else { 654 * fractBits >>>= expShift+1-nFractBits; 655 * fractBits *= long5pow[ nTinyBits ]; 656 * halfULP = long5pow[ nTinyBits ] >> (1+nSignificantBits-nFractBits); 657 * developLongDigits( -nTinyBits, fractBits, halfULP ); 658 * return; 659 * } 660 */ 661 } 662 } 663 /* 664 * This is the hard case. We are going to compute large positive 665 * integers B and S and integer decExp, s.t. 666 * d = ( B / S ) * 10^decExp 667 * 1 <= B / S < 10 668 * Obvious choices are: 669 * decExp = floor( log10(d) ) 670 * B = d * 2^nTinyBits * 10^max( 0, -decExp ) 671 * S = 10^max( 0, decExp) * 2^nTinyBits 672 * (noting that nTinyBits has already been forced to non-negative) 673 * I am also going to compute a large positive integer 674 * M = (1/2^nSignificantBits) * 2^nTinyBits * 10^max( 0, -decExp ) 675 * i.e. M is (1/2) of the ULP of d, scaled like B. 676 * When we iterate through dividing B/S and picking off the 677 * quotient bits, we will know when to stop when the remainder 678 * is <= M. 679 * 680 * We keep track of powers of 2 and powers of 5. 681 */ 682 683 /* 684 * Estimate decimal exponent. (If it is small-ish, 685 * we could double-check.) 686 * 687 * First, scale the mantissa bits such that 1 <= d2 < 2. 688 * We are then going to estimate 689 * log10(d2) ~=~ (d2-1.5)/1.5 + log(1.5) 690 * and so we can estimate 691 * log10(d) ~=~ log10(d2) + binExp * log10(2) 692 * take the floor and call it decExp. 693 * FIXME -- use more precise constants here. It costs no more. 694 */ 695 double d2 = Double.longBitsToDouble( 696 expOne | ( fractBits &~ fractHOB ) ); 697 decExp = (int)Math.floor( 698 (d2-1.5D)*0.289529654D + 0.176091259 + (double)binExp * 0.301029995663981 ); 699 int B2, B5; // powers of 2 and powers of 5, respectively, in B 700 int S2, S5; // powers of 2 and powers of 5, respectively, in S 701 int M2, M5; // powers of 2 and powers of 5, respectively, in M 702 int Bbits; // binary digits needed to represent B, approx. 703 int tenSbits; // binary digits needed to represent 10*S, approx. 704 FDBigInt Sval, Bval, Mval; 705 706 B5 = Math.max( 0, -decExp ); 707 B2 = B5 + nTinyBits + binExp; 708 709 S5 = Math.max( 0, decExp ); 710 S2 = S5 + nTinyBits; 711 712 M5 = B5; 713 M2 = B2 - nSignificantBits; 714 715 /* 716 * the long integer fractBits contains the (nFractBits) interesting 717 * bits from the mantissa of d ( hidden 1 added if necessary) followed 718 * by (expShift+1-nFractBits) zeros. In the interest of compactness, 719 * I will shift out those zeros before turning fractBits into a 720 * FDBigInt. The resulting whole number will be 721 * d * 2^(nFractBits-1-binExp). 722 */ 723 fractBits >>>= (expShift+1-nFractBits); 724 B2 -= nFractBits-1; 725 int common2factor = Math.min( B2, S2 ); 726 B2 -= common2factor; 727 S2 -= common2factor; 728 M2 -= common2factor; 729 730 /* 731 * HACK!! For exact powers of two, the next smallest number 732 * is only half as far away as we think (because the meaning of 733 * ULP changes at power-of-two bounds) for this reason, we 734 * hack M2. Hope this works. 735 */ 736 if ( nFractBits == 1 ) 737 M2 -= 1; 738 739 if ( M2 < 0 ){ 740 // oops. 741 // since we cannot scale M down far enough, 742 // we must scale the other values up. 743 B2 -= M2; 744 S2 -= M2; 745 M2 = 0; 746 } 747 /* 748 * Construct, Scale, iterate. 749 * Some day, we'll write a stopping test that takes 750 * account of the assymetry of the spacing of floating-point 751 * numbers below perfect powers of 2 752 * 26 Sept 96 is not that day. 753 * So we use a symmetric test. 754 */ 755 char digits[] = this.digits = new char[18]; 756 int ndigit = 0; 757 boolean low, high; 758 long lowDigitDifference; 759 int q; 760 761 /* 762 * Detect the special cases where all the numbers we are about 763 * to compute will fit in int or long integers. 764 * In these cases, we will avoid doing FDBigInt arithmetic. 765 * We use the same algorithms, except that we "normalize" 766 * our FDBigInts before iterating. This is to make division easier, 767 * as it makes our fist guess (quotient of high-order words) 768 * more accurate! 769 * 770 * Some day, we'll write a stopping test that takes 771 * account of the assymetry of the spacing of floating-point 772 * numbers below perfect powers of 2 773 * 26 Sept 96 is not that day. 774 * So we use a symmetric test. 775 */ 776 Bbits = nFractBits + B2 + (( B5 < n5bits.length )? n5bits[B5] : ( B5*3 )); 777 tenSbits = S2+1 + (( (S5+1) < n5bits.length )? n5bits[(S5+1)] : ( (S5+1)*3 )); 778 if ( Bbits < 64 && tenSbits < 64){ 779 if ( Bbits < 32 && tenSbits < 32){ 780 // wa-hoo! They're all ints! 781 int b = ((int)fractBits * small5pow[B5] ) << B2; 782 int s = small5pow[S5] << S2; 783 int m = small5pow[M5] << M2; 784 int tens = s * 10; 785 /* 786 * Unroll the first iteration. If our decExp estimate 787 * was too high, our first quotient will be zero. In this 788 * case, we discard it and decrement decExp. 789 */ 790 ndigit = 0; 791 q = b / s; 792 b = 10 * ( b % s ); 793 m *= 10; 794 low = (b < m ); 795 high = (b+m > tens ); 796 assert q < 10 : q; // excessively large digit 797 if ( (q == 0) && ! high ){ 798 // oops. Usually ignore leading zero. 799 decExp--; 800 } else { 801 digits[ndigit++] = (char)('0' + q); 802 } 803 /* 804 * HACK! Java spec sez that we always have at least 805 * one digit after the . in either F- or E-form output. 806 * Thus we will need more than one digit if we're using 807 * E-form 808 */ 809 if (! (form == Form.COMPATIBLE && -3 < decExp && decExp < 8)) { 810 high = low = false; 811 } 812 while( ! low && ! high ){ 813 q = b / s; 814 b = 10 * ( b % s ); 815 m *= 10; 816 assert q < 10 : q; // excessively large digit 817 if ( m > 0L ){ 818 low = (b < m ); 819 high = (b+m > tens ); 820 } else { 821 // hack -- m might overflow! 822 // in this case, it is certainly > b, 823 // which won't 824 // and b+m > tens, too, since that has overflowed 825 // either! 826 low = true; 827 high = true; 828 } 829 digits[ndigit++] = (char)('0' + q); 830 } 831 lowDigitDifference = (b<<1) - tens; 832 } else { 833 // still good! they're all longs! 834 long b = (fractBits * long5pow[B5] ) << B2; 835 long s = long5pow[S5] << S2; 836 long m = long5pow[M5] << M2; 837 long tens = s * 10L; 838 /* 839 * Unroll the first iteration. If our decExp estimate 840 * was too high, our first quotient will be zero. In this 841 * case, we discard it and decrement decExp. 842 */ 843 ndigit = 0; 844 q = (int) ( b / s ); 845 b = 10L * ( b % s ); 846 m *= 10L; 847 low = (b < m ); 848 high = (b+m > tens ); 849 assert q < 10 : q; // excessively large digit 850 if ( (q == 0) && ! high ){ 851 // oops. Usually ignore leading zero. 852 decExp--; 853 } else { 854 digits[ndigit++] = (char)('0' + q); 855 } 856 /* 857 * HACK! Java spec sez that we always have at least 858 * one digit after the . in either F- or E-form output. 859 * Thus we will need more than one digit if we're using 860 * E-form 861 */ 862 if (! (form == Form.COMPATIBLE && -3 < decExp && decExp < 8)) { 863 high = low = false; 864 } 865 while( ! low && ! high ){ 866 q = (int) ( b / s ); 867 b = 10 * ( b % s ); 868 m *= 10; 869 assert q < 10 : q; // excessively large digit 870 if ( m > 0L ){ 871 low = (b < m ); 872 high = (b+m > tens ); 873 } else { 874 // hack -- m might overflow! 875 // in this case, it is certainly > b, 876 // which won't 877 // and b+m > tens, too, since that has overflowed 878 // either! 879 low = true; 880 high = true; 881 } 882 digits[ndigit++] = (char)('0' + q); 883 } 884 lowDigitDifference = (b<<1) - tens; 885 } 886 } else { 887 FDBigInt tenSval; 888 int shiftBias; 889 890 /* 891 * We really must do FDBigInt arithmetic. 892 * Fist, construct our FDBigInt initial values. 893 */ 894 Bval = multPow52( new FDBigInt( fractBits ), B5, B2 ); 895 Sval = constructPow52( S5, S2 ); 896 Mval = constructPow52( M5, M2 ); 897 898 899 // normalize so that division works better 900 Bval.lshiftMe( shiftBias = Sval.normalizeMe() ); 901 Mval.lshiftMe( shiftBias ); 902 tenSval = Sval.mult( 10 ); 903 /* 904 * Unroll the first iteration. If our decExp estimate 905 * was too high, our first quotient will be zero. In this 906 * case, we discard it and decrement decExp. 907 */ 908 ndigit = 0; 909 q = Bval.quoRemIteration( Sval ); 910 Mval = Mval.mult( 10 ); 911 low = (Bval.cmp( Mval ) < 0); 912 high = (Bval.add( Mval ).cmp( tenSval ) > 0 ); 913 assert q < 10 : q; // excessively large digit 914 if ( (q == 0) && ! high ){ 915 // oops. Usually ignore leading zero. 916 decExp--; 917 } else { 918 digits[ndigit++] = (char)('0' + q); 919 } 920 /* 921 * HACK! Java spec sez that we always have at least 922 * one digit after the . in either F- or E-form output. 923 * Thus we will need more than one digit if we're using 924 * E-form 925 */ 926 if (! (form == Form.COMPATIBLE && -3 < decExp && decExp < 8)) { 927 high = low = false; 928 } 929 while( ! low && ! high ){ 930 q = Bval.quoRemIteration( Sval ); 931 Mval = Mval.mult( 10 ); 932 assert q < 10 : q; // excessively large digit 933 low = (Bval.cmp( Mval ) < 0); 934 high = (Bval.add( Mval ).cmp( tenSval ) > 0 ); 935 digits[ndigit++] = (char)('0' + q); 936 } 937 if ( high && low ){ 938 Bval.lshiftMe(1); 939 lowDigitDifference = Bval.cmp(tenSval); 940 } else 941 lowDigitDifference = 0L; // this here only for flow analysis! 942 } 943 this.decExponent = decExp+1; 944 this.digits = digits; 945 this.nDigits = ndigit; 946 /* 947 * Last digit gets rounded based on stopping condition. 948 */ 949 if ( high ){ 950 if ( low ){ 951 if ( lowDigitDifference == 0L ){ 952 // it's a tie! 953 // choose based on which digits we like. 954 if ( (digits[nDigits-1]&1) != 0 ) roundup(); 955 } else if ( lowDigitDifference > 0 ){ 956 roundup(); 957 } 958 } else { 959 roundup(); 960 } 961 } 962 } 963 964 public String 965 toString(){ 966 // most brain-dead version 967 StringBuffer result = new StringBuffer( nDigits+8 ); 968 if ( isNegative ){ result.append( '-' ); } 969 if ( isExceptional ){ 970 result.append( digits, 0, nDigits ); 971 } else { 972 result.append( "0."); 973 result.append( digits, 0, nDigits ); 974 result.append('e'); 975 result.append( decExponent ); 976 } 977 return new String(result); 978 } 979 980 // returns the exponent before rounding 981 public int getExponent() { 982 return decExponent - 1; 983 } 984 985 // returns the exponent after rounding has been done by applyPrecision 986 public int getExponentRounded() { 987 return decExponentRounded - 1; 988 } 989 990 public int getChars(char[] result) { 991 assert nDigits <= 19 : nDigits; // generous bound on size of nDigits 992 int i = 0; 993 if (isNegative) { result[0] = '-'; i = 1; } 994 if (isExceptional) { 995 System.arraycopy(digits, 0, result, i, nDigits); 996 i += nDigits; 997 } else { 998 char digits [] = this.digits; 999 int exp = decExponent; 1000 switch (form) { 1001 case COMPATIBLE: 1002 break; 1003 case DECIMAL_FLOAT: 1004 exp = checkExponent(decExponent + precision); 1005 digits = applyPrecision(decExponent + precision); 1006 break; 1007 case SCIENTIFIC: 1008 exp = checkExponent(precision + 1); 1009 digits = applyPrecision(precision + 1); 1010 break; 1011 case GENERAL: 1012 exp = checkExponent(precision); 1013 digits = applyPrecision(precision); 1014 // adjust precision to be the number of digits to right of decimal 1015 // the real exponent to be output is actually exp - 1, not exp 1016 if (exp - 1 < -4 || exp - 1 >= precision) { 1017 form = Form.SCIENTIFIC; 1018 precision--; 1019 } else { 1020 form = Form.DECIMAL_FLOAT; 1021 precision = precision - exp; 1022 } 1023 break; 1024 default: 1025 assert false; 1026 } 1027 decExponentRounded = exp; 1028 1029 if (exp > 0 1030 && ((form == Form.COMPATIBLE && (exp < 8)) 1031 || (form == Form.DECIMAL_FLOAT))) 1032 { 1033 // print digits.digits. 1034 int charLength = Math.min(nDigits, exp); 1035 System.arraycopy(digits, 0, result, i, charLength); 1036 i += charLength; 1037 if (charLength < exp) { 1038 charLength = exp-charLength; 1039 for (int nz = 0; nz < charLength; nz++) 1040 result[i++] = '0'; 1041 // Do not append ".0" for formatted floats since the user 1042 // may request that it be omitted. It is added as necessary 1043 // by the Formatter. 1044 if (form == Form.COMPATIBLE) { 1045 result[i++] = '.'; 1046 result[i++] = '0'; 1047 } 1048 } else { 1049 // Do not append ".0" for formatted floats since the user 1050 // may request that it be omitted. It is added as necessary 1051 // by the Formatter. 1052 if (form == Form.COMPATIBLE) { 1053 result[i++] = '.'; 1054 if (charLength < nDigits) { 1055 int t = Math.min(nDigits - charLength, precision); 1056 System.arraycopy(digits, charLength, result, i, t); 1057 i += t; 1058 } else { 1059 result[i++] = '0'; 1060 } 1061 } else { 1062 int t = Math.min(nDigits - charLength, precision); 1063 if (t > 0) { 1064 result[i++] = '.'; 1065 System.arraycopy(digits, charLength, result, i, t); 1066 i += t; 1067 } 1068 } 1069 } 1070 } else if (exp <= 0 1071 && ((form == Form.COMPATIBLE && exp > -3) 1072 || (form == Form.DECIMAL_FLOAT))) 1073 { 1074 // print 0.0* digits 1075 result[i++] = '0'; 1076 if (exp != 0) { 1077 // write '0' s before the significant digits 1078 int t = Math.min(-exp, precision); 1079 if (t > 0) { 1080 result[i++] = '.'; 1081 for (int nz = 0; nz < t; nz++) 1082 result[i++] = '0'; 1083 } 1084 } 1085 int t = Math.min(digits.length, precision + exp); 1086 if (t > 0) { 1087 if (i == 1) 1088 result[i++] = '.'; 1089 // copy only when significant digits are within the precision 1090 System.arraycopy(digits, 0, result, i, t); 1091 i += t; 1092 } 1093 } else { 1094 result[i++] = digits[0]; 1095 if (form == Form.COMPATIBLE) { 1096 result[i++] = '.'; 1097 if (nDigits > 1) { 1098 System.arraycopy(digits, 1, result, i, nDigits-1); 1099 i += nDigits-1; 1100 } else { 1101 result[i++] = '0'; 1102 } 1103 result[i++] = 'E'; 1104 } else { 1105 if (nDigits > 1) { 1106 int t = Math.min(nDigits -1, precision); 1107 if (t > 0) { 1108 result[i++] = '.'; 1109 System.arraycopy(digits, 1, result, i, t); 1110 i += t; 1111 } 1112 } 1113 result[i++] = 'e'; 1114 } 1115 int e; 1116 if (exp <= 0) { 1117 result[i++] = '-'; 1118 e = -exp+1; 1119 } else { 1120 if (form != Form.COMPATIBLE) 1121 result[i++] = '+'; 1122 e = exp-1; 1123 } 1124 // decExponent has 1, 2, or 3, digits 1125 if (e <= 9) { 1126 if (form != Form.COMPATIBLE) 1127 result[i++] = '0'; 1128 result[i++] = (char)(e+'0'); 1129 } else if (e <= 99) { 1130 result[i++] = (char)(e/10 +'0'); 1131 result[i++] = (char)(e%10 + '0'); 1132 } else { 1133 result[i++] = (char)(e/100+'0'); 1134 e %= 100; 1135 result[i++] = (char)(e/10+'0'); 1136 result[i++] = (char)(e%10 + '0'); 1137 } 1138 } 1139 } 1140 return i; 1141 } 1142 1143 // Per-thread buffer for string/stringbuffer conversion 1144 private static ThreadLocal perThreadBuffer = new ThreadLocal() { 1145 protected synchronized Object initialValue() { 1146 return new char[26]; 1147 } 1148 }; 1149 1150 /* 1151 * Take a FormattedFloatingDecimal, which we presumably just scanned in, 1152 * and find out what its value is, as a double. 1153 * 1154 * AS A SIDE EFFECT, SET roundDir TO INDICATE PREFERRED 1155 * ROUNDING DIRECTION in case the result is really destined 1156 * for a single-precision float. 1157 */ 1158 1159 public strictfp double doubleValue(){ 1160 int kDigits = Math.min( nDigits, maxDecimalDigits+1 ); 1161 long lValue; 1162 double dValue; 1163 double rValue, tValue; 1164 1165 // First, check for NaN and Infinity values 1166 if(digits == infinity || digits == notANumber) { 1167 if(digits == notANumber) 1168 return Double.NaN; 1169 else 1170 return (isNegative?Double.NEGATIVE_INFINITY:Double.POSITIVE_INFINITY); 1171 } 1172 else { 1173 if (mustSetRoundDir) { 1174 roundDir = 0; 1175 } 1176 /* 1177 * convert the lead kDigits to a long integer. 1178 */ 1179 // (special performance hack: start to do it using int) 1180 int iValue = (int)digits[0]-(int)'0'; 1181 int iDigits = Math.min( kDigits, intDecimalDigits ); 1182 for ( int i=1; i < iDigits; i++ ){ 1183 iValue = iValue*10 + (int)digits[i]-(int)'0'; 1184 } 1185 lValue = (long)iValue; 1186 for ( int i=iDigits; i < kDigits; i++ ){ 1187 lValue = lValue*10L + (long)((int)digits[i]-(int)'0'); 1188 } 1189 dValue = (double)lValue; 1190 int exp = decExponent-kDigits; 1191 /* 1192 * lValue now contains a long integer with the value of 1193 * the first kDigits digits of the number. 1194 * dValue contains the (double) of the same. 1195 */ 1196 1197 if ( nDigits <= maxDecimalDigits ){ 1198 /* 1199 * possibly an easy case. 1200 * We know that the digits can be represented 1201 * exactly. And if the exponent isn't too outrageous, 1202 * the whole thing can be done with one operation, 1203 * thus one rounding error. 1204 * Note that all our constructors trim all leading and 1205 * trailing zeros, so simple values (including zero) 1206 * will always end up here 1207 */ 1208 if (exp == 0 || dValue == 0.0) 1209 return (isNegative)? -dValue : dValue; // small floating integer 1210 else if ( exp >= 0 ){ 1211 if ( exp <= maxSmallTen ){ 1212 /* 1213 * Can get the answer with one operation, 1214 * thus one roundoff. 1215 */ 1216 rValue = dValue * small10pow[exp]; 1217 if ( mustSetRoundDir ){ 1218 tValue = rValue / small10pow[exp]; 1219 roundDir = ( tValue == dValue ) ? 0 1220 :( tValue < dValue ) ? 1 1221 : -1; 1222 } 1223 return (isNegative)? -rValue : rValue; 1224 } 1225 int slop = maxDecimalDigits - kDigits; 1226 if ( exp <= maxSmallTen+slop ){ 1227 /* 1228 * We can multiply dValue by 10^(slop) 1229 * and it is still "small" and exact. 1230 * Then we can multiply by 10^(exp-slop) 1231 * with one rounding. 1232 */ 1233 dValue *= small10pow[slop]; 1234 rValue = dValue * small10pow[exp-slop]; 1235 1236 if ( mustSetRoundDir ){ 1237 tValue = rValue / small10pow[exp-slop]; 1238 roundDir = ( tValue == dValue ) ? 0 1239 :( tValue < dValue ) ? 1 1240 : -1; 1241 } 1242 return (isNegative)? -rValue : rValue; 1243 } 1244 /* 1245 * Else we have a hard case with a positive exp. 1246 */ 1247 } else { 1248 if ( exp >= -maxSmallTen ){ 1249 /* 1250 * Can get the answer in one division. 1251 */ 1252 rValue = dValue / small10pow[-exp]; 1253 tValue = rValue * small10pow[-exp]; 1254 if ( mustSetRoundDir ){ 1255 roundDir = ( tValue == dValue ) ? 0 1256 :( tValue < dValue ) ? 1 1257 : -1; 1258 } 1259 return (isNegative)? -rValue : rValue; 1260 } 1261 /* 1262 * Else we have a hard case with a negative exp. 1263 */ 1264 } 1265 } 1266 1267 /* 1268 * Harder cases: 1269 * The sum of digits plus exponent is greater than 1270 * what we think we can do with one error. 1271 * 1272 * Start by approximating the right answer by, 1273 * naively, scaling by powers of 10. 1274 */ 1275 if ( exp > 0 ){ 1276 if ( decExponent > maxDecimalExponent+1 ){ 1277 /* 1278 * Lets face it. This is going to be 1279 * Infinity. Cut to the chase. 1280 */ 1281 return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY; 1282 } 1283 if ( (exp&15) != 0 ){ 1284 dValue *= small10pow[exp&15]; 1285 } 1286 if ( (exp>>=4) != 0 ){ 1287 int j; 1288 for( j = 0; exp > 1; j++, exp>>=1 ){ 1289 if ( (exp&1)!=0) 1290 dValue *= big10pow[j]; 1291 } 1292 /* 1293 * The reason for the weird exp > 1 condition 1294 * in the above loop was so that the last multiply 1295 * would get unrolled. We handle it here. 1296 * It could overflow. 1297 */ 1298 double t = dValue * big10pow[j]; 1299 if ( Double.isInfinite( t ) ){ 1300 /* 1301 * It did overflow. 1302 * Look more closely at the result. 1303 * If the exponent is just one too large, 1304 * then use the maximum finite as our estimate 1305 * value. Else call the result infinity 1306 * and punt it. 1307 * ( I presume this could happen because 1308 * rounding forces the result here to be 1309 * an ULP or two larger than 1310 * Double.MAX_VALUE ). 1311 */ 1312 t = dValue / 2.0; 1313 t *= big10pow[j]; 1314 if ( Double.isInfinite( t ) ){ 1315 return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY; 1316 } 1317 t = Double.MAX_VALUE; 1318 } 1319 dValue = t; 1320 } 1321 } else if ( exp < 0 ){ 1322 exp = -exp; 1323 if ( decExponent < minDecimalExponent-1 ){ 1324 /* 1325 * Lets face it. This is going to be 1326 * zero. Cut to the chase. 1327 */ 1328 return (isNegative)? -0.0 : 0.0; 1329 } 1330 if ( (exp&15) != 0 ){ 1331 dValue /= small10pow[exp&15]; 1332 } 1333 if ( (exp>>=4) != 0 ){ 1334 int j; 1335 for( j = 0; exp > 1; j++, exp>>=1 ){ 1336 if ( (exp&1)!=0) 1337 dValue *= tiny10pow[j]; 1338 } 1339 /* 1340 * The reason for the weird exp > 1 condition 1341 * in the above loop was so that the last multiply 1342 * would get unrolled. We handle it here. 1343 * It could underflow. 1344 */ 1345 double t = dValue * tiny10pow[j]; 1346 if ( t == 0.0 ){ 1347 /* 1348 * It did underflow. 1349 * Look more closely at the result. 1350 * If the exponent is just one too small, 1351 * then use the minimum finite as our estimate 1352 * value. Else call the result 0.0 1353 * and punt it. 1354 * ( I presume this could happen because 1355 * rounding forces the result here to be 1356 * an ULP or two less than 1357 * Double.MIN_VALUE ). 1358 */ 1359 t = dValue * 2.0; 1360 t *= tiny10pow[j]; 1361 if ( t == 0.0 ){ 1362 return (isNegative)? -0.0 : 0.0; 1363 } 1364 t = Double.MIN_VALUE; 1365 } 1366 dValue = t; 1367 } 1368 } 1369 1370 /* 1371 * dValue is now approximately the result. 1372 * The hard part is adjusting it, by comparison 1373 * with FDBigInt arithmetic. 1374 * Formulate the EXACT big-number result as 1375 * bigD0 * 10^exp 1376 */ 1377 FDBigInt bigD0 = new FDBigInt( lValue, digits, kDigits, nDigits ); 1378 exp = decExponent - nDigits; 1379 1380 correctionLoop: 1381 while(true){ 1382 /* AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES 1383 * bigIntExp and bigIntNBits 1384 */ 1385 FDBigInt bigB = doubleToBigInt( dValue ); 1386 1387 /* 1388 * Scale bigD, bigB appropriately for 1389 * big-integer operations. 1390 * Naively, we multipy by powers of ten 1391 * and powers of two. What we actually do 1392 * is keep track of the powers of 5 and 1393 * powers of 2 we would use, then factor out 1394 * common divisors before doing the work. 1395 */ 1396 int B2, B5; // powers of 2, 5 in bigB 1397 int D2, D5; // powers of 2, 5 in bigD 1398 int Ulp2; // powers of 2 in halfUlp. 1399 if ( exp >= 0 ){ 1400 B2 = B5 = 0; 1401 D2 = D5 = exp; 1402 } else { 1403 B2 = B5 = -exp; 1404 D2 = D5 = 0; 1405 } 1406 if ( bigIntExp >= 0 ){ 1407 B2 += bigIntExp; 1408 } else { 1409 D2 -= bigIntExp; 1410 } 1411 Ulp2 = B2; 1412 // shift bigB and bigD left by a number s. t. 1413 // halfUlp is still an integer. 1414 int hulpbias; 1415 if ( bigIntExp+bigIntNBits <= -expBias+1 ){ 1416 // This is going to be a denormalized number 1417 // (if not actually zero). 1418 // half an ULP is at 2^-(expBias+expShift+1) 1419 hulpbias = bigIntExp+ expBias + expShift; 1420 } else { 1421 hulpbias = expShift + 2 - bigIntNBits; 1422 } 1423 B2 += hulpbias; 1424 D2 += hulpbias; 1425 // if there are common factors of 2, we might just as well 1426 // factor them out, as they add nothing useful. 1427 int common2 = Math.min( B2, Math.min( D2, Ulp2 ) ); 1428 B2 -= common2; 1429 D2 -= common2; 1430 Ulp2 -= common2; 1431 // do multiplications by powers of 5 and 2 1432 bigB = multPow52( bigB, B5, B2 ); 1433 FDBigInt bigD = multPow52( new FDBigInt( bigD0 ), D5, D2 ); 1434 // 1435 // to recap: 1436 // bigB is the scaled-big-int version of our floating-point 1437 // candidate. 1438 // bigD is the scaled-big-int version of the exact value 1439 // as we understand it. 1440 // halfUlp is 1/2 an ulp of bigB, except for special cases 1441 // of exact powers of 2 1442 // 1443 // the plan is to compare bigB with bigD, and if the difference 1444 // is less than halfUlp, then we're satisfied. Otherwise, 1445 // use the ratio of difference to halfUlp to calculate a fudge 1446 // factor to add to the floating value, then go 'round again. 1447 // 1448 FDBigInt diff; 1449 int cmpResult; 1450 boolean overvalue; 1451 if ( (cmpResult = bigB.cmp( bigD ) ) > 0 ){ 1452 overvalue = true; // our candidate is too big. 1453 diff = bigB.sub( bigD ); 1454 if ( (bigIntNBits == 1) && (bigIntExp > -expBias) ){ 1455 // candidate is a normalized exact power of 2 and 1456 // is too big. We will be subtracting. 1457 // For our purposes, ulp is the ulp of the 1458 // next smaller range. 1459 Ulp2 -= 1; 1460 if ( Ulp2 < 0 ){ 1461 // rats. Cannot de-scale ulp this far. 1462 // must scale diff in other direction. 1463 Ulp2 = 0; 1464 diff.lshiftMe( 1 ); 1465 } 1466 } 1467 } else if ( cmpResult < 0 ){ 1468 overvalue = false; // our candidate is too small. 1469 diff = bigD.sub( bigB ); 1470 } else { 1471 // the candidate is exactly right! 1472 // this happens with surprising fequency 1473 break correctionLoop; 1474 } 1475 FDBigInt halfUlp = constructPow52( B5, Ulp2 ); 1476 if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){ 1477 // difference is small. 1478 // this is close enough 1479 if (mustSetRoundDir) { 1480 roundDir = overvalue ? -1 : 1; 1481 } 1482 break correctionLoop; 1483 } else if ( cmpResult == 0 ){ 1484 // difference is exactly half an ULP 1485 // round to some other value maybe, then finish 1486 dValue += 0.5*ulp( dValue, overvalue ); 1487 // should check for bigIntNBits == 1 here?? 1488 if (mustSetRoundDir) { 1489 roundDir = overvalue ? -1 : 1; 1490 } 1491 break correctionLoop; 1492 } else { 1493 // difference is non-trivial. 1494 // could scale addend by ratio of difference to 1495 // halfUlp here, if we bothered to compute that difference. 1496 // Most of the time ( I hope ) it is about 1 anyway. 1497 dValue += ulp( dValue, overvalue ); 1498 if ( dValue == 0.0 || dValue == Double.POSITIVE_INFINITY ) 1499 break correctionLoop; // oops. Fell off end of range. 1500 continue; // try again. 1501 } 1502 1503 } 1504 return (isNegative)? -dValue : dValue; 1505 } 1506 } 1507 1508 /* 1509 * Take a FormattedFloatingDecimal, which we presumably just scanned in, 1510 * and find out what its value is, as a float. 1511 * This is distinct from doubleValue() to avoid the extremely 1512 * unlikely case of a double rounding error, wherein the converstion 1513 * to double has one rounding error, and the conversion of that double 1514 * to a float has another rounding error, IN THE WRONG DIRECTION, 1515 * ( because of the preference to a zero low-order bit ). 1516 */ 1517 1518 public strictfp float floatValue(){ 1519 int kDigits = Math.min( nDigits, singleMaxDecimalDigits+1 ); 1520 int iValue; 1521 float fValue; 1522 1523 // First, check for NaN and Infinity values 1524 if(digits == infinity || digits == notANumber) { 1525 if(digits == notANumber) 1526 return Float.NaN; 1527 else 1528 return (isNegative?Float.NEGATIVE_INFINITY:Float.POSITIVE_INFINITY); 1529 } 1530 else { 1531 /* 1532 * convert the lead kDigits to an integer. 1533 */ 1534 iValue = (int)digits[0]-(int)'0'; 1535 for ( int i=1; i < kDigits; i++ ){ 1536 iValue = iValue*10 + (int)digits[i]-(int)'0'; 1537 } 1538 fValue = (float)iValue; 1539 int exp = decExponent-kDigits; 1540 /* 1541 * iValue now contains an integer with the value of 1542 * the first kDigits digits of the number. 1543 * fValue contains the (float) of the same. 1544 */ 1545 1546 if ( nDigits <= singleMaxDecimalDigits ){ 1547 /* 1548 * possibly an easy case. 1549 * We know that the digits can be represented 1550 * exactly. And if the exponent isn't too outrageous, 1551 * the whole thing can be done with one operation, 1552 * thus one rounding error. 1553 * Note that all our constructors trim all leading and 1554 * trailing zeros, so simple values (including zero) 1555 * will always end up here. 1556 */ 1557 if (exp == 0 || fValue == 0.0f) 1558 return (isNegative)? -fValue : fValue; // small floating integer 1559 else if ( exp >= 0 ){ 1560 if ( exp <= singleMaxSmallTen ){ 1561 /* 1562 * Can get the answer with one operation, 1563 * thus one roundoff. 1564 */ 1565 fValue *= singleSmall10pow[exp]; 1566 return (isNegative)? -fValue : fValue; 1567 } 1568 int slop = singleMaxDecimalDigits - kDigits; 1569 if ( exp <= singleMaxSmallTen+slop ){ 1570 /* 1571 * We can multiply dValue by 10^(slop) 1572 * and it is still "small" and exact. 1573 * Then we can multiply by 10^(exp-slop) 1574 * with one rounding. 1575 */ 1576 fValue *= singleSmall10pow[slop]; 1577 fValue *= singleSmall10pow[exp-slop]; 1578 return (isNegative)? -fValue : fValue; 1579 } 1580 /* 1581 * Else we have a hard case with a positive exp. 1582 */ 1583 } else { 1584 if ( exp >= -singleMaxSmallTen ){ 1585 /* 1586 * Can get the answer in one division. 1587 */ 1588 fValue /= singleSmall10pow[-exp]; 1589 return (isNegative)? -fValue : fValue; 1590 } 1591 /* 1592 * Else we have a hard case with a negative exp. 1593 */ 1594 } 1595 } else if ( (decExponent >= nDigits) && (nDigits+decExponent <= maxDecimalDigits) ){ 1596 /* 1597 * In double-precision, this is an exact floating integer. 1598 * So we can compute to double, then shorten to float 1599 * with one round, and get the right answer. 1600 * 1601 * First, finish accumulating digits. 1602 * Then convert that integer to a double, multiply 1603 * by the appropriate power of ten, and convert to float. 1604 */ 1605 long lValue = (long)iValue; 1606 for ( int i=kDigits; i < nDigits; i++ ){ 1607 lValue = lValue*10L + (long)((int)digits[i]-(int)'0'); 1608 } 1609 double dValue = (double)lValue; 1610 exp = decExponent-nDigits; 1611 dValue *= small10pow[exp]; 1612 fValue = (float)dValue; 1613 return (isNegative)? -fValue : fValue; 1614 1615 } 1616 /* 1617 * Harder cases: 1618 * The sum of digits plus exponent is greater than 1619 * what we think we can do with one error. 1620 * 1621 * Start by weeding out obviously out-of-range 1622 * results, then convert to double and go to 1623 * common hard-case code. 1624 */ 1625 if ( decExponent > singleMaxDecimalExponent+1 ){ 1626 /* 1627 * Lets face it. This is going to be 1628 * Infinity. Cut to the chase. 1629 */ 1630 return (isNegative)? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY; 1631 } else if ( decExponent < singleMinDecimalExponent-1 ){ 1632 /* 1633 * Lets face it. This is going to be 1634 * zero. Cut to the chase. 1635 */ 1636 return (isNegative)? -0.0f : 0.0f; 1637 } 1638 1639 /* 1640 * Here, we do 'way too much work, but throwing away 1641 * our partial results, and going and doing the whole 1642 * thing as double, then throwing away half the bits that computes 1643 * when we convert back to float. 1644 * 1645 * The alternative is to reproduce the whole multiple-precision 1646 * algorythm for float precision, or to try to parameterize it 1647 * for common usage. The former will take about 400 lines of code, 1648 * and the latter I tried without success. Thus the semi-hack 1649 * answer here. 1650 */ 1651 mustSetRoundDir = !fromHex; 1652 double dValue = doubleValue(); 1653 return stickyRound( dValue ); 1654 } 1655 } 1656 1657 1658 /* 1659 * All the positive powers of 10 that can be 1660 * represented exactly in double/float. 1661 */ 1662 private static final double small10pow[] = { 1663 1.0e0, 1664 1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5, 1665 1.0e6, 1.0e7, 1.0e8, 1.0e9, 1.0e10, 1666 1.0e11, 1.0e12, 1.0e13, 1.0e14, 1.0e15, 1667 1.0e16, 1.0e17, 1.0e18, 1.0e19, 1.0e20, 1668 1.0e21, 1.0e22 1669 }; 1670 1671 private static final float singleSmall10pow[] = { 1672 1.0e0f, 1673 1.0e1f, 1.0e2f, 1.0e3f, 1.0e4f, 1.0e5f, 1674 1.0e6f, 1.0e7f, 1.0e8f, 1.0e9f, 1.0e10f 1675 }; 1676 1677 private static final double big10pow[] = { 1678 1e16, 1e32, 1e64, 1e128, 1e256 }; 1679 private static final double tiny10pow[] = { 1680 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 }; 1681 1682 private static final int maxSmallTen = small10pow.length-1; 1683 private static final int singleMaxSmallTen = singleSmall10pow.length-1; 1684 1685 private static final int small5pow[] = { 1686 1, 1687 5, 1688 5*5, 1689 5*5*5, 1690 5*5*5*5, 1691 5*5*5*5*5, 1692 5*5*5*5*5*5, 1693 5*5*5*5*5*5*5, 1694 5*5*5*5*5*5*5*5, 1695 5*5*5*5*5*5*5*5*5, 1696 5*5*5*5*5*5*5*5*5*5, 1697 5*5*5*5*5*5*5*5*5*5*5, 1698 5*5*5*5*5*5*5*5*5*5*5*5, 1699 5*5*5*5*5*5*5*5*5*5*5*5*5 1700 }; 1701 1702 1703 private static final long long5pow[] = { 1704 1L, 1705 5L, 1706 5L*5, 1707 5L*5*5, 1708 5L*5*5*5, 1709 5L*5*5*5*5, 1710 5L*5*5*5*5*5, 1711 5L*5*5*5*5*5*5, 1712 5L*5*5*5*5*5*5*5, 1713 5L*5*5*5*5*5*5*5*5, 1714 5L*5*5*5*5*5*5*5*5*5, 1715 5L*5*5*5*5*5*5*5*5*5*5, 1716 5L*5*5*5*5*5*5*5*5*5*5*5, 1717 5L*5*5*5*5*5*5*5*5*5*5*5*5, 1718 5L*5*5*5*5*5*5*5*5*5*5*5*5*5, 1719 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1720 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1721 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1722 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1723 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1724 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1725 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1726 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1727 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1728 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1729 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1730 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1731 }; 1732 1733 // approximately ceil( log2( long5pow[i] ) ) 1734 private static final int n5bits[] = { 1735 0, 1736 3, 1737 5, 1738 7, 1739 10, 1740 12, 1741 14, 1742 17, 1743 19, 1744 21, 1745 24, 1746 26, 1747 28, 1748 31, 1749 33, 1750 35, 1751 38, 1752 40, 1753 42, 1754 45, 1755 47, 1756 49, 1757 52, 1758 54, 1759 56, 1760 59, 1761 61, 1762 }; 1763 1764 private static final char infinity[] = { 'I', 'n', 'f', 'i', 'n', 'i', 't', 'y' }; 1765 private static final char notANumber[] = { 'N', 'a', 'N' }; 1766 private static final char zero[] = { '0', '0', '0', '0', '0', '0', '0', '0' }; 1767 } 1768