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