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