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