1 /*
2 * Copyright (c) 2003, 2011, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
23 * questions.
24 */
25
26 package sun.misc;
27
28 import sun.misc.DoubleConsts;
29 import sun.misc.FloatConsts;
30 import java.util.regex.*;
31
32 public class FormattedFloatingDecimal{
33 boolean isExceptional;
34 boolean isNegative;
35 int decExponent; // value set at construction, then immutable
36 int decExponentRounded;
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 int precision; // number of digits to the right of decimal
45
46 public enum Form { SCIENTIFIC, COMPATIBLE, DECIMAL_FLOAT, GENERAL };
47
48 private Form form;
49
50 private FormattedFloatingDecimal( boolean negSign, int decExponent, char []digits, int n, boolean e, int precision, Form form )
51 {
52 isNegative = negSign;
53 isExceptional = e;
54 this.decExponent = decExponent;
55 this.digits = digits;
56 this.nDigits = n;
57 this.precision = precision;
58 this.form = form;
59 }
60
61 /*
62 * Constants of the implementation
63 * Most are IEEE-754 related.
64 * (There are more really boring constants at the end.)
65 */
66 static final long signMask = 0x8000000000000000L;
67 static final long expMask = 0x7ff0000000000000L;
68 static final long fractMask= ~(signMask|expMask);
69 static final int expShift = 52;
70 static final int expBias = 1023;
71 static final long fractHOB = ( 1L<<expShift ); // assumed High-Order bit
72 static final long expOne = ((long)expBias)<<expShift; // exponent of 1.0
73 static final int maxSmallBinExp = 62;
74 static final int minSmallBinExp = -( 63 / 3 );
75 static final int maxDecimalDigits = 15;
76 static final int maxDecimalExponent = 308;
77 static final int minDecimalExponent = -324;
78 static final int bigDecimalExponent = 324; // i.e. abs(minDecimalExponent)
79
80 static final long highbyte = 0xff00000000000000L;
81 static final long highbit = 0x8000000000000000L;
82 static final long lowbytes = ~highbyte;
83
84 static final int singleSignMask = 0x80000000;
85 static final int singleExpMask = 0x7f800000;
86 static final int singleFractMask = ~(singleSignMask|singleExpMask);
87 static final int singleExpShift = 23;
88 static final int singleFractHOB = 1<<singleExpShift;
89 static final int singleExpBias = 127;
90 static final int singleMaxDecimalDigits = 7;
91 static final int singleMaxDecimalExponent = 38;
92 static final int singleMinDecimalExponent = -45;
93
94 static final int intDecimalDigits = 9;
95
96
97 /*
98 * count number of bits from high-order 1 bit to low-order 1 bit,
99 * inclusive.
100 */
101 private static int
102 countBits( long v ){
103 //
104 // the strategy is to shift until we get a non-zero sign bit
105 // then shift until we have no bits left, counting the difference.
106 // we do byte shifting as a hack. Hope it helps.
107 //
108 if ( v == 0L ) return 0;
109
110 while ( ( v & highbyte ) == 0L ){
111 v <<= 8;
112 }
113 while ( v > 0L ) { // i.e. while ((v&highbit) == 0L )
114 v <<= 1;
115 }
116
117 int n = 0;
118 while (( v & lowbytes ) != 0L ){
119 v <<= 8;
120 n += 8;
121 }
122 while ( v != 0L ){
123 v <<= 1;
124 n += 1;
125 }
126 return n;
127 }
128
129 /*
130 * Keep big powers of 5 handy for future reference.
131 */
132 private static FDBigInt b5p[];
133
134 private static synchronized FDBigInt
135 big5pow( int p ){
136 assert p >= 0 : p; // negative power of 5
137 if ( b5p == null ){
138 b5p = new FDBigInt[ p+1 ];
139 }else if (b5p.length <= p ){
140 FDBigInt t[] = new FDBigInt[ p+1 ];
141 System.arraycopy( b5p, 0, t, 0, b5p.length );
142 b5p = t;
143 }
144 if ( b5p[p] != null )
145 return b5p[p];
146 else if ( p < small5pow.length )
147 return b5p[p] = new FDBigInt( small5pow[p] );
148 else if ( p < long5pow.length )
149 return b5p[p] = new FDBigInt( long5pow[p] );
150 else {
151 // construct the value.
152 // recursively.
153 int q, r;
154 // in order to compute 5^p,
155 // compute its square root, 5^(p/2) and square.
156 // or, let q = p / 2, r = p -q, then
157 // 5^p = 5^(q+r) = 5^q * 5^r
158 q = p >> 1;
159 r = p - q;
160 FDBigInt bigq = b5p[q];
161 if ( bigq == null )
162 bigq = big5pow ( q );
163 if ( r < small5pow.length ){
164 return (b5p[p] = bigq.mult( small5pow[r] ) );
165 }else{
166 FDBigInt bigr = b5p[ r ];
167 if ( bigr == null )
168 bigr = big5pow( r );
169 return (b5p[p] = bigq.mult( bigr ) );
170 }
171 }
172 }
173
174 //
175 // a common operation
176 //
177 private static FDBigInt
178 multPow52( FDBigInt v, int p5, int p2 ){
179 if ( p5 != 0 ){
180 if ( p5 < small5pow.length ){
181 v = v.mult( small5pow[p5] );
182 } else {
183 v = v.mult( big5pow( p5 ) );
184 }
185 }
186 if ( p2 != 0 ){
187 v.lshiftMe( p2 );
188 }
189 return v;
190 }
191
192 //
193 // another common operation
194 //
195 private static FDBigInt
196 constructPow52( int p5, int p2 ){
197 FDBigInt v = new FDBigInt( big5pow( p5 ) );
198 if ( p2 != 0 ){
199 v.lshiftMe( p2 );
200 }
201 return v;
202 }
203
204 /*
205 * Make a floating double into a FDBigInt.
206 * This could also be structured as a FDBigInt
207 * constructor, but we'd have to build a lot of knowledge
208 * about floating-point representation into it, and we don't want to.
209 *
210 * AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
211 * bigIntExp and bigIntNBits
212 *
213 */
214 private FDBigInt
215 doubleToBigInt( double dval ){
216 long lbits = Double.doubleToLongBits( dval ) & ~signMask;
217 int binexp = (int)(lbits >>> expShift);
218 lbits &= fractMask;
219 if ( binexp > 0 ){
220 lbits |= fractHOB;
221 } else {
222 assert lbits != 0L : lbits; // doubleToBigInt(0.0)
223 binexp +=1;
224 while ( (lbits & fractHOB ) == 0L){
225 lbits <<= 1;
226 binexp -= 1;
227 }
228 }
229 binexp -= expBias;
230 int nbits = countBits( lbits );
231 /*
232 * We now know where the high-order 1 bit is,
233 * and we know how many there are.
234 */
235 int lowOrderZeros = expShift+1-nbits;
236 lbits >>>= lowOrderZeros;
237
238 bigIntExp = binexp+1-nbits;
239 bigIntNBits = nbits;
240 return new FDBigInt( lbits );
241 }
242
243 /*
244 * Compute a number that is the ULP of the given value,
245 * for purposes of addition/subtraction. Generally easy.
246 * More difficult if subtracting and the argument
247 * is a normalized a power of 2, as the ULP changes at these points.
248 */
249 private static double ulp( double dval, boolean subtracting ){
250 long lbits = Double.doubleToLongBits( dval ) & ~signMask;
251 int binexp = (int)(lbits >>> expShift);
252 double ulpval;
253 if ( subtracting && ( binexp >= expShift ) && ((lbits&fractMask) == 0L) ){
254 // for subtraction from normalized, powers of 2,
255 // use next-smaller exponent
256 binexp -= 1;
257 }
258 if ( binexp > expShift ){
259 ulpval = Double.longBitsToDouble( ((long)(binexp-expShift))<<expShift );
260 } else if ( binexp == 0 ){
261 ulpval = Double.MIN_VALUE;
262 } else {
263 ulpval = Double.longBitsToDouble( 1L<<(binexp-1) );
264 }
265 if ( subtracting ) ulpval = - ulpval;
266
267 return ulpval;
268 }
269
270 /*
271 * Round a double to a float.
272 * In addition to the fraction bits of the double,
273 * look at the class instance variable roundDir,
274 * which should help us avoid double-rounding error.
275 * roundDir was set in hardValueOf if the estimate was
276 * close enough, but not exact. It tells us which direction
277 * of rounding is preferred.
278 */
279 float
280 stickyRound( double dval ){
281 long lbits = Double.doubleToLongBits( dval );
282 long binexp = lbits & expMask;
283 if ( binexp == 0L || binexp == expMask ){
284 // what we have here is special.
285 // don't worry, the right thing will happen.
286 return (float) dval;
287 }
288 lbits += (long)roundDir; // hack-o-matic.
289 return (float)Double.longBitsToDouble( lbits );
290 }
291
292
293 /*
294 * This is the easy subcase --
295 * all the significant bits, after scaling, are held in lvalue.
296 * negSign and decExponent tell us what processing and scaling
297 * has already been done. Exceptional cases have already been
298 * stripped out.
299 * In particular:
300 * lvalue is a finite number (not Inf, nor NaN)
301 * lvalue > 0L (not zero, nor negative).
302 *
303 * The only reason that we develop the digits here, rather than
304 * calling on Long.toString() is that we can do it a little faster,
305 * and besides want to treat trailing 0s specially. If Long.toString
306 * changes, we should re-evaluate this strategy!
307 */
308 private void
309 developLongDigits( int decExponent, long lvalue, long insignificant ){
310 char digits[];
311 int ndigits;
312 int digitno;
313 int c;
314 //
315 // Discard non-significant low-order bits, while rounding,
316 // up to insignificant value.
317 int i;
318 for ( i = 0; insignificant >= 10L; i++ )
319 insignificant /= 10L;
320 if ( i != 0 ){
321 long pow10 = long5pow[i] << i; // 10^i == 5^i * 2^i;
322 long residue = lvalue % pow10;
323 lvalue /= pow10;
324 decExponent += i;
325 if ( residue >= (pow10>>1) ){
326 // round up based on the low-order bits we're discarding
327 lvalue++;
328 }
329 }
330 if ( lvalue <= Integer.MAX_VALUE ){
331 assert lvalue > 0L : lvalue; // lvalue <= 0
332 // even easier subcase!
333 // can do int arithmetic rather than long!
334 int ivalue = (int)lvalue;
335 ndigits = 10;
336 digits = (char[])(perThreadBuffer.get());
337 digitno = ndigits-1;
338 c = ivalue%10;
339 ivalue /= 10;
340 while ( c == 0 ){
341 decExponent++;
342 c = ivalue%10;
343 ivalue /= 10;
344 }
345 while ( ivalue != 0){
346 digits[digitno--] = (char)(c+'0');
347 decExponent++;
348 c = ivalue%10;
349 ivalue /= 10;
350 }
351 digits[digitno] = (char)(c+'0');
352 } else {
353 // same algorithm as above (same bugs, too )
354 // but using long arithmetic.
355 ndigits = 20;
356 digits = (char[])(perThreadBuffer.get());
357 digitno = ndigits-1;
358 c = (int)(lvalue%10L);
359 lvalue /= 10L;
360 while ( c == 0 ){
361 decExponent++;
362 c = (int)(lvalue%10L);
363 lvalue /= 10L;
364 }
365 while ( lvalue != 0L ){
366 digits[digitno--] = (char)(c+'0');
367 decExponent++;
368 c = (int)(lvalue%10L);
369 lvalue /= 10;
370 }
371 digits[digitno] = (char)(c+'0');
372 }
373 char result [];
374 ndigits -= digitno;
375 result = new char[ ndigits ];
376 System.arraycopy( digits, digitno, result, 0, ndigits );
377 this.digits = result;
378 this.decExponent = decExponent+1;
379 this.nDigits = ndigits;
380 }
381
382 //
383 // add one to the least significant digit.
384 // in the unlikely event there is a carry out,
385 // deal with it.
386 // assert that this will only happen where there
387 // is only one digit, e.g. (float)1e-44 seems to do it.
388 //
389 private void
390 roundup(){
391 int i;
392 int q = digits[ i = (nDigits-1)];
393 if ( q == '9' ){
394 while ( q == '9' && i > 0 ){
395 digits[i] = '0';
396 q = digits[--i];
397 }
398 if ( q == '9' ){
399 // carryout! High-order 1, rest 0s, larger exp.
400 decExponent += 1;
401 digits[0] = '1';
402 return;
403 }
404 // else fall through.
405 }
406 digits[i] = (char)(q+1);
407 }
408
409 // Given the desired number of digits predict the result's exponent.
410 private int checkExponent(int length) {
411 if (length >= nDigits || length < 0)
412 return decExponent;
413
414 for (int i = 0; i < length; i++)
415 if (digits[i] != '9')
416 // a '9' anywhere in digits will absorb the round
417 return decExponent;
418 return decExponent + (digits[length] >= '5' ? 1 : 0);
419 }
420
421 // Unlike roundup(), this method does not modify digits. It also
422 // rounds at a particular precision.
423 private char [] applyPrecision(int length) {
424 char [] result = new char[nDigits];
425 for (int i = 0; i < result.length; i++) result[i] = '0';
426
427 if (length >= nDigits || length < 0) {
428 // no rounding necessary
429 System.arraycopy(digits, 0, result, 0, nDigits);
430 return result;
431 }
432 if (length == 0) {
433 // only one digit (0 or 1) is returned because the precision
434 // excludes all significant digits
435 if (digits[0] >= '5') {
436 result[0] = '1';
437 }
438 return result;
439 }
440
441 int i = length;
442 int q = digits[i];
443 if (q >= '5' && i > 0) {
444 q = digits[--i];
445 if ( q == '9' ) {
446 while ( q == '9' && i > 0 ){
447 q = digits[--i];
448 }
449 if ( q == '9' ){
450 // carryout! High-order 1, rest 0s, larger exp.
451 result[0] = '1';
452 return result;
453 }
454 }
455 result[i] = (char)(q + 1);
456 }
457 while (--i >= 0) {
458 result[i] = digits[i];
459 }
460 return result;
461 }
462
463 /*
464 * FIRST IMPORTANT CONSTRUCTOR: DOUBLE
465 */
466 public FormattedFloatingDecimal( double d )
467 {
468 this(d, Integer.MAX_VALUE, Form.COMPATIBLE);
469 }
470
471 public FormattedFloatingDecimal( double d, int precision, Form form )
472 {
473 long dBits = Double.doubleToLongBits( d );
474 long fractBits;
475 int binExp;
476 int nSignificantBits;
477
478 this.precision = precision;
479 this.form = form;
480
481 // discover and delete sign
482 if ( (dBits&signMask) != 0 ){
483 isNegative = true;
484 dBits ^= signMask;
485 } else {
486 isNegative = false;
487 }
488 // Begin to unpack
489 // Discover obvious special cases of NaN and Infinity.
490 binExp = (int)( (dBits&expMask) >> expShift );
491 fractBits = dBits&fractMask;
492 if ( binExp == (int)(expMask>>expShift) ) {
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 == 0L ){
510 // not a denorm, just a 0!
511 decExponent = 0;
512 digits = zero;
513 nDigits = 1;
514 return;
515 }
516 while ( (fractBits&fractHOB) == 0L ){
517 fractBits <<= 1;
518 binExp -= 1;
519 }
520 nSignificantBits = expShift + binExp +1; // recall binExp is - shift count.
521 binExp += 1;
522 } else {
523 fractBits |= fractHOB;
524 nSignificantBits = expShift+1;
525 }
526 binExp -= expBias;
527 // call the routine that actually does all the hard work.
528 dtoa( binExp, fractBits, nSignificantBits );
529 }
530
531 /*
532 * SECOND IMPORTANT CONSTRUCTOR: SINGLE
533 */
534 public FormattedFloatingDecimal( float f )
535 {
536 this(f, Integer.MAX_VALUE, Form.COMPATIBLE);
537 }
538 public FormattedFloatingDecimal( float f, int precision, Form form )
539 {
540 int fBits = Float.floatToIntBits( f );
541 int fractBits;
542 int binExp;
543 int nSignificantBits;
544
545 this.precision = precision;
546 this.form = form;
547
548 // discover and delete sign
549 if ( (fBits&singleSignMask) != 0 ){
550 isNegative = true;
551 fBits ^= singleSignMask;
552 } else {
553 isNegative = false;
554 }
555 // Begin to unpack
556 // Discover obvious special cases of NaN and Infinity.
557 binExp = (int)( (fBits&singleExpMask) >> singleExpShift );
558 fractBits = fBits&singleFractMask;
559 if ( binExp == (int)(singleExpMask>>singleExpShift) ) {
560 isExceptional = true;
561 if ( fractBits == 0L ){
562 digits = infinity;
563 } else {
564 digits = notANumber;
565 isNegative = false; // NaN has no sign!
566 }
567 nDigits = digits.length;
568 return;
569 }
570 isExceptional = false;
571 // Finish unpacking
572 // Normalize denormalized numbers.
573 // Insert assumed high-order bit for normalized numbers.
574 // Subtract exponent bias.
575 if ( binExp == 0 ){
576 if ( fractBits == 0 ){
577 // not a denorm, just a 0!
578 decExponent = 0;
579 digits = zero;
580 nDigits = 1;
581 return;
582 }
583 while ( (fractBits&singleFractHOB) == 0 ){
584 fractBits <<= 1;
585 binExp -= 1;
586 }
587 nSignificantBits = singleExpShift + binExp +1; // recall binExp is - shift count.
588 binExp += 1;
589 } else {
590 fractBits |= singleFractHOB;
591 nSignificantBits = singleExpShift+1;
592 }
593 binExp -= singleExpBias;
594 // call the routine that actually does all the hard work.
595 dtoa( binExp, ((long)fractBits)<<(expShift-singleExpShift), nSignificantBits );
596 }
597
598 private void
599 dtoa( int binExp, long fractBits, int nSignificantBits )
600 {
601 int nFractBits; // number of significant bits of fractBits;
602 int nTinyBits; // number of these to the right of the point.
603 int decExp;
604
605 // Examine number. Determine if it is an easy case,
606 // which we can do pretty trivially using float/long conversion,
607 // or whether we must do real work.
608 nFractBits = countBits( fractBits );
609 nTinyBits = Math.max( 0, nFractBits - binExp - 1 );
610 if ( binExp <= maxSmallBinExp && binExp >= minSmallBinExp ){
611 // Look more closely at the number to decide if,
612 // with scaling by 10^nTinyBits, the result will fit in
613 // a long.
614 if ( (nTinyBits < long5pow.length) && ((nFractBits + n5bits[nTinyBits]) < 64 ) ){
615 /*
616 * We can do this:
617 * take the fraction bits, which are normalized.
618 * (a) nTinyBits == 0: Shift left or right appropriately
619 * to align the binary point at the extreme right, i.e.
620 * where a long int point is expected to be. The integer
621 * result is easily converted to a string.
622 * (b) nTinyBits > 0: Shift right by expShift-nFractBits,
623 * which effectively converts to long and scales by
624 * 2^nTinyBits. Then multiply by 5^nTinyBits to
625 * complete the scaling. We know this won't overflow
626 * because we just counted the number of bits necessary
627 * in the result. The integer you get from this can
628 * then be converted to a string pretty easily.
629 */
630 long halfULP;
631 if ( nTinyBits == 0 ) {
632 if ( binExp > nSignificantBits ){
633 halfULP = 1L << ( binExp-nSignificantBits-1);
634 } else {
635 halfULP = 0L;
636 }
637 if ( binExp >= expShift ){
638 fractBits <<= (binExp-expShift);
639 } else {
640 fractBits >>>= (expShift-binExp) ;
641 }
642 developLongDigits( 0, fractBits, halfULP );
643 return;
644 }
645 /*
646 * The following causes excess digits to be printed
647 * out in the single-float case. Our manipulation of
648 * halfULP here is apparently not correct. If we
649 * better understand how this works, perhaps we can
650 * use this special case again. But for the time being,
651 * we do not.
652 * else {
653 * fractBits >>>= expShift+1-nFractBits;
654 * fractBits *= long5pow[ nTinyBits ];
655 * halfULP = long5pow[ nTinyBits ] >> (1+nSignificantBits-nFractBits);
656 * developLongDigits( -nTinyBits, fractBits, halfULP );
657 * return;
658 * }
659 */
660 }
661 }
662 /*
663 * This is the hard case. We are going to compute large positive
664 * integers B and S and integer decExp, s.t.
665 * d = ( B / S ) * 10^decExp
666 * 1 <= B / S < 10
667 * Obvious choices are:
668 * decExp = floor( log10(d) )
669 * B = d * 2^nTinyBits * 10^max( 0, -decExp )
670 * S = 10^max( 0, decExp) * 2^nTinyBits
671 * (noting that nTinyBits has already been forced to non-negative)
672 * I am also going to compute a large positive integer
673 * M = (1/2^nSignificantBits) * 2^nTinyBits * 10^max( 0, -decExp )
674 * i.e. M is (1/2) of the ULP of d, scaled like B.
675 * When we iterate through dividing B/S and picking off the
676 * quotient bits, we will know when to stop when the remainder
677 * is <= M.
678 *
679 * We keep track of powers of 2 and powers of 5.
680 */
681
682 /*
683 * Estimate decimal exponent. (If it is small-ish,
684 * we could double-check.)
685 *
686 * First, scale the mantissa bits such that 1 <= d2 < 2.
687 * We are then going to estimate
688 * log10(d2) ~=~ (d2-1.5)/1.5 + log(1.5)
689 * and so we can estimate
690 * log10(d) ~=~ log10(d2) + binExp * log10(2)
691 * take the floor and call it decExp.
692 * FIXME -- use more precise constants here. It costs no more.
693 */
694 double d2 = Double.longBitsToDouble(
695 expOne | ( fractBits &~ fractHOB ) );
696 decExp = (int)Math.floor(
697 (d2-1.5D)*0.289529654D + 0.176091259 + (double)binExp * 0.301029995663981 );
698 int B2, B5; // powers of 2 and powers of 5, respectively, in B
699 int S2, S5; // powers of 2 and powers of 5, respectively, in S
700 int M2, M5; // powers of 2 and powers of 5, respectively, in M
701 int Bbits; // binary digits needed to represent B, approx.
702 int tenSbits; // binary digits needed to represent 10*S, approx.
703 FDBigInt Sval, Bval, Mval;
704
705 B5 = Math.max( 0, -decExp );
706 B2 = B5 + nTinyBits + binExp;
707
708 S5 = Math.max( 0, decExp );
709 S2 = S5 + nTinyBits;
710
711 M5 = B5;
712 M2 = B2 - nSignificantBits;
713
714 /*
715 * the long integer fractBits contains the (nFractBits) interesting
716 * bits from the mantissa of d ( hidden 1 added if necessary) followed
717 * by (expShift+1-nFractBits) zeros. In the interest of compactness,
718 * I will shift out those zeros before turning fractBits into a
719 * FDBigInt. The resulting whole number will be
720 * d * 2^(nFractBits-1-binExp).
721 */
722 fractBits >>>= (expShift+1-nFractBits);
723 B2 -= nFractBits-1;
724 int common2factor = Math.min( B2, S2 );
725 B2 -= common2factor;
726 S2 -= common2factor;
727 M2 -= common2factor;
728
729 /*
730 * HACK!! For exact powers of two, the next smallest number
731 * is only half as far away as we think (because the meaning of
732 * ULP changes at power-of-two bounds) for this reason, we
733 * hack M2. Hope this works.
734 */
735 if ( nFractBits == 1 )
736 M2 -= 1;
737
738 if ( M2 < 0 ){
739 // oops.
740 // since we cannot scale M down far enough,
741 // we must scale the other values up.
742 B2 -= M2;
743 S2 -= M2;
744 M2 = 0;
745 }
746 /*
747 * Construct, Scale, iterate.
748 * Some day, we'll write a stopping test that takes
749 * account of the assymetry of the spacing of floating-point
750 * numbers below perfect powers of 2
751 * 26 Sept 96 is not that day.
752 * So we use a symmetric test.
753 */
754 char digits[] = this.digits = new char[18];
755 int ndigit = 0;
756 boolean low, high;
757 long lowDigitDifference;
758 int q;
759
760 /*
761 * Detect the special cases where all the numbers we are about
762 * to compute will fit in int or long integers.
763 * In these cases, we will avoid doing FDBigInt arithmetic.
764 * We use the same algorithms, except that we "normalize"
765 * our FDBigInts before iterating. This is to make division easier,
766 * as it makes our fist guess (quotient of high-order words)
767 * more accurate!
768 *
769 * Some day, we'll write a stopping test that takes
770 * account of the assymetry of the spacing of floating-point
771 * numbers below perfect powers of 2
772 * 26 Sept 96 is not that day.
773 * So we use a symmetric test.
774 */
775 Bbits = nFractBits + B2 + (( B5 < n5bits.length )? n5bits[B5] : ( B5*3 ));
776 tenSbits = S2+1 + (( (S5+1) < n5bits.length )? n5bits[(S5+1)] : ( (S5+1)*3 ));
777 if ( Bbits < 64 && tenSbits < 64){
778 if ( Bbits < 32 && tenSbits < 32){
779 // wa-hoo! They're all ints!
780 int b = ((int)fractBits * small5pow[B5] ) << B2;
781 int s = small5pow[S5] << S2;
782 int m = small5pow[M5] << M2;
783 int tens = s * 10;
784 /*
785 * Unroll the first iteration. If our decExp estimate
786 * was too high, our first quotient will be zero. In this
787 * case, we discard it and decrement decExp.
788 */
789 ndigit = 0;
790 q = b / s;
791 b = 10 * ( b % s );
792 m *= 10;
793 low = (b < m );
794 high = (b+m > tens );
795 assert q < 10 : q; // excessively large digit
796 if ( (q == 0) && ! high ){
797 // oops. Usually ignore leading zero.
798 decExp--;
799 } else {
800 digits[ndigit++] = (char)('0' + q);
801 }
802 /*
803 * HACK! Java spec sez that we always have at least
804 * one digit after the . in either F- or E-form output.
805 * Thus we will need more than one digit if we're using
806 * E-form
807 */
808 if (! (form == Form.COMPATIBLE && -3 < decExp && decExp < 8)) {
809 high = low = false;
810 }
811 while( ! low && ! high ){
812 q = b / s;
813 b = 10 * ( b % s );
814 m *= 10;
815 assert q < 10 : q; // excessively large digit
816 if ( m > 0L ){
817 low = (b < m );
818 high = (b+m > tens );
819 } else {
820 // hack -- m might overflow!
821 // in this case, it is certainly > b,
822 // which won't
823 // and b+m > tens, too, since that has overflowed
824 // either!
825 low = true;
826 high = true;
827 }
828 digits[ndigit++] = (char)('0' + q);
829 }
830 lowDigitDifference = (b<<1) - tens;
831 } else {
832 // still good! they're all longs!
833 long b = (fractBits * long5pow[B5] ) << B2;
834 long s = long5pow[S5] << S2;
835 long m = long5pow[M5] << M2;
836 long tens = s * 10L;
837 /*
838 * Unroll the first iteration. If our decExp estimate
839 * was too high, our first quotient will be zero. In this
840 * case, we discard it and decrement decExp.
841 */
842 ndigit = 0;
843 q = (int) ( b / s );
844 b = 10L * ( b % s );
845 m *= 10L;
846 low = (b < m );
847 high = (b+m > tens );
848 assert q < 10 : q; // excessively large digit
849 if ( (q == 0) && ! high ){
850 // oops. Usually ignore leading zero.
851 decExp--;
852 } else {
853 digits[ndigit++] = (char)('0' + q);
854 }
855 /*
856 * HACK! Java spec sez that we always have at least
857 * one digit after the . in either F- or E-form output.
858 * Thus we will need more than one digit if we're using
859 * E-form
860 */
861 if (! (form == Form.COMPATIBLE && -3 < decExp && decExp < 8)) {
862 high = low = false;
863 }
864 while( ! low && ! high ){
865 q = (int) ( b / s );
866 b = 10 * ( b % s );
867 m *= 10;
868 assert q < 10 : q; // excessively large digit
869 if ( m > 0L ){
870 low = (b < m );
871 high = (b+m > tens );
872 } else {
873 // hack -- m might overflow!
874 // in this case, it is certainly > b,
875 // which won't
876 // and b+m > tens, too, since that has overflowed
877 // either!
878 low = true;
879 high = true;
880 }
881 digits[ndigit++] = (char)('0' + q);
882 }
883 lowDigitDifference = (b<<1) - tens;
884 }
885 } else {
886 FDBigInt tenSval;
887 int shiftBias;
888
889 /*
890 * We really must do FDBigInt arithmetic.
891 * Fist, construct our FDBigInt initial values.
892 */
893 Bval = multPow52( new FDBigInt( fractBits ), B5, B2 );
894 Sval = constructPow52( S5, S2 );
895 Mval = constructPow52( M5, M2 );
896
897
898 // normalize so that division works better
899 Bval.lshiftMe( shiftBias = Sval.normalizeMe() );
900 Mval.lshiftMe( shiftBias );
901 tenSval = Sval.mult( 10 );
902 /*
903 * Unroll the first iteration. If our decExp estimate
904 * was too high, our first quotient will be zero. In this
905 * case, we discard it and decrement decExp.
906 */
907 ndigit = 0;
908 q = Bval.quoRemIteration( Sval );
909 Mval = Mval.mult( 10 );
910 low = (Bval.cmp( Mval ) < 0);
911 high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
912 assert q < 10 : q; // excessively large digit
913 if ( (q == 0) && ! high ){
914 // oops. Usually ignore leading zero.
915 decExp--;
916 } else {
917 digits[ndigit++] = (char)('0' + q);
918 }
919 /*
920 * HACK! Java spec sez that we always have at least
921 * one digit after the . in either F- or E-form output.
922 * Thus we will need more than one digit if we're using
923 * E-form
924 */
925 if (! (form == Form.COMPATIBLE && -3 < decExp && decExp < 8)) {
926 high = low = false;
927 }
928 while( ! low && ! high ){
929 q = Bval.quoRemIteration( Sval );
930 Mval = Mval.mult( 10 );
931 assert q < 10 : q; // excessively large digit
932 low = (Bval.cmp( Mval ) < 0);
933 high = (Bval.add( Mval ).cmp( tenSval ) > 0 );
934 digits[ndigit++] = (char)('0' + q);
935 }
936 if ( high && low ){
937 Bval.lshiftMe(1);
938 lowDigitDifference = Bval.cmp(tenSval);
939 } else
940 lowDigitDifference = 0L; // this here only for flow analysis!
941 }
942 this.decExponent = decExp+1;
943 this.digits = digits;
944 this.nDigits = ndigit;
945 /*
946 * Last digit gets rounded based on stopping condition.
947 */
948 if ( high ){
949 if ( low ){
950 if ( lowDigitDifference == 0L ){
951 // it's a tie!
952 // choose based on which digits we like.
953 if ( (digits[nDigits-1]&1) != 0 ) roundup();
954 } else if ( lowDigitDifference > 0 ){
955 roundup();
956 }
957 } else {
958 roundup();
959 }
960 }
961 }
962
963 public String
964 toString(){
965 // most brain-dead version
966 StringBuffer result = new StringBuffer( nDigits+8 );
967 if ( isNegative ){ result.append( '-' ); }
968 if ( isExceptional ){
969 result.append( digits, 0, nDigits );
970 } else {
971 result.append( "0.");
972 result.append( digits, 0, nDigits );
973 result.append('e');
974 result.append( decExponent );
975 }
976 return new String(result);
977 }
978
979 // returns the exponent before rounding
980 public int getExponent() {
981 return decExponent - 1;
982 }
983
984 // returns the exponent after rounding has been done by applyPrecision
985 public int getExponentRounded() {
986 return decExponentRounded - 1;
987 }
988
989 public int getChars(char[] result) {
990 assert nDigits <= 19 : nDigits; // generous bound on size of nDigits
991 int i = 0;
992 if (isNegative) { result[0] = '-'; i = 1; }
993 if (isExceptional) {
994 System.arraycopy(digits, 0, result, i, nDigits);
995 i += nDigits;
996 } else {
997 char digits [] = this.digits;
998 int exp = decExponent;
999 switch (form) {
1000 case COMPATIBLE:
1001 break;
1002 case DECIMAL_FLOAT:
1003 exp = checkExponent(decExponent + precision);
1004 digits = applyPrecision(decExponent + precision);
1005 break;
1006 case SCIENTIFIC:
1007 exp = checkExponent(precision + 1);
1008 digits = applyPrecision(precision + 1);
1009 break;
1010 case GENERAL:
1011 exp = checkExponent(precision);
1012 digits = applyPrecision(precision);
1013 // adjust precision to be the number of digits to right of decimal
1014 // the real exponent to be output is actually exp - 1, not exp
1015 if (exp - 1 < -4 || exp - 1 >= precision) {
1016 form = Form.SCIENTIFIC;
1017 precision--;
1018 } else {
1019 form = Form.DECIMAL_FLOAT;
1020 precision = precision - exp;
1021 }
1022 break;
1023 default:
1024 assert false;
1025 }
1026 decExponentRounded = exp;
1027
1028 if (exp > 0
1029 && ((form == Form.COMPATIBLE && (exp < 8))
1030 || (form == Form.DECIMAL_FLOAT)))
1031 {
1032 // print digits.digits.
1033 int charLength = Math.min(nDigits, exp);
1034 System.arraycopy(digits, 0, result, i, charLength);
1035 i += charLength;
1036 if (charLength < exp) {
1037 charLength = exp-charLength;
1038 for (int nz = 0; nz < charLength; nz++)
1039 result[i++] = '0';
1040 // Do not append ".0" for formatted floats since the user
1041 // may request that it be omitted. It is added as necessary
1042 // by the Formatter.
1043 if (form == Form.COMPATIBLE) {
1044 result[i++] = '.';
1045 result[i++] = '0';
1046 }
1047 } else {
1048 // Do not append ".0" for formatted floats since the user
1049 // may request that it be omitted. It is added as necessary
1050 // by the Formatter.
1051 if (form == Form.COMPATIBLE) {
1052 result[i++] = '.';
1053 if (charLength < nDigits) {
1054 int t = Math.min(nDigits - charLength, precision);
1055 System.arraycopy(digits, charLength, result, i, t);
1056 i += t;
1057 } else {
1058 result[i++] = '0';
1059 }
1060 } else {
1061 int t = Math.min(nDigits - charLength, precision);
1062 if (t > 0) {
1063 result[i++] = '.';
1064 System.arraycopy(digits, charLength, result, i, t);
1065 i += t;
1066 }
1067 }
1068 }
1069 } else if (exp <= 0
1070 && ((form == Form.COMPATIBLE && exp > -3)
1071 || (form == Form.DECIMAL_FLOAT)))
1072 {
1073 // print 0.0* digits
1074 result[i++] = '0';
1075 if (exp != 0) {
1076 // write '0' s before the significant digits
1077 int t = Math.min(-exp, precision);
1078 if (t > 0) {
1079 result[i++] = '.';
1080 for (int nz = 0; nz < t; nz++)
1081 result[i++] = '0';
1082 }
1083 }
1084 int t = Math.min(digits.length, precision + exp);
1085 if (t > 0) {
1086 if (i == 1)
1087 result[i++] = '.';
1088 // copy only when significant digits are within the precision
1089 System.arraycopy(digits, 0, result, i, t);
1090 i += t;
1091 }
1092 } else {
1093 result[i++] = digits[0];
1094 if (form == Form.COMPATIBLE) {
1095 result[i++] = '.';
1096 if (nDigits > 1) {
1097 System.arraycopy(digits, 1, result, i, nDigits-1);
1098 i += nDigits-1;
1099 } else {
1100 result[i++] = '0';
1101 }
1102 result[i++] = 'E';
1103 } else {
1104 if (nDigits > 1) {
1105 int t = Math.min(nDigits -1, precision);
1106 if (t > 0) {
1107 result[i++] = '.';
1108 System.arraycopy(digits, 1, result, i, t);
1109 i += t;
1110 }
1111 }
1112 result[i++] = 'e';
1113 }
1114 int e;
1115 if (exp <= 0) {
1116 result[i++] = '-';
1117 e = -exp+1;
1118 } else {
1119 if (form != Form.COMPATIBLE)
1120 result[i++] = '+';
1121 e = exp-1;
1122 }
1123 // decExponent has 1, 2, or 3, digits
1124 if (e <= 9) {
1125 if (form != Form.COMPATIBLE)
1126 result[i++] = '0';
1127 result[i++] = (char)(e+'0');
1128 } else if (e <= 99) {
1129 result[i++] = (char)(e/10 +'0');
1130 result[i++] = (char)(e%10 + '0');
1131 } else {
1132 result[i++] = (char)(e/100+'0');
1133 e %= 100;
1134 result[i++] = (char)(e/10+'0');
1135 result[i++] = (char)(e%10 + '0');
1136 }
1137 }
1138 }
1139 return i;
1140 }
1141
1142 // Per-thread buffer for string/stringbuffer conversion
1143 private static ThreadLocal perThreadBuffer = new ThreadLocal() {
1144 protected synchronized Object initialValue() {
1145 return new char[26];
1146 }
1147 };
1148
1149 /*
1150 * Take a FormattedFloatingDecimal, which we presumably just scanned in,
1151 * and find out what its value is, as a double.
1152 *
1153 * AS A SIDE EFFECT, SET roundDir TO INDICATE PREFERRED
1154 * ROUNDING DIRECTION in case the result is really destined
1155 * for a single-precision float.
1156 */
1157
1158 public strictfp double doubleValue(){
1159 int kDigits = Math.min( nDigits, maxDecimalDigits+1 );
1160 long lValue;
1161 double dValue;
1162 double rValue, tValue;
1163
1164 // First, check for NaN and Infinity values
1165 if(digits == infinity || digits == notANumber) {
1166 if(digits == notANumber)
1167 return Double.NaN;
1168 else
1169 return (isNegative?Double.NEGATIVE_INFINITY:Double.POSITIVE_INFINITY);
1170 }
1171 else {
1172 if (mustSetRoundDir) {
1173 roundDir = 0;
1174 }
1175 /*
1176 * convert the lead kDigits to a long integer.
1177 */
1178 // (special performance hack: start to do it using int)
1179 int iValue = (int)digits[0]-(int)'0';
1180 int iDigits = Math.min( kDigits, intDecimalDigits );
1181 for ( int i=1; i < iDigits; i++ ){
1182 iValue = iValue*10 + (int)digits[i]-(int)'0';
1183 }
1184 lValue = (long)iValue;
1185 for ( int i=iDigits; i < kDigits; i++ ){
1186 lValue = lValue*10L + (long)((int)digits[i]-(int)'0');
1187 }
1188 dValue = (double)lValue;
1189 int exp = decExponent-kDigits;
1190 /*
1191 * lValue now contains a long integer with the value of
1192 * the first kDigits digits of the number.
1193 * dValue contains the (double) of the same.
1194 */
1195
1196 if ( nDigits <= maxDecimalDigits ){
1197 /*
1198 * possibly an easy case.
1199 * We know that the digits can be represented
1200 * exactly. And if the exponent isn't too outrageous,
1201 * the whole thing can be done with one operation,
1202 * thus one rounding error.
1203 * Note that all our constructors trim all leading and
1204 * trailing zeros, so simple values (including zero)
1205 * will always end up here
1206 */
1207 if (exp == 0 || dValue == 0.0)
1208 return (isNegative)? -dValue : dValue; // small floating integer
1209 else if ( exp >= 0 ){
1210 if ( exp <= maxSmallTen ){
1211 /*
1212 * Can get the answer with one operation,
1213 * thus one roundoff.
1214 */
1215 rValue = dValue * small10pow[exp];
1216 if ( mustSetRoundDir ){
1217 tValue = rValue / small10pow[exp];
1218 roundDir = ( tValue == dValue ) ? 0
1219 :( tValue < dValue ) ? 1
1220 : -1;
1221 }
1222 return (isNegative)? -rValue : rValue;
1223 }
1224 int slop = maxDecimalDigits - kDigits;
1225 if ( exp <= maxSmallTen+slop ){
1226 /*
1227 * We can multiply dValue by 10^(slop)
1228 * and it is still "small" and exact.
1229 * Then we can multiply by 10^(exp-slop)
1230 * with one rounding.
1231 */
1232 dValue *= small10pow[slop];
1233 rValue = dValue * small10pow[exp-slop];
1234
1235 if ( mustSetRoundDir ){
1236 tValue = rValue / small10pow[exp-slop];
1237 roundDir = ( tValue == dValue ) ? 0
1238 :( tValue < dValue ) ? 1
1239 : -1;
1240 }
1241 return (isNegative)? -rValue : rValue;
1242 }
1243 /*
1244 * Else we have a hard case with a positive exp.
1245 */
1246 } else {
1247 if ( exp >= -maxSmallTen ){
1248 /*
1249 * Can get the answer in one division.
1250 */
1251 rValue = dValue / small10pow[-exp];
1252 tValue = rValue * small10pow[-exp];
1253 if ( mustSetRoundDir ){
1254 roundDir = ( tValue == dValue ) ? 0
1255 :( tValue < dValue ) ? 1
1256 : -1;
1257 }
1258 return (isNegative)? -rValue : rValue;
1259 }
1260 /*
1261 * Else we have a hard case with a negative exp.
1262 */
1263 }
1264 }
1265
1266 /*
1267 * Harder cases:
1268 * The sum of digits plus exponent is greater than
1269 * what we think we can do with one error.
1270 *
1271 * Start by approximating the right answer by,
1272 * naively, scaling by powers of 10.
1273 */
1274 if ( exp > 0 ){
1275 if ( decExponent > maxDecimalExponent+1 ){
1276 /*
1277 * Lets face it. This is going to be
1278 * Infinity. Cut to the chase.
1279 */
1280 return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1281 }
1282 if ( (exp&15) != 0 ){
1283 dValue *= small10pow[exp&15];
1284 }
1285 if ( (exp>>=4) != 0 ){
1286 int j;
1287 for( j = 0; exp > 1; j++, exp>>=1 ){
1288 if ( (exp&1)!=0)
1289 dValue *= big10pow[j];
1290 }
1291 /*
1292 * The reason for the weird exp > 1 condition
1293 * in the above loop was so that the last multiply
1294 * would get unrolled. We handle it here.
1295 * It could overflow.
1296 */
1297 double t = dValue * big10pow[j];
1298 if ( Double.isInfinite( t ) ){
1299 /*
1300 * It did overflow.
1301 * Look more closely at the result.
1302 * If the exponent is just one too large,
1303 * then use the maximum finite as our estimate
1304 * value. Else call the result infinity
1305 * and punt it.
1306 * ( I presume this could happen because
1307 * rounding forces the result here to be
1308 * an ULP or two larger than
1309 * Double.MAX_VALUE ).
1310 */
1311 t = dValue / 2.0;
1312 t *= big10pow[j];
1313 if ( Double.isInfinite( t ) ){
1314 return (isNegative)? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1315 }
1316 t = Double.MAX_VALUE;
1317 }
1318 dValue = t;
1319 }
1320 } else if ( exp < 0 ){
1321 exp = -exp;
1322 if ( decExponent < minDecimalExponent-1 ){
1323 /*
1324 * Lets face it. This is going to be
1325 * zero. Cut to the chase.
1326 */
1327 return (isNegative)? -0.0 : 0.0;
1328 }
1329 if ( (exp&15) != 0 ){
1330 dValue /= small10pow[exp&15];
1331 }
1332 if ( (exp>>=4) != 0 ){
1333 int j;
1334 for( j = 0; exp > 1; j++, exp>>=1 ){
1335 if ( (exp&1)!=0)
1336 dValue *= tiny10pow[j];
1337 }
1338 /*
1339 * The reason for the weird exp > 1 condition
1340 * in the above loop was so that the last multiply
1341 * would get unrolled. We handle it here.
1342 * It could underflow.
1343 */
1344 double t = dValue * tiny10pow[j];
1345 if ( t == 0.0 ){
1346 /*
1347 * It did underflow.
1348 * Look more closely at the result.
1349 * If the exponent is just one too small,
1350 * then use the minimum finite as our estimate
1351 * value. Else call the result 0.0
1352 * and punt it.
1353 * ( I presume this could happen because
1354 * rounding forces the result here to be
1355 * an ULP or two less than
1356 * Double.MIN_VALUE ).
1357 */
1358 t = dValue * 2.0;
1359 t *= tiny10pow[j];
1360 if ( t == 0.0 ){
1361 return (isNegative)? -0.0 : 0.0;
1362 }
1363 t = Double.MIN_VALUE;
1364 }
1365 dValue = t;
1366 }
1367 }
1368
1369 /*
1370 * dValue is now approximately the result.
1371 * The hard part is adjusting it, by comparison
1372 * with FDBigInt arithmetic.
1373 * Formulate the EXACT big-number result as
1374 * bigD0 * 10^exp
1375 */
1376 FDBigInt bigD0 = new FDBigInt( lValue, digits, kDigits, nDigits );
1377 exp = decExponent - nDigits;
1378
1379 correctionLoop:
1380 while(true){
1381 /* AS A SIDE EFFECT, THIS METHOD WILL SET THE INSTANCE VARIABLES
1382 * bigIntExp and bigIntNBits
1383 */
1384 FDBigInt bigB = doubleToBigInt( dValue );
1385
1386 /*
1387 * Scale bigD, bigB appropriately for
1388 * big-integer operations.
1389 * Naively, we multipy by powers of ten
1390 * and powers of two. What we actually do
1391 * is keep track of the powers of 5 and
1392 * powers of 2 we would use, then factor out
1393 * common divisors before doing the work.
1394 */
1395 int B2, B5; // powers of 2, 5 in bigB
1396 int D2, D5; // powers of 2, 5 in bigD
1397 int Ulp2; // powers of 2 in halfUlp.
1398 if ( exp >= 0 ){
1399 B2 = B5 = 0;
1400 D2 = D5 = exp;
1401 } else {
1402 B2 = B5 = -exp;
1403 D2 = D5 = 0;
1404 }
1405 if ( bigIntExp >= 0 ){
1406 B2 += bigIntExp;
1407 } else {
1408 D2 -= bigIntExp;
1409 }
1410 Ulp2 = B2;
1411 // shift bigB and bigD left by a number s. t.
1412 // halfUlp is still an integer.
1413 int hulpbias;
1414 if ( bigIntExp+bigIntNBits <= -expBias+1 ){
1415 // This is going to be a denormalized number
1416 // (if not actually zero).
1417 // half an ULP is at 2^-(expBias+expShift+1)
1418 hulpbias = bigIntExp+ expBias + expShift;
1419 } else {
1420 hulpbias = expShift + 2 - bigIntNBits;
1421 }
1422 B2 += hulpbias;
1423 D2 += hulpbias;
1424 // if there are common factors of 2, we might just as well
1425 // factor them out, as they add nothing useful.
1426 int common2 = Math.min( B2, Math.min( D2, Ulp2 ) );
1427 B2 -= common2;
1428 D2 -= common2;
1429 Ulp2 -= common2;
1430 // do multiplications by powers of 5 and 2
1431 bigB = multPow52( bigB, B5, B2 );
1432 FDBigInt bigD = multPow52( new FDBigInt( bigD0 ), D5, D2 );
1433 //
1434 // to recap:
1435 // bigB is the scaled-big-int version of our floating-point
1436 // candidate.
1437 // bigD is the scaled-big-int version of the exact value
1438 // as we understand it.
1439 // halfUlp is 1/2 an ulp of bigB, except for special cases
1440 // of exact powers of 2
1441 //
1442 // the plan is to compare bigB with bigD, and if the difference
1443 // is less than halfUlp, then we're satisfied. Otherwise,
1444 // use the ratio of difference to halfUlp to calculate a fudge
1445 // factor to add to the floating value, then go 'round again.
1446 //
1447 FDBigInt diff;
1448 int cmpResult;
1449 boolean overvalue;
1450 if ( (cmpResult = bigB.cmp( bigD ) ) > 0 ){
1451 overvalue = true; // our candidate is too big.
1452 diff = bigB.sub( bigD );
1453 if ( (bigIntNBits == 1) && (bigIntExp > -expBias) ){
1454 // candidate is a normalized exact power of 2 and
1455 // is too big. We will be subtracting.
1456 // For our purposes, ulp is the ulp of the
1457 // next smaller range.
1458 Ulp2 -= 1;
1459 if ( Ulp2 < 0 ){
1460 // rats. Cannot de-scale ulp this far.
1461 // must scale diff in other direction.
1462 Ulp2 = 0;
1463 diff.lshiftMe( 1 );
1464 }
1465 }
1466 } else if ( cmpResult < 0 ){
1467 overvalue = false; // our candidate is too small.
1468 diff = bigD.sub( bigB );
1469 } else {
1470 // the candidate is exactly right!
1471 // this happens with surprising fequency
1472 break correctionLoop;
1473 }
1474 FDBigInt halfUlp = constructPow52( B5, Ulp2 );
1475 if ( (cmpResult = diff.cmp( halfUlp ) ) < 0 ){
1476 // difference is small.
1477 // this is close enough
1478 if (mustSetRoundDir) {
1479 roundDir = overvalue ? -1 : 1;
1480 }
1481 break correctionLoop;
1482 } else if ( cmpResult == 0 ){
1483 // difference is exactly half an ULP
1484 // round to some other value maybe, then finish
1485 dValue += 0.5*ulp( dValue, overvalue );
1486 // should check for bigIntNBits == 1 here??
1487 if (mustSetRoundDir) {
1488 roundDir = overvalue ? -1 : 1;
1489 }
1490 break correctionLoop;
1491 } else {
1492 // difference is non-trivial.
1493 // could scale addend by ratio of difference to
1494 // halfUlp here, if we bothered to compute that difference.
1495 // Most of the time ( I hope ) it is about 1 anyway.
1496 dValue += ulp( dValue, overvalue );
1497 if ( dValue == 0.0 || dValue == Double.POSITIVE_INFINITY )
1498 break correctionLoop; // oops. Fell off end of range.
1499 continue; // try again.
1500 }
1501
1502 }
1503 return (isNegative)? -dValue : dValue;
1504 }
1505 }
1506
1507 /*
1508 * Take a FormattedFloatingDecimal, which we presumably just scanned in,
1509 * and find out what its value is, as a float.
1510 * This is distinct from doubleValue() to avoid the extremely
1511 * unlikely case of a double rounding error, wherein the converstion
1512 * to double has one rounding error, and the conversion of that double
1513 * to a float has another rounding error, IN THE WRONG DIRECTION,
1514 * ( because of the preference to a zero low-order bit ).
1515 */
1516
1517 public strictfp float floatValue(){
1518 int kDigits = Math.min( nDigits, singleMaxDecimalDigits+1 );
1519 int iValue;
1520 float fValue;
1521
1522 // First, check for NaN and Infinity values
1523 if(digits == infinity || digits == notANumber) {
1524 if(digits == notANumber)
1525 return Float.NaN;
1526 else
1527 return (isNegative?Float.NEGATIVE_INFINITY:Float.POSITIVE_INFINITY);
1528 }
1529 else {
1530 /*
1531 * convert the lead kDigits to an integer.
1532 */
1533 iValue = (int)digits[0]-(int)'0';
1534 for ( int i=1; i < kDigits; i++ ){
1535 iValue = iValue*10 + (int)digits[i]-(int)'0';
1536 }
1537 fValue = (float)iValue;
1538 int exp = decExponent-kDigits;
1539 /*
1540 * iValue now contains an integer with the value of
1541 * the first kDigits digits of the number.
1542 * fValue contains the (float) of the same.
1543 */
1544
1545 if ( nDigits <= singleMaxDecimalDigits ){
1546 /*
1547 * possibly an easy case.
1548 * We know that the digits can be represented
1549 * exactly. And if the exponent isn't too outrageous,
1550 * the whole thing can be done with one operation,
1551 * thus one rounding error.
1552 * Note that all our constructors trim all leading and
1553 * trailing zeros, so simple values (including zero)
1554 * will always end up here.
1555 */
1556 if (exp == 0 || fValue == 0.0f)
1557 return (isNegative)? -fValue : fValue; // small floating integer
1558 else if ( exp >= 0 ){
1559 if ( exp <= singleMaxSmallTen ){
1560 /*
1561 * Can get the answer with one operation,
1562 * thus one roundoff.
1563 */
1564 fValue *= singleSmall10pow[exp];
1565 return (isNegative)? -fValue : fValue;
1566 }
1567 int slop = singleMaxDecimalDigits - kDigits;
1568 if ( exp <= singleMaxSmallTen+slop ){
1569 /*
1570 * We can multiply dValue by 10^(slop)
1571 * and it is still "small" and exact.
1572 * Then we can multiply by 10^(exp-slop)
1573 * with one rounding.
1574 */
1575 fValue *= singleSmall10pow[slop];
1576 fValue *= singleSmall10pow[exp-slop];
1577 return (isNegative)? -fValue : fValue;
1578 }
1579 /*
1580 * Else we have a hard case with a positive exp.
1581 */
1582 } else {
1583 if ( exp >= -singleMaxSmallTen ){
1584 /*
1585 * Can get the answer in one division.
1586 */
1587 fValue /= singleSmall10pow[-exp];
1588 return (isNegative)? -fValue : fValue;
1589 }
1590 /*
1591 * Else we have a hard case with a negative exp.
1592 */
1593 }
1594 } else if ( (decExponent >= nDigits) && (nDigits+decExponent <= maxDecimalDigits) ){
1595 /*
1596 * In double-precision, this is an exact floating integer.
1597 * So we can compute to double, then shorten to float
1598 * with one round, and get the right answer.
1599 *
1600 * First, finish accumulating digits.
1601 * Then convert that integer to a double, multiply
1602 * by the appropriate power of ten, and convert to float.
1603 */
1604 long lValue = (long)iValue;
1605 for ( int i=kDigits; i < nDigits; i++ ){
1606 lValue = lValue*10L + (long)((int)digits[i]-(int)'0');
1607 }
1608 double dValue = (double)lValue;
1609 exp = decExponent-nDigits;
1610 dValue *= small10pow[exp];
1611 fValue = (float)dValue;
1612 return (isNegative)? -fValue : fValue;
1613
1614 }
1615 /*
1616 * Harder cases:
1617 * The sum of digits plus exponent is greater than
1618 * what we think we can do with one error.
1619 *
1620 * Start by weeding out obviously out-of-range
1621 * results, then convert to double and go to
1622 * common hard-case code.
1623 */
1624 if ( decExponent > singleMaxDecimalExponent+1 ){
1625 /*
1626 * Lets face it. This is going to be
1627 * Infinity. Cut to the chase.
1628 */
1629 return (isNegative)? Float.NEGATIVE_INFINITY : Float.POSITIVE_INFINITY;
1630 } else if ( decExponent < singleMinDecimalExponent-1 ){
1631 /*
1632 * Lets face it. This is going to be
1633 * zero. Cut to the chase.
1634 */
1635 return (isNegative)? -0.0f : 0.0f;
1636 }
1637
1638 /*
1639 * Here, we do 'way too much work, but throwing away
1640 * our partial results, and going and doing the whole
1641 * thing as double, then throwing away half the bits that computes
1642 * when we convert back to float.
1643 *
1644 * The alternative is to reproduce the whole multiple-precision
1645 * algorythm for float precision, or to try to parameterize it
1646 * for common usage. The former will take about 400 lines of code,
1647 * and the latter I tried without success. Thus the semi-hack
1648 * answer here.
1649 */
1650 mustSetRoundDir = !fromHex;
1651 double dValue = doubleValue();
1652 return stickyRound( dValue );
1653 }
1654 }
1655
1656
1657 /*
1658 * All the positive powers of 10 that can be
1659 * represented exactly in double/float.
1660 */
1661 private static final double small10pow[] = {
1662 1.0e0,
1663 1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5,
1664 1.0e6, 1.0e7, 1.0e8, 1.0e9, 1.0e10,
1665 1.0e11, 1.0e12, 1.0e13, 1.0e14, 1.0e15,
1666 1.0e16, 1.0e17, 1.0e18, 1.0e19, 1.0e20,
1667 1.0e21, 1.0e22
1668 };
1669
1670 private static final float singleSmall10pow[] = {
1671 1.0e0f,
1672 1.0e1f, 1.0e2f, 1.0e3f, 1.0e4f, 1.0e5f,
1673 1.0e6f, 1.0e7f, 1.0e8f, 1.0e9f, 1.0e10f
1674 };
1675
1676 private static final double big10pow[] = {
1677 1e16, 1e32, 1e64, 1e128, 1e256 };
1678 private static final double tiny10pow[] = {
1679 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1680
1681 private static final int maxSmallTen = small10pow.length-1;
1682 private static final int singleMaxSmallTen = singleSmall10pow.length-1;
1683
1684 private static final int small5pow[] = {
1685 1,
1686 5,
1687 5*5,
1688 5*5*5,
1689 5*5*5*5,
1690 5*5*5*5*5,
1691 5*5*5*5*5*5,
1692 5*5*5*5*5*5*5,
1693 5*5*5*5*5*5*5*5,
1694 5*5*5*5*5*5*5*5*5,
1695 5*5*5*5*5*5*5*5*5*5,
1696 5*5*5*5*5*5*5*5*5*5*5,
1697 5*5*5*5*5*5*5*5*5*5*5*5,
1698 5*5*5*5*5*5*5*5*5*5*5*5*5
1699 };
1700
1701
1702 private static final long long5pow[] = {
1703 1L,
1704 5L,
1705 5L*5,
1706 5L*5*5,
1707 5L*5*5*5,
1708 5L*5*5*5*5,
1709 5L*5*5*5*5*5,
1710 5L*5*5*5*5*5*5,
1711 5L*5*5*5*5*5*5*5,
1712 5L*5*5*5*5*5*5*5*5,
1713 5L*5*5*5*5*5*5*5*5*5,
1714 5L*5*5*5*5*5*5*5*5*5*5,
1715 5L*5*5*5*5*5*5*5*5*5*5*5,
1716 5L*5*5*5*5*5*5*5*5*5*5*5*5,
1717 5L*5*5*5*5*5*5*5*5*5*5*5*5*5,
1718 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1719 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1720 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1721 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1722 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1723 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1724 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1725 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1726 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1727 5L*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5*5,
1728 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,
1729 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,
1730 };
1731
1732 // approximately ceil( log2( long5pow[i] ) )
1733 private static final int n5bits[] = {
1734 0,
1735 3,
1736 5,
1737 7,
1738 10,
1739 12,
1740 14,
1741 17,
1742 19,
1743 21,
1744 24,
1745 26,
1746 28,
1747 31,
1748 33,
1749 35,
1750 38,
1751 40,
1752 42,
1753 45,
1754 47,
1755 49,
1756 52,
1757 54,
1758 56,
1759 59,
1760 61,
1761 };
1762
1763 private static final char infinity[] = { 'I', 'n', 'f', 'i', 'n', 'i', 't', 'y' };
1764 private static final char notANumber[] = { 'N', 'a', 'N' };
1765 private static final char zero[] = { '0', '0', '0', '0', '0', '0', '0', '0' };
1766 }