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