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