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