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