1 | /* (C) Copyright 2001 William Rossi |
2 | */ |
3 | |
4 | package rossi.dfp; |
5 | |
6 | /**<pre> |
7 | * Decimal floating point library for Java |
8 | * |
9 | * Another floating point class. This one is built using radix 10000 |
10 | * which is 10^4, so its almost decimal. |
11 | * |
12 | * The design goals here are - |
13 | * 1.) decimal math, or close to it. |
14 | * 2.) Compile-time settable precision |
15 | * 3.) Portability. Code should be keep as portable as possible. |
16 | * 4.) Performance |
17 | * 5.) Accuracy - Results should always be +/- 1 ULP for basic |
18 | * algebraic operation |
19 | * 6.) Comply with IEEE 854-1987 as much as possible. |
20 | * (See IEEE 854-1987 notes below) |
21 | * |
22 | * The trade offs - |
23 | * 1.) Memory foot print. I'm using more memory than necessary to |
24 | * represent numbers to get better performance. |
25 | * 2.) Digits are bigger, so rounding is a greater loss. So, if you |
26 | * really need 12 decimal digits, better use 4 base 10000 digits |
27 | * there can be one partially filled. |
28 | * |
29 | * Numbers are represented in the following form: |
30 | * |
31 | * n = sign * mant * (radix) ^ exp; |
32 | * |
33 | * where sign is +/- 1, mantissa represents a fractional number between |
34 | * zero and one. mant[0] is the least significant digit. |
35 | * exp is in the range of -32767 to 32768 |
36 | * |
37 | * dfp objects are immuatable via their public interface. |
38 | * |
39 | * IEEE 854-1987 Notes and differences |
40 | * |
41 | * IEEE 854 requires the radix to be either 2 or 10. The radix here is |
42 | * 10000, so that requirement is not met, but it is possible that a |
43 | * subclassed can be made to make it behave as a radix 10 |
44 | * number. It is my opinion that if it looks and behaves as a radix |
45 | * 10 number then it is one and that requirement would be met. |
46 | * |
47 | * The radix of 10000 was chosen because it should be faster to operate |
48 | * on 4 decimal digits at once intead of one at a time. Radix 10 behaviour |
49 | * can be realized by add an additional rounding step to ensure that |
50 | * the number of decimal digits represented is constant. |
51 | * |
52 | * The IEEE standard specifically leaves out internal data encoding, |
53 | * so it is reasonable to conclude that such a subclass of this radix |
54 | * 10000 system is merely an encoding of a radix 10 system. |
55 | * |
56 | * IEEE 854 also specifies the existance of "sub-normal" numbers. This |
57 | * class does not contain any such entities. The most significant radix |
58 | * 10000 digit is always non-zero. Instead, we support "gradual underflow" |
59 | * by raising the underflow flag for numbers less with exponent less than |
60 | * expMin, but dont flush to zero until the exponent reaches expMin-DIGITS. |
61 | * Thus the smallest number we can represent would be: |
62 | * 1E(-(minExp-DIGITS-1)*4), eg, for DIGITS=5, minExp=-32767, that would |
63 | * be 1e-131092. |
64 | * |
65 | * IEEE 854 defines that the implied radix point lies just to the right |
66 | * of the most significant digit and to the left of the remaining digits. |
67 | * This implementation puts the implied radix point to the left of all |
68 | * digits including the most significant one. The most significant digit |
69 | * here is the one just to the right of the radix point. This is a fine |
70 | * detail and is really only a matter of definition. Any side effects of |
71 | * this can be rendered invisible by a subclass. |
72 | * </pre> |
73 | */ |
74 | |
75 | public class dfp |
76 | { |
77 | protected int[] mant; // the mantissa |
78 | protected byte sign; // the sign bit. 1 for positive, -1 for negative |
79 | protected int exp; // the exponent. |
80 | protected byte nans; // Indicates non-finite / non-number values |
81 | |
82 | protected static int rMode = 4; // Current rounding mode |
83 | protected static int ieeeFlags = 0; // IEEE 854-1987 signals |
84 | |
85 | /** The number of digits. note these are radix 10000 digits, so each |
86 | * one is equivilent to 4 decimal digits */ |
87 | public final static int DIGITS = 5; // each digit yeilds 4 decimal digits |
88 | |
89 | /** The radix, or base of this system. Set to 10000 */ |
90 | public final static int radix = 10000; |
91 | |
92 | /** The minium exponent before underflow is signaled. Flush to zero |
93 | * occurs at minExp-DIGITS */ |
94 | public final static int minExp = -32767; |
95 | |
96 | /** The maximum exponent before overflow is signaled and results flushed |
97 | * to infinity */ |
98 | public final static int maxExp = 32768; |
99 | |
100 | /** The amount under/overflows are scaled by before going to trap handler */ |
101 | public final static int errScale = 32760; |
102 | |
103 | /**** Rounding modes *****/ |
104 | |
105 | /** Rounds toward zero. I.E. truncation */ |
106 | public final static int ROUND_DOWN = 0; |
107 | |
108 | /** Rounds away from zero if discarded digit is non-zero */ |
109 | public final static int ROUND_UP = 1; |
110 | |
111 | /** Rounds towards nearest unless both are equidistant in which case |
112 | * it rounds away from zero */ |
113 | public final static int ROUND_HALF_UP = 2; |
114 | |
115 | /** Rounds towards nearest unless both are equidistant in which case |
116 | * it rounds toward zero */ |
117 | public final static int ROUND_HALF_DOWN = 3; |
118 | |
119 | /** Rounds towards nearest unless both are equidistant in which case |
120 | * it rounds toward the even neighbor. This is the default as |
121 | * specified by IEEE 854-1987 */ |
122 | public final static int ROUND_HALF_EVEN = 4; |
123 | |
124 | /** Rounds towards nearest unless both are equidistant in which case |
125 | * it rounds toward the odd neighbor. */ |
126 | public final static int ROUND_HALF_ODD = 5; |
127 | |
128 | /** Rounds towards positive infinity */ |
129 | public final static int ROUND_CEIL = 6; |
130 | |
131 | /** Rounds towards negative infinity */ |
132 | public final static int ROUND_FLOOR = 7; |
133 | |
134 | /* non-numbers per IEEE 854-1987 */ |
135 | public final static byte FINITE = 0; // Normal finite numbers |
136 | public final static byte INFINITE = 1; // Infinity |
137 | public final static byte SNAN = 2; // Signaling NaN |
138 | public final static byte QNAN = 3; // Quiet NaN |
139 | |
140 | /* Flags */ |
141 | public final static int FLAG_INVALID = 1; // Invalid operation |
142 | public final static int FLAG_DIV_ZERO = 2; // Division by zero |
143 | public final static int FLAG_OVERFLOW = 4; // Overflow |
144 | public final static int FLAG_UNDERFLOW = 8; // Underflow |
145 | public final static int FLAG_INEXACT = 16; // Inexact |
146 | |
147 | /* Handy constants */ |
148 | public final static dfp zero = new dfp(); |
149 | public final static dfp one = new dfp("1"); |
150 | public final static dfp two = new dfp("2"); |
151 | |
152 | /** Default constructor. Makes a dfp with a value of zero */ |
153 | public dfp() |
154 | { |
155 | mant = new int[DIGITS]; |
156 | for (int i=DIGITS-1; i>=0; i--) |
157 | mant[i] = 0; |
158 | sign = 1; |
159 | exp = 0; |
160 | nans = FINITE; |
161 | } |
162 | |
163 | /** Copy constructor. Creates a copy of the supplied dfp */ |
164 | public dfp(dfp d) |
165 | { |
166 | mant = new int[DIGITS]; |
167 | |
168 | for (int i=DIGITS-1; i>=0; i--) |
169 | mant[i] = d.mant[i]; |
170 | sign = d.sign; |
171 | exp = d.exp; |
172 | nans = d.nans; |
173 | } |
174 | |
175 | /** Create a dfp given a String representation */ |
176 | public dfp(String s) |
177 | { |
178 | dfp r = string2dfp(s); |
179 | this.mant = r.mant; |
180 | this.exp = r.exp; |
181 | this.sign = r.sign; |
182 | this.nans = r.nans; |
183 | } |
184 | |
185 | /** Create a dfp. Use this internally in preferenct to constructors to facilitate |
186 | subclasses. */ |
187 | public dfp newInstance(dfp d) |
188 | { |
189 | return new dfp(d); |
190 | } |
191 | |
192 | /** Create a dfp. Use this internally in preferenct to constructors to facilitate |
193 | subclasses. */ |
194 | public dfp newInstance(String s) |
195 | { |
196 | return new dfp(s); |
197 | } |
198 | |
199 | /** Shift the mantissa left, and adjust the exponent to compensate */ |
200 | protected void shiftLeft() |
201 | { |
202 | for (int i=DIGITS-1; i>0; i--) |
203 | mant[i] = mant[i-1]; |
204 | mant[0] = 0; |
205 | exp--; |
206 | } |
207 | |
208 | /* Note that shiftRight() does not call round() as that round() itself |
209 | uses shiftRight() */ |
210 | /** Shift the mantissa right, and adjust the exponent to compensate */ |
211 | protected void shiftRight() |
212 | { |
213 | for (int i=0; i<DIGITS-1; i++) |
214 | mant[i] = mant[i+1]; |
215 | mant[DIGITS-1] = 0; |
216 | exp++; |
217 | } |
218 | |
219 | /** Make our exp equal to the supplied one. This may cause rounding. |
220 | * Also causes de-normalized numbers. These numbers are generally |
221 | * dangerous because most routines assume normalized numbers. |
222 | * Align doesn't round, so it will return the last digit destroyed |
223 | * by shifting right. |
224 | */ |
225 | protected int align(int e) |
226 | { |
227 | int diff; |
228 | int adiff; |
229 | int lostdigit = 0; |
230 | boolean inexact = false; |
231 | |
232 | diff = exp - e; |
233 | |
234 | adiff = diff; |
235 | if (adiff < 0) |
236 | adiff = -adiff; |
237 | |
238 | if (diff == 0) |
239 | return 0; |
240 | |
241 | if (adiff > (DIGITS+1)) // Special case |
242 | { |
243 | for (int i=DIGITS-1; i>=0; i--) |
244 | mant[i] = 0; |
245 | exp = e; |
246 | |
247 | ieeeFlags |= FLAG_INEXACT; |
248 | dotrap(FLAG_INEXACT, "align", this, this); |
249 | |
250 | return 0; |
251 | } |
252 | |
253 | for (int i=0; i<adiff; i++) |
254 | { |
255 | if (diff < 0) |
256 | { |
257 | /* Keep track of loss -- only signal inexact after losing 2 digits. |
258 | * the first lost digit is returned to add() and may be incorporated |
259 | * into the result. |
260 | */ |
261 | if (lostdigit != 0) |
262 | inexact = true; |
263 | |
264 | lostdigit = mant[0]; |
265 | |
266 | shiftRight(); |
267 | } |
268 | else |
269 | shiftLeft(); |
270 | } |
271 | |
272 | if (inexact) |
273 | { |
274 | ieeeFlags |= FLAG_INEXACT; |
275 | dotrap(FLAG_INEXACT, "align", this, this); |
276 | } |
277 | |
278 | return lostdigit; |
279 | } |
280 | |
281 | /** returns true if this is less than x. |
282 | * returns false if this or x is NaN */ |
283 | public boolean lessThan(dfp x) |
284 | { |
285 | /* if a nan is involved, signal invalid and return false */ |
286 | if (nans == SNAN || nans == QNAN || x.nans == SNAN || x.nans == QNAN) |
287 | { |
288 | ieeeFlags |= FLAG_INVALID; |
289 | dotrap(FLAG_INVALID, "lessThan", x, newInstance(zero)); |
290 | return false; |
291 | } |
292 | |
293 | return (compare(this, x) < 0); |
294 | } |
295 | |
296 | /** returns true if this is greater than x. |
297 | * returns false if this or x is NaN */ |
298 | public boolean greaterThan(dfp x) |
299 | { |
300 | /* if a nan is involved, signal invalid and return false */ |
301 | if (nans == SNAN || nans == QNAN || x.nans == SNAN || x.nans == QNAN) |
302 | { |
303 | ieeeFlags |= FLAG_INVALID; |
304 | dotrap(FLAG_INVALID, "lessThan", x, newInstance(zero)); |
305 | return false; |
306 | } |
307 | |
308 | return (compare(this, x) > 0); |
309 | } |
310 | |
311 | /** returns true if this is equal to x. |
312 | * returns false if this or x is NaN */ |
313 | public boolean equal(dfp x) |
314 | { |
315 | if (nans == SNAN || nans == QNAN || x.nans == SNAN || x.nans == QNAN) |
316 | return false; |
317 | |
318 | return (compare(this, x) == 0); |
319 | } |
320 | |
321 | /** returns true if this is not equal to x. |
322 | * different from !equal(x) in the way NaNs are handled. |
323 | */ |
324 | public boolean unequal(dfp x) |
325 | { |
326 | if (nans == SNAN || nans == QNAN || x.nans == SNAN || x.nans == QNAN) |
327 | return false; |
328 | |
329 | return (greaterThan(x) || lessThan(x)); |
330 | } |
331 | |
332 | /** compare a and b. return -1 if a<b, 1 if a>b and 0 if a==b |
333 | * Note this method does not properly handle NaNs. */ |
334 | protected static int compare(dfp a, dfp b) |
335 | { |
336 | /* Ignore the sign of zero */ |
337 | if (a.mant[DIGITS-1] == 0 && b.mant[DIGITS-1] == 0 && a.nans == FINITE && b.nans == FINITE) |
338 | return 0; |
339 | |
340 | if (a.sign != b.sign) |
341 | { |
342 | if (a.sign == -1) |
343 | return -1; |
344 | else |
345 | return 1; |
346 | } |
347 | |
348 | /* deal with the infinities */ |
349 | if (a.nans == INFINITE && b.nans == FINITE) |
350 | return a.sign; |
351 | |
352 | if (a.nans == FINITE && b.nans == INFINITE) |
353 | return -b.sign; |
354 | |
355 | if (a.nans == INFINITE && b.nans == INFINITE) |
356 | return 0; |
357 | |
358 | /* Handle special case when a or b is zero, by ignoring the exponents */ |
359 | if (b.mant[DIGITS-1] != 0 && a.mant[DIGITS-1] != 0) |
360 | { |
361 | if (a.exp < b.exp) |
362 | return -a.sign; |
363 | |
364 | if (a.exp > b.exp) |
365 | return a.sign; |
366 | } |
367 | |
368 | /* compare the mantissas */ |
369 | for (int i=DIGITS-1; i>=0; i--) |
370 | { |
371 | if (a.mant[i] > b.mant[i]) |
372 | return a.sign; |
373 | |
374 | if (a.mant[i] < b.mant[i]) |
375 | return -a.sign; |
376 | } |
377 | |
378 | return 0; |
379 | } |
380 | |
381 | /** Round to nearest integer using the round-half-even method. |
382 | * That is round to nearest integer unless both are equidistant. |
383 | * In which case round to the even one. |
384 | */ |
385 | public dfp rint() |
386 | { |
387 | return trunc(ROUND_HALF_EVEN); |
388 | } |
389 | |
390 | /** Round to an integer using the round floor mode. That is, |
391 | * round toward -Infinity */ |
392 | public dfp floor() |
393 | { |
394 | return trunc(ROUND_FLOOR); |
395 | } |
396 | |
397 | /** Round to an integer using the ceil floor mode. That is, |
398 | * round toward +Infinity */ |
399 | public dfp ceil() |
400 | { |
401 | return trunc(ROUND_CEIL); |
402 | } |
403 | |
404 | /** Returns the IEEE remainder. That is the result of this less |
405 | * n times d, where n is the integer closest to this/d. |
406 | */ |
407 | public dfp remainder(dfp d) |
408 | { |
409 | dfp q = this.divide(d); |
410 | dfp result; |
411 | |
412 | result = this.subtract(this.divide(d).rint().multiply(d)); |
413 | |
414 | /* IEEE 854-1987 says that if the result is zero, then it |
415 | carries the sign of this */ |
416 | |
417 | if (result.mant[DIGITS-1] == 0) |
418 | result.sign = sign; |
419 | |
420 | return result; |
421 | } |
422 | |
423 | /** Does the integer conversions with the spec rounding. */ |
424 | protected dfp trunc(int rmode) |
425 | { |
426 | dfp result, a, half; |
427 | boolean changed = false; |
428 | |
429 | if (nans == SNAN || nans == QNAN) |
430 | return newInstance(this); |
431 | |
432 | if (nans == INFINITE) |
433 | return newInstance(this); |
434 | |
435 | if (mant[DIGITS-1] == 0) // a is zero |
436 | return newInstance(this); |
437 | |
438 | /* If the exponent is less than zero then we can certainly |
439 | * return zero */ |
440 | if (exp < 0) |
441 | { |
442 | ieeeFlags |= FLAG_INEXACT; |
443 | result = newInstance(zero); |
444 | result = dotrap(FLAG_INEXACT, "trunc", this, result); |
445 | return result; |
446 | } |
447 | |
448 | /* If the exponent is greater than or equal to digits, then it |
449 | * must already be an integer since there is no precision left |
450 | * for any fractional part */ |
451 | |
452 | if (exp >= DIGITS) |
453 | return newInstance(this); |
454 | |
455 | /* General case: create another dfp, result, that contains the |
456 | * a with the fractional part lopped off. */ |
457 | |
458 | result = newInstance(this); |
459 | for (int i=0; i<(DIGITS-result.exp); i++) |
460 | { |
461 | changed |= (result.mant[i] != 0); |
462 | result.mant[i] = 0; |
463 | } |
464 | |
465 | if (changed) |
466 | { |
467 | switch (rmode) |
468 | { |
469 | case ROUND_FLOOR: |
470 | if (result.sign == -1) // then we must increment the mantissa by one |
471 | result = result.add(newInstance("-1")); |
472 | break; |
473 | |
474 | case ROUND_CEIL: |
475 | if (result.sign == 1) // then we must increment the mantissa by one |
476 | result = result.add(one); |
477 | break; |
478 | |
479 | case ROUND_HALF_EVEN: |
480 | default: |
481 | half = newInstance("0.5"); |
482 | a = subtract(result); // difference between this and result |
483 | a.sign = 1; // force positive (take abs) |
484 | if (a.greaterThan(half)) |
485 | { |
486 | a = newInstance(one); |
487 | a.sign = sign; |
488 | result = result.add(a); |
489 | } |
490 | |
491 | /** If exactly equal to 1/2 and odd then increment */ |
492 | if (a.equal(half) && result.exp > 0 && (result.mant[DIGITS-result.exp]&1) != 0) |
493 | { |
494 | a = newInstance(one); |
495 | a.sign = sign; |
496 | result = result.add(a); |
497 | } |
498 | break; |
499 | } |
500 | |
501 | ieeeFlags |= FLAG_INEXACT; // signal inexact |
502 | result = dotrap(FLAG_INEXACT, "trunc", this, result); |
503 | return result; |
504 | } |
505 | |
506 | return result; |
507 | } |
508 | |
509 | /** Convert this to an integer. If greater than 2147483647, it returns |
510 | * 2147483647. If less than -2147483648 it returns -2147483648. |
511 | */ |
512 | public int intValue() |
513 | { |
514 | dfp rounded; |
515 | int result = 0; |
516 | |
517 | rounded = rint(); |
518 | |
519 | if (rounded.greaterThan(newInstance("2147483647"))) |
520 | return 2147483647; |
521 | |
522 | if (rounded.lessThan(newInstance("-2147483648"))) |
523 | return -2147483648; |
524 | |
525 | for (int i=DIGITS-1; i>=DIGITS-rounded.exp; i--) |
526 | result = result*radix+rounded.mant[i]; |
527 | |
528 | if (rounded.sign == -1) |
529 | result = -result; |
530 | |
531 | return result; |
532 | } |
533 | |
534 | /** Returns the exponent of the greatest power of 10000 that is |
535 | * less than or equal to the absolute value of this. I.E. if |
536 | * this is 10e6 then log10K would return 1. |
537 | */ |
538 | public int log10K() |
539 | { |
540 | return exp-1; |
541 | } |
542 | |
543 | /** Return the specified power of 10000 */ |
544 | public dfp power10K(int e) |
545 | { |
546 | dfp d = newInstance(one); |
547 | d.exp = e + 1; |
548 | return d; |
549 | } |
550 | |
551 | /** |
552 | * Return the exponent of the greatest power of 10 that is less than |
553 | * or equal to than abs(this). |
554 | */ |
555 | public int log10() |
556 | { |
557 | if (mant[DIGITS-1] > 1000) return exp * 4 - 1; |
558 | if (mant[DIGITS-1] > 100) return exp * 4 - 2; |
559 | if (mant[DIGITS-1] > 10) return exp * 4 - 3; |
560 | return exp * 4 - 4; |
561 | } |
562 | |
563 | /** Return the specified power of 10 */ |
564 | public dfp power10(int e) |
565 | { |
566 | dfp d = newInstance(one); |
567 | |
568 | if (e >= 0) |
569 | d.exp = e/4 + 1; |
570 | else |
571 | d.exp = (e+1)/4; |
572 | |
573 | switch ((e%4+4)%4) |
574 | { |
575 | case 0: break; |
576 | case 1: d = d.multiply(10); break; |
577 | case 2: d = d.multiply(100); break; |
578 | case 3: d = d.multiply(1000); break; |
579 | } |
580 | |
581 | return d; |
582 | } |
583 | |
584 | /** Negate the mantissa of this by computing the complement. |
585 | * Leaves the sign bit unchanged, used internally by add. |
586 | * Denormalized numbers are handled properly here. */ |
587 | protected int complement(int extra) |
588 | { |
589 | int r, rl, rh; |
590 | |
591 | extra = radix-extra; |
592 | for (int i=0; i<DIGITS; i++) |
593 | mant[i] = radix-mant[i]-1; |
594 | |
595 | rh = extra / radix; |
596 | extra = extra % radix; |
597 | for (int i=0; i<DIGITS; i++) |
598 | { |
599 | r = mant[i]+rh; |
600 | rl = r % radix; |
601 | rh = r / radix; |
602 | mant[i] = rl; |
603 | } |
604 | |
605 | return extra; |
606 | } |
607 | |
608 | /** Add x to this and return the result */ |
609 | public dfp add(dfp x) |
610 | { |
611 | int r, rh, rl; |
612 | dfp a, b, result; |
613 | byte asign, bsign, rsign; |
614 | int aextradigit = 0, bextradigit = 0; |
615 | |
616 | /* handle special cases */ |
617 | if (nans != FINITE || x.nans != FINITE) |
618 | { |
619 | if (nans == QNAN || nans == SNAN) |
620 | return this; |
621 | |
622 | if (x.nans == QNAN || x.nans == SNAN) |
623 | return x; |
624 | |
625 | if (nans == INFINITE && x.nans == FINITE) |
626 | return this; |
627 | |
628 | if (x.nans == INFINITE && nans == FINITE) |
629 | return x; |
630 | |
631 | if (x.nans == INFINITE && nans == INFINITE && sign == x.sign) |
632 | return x; |
633 | |
634 | if (x.nans == INFINITE && nans == INFINITE && sign != x.sign) |
635 | { |
636 | ieeeFlags |= FLAG_INVALID; |
637 | result = newInstance(zero); |
638 | result.nans = QNAN; |
639 | result = dotrap(FLAG_INVALID, "add", x, result); |
640 | return result; |
641 | } |
642 | } |
643 | |
644 | /* copy this and the arg */ |
645 | a = newInstance(this); |
646 | b = newInstance(x); |
647 | |
648 | /* initialize the result object */ |
649 | result = newInstance(zero); |
650 | |
651 | /* Make all numbers positive, but remember their sign */ |
652 | asign = a.sign; |
653 | bsign = b.sign; |
654 | |
655 | a.sign = 1; |
656 | b.sign = 1; |
657 | |
658 | /* The result will be signed like the arg with greatest magnitude */ |
659 | rsign = bsign; |
660 | if (compare(a, b) > 0) |
661 | rsign = asign; |
662 | |
663 | /* Handle special case when a or b is zero, by setting the exponent |
664 | of the zero number equal to the other one. This avoids an alignment |
665 | which would cause catastropic loss of precision */ |
666 | if (b.mant[DIGITS-1] == 0) |
667 | b.exp = a.exp; |
668 | |
669 | if (a.mant[DIGITS-1] == 0) |
670 | a.exp = b.exp; |
671 | |
672 | /* align number with the smaller exponent */ |
673 | if (a.exp < b.exp) |
674 | aextradigit = a.align(b.exp); |
675 | else |
676 | bextradigit = b.align(a.exp); |
677 | |
678 | /* complement the smaller of the two if the signs are different */ |
679 | if (asign != bsign) |
680 | { |
681 | if (asign == rsign) |
682 | bextradigit = b.complement(bextradigit); |
683 | else |
684 | aextradigit = a.complement(aextradigit); |
685 | } |
686 | |
687 | /* add the mantissas */ |
688 | rh = 0; /* acts as a carry */ |
689 | for (int i=0; i<DIGITS; i++) |
690 | { |
691 | r = a.mant[i]+b.mant[i]+rh; |
692 | rl = r % radix; |
693 | rh = r / radix; |
694 | result.mant[i] = rl; |
695 | } |
696 | result.exp = a.exp; |
697 | result.sign = rsign; |
698 | |
699 | //System.out.println("a = "+a.mant[2]+" "+a.mant[1]+" "+a.mant[0]); |
700 | //System.out.println("b = "+b.mant[2]+" "+b.mant[1]+" "+b.mant[0]); |
701 | //System.out.println("result = "+result.mant[2]+" "+result.mant[1]+" "+result.mant[0]); |
702 | |
703 | /* handle overflow -- note, when asign!=bsign an overflow is |
704 | * normal and should be ignored. */ |
705 | |
706 | if (rh != 0 && (asign == bsign)) |
707 | { |
708 | int lostdigit = result.mant[0]; |
709 | result.shiftRight(); |
710 | result.mant[DIGITS-1] = rh; |
711 | int excp = result.round(lostdigit); |
712 | if (excp != 0) |
713 | result = dotrap(excp, "add", x, result); |
714 | } |
715 | |
716 | /* normalize the result */ |
717 | for (int i=0; i<DIGITS; i++) |
718 | { |
719 | if (result.mant[DIGITS-1] != 0) |
720 | break; |
721 | result.shiftLeft(); |
722 | if (i == 0) |
723 | { |
724 | result.mant[0] = aextradigit+bextradigit; |
725 | aextradigit = 0; |
726 | bextradigit = 0; |
727 | } |
728 | } |
729 | |
730 | /* result is zero if after normalization the most sig. digit is zero */ |
731 | if (result.mant[DIGITS-1] == 0) |
732 | { |
733 | result.exp = 0; |
734 | |
735 | if (asign != bsign) // Unless adding 2 negative zeros, sign is positive |
736 | result.sign = 1; // Per IEEE 854-1987 Section 6.3 |
737 | } |
738 | |
739 | /* Call round to test for over/under flows */ |
740 | int excp = result.round((aextradigit+bextradigit)); |
741 | if (excp != 0) |
742 | result = dotrap(excp, "add", x, result); |
743 | |
744 | return result; |
745 | } |
746 | |
747 | /** Returns a number that is this number with the sign bit reversed */ |
748 | public dfp negate() |
749 | { |
750 | dfp result = newInstance(this); |
751 | result.sign = (byte) -result.sign; |
752 | return result; |
753 | } |
754 | |
755 | /** Subtract a from this */ |
756 | public dfp subtract(dfp a) |
757 | { |
758 | return add(a.negate()); |
759 | } |
760 | |
761 | /** round this given the next digit n using the current rounding mode |
762 | * returns a flag if an exception occured |
763 | */ |
764 | protected int round(int n) |
765 | { |
766 | int r, rh, rl; |
767 | boolean inc=false; |
768 | |
769 | switch (rMode) |
770 | { |
771 | case ROUND_DOWN: |
772 | inc = false; |
773 | break; |
774 | |
775 | case ROUND_UP: |
776 | inc = (n!=0); // round up if n!=0 |
777 | break; |
778 | |
779 | case ROUND_HALF_UP: |
780 | inc = (n >= 5000); // round half up |
781 | break; |
782 | |
783 | case ROUND_HALF_DOWN: |
784 | inc = (n > 5000); // round half down |
785 | break; |
786 | |
787 | case ROUND_HALF_EVEN: |
788 | inc = (n > 5000 || (n == 5000 && (mant[0]&1)==1)); // round half-even |
789 | break; |
790 | |
791 | case ROUND_HALF_ODD: |
792 | inc = (n > 5000 || (n == 5000 && (mant[0]&1)==0)); // round half-odd |
793 | break; |
794 | |
795 | case ROUND_CEIL: |
796 | inc = (sign == 1 && n != 0); // round ceil |
797 | break; |
798 | |
799 | case ROUND_FLOOR: |
800 | inc = (sign == -1 && n != 0); // round floor |
801 | break; |
802 | } |
803 | |
804 | if (inc) // increment if necessary |
805 | { |
806 | rh = 1; |
807 | for (int i=0; i<DIGITS; i++) |
808 | { |
809 | r = mant[i] + rh; |
810 | rh = r / radix; |
811 | rl = r % radix; |
812 | mant[i] = rl; |
813 | } |
814 | |
815 | if (rh != 0) |
816 | { |
817 | shiftRight(); |
818 | mant[DIGITS-1]=rh; |
819 | } |
820 | } |
821 | |
822 | /* Check for exceptional cases and raise signals if necessary */ |
823 | if (exp < minExp) // Gradual Underflow |
824 | { |
825 | ieeeFlags |= FLAG_UNDERFLOW; |
826 | return FLAG_UNDERFLOW; |
827 | } |
828 | |
829 | if (exp > maxExp) // Overflow |
830 | { |
831 | ieeeFlags |= FLAG_OVERFLOW; |
832 | return FLAG_OVERFLOW; |
833 | } |
834 | |
835 | if (n != 0) // Inexact |
836 | { |
837 | ieeeFlags |= FLAG_INEXACT; |
838 | return FLAG_INEXACT; |
839 | } |
840 | return 0; |
841 | } |
842 | |
843 | /** Multiply this by x */ |
844 | public dfp multiply(dfp x) |
845 | { |
846 | int product[]; // product array |
847 | int r, rh, rl; // working registers |
848 | int md=0; // most sig digit in result |
849 | int excp; // exception code if any |
850 | dfp result = newInstance(zero); |
851 | |
852 | /* handle special cases */ |
853 | if (nans != FINITE || x.nans != FINITE) |
854 | { |
855 | if (nans == QNAN || nans == SNAN) |
856 | return this; |
857 | |
858 | if (x.nans == QNAN || x.nans == SNAN) |
859 | return x; |
860 | |
861 | if (nans == INFINITE && x.nans == FINITE && x.mant[DIGITS-1] != 0) |
862 | { |
863 | result = newInstance(this); |
864 | result.sign = (byte) (sign * x.sign); |
865 | return result; |
866 | } |
867 | |
868 | if (x.nans == INFINITE && nans == FINITE && mant[DIGITS-1] != 0) |
869 | { |
870 | result = newInstance(x); |
871 | result.sign = (byte) (sign * x.sign); |
872 | return result; |
873 | } |
874 | |
875 | if (x.nans == INFINITE && nans == INFINITE) |
876 | { |
877 | result = newInstance(this); |
878 | result.sign = (byte) (sign * x.sign); |
879 | return result; |
880 | } |
881 | |
882 | if ( (x.nans == INFINITE && nans == FINITE && mant[DIGITS-1] == 0) || |
883 | (nans == INFINITE && x.nans == FINITE && x.mant[DIGITS-1] == 0) ) |
884 | { |
885 | ieeeFlags |= FLAG_INVALID; |
886 | result = newInstance(zero); |
887 | result.nans = QNAN; |
888 | result = dotrap(FLAG_INVALID, "multiply", x, result); |
889 | return result; |
890 | } |
891 | } |
892 | |
893 | product = new int[DIGITS*2]; // Big enough to hold even the largest result |
894 | |
895 | for (int i=0; i<DIGITS*2; i++) |
896 | product[i] = 0; |
897 | |
898 | for (int i=0; i<DIGITS; i++) |
899 | { |
900 | rh = 0; // acts as a carry |
901 | for (int j=0; j<DIGITS; j++) |
902 | { |
903 | r = mant[i] * x.mant[j]; // multiply the 2 digits |
904 | r = r + product[i+j] + rh; // add to the product digit with carry in |
905 | rl = r % radix; |
906 | rh = r / radix; |
907 | product[i+j] = rl; |
908 | } |
909 | product[i+DIGITS] = rh; |
910 | } |
911 | |
912 | /* Find the most sig digit */ |
913 | md = DIGITS*2-1; // default, in case result is zero |
914 | for (int i=DIGITS*2-1; i>=0; i--) |
915 | { |
916 | if (product[i] != 0) |
917 | { |
918 | md = i; |
919 | break; |
920 | } |
921 | } |
922 | |
923 | /* Copy the digits into the result */ |
924 | for (int i=0; i<DIGITS; i++) |
925 | result.mant[DIGITS-i-1] = product[md-i]; |
926 | |
927 | /* Fixup the exponent. */ |
928 | result.exp = (exp + x.exp + md - 2*DIGITS + 1); |
929 | result.sign = (byte)((sign == x.sign)?1:-1); |
930 | |
931 | if (result.mant[DIGITS-1] == 0) // if result is zero, set exp to zero |
932 | result.exp = 0; |
933 | |
934 | if (md > (DIGITS-1)) |
935 | excp = result.round(product[md-DIGITS]); |
936 | else |
937 | excp = result.round(0); // has no effect except to check status |
938 | |
939 | if (excp != 0) |
940 | result = dotrap(excp, "multiply", x, result); |
941 | |
942 | return result; |
943 | } |
944 | |
945 | /** Multiply this by a single digit 0<=x<radix. There are speed |
946 | * advantages in this special case */ |
947 | public dfp multiply(int x) |
948 | { |
949 | int product[]; // product array |
950 | int r, rh, rl; // working registers |
951 | int excp; // exception code if any |
952 | int lostdigit; // rounded off digit |
953 | dfp result = newInstance(this); |
954 | |
955 | /* handle special cases */ |
956 | if (nans != FINITE) |
957 | { |
958 | if (nans == QNAN || nans == SNAN) |
959 | return this; |
960 | |
961 | if (nans == INFINITE && x != 0) |
962 | { |
963 | result = newInstance(this); |
964 | return result; |
965 | } |
966 | |
967 | if (nans == INFINITE && x == 0) |
968 | { |
969 | ieeeFlags |= FLAG_INVALID; |
970 | result = newInstance(zero); |
971 | result.nans = QNAN; |
972 | result = dotrap(FLAG_INVALID, "multiply", newInstance(zero), result); |
973 | return result; |
974 | } |
975 | } |
976 | |
977 | /* range check x */ |
978 | if (x < 0 || x >= radix) |
979 | { |
980 | ieeeFlags |= FLAG_INVALID; |
981 | result = newInstance(zero); |
982 | result.nans = QNAN; |
983 | result = dotrap(FLAG_INVALID, "multiply", result, result); |
984 | return result; |
985 | } |
986 | |
987 | rh = 0; |
988 | for (int i=0; i<DIGITS; i++) |
989 | { |
990 | r = mant[i] * x + rh; |
991 | rl = r % radix; |
992 | rh = r / radix; |
993 | result.mant[i] = rl; |
994 | } |
995 | |
996 | lostdigit = 0; |
997 | if (rh != 0) |
998 | { |
999 | lostdigit = result.mant[0]; |
1000 | result.shiftRight(); |
1001 | result.mant[DIGITS-1] = rh; |
1002 | } |
1003 | |
1004 | if (result.mant[DIGITS-1] == 0) // if result is zero, set exp to zero |
1005 | result.exp = 0; |
1006 | |
1007 | excp = result.round(lostdigit); |
1008 | |
1009 | if (excp != 0) |
1010 | result = dotrap(excp, "multiply", result, result); |
1011 | |
1012 | return result; |
1013 | } |
1014 | |
1015 | /** Divide this by divisor */ |
1016 | public dfp divide(dfp divisor) |
1017 | { |
1018 | int dividend[]; // current status of the dividend |
1019 | int quotient[]; // quotient |
1020 | int remainder[];// remainder |
1021 | int qd; // current quotient digit we're working with |
1022 | int nsqd; // number of significant quotient digits we have |
1023 | int trial=0; // trial quotient digit |
1024 | int min, max; // quotient digit limits |
1025 | int minadj; // minimum adjustment |
1026 | boolean trialgood; // Flag to indicate a good trail digit |
1027 | int r, rh, rl; // working registers |
1028 | int md=0; // most sig digit in result |
1029 | int excp; // exceptions |
1030 | dfp result = newInstance(zero); |
1031 | |
1032 | /* handle special cases */ |
1033 | if (nans != FINITE || divisor.nans != FINITE) |
1034 | { |
1035 | if (nans == QNAN || nans == SNAN) |
1036 | return this; |
1037 | |
1038 | if (divisor.nans == QNAN || divisor.nans == SNAN) |
1039 | return divisor; |
1040 | |
1041 | if (nans == INFINITE && divisor.nans == FINITE) |
1042 | { |
1043 | result = newInstance(this); |
1044 | result.sign = (byte) (sign * divisor.sign); |
1045 | return result; |
1046 | } |
1047 | |
1048 | if (divisor.nans == INFINITE && nans == FINITE) |
1049 | { |
1050 | result = newInstance(zero); |
1051 | result.sign = (byte) (sign * divisor.sign); |
1052 | return result; |
1053 | } |
1054 | |
1055 | if (divisor.nans == INFINITE && nans == INFINITE) |
1056 | { |
1057 | ieeeFlags |= FLAG_INVALID; |
1058 | result = newInstance(zero); |
1059 | result.nans = QNAN; |
1060 | result = dotrap(FLAG_INVALID, "divide", divisor, result); |
1061 | return result; |
1062 | } |
1063 | } |
1064 | |
1065 | /* Test for divide by zero */ |
1066 | if (divisor.mant[DIGITS-1] == 0) |
1067 | { |
1068 | ieeeFlags |= FLAG_DIV_ZERO; |
1069 | result = newInstance(zero); |
1070 | result.sign = (byte) (sign * divisor.sign); |
1071 | result.nans = INFINITE; |
1072 | result = dotrap(FLAG_DIV_ZERO, "divide", divisor, result); |
1073 | return result; |
1074 | } |
1075 | |
1076 | dividend = new int[DIGITS+1]; // one extra digit needed |
1077 | quotient = new int[DIGITS+2]; // two extra digits needed 1 for overflow, 1 for rounding |
1078 | remainder = new int[DIGITS+1]; // one extra digit needed |
1079 | |
1080 | /* Initialize our most significat digits to zero */ |
1081 | |
1082 | dividend[DIGITS] = 0; |
1083 | quotient[DIGITS] = 0; |
1084 | quotient[DIGITS+1] = 0; |
1085 | remainder[DIGITS] = 0; |
1086 | |
1087 | /* copy our mantissa into the dividend, initialize the |
1088 | quotient while we are at it */ |
1089 | |
1090 | for (int i=0; i<DIGITS; i++) |
1091 | { |
1092 | dividend[i] = mant[i]; |
1093 | quotient[i] = 0; |
1094 | remainder[i] = 0; |
1095 | } |
1096 | |
1097 | /* outer loop. Once per quotient digit */ |
1098 | nsqd = 0; |
1099 | for (qd = DIGITS+1; qd >= 0; qd--) |
1100 | { |
1101 | /* Determine outer limits of our quotient digit */ |
1102 | |
1103 | // r = most sig 2 digits of dividend |
1104 | r = dividend[DIGITS]*radix+dividend[DIGITS-1]; |
1105 | min = r / (divisor.mant[DIGITS-1]+1); |
1106 | max = (r+1) / divisor.mant[DIGITS-1]; |
1107 | |
1108 | trialgood = false; |
1109 | |
1110 | while (!trialgood) |
1111 | { |
1112 | // try the mean |
1113 | trial = (min+max)/2; |
1114 | |
1115 | //System.out.println("dividend = "+dividend[2]+" "+dividend[1]+" "+dividend[0]); |
1116 | //System.out.println("divisor = "+divisor.mant[1]+" "+divisor.mant[0]); |
1117 | //System.out.println("min = "+min+" max = "+max+" trial = "+trial); |
1118 | |
1119 | /* Multiply by divisor and store as remainder */ |
1120 | rh = 0; |
1121 | for (int i=0; i<(DIGITS+1); i++) |
1122 | { |
1123 | int dm = (i<DIGITS)?divisor.mant[i]:0; |
1124 | r = (dm * trial) + rh; |
1125 | |
1126 | rh = r / radix; |
1127 | rl = r % radix; |
1128 | remainder[i] = rl; |
1129 | } |
1130 | |
1131 | //System.out.println(" *remainder = "+remainder[1]+" "+remainder[0]); |
1132 | |
1133 | /* subtract the remainder from the dividend */ |
1134 | rh = 1; // carry in to aid the subtraction |
1135 | for (int i=0; i<(DIGITS+1); i++) |
1136 | { |
1137 | r = ((radix-1) - remainder[i]) + dividend[i] + rh; |
1138 | rh = r / radix; |
1139 | rl = r % radix; |
1140 | remainder[i] = rl; |
1141 | } |
1142 | //System.out.println(" +remainder = "+remainder[1]+" "+remainder[0]); |
1143 | |
1144 | /* Lets analyse what we have here */ |
1145 | if (rh == 0) // trial is too big -- negative remainder |
1146 | { |
1147 | max = trial-1; |
1148 | //System.out.println("neg remainder"); |
1149 | //System.out.println(" remainder = "+remainder[1]+" "+remainder[0]); |
1150 | //System.out.println(" rh = "+rh); |
1151 | continue; |
1152 | } |
1153 | |
1154 | /* find out how far off the remainder is telling us we are */ |
1155 | minadj = (remainder[DIGITS] * radix)+remainder[DIGITS-1]; |
1156 | minadj = minadj / (divisor.mant[DIGITS-1]+1); |
1157 | |
1158 | //System.out.println("minadj = "+minadj); |
1159 | |
1160 | if (minadj >= 2) |
1161 | { |
1162 | min = trial+minadj; // update the minium |
1163 | continue; |
1164 | } |
1165 | |
1166 | /* May have a good one here, check more thoughly. Basically |
1167 | its a good one if it is less than the divisor */ |
1168 | trialgood = false; // assume false |
1169 | for (int i=(DIGITS-1); i >=0; i--) |
1170 | { |
1171 | if (divisor.mant[i] > remainder[i]) |
1172 | trialgood = true; |
1173 | if (divisor.mant[i] < remainder[i]) |
1174 | break; |
1175 | } |
1176 | |
1177 | //System.out.println("remainder = "+remainder[1]+" "+remainder[0]); |
1178 | if (remainder[DIGITS] != 0) |
1179 | trialgood = false; |
1180 | |
1181 | if (trialgood == false) |
1182 | min = trial+1; |
1183 | } |
1184 | |
1185 | /* Great we have a digit! */ |
1186 | quotient[qd] = trial; |
1187 | if (trial != 0 || nsqd != 0) |
1188 | nsqd++; |
1189 | |
1190 | if (rMode == ROUND_DOWN && nsqd == DIGITS) // We have enough for this mode |
1191 | break; |
1192 | |
1193 | if (nsqd > DIGITS) // We have enough digits |
1194 | break; |
1195 | |
1196 | /* move the remainder into the dividend while left shifting */ |
1197 | dividend[0] = 0; |
1198 | for (int i=0; i<DIGITS; i++) |
1199 | dividend[i+1] = remainder[i]; |
1200 | } |
1201 | |
1202 | /* Find the most sig digit */ |
1203 | md = DIGITS; // default |
1204 | for (int i=DIGITS+1; i>=0; i--) |
1205 | { |
1206 | if (quotient[i] != 0) |
1207 | { |
1208 | md = i; |
1209 | break; |
1210 | } |
1211 | } |
1212 | |
1213 | //System.out.println("quotient = "+quotient[2]+" "+quotient[1]+" "+quotient[0]); |
1214 | |
1215 | /* Copy the digits into the result */ |
1216 | for (int i=0; i<DIGITS; i++) |
1217 | result.mant[DIGITS-i-1] = quotient[md-i]; |
1218 | |
1219 | /* Fixup the exponent. */ |
1220 | result.exp = (exp - divisor.exp + md - DIGITS + 1 - 1); |
1221 | result.sign = (byte)((sign == divisor.sign)?1:-1); |
1222 | |
1223 | if (result.mant[DIGITS-1] == 0) // if result is zero, set exp to zero |
1224 | result.exp = 0; |
1225 | |
1226 | if (md > (DIGITS-1)) |
1227 | excp = result.round(quotient[md-DIGITS]); |
1228 | else |
1229 | excp = result.round(0); |
1230 | |
1231 | if (excp != 0) |
1232 | result = dotrap(excp, "divide", divisor, result); |
1233 | |
1234 | return result; |
1235 | } |
1236 | |
1237 | /** Divide by a single digit less than radix. |
1238 | * Special case, so there are speed advantages. |
1239 | * 0 <= divisor < radix */ |
1240 | public dfp divide(int divisor) |
1241 | { |
1242 | dfp result; |
1243 | int r, rh, rl; |
1244 | int excp; |
1245 | |
1246 | /* handle special cases */ |
1247 | if (nans != FINITE) |
1248 | { |
1249 | if (nans == QNAN || nans == SNAN) |
1250 | return this; |
1251 | |
1252 | if (nans == INFINITE) |
1253 | { |
1254 | result = newInstance(this); |
1255 | return result; |
1256 | } |
1257 | } |
1258 | |
1259 | /* Test for divide by zero */ |
1260 | if (divisor == 0) |
1261 | { |
1262 | ieeeFlags |= FLAG_DIV_ZERO; |
1263 | result = newInstance(zero); |
1264 | result.sign = sign; |
1265 | result.nans = INFINITE; |
1266 | result = dotrap(FLAG_DIV_ZERO, "divide", zero, result); |
1267 | return result; |
1268 | } |
1269 | |
1270 | /* range check divisor */ |
1271 | if (divisor < 0 || divisor >= radix) |
1272 | { |
1273 | ieeeFlags |= FLAG_INVALID; |
1274 | result = newInstance(zero); |
1275 | result.nans = QNAN; |
1276 | result = dotrap(FLAG_INVALID, "divide", result, result); |
1277 | return result; |
1278 | } |
1279 | |
1280 | result = newInstance(this); |
1281 | |
1282 | rl = 0; |
1283 | for (int i=DIGITS-1; i>=0; i--) |
1284 | { |
1285 | r = rl*radix + result.mant[i]; |
1286 | rh = r / divisor; |
1287 | rl = r % divisor; |
1288 | result.mant[i] = rh; |
1289 | } |
1290 | |
1291 | if (result.mant[DIGITS-1] == 0) // normalize |
1292 | { |
1293 | result.shiftLeft(); |
1294 | r = rl*radix; // compute the next digit and put it in |
1295 | rh = r / divisor; |
1296 | rl = r % divisor; |
1297 | result.mant[0] = rh; |
1298 | } |
1299 | |
1300 | excp = result.round(rl*radix/divisor); // do the rounding |
1301 | |
1302 | if (excp != 0) |
1303 | result = dotrap(excp, "divide", result, result); |
1304 | |
1305 | return result; |
1306 | } |
1307 | |
1308 | public dfp sqrt() /* returns the square root of this */ |
1309 | { |
1310 | dfp x, dx, px; |
1311 | |
1312 | /* check for unusual cases */ |
1313 | if (nans == FINITE && mant[DIGITS-1] == 0) // if zero |
1314 | return newInstance(this); |
1315 | |
1316 | if (nans != FINITE) |
1317 | { |
1318 | if (nans == INFINITE && sign == 1) // if positive infinity |
1319 | return newInstance(this); |
1320 | |
1321 | if (nans == QNAN) |
1322 | return newInstance(this); |
1323 | |
1324 | if (nans == SNAN) |
1325 | { |
1326 | dfp result; |
1327 | |
1328 | ieeeFlags |= FLAG_INVALID; |
1329 | result = newInstance(this); |
1330 | result = dotrap(FLAG_INVALID, "sqrt", null, result); |
1331 | return result; |
1332 | } |
1333 | } |
1334 | |
1335 | if (sign == -1) // if negative |
1336 | { |
1337 | dfp result; |
1338 | |
1339 | ieeeFlags |= FLAG_INVALID; |
1340 | result = newInstance(this); |
1341 | result.nans = QNAN; |
1342 | result = dotrap(FLAG_INVALID, "sqrt", null, result); |
1343 | return result; |
1344 | } |
1345 | |
1346 | x = newInstance(this); |
1347 | |
1348 | /* Lets make a reasonable guess as to the size of the square root */ |
1349 | if (x.exp < -1 || x.exp > 1) |
1350 | x.exp = (this.exp/2); |
1351 | |
1352 | /* Coarsely estimate the mantissa */ |
1353 | switch (x.mant[DIGITS-1] / 2000) |
1354 | { |
1355 | case 0: x.mant[DIGITS-1] = x.mant[DIGITS-1]/2+1; break; |
1356 | case 2: x.mant[DIGITS-1] = 1500; break; |
1357 | case 3: x.mant[DIGITS-1] = 2200; break; |
1358 | case 4: x.mant[DIGITS-1] = 3000; break; |
1359 | } |
1360 | |
1361 | dx = newInstance(x); |
1362 | |
1363 | /* Now that we have the first pass estimiate, compute the rest |
1364 | by the formula dx = (y - x*x) / (2x); */ |
1365 | |
1366 | px = zero; |
1367 | while (x.unequal(px)) |
1368 | { |
1369 | dx = newInstance(x); |
1370 | dx.sign = -1; |
1371 | dx = dx.add(this.divide(x)); |
1372 | dx = dx.divide(2); |
1373 | px = x; |
1374 | x = x.add(dx); |
1375 | |
1376 | // if dx is zero, break. Note testing the most sig digit |
1377 | // is a sufficient test since dx is normalized |
1378 | if (dx.mant[DIGITS-1] == 0) |
1379 | break; |
1380 | } |
1381 | |
1382 | return x; |
1383 | } |
1384 | |
1385 | public String toString() |
1386 | { |
1387 | if (nans != FINITE) // if non-finite exceptional cases |
1388 | { |
1389 | switch (sign*nans) |
1390 | { |
1391 | case INFINITE: return "Infinity"; |
1392 | case -INFINITE: return "-Infinity"; |
1393 | case QNAN: return "NaN"; |
1394 | case -QNAN: return "NaN"; |
1395 | case SNAN: return "NaN"; |
1396 | case -SNAN: return "NaN"; |
1397 | } |
1398 | } |
1399 | |
1400 | if (exp > DIGITS || exp < -1) |
1401 | return dfp2sci(this); |
1402 | |
1403 | return dfp2string(this); |
1404 | } |
1405 | |
1406 | /* Convert a dfp to a string using scientific notation */ |
1407 | protected String dfp2sci(dfp a) |
1408 | { |
1409 | char rawdigits[] = new char[DIGITS*4]; |
1410 | char outputbuffer[] = new char[DIGITS*4 + 20]; |
1411 | int p; |
1412 | int q; |
1413 | int e; |
1414 | int ae; |
1415 | int shf; |
1416 | |
1417 | |
1418 | /* Get all the digits */ |
1419 | p = 0; |
1420 | for (int i=DIGITS-1; i>=0; i--) |
1421 | { |
1422 | rawdigits[p++] = (char) ((a.mant[i]/1000) + '0'); |
1423 | rawdigits[p++] = (char) (((a.mant[i]/100)%10) + '0'); |
1424 | rawdigits[p++] = (char) (((a.mant[i]/10)%10) + '0'); |
1425 | rawdigits[p++] = (char) (((a.mant[i])%10) + '0'); |
1426 | } |
1427 | |
1428 | /* find the first non-zero one */ |
1429 | for (p=0; p<rawdigits.length; p++) |
1430 | if (rawdigits[p] != '0') |
1431 | break; |
1432 | shf = p; |
1433 | |
1434 | /* Now do the conversion */ |
1435 | q = 0; |
1436 | if (a.sign == -1) |
1437 | outputbuffer[q++] = '-'; |
1438 | |
1439 | if (p != rawdigits.length) // there are non zero digits... |
1440 | { |
1441 | outputbuffer[q++] = rawdigits[p++]; |
1442 | outputbuffer[q++] = '.'; |
1443 | |
1444 | while (p<rawdigits.length) |
1445 | outputbuffer[q++] = rawdigits[p++]; |
1446 | } |
1447 | else |
1448 | { |
1449 | outputbuffer[q++] = '0'; |
1450 | outputbuffer[q++] = '.'; |
1451 | outputbuffer[q++] = '0'; |
1452 | outputbuffer[q++] = 'e'; |
1453 | outputbuffer[q++] = '0'; |
1454 | return new String(outputbuffer, 0, 5); |
1455 | } |
1456 | |
1457 | outputbuffer[q++] = 'e'; |
1458 | |
1459 | /* find the msd of the exponent */ |
1460 | |
1461 | e = a.exp * 4 - shf - 1; |
1462 | ae = e; |
1463 | if (e < 0) |
1464 | ae = -e; |
1465 | |
1466 | /* Find the largest p such that p < e */ |
1467 | for (p=1000000000; p>ae; p /= 10); |
1468 | |
1469 | if (e < 0) |
1470 | outputbuffer[q++] = '-'; |
1471 | |
1472 | while (p > 0) |
1473 | { |
1474 | outputbuffer[q++] = (char)(ae/p + '0'); |
1475 | ae = ae % p; |
1476 | p = p / 10; |
1477 | } |
1478 | |
1479 | return new String(outputbuffer, 0, q); |
1480 | } |
1481 | |
1482 | /* converts a dfp to a string handling the normal case */ |
1483 | protected String dfp2string(dfp a) |
1484 | { |
1485 | char buffer[] = new char[DIGITS*4 + 20]; |
1486 | int p = 1; |
1487 | int q; |
1488 | int e = a.exp; |
1489 | boolean pointInserted = false; |
1490 | |
1491 | buffer[0] = ' '; |
1492 | |
1493 | if (e <= 0) |
1494 | { |
1495 | buffer[p++] = '0'; |
1496 | buffer[p++] = '.'; |
1497 | pointInserted = true; |
1498 | } |
1499 | |
1500 | while (e < 0) |
1501 | { |
1502 | buffer[p++] = '0'; |
1503 | buffer[p++] = '0'; |
1504 | buffer[p++] = '0'; |
1505 | buffer[p++] = '0'; |
1506 | e++; |
1507 | } |
1508 | |
1509 | for (int i=DIGITS-1; i>=0; i--) |
1510 | { |
1511 | buffer[p++] = (char) ((a.mant[i]/1000) + '0'); |
1512 | buffer[p++] = (char) (((a.mant[i]/100)%10) + '0'); |
1513 | buffer[p++] = (char) (((a.mant[i]/10)%10) + '0'); |
1514 | buffer[p++] = (char) (((a.mant[i])%10) + '0'); |
1515 | if (--e == 0) |
1516 | { |
1517 | buffer[p++] = '.'; |
1518 | pointInserted = true; |
1519 | } |
1520 | } |
1521 | |
1522 | while (e > 0) |
1523 | { |
1524 | buffer[p++] = '0'; |
1525 | buffer[p++] = '0'; |
1526 | buffer[p++] = '0'; |
1527 | buffer[p++] = '0'; |
1528 | e--; |
1529 | } |
1530 | |
1531 | if (!pointInserted) /* Ensure we have a radix point! */ |
1532 | buffer[p++] = '.'; |
1533 | |
1534 | /* Supress leading zeros */ |
1535 | q = 1; |
1536 | while (buffer[q] == '0') |
1537 | q++; |
1538 | if (buffer[q] == '.') |
1539 | q--; |
1540 | |
1541 | /* Suppress trailing zeros */ |
1542 | while (buffer[p-1] == '0') |
1543 | p--; |
1544 | |
1545 | /* Insert sign */ |
1546 | if (a.sign < 0) |
1547 | buffer[--q] = '-'; |
1548 | |
1549 | return new String(buffer, q, p-q); |
1550 | } |
1551 | |
1552 | protected dfp string2dfp(String fpin) |
1553 | { |
1554 | String fpdecimal; |
1555 | int trailing_zeros; |
1556 | int significant_digits; |
1557 | int decimal; |
1558 | int p, q; |
1559 | char Striped[]; |
1560 | dfp result; |
1561 | boolean decimalFound = false; |
1562 | int decimalPos = 0; // position of the decimal. |
1563 | final int rsize = 4; // size of radix in decimal digits |
1564 | final int offset = 4; // Starting offset into Striped |
1565 | int sciexp = 0; |
1566 | int i; |
1567 | |
1568 | Striped = new char[DIGITS*rsize+offset*2]; |
1569 | |
1570 | /* Check some special cases */ |
1571 | if (fpin.equals("Infinite")) |
1572 | return create((byte) 1, (byte) INFINITE); |
1573 | |
1574 | if (fpin.equals("-Infinite")) |
1575 | return create((byte) -1, (byte) INFINITE); |
1576 | |
1577 | if (fpin.equals("NaN")) |
1578 | return create((byte) 1, (byte) QNAN); |
1579 | |
1580 | result = newInstance(zero); |
1581 | |
1582 | /* Check for scientific notation */ |
1583 | p = fpin.indexOf("e"); |
1584 | if (p == -1) // try upper case? |
1585 | p = fpin.indexOf("E"); |
1586 | |
1587 | if (p != -1) // scientific notation |
1588 | { |
1589 | fpdecimal = fpin.substring(0, p); |
1590 | String fpexp = fpin.substring(p+1); |
1591 | boolean negative = false; |
1592 | |
1593 | sciexp = 0; |
1594 | for (i=0; i<fpexp.length(); i++) |
1595 | { |
1596 | if (fpexp.charAt(i) == '-') |
1597 | { |
1598 | negative = true; |
1599 | continue; |
1600 | } |
1601 | if (fpexp.charAt(i) >= '0' && fpexp.charAt(i) <= '9') |
1602 | sciexp = sciexp * 10 + fpexp.charAt(i) - '0'; |
1603 | } |
1604 | |
1605 | if (negative) |
1606 | sciexp = -sciexp; |
1607 | } |
1608 | else // normal case |
1609 | { |
1610 | fpdecimal = fpin; |
1611 | } |
1612 | |
1613 | /* If there is a minus sign in the number then it is negative */ |
1614 | |
1615 | if (fpdecimal.indexOf("-") != -1) |
1616 | result.sign = -1; |
1617 | |
1618 | /* First off, find all of the leading zeros, trailing zeros, and |
1619 | siginificant digits */ |
1620 | |
1621 | p = 0; |
1622 | |
1623 | /* Move p to first significant digit */ |
1624 | |
1625 | for(;;) |
1626 | { |
1627 | if (fpdecimal.charAt(p) >= '1' && fpdecimal.charAt(p) <= '9') |
1628 | break; |
1629 | |
1630 | if (decimalFound && fpdecimal.charAt(p) == '0') |
1631 | decimalPos--; |
1632 | |
1633 | if (fpdecimal.charAt(p) == '.') |
1634 | decimalFound = true; |
1635 | |
1636 | p++; |
1637 | |
1638 | if (p == fpdecimal.length()) |
1639 | break; |
1640 | } |
1641 | |
1642 | /* Copy the string onto Stripped */ |
1643 | |
1644 | q = offset; |
1645 | Striped[0] = '0'; |
1646 | Striped[1] = '0'; |
1647 | Striped[2] = '0'; |
1648 | Striped[3] = '0'; |
1649 | significant_digits=0; |
1650 | for(;;) |
1651 | { |
1652 | if (p == (fpdecimal.length())) |
1653 | break; |
1654 | |
1655 | // Dont want to run pass the end of the array |
1656 | if (q == DIGITS*rsize+offset+1) |
1657 | break; |
1658 | |
1659 | if (fpdecimal.charAt(p) == '.') |
1660 | { |
1661 | decimalFound = true; |
1662 | decimalPos = significant_digits; |
1663 | p++; |
1664 | continue; |
1665 | } |
1666 | |
1667 | if (fpdecimal.charAt(p) < '0' || fpdecimal.charAt(p) > '9') |
1668 | { |
1669 | p++; |
1670 | continue; |
1671 | } |
1672 | |
1673 | Striped[q] = fpdecimal.charAt(p); |
1674 | q++; |
1675 | p++; |
1676 | significant_digits++; |
1677 | } |
1678 | |
1679 | |
1680 | // If the decimal point has been found then get rid of trailing zeros. |
1681 | if (decimalFound && q != offset) |
1682 | { |
1683 | for (;;) |
1684 | { |
1685 | q--; |
1686 | if (q == offset) |
1687 | break; |
1688 | if (Striped[q] == '0') |
1689 | { |
1690 | significant_digits--; |
1691 | } |
1692 | else |
1693 | { |
1694 | break; |
1695 | } |
1696 | } |
1697 | } |
1698 | |
1699 | // special case of numbers like "0.00000" |
1700 | if (decimalFound && significant_digits == 0) |
1701 | decimalPos = 0; |
1702 | |
1703 | // Implicit decimal point at end of number if not present |
1704 | if (!decimalFound) |
1705 | decimalPos = q-offset; |
1706 | |
1707 | /* Find the number of significant trailing zeros */ |
1708 | |
1709 | q = offset; // set q to point to first sig digit |
1710 | p = significant_digits-1+offset; |
1711 | |
1712 | trailing_zeros = 0; |
1713 | while (p>q) |
1714 | { |
1715 | if (Striped[p] != '0') |
1716 | break; |
1717 | trailing_zeros++; |
1718 | p--; |
1719 | } |
1720 | |
1721 | /* Make sure the decimal is on a mod 10000 boundary */ |
1722 | i = (((rsize*100) - decimalPos - sciexp%rsize) % rsize); |
1723 | q -= i; |
1724 | decimalPos += i; |
1725 | |
1726 | /* Make the mantissa length right by adding zeros at the end if |
1727 | necessary */ |
1728 | |
1729 | while ( (p-q) < (DIGITS*rsize) ) |
1730 | { |
1731 | for (i=0; i<rsize; i++) |
1732 | Striped[++p] = '0'; |
1733 | } |
1734 | |
1735 | /* Ok, now we know how many trailing zeros there are, and |
1736 | where the least significant digit is. */ |
1737 | |
1738 | for (i=(DIGITS-1); i>=0; i--) |
1739 | { |
1740 | result.mant[i] = (Striped[q] - '0')*1000 + |
1741 | (Striped[q+1] - '0')*100 + |
1742 | (Striped[q+2] - '0')*10 + |
1743 | (Striped[q+3] - '0'); |
1744 | q += 4; |
1745 | } |
1746 | |
1747 | |
1748 | result.exp = ((decimalPos+sciexp) / rsize); |
1749 | |
1750 | if (q < Striped.length) // Is there possible another digit? |
1751 | result.round((Striped[q] - '0')*1000); |
1752 | |
1753 | return result; |
1754 | } |
1755 | |
1756 | /** Set the rounding mode to be one of the following values: |
1757 | * ROUND_UP, ROUND_DOWN, ROUND_HALF_UP, ROUND_HALF_DOWN, |
1758 | * ROUND_HALF_EVEN, ROUND_HALF_ODD, ROUND_CEIL, ROUND_FLOOR. |
1759 | * |
1760 | * Default is ROUND_HALF_EVEN |
1761 | * |
1762 | * Note that the rounding mode is common to all instances |
1763 | * in the system and will effect all future calculations. |
1764 | */ |
1765 | public static void setRoundingMode(int mode) |
1766 | { |
1767 | rMode = mode; |
1768 | } |
1769 | |
1770 | /** Returns the current rounding mode */ |
1771 | public static int getRoundingMode() |
1772 | { |
1773 | return rMode; |
1774 | } |
1775 | |
1776 | /** Raises a trap. This does not set the corresponding flag however. |
1777 | * @param type the trap type |
1778 | * @param what - name of routine trap occured in |
1779 | * @param oper - input operator to function |
1780 | * @param result - the result computed prior to the trap |
1781 | * @return The suggested return value from the trap handler |
1782 | */ |
1783 | public dfp dotrap(int type, String what, dfp oper, dfp result) |
1784 | { |
1785 | dfp def = result; |
1786 | |
1787 | switch (type) |
1788 | { |
1789 | case FLAG_INVALID: |
1790 | def = newInstance(zero); |
1791 | def.sign = result.sign; |
1792 | def.nans = QNAN; |
1793 | break; |
1794 | |
1795 | case FLAG_DIV_ZERO: |
1796 | if (nans == FINITE && mant[DIGITS-1] != 0) // normal case, we are finite, non-zero |
1797 | { |
1798 | def = newInstance(zero); |
1799 | def.sign = (byte)(sign*oper.sign); |
1800 | def.nans = INFINITE; |
1801 | } |
1802 | |
1803 | if (nans == FINITE && mant[DIGITS-1] == 0) // 0/0 |
1804 | { |
1805 | def = newInstance(zero); |
1806 | def.nans = QNAN; |
1807 | } |
1808 | |
1809 | if (nans == INFINITE || nans == QNAN) |
1810 | { |
1811 | def = newInstance(zero); |
1812 | def.nans = QNAN; |
1813 | } |
1814 | |
1815 | if (nans == INFINITE || nans == SNAN) |
1816 | { |
1817 | def = newInstance(zero); |
1818 | def.nans = QNAN; |
1819 | } |
1820 | break; |
1821 | |
1822 | case FLAG_UNDERFLOW: |
1823 | if ( (result.exp+DIGITS) < minExp) |
1824 | { |
1825 | def = newInstance(zero); |
1826 | def.sign = result.sign; |
1827 | } |
1828 | else |
1829 | { |
1830 | def = newInstance(result); // gradual underflow |
1831 | } |
1832 | result.exp = result.exp + errScale; |
1833 | break; |
1834 | |
1835 | case FLAG_OVERFLOW: |
1836 | result.exp = result.exp - errScale; |
1837 | def = newInstance(zero); |
1838 | def.sign = result.sign; |
1839 | def.nans = INFINITE; |
1840 | break; |
1841 | |
1842 | default: def = result; break; |
1843 | } |
1844 | |
1845 | return trap(type, what, oper, def, result); |
1846 | } |
1847 | |
1848 | /** Trap handler. Subclasses may override this to provide trap |
1849 | * functionality per IEEE 854-1987. |
1850 | * |
1851 | * @param type The exception type - e.g. FLAG_OVERFLOW |
1852 | * @param what The name of the routine we were in e.g. divide() |
1853 | * @param oper An operand to this function if any |
1854 | * @param def The default return value if trap not enabled |
1855 | * @param result The result that is spcefied to be delivered per |
1856 | * IEEE 854, if any |
1857 | */ |
1858 | protected dfp trap(int type, String what, dfp oper, |
1859 | dfp def, dfp result) |
1860 | { |
1861 | return def; |
1862 | } |
1863 | |
1864 | /** Returns the IEEE 854 status flags */ |
1865 | public static int getIEEEFlags() |
1866 | { |
1867 | return ieeeFlags; |
1868 | } |
1869 | |
1870 | /** Clears the IEEE 854 status flags */ |
1871 | public static void clearIEEEFlags() |
1872 | { |
1873 | ieeeFlags = 0; |
1874 | } |
1875 | |
1876 | /** Sets the IEEE 854 status flags */ |
1877 | public static void setIEEEFlags(int flags) |
1878 | { |
1879 | ieeeFlags = flags; |
1880 | } |
1881 | |
1882 | /** Returns the type - one of FINITE, INFINITE, SNAN, QNAN */ |
1883 | public int classify() |
1884 | { |
1885 | return nans; |
1886 | } |
1887 | |
1888 | /** Creates a dfp with a non-finite value */ |
1889 | public static dfp create(byte sign, byte nans) |
1890 | { |
1891 | dfp result = new dfp(); |
1892 | result.sign = sign; |
1893 | result.nans = nans; |
1894 | |
1895 | return result; |
1896 | } |
1897 | |
1898 | /** Creates a dfp that is the same as x except that it has |
1899 | * the sign of y. abs(x) = dfp.copysign(x, dfp.one) */ |
1900 | public static dfp copysign(dfp x, dfp y) |
1901 | { |
1902 | dfp result = x.newInstance(x); |
1903 | result.sign = y.sign; |
1904 | return result; |
1905 | } |
1906 | |
1907 | /** Returns the next number greater than this one in the direction |
1908 | * of x. If this==x then simply returns this. */ |
1909 | |
1910 | public dfp nextAfter(dfp x) |
1911 | { |
1912 | boolean up = false; |
1913 | dfp result, inc; |
1914 | |
1915 | // if this is greater than x |
1916 | if (this.lessThan(x)) |
1917 | up = true; |
1918 | |
1919 | if (compare(this, x) == 0) |
1920 | return newInstance(x); |
1921 | |
1922 | if (lessThan(zero)) |
1923 | up = !up; |
1924 | |
1925 | if (up) |
1926 | { |
1927 | inc = newInstance(one); |
1928 | inc.exp = this.exp-DIGITS+1; |
1929 | inc.sign = this.sign; |
1930 | |
1931 | if (this.equal(zero)) |
1932 | inc.exp = minExp-DIGITS; |
1933 | |
1934 | result = add(inc); |
1935 | } |
1936 | else |
1937 | { |
1938 | inc = newInstance(one); |
1939 | inc.exp = this.exp; |
1940 | inc.sign = this.sign; |
1941 | |
1942 | if (this.equal(inc)) |
1943 | inc.exp = this.exp-DIGITS; |
1944 | else |
1945 | inc.exp = this.exp-DIGITS+1; |
1946 | |
1947 | if (this.equal(zero)) |
1948 | inc.exp = minExp-DIGITS; |
1949 | |
1950 | result = this.subtract(inc); |
1951 | } |
1952 | |
1953 | if (result.classify() == INFINITE && this.classify() != INFINITE) |
1954 | { |
1955 | ieeeFlags |= FLAG_INEXACT; |
1956 | result = dotrap(FLAG_INEXACT, "nextAfter", x, result); |
1957 | } |
1958 | |
1959 | if (result.equal(zero) && this.equal(zero) == false) |
1960 | { |
1961 | ieeeFlags |= FLAG_INEXACT; |
1962 | result = dotrap(FLAG_INEXACT, "nextAfter", x, result); |
1963 | } |
1964 | |
1965 | return result; |
1966 | } |
1967 | } |