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