History

My first initial commit with all of mine experimental changes.

Dmitry Kabak (userdima committed 13 years ago

Commit d37b88d1 parent 3f5b46618352 files · +2598093 −-8690Browse files

Showing the changes to java/pkg/sun/misc/FloatingDecimal.java (1 of 18352 files). Show all

Unified · Split

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