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