1 /* 2 * Copyright (c) 1996, 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 FloatingDecimal{ 34 boolean isExceptional; 35 boolean isNegative; 36 int decExponent; 37 char digits[]; 38 int nDigits; 39 int bigIntExp; 40 int bigIntNBits; 41 boolean mustSetRoundDir = false; 42 boolean fromHex = false; 43 int roundDir = 0; // set by doubleValue 44 45 private FloatingDecimal( boolean negSign, int decExponent, char []digits, int n, boolean e ) 46 { 47 isNegative = negSign; 48 isExceptional = e; 49 this.decExponent = decExponent; 50 this.digits = digits; 51 this.nDigits = n; 52 } 53 54 /* 55 * Constants of the implementation 56 * Most are IEEE-754 related. 57 * (There are more really boring constants at the end.) 58 */ 59 static final long signMask = 0x8000000000000000L; 60 static final long expMask = 0x7ff0000000000000L; 61 static final long fractMask= ~(signMask|expMask); 62 static final int expShift = 52; 63 static final int expBias = 1023; 64 static final long fractHOB = ( 1L<<expShift ); // assumed High-Order bit 65 static final long expOne = ((long)expBias)<<expShift; // exponent of 1.0 66 static final int maxSmallBinExp = 62; 67 static final int minSmallBinExp = -( 63 / 3 ); 68 static final int maxDecimalDigits = 15; 69 static final int maxDecimalExponent = 308; 70 static final int minDecimalExponent = -324; 71 static final int bigDecimalExponent = 324; // i.e. abs(minDecimalExponent) 72 73 static final long highbyte = 0xff00000000000000L; 74 static final long highbit = 0x8000000000000000L; 75 static final long lowbytes = ~highbyte; 76 77 static final int singleSignMask = 0x80000000; 78 static final int singleExpMask = 0x7f800000; 79 static final int singleFractMask = ~(singleSignMask|singleExpMask); 80 static final int singleExpShift = 23; 81 static final int singleFractHOB = 1<<singleExpShift; 82 static final int singleExpBias = 127; 83 static final int singleMaxDecimalDigits = 7; 84 static final int singleMaxDecimalExponent = 38; 85 static final int singleMinDecimalExponent = -45; 86 87 static final int intDecimalDigits = 9; 88 89 90 /* 91 * count number of bits from high-order 1 bit to low-order 1 bit, 92 * inclusive. 93 */ 94 private static int 95 countBits( long v ){ 96 // 97 // the strategy is to shift until we get a non-zero sign bit 98 // then shift until we have no bits left, counting the difference. 99 // we do byte shifting as a hack. Hope it helps. 100 // 101 if ( v == 0L ) return 0; 102 103 while ( ( v & highbyte ) == 0L ){ 104 v <<= 8; 105 } 106 while ( v > 0L ) { // i.e. while ((v&highbit) == 0L ) 107 v <<= 1; 108 } 109 110 int n = 0; 111 while (( v & lowbytes ) != 0L ){ 112 v <<= 8; 113 n += 8; 114 } 115 while ( v != 0L ){ 116 v <<= 1; 117 n += 1; 118 } 119 return n; 120 } 121 122 /* 123 * Keep big powers of 5 handy for future reference. 124 */ 125 private static FDBigInt b5p[]; 126 127 private static synchronized FDBigInt 128 big5pow( int p ){ 129 assert p >= 0 : p; // negative power of 5 130 if ( b5p == null ){ 131 b5p = new FDBigInt[ p+1 ]; 132 }else if (b5p.length <= p ){ 133 FDBigInt t[] = new FDBigInt[ p+1 ]; 134 System.arraycopy( b5p, 0, t, 0, b5p.length ); 135 b5p = t; 136 } 137 if ( b5p[p] != null ) 138 return b5p[p]; 139 else if ( p < small5pow.length ) 140 return b5p[p] = new FDBigInt( small5pow[p] ); 141 else if ( p < long5pow.length ) 142 return b5p[p] = new FDBigInt( long5pow[p] ); 143 else { 144 // construct the value. 145 // recursively. 146 int q, r; 147 // in order to compute 5^p, 148 // compute its square root, 5^(p/2) and square. 149 // or, let q = p / 2, r = p -q, then 150 // 5^p = 5^(q+r) = 5^q * 5^r 151 q = p >> 1; 152 r = p - q; 153 FDBigInt bigq = b5p[q]; 154 if ( bigq == null ) 155 bigq = big5pow ( q ); 156 if ( r < small5pow.length ){ 157 return (b5p[p] = bigq.mult( small5pow[r] ) ); 158 }else{ 159 FDBigInt bigr = b5p[ r ]; 160 if ( bigr == null ) 161 bigr = big5pow( r ); 162 return (b5p[p] = bigq.mult( bigr ) ); 163 } 164 } 165 } 166 167 // 168 // a common operation 169 // 170 private static FDBigInt 171 multPow52( FDBigInt v, int p5, int p2 ){ 172 if ( p5 != 0 ){ 173 if ( p5 < small5pow.length ){ 174 v = v.mult( small5pow[p5] ); 175 } else { 176 v = v.mult( big5pow( p5 ) ); 177 } 178 } 179 if ( p2 != 0 ){ 180 v.lshiftMe( p2 ); 181 } 182 return v; 183 } 184 185 // 186 // another common operation 187 // 188 private static FDBigInt 189 constructPow52( int p5, int p2 ){ 190 FDBigInt v = new FDBigInt( big5pow( p5 ) ); 191 if ( p2 != 0 ){ 192 v.lshiftMe( p2 ); 193 } 194 return v; 195 } 196 197 /* 198 * Make a floating double into a FDBigInt. 199 * This could also be structured as a FDBigInt 200 * constructor, but we'd have to build a lot of knowledge 201 * about floating-point representation into it, and we don't want to. 202 * 203 * AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES 204 * bigIntExp and bigIntNBits 205 * 206 */ 207 private FDBigInt 208 doubleToBigInt( double dval ){ 209 long lbits = Double.doubleToLongBits( dval ) & ~signMask; 210 int binexp = (int)(lbits >>> expShift); 211 lbits &= fractMask; 212 if ( binexp > 0 ){ 213 lbits |= fractHOB; 214 } else { 215 assert lbits != 0L : lbits; // doubleToBigInt(0.0) 216 binexp +=1; 217 while ( (lbits & fractHOB ) == 0L){ 218 lbits <<= 1; 219 binexp -= 1; 220 } 221 } 222 binexp -= expBias; 223 int nbits = countBits( lbits ); 224 /* 225 * We now know where the high-order 1 bit is, 226 * and we know how many there are. 227 */ 228 int lowOrderZeros = expShift+1-nbits; 229 lbits >>>= lowOrderZeros; 230 231 bigIntExp = binexp+1-nbits; 232 bigIntNBits = nbits; 233 return new FDBigInt( lbits ); 234 } 235 236 /* 237 * Compute a number that is the ULP of the given value, 238 * for purposes of addition/subtraction. Generally easy. 239 * More difficult if subtracting and the argument 240 * is a normalized a power of 2, as the ULP changes at these points. 241 */ 242 private static double ulp( double dval, boolean subtracting ){ 243 long lbits = Double.doubleToLongBits( dval ) & ~signMask; 244 int binexp = (int)(lbits >>> expShift); 245 double ulpval; 246 if ( subtracting && ( binexp >= expShift ) && ((lbits&fractMask) == 0L) ){ 247 // for subtraction from normalized, powers of 2, 248 // use next-smaller exponent 249 binexp -= 1; 250 } 251 if ( binexp > expShift ){ 252 ulpval = Double.longBitsToDouble( ((long)(binexp-expShift))<<expShift ); 253 } else if ( binexp == 0 ){ 254 ulpval = Double.MIN_VALUE; 255 } else { 256 ulpval = Double.longBitsToDouble( 1L<<(binexp-1) ); 257 } 258 if ( subtracting ) ulpval = - ulpval; 259 260 return ulpval; 261 } 262 263 /* 264 * Round a double to a float. 265 * In addition to the fraction bits of the double, 266 * look at the class instance variable roundDir, 267 * which should help us avoid double-rounding error. 268 * roundDir was set in hardValueOf if the estimate was 269 * close enough, but not exact. It tells us which direction 270 * of rounding is preferred. 271 */ 272 float 273 stickyRound( double dval ){ 274 long lbits = Double.doubleToLongBits( dval ); 275 long binexp = lbits & expMask; 276 if ( binexp == 0L || binexp == expMask ){ 277 // what we have here is special. 278 // don't worry, the right thing will happen. 279 return (float) dval; 280 } 281 lbits += (long)roundDir; // hack-o-matic. 282 return (float)Double.longBitsToDouble( lbits ); 283 } 284 285 286 /* 287 * This is the easy subcase -- 288 * all the significant bits, after scaling, are held in lvalue. 289 * negSign and decExponent tell us what processing and scaling 290 * has already been done. Exceptional cases have already been 291 * stripped out. 292 * In particular: 293 * lvalue is a finite number (not Inf, nor NaN) 294 * lvalue > 0L (not zero, nor negative). 295 * 296 * The only reason that we develop the digits here, rather than 297 * calling on Long.toString() is that we can do it a little faster, 298 * and besides want to treat trailing 0s specially. If Long.toString 299 * changes, we should re-evaluate this strategy! 300 */ 301 private void 302 developLongDigits( int decExponent, long lvalue, long insignificant ){ 303 char digits[]; 304 int ndigits; 305 int digitno; 306 int c; 307 // 308 // Discard non-significant low-order bits, while rounding, 309 // up to insignificant value. 310 int i; 311 for ( i = 0; insignificant >= 10L; i++ ) 312 insignificant /= 10L; 313 if ( i != 0 ){ 314 long pow10 = long5pow[i] << i; // 10^i == 5^i * 2^i; 315 long residue = lvalue % pow10; 316 lvalue /= pow10; 317 decExponent += i; 318 if ( residue >= (pow10>>1) ){ 319 // round up based on the low-order bits we're discarding 320 lvalue++; 321 } 322 } 323 if ( lvalue <= Integer.MAX_VALUE ){ 324 assert lvalue > 0L : lvalue; // lvalue <= 0 325 // even easier subcase! 326 // can do int arithmetic rather than long! 327 int ivalue = (int)lvalue; 328 ndigits = 10; 329 digits = (char[])(perThreadBuffer.get()); 330 digitno = ndigits-1; 331 c = ivalue%10; 332 ivalue /= 10; 333 while ( c == 0 ){ 334 decExponent++; 335 c = ivalue%10; 336 ivalue /= 10; 337 } 338 while ( ivalue != 0){ 339 digits[digitno--] = (char)(c+'0'); 340 decExponent++; 341 c = ivalue%10; 342 ivalue /= 10; 343 } 344 digits[digitno] = (char)(c+'0'); 345 } else { 346 // same algorithm as above (same bugs, too ) 347 // but using long arithmetic. 348 ndigits = 20; 349 digits = (char[])(perThreadBuffer.get()); 350 digitno = ndigits-1; 351 c = (int)(lvalue%10L); 352 lvalue /= 10L; 353 while ( c == 0 ){ 354 decExponent++; 355 c = (int)(lvalue%10L); 356 lvalue /= 10L; 357 } 358 while ( lvalue != 0L ){ 359 digits[digitno--] = (char)(c+'0'); 360 decExponent++; 361 c = (int)(lvalue%10L); 362 lvalue /= 10; 363 } 364 digits[digitno] = (char)(c+'0'); 365 } 366 char result []; 367 ndigits -= digitno; 368 result = new char[ ndigits ]; 369 System.arraycopy( digits, digitno, result, 0, ndigits ); 370 this.digits = result; 371 this.decExponent = decExponent+1; 372 this.nDigits = ndigits; 373 } 374 375 // 376 // add one to the least significant digit. 377 // in the unlikely event there is a carry out, 378 // deal with it. 379 // assert that this will only happen where there 380 // is only one digit, e.g. (float)1e-44 seems to do it. 381 // 382 private void 383 roundup(){ 384 int i; 385 int q = digits[ i = (nDigits-1)]; 386 if ( q == '9' ){ 387 while ( q == '9' && i > 0 ){ 388 digits[i] = '0'; 389 q = digits[--i]; 390 } 391 if ( q == '9' ){ 392 // carryout! High-order 1, rest 0s, larger exp. 393 decExponent += 1; 394 digits[0] = '1'; 395 return; 396 } 397 // else fall through. 398 } 399 digits[i] = (char)(q+1); 400 } 401 402 /* 403 * FIRST IMPORTANT CONSTRUCTOR: DOUBLE 404 */ 405 public FloatingDecimal( double d ) 406 { 407 long dBits = Double.doubleToLongBits( d ); 408 long fractBits; 409 int binExp; 410 int nSignificantBits; 411 412 // discover and delete sign 413 if ( (dBits&signMask) != 0 ){ 414 isNegative = true; 415 dBits ^= signMask; 416 } else { 417 isNegative = false; 418 } 419 // Begin to unpack 420 // Discover obvious special cases of NaN and Infinity. 421 binExp = (int)( (dBits&expMask) >> expShift ); 422 fractBits = dBits&fractMask; 423 if ( binExp == (int)(expMask>>expShift) ) { 424 isExceptional = true; 425 if ( fractBits == 0L ){ 426 digits = infinity; 427 } else { 428 digits = notANumber; 429 isNegative = false; // NaN has no sign! 430 } 431 nDigits = digits.length; 432 return; 433 } 434 isExceptional = false; 435 // Finish unpacking 436 // Normalize denormalized numbers. 437 // Insert assumed high-order bit for normalized numbers. 438 // Subtract exponent bias. 439 if ( binExp == 0 ){ 440 if ( fractBits == 0L ){ 441 // not a denorm, just a 0! 442 decExponent = 0; 443 digits = zero; 444 nDigits = 1; 445 return; 446 } 447 while ( (fractBits&fractHOB) == 0L ){ 448 fractBits <<= 1; 449 binExp -= 1; 450 } 451 nSignificantBits = expShift + binExp +1; // recall binExp is - shift count. 452 binExp += 1; 453 } else { 454 fractBits |= fractHOB; 455 nSignificantBits = expShift+1; 456 } 457 binExp -= expBias; 458 // call the routine that actually does all the hard work. 459 dtoa( binExp, fractBits, nSignificantBits ); 460 } 461 462 /* 463 * SECOND IMPORTANT CONSTRUCTOR: SINGLE 464 */ 465 public FloatingDecimal( float f ) 466 { 467 int fBits = Float.floatToIntBits( f ); 468 int fractBits; 469 int binExp; 470 int nSignificantBits; 471 472 // discover and delete sign 473 if ( (fBits&singleSignMask) != 0 ){ 474 isNegative = true; 475 fBits ^= singleSignMask; 476 } else { 477 isNegative = false; 478 } 479 // Begin to unpack 480 // Discover obvious special cases of NaN and Infinity. 481 binExp = (int)( (fBits&singleExpMask) >> singleExpShift ); 482 fractBits = fBits&singleFractMask; 483 if ( binExp == (int)(singleExpMask>>singleExpShift) ) { 484 isExceptional = true; 485 if ( fractBits == 0L ){ 486 digits = infinity; 487 } else { 488 digits = notANumber; 489 isNegative = false; // NaN has no sign! 490 } 491 nDigits = digits.length; 492 return; 493 } 494 isExceptional = false; 495 // Finish unpacking 496 // Normalize denormalized numbers. 497 // Insert assumed high-order bit for normalized numbers. 498 // Subtract exponent bias. 499 if ( binExp == 0 ){ 500 if ( fractBits == 0 ){ 501 // not a denorm, just a 0! 502 decExponent = 0; 503 digits = zero; 504 nDigits = 1; 505 return; 506 } 507 while ( (fractBits&singleFractHOB) == 0 ){ 508 fractBits <<= 1; 509 binExp -= 1; 510 } 511 nSignificantBits = singleExpShift + binExp +1; // recall binExp is - shift count. 512 binExp += 1; 513 } else { 514 fractBits |= singleFractHOB; 515 nSignificantBits = singleExpShift+1; 516 } 517 binExp -= singleExpBias; 518 // call the routine that actually does all the hard work. 519 dtoa( binExp, ((long)fractBits)<<(expShift-singleExpShift), nSignificantBits ); 520 } 521 522 private void 523 dtoa( int binExp, long fractBits, int nSignificantBits ) 524 { 525 int nFractBits; // number of significant bits of fractBits; 526 int nTinyBits; // number of these to the right of the point. 527 int decExp; 528 529 // Examine number. Determine if it is an easy case, 530 // which we can do pretty trivially using float/long conversion, 531 // or whether we must do real work. 532 nFractBits = countBits( fractBits ); 533 nTinyBits = Math.max( 0, nFractBits - binExp - 1 ); 534 if ( binExp <= maxSmallBinExp && binExp >= minSmallBinExp ){ 535 // Look more closely at the number to decide if, 536 // with scaling by 10^nTinyBits, the result will fit in 537 // a long. 538 if ( (nTinyBits < long5pow.length) && ((nFractBits + n5bits[nTinyBits]) < 64 ) ){ 539 /* 540 * We can do this: 541 * take the fraction bits, which are normalized. 542 * (a) nTinyBits == 0: Shift left or right appropriately 543 * to align the binary point at the extreme right, i.e. 544 * where a long int point is expected to be. The integer 545 * result is easily converted to a string. 546 * (b) nTinyBits > 0: Shift right by expShift-nFractBits, 547 * which effectively converts to long and scales by 548 * 2^nTinyBits. Then multiply by 5^nTinyBits to 549 * complete the scaling. We know this won't overflow 550 * because we just counted the number of bits necessary 551 * in the result. The integer you get from this can 552 * then be converted to a string pretty easily. 553 */ 554 long halfULP; 555 if ( nTinyBits == 0 ) { 556 if ( binExp > nSignificantBits ){ 557 halfULP = 1L << ( binExp-nSignificantBits-1); 558 } else { 559 halfULP = 0L; 560 } 561 if ( binExp >= expShift ){ 562 fractBits <<= (binExp-expShift); 563 } else { 564 fractBits >>>= (expShift-binExp) ; 565 } 566 developLongDigits( 0, fractBits, halfULP ); 567 return; 568 } 569 /* 570 * The following causes excess digits to be printed 571 * out in the single-float case. Our manipulation of 572 * halfULP here is apparently not correct. If we 573 * better understand how this works, perhaps we can 574 * use this special case again. But for the time being, 575 * we do not. 576 * else { 577 * fractBits >>>= expShift+1-nFractBits; 578 * fractBits *= long5pow[ nTinyBits ]; 579 * halfULP = long5pow[ nTinyBits ] >> (1+nSignificantBits-nFractBits); 580 * developLongDigits( -nTinyBits, fractBits, halfULP ); 581 * return; 582 * } 583 */ 584 } 585 } 586 /* 587 * This is the hard case. We are going to compute large positive 588 * integers B and S and integer decExp, s.t. 589 * d = ( B / S ) * 10^decExp 590 * 1 <= B / S < 10 591 * Obvious choices are: 592 * decExp = floor( log10(d) ) 593 * B = d * 2^nTinyBits * 10^max( 0, -decExp ) 594 * S = 10^max( 0, decExp) * 2^nTinyBits 595 * (noting that nTinyBits has already been forced to non-negative) 596 * I am also going to compute a large positive integer 597 * M = (1/2^nSignificantBits) * 2^nTinyBits * 10^max( 0, -decExp ) 598 * i.e. M is (1/2) of the ULP of d, scaled like B. 599 * When we iterate through dividing B/S and picking off the 600 * quotient bits, we will know when to stop when the remainder 601 * is <= M. 602 * 603 * We keep track of powers of 2 and powers of 5. 604 */ 605 606 /* 607 * Estimate decimal exponent. (If it is small-ish, 608 * we could double-check.) 609 * 610 * First, scale the mantissa bits such that 1 <= d2 < 2. 611 * We are then going to estimate 612 * log10(d2) ~=~ (d2-1.5)/1.5 + log(1.5) 613 * and so we can estimate 614 * log10(d) ~=~ log10(d2) + binExp * log10(2) 615 * take the floor and call it decExp. 616 * FIXME -- use more precise constants here. It costs no more. 617 */ 618 double d2 = Double.longBitsToDouble( 619 expOne | ( fractBits &~ fractHOB ) ); 620 decExp = (int)Math.floor( 621 (d2-1.5D)*0.289529654D + 0.176091259 + (double)binExp * 0.301029995663981 ); 622 int B2, B5; // powers of 2 and powers of 5, respectively, in B 623 int S2, S5; // powers of 2 and powers of 5, respectively, in S 624 int M2, M5; // powers of 2 and powers of 5, respectively, in M 625 int Bbits; // binary digits needed to represent B, approx. 626 int tenSbits; // binary digits needed to represent 10*S, approx. 627 FDBigInt Sval, Bval, Mval; 628 629 B5 = Math.max( 0, -decExp ); 630 B2 = B5 + nTinyBits + binExp; 631 632 S5 = Math.max( 0, decExp ); 633 S2 = S5 + nTinyBits; 634 635 M5 = B5; 636 M2 = B2 - nSignificantBits; 637 638 /* 639 * the long integer fractBits contains the (nFractBits) interesting 640 * bits from the mantissa of d ( hidden 1 added if necessary) followed 641 * by (expShift+1-nFractBits) zeros. In the interest of compactness, 642 * I will shift out those zeros before turning fractBits into a 643 * FDBigInt. The resulting whole number will be 644 * d * 2^(nFractBits-1-binExp). 645 */ 646 fractBits >>>= (expShift+1-nFractBits); 647 B2 -= nFractBits-1; 648 int common2factor = Math.min( B2, S2 ); 649 B2 -= common2factor; 650 S2 -= common2factor; 651 M2 -= common2factor; 652 653 /* 654 * HACK!! For exact powers of two, the next smallest number 655 * is only half as far away as we think (because the meaning of 656 * ULP changes at power-of-two bounds) for this reason, we 657 * hack M2. Hope this works. 658 */ 659 if ( nFractBits == 1 ) 660 M2 -= 1; 661 662 if ( M2 < 0 ){ 663 // oops. 664 // since we cannot scale M down far enough, 665 // we must scale the other values up. 666 B2 -= M2; 667 S2 -= M2; 668 M2 = 0; 669 } 670 /* 671 * Construct, Scale, iterate. 672 * Some day, we'll write a stopping test that takes 673 * account of the asymmetry of the spacing of floating-point 674 * numbers below perfect powers of 2 675 * 26 Sept 96 is not that day. 676 * So we use a symmetric test. 677 */ 678 char digits[] = this.digits = new char[18]; 679 int ndigit = 0; 680 boolean low, high; 681 long lowDigitDifference; 682 int q; 683 684 /* 685 * Detect the special cases where all the numbers we are about 686 * to compute will fit in int or long integers. 687 * In these cases, we will avoid doing FDBigInt arithmetic. 688 * We use the same algorithms, except that we "normalize" 689 * our FDBigInts before iterating. This is to make division easier, 690 * as it makes our fist guess (quotient of high-order words) 691 * more accurate! 692 * 693 * Some day, we'll write a stopping test that takes 694 * account of the asymmetry of the spacing of floating-point 695 * numbers below perfect powers of 2 696 * 26 Sept 96 is not that day. 697 * So we use a symmetric test. 698 */ 699 Bbits = nFractBits + B2 + (( B5 < n5bits.length )? n5bits[B5] : ( B5*3 )); 700 tenSbits = S2+1 + (( (S5+1) < n5bits.length )? n5bits[(S5+1)] : ( (S5+1)*3 )); 701 if ( Bbits < 64 && tenSbits < 64){ 702 if ( Bbits < 32 && tenSbits < 32){ 703 // wa-hoo! They're all ints! 704 int b = ((int)fractBits * small5pow[B5] ) << B2; 705 int s = small5pow[S5] << S2; 706 int m = small5pow[M5] << M2; 707 int tens = s * 10; 708 /* 709 * Unroll the first iteration. If our decExp estimate 710 * was too high, our first quotient will be zero. In this 711 * case, we discard it and decrement decExp. 712 */ 713 ndigit = 0; 714 q = b / s; 715 b = 10 * ( b % s ); 716 m *= 10; 717 low = (b < m ); 718 high = (b+m > tens ); 719 assert q < 10 : q; // excessively large digit 720 if ( (q == 0) && ! high ){ 721 // oops. Usually ignore leading zero. 722 decExp--; 723 } else { 724 digits[ndigit++] = (char)('0' + q); 725 } 726 /* 727 * HACK! Java spec sez that we always have at least 728 * one digit after the . in either F- or E-form output. 729 * Thus we will need more than one digit if we're using 730 * E-form 731 */ 732 if ( decExp < -3 || decExp >= 8 ){ 733 high = low = false; 734 } 735 while( ! low && ! high ){ 736 q = b / s; 737 b = 10 * ( b % s ); 738 m *= 10; 739 assert q < 10 : q; // excessively large digit 740 if ( m > 0L ){ 741 low = (b < m ); 742 high = (b+m > tens ); 743 } else { 744 // hack -- m might overflow! 745 // in this case, it is certainly > b, 746 // which won't 747 // and b+m > tens, too, since that has overflowed 748 // either! 749 low = true; 750 high = true; 751 } 752 digits[ndigit++] = (char)('0' + q); 753 } 754 lowDigitDifference = (b<<1) - tens; 755 } else { 756 // still good! they're all longs! 757 long b = (fractBits * long5pow[B5] ) << B2; 758 long s = long5pow[S5] << S2; 759 long m = long5pow[M5] << M2; 760 long tens = s * 10L; 761 /* 762 * Unroll the first iteration. If our decExp estimate 763 * was too high, our first quotient will be zero. In this 764 * case, we discard it and decrement decExp. 765 */ 766 ndigit = 0; 767 q = (int) ( b / s ); 768 b = 10L * ( b % s ); 769 m *= 10L; 770 low = (b < m ); 771 high = (b+m > tens ); 772 assert q < 10 : q; // excessively large digit 773 if ( (q == 0) && ! high ){ 774 // oops. Usually ignore leading zero. 775 decExp--; 776 } else { 777 digits[ndigit++] = (char)('0' + q); 778 } 779 /* 780 * HACK! Java spec sez that we always have at least 781 * one digit after the . in either F- or E-form output. 782 * Thus we will need more than one digit if we're using 783 * E-form 784 */ 785 if ( decExp < -3 || decExp >= 8 ){ 786 high = low = false; 787 } 788 while( ! low && ! high ){ 789 q = (int) ( b / s ); 790 b = 10 * ( b % s ); 791 m *= 10; 792 assert q < 10 : q; // excessively large digit 793 if ( m > 0L ){ 794 low = (b < m ); 795 high = (b+m > tens ); 796 } else { 797 // hack -- m might overflow! 798 // in this case, it is certainly > b, 799 // which won't 800 // and b+m > tens, too, since that has overflowed 801 // either! 802 low = true; 803 high = true; 804 } 805 digits[ndigit++] = (char)('0' + q); 806 } 807 lowDigitDifference = (b<<1) - tens; 808 } 809 } else { 810 FDBigInt tenSval; 811 int shiftBias; 812 813 /* 814 * We really must do FDBigInt arithmetic. 815 * Fist, construct our FDBigInt initial values. 816 */ 817 Bval = multPow52( new FDBigInt( fractBits ), B5, B2 ); 818 Sval = constructPow52( S5, S2 ); 819 Mval = constructPow52( M5, M2 ); 820 821 822 // normalize so that division works better 823 Bval.lshiftMe( shiftBias = Sval.normalizeMe() ); 824 Mval.lshiftMe( shiftBias ); 825 tenSval = Sval.mult( 10 ); 826 /* 827 * Unroll the first iteration. If our decExp estimate 828 * was too high, our first quotient will be zero. In this 829 * case, we discard it and decrement decExp. 830 */ 831 ndigit = 0; 832 q = Bval.quoRemIteration( Sval ); 833 Mval = Mval.mult( 10 ); 834 low = (Bval.cmp( Mval ) < 0); 835 high = (Bval.add( Mval ).cmp( tenSval ) > 0 ); 836 assert q < 10 : q; // excessively large digit 837 if ( (q == 0) && ! high ){ 838 // oops. Usually ignore leading zero. 839 decExp--; 840 } else { 841 digits[ndigit++] = (char)('0' + q); 842 } 843 /* 844 * HACK! Java spec sez that we always have at least 845 * one digit after the . in either F- or E-form output. 846 * Thus we will need more than one digit if we're using 847 * E-form 848 */ 849 if ( decExp < -3 || decExp >= 8 ){ 850 high = low = false; 851 } 852 while( ! low && ! high ){ 853 q = Bval.quoRemIteration( Sval ); 854 Mval = Mval.mult( 10 ); 855 assert q < 10 : q; // excessively large digit 856 low = (Bval.cmp( Mval ) < 0); 857 high = (Bval.add( Mval ).cmp( tenSval ) > 0 ); 858 digits[ndigit++] = (char)('0' + q); 859 } 860 if ( high && low ){ 861 Bval.lshiftMe(1); 862 lowDigitDifference = Bval.cmp(tenSval); 863 } else 864 lowDigitDifference = 0L; // this here only for flow analysis! 865 } 866 this.decExponent = decExp+1; 867 this.digits = digits; 868 this.nDigits = ndigit; 869 /* 870 * Last digit gets rounded based on stopping condition. 871 */ 872 if ( high ){ 873 if ( low ){ 874 if ( lowDigitDifference == 0L ){ 875 // it's a tie! 876 // choose based on which digits we like. 877 if ( (digits[nDigits-1]&1) != 0 ) roundup(); 878 } else if ( lowDigitDifference > 0 ){ 879 roundup(); 880 } 881 } else { 882 roundup(); 883 } 884 } 885 } 886 887 public String 888 toString(){ 889 // most brain-dead version 890 StringBuffer result = new StringBuffer( nDigits+8 ); 891 if ( isNegative ){ result.append( '-' ); } 892 if ( isExceptional ){ 893 result.append( digits, 0, nDigits ); 894 } else { 895 result.append( "0."); 896 result.append( digits, 0, nDigits ); 897 result.append('e'); 898 result.append( decExponent ); 899 } 900 return new String(result); 901 } 902 903 public String toJavaFormatString() { 904 char result[] = (char[])(perThreadBuffer.get()); 905 int i = getChars(result); 906 return new String(result, 0, i); 907 } 908 909 private int getChars(char[] result) { 910 assert nDigits <= 19 : nDigits; // generous bound on size of nDigits 911 int i = 0; 912 if (isNegative) { result[0] = '-'; i = 1; } 913 if (isExceptional) { 914 System.arraycopy(digits, 0, result, i, nDigits); 915 i += nDigits; 916 } else { 917 if (decExponent > 0 && decExponent < 8) { 918 // print digits.digits. 919 int charLength = Math.min(nDigits, decExponent); 920 System.arraycopy(digits, 0, result, i, charLength); 921 i += charLength; 922 if (charLength < decExponent) { 923 charLength = decExponent-charLength; 924 System.arraycopy(zero, 0, result, i, charLength); 925 i += charLength; 926 result[i++] = '.'; 927 result[i++] = '0'; 928 } else { 929 result[i++] = '.'; 930 if (charLength < nDigits) { 931 int t = nDigits - charLength; 932 System.arraycopy(digits, charLength, result, i, t); 933 i += t; 934 } else { 935 result[i++] = '0'; 936 } 937 } 938 } else if (decExponent <=0 && decExponent > -3) { 939 result[i++] = '0'; 940 result[i++] = '.'; 941 if (decExponent != 0) { 942 System.arraycopy(zero, 0, result, i, -decExponent); 943 i -= decExponent; 944 } 945 System.arraycopy(digits, 0, result, i, nDigits); 946 i += nDigits; 947 } else { 948 result[i++] = digits[0]; 949 result[i++] = '.'; 950 if (nDigits > 1) { 951 System.arraycopy(digits, 1, result, i, nDigits-1); 952 i += nDigits-1; 953 } else { 954 result[i++] = '0'; 955 } 956 result[i++] = 'E'; 957 int e; 958 if (decExponent <= 0) { 959 result[i++] = '-'; 960 e = -decExponent+1; 961 } else { 962 e = decExponent-1; 963 } 964 // decExponent has 1, 2, or 3, digits 965 if (e <= 9) { 966 result[i++] = (char)(e+'0'); 967 } else if (e <= 99) { 968 result[i++] = (char)(e/10 +'0'); 969 result[i++] = (char)(e%10 + '0'); 970 } else { 971 result[i++] = (char)(e/100+'0'); 972 e %= 100; 973 result[i++] = (char)(e/10+'0'); 974 result[i++] = (char)(e%10 + '0'); 975 } 976 } 977 } 978 return i; 979 } 980 981 // Per-thread buffer for string/stringbuffer conversion 982 private static ThreadLocal perThreadBuffer = new ThreadLocal() { 983 protected synchronized Object initialValue() { 984 return new char[26]; 985 } 986 }; 987 988 public void appendTo(Appendable buf) { 989 char result[] = (char[])(perThreadBuffer.get()); 990 int i = getChars(result); 991 if (buf instanceof StringBuilder) 992 ((StringBuilder) buf).append(result, 0, i); 993 else if (buf instanceof StringBuffer) 994 ((StringBuffer) buf).append(result, 0, i); 995 else 996 assert false; 997 } 998 999 public static FloatingDecimal 1000 readJavaFormatString( String in ) throws NumberFormatException { 1001 boolean isNegative = false; 1002 boolean signSeen = false; 1003 int decExp; 1004 char c; 1005 1006 parseNumber: 1007 try{ 1008 in = in.trim(); // don't fool around with white space. 1009 // throws NullPointerException if null 1010 int l = in.length(); 1011 if ( l == 0 ) throw new NumberFormatException("empty String"); 1012 int i = 0; 1013 switch ( c = in.charAt( i ) ){ 1014 case '-': 1015 isNegative = true; 1016 //FALLTHROUGH 1017 case '+': 1018 i++; 1019 signSeen = true; 1020 } 1021 1022 // Check for NaN and Infinity strings 1023 c = in.charAt(i); 1024 if(c == 'N' || c == 'I') { // possible NaN or infinity 1025 boolean potentialNaN = false; 1026 char targetChars[] = null; // char array of "NaN" or "Infinity" 1027 1028 if(c == 'N') { 1029 targetChars = notANumber; 1030 potentialNaN = true; 1031 } else { 1032 targetChars = infinity; 1033 } 1034 1035 // compare Input string to "NaN" or "Infinity" 1036 int j = 0; 1037 while(i < l && j < targetChars.length) { 1038 if(in.charAt(i) == targetChars[j]) { 1039 i++; j++; 1040 } 1041 else // something is amiss, throw exception 1042 break parseNumber; 1043 } 1044 1045 // For the candidate string to be a NaN or infinity, 1046 // all characters in input string and target char[] 1047 // must be matched ==> j must equal targetChars.length 1048 // and i must equal l 1049 if( (j == targetChars.length) && (i == l) ) { // return NaN or infinity 1050 return (potentialNaN ? new FloatingDecimal(Double.NaN) // NaN has no sign 1051 : new FloatingDecimal(isNegative? 1052 Double.NEGATIVE_INFINITY: 1053 Double.POSITIVE_INFINITY)) ; 1054 } 1055 else { // something went wrong, throw exception 1056 break parseNumber; 1057 } 1058 1059 } else if (c == '0') { // check for hexadecimal floating-point number 1060 if (l > i+1 ) { 1061 char ch = in.charAt(i+1); 1062 if (ch == 'x' || ch == 'X' ) // possible hex string 1063 return parseHexString(in); 1064 } 1065 } // look for and process decimal floating-point string 1066 1067 char[] digits = new char[ l ]; 1068 int nDigits= 0; 1069 boolean decSeen = false; 1070 int decPt = 0; 1071 int nLeadZero = 0; 1072 int nTrailZero= 0; 1073 digitLoop: 1074 while ( i < l ){ 1075 switch ( c = in.charAt( i ) ){ 1076 case '0': 1077 if ( nDigits > 0 ){ 1078 nTrailZero += 1; 1079 } else { 1080 nLeadZero += 1; 1081 } 1082 break; // out of switch. 1083 case '1': 1084 case '2': 1085 case '3': 1086 case '4': 1087 case '5': 1088 case '6': 1089 case '7': 1090 case '8': 1091 case '9': 1092 while ( nTrailZero > 0 ){ 1093 digits[nDigits++] = '0'; 1094 nTrailZero -= 1; 1095 } 1096 digits[nDigits++] = c; 1097 break; // out of switch. 1098 case '.': 1099 if ( decSeen ){ 1100 // already saw one ., this is the 2nd. 1101 throw new NumberFormatException("multiple points"); 1102 } 1103 decPt = i; 1104 if ( signSeen ){ 1105 decPt -= 1; 1106 } 1107 decSeen = true; 1108 break; // out of switch. 1109 default: 1110 break digitLoop; 1111 } 1112 i++; 1113 } 1114 /* 1115 * At this point, we've scanned all the digits and decimal 1116 * point we're going to see. Trim off leading and trailing 1117 * zeros, which will just confuse us later, and adjust 1118 * our initial decimal exponent accordingly. 1119 * To review: 1120 * we have seen i total characters. 1121 * nLeadZero of them were zeros before any other digits. 1122 * nTrailZero of them were zeros after any other digits. 1123 * if ( decSeen ), then a . was seen after decPt characters 1124 * ( including leading zeros which have been discarded ) 1125 * nDigits characters were neither lead nor trailing 1126 * zeros, nor point 1127 */ 1128 /* 1129 * special hack: if we saw no non-zero digits, then the 1130 * answer is zero! 1131 * Unfortunately, we feel honor-bound to keep parsing! 1132 */ 1133 if ( nDigits == 0 ){ 1134 digits = zero; 1135 nDigits = 1; 1136 if ( nLeadZero == 0 ){ 1137 // we saw NO DIGITS AT ALL, 1138 // not even a crummy 0! 1139 // this is not allowed. 1140 break parseNumber; // go throw exception 1141 } 1142 1143 } 1144 1145 /* Our initial exponent is decPt, adjusted by the number of 1146 * discarded zeros. Or, if there was no decPt, 1147 * then its just nDigits adjusted by discarded trailing zeros. 1148 */ 1149 if ( decSeen ){ 1150 decExp = decPt - nLeadZero; 1151 } else { 1152 decExp = nDigits+nTrailZero; 1153 } 1154 1155 /* 1156 * Look for 'e' or 'E' and an optionally signed integer. 1157 */ 1158 if ( (i < l) && (((c = in.charAt(i) )=='e') || (c == 'E') ) ){ 1159 int expSign = 1; 1160 int expVal = 0; 1161 int reallyBig = Integer.MAX_VALUE / 10; 1162 boolean expOverflow = false; 1163 switch( in.charAt(++i) ){ 1164 case '-': 1165 expSign = -1; 1166 //FALLTHROUGH 1167 case '+': 1168 i++; 1169 } 1170 int expAt = i; 1171 expLoop: 1172 while ( i < l ){ 1173 if ( expVal >= reallyBig ){ 1174 // the next character will cause integer 1175 // overflow. 1176 expOverflow = true; 1177 } 1178 switch ( c = in.charAt(i++) ){ 1179 case '0': 1180 case '1': 1181 case '2': 1182 case '3': 1183 case '4': 1184 case '5': 1185 case '6': 1186 case '7': 1187 case '8': 1188 case '9': 1189 expVal = expVal*10 + ( (int)c - (int)'0' ); 1190 continue; 1191 default: 1192 i--; // back up. 1193 break expLoop; // stop parsing exponent. 1194 } 1195 } 1196 int expLimit = bigDecimalExponent+nDigits+nTrailZero; 1197 if ( expOverflow || ( expVal > expLimit ) ){ 1198 // 1199 // The intent here is to end up with 1200 // infinity or zero, as appropriate. 1201 // The reason for yielding such a small decExponent, 1202 // rather than something intuitive such as 1203 // expSign*Integer.MAX_VALUE, is that this value 1204 // is subject to further manipulation in 1205 // doubleValue() and floatValue(), and I don't want 1206 // it to be able to cause overflow there! 1207 // (The only way we can get into trouble here is for 1208 // really outrageous nDigits+nTrailZero, such as 2 billion. ) 1209 // 1210 decExp = expSign*expLimit; 1211 } else { 1212 // this should not overflow, since we tested 1213 // for expVal > (MAX+N), where N >= abs(decExp) 1214 decExp = decExp + expSign*expVal; 1215 } 1216 1217 // if we saw something not a digit ( or end of string ) 1218 // after the [Ee][+-], without seeing any digits at all 1219 // this is certainly an error. If we saw some digits, 1220 // but then some trailing garbage, that might be ok. 1221 // so we just fall through in that case. 1222 // HUMBUG 1223 if ( i == expAt ) 1224 break parseNumber; // certainly bad 1225 } 1226 /* 1227 * We parsed everything we could. 1228 * If there are leftovers, then this is not good input! 1229 */ 1230 if ( i < l && 1231 ((i != l - 1) || 1232 (in.charAt(i) != 'f' && 1233 in.charAt(i) != 'F' && 1234 in.charAt(i) != 'd' && 1235 in.charAt(i) != 'D'))) { 1236 break parseNumber; // go throw exception 1237 } 1238 1239 return new FloatingDecimal( isNegative, decExp, digits, nDigits, false ); 1240 } catch ( StringIndexOutOfBoundsException e ){ } 1241 throw new NumberFormatException("For input string: \"" + in + "\""); 1242 } 1243 1244 /* 1245 * Take a FloatingDecimal, which we presumably just scanned in, 1246 * and find out what its value is, as a double. 1247 * 1248 * AS A SIDE EFFECT, SET roundDir TO INDICATE PREFERRED 1249 * ROUNDING DIRECTION in case the result is really destined 1250 * for a single-precision float. 1251 */ 1252 1253 public strictfp double doubleValue(){ 1254 int kDigits = Math.min( nDigits, maxDecimalDigits+1 ); 1255 long lValue; 1256 double dValue; 1257 double rValue, tValue; 1258 1259 // First, check for NaN and Infinity values 1260 if(digits == infinity || digits == notANumber) { 1261 if(digits == notANumber) 1262 return Double.NaN; 1263 else 1264 return (isNegative?Double.NEGATIVE_INFINITY:Double.POSITIVE_INFINITY); 1265 } 1266 else { 1267 if (mustSetRoundDir) { 1268 roundDir = 0; 1269 } 1270 /* 1271 * convert the lead kDigits to a long integer. 1272 */ 1273 // (special performance hack: start to do it using int) 1274 int iValue = (int)digits[0]-(int)'0'; 1275 int iDigits = Math.min( kDigits, intDecimalDigits ); 1276 for ( int i=1; i < iDigits; i++ ){ 1277 iValue = iValue*10 + (int)digits[i]-(int)'0'; 1278 } 1279 lValue = (long)iValue; 1280 for ( int i=iDigits; i < kDigits; i++ ){ 1281 lValue = lValue*10L + (long)((int)digits[i]-(int)'0'); 1282 } 1283 dValue = (double)lValue; 1284 int exp = decExponent-kDigits; 1285 /* 1286 * lValue now contains a long integer with the value of 1287 * the first kDigits digits of the number. 1288 * dValue contains the (double) of the same. 1289 */ 1290 1291 if ( nDigits <= maxDecimalDigits ){ 1292 /* 1293 * possibly an easy case. 1294 * We know that the digits can be represented 1295 * exactly. And if the exponent isn't too outrageous, 1296 * the whole thing can be done with one operation, 1297 * thus one rounding error. 1298 * Note that all our constructors trim all leading and 1299 * trailing zeros, so simple values (including zero) 1300 * will always end up here 1301 */ 1302 if (exp == 0 || dValue == 0.0) 1303 return (isNegative)? -dValue : dValue; // small floating integer 1304 else if ( exp >= 0 ){ 1305 if ( exp <= maxSmallTen ){ 1306 /* 1307 * Can get the answer with one operation, 1308 * thus one roundoff. 1309 */ 1310 rValue = dValue * small10pow[exp]; 1311 if ( mustSetRoundDir ){ 1312 tValue = rValue / small10pow[exp]; 1313 roundDir = ( tValue == dValue ) ? 0 1314 :( tValue < dValue ) ? 1 1315 : -1; 1316 } 1317 return (isNegative)? -rValue : rValue; 1318 } 1319 int slop = maxDecimalDigits - kDigits; 1320 if ( exp <= maxSmallTen+slop ){ 1321 /* 1322 * We can multiply dValue by 10^(slop) 1323 * and it is still "small" and exact. 1324 * Then we can multiply by 10^(exp-slop) 1325 * with one rounding. 1326 */ 1327 dValue *= small10pow[slop]; 1328 rValue = dValue * small10pow[exp-slop]; 1329 1330 if ( mustSetRoundDir ){ 1331 tValue = rValue / small10pow[exp-slop]; 1332 roundDir = ( tValue == dValue ) ? 0 1333 :( tValue < dValue ) ? 1 1334 : -1; 1335 } 1336 return (isNegative)? -rValue : rValue; 1337 } 1338 /* 1339 * Else we have a hard case with a positive exp. 1340 */ 1341 } else { 1342 if ( exp >= -maxSmallTen ){ 1343 /* 1344 * Can get the answer in one division. 1345 */ 1346 rValue = dValue / small10pow[-exp]; 1347 tValue = rValue * small10pow[-exp]; 1348 if ( mustSetRoundDir ){ 1349 roundDir = ( tValue == dValue ) ? 0 1350 :( tValue < dValue ) ? 1 1351 : -1; 1352 } 1353 return (isNegative)? -rValue : rValue; 1354 } 1355 /* 1356 * Else we have a hard case with a negative exp. 1357 */ 1358 } 1359 } 1360 1361 /* 1362 * Harder cases: 1363 * The sum of digits plus exponent is greater than 1364 * what we think we can do with one error. 1365 * 1366 * Start by approximating the right answer by, 1367 * naively, scaling by powers of 10. 1368 */ 1369 if ( exp > 0 ){ 1370 if ( decExponent > maxDecimalExponent+1 ){ 1371 /* 1372 * Lets face it. This is going to be 1373 * Infinity. Cut to the chase. 1374 */ 1375 return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY; 1376 } 1377 if ( (exp&15) != 0 ){ 1378 dValue *= small10pow[exp&15]; 1379 } 1380 if ( (exp>>=4) != 0 ){ 1381 int j; 1382 for( j = 0; exp > 1; j++, exp>>=1 ){ 1383 if ( (exp&1)!=0) 1384 dValue *= big10pow[j]; 1385 } 1386 /* 1387 * The reason for the weird exp > 1 condition 1388 * in the above loop was so that the last multiply 1389 * would get unrolled. We handle it here. 1390 * It could overflow. 1391 */ 1392 double t = dValue * big10pow[j]; 1393 if ( Double.isInfinite( t ) ){ 1394 /* 1395 * It did overflow. 1396 * Look more closely at the result. 1397 * If the exponent is just one too large, 1398 * then use the maximum finite as our estimate 1399 * value. Else call the result infinity 1400 * and punt it. 1401 * ( I presume this could happen because 1402 * rounding forces the result here to be 1403 * an ULP or two larger than 1404 * Double.MAX_VALUE ). 1405 */ 1406 t = dValue / 2.0; 1407 t *= big10pow[j]; 1408 if ( Double.isInfinite( t ) ){ 1409 return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY; 1410 } 1411 t = Double.MAX_VALUE; 1412 } 1413 dValue = t; 1414 } 1415 } else if ( exp < 0 ){ 1416 exp = -exp; 1417 if ( decExponent < minDecimalExponent-1 ){ 1418 /* 1419 * Lets face it. This is going to be 1420 * zero. Cut to the chase. 1421 */ 1422 return (isNegative)? -0.0 : 0.0; 1423 } 1424 if ( (exp&15) != 0 ){ 1425 dValue /= small10pow[exp&15]; 1426 } 1427 if ( (exp>>=4) != 0 ){ 1428 int j; 1429 for( j = 0; exp > 1; j++, exp>>=1 ){ 1430 if ( (exp&1)!=0) 1431 dValue *= tiny10pow[j]; 1432 } 1433 /* 1434 * The reason for the weird exp > 1 condition 1435 * in the above loop was so that the last multiply 1436 * would get unrolled. We handle it here. 1437 * It could underflow. 1438 */ 1439 double t = dValue * tiny10pow[j]; 1440 if ( t == 0.0 ){ 1441 /* 1442 * It did underflow. 1443 * Look more closely at the result. 1444 * If the exponent is just one too small, 1445 * then use the minimum finite as our estimate 1446 * value. Else call the result 0.0 1447 * and punt it. 1448 * ( I presume this could happen because 1449 * rounding forces the result here to be 1450 * an ULP or two less than 1451 * Double.MIN_VALUE ). 1452 */ 1453 t = dValue * 2.0; 1454 t *= tiny10pow[j]; 1455 if ( t == 0.0 ){ 1456 return (isNegative)? -0.0 : 0.0; 1457 } 1458 t = Double.MIN_VALUE; 1459 } 1460 dValue = t; 1461 } 1462 } 1463 1464 /* 1465 * dValue is now approximately the result. 1466 * The hard part is adjusting it, by comparison 1467 * with FDBigInt arithmetic. 1468 * Formulate the EXACT big-number result as 1469 * bigD0 * 10^exp 1470 */ 1471 FDBigInt bigD0 = new FDBigInt( lValue, digits, kDigits, nDigits ); 1472 exp = decExponent - nDigits; 1473 1474 correctionLoop: 1475 while(true){ 1476 /* AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES 1477 * bigIntExp and bigIntNBits 1478 */ 1479 FDBigInt bigB = doubleToBigInt( dValue ); 1480 1481 /* 1482 * Scale bigD, bigB appropriately for 1483 * big-integer operations. 1484 * Naively, we multiply by powers of ten 1485 * and powers of two. What we actually do 1486 * is keep track of the powers of 5 and 1487 * powers of 2 we would use, then factor out 1488 * common divisors before doing the work. 1489 */ 1490 int B2, B5; // powers of 2, 5 in bigB 1491 int D2, D5; // powers of 2, 5 in bigD 1492 int Ulp2; // powers of 2 in halfUlp. 1493 if ( exp >= 0 ){ 1494 B2 = B5 = 0; 1495 D2 = D5 = exp; 1496 } else { 1497 B2 = B5 = -exp; 1498 D2 = D5 = 0; 1499 } 1500 if ( bigIntExp >= 0 ){ 1501 B2 += bigIntExp; 1502 } else { 1503 D2 -= bigIntExp; 1504 } 1505 Ulp2 = B2; 1506 // shift bigB and bigD left by a number s. t. 1507 // halfUlp is still an integer. 1508 int hulpbias; 1509 if ( bigIntExp+bigIntNBits <= -expBias+1 ){ 1510 // This is going to be a denormalized number 1511 // (if not actually zero). 1512 // half an ULP is at 2^-(expBias+expShift+1) 1513 hulpbias = bigIntExp+ expBias + expShift; 1514 } else { 1515 hulpbias = expShift + 2 - bigIntNBits; 1516 } 1517 B2 += hulpbias; 1518 D2 += hulpbias; 1519 // if there are common factors of 2, we might just as well 1520 // factor them out, as they add nothing useful. 1521 int common2 = Math.min( B2, Math.min( D2, Ulp2 ) ); 1522 B2 -= common2; 1523 D2 -= common2; 1524 Ulp2 -= common2; 1525 // do multiplications by powers of 5 and 2 1526 bigB = multPow52( bigB, B5, B2 ); 1527 FDBigInt bigD = multPow52( new FDBigInt( bigD0 ), D5, D2 ); 1528 // 1529 // to recap: 1530 // bigB is the scaled-big-int version of our floating-point 1531 // candidate. 1532 // bigD is the scaled-big-int version of the exact value 1533 // as we understand it. 1534 // halfUlp is 1/2 an ulp of bigB, except for special cases 1535 // of exact powers of 2 1536 // 1537 // the plan is to compare bigB with bigD, and if the difference 1538 // is less than halfUlp, then we're satisfied. Otherwise, 1539 // use the ratio of difference to halfUlp to calculate a fudge 1540 // factor to add to the floating value, then go 'round again. 1541 // 1542 FDBigInt diff; 1543 int cmpResult; 1544 boolean overvalue; 1545 if ( (cmpResult = bigB.cmp( bigD ) ) > 0 ){ 1546 overvalue = true; // our candidate is too big. 1547 diff = bigB.sub( bigD ); 1548 if ( (bigIntNBits == 1) && (bigIntExp > -expBias+1) ){ 1549 // candidate is a normalized exact power of 2 and 1550 // is too big. We will be subtracting. 1551 // For our purposes, ulp is the ulp of the 1552 // next smaller range. 1553 Ulp2 -= 1; 1554 if ( Ulp2 < 0 ){ 1555 // rats. Cannot de-scale ulp this far. 1556 // must scale diff in other direction. 1557 Ulp2 = 0; 1558 diff.lshiftMe( 1 ); 1559 } 1560 } 1561 } else if ( cmpResult < 0 ){ 1562 overvalue = false; // our candidate is too small. 1563 diff = bigD.sub( bigB ); 1564 } else { 1565 // the candidate is exactly right! 1566 // this happens with surprising frequency 1567 break correctionLoop; 1568 } 1569 FDBigInt halfUlp = constructPow52( B5, Ulp2 ); 1570 if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){ 1571 // difference is small. 1572 // this is close enough 1573 if (mustSetRoundDir) { 1574 roundDir = overvalue ? -1 : 1; 1575 } 1576 break correctionLoop; 1577 } else if ( cmpResult == 0 ){ 1578 // difference is exactly half an ULP 1579 // round to some other value maybe, then finish 1580 dValue += 0.5*ulp( dValue, overvalue ); 1581 // should check for bigIntNBits == 1 here?? 1582 if (mustSetRoundDir) { 1583 roundDir = overvalue ? -1 : 1; 1584 } 1585 break correctionLoop; 1586 } else { 1587 // difference is non-trivial. 1588 // could scale addend by ratio of difference to 1589 // halfUlp here, if we bothered to compute that difference. 1590 // Most of the time ( I hope ) it is about 1 anyway. 1591 dValue += ulp( dValue, overvalue ); 1592 if ( dValue == 0.0 || dValue == Double.POSITIVE_INFINITY ) 1593 break correctionLoop; // oops. Fell off end of range. 1594 continue; // try again. 1595 } 1596 1597 } 1598 return (isNegative)? -dValue : dValue; 1599 } 1600 } 1601 1602 /* 1603 * Take a FloatingDecimal, which we presumably just scanned in, 1604 * and find out what its value is, as a float. 1605 * This is distinct from doubleValue() to avoid the extremely 1606 * unlikely case of a double rounding error, wherein the conversion 1607 * to double has one rounding error, and the conversion of that double 1608 * to a float has another rounding error, IN THE WRONG DIRECTION, 1609 * ( because of the preference to a zero low-order bit ). 1610 */ 1611 1612 public strictfp float floatValue(){ 1613 int kDigits = Math.min( nDigits, singleMaxDecimalDigits+1 ); 1614 int iValue; 1615 float fValue; 1616 1617 // First, check for NaN and Infinity values 1618 if(digits == infinity || digits == notANumber) { 1619 if(digits == notANumber) 1620 return Float.NaN; 1621 else 1622 return (isNegative?Float.NEGATIVE_INFINITY:Float.POSITIVE_INFINITY); 1623 } 1624 else { 1625 /* 1626 * convert the lead kDigits to an integer. 1627 */ 1628 iValue = (int)digits[0]-(int)'0'; 1629 for ( int i=1; i < kDigits; i++ ){ 1630 iValue = iValue*10 + (int)digits[i]-(int)'0'; 1631 } 1632 fValue = (float)iValue; 1633 int exp = decExponent-kDigits; 1634 /* 1635 * iValue now contains an integer with the value of 1636 * the first kDigits digits of the number. 1637 * fValue contains the (float) of the same. 1638 */ 1639 1640 if ( nDigits <= singleMaxDecimalDigits ){ 1641 /* 1642 * possibly an easy case. 1643 * We know that the digits can be represented 1644 * exactly. And if the exponent isn't too outrageous, 1645 * the whole thing can be done with one operation, 1646 * thus one rounding error. 1647 * Note that all our constructors trim all leading and 1648 * trailing zeros, so simple values (including zero) 1649 * will always end up here. 1650 */ 1651 if (exp == 0 || fValue == 0.0f) 1652 return (isNegative)? -fValue : fValue; // small floating integer 1653 else if ( exp >= 0 ){ 1654 if ( exp <= singleMaxSmallTen ){ 1655 /* 1656 * Can get the answer with one operation, 1657 * thus one roundoff. 1658 */ 1659 fValue *= singleSmall10pow[exp]; 1660 return (isNegative)? -fValue : fValue; 1661 } 1662 int slop = singleMaxDecimalDigits - kDigits; 1663 if ( exp <= singleMaxSmallTen+slop ){ 1664 /* 1665 * We can multiply dValue by 10^(slop) 1666 * and it is still "small" and exact. 1667 * Then we can multiply by 10^(exp-slop) 1668 * with one rounding. 1669 */ 1670 fValue *= singleSmall10pow[slop]; 1671 fValue *= singleSmall10pow[exp-slop]; 1672 return (isNegative)? -fValue : fValue; 1673 } 1674 /* 1675 * Else we have a hard case with a positive exp. 1676 */ 1677 } else { 1678 if ( exp >= -singleMaxSmallTen ){ 1679 /* 1680 * Can get the answer in one division. 1681 */ 1682 fValue /= singleSmall10pow[-exp]; 1683 return (isNegative)? -fValue : fValue; 1684 } 1685 /* 1686 * Else we have a hard case with a negative exp. 1687 */ 1688 } 1689 } else if ( (decExponent >= nDigits) && (nDigits+decExponent <= maxDecimalDigits) ){ 1690 /* 1691 * In double-precision, this is an exact floating integer. 1692 * So we can compute to double, then shorten to float 1693 * with one round, and get the right answer. 1694 * 1695 * First, finish accumulating digits. 1696 * Then convert that integer to a double, multiply 1697 * by the appropriate power of ten, and convert to float. 1698 */ 1699 long lValue = (long)iValue; 1700 for ( int i=kDigits; i < nDigits; i++ ){ 1701 lValue = lValue*10L + (long)((int)digits[i]-(int)'0'); 1702 } 1703 double dValue = (double)lValue; 1704 exp = decExponent-nDigits; 1705 dValue *= small10pow[exp]; 1706 fValue = (float)dValue; 1707 return (isNegative)? -fValue : fValue; 1708 1709 } 1710 /* 1711 * Harder cases: 1712 * The sum of digits plus exponent is greater than 1713 * what we think we can do with one error. 1714 * 1715 * Start by weeding out obviously out-of-range 1716 * results, then convert to double and go to 1717 * common hard-case code. 1718 */ 1719 if ( decExponent > singleMaxDecimalExponent+1 ){ 1720 /* 1721 * Lets face it. This is going to be 1722 * Infinity. Cut to the chase. 1723 */ 1724 return (isNegative)? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY; 1725 } else if ( decExponent < singleMinDecimalExponent-1 ){ 1726 /* 1727 * Lets face it. This is going to be 1728 * zero. Cut to the chase. 1729 */ 1730 return (isNegative)? -0.0f : 0.0f; 1731 } 1732 1733 /* 1734 * Here, we do 'way too much work, but throwing away 1735 * our partial results, and going and doing the whole 1736 * thing as double, then throwing away half the bits that computes 1737 * when we convert back to float. 1738 * 1739 * The alternative is to reproduce the whole multiple-precision 1740 * algorithm for float precision, or to try to parameterize it 1741 * for common usage. The former will take about 400 lines of code, 1742 * and the latter I tried without success. Thus the semi-hack 1743 * answer here. 1744 */ 1745 mustSetRoundDir = !fromHex; 1746 double dValue = doubleValue(); 1747 return stickyRound( dValue ); 1748 } 1749 } 1750 1751 1752 /* 1753 * All the positive powers of 10 that can be 1754 * represented exactly in double/float. 1755 */ 1756 private static final double small10pow[] = { 1757 1.0e0, 1758 1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5, 1759 1.0e6, 1.0e7, 1.0e8, 1.0e9, 1.0e10, 1760 1.0e11, 1.0e12, 1.0e13, 1.0e14, 1.0e15, 1761 1.0e16, 1.0e17, 1.0e18, 1.0e19, 1.0e20, 1762 1.0e21, 1.0e22 1763 }; 1764 1765 private static final float singleSmall10pow[] = { 1766 1.0e0f, 1767 1.0e1f, 1.0e2f, 1.0e3f, 1.0e4f, 1.0e5f, 1768 1.0e6f, 1.0e7f, 1.0e8f, 1.0e9f, 1.0e10f 1769 }; 1770 1771 private static final double big10pow[] = { 1772 1e16, 1e32, 1e64, 1e128, 1e256 }; 1773 private static final double tiny10pow[] = { 1774 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 }; 1775 1776 private static final int maxSmallTen = small10pow.length-1; 1777 private static final int singleMaxSmallTen = singleSmall10pow.length-1; 1778 1779 private static final int small5pow[] = { 1780 1, 1781 5, 1782 5*5, 1783 5*5*5, 1784 5*5*5*5, 1785 5*5*5*5*5, 1786 5*5*5*5*5*5, 1787 5*5*5*5*5*5*5, 1788 5*5*5*5*5*5*5*5, 1789 5*5*5*5*5*5*5*5*5, 1790 5*5*5*5*5*5*5*5*5*5, 1791 5*5*5*5*5*5*5*5*5*5*5, 1792 5*5*5*5*5*5*5*5*5*5*5*5, 1793 5*5*5*5*5*5*5*5*5*5*5*5*5 1794 }; 1795 1796 1797 private static final long long5pow[] = { 1798 1L, 1799 5L, 1800 5L*5, 1801 5L*5*5, 1802 5L*5*5*5, 1803 5L*5*5*5*5, 1804 5L*5*5*5*5*5, 1805 5L*5*5*5*5*5*5, 1806 5L*5*5*5*5*5*5*5, 1807 5L*5*5*5*5*5*5*5*5, 1808 5L*5*5*5*5*5*5*5*5*5, 1809 5L*5*5*5*5*5*5*5*5*5*5, 1810 5L*5*5*5*5*5*5*5*5*5*5*5, 1811 5L*5*5*5*5*5*5*5*5*5*5*5*5, 1812 5L*5*5*5*5*5*5*5*5*5*5*5*5*5, 1813 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1814 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1815 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1816 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1817 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1818 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1819 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1820 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1821 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1822 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5, 1823 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, 1824 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, 1825 }; 1826 1827 // approximately ceil( log2( long5pow[i] ) ) 1828 private static final int n5bits[] = { 1829 0, 1830 3, 1831 5, 1832 7, 1833 10, 1834 12, 1835 14, 1836 17, 1837 19, 1838 21, 1839 24, 1840 26, 1841 28, 1842 31, 1843 33, 1844 35, 1845 38, 1846 40, 1847 42, 1848 45, 1849 47, 1850 49, 1851 52, 1852 54, 1853 56, 1854 59, 1855 61, 1856 }; 1857 1858 private static final char infinity[] = { 'I', 'n', 'f', 'i', 'n', 'i', 't', 'y' }; 1859 private static final char notANumber[] = { 'N', 'a', 'N' }; 1860 private static final char zero[] = { '0', '0', '0', '0', '0', '0', '0', '0' }; 1861 1862 1863 /* 1864 * Grammar is compatible with hexadecimal floating-point constants 1865 * described in section 6.4.4.2 of the C99 specification. 1866 */ 1867 private static Pattern hexFloatPattern = null; 1868 private static synchronized Pattern getHexFloatPattern() { 1869 if (hexFloatPattern == null) { 1870 hexFloatPattern = Pattern.compile( 1871 //1 234 56 7 8 9 1872 "([-+])?0[xX](((\\p{XDigit}+)\\.?)|((\\p{XDigit}*)\\.(\\p{XDigit}+)))[pP]([-+])?(\\p{Digit}+)[fFdD]?" 1873 ); 1874 } 1875 return hexFloatPattern; 1876 } 1877 1878 /* 1879 * Convert string s to a suitable floating decimal; uses the 1880 * double constructor and set the roundDir variable appropriately 1881 * in case the value is later converted to a float. 1882 */ 1883 static FloatingDecimal parseHexString(String s) { 1884 // Verify string is a member of the hexadecimal floating-point 1885 // string language. 1886 Matcher m = getHexFloatPattern().matcher(s); 1887 boolean validInput = m.matches(); 1888 1889 if (!validInput) { 1890 // Input does not match pattern 1891 throw new NumberFormatException("For input string: \"" + s + "\""); 1892 } else { // validInput 1893 /* 1894 * We must isolate the sign, significand, and exponent 1895 * fields. The sign value is straightforward. Since 1896 * floating-point numbers are stored with a normalized 1897 * representation, the significand and exponent are 1898 * interrelated. 1899 * 1900 * After extracting the sign, we normalized the 1901 * significand as a hexadecimal value, calculating an 1902 * exponent adjust for any shifts made during 1903 * normalization. If the significand is zero, the 1904 * exponent doesn't need to be examined since the output 1905 * will be zero. 1906 * 1907 * Next the exponent in the input string is extracted. 1908 * Afterwards, the significand is normalized as a *binary* 1909 * value and the input value's normalized exponent can be 1910 * computed. The significand bits are copied into a 1911 * double significand; if the string has more logical bits 1912 * than can fit in a double, the extra bits affect the 1913 * round and sticky bits which are used to round the final 1914 * value. 1915 */ 1916 1917 // Extract significand sign 1918 String group1 = m.group(1); 1919 double sign = (( group1 == null ) || group1.equals("+"))? 1.0 : -1.0; 1920 1921 1922 // Extract Significand magnitude 1923 /* 1924 * Based on the form of the significand, calculate how the 1925 * binary exponent needs to be adjusted to create a 1926 * normalized *hexadecimal* floating-point number; that 1927 * is, a number where there is one nonzero hex digit to 1928 * the left of the (hexa)decimal point. Since we are 1929 * adjusting a binary, not hexadecimal exponent, the 1930 * exponent is adjusted by a multiple of 4. 1931 * 1932 * There are a number of significand scenarios to consider; 1933 * letters are used in indicate nonzero digits: 1934 * 1935 * 1. 000xxxx => x.xxx normalized 1936 * increase exponent by (number of x's - 1)*4 1937 * 1938 * 2. 000xxx.yyyy => x.xxyyyy normalized 1939 * increase exponent by (number of x's - 1)*4 1940 * 1941 * 3. .000yyy => y.yy normalized 1942 * decrease exponent by (number of zeros + 1)*4 1943 * 1944 * 4. 000.00000yyy => y.yy normalized 1945 * decrease exponent by (number of zeros to right of point + 1)*4 1946 * 1947 * If the significand is exactly zero, return a properly 1948 * signed zero. 1949 */ 1950 1951 String significandString =null; 1952 int signifLength = 0; 1953 int exponentAdjust = 0; 1954 { 1955 int leftDigits = 0; // number of meaningful digits to 1956 // left of "decimal" point 1957 // (leading zeros stripped) 1958 int rightDigits = 0; // number of digits to right of 1959 // "decimal" point; leading zeros 1960 // must always be accounted for 1961 /* 1962 * The significand is made up of either 1963 * 1964 * 1. group 4 entirely (integer portion only) 1965 * 1966 * OR 1967 * 1968 * 2. the fractional portion from group 7 plus any 1969 * (optional) integer portions from group 6. 1970 */ 1971 String group4; 1972 if( (group4 = m.group(4)) != null) { // Integer-only significand 1973 // Leading zeros never matter on the integer portion 1974 significandString = stripLeadingZeros(group4); 1975 leftDigits = significandString.length(); 1976 } 1977 else { 1978 // Group 6 is the optional integer; leading zeros 1979 // never matter on the integer portion 1980 String group6 = stripLeadingZeros(m.group(6)); 1981 leftDigits = group6.length(); 1982 1983 // fraction 1984 String group7 = m.group(7); 1985 rightDigits = group7.length(); 1986 1987 // Turn "integer.fraction" into "integer"+"fraction" 1988 significandString = 1989 ((group6 == null)?"":group6) + // is the null 1990 // check necessary? 1991 group7; 1992 } 1993 1994 significandString = stripLeadingZeros(significandString); 1995 signifLength = significandString.length(); 1996 1997 /* 1998 * Adjust exponent as described above 1999 */ 2000 if (leftDigits >= 1) { // Cases 1 and 2 2001 exponentAdjust = 4*(leftDigits - 1); 2002 } else { // Cases 3 and 4 2003 exponentAdjust = -4*( rightDigits - signifLength + 1); 2004 } 2005 2006 // If the significand is zero, the exponent doesn't 2007 // matter; return a properly signed zero. 2008 2009 if (signifLength == 0) { // Only zeros in input 2010 return new FloatingDecimal(sign * 0.0); 2011 } 2012 } 2013 2014 // Extract Exponent 2015 /* 2016 * Use an int to read in the exponent value; this should 2017 * provide more than sufficient range for non-contrived 2018 * inputs. If reading the exponent in as an int does 2019 * overflow, examine the sign of the exponent and 2020 * significand to determine what to do. 2021 */ 2022 String group8 = m.group(8); 2023 boolean positiveExponent = ( group8 == null ) || group8.equals("+"); 2024 long unsignedRawExponent; 2025 try { 2026 unsignedRawExponent = Integer.parseInt(m.group(9)); 2027 } 2028 catch (NumberFormatException e) { 2029 // At this point, we know the exponent is 2030 // syntactically well-formed as a sequence of 2031 // digits. Therefore, if an NumberFormatException 2032 // is thrown, it must be due to overflowing int's 2033 // range. Also, at this point, we have already 2034 // checked for a zero significand. Thus the signs 2035 // of the exponent and significand determine the 2036 // final result: 2037 // 2038 // significand 2039 // + - 2040 // exponent + +infinity -infinity 2041 // - +0.0 -0.0 2042 return new FloatingDecimal(sign * (positiveExponent ? 2043 Double.POSITIVE_INFINITY : 0.0)); 2044 } 2045 2046 long rawExponent = 2047 (positiveExponent ? 1L : -1L) * // exponent sign 2048 unsignedRawExponent; // exponent magnitude 2049 2050 // Calculate partially adjusted exponent 2051 long exponent = rawExponent + exponentAdjust ; 2052 2053 // Starting copying non-zero bits into proper position in 2054 // a long; copy explicit bit too; this will be masked 2055 // later for normal values. 2056 2057 boolean round = false; 2058 boolean sticky = false; 2059 int bitsCopied=0; 2060 int nextShift=0; 2061 long significand=0L; 2062 // First iteration is different, since we only copy 2063 // from the leading significand bit; one more exponent 2064 // adjust will be needed... 2065 2066 // IMPORTANT: make leadingDigit a long to avoid 2067 // surprising shift semantics! 2068 long leadingDigit = getHexDigit(significandString, 0); 2069 2070 /* 2071 * Left shift the leading digit (53 - (bit position of 2072 * leading 1 in digit)); this sets the top bit of the 2073 * significand to 1. The nextShift value is adjusted 2074 * to take into account the number of bit positions of 2075 * the leadingDigit actually used. Finally, the 2076 * exponent is adjusted to normalize the significand 2077 * as a binary value, not just a hex value. 2078 */ 2079 if (leadingDigit == 1) { 2080 significand |= leadingDigit << 52; 2081 nextShift = 52 - 4; 2082 /* exponent += 0 */ } 2083 else if (leadingDigit <= 3) { // [2, 3] 2084 significand |= leadingDigit << 51; 2085 nextShift = 52 - 5; 2086 exponent += 1; 2087 } 2088 else if (leadingDigit <= 7) { // [4, 7] 2089 significand |= leadingDigit << 50; 2090 nextShift = 52 - 6; 2091 exponent += 2; 2092 } 2093 else if (leadingDigit <= 15) { // [8, f] 2094 significand |= leadingDigit << 49; 2095 nextShift = 52 - 7; 2096 exponent += 3; 2097 } else { 2098 throw new AssertionError("Result from digit conversion too large!"); 2099 } 2100 // The preceding if-else could be replaced by a single 2101 // code block based on the high-order bit set in 2102 // leadingDigit. Given leadingOnePosition, 2103 2104 // significand |= leadingDigit << (SIGNIFICAND_WIDTH - leadingOnePosition); 2105 // nextShift = 52 - (3 + leadingOnePosition); 2106 // exponent += (leadingOnePosition-1); 2107 2108 2109 /* 2110 * Now the exponent variable is equal to the normalized 2111 * binary exponent. Code below will make representation 2112 * adjustments if the exponent is incremented after 2113 * rounding (includes overflows to infinity) or if the 2114 * result is subnormal. 2115 */ 2116 2117 // Copy digit into significand until the significand can't 2118 // hold another full hex digit or there are no more input 2119 // hex digits. 2120 int i = 0; 2121 for(i = 1; 2122 i < signifLength && nextShift >= 0; 2123 i++) { 2124 long currentDigit = getHexDigit(significandString, i); 2125 significand |= (currentDigit << nextShift); 2126 nextShift-=4; 2127 } 2128 2129 // After the above loop, the bulk of the string is copied. 2130 // Now, we must copy any partial hex digits into the 2131 // significand AND compute the round bit and start computing 2132 // sticky bit. 2133 2134 if ( i < signifLength ) { // at least one hex input digit exists 2135 long currentDigit = getHexDigit(significandString, i); 2136 2137 // from nextShift, figure out how many bits need 2138 // to be copied, if any 2139 switch(nextShift) { // must be negative 2140 case -1: 2141 // three bits need to be copied in; can 2142 // set round bit 2143 significand |= ((currentDigit & 0xEL) >> 1); 2144 round = (currentDigit & 0x1L) != 0L; 2145 break; 2146 2147 case -2: 2148 // two bits need to be copied in; can 2149 // set round and start sticky 2150 significand |= ((currentDigit & 0xCL) >> 2); 2151 round = (currentDigit &0x2L) != 0L; 2152 sticky = (currentDigit & 0x1L) != 0; 2153 break; 2154 2155 case -3: 2156 // one bit needs to be copied in 2157 significand |= ((currentDigit & 0x8L)>>3); 2158 // Now set round and start sticky, if possible 2159 round = (currentDigit &0x4L) != 0L; 2160 sticky = (currentDigit & 0x3L) != 0; 2161 break; 2162 2163 case -4: 2164 // all bits copied into significand; set 2165 // round and start sticky 2166 round = ((currentDigit & 0x8L) != 0); // is top bit set? 2167 // nonzeros in three low order bits? 2168 sticky = (currentDigit & 0x7L) != 0; 2169 break; 2170 2171 default: 2172 throw new AssertionError("Unexpected shift distance remainder."); 2173 // break; 2174 } 2175 2176 // Round is set; sticky might be set. 2177 2178 // For the sticky bit, it suffices to check the 2179 // current digit and test for any nonzero digits in 2180 // the remaining unprocessed input. 2181 i++; 2182 while(i < signifLength && !sticky) { 2183 currentDigit = getHexDigit(significandString,i); 2184 sticky = sticky || (currentDigit != 0); 2185 i++; 2186 } 2187 2188 } 2189 // else all of string was seen, round and sticky are 2190 // correct as false. 2191 2192 2193 // Check for overflow and update exponent accordingly. 2194 2195 if (exponent > DoubleConsts.MAX_EXPONENT) { // Infinite result 2196 // overflow to properly signed infinity 2197 return new FloatingDecimal(sign * Double.POSITIVE_INFINITY); 2198 } else { // Finite return value 2199 if (exponent <= DoubleConsts.MAX_EXPONENT && // (Usually) normal result 2200 exponent >= DoubleConsts.MIN_EXPONENT) { 2201 2202 // The result returned in this block cannot be a 2203 // zero or subnormal; however after the 2204 // significand is adjusted from rounding, we could 2205 // still overflow in infinity. 2206 2207 // AND exponent bits into significand; if the 2208 // significand is incremented and overflows from 2209 // rounding, this combination will update the 2210 // exponent correctly, even in the case of 2211 // Double.MAX_VALUE overflowing to infinity. 2212 2213 significand = (( ((long)exponent + 2214 (long)DoubleConsts.EXP_BIAS) << 2215 (DoubleConsts.SIGNIFICAND_WIDTH-1)) 2216 & DoubleConsts.EXP_BIT_MASK) | 2217 (DoubleConsts.SIGNIF_BIT_MASK & significand); 2218 2219 } else { // Subnormal or zero 2220 // (exponent < DoubleConsts.MIN_EXPONENT) 2221 2222 if (exponent < (DoubleConsts.MIN_SUB_EXPONENT -1 )) { 2223 // No way to round back to nonzero value 2224 // regardless of significand if the exponent is 2225 // less than -1075. 2226 return new FloatingDecimal(sign * 0.0); 2227 } else { // -1075 <= exponent <= MIN_EXPONENT -1 = -1023 2228 /* 2229 * Find bit position to round to; recompute 2230 * round and sticky bits, and shift 2231 * significand right appropriately. 2232 */ 2233 2234 sticky = sticky || round; 2235 round = false; 2236 2237 // Number of bits of significand to preserve is 2238 // exponent - abs_min_exp +1 2239 // check: 2240 // -1075 +1074 + 1 = 0 2241 // -1023 +1074 + 1 = 52 2242 2243 int bitsDiscarded = 53 - 2244 ((int)exponent - DoubleConsts.MIN_SUB_EXPONENT + 1); 2245 assert bitsDiscarded >= 1 && bitsDiscarded <= 53; 2246 2247 // What to do here: 2248 // First, isolate the new round bit 2249 round = (significand & (1L << (bitsDiscarded -1))) != 0L; 2250 if (bitsDiscarded > 1) { 2251 // create mask to update sticky bits; low 2252 // order bitsDiscarded bits should be 1 2253 long mask = ~((~0L) << (bitsDiscarded -1)); 2254 sticky = sticky || ((significand & mask) != 0L ) ; 2255 } 2256 2257 // Now, discard the bits 2258 significand = significand >> bitsDiscarded; 2259 2260 significand = (( ((long)(DoubleConsts.MIN_EXPONENT -1) + // subnorm exp. 2261 (long)DoubleConsts.EXP_BIAS) << 2262 (DoubleConsts.SIGNIFICAND_WIDTH-1)) 2263 & DoubleConsts.EXP_BIT_MASK) | 2264 (DoubleConsts.SIGNIF_BIT_MASK & significand); 2265 } 2266 } 2267 2268 // The significand variable now contains the currently 2269 // appropriate exponent bits too. 2270 2271 /* 2272 * Determine if significand should be incremented; 2273 * making this determination depends on the least 2274 * significant bit and the round and sticky bits. 2275 * 2276 * Round to nearest even rounding table, adapted from 2277 * table 4.7 in "Computer Arithmetic" by IsraelKoren. 2278 * The digit to the left of the "decimal" point is the 2279 * least significant bit, the digits to the right of 2280 * the point are the round and sticky bits 2281 * 2282 * Number Round(x) 2283 * x0.00 x0. 2284 * x0.01 x0. 2285 * x0.10 x0. 2286 * x0.11 x1. = x0. +1 2287 * x1.00 x1. 2288 * x1.01 x1. 2289 * x1.10 x1. + 1 2290 * x1.11 x1. + 1 2291 */ 2292 boolean incremented = false; 2293 boolean leastZero = ((significand & 1L) == 0L); 2294 if( ( leastZero && round && sticky ) || 2295 ((!leastZero) && round )) { 2296 incremented = true; 2297 significand++; 2298 } 2299 2300 FloatingDecimal fd = new FloatingDecimal(FpUtils.rawCopySign( 2301 Double.longBitsToDouble(significand), 2302 sign)); 2303 2304 /* 2305 * Set roundingDir variable field of fd properly so 2306 * that the input string can be properly rounded to a 2307 * float value. There are two cases to consider: 2308 * 2309 * 1. rounding to double discards sticky bit 2310 * information that would change the result of a float 2311 * rounding (near halfway case between two floats) 2312 * 2313 * 2. rounding to double rounds up when rounding up 2314 * would not occur when rounding to float. 2315 * 2316 * For former case only needs to be considered when 2317 * the bits rounded away when casting to float are all 2318 * zero; otherwise, float round bit is properly set 2319 * and sticky will already be true. 2320 * 2321 * The lower exponent bound for the code below is the 2322 * minimum (normalized) subnormal exponent - 1 since a 2323 * value with that exponent can round up to the 2324 * minimum subnormal value and the sticky bit 2325 * information must be preserved (i.e. case 1). 2326 */ 2327 if ((exponent >= FloatConsts.MIN_SUB_EXPONENT-1) && 2328 (exponent <= FloatConsts.MAX_EXPONENT ) ){ 2329 // Outside above exponent range, the float value 2330 // will be zero or infinity. 2331 2332 /* 2333 * If the low-order 28 bits of a rounded double 2334 * significand are 0, the double could be a 2335 * half-way case for a rounding to float. If the 2336 * double value is a half-way case, the double 2337 * significand may have to be modified to round 2338 * the the right float value (see the stickyRound 2339 * method). If the rounding to double has lost 2340 * what would be float sticky bit information, the 2341 * double significand must be incremented. If the 2342 * double value's significand was itself 2343 * incremented, the float value may end up too 2344 * large so the increment should be undone. 2345 */ 2346 if ((significand & 0xfffffffL) == 0x0L) { 2347 // For negative values, the sign of the 2348 // roundDir is the same as for positive values 2349 // since adding 1 increasing the significand's 2350 // magnitude and subtracting 1 decreases the 2351 // significand's magnitude. If neither round 2352 // nor sticky is true, the double value is 2353 // exact and no adjustment is required for a 2354 // proper float rounding. 2355 if( round || sticky) { 2356 if (leastZero) { // prerounding lsb is 0 2357 // If round and sticky were both true, 2358 // and the least significant 2359 // significand bit were 0, the rounded 2360 // significand would not have its 2361 // low-order bits be zero. Therefore, 2362 // we only need to adjust the 2363 // significand if round XOR sticky is 2364 // true. 2365 if (round ^ sticky) { 2366 fd.roundDir = 1; 2367 } 2368 } 2369 else { // prerounding lsb is 1 2370 // If the prerounding lsb is 1 and the 2371 // resulting significand has its 2372 // low-order bits zero, the significand 2373 // was incremented. Here, we undo the 2374 // increment, which will ensure the 2375 // right guard and sticky bits for the 2376 // float rounding. 2377 if (round) 2378 fd.roundDir = -1; 2379 } 2380 } 2381 } 2382 } 2383 2384 fd.fromHex = true; 2385 return fd; 2386 } 2387 } 2388 } 2389 2390 /** 2391 * Return <code>s</code> with any leading zeros removed. 2392 */ 2393 static String stripLeadingZeros(String s) { 2394 return s.replaceFirst("^0+", ""); 2395 } 2396 2397 /** 2398 * Extract a hexadecimal digit from position <code>position</code> 2399 * of string <code>s</code>. 2400 */ 2401 static int getHexDigit(String s, int position) { 2402 int value = Character.digit(s.charAt(position), 16); 2403 if (value <= -1 || value >= 16) { 2404 throw new AssertionError("Unexpected failure of digit conversion of " + 2405 s.charAt(position)); 2406 } 2407 return value; 2408 } 2409 2410 2411 } 2412 2413 /* 2414 * A really, really simple bigint package 2415 * tailored to the needs of floating base conversion. 2416 */ 2417 class FDBigInt { 2418 int nWords; // number of words used 2419 int data[]; // value: data[0] is least significant 2420 2421 2422 public FDBigInt( int v ){ 2423 nWords = 1; 2424 data = new int[1]; 2425 data[0] = v; 2426 } 2427 2428 public FDBigInt( long v ){ 2429 data = new int[2]; 2430 data[0] = (int)v; 2431 data[1] = (int)(v>>>32); 2432 nWords = (data[1]==0) ? 1 : 2; 2433 } 2434 2435 public FDBigInt( FDBigInt other ){ 2436 data = new int[nWords = other.nWords]; 2437 System.arraycopy( other.data, 0, data, 0, nWords ); 2438 } 2439 2440 private FDBigInt( int [] d, int n ){ 2441 data = d; 2442 nWords = n; 2443 } 2444 2445 public FDBigInt( long seed, char digit[], int nd0, int nd ){ 2446 int n= (nd+8)/9; // estimate size needed. 2447 if ( n < 2 ) n = 2; 2448 data = new int[n]; // allocate enough space 2449 data[0] = (int)seed; // starting value 2450 data[1] = (int)(seed>>>32); 2451 nWords = (data[1]==0) ? 1 : 2; 2452 int i = nd0; 2453 int limit = nd-5; // slurp digits 5 at a time. 2454 int v; 2455 while ( i < limit ){ 2456 int ilim = i+5; 2457 v = (int)digit[i++]-(int)'0'; 2458 while( i <ilim ){ 2459 v = 10*v + (int)digit[i++]-(int)'0'; 2460 } 2461 multaddMe( 100000, v); // ... where 100000 is 10^5. 2462 } 2463 int factor = 1; 2464 v = 0; 2465 while ( i < nd ){ 2466 v = 10*v + (int)digit[i++]-(int)'0'; 2467 factor *= 10; 2468 } 2469 if ( factor != 1 ){ 2470 multaddMe( factor, v ); 2471 } 2472 } 2473 2474 /* 2475 * Left shift by c bits. 2476 * Shifts this in place. 2477 */ 2478 public void 2479 lshiftMe( int c )throws IllegalArgumentException { 2480 if ( c <= 0 ){ 2481 if ( c == 0 ) 2482 return; // silly. 2483 else 2484 throw new IllegalArgumentException("negative shift count"); 2485 } 2486 int wordcount = c>>5; 2487 int bitcount = c & 0x1f; 2488 int anticount = 32-bitcount; 2489 int t[] = data; 2490 int s[] = data; 2491 if ( nWords+wordcount+1 > t.length ){ 2492 // reallocate. 2493 t = new int[ nWords+wordcount+1 ]; 2494 } 2495 int target = nWords+wordcount; 2496 int src = nWords-1; 2497 if ( bitcount == 0 ){ 2498 // special hack, since an anticount of 32 won't go! 2499 System.arraycopy( s, 0, t, wordcount, nWords ); 2500 target = wordcount-1; 2501 } else { 2502 t[target--] = s[src]>>>anticount; 2503 while ( src >= 1 ){ 2504 t[target--] = (s[src]<<bitcount) | (s[--src]>>>anticount); 2505 } 2506 t[target--] = s[src]<<bitcount; 2507 } 2508 while( target >= 0 ){ 2509 t[target--] = 0; 2510 } 2511 data = t; 2512 nWords += wordcount + 1; 2513 // may have constructed high-order word of 0. 2514 // if so, trim it 2515 while ( nWords > 1 && data[nWords-1] == 0 ) 2516 nWords--; 2517 } 2518 2519 /* 2520 * normalize this number by shifting until 2521 * the MSB of the number is at 0x08000000. 2522 * This is in preparation for quoRemIteration, below. 2523 * The idea is that, to make division easier, we want the 2524 * divisor to be "normalized" -- usually this means shifting 2525 * the MSB into the high words sign bit. But because we know that 2526 * the quotient will be 0 < q < 10, we would like to arrange that 2527 * the dividend not span up into another word of precision. 2528 * (This needs to be explained more clearly!) 2529 */ 2530 public int 2531 normalizeMe() throws IllegalArgumentException { 2532 int src; 2533 int wordcount = 0; 2534 int bitcount = 0; 2535 int v = 0; 2536 for ( src= nWords-1 ; src >= 0 && (v=data[src]) == 0 ; src--){ 2537 wordcount += 1; 2538 } 2539 if ( src < 0 ){ 2540 // oops. Value is zero. Cannot normalize it! 2541 throw new IllegalArgumentException("zero value"); 2542 } 2543 /* 2544 * In most cases, we assume that wordcount is zero. This only 2545 * makes sense, as we try not to maintain any high-order 2546 * words full of zeros. In fact, if there are zeros, we will 2547 * simply SHORTEN our number at this point. Watch closely... 2548 */ 2549 nWords -= wordcount; 2550 /* 2551 * Compute how far left we have to shift v s.t. its highest- 2552 * order bit is in the right place. Then call lshiftMe to 2553 * do the work. 2554 */ 2555 if ( (v & 0xf0000000) != 0 ){ 2556 // will have to shift up into the next word. 2557 // too bad. 2558 for( bitcount = 32 ; (v & 0xf0000000) != 0 ; bitcount-- ) 2559 v >>>= 1; 2560 } else { 2561 while ( v <= 0x000fffff ){ 2562 // hack: byte-at-a-time shifting 2563 v <<= 8; 2564 bitcount += 8; 2565 } 2566 while ( v <= 0x07ffffff ){ 2567 v <<= 1; 2568 bitcount += 1; 2569 } 2570 } 2571 if ( bitcount != 0 ) 2572 lshiftMe( bitcount ); 2573 return bitcount; 2574 } 2575 2576 /* 2577 * Multiply a FDBigInt by an int. 2578 * Result is a new FDBigInt. 2579 */ 2580 public FDBigInt 2581 mult( int iv ) { 2582 long v = iv; 2583 int r[]; 2584 long p; 2585 2586 // guess adequate size of r. 2587 r = new int[ ( v * ((long)data[nWords-1]&0xffffffffL) > 0xfffffffL ) ? nWords+1 : nWords ]; 2588 p = 0L; 2589 for( int i=0; i < nWords; i++ ) { 2590 p += v * ((long)data[i]&0xffffffffL); 2591 r[i] = (int)p; 2592 p >>>= 32; 2593 } 2594 if ( p == 0L){ 2595 return new FDBigInt( r, nWords ); 2596 } else { 2597 r[nWords] = (int)p; 2598 return new FDBigInt( r, nWords+1 ); 2599 } 2600 } 2601 2602 /* 2603 * Multiply a FDBigInt by an int and add another int. 2604 * Result is computed in place. 2605 * Hope it fits! 2606 */ 2607 public void 2608 multaddMe( int iv, int addend ) { 2609 long v = iv; 2610 long p; 2611 2612 // unroll 0th iteration, doing addition. 2613 p = v * ((long)data[0]&0xffffffffL) + ((long)addend&0xffffffffL); 2614 data[0] = (int)p; 2615 p >>>= 32; 2616 for( int i=1; i < nWords; i++ ) { 2617 p += v * ((long)data[i]&0xffffffffL); 2618 data[i] = (int)p; 2619 p >>>= 32; 2620 } 2621 if ( p != 0L){ 2622 data[nWords] = (int)p; // will fail noisily if illegal! 2623 nWords++; 2624 } 2625 } 2626 2627 /* 2628 * Multiply a FDBigInt by another FDBigInt. 2629 * Result is a new FDBigInt. 2630 */ 2631 public FDBigInt 2632 mult( FDBigInt other ){ 2633 // crudely guess adequate size for r 2634 int r[] = new int[ nWords + other.nWords ]; 2635 int i; 2636 // I think I am promised zeros... 2637 2638 for( i = 0; i < this.nWords; i++ ){ 2639 long v = (long)this.data[i] & 0xffffffffL; // UNSIGNED CONVERSION 2640 long p = 0L; 2641 int j; 2642 for( j = 0; j < other.nWords; j++ ){ 2643 p += ((long)r[i+j]&0xffffffffL) + v*((long)other.data[j]&0xffffffffL); // UNSIGNED CONVERSIONS ALL 'ROUND. 2644 r[i+j] = (int)p; 2645 p >>>= 32; 2646 } 2647 r[i+j] = (int)p; 2648 } 2649 // compute how much of r we actually needed for all that. 2650 for ( i = r.length-1; i> 0; i--) 2651 if ( r[i] != 0 ) 2652 break; 2653 return new FDBigInt( r, i+1 ); 2654 } 2655 2656 /* 2657 * Add one FDBigInt to another. Return a FDBigInt 2658 */ 2659 public FDBigInt 2660 add( FDBigInt other ){ 2661 int i; 2662 int a[], b[]; 2663 int n, m; 2664 long c = 0L; 2665 // arrange such that a.nWords >= b.nWords; 2666 // n = a.nWords, m = b.nWords 2667 if ( this.nWords >= other.nWords ){ 2668 a = this.data; 2669 n = this.nWords; 2670 b = other.data; 2671 m = other.nWords; 2672 } else { 2673 a = other.data; 2674 n = other.nWords; 2675 b = this.data; 2676 m = this.nWords; 2677 } 2678 int r[] = new int[ n ]; 2679 for ( i = 0; i < n; i++ ){ 2680 c += (long)a[i] & 0xffffffffL; 2681 if ( i < m ){ 2682 c += (long)b[i] & 0xffffffffL; 2683 } 2684 r[i] = (int) c; 2685 c >>= 32; // signed shift. 2686 } 2687 if ( c != 0L ){ 2688 // oops -- carry out -- need longer result. 2689 int s[] = new int[ r.length+1 ]; 2690 System.arraycopy( r, 0, s, 0, r.length ); 2691 s[i++] = (int)c; 2692 return new FDBigInt( s, i ); 2693 } 2694 return new FDBigInt( r, i ); 2695 } 2696 2697 /* 2698 * Subtract one FDBigInt from another. Return a FDBigInt 2699 * Assert that the result is positive. 2700 */ 2701 public FDBigInt 2702 sub( FDBigInt other ){ 2703 int r[] = new int[ this.nWords ]; 2704 int i; 2705 int n = this.nWords; 2706 int m = other.nWords; 2707 int nzeros = 0; 2708 long c = 0L; 2709 for ( i = 0; i < n; i++ ){ 2710 c += (long)this.data[i] & 0xffffffffL; 2711 if ( i < m ){ 2712 c -= (long)other.data[i] & 0xffffffffL; 2713 } 2714 if ( ( r[i] = (int) c ) == 0 ) 2715 nzeros++; 2716 else 2717 nzeros = 0; 2718 c >>= 32; // signed shift 2719 } 2720 assert c == 0L : c; // borrow out of subtract 2721 assert dataInRangeIsZero(i, m, other); // negative result of subtract 2722 return new FDBigInt( r, n-nzeros ); 2723 } 2724 2725 private static boolean dataInRangeIsZero(int i, int m, FDBigInt other) { 2726 while ( i < m ) 2727 if (other.data[i++] != 0) 2728 return false; 2729 return true; 2730 } 2731 2732 /* 2733 * Compare FDBigInt with another FDBigInt. Return an integer 2734 * >0: this > other 2735 * 0: this == other 2736 * <0: this < other 2737 */ 2738 public int 2739 cmp( FDBigInt other ){ 2740 int i; 2741 if ( this.nWords > other.nWords ){ 2742 // if any of my high-order words is non-zero, 2743 // then the answer is evident 2744 int j = other.nWords-1; 2745 for ( i = this.nWords-1; i > j ; i-- ) 2746 if ( this.data[i] != 0 ) return 1; 2747 }else if ( this.nWords < other.nWords ){ 2748 // if any of other's high-order words is non-zero, 2749 // then the answer is evident 2750 int j = this.nWords-1; 2751 for ( i = other.nWords-1; i > j ; i-- ) 2752 if ( other.data[i] != 0 ) return -1; 2753 } else{ 2754 i = this.nWords-1; 2755 } 2756 for ( ; i > 0 ; i-- ) 2757 if ( this.data[i] != other.data[i] ) 2758 break; 2759 // careful! want unsigned compare! 2760 // use brute force here. 2761 int a = this.data[i]; 2762 int b = other.data[i]; 2763 if ( a < 0 ){ 2764 // a is really big, unsigned 2765 if ( b < 0 ){ 2766 return a-b; // both big, negative 2767 } else { 2768 return 1; // b not big, answer is obvious; 2769 } 2770 } else { 2771 // a is not really big 2772 if ( b < 0 ) { 2773 // but b is really big 2774 return -1; 2775 } else { 2776 return a - b; 2777 } 2778 } 2779 } 2780 2781 /* 2782 * Compute 2783 * q = (int)( this / S ) 2784 * this = 10 * ( this mod S ) 2785 * Return q. 2786 * This is the iteration step of digit development for output. 2787 * We assume that S has been normalized, as above, and that 2788 * "this" has been lshift'ed accordingly. 2789 * Also assume, of course, that the result, q, can be expressed 2790 * as an integer, 0 <= q < 10. 2791 */ 2792 public int 2793 quoRemIteration( FDBigInt S )throws IllegalArgumentException { 2794 // ensure that this and S have the same number of 2795 // digits. If S is properly normalized and q < 10 then 2796 // this must be so. 2797 if ( nWords != S.nWords ){ 2798 throw new IllegalArgumentException("disparate values"); 2799 } 2800 // estimate q the obvious way. We will usually be 2801 // right. If not, then we're only off by a little and 2802 // will re-add. 2803 int n = nWords-1; 2804 long q = ((long)data[n]&0xffffffffL) / (long)S.data[n]; 2805 long diff = 0L; 2806 for ( int i = 0; i <= n ; i++ ){ 2807 diff += ((long)data[i]&0xffffffffL) - q*((long)S.data[i]&0xffffffffL); 2808 data[i] = (int)diff; 2809 diff >>= 32; // N.B. SIGNED shift. 2810 } 2811 if ( diff != 0L ) { 2812 // damn, damn, damn. q is too big. 2813 // add S back in until this turns +. This should 2814 // not be very many times! 2815 long sum = 0L; 2816 while ( sum == 0L ){ 2817 sum = 0L; 2818 for ( int i = 0; i <= n; i++ ){ 2819 sum += ((long)data[i]&0xffffffffL) + ((long)S.data[i]&0xffffffffL); 2820 data[i] = (int) sum; 2821 sum >>= 32; // Signed or unsigned, answer is 0 or 1 2822 } 2823 /* 2824 * Originally the following line read 2825 * "if ( sum !=0 && sum != -1 )" 2826 * but that would be wrong, because of the 2827 * treatment of the two values as entirely unsigned, 2828 * it would be impossible for a carry-out to be interpreted 2829 * as -1 -- it would have to be a single-bit carry-out, or 2830 * +1. 2831 */ 2832 assert sum == 0 || sum == 1 : sum; // carry out of division correction 2833 q -= 1; 2834 } 2835 } 2836 // finally, we can multiply this by 10. 2837 // it cannot overflow, right, as the high-order word has 2838 // at least 4 high-order zeros! 2839 long p = 0L; 2840 for ( int i = 0; i <= n; i++ ){ 2841 p += 10*((long)data[i]&0xffffffffL); 2842 data[i] = (int)p; 2843 p >>= 32; // SIGNED shift. 2844 } 2845 assert p == 0L : p; // Carry out of *10 2846 return (int)q; 2847 } 2848 2849 public long 2850 longValue(){ 2851 // if this can be represented as a long, return the value 2852 assert this.nWords > 0 : this.nWords; // longValue confused 2853 2854 if (this.nWords == 1) 2855 return ((long)data[0]&0xffffffffL); 2856 2857 assert dataInRangeIsZero(2, this.nWords, this); // value too big 2858 assert data[1] >= 0; // value too big 2859 return ((long)(data[1]) << 32) | ((long)data[0]&0xffffffffL); 2860 } 2861 2862 public String 2863 toString() { 2864 StringBuffer r = new StringBuffer(30); 2865 r.append('['); 2866 int i = Math.min( nWords-1, data.length-1) ; 2867 if ( nWords > data.length ){ 2868 r.append( "("+data.length+"<"+nWords+"!)" ); 2869 } 2870 for( ; i> 0 ; i-- ){ 2871 r.append( Integer.toHexString( data[i] ) ); 2872 r.append(' '); 2873 } 2874 r.append( Integer.toHexString( data[0] ) ); 2875 r.append(']'); 2876 return new String( r ); 2877 } 2878 }