001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math3.dfp;
019
020 import java.util.Arrays;
021
022 import org.apache.commons.math3.FieldElement;
023
024 /**
025 * Decimal floating point library for Java
026 *
027 * <p>Another floating point class. This one is built using radix 10000
028 * which is 10<sup>4</sup>, so its almost decimal.</p>
029 *
030 * <p>The design goals here are:
031 * <ol>
032 * <li>Decimal math, or close to it</li>
033 * <li>Settable precision (but no mix between numbers using different settings)</li>
034 * <li>Portability. Code should be kept as portable as possible.</li>
035 * <li>Performance</li>
036 * <li>Accuracy - Results should always be +/- 1 ULP for basic
037 * algebraic operation</li>
038 * <li>Comply with IEEE 854-1987 as much as possible.
039 * (See IEEE 854-1987 notes below)</li>
040 * </ol></p>
041 *
042 * <p>Trade offs:
043 * <ol>
044 * <li>Memory foot print. I'm using more memory than necessary to
045 * represent numbers to get better performance.</li>
046 * <li>Digits are bigger, so rounding is a greater loss. So, if you
047 * really need 12 decimal digits, better use 4 base 10000 digits
048 * there can be one partially filled.</li>
049 * </ol></p>
050 *
051 * <p>Numbers are represented in the following form:
052 * <pre>
053 * n = sign × mant × (radix)<sup>exp</sup>;</p>
054 * </pre>
055 * where sign is ±1, mantissa represents a fractional number between
056 * zero and one. mant[0] is the least significant digit.
057 * exp is in the range of -32767 to 32768</p>
058 *
059 * <p>IEEE 854-1987 Notes and differences</p>
060 *
061 * <p>IEEE 854 requires the radix to be either 2 or 10. The radix here is
062 * 10000, so that requirement is not met, but it is possible that a
063 * subclassed can be made to make it behave as a radix 10
064 * number. It is my opinion that if it looks and behaves as a radix
065 * 10 number then it is one and that requirement would be met.</p>
066 *
067 * <p>The radix of 10000 was chosen because it should be faster to operate
068 * on 4 decimal digits at once instead of one at a time. Radix 10 behavior
069 * can be realized by adding an additional rounding step to ensure that
070 * the number of decimal digits represented is constant.</p>
071 *
072 * <p>The IEEE standard specifically leaves out internal data encoding,
073 * so it is reasonable to conclude that such a subclass of this radix
074 * 10000 system is merely an encoding of a radix 10 system.</p>
075 *
076 * <p>IEEE 854 also specifies the existence of "sub-normal" numbers. This
077 * class does not contain any such entities. The most significant radix
078 * 10000 digit is always non-zero. Instead, we support "gradual underflow"
079 * by raising the underflow flag for numbers less with exponent less than
080 * expMin, but don't flush to zero until the exponent reaches MIN_EXP-digits.
081 * Thus the smallest number we can represent would be:
082 * 1E(-(MIN_EXP-digits-1)*4), eg, for digits=5, MIN_EXP=-32767, that would
083 * be 1e-131092.</p>
084 *
085 * <p>IEEE 854 defines that the implied radix point lies just to the right
086 * of the most significant digit and to the left of the remaining digits.
087 * This implementation puts the implied radix point to the left of all
088 * digits including the most significant one. The most significant digit
089 * here is the one just to the right of the radix point. This is a fine
090 * detail and is really only a matter of definition. Any side effects of
091 * this can be rendered invisible by a subclass.</p>
092 * @see DfpField
093 * @version $Id: Dfp.java 1416643 2012-12-03 19:37:14Z tn $
094 * @since 2.2
095 */
096 public class Dfp implements FieldElement<Dfp> {
097
098 /** The radix, or base of this system. Set to 10000 */
099 public static final int RADIX = 10000;
100
101 /** The minimum exponent before underflow is signaled. Flush to zero
102 * occurs at minExp-DIGITS */
103 public static final int MIN_EXP = -32767;
104
105 /** The maximum exponent before overflow is signaled and results flushed
106 * to infinity */
107 public static final int MAX_EXP = 32768;
108
109 /** The amount under/overflows are scaled by before going to trap handler */
110 public static final int ERR_SCALE = 32760;
111
112 /** Indicator value for normal finite numbers. */
113 public static final byte FINITE = 0;
114
115 /** Indicator value for Infinity. */
116 public static final byte INFINITE = 1;
117
118 /** Indicator value for signaling NaN. */
119 public static final byte SNAN = 2;
120
121 /** Indicator value for quiet NaN. */
122 public static final byte QNAN = 3;
123
124 /** String for NaN representation. */
125 private static final String NAN_STRING = "NaN";
126
127 /** String for positive infinity representation. */
128 private static final String POS_INFINITY_STRING = "Infinity";
129
130 /** String for negative infinity representation. */
131 private static final String NEG_INFINITY_STRING = "-Infinity";
132
133 /** Name for traps triggered by addition. */
134 private static final String ADD_TRAP = "add";
135
136 /** Name for traps triggered by multiplication. */
137 private static final String MULTIPLY_TRAP = "multiply";
138
139 /** Name for traps triggered by division. */
140 private static final String DIVIDE_TRAP = "divide";
141
142 /** Name for traps triggered by square root. */
143 private static final String SQRT_TRAP = "sqrt";
144
145 /** Name for traps triggered by alignment. */
146 private static final String ALIGN_TRAP = "align";
147
148 /** Name for traps triggered by truncation. */
149 private static final String TRUNC_TRAP = "trunc";
150
151 /** Name for traps triggered by nextAfter. */
152 private static final String NEXT_AFTER_TRAP = "nextAfter";
153
154 /** Name for traps triggered by lessThan. */
155 private static final String LESS_THAN_TRAP = "lessThan";
156
157 /** Name for traps triggered by greaterThan. */
158 private static final String GREATER_THAN_TRAP = "greaterThan";
159
160 /** Name for traps triggered by newInstance. */
161 private static final String NEW_INSTANCE_TRAP = "newInstance";
162
163 /** Mantissa. */
164 protected int[] mant;
165
166 /** Sign bit: 1 for positive, -1 for negative. */
167 protected byte sign;
168
169 /** Exponent. */
170 protected int exp;
171
172 /** Indicator for non-finite / non-number values. */
173 protected byte nans;
174
175 /** Factory building similar Dfp's. */
176 private final DfpField field;
177
178 /** Makes an instance with a value of zero.
179 * @param field field to which this instance belongs
180 */
181 protected Dfp(final DfpField field) {
182 mant = new int[field.getRadixDigits()];
183 sign = 1;
184 exp = 0;
185 nans = FINITE;
186 this.field = field;
187 }
188
189 /** Create an instance from a byte value.
190 * @param field field to which this instance belongs
191 * @param x value to convert to an instance
192 */
193 protected Dfp(final DfpField field, byte x) {
194 this(field, (long) x);
195 }
196
197 /** Create an instance from an int value.
198 * @param field field to which this instance belongs
199 * @param x value to convert to an instance
200 */
201 protected Dfp(final DfpField field, int x) {
202 this(field, (long) x);
203 }
204
205 /** Create an instance from a long value.
206 * @param field field to which this instance belongs
207 * @param x value to convert to an instance
208 */
209 protected Dfp(final DfpField field, long x) {
210
211 // initialize as if 0
212 mant = new int[field.getRadixDigits()];
213 nans = FINITE;
214 this.field = field;
215
216 boolean isLongMin = false;
217 if (x == Long.MIN_VALUE) {
218 // special case for Long.MIN_VALUE (-9223372036854775808)
219 // we must shift it before taking its absolute value
220 isLongMin = true;
221 ++x;
222 }
223
224 // set the sign
225 if (x < 0) {
226 sign = -1;
227 x = -x;
228 } else {
229 sign = 1;
230 }
231
232 exp = 0;
233 while (x != 0) {
234 System.arraycopy(mant, mant.length - exp, mant, mant.length - 1 - exp, exp);
235 mant[mant.length - 1] = (int) (x % RADIX);
236 x /= RADIX;
237 exp++;
238 }
239
240 if (isLongMin) {
241 // remove the shift added for Long.MIN_VALUE
242 // we know in this case that fixing the last digit is sufficient
243 for (int i = 0; i < mant.length - 1; i++) {
244 if (mant[i] != 0) {
245 mant[i]++;
246 break;
247 }
248 }
249 }
250 }
251
252 /** Create an instance from a double value.
253 * @param field field to which this instance belongs
254 * @param x value to convert to an instance
255 */
256 protected Dfp(final DfpField field, double x) {
257
258 // initialize as if 0
259 mant = new int[field.getRadixDigits()];
260 sign = 1;
261 exp = 0;
262 nans = FINITE;
263 this.field = field;
264
265 long bits = Double.doubleToLongBits(x);
266 long mantissa = bits & 0x000fffffffffffffL;
267 int exponent = (int) ((bits & 0x7ff0000000000000L) >> 52) - 1023;
268
269 if (exponent == -1023) {
270 // Zero or sub-normal
271 if (x == 0) {
272 // make sure 0 has the right sign
273 if ((bits & 0x8000000000000000L) != 0) {
274 sign = -1;
275 }
276 return;
277 }
278
279 exponent++;
280
281 // Normalize the subnormal number
282 while ( (mantissa & 0x0010000000000000L) == 0) {
283 exponent--;
284 mantissa <<= 1;
285 }
286 mantissa &= 0x000fffffffffffffL;
287 }
288
289 if (exponent == 1024) {
290 // infinity or NAN
291 if (x != x) {
292 sign = (byte) 1;
293 nans = QNAN;
294 } else if (x < 0) {
295 sign = (byte) -1;
296 nans = INFINITE;
297 } else {
298 sign = (byte) 1;
299 nans = INFINITE;
300 }
301 return;
302 }
303
304 Dfp xdfp = new Dfp(field, mantissa);
305 xdfp = xdfp.divide(new Dfp(field, 4503599627370496l)).add(field.getOne()); // Divide by 2^52, then add one
306 xdfp = xdfp.multiply(DfpMath.pow(field.getTwo(), exponent));
307
308 if ((bits & 0x8000000000000000L) != 0) {
309 xdfp = xdfp.negate();
310 }
311
312 System.arraycopy(xdfp.mant, 0, mant, 0, mant.length);
313 sign = xdfp.sign;
314 exp = xdfp.exp;
315 nans = xdfp.nans;
316
317 }
318
319 /** Copy constructor.
320 * @param d instance to copy
321 */
322 public Dfp(final Dfp d) {
323 mant = d.mant.clone();
324 sign = d.sign;
325 exp = d.exp;
326 nans = d.nans;
327 field = d.field;
328 }
329
330 /** Create an instance from a String representation.
331 * @param field field to which this instance belongs
332 * @param s string representation of the instance
333 */
334 protected Dfp(final DfpField field, final String s) {
335
336 // initialize as if 0
337 mant = new int[field.getRadixDigits()];
338 sign = 1;
339 exp = 0;
340 nans = FINITE;
341 this.field = field;
342
343 boolean decimalFound = false;
344 final int rsize = 4; // size of radix in decimal digits
345 final int offset = 4; // Starting offset into Striped
346 final char[] striped = new char[getRadixDigits() * rsize + offset * 2];
347
348 // Check some special cases
349 if (s.equals(POS_INFINITY_STRING)) {
350 sign = (byte) 1;
351 nans = INFINITE;
352 return;
353 }
354
355 if (s.equals(NEG_INFINITY_STRING)) {
356 sign = (byte) -1;
357 nans = INFINITE;
358 return;
359 }
360
361 if (s.equals(NAN_STRING)) {
362 sign = (byte) 1;
363 nans = QNAN;
364 return;
365 }
366
367 // Check for scientific notation
368 int p = s.indexOf("e");
369 if (p == -1) { // try upper case?
370 p = s.indexOf("E");
371 }
372
373 final String fpdecimal;
374 int sciexp = 0;
375 if (p != -1) {
376 // scientific notation
377 fpdecimal = s.substring(0, p);
378 String fpexp = s.substring(p+1);
379 boolean negative = false;
380
381 for (int i=0; i<fpexp.length(); i++)
382 {
383 if (fpexp.charAt(i) == '-')
384 {
385 negative = true;
386 continue;
387 }
388 if (fpexp.charAt(i) >= '0' && fpexp.charAt(i) <= '9') {
389 sciexp = sciexp * 10 + fpexp.charAt(i) - '0';
390 }
391 }
392
393 if (negative) {
394 sciexp = -sciexp;
395 }
396 } else {
397 // normal case
398 fpdecimal = s;
399 }
400
401 // If there is a minus sign in the number then it is negative
402 if (fpdecimal.indexOf("-") != -1) {
403 sign = -1;
404 }
405
406 // First off, find all of the leading zeros, trailing zeros, and significant digits
407 p = 0;
408
409 // Move p to first significant digit
410 int decimalPos = 0;
411 for (;;) {
412 if (fpdecimal.charAt(p) >= '1' && fpdecimal.charAt(p) <= '9') {
413 break;
414 }
415
416 if (decimalFound && fpdecimal.charAt(p) == '0') {
417 decimalPos--;
418 }
419
420 if (fpdecimal.charAt(p) == '.') {
421 decimalFound = true;
422 }
423
424 p++;
425
426 if (p == fpdecimal.length()) {
427 break;
428 }
429 }
430
431 // Copy the string onto Stripped
432 int q = offset;
433 striped[0] = '0';
434 striped[1] = '0';
435 striped[2] = '0';
436 striped[3] = '0';
437 int significantDigits=0;
438 for(;;) {
439 if (p == (fpdecimal.length())) {
440 break;
441 }
442
443 // Don't want to run pass the end of the array
444 if (q == mant.length*rsize+offset+1) {
445 break;
446 }
447
448 if (fpdecimal.charAt(p) == '.') {
449 decimalFound = true;
450 decimalPos = significantDigits;
451 p++;
452 continue;
453 }
454
455 if (fpdecimal.charAt(p) < '0' || fpdecimal.charAt(p) > '9') {
456 p++;
457 continue;
458 }
459
460 striped[q] = fpdecimal.charAt(p);
461 q++;
462 p++;
463 significantDigits++;
464 }
465
466
467 // If the decimal point has been found then get rid of trailing zeros.
468 if (decimalFound && q != offset) {
469 for (;;) {
470 q--;
471 if (q == offset) {
472 break;
473 }
474 if (striped[q] == '0') {
475 significantDigits--;
476 } else {
477 break;
478 }
479 }
480 }
481
482 // special case of numbers like "0.00000"
483 if (decimalFound && significantDigits == 0) {
484 decimalPos = 0;
485 }
486
487 // Implicit decimal point at end of number if not present
488 if (!decimalFound) {
489 decimalPos = q-offset;
490 }
491
492 // Find the number of significant trailing zeros
493 q = offset; // set q to point to first sig digit
494 p = significantDigits-1+offset;
495
496 while (p > q) {
497 if (striped[p] != '0') {
498 break;
499 }
500 p--;
501 }
502
503 // Make sure the decimal is on a mod 10000 boundary
504 int i = ((rsize * 100) - decimalPos - sciexp % rsize) % rsize;
505 q -= i;
506 decimalPos += i;
507
508 // Make the mantissa length right by adding zeros at the end if necessary
509 while ((p - q) < (mant.length * rsize)) {
510 for (i = 0; i < rsize; i++) {
511 striped[++p] = '0';
512 }
513 }
514
515 // Ok, now we know how many trailing zeros there are,
516 // and where the least significant digit is
517 for (i = mant.length - 1; i >= 0; i--) {
518 mant[i] = (striped[q] - '0') * 1000 +
519 (striped[q+1] - '0') * 100 +
520 (striped[q+2] - '0') * 10 +
521 (striped[q+3] - '0');
522 q += 4;
523 }
524
525
526 exp = (decimalPos+sciexp) / rsize;
527
528 if (q < striped.length) {
529 // Is there possible another digit?
530 round((striped[q] - '0')*1000);
531 }
532
533 }
534
535 /** Creates an instance with a non-finite value.
536 * @param field field to which this instance belongs
537 * @param sign sign of the Dfp to create
538 * @param nans code of the value, must be one of {@link #INFINITE},
539 * {@link #SNAN}, {@link #QNAN}
540 */
541 protected Dfp(final DfpField field, final byte sign, final byte nans) {
542 this.field = field;
543 this.mant = new int[field.getRadixDigits()];
544 this.sign = sign;
545 this.exp = 0;
546 this.nans = nans;
547 }
548
549 /** Create an instance with a value of 0.
550 * Use this internally in preference to constructors to facilitate subclasses
551 * @return a new instance with a value of 0
552 */
553 public Dfp newInstance() {
554 return new Dfp(getField());
555 }
556
557 /** Create an instance from a byte value.
558 * @param x value to convert to an instance
559 * @return a new instance with value x
560 */
561 public Dfp newInstance(final byte x) {
562 return new Dfp(getField(), x);
563 }
564
565 /** Create an instance from an int value.
566 * @param x value to convert to an instance
567 * @return a new instance with value x
568 */
569 public Dfp newInstance(final int x) {
570 return new Dfp(getField(), x);
571 }
572
573 /** Create an instance from a long value.
574 * @param x value to convert to an instance
575 * @return a new instance with value x
576 */
577 public Dfp newInstance(final long x) {
578 return new Dfp(getField(), x);
579 }
580
581 /** Create an instance from a double value.
582 * @param x value to convert to an instance
583 * @return a new instance with value x
584 */
585 public Dfp newInstance(final double x) {
586 return new Dfp(getField(), x);
587 }
588
589 /** Create an instance by copying an existing one.
590 * Use this internally in preference to constructors to facilitate subclasses.
591 * @param d instance to copy
592 * @return a new instance with the same value as d
593 */
594 public Dfp newInstance(final Dfp d) {
595
596 // make sure we don't mix number with different precision
597 if (field.getRadixDigits() != d.field.getRadixDigits()) {
598 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
599 final Dfp result = newInstance(getZero());
600 result.nans = QNAN;
601 return dotrap(DfpField.FLAG_INVALID, NEW_INSTANCE_TRAP, d, result);
602 }
603
604 return new Dfp(d);
605
606 }
607
608 /** Create an instance from a String representation.
609 * Use this internally in preference to constructors to facilitate subclasses.
610 * @param s string representation of the instance
611 * @return a new instance parsed from specified string
612 */
613 public Dfp newInstance(final String s) {
614 return new Dfp(field, s);
615 }
616
617 /** Creates an instance with a non-finite value.
618 * @param sig sign of the Dfp to create
619 * @param code code of the value, must be one of {@link #INFINITE},
620 * {@link #SNAN}, {@link #QNAN}
621 * @return a new instance with a non-finite value
622 */
623 public Dfp newInstance(final byte sig, final byte code) {
624 return field.newDfp(sig, code);
625 }
626
627 /** Get the {@link org.apache.commons.math3.Field Field} (really a {@link DfpField}) to which the instance belongs.
628 * <p>
629 * The field is linked to the number of digits and acts as a factory
630 * for {@link Dfp} instances.
631 * </p>
632 * @return {@link org.apache.commons.math3.Field Field} (really a {@link DfpField}) to which the instance belongs
633 */
634 public DfpField getField() {
635 return field;
636 }
637
638 /** Get the number of radix digits of the instance.
639 * @return number of radix digits
640 */
641 public int getRadixDigits() {
642 return field.getRadixDigits();
643 }
644
645 /** Get the constant 0.
646 * @return a Dfp with value zero
647 */
648 public Dfp getZero() {
649 return field.getZero();
650 }
651
652 /** Get the constant 1.
653 * @return a Dfp with value one
654 */
655 public Dfp getOne() {
656 return field.getOne();
657 }
658
659 /** Get the constant 2.
660 * @return a Dfp with value two
661 */
662 public Dfp getTwo() {
663 return field.getTwo();
664 }
665
666 /** Shift the mantissa left, and adjust the exponent to compensate.
667 */
668 protected void shiftLeft() {
669 for (int i = mant.length - 1; i > 0; i--) {
670 mant[i] = mant[i-1];
671 }
672 mant[0] = 0;
673 exp--;
674 }
675
676 /* Note that shiftRight() does not call round() as that round() itself
677 uses shiftRight() */
678 /** Shift the mantissa right, and adjust the exponent to compensate.
679 */
680 protected void shiftRight() {
681 for (int i = 0; i < mant.length - 1; i++) {
682 mant[i] = mant[i+1];
683 }
684 mant[mant.length - 1] = 0;
685 exp++;
686 }
687
688 /** Make our exp equal to the supplied one, this may cause rounding.
689 * Also causes de-normalized numbers. These numbers are generally
690 * dangerous because most routines assume normalized numbers.
691 * Align doesn't round, so it will return the last digit destroyed
692 * by shifting right.
693 * @param e desired exponent
694 * @return last digit destroyed by shifting right
695 */
696 protected int align(int e) {
697 int lostdigit = 0;
698 boolean inexact = false;
699
700 int diff = exp - e;
701
702 int adiff = diff;
703 if (adiff < 0) {
704 adiff = -adiff;
705 }
706
707 if (diff == 0) {
708 return 0;
709 }
710
711 if (adiff > (mant.length + 1)) {
712 // Special case
713 Arrays.fill(mant, 0);
714 exp = e;
715
716 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
717 dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
718
719 return 0;
720 }
721
722 for (int i = 0; i < adiff; i++) {
723 if (diff < 0) {
724 /* Keep track of loss -- only signal inexact after losing 2 digits.
725 * the first lost digit is returned to add() and may be incorporated
726 * into the result.
727 */
728 if (lostdigit != 0) {
729 inexact = true;
730 }
731
732 lostdigit = mant[0];
733
734 shiftRight();
735 } else {
736 shiftLeft();
737 }
738 }
739
740 if (inexact) {
741 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
742 dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
743 }
744
745 return lostdigit;
746
747 }
748
749 /** Check if instance is less than x.
750 * @param x number to check instance against
751 * @return true if instance is less than x and neither are NaN, false otherwise
752 */
753 public boolean lessThan(final Dfp x) {
754
755 // make sure we don't mix number with different precision
756 if (field.getRadixDigits() != x.field.getRadixDigits()) {
757 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
758 final Dfp result = newInstance(getZero());
759 result.nans = QNAN;
760 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, result);
761 return false;
762 }
763
764 /* if a nan is involved, signal invalid and return false */
765 if (isNaN() || x.isNaN()) {
766 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
767 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, newInstance(getZero()));
768 return false;
769 }
770
771 return compare(this, x) < 0;
772 }
773
774 /** Check if instance is greater than x.
775 * @param x number to check instance against
776 * @return true if instance is greater than x and neither are NaN, false otherwise
777 */
778 public boolean greaterThan(final Dfp x) {
779
780 // make sure we don't mix number with different precision
781 if (field.getRadixDigits() != x.field.getRadixDigits()) {
782 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
783 final Dfp result = newInstance(getZero());
784 result.nans = QNAN;
785 dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, result);
786 return false;
787 }
788
789 /* if a nan is involved, signal invalid and return false */
790 if (isNaN() || x.isNaN()) {
791 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
792 dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, newInstance(getZero()));
793 return false;
794 }
795
796 return compare(this, x) > 0;
797 }
798
799 /** Check if instance is less than or equal to 0.
800 * @return true if instance is not NaN and less than or equal to 0, false otherwise
801 */
802 public boolean negativeOrNull() {
803
804 if (isNaN()) {
805 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
806 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
807 return false;
808 }
809
810 return (sign < 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
811
812 }
813
814 /** Check if instance is strictly less than 0.
815 * @return true if instance is not NaN and less than or equal to 0, false otherwise
816 */
817 public boolean strictlyNegative() {
818
819 if (isNaN()) {
820 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
821 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
822 return false;
823 }
824
825 return (sign < 0) && ((mant[mant.length - 1] != 0) || isInfinite());
826
827 }
828
829 /** Check if instance is greater than or equal to 0.
830 * @return true if instance is not NaN and greater than or equal to 0, false otherwise
831 */
832 public boolean positiveOrNull() {
833
834 if (isNaN()) {
835 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
836 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
837 return false;
838 }
839
840 return (sign > 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
841
842 }
843
844 /** Check if instance is strictly greater than 0.
845 * @return true if instance is not NaN and greater than or equal to 0, false otherwise
846 */
847 public boolean strictlyPositive() {
848
849 if (isNaN()) {
850 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
851 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
852 return false;
853 }
854
855 return (sign > 0) && ((mant[mant.length - 1] != 0) || isInfinite());
856
857 }
858
859 /** Get the absolute value of instance.
860 * @return absolute value of instance
861 */
862 public Dfp abs() {
863 Dfp result = newInstance(this);
864 result.sign = 1;
865 return result;
866 }
867
868 /** Check if instance is infinite.
869 * @return true if instance is infinite
870 */
871 public boolean isInfinite() {
872 return nans == INFINITE;
873 }
874
875 /** Check if instance is not a number.
876 * @return true if instance is not a number
877 */
878 public boolean isNaN() {
879 return (nans == QNAN) || (nans == SNAN);
880 }
881
882 /** Check if instance is equal to zero.
883 * @return true if instance is equal to zero
884 */
885 public boolean isZero() {
886
887 if (isNaN()) {
888 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
889 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
890 return false;
891 }
892
893 return (mant[mant.length - 1] == 0) && !isInfinite();
894
895 }
896
897 /** Check if instance is equal to x.
898 * @param other object to check instance against
899 * @return true if instance is equal to x and neither are NaN, false otherwise
900 */
901 @Override
902 public boolean equals(final Object other) {
903
904 if (other instanceof Dfp) {
905 final Dfp x = (Dfp) other;
906 if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
907 return false;
908 }
909
910 return compare(this, x) == 0;
911 }
912
913 return false;
914
915 }
916
917 /**
918 * Gets a hashCode for the instance.
919 * @return a hash code value for this object
920 */
921 @Override
922 public int hashCode() {
923 return 17 + (sign << 8) + (nans << 16) + exp + Arrays.hashCode(mant);
924 }
925
926 /** Check if instance is not equal to x.
927 * @param x number to check instance against
928 * @return true if instance is not equal to x and neither are NaN, false otherwise
929 */
930 public boolean unequal(final Dfp x) {
931 if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
932 return false;
933 }
934
935 return greaterThan(x) || lessThan(x);
936 }
937
938 /** Compare two instances.
939 * @param a first instance in comparison
940 * @param b second instance in comparison
941 * @return -1 if a<b, 1 if a>b and 0 if a==b
942 * Note this method does not properly handle NaNs or numbers with different precision.
943 */
944 private static int compare(final Dfp a, final Dfp b) {
945 // Ignore the sign of zero
946 if (a.mant[a.mant.length - 1] == 0 && b.mant[b.mant.length - 1] == 0 &&
947 a.nans == FINITE && b.nans == FINITE) {
948 return 0;
949 }
950
951 if (a.sign != b.sign) {
952 if (a.sign == -1) {
953 return -1;
954 } else {
955 return 1;
956 }
957 }
958
959 // deal with the infinities
960 if (a.nans == INFINITE && b.nans == FINITE) {
961 return a.sign;
962 }
963
964 if (a.nans == FINITE && b.nans == INFINITE) {
965 return -b.sign;
966 }
967
968 if (a.nans == INFINITE && b.nans == INFINITE) {
969 return 0;
970 }
971
972 // Handle special case when a or b is zero, by ignoring the exponents
973 if (b.mant[b.mant.length-1] != 0 && a.mant[b.mant.length-1] != 0) {
974 if (a.exp < b.exp) {
975 return -a.sign;
976 }
977
978 if (a.exp > b.exp) {
979 return a.sign;
980 }
981 }
982
983 // compare the mantissas
984 for (int i = a.mant.length - 1; i >= 0; i--) {
985 if (a.mant[i] > b.mant[i]) {
986 return a.sign;
987 }
988
989 if (a.mant[i] < b.mant[i]) {
990 return -a.sign;
991 }
992 }
993
994 return 0;
995
996 }
997
998 /** Round to nearest integer using the round-half-even method.
999 * That is round to nearest integer unless both are equidistant.
1000 * In which case round to the even one.
1001 * @return rounded value
1002 */
1003 public Dfp rint() {
1004 return trunc(DfpField.RoundingMode.ROUND_HALF_EVEN);
1005 }
1006
1007 /** Round to an integer using the round floor mode.
1008 * That is, round toward -Infinity
1009 * @return rounded value
1010 */
1011 public Dfp floor() {
1012 return trunc(DfpField.RoundingMode.ROUND_FLOOR);
1013 }
1014
1015 /** Round to an integer using the round ceil mode.
1016 * That is, round toward +Infinity
1017 * @return rounded value
1018 */
1019 public Dfp ceil() {
1020 return trunc(DfpField.RoundingMode.ROUND_CEIL);
1021 }
1022
1023 /** Returns the IEEE remainder.
1024 * @param d divisor
1025 * @return this less n × d, where n is the integer closest to this/d
1026 */
1027 public Dfp remainder(final Dfp d) {
1028
1029 final Dfp result = this.subtract(this.divide(d).rint().multiply(d));
1030
1031 // IEEE 854-1987 says that if the result is zero, then it carries the sign of this
1032 if (result.mant[mant.length-1] == 0) {
1033 result.sign = sign;
1034 }
1035
1036 return result;
1037
1038 }
1039
1040 /** Does the integer conversions with the specified rounding.
1041 * @param rmode rounding mode to use
1042 * @return truncated value
1043 */
1044 protected Dfp trunc(final DfpField.RoundingMode rmode) {
1045 boolean changed = false;
1046
1047 if (isNaN()) {
1048 return newInstance(this);
1049 }
1050
1051 if (nans == INFINITE) {
1052 return newInstance(this);
1053 }
1054
1055 if (mant[mant.length-1] == 0) {
1056 // a is zero
1057 return newInstance(this);
1058 }
1059
1060 /* If the exponent is less than zero then we can certainly
1061 * return zero */
1062 if (exp < 0) {
1063 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
1064 Dfp result = newInstance(getZero());
1065 result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
1066 return result;
1067 }
1068
1069 /* If the exponent is greater than or equal to digits, then it
1070 * must already be an integer since there is no precision left
1071 * for any fractional part */
1072
1073 if (exp >= mant.length) {
1074 return newInstance(this);
1075 }
1076
1077 /* General case: create another dfp, result, that contains the
1078 * a with the fractional part lopped off. */
1079
1080 Dfp result = newInstance(this);
1081 for (int i = 0; i < mant.length-result.exp; i++) {
1082 changed |= result.mant[i] != 0;
1083 result.mant[i] = 0;
1084 }
1085
1086 if (changed) {
1087 switch (rmode) {
1088 case ROUND_FLOOR:
1089 if (result.sign == -1) {
1090 // then we must increment the mantissa by one
1091 result = result.add(newInstance(-1));
1092 }
1093 break;
1094
1095 case ROUND_CEIL:
1096 if (result.sign == 1) {
1097 // then we must increment the mantissa by one
1098 result = result.add(getOne());
1099 }
1100 break;
1101
1102 case ROUND_HALF_EVEN:
1103 default:
1104 final Dfp half = newInstance("0.5");
1105 Dfp a = subtract(result); // difference between this and result
1106 a.sign = 1; // force positive (take abs)
1107 if (a.greaterThan(half)) {
1108 a = newInstance(getOne());
1109 a.sign = sign;
1110 result = result.add(a);
1111 }
1112
1113 /** If exactly equal to 1/2 and odd then increment */
1114 if (a.equals(half) && result.exp > 0 && (result.mant[mant.length-result.exp]&1) != 0) {
1115 a = newInstance(getOne());
1116 a.sign = sign;
1117 result = result.add(a);
1118 }
1119 break;
1120 }
1121
1122 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT); // signal inexact
1123 result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
1124 return result;
1125 }
1126
1127 return result;
1128 }
1129
1130 /** Convert this to an integer.
1131 * If greater than 2147483647, it returns 2147483647. If less than -2147483648 it returns -2147483648.
1132 * @return converted number
1133 */
1134 public int intValue() {
1135 Dfp rounded;
1136 int result = 0;
1137
1138 rounded = rint();
1139
1140 if (rounded.greaterThan(newInstance(2147483647))) {
1141 return 2147483647;
1142 }
1143
1144 if (rounded.lessThan(newInstance(-2147483648))) {
1145 return -2147483648;
1146 }
1147
1148 for (int i = mant.length - 1; i >= mant.length - rounded.exp; i--) {
1149 result = result * RADIX + rounded.mant[i];
1150 }
1151
1152 if (rounded.sign == -1) {
1153 result = -result;
1154 }
1155
1156 return result;
1157 }
1158
1159 /** Get the exponent of the greatest power of 10000 that is
1160 * less than or equal to the absolute value of this. I.E. if
1161 * this is 10<sup>6</sup> then log10K would return 1.
1162 * @return integer base 10000 logarithm
1163 */
1164 public int log10K() {
1165 return exp - 1;
1166 }
1167
1168 /** Get the specified power of 10000.
1169 * @param e desired power
1170 * @return 10000<sup>e</sup>
1171 */
1172 public Dfp power10K(final int e) {
1173 Dfp d = newInstance(getOne());
1174 d.exp = e + 1;
1175 return d;
1176 }
1177
1178 /** Get the exponent of the greatest power of 10 that is less than or equal to abs(this).
1179 * @return integer base 10 logarithm
1180 */
1181 public int log10() {
1182 if (mant[mant.length-1] > 1000) {
1183 return exp * 4 - 1;
1184 }
1185 if (mant[mant.length-1] > 100) {
1186 return exp * 4 - 2;
1187 }
1188 if (mant[mant.length-1] > 10) {
1189 return exp * 4 - 3;
1190 }
1191 return exp * 4 - 4;
1192 }
1193
1194 /** Return the specified power of 10.
1195 * @param e desired power
1196 * @return 10<sup>e</sup>
1197 */
1198 public Dfp power10(final int e) {
1199 Dfp d = newInstance(getOne());
1200
1201 if (e >= 0) {
1202 d.exp = e / 4 + 1;
1203 } else {
1204 d.exp = (e + 1) / 4;
1205 }
1206
1207 switch ((e % 4 + 4) % 4) {
1208 case 0:
1209 break;
1210 case 1:
1211 d = d.multiply(10);
1212 break;
1213 case 2:
1214 d = d.multiply(100);
1215 break;
1216 default:
1217 d = d.multiply(1000);
1218 }
1219
1220 return d;
1221 }
1222
1223 /** Negate the mantissa of this by computing the complement.
1224 * Leaves the sign bit unchanged, used internally by add.
1225 * Denormalized numbers are handled properly here.
1226 * @param extra ???
1227 * @return ???
1228 */
1229 protected int complement(int extra) {
1230
1231 extra = RADIX-extra;
1232 for (int i = 0; i < mant.length; i++) {
1233 mant[i] = RADIX-mant[i]-1;
1234 }
1235
1236 int rh = extra / RADIX;
1237 extra = extra - rh * RADIX;
1238 for (int i = 0; i < mant.length; i++) {
1239 final int r = mant[i] + rh;
1240 rh = r / RADIX;
1241 mant[i] = r - rh * RADIX;
1242 }
1243
1244 return extra;
1245 }
1246
1247 /** Add x to this.
1248 * @param x number to add
1249 * @return sum of this and x
1250 */
1251 public Dfp add(final Dfp x) {
1252
1253 // make sure we don't mix number with different precision
1254 if (field.getRadixDigits() != x.field.getRadixDigits()) {
1255 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1256 final Dfp result = newInstance(getZero());
1257 result.nans = QNAN;
1258 return dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
1259 }
1260
1261 /* handle special cases */
1262 if (nans != FINITE || x.nans != FINITE) {
1263 if (isNaN()) {
1264 return this;
1265 }
1266
1267 if (x.isNaN()) {
1268 return x;
1269 }
1270
1271 if (nans == INFINITE && x.nans == FINITE) {
1272 return this;
1273 }
1274
1275 if (x.nans == INFINITE && nans == FINITE) {
1276 return x;
1277 }
1278
1279 if (x.nans == INFINITE && nans == INFINITE && sign == x.sign) {
1280 return x;
1281 }
1282
1283 if (x.nans == INFINITE && nans == INFINITE && sign != x.sign) {
1284 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1285 Dfp result = newInstance(getZero());
1286 result.nans = QNAN;
1287 result = dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
1288 return result;
1289 }
1290 }
1291
1292 /* copy this and the arg */
1293 Dfp a = newInstance(this);
1294 Dfp b = newInstance(x);
1295
1296 /* initialize the result object */
1297 Dfp result = newInstance(getZero());
1298
1299 /* Make all numbers positive, but remember their sign */
1300 final byte asign = a.sign;
1301 final byte bsign = b.sign;
1302
1303 a.sign = 1;
1304 b.sign = 1;
1305
1306 /* The result will be signed like the arg with greatest magnitude */
1307 byte rsign = bsign;
1308 if (compare(a, b) > 0) {
1309 rsign = asign;
1310 }
1311
1312 /* Handle special case when a or b is zero, by setting the exponent
1313 of the zero number equal to the other one. This avoids an alignment
1314 which would cause catastropic loss of precision */
1315 if (b.mant[mant.length-1] == 0) {
1316 b.exp = a.exp;
1317 }
1318
1319 if (a.mant[mant.length-1] == 0) {
1320 a.exp = b.exp;
1321 }
1322
1323 /* align number with the smaller exponent */
1324 int aextradigit = 0;
1325 int bextradigit = 0;
1326 if (a.exp < b.exp) {
1327 aextradigit = a.align(b.exp);
1328 } else {
1329 bextradigit = b.align(a.exp);
1330 }
1331
1332 /* complement the smaller of the two if the signs are different */
1333 if (asign != bsign) {
1334 if (asign == rsign) {
1335 bextradigit = b.complement(bextradigit);
1336 } else {
1337 aextradigit = a.complement(aextradigit);
1338 }
1339 }
1340
1341 /* add the mantissas */
1342 int rh = 0; /* acts as a carry */
1343 for (int i = 0; i < mant.length; i++) {
1344 final int r = a.mant[i]+b.mant[i]+rh;
1345 rh = r / RADIX;
1346 result.mant[i] = r - rh * RADIX;
1347 }
1348 result.exp = a.exp;
1349 result.sign = rsign;
1350
1351 /* handle overflow -- note, when asign!=bsign an overflow is
1352 * normal and should be ignored. */
1353
1354 if (rh != 0 && (asign == bsign)) {
1355 final int lostdigit = result.mant[0];
1356 result.shiftRight();
1357 result.mant[mant.length-1] = rh;
1358 final int excp = result.round(lostdigit);
1359 if (excp != 0) {
1360 result = dotrap(excp, ADD_TRAP, x, result);
1361 }
1362 }
1363
1364 /* normalize the result */
1365 for (int i = 0; i < mant.length; i++) {
1366 if (result.mant[mant.length-1] != 0) {
1367 break;
1368 }
1369 result.shiftLeft();
1370 if (i == 0) {
1371 result.mant[0] = aextradigit+bextradigit;
1372 aextradigit = 0;
1373 bextradigit = 0;
1374 }
1375 }
1376
1377 /* result is zero if after normalization the most sig. digit is zero */
1378 if (result.mant[mant.length-1] == 0) {
1379 result.exp = 0;
1380
1381 if (asign != bsign) {
1382 // Unless adding 2 negative zeros, sign is positive
1383 result.sign = 1; // Per IEEE 854-1987 Section 6.3
1384 }
1385 }
1386
1387 /* Call round to test for over/under flows */
1388 final int excp = result.round(aextradigit + bextradigit);
1389 if (excp != 0) {
1390 result = dotrap(excp, ADD_TRAP, x, result);
1391 }
1392
1393 return result;
1394 }
1395
1396 /** Returns a number that is this number with the sign bit reversed.
1397 * @return the opposite of this
1398 */
1399 public Dfp negate() {
1400 Dfp result = newInstance(this);
1401 result.sign = (byte) - result.sign;
1402 return result;
1403 }
1404
1405 /** Subtract x from this.
1406 * @param x number to subtract
1407 * @return difference of this and a
1408 */
1409 public Dfp subtract(final Dfp x) {
1410 return add(x.negate());
1411 }
1412
1413 /** Round this given the next digit n using the current rounding mode.
1414 * @param n ???
1415 * @return the IEEE flag if an exception occurred
1416 */
1417 protected int round(int n) {
1418 boolean inc = false;
1419 switch (field.getRoundingMode()) {
1420 case ROUND_DOWN:
1421 inc = false;
1422 break;
1423
1424 case ROUND_UP:
1425 inc = n != 0; // round up if n!=0
1426 break;
1427
1428 case ROUND_HALF_UP:
1429 inc = n >= 5000; // round half up
1430 break;
1431
1432 case ROUND_HALF_DOWN:
1433 inc = n > 5000; // round half down
1434 break;
1435
1436 case ROUND_HALF_EVEN:
1437 inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 1); // round half-even
1438 break;
1439
1440 case ROUND_HALF_ODD:
1441 inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 0); // round half-odd
1442 break;
1443
1444 case ROUND_CEIL:
1445 inc = sign == 1 && n != 0; // round ceil
1446 break;
1447
1448 case ROUND_FLOOR:
1449 default:
1450 inc = sign == -1 && n != 0; // round floor
1451 break;
1452 }
1453
1454 if (inc) {
1455 // increment if necessary
1456 int rh = 1;
1457 for (int i = 0; i < mant.length; i++) {
1458 final int r = mant[i] + rh;
1459 rh = r / RADIX;
1460 mant[i] = r - rh * RADIX;
1461 }
1462
1463 if (rh != 0) {
1464 shiftRight();
1465 mant[mant.length-1] = rh;
1466 }
1467 }
1468
1469 // check for exceptional cases and raise signals if necessary
1470 if (exp < MIN_EXP) {
1471 // Gradual Underflow
1472 field.setIEEEFlagsBits(DfpField.FLAG_UNDERFLOW);
1473 return DfpField.FLAG_UNDERFLOW;
1474 }
1475
1476 if (exp > MAX_EXP) {
1477 // Overflow
1478 field.setIEEEFlagsBits(DfpField.FLAG_OVERFLOW);
1479 return DfpField.FLAG_OVERFLOW;
1480 }
1481
1482 if (n != 0) {
1483 // Inexact
1484 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
1485 return DfpField.FLAG_INEXACT;
1486 }
1487
1488 return 0;
1489
1490 }
1491
1492 /** Multiply this by x.
1493 * @param x multiplicand
1494 * @return product of this and x
1495 */
1496 public Dfp multiply(final Dfp x) {
1497
1498 // make sure we don't mix number with different precision
1499 if (field.getRadixDigits() != x.field.getRadixDigits()) {
1500 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1501 final Dfp result = newInstance(getZero());
1502 result.nans = QNAN;
1503 return dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
1504 }
1505
1506 Dfp result = newInstance(getZero());
1507
1508 /* handle special cases */
1509 if (nans != FINITE || x.nans != FINITE) {
1510 if (isNaN()) {
1511 return this;
1512 }
1513
1514 if (x.isNaN()) {
1515 return x;
1516 }
1517
1518 if (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] != 0) {
1519 result = newInstance(this);
1520 result.sign = (byte) (sign * x.sign);
1521 return result;
1522 }
1523
1524 if (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] != 0) {
1525 result = newInstance(x);
1526 result.sign = (byte) (sign * x.sign);
1527 return result;
1528 }
1529
1530 if (x.nans == INFINITE && nans == INFINITE) {
1531 result = newInstance(this);
1532 result.sign = (byte) (sign * x.sign);
1533 return result;
1534 }
1535
1536 if ( (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] == 0) ||
1537 (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] == 0) ) {
1538 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1539 result = newInstance(getZero());
1540 result.nans = QNAN;
1541 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
1542 return result;
1543 }
1544 }
1545
1546 int[] product = new int[mant.length*2]; // Big enough to hold even the largest result
1547
1548 for (int i = 0; i < mant.length; i++) {
1549 int rh = 0; // acts as a carry
1550 for (int j=0; j<mant.length; j++) {
1551 int r = mant[i] * x.mant[j]; // multiply the 2 digits
1552 r = r + product[i+j] + rh; // add to the product digit with carry in
1553
1554 rh = r / RADIX;
1555 product[i+j] = r - rh * RADIX;
1556 }
1557 product[i+mant.length] = rh;
1558 }
1559
1560 // Find the most sig digit
1561 int md = mant.length * 2 - 1; // default, in case result is zero
1562 for (int i = mant.length * 2 - 1; i >= 0; i--) {
1563 if (product[i] != 0) {
1564 md = i;
1565 break;
1566 }
1567 }
1568
1569 // Copy the digits into the result
1570 for (int i = 0; i < mant.length; i++) {
1571 result.mant[mant.length - i - 1] = product[md - i];
1572 }
1573
1574 // Fixup the exponent.
1575 result.exp = exp + x.exp + md - 2 * mant.length + 1;
1576 result.sign = (byte)((sign == x.sign)?1:-1);
1577
1578 if (result.mant[mant.length-1] == 0) {
1579 // if result is zero, set exp to zero
1580 result.exp = 0;
1581 }
1582
1583 final int excp;
1584 if (md > (mant.length-1)) {
1585 excp = result.round(product[md-mant.length]);
1586 } else {
1587 excp = result.round(0); // has no effect except to check status
1588 }
1589
1590 if (excp != 0) {
1591 result = dotrap(excp, MULTIPLY_TRAP, x, result);
1592 }
1593
1594 return result;
1595
1596 }
1597
1598 /** Multiply this by a single digit x.
1599 * @param x multiplicand
1600 * @return product of this and x
1601 */
1602 public Dfp multiply(final int x) {
1603 if (x >= 0 && x < RADIX) {
1604 return multiplyFast(x);
1605 } else {
1606 return multiply(newInstance(x));
1607 }
1608 }
1609
1610 /** Multiply this by a single digit 0<=x<radix.
1611 * There are speed advantages in this special case.
1612 * @param x multiplicand
1613 * @return product of this and x
1614 */
1615 private Dfp multiplyFast(final int x) {
1616 Dfp result = newInstance(this);
1617
1618 /* handle special cases */
1619 if (nans != FINITE) {
1620 if (isNaN()) {
1621 return this;
1622 }
1623
1624 if (nans == INFINITE && x != 0) {
1625 result = newInstance(this);
1626 return result;
1627 }
1628
1629 if (nans == INFINITE && x == 0) {
1630 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1631 result = newInstance(getZero());
1632 result.nans = QNAN;
1633 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, newInstance(getZero()), result);
1634 return result;
1635 }
1636 }
1637
1638 /* range check x */
1639 if (x < 0 || x >= RADIX) {
1640 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1641 result = newInstance(getZero());
1642 result.nans = QNAN;
1643 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, result, result);
1644 return result;
1645 }
1646
1647 int rh = 0;
1648 for (int i = 0; i < mant.length; i++) {
1649 final int r = mant[i] * x + rh;
1650 rh = r / RADIX;
1651 result.mant[i] = r - rh * RADIX;
1652 }
1653
1654 int lostdigit = 0;
1655 if (rh != 0) {
1656 lostdigit = result.mant[0];
1657 result.shiftRight();
1658 result.mant[mant.length-1] = rh;
1659 }
1660
1661 if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
1662 result.exp = 0;
1663 }
1664
1665 final int excp = result.round(lostdigit);
1666 if (excp != 0) {
1667 result = dotrap(excp, MULTIPLY_TRAP, result, result);
1668 }
1669
1670 return result;
1671 }
1672
1673 /** Divide this by divisor.
1674 * @param divisor divisor
1675 * @return quotient of this by divisor
1676 */
1677 public Dfp divide(Dfp divisor) {
1678 int dividend[]; // current status of the dividend
1679 int quotient[]; // quotient
1680 int remainder[];// remainder
1681 int qd; // current quotient digit we're working with
1682 int nsqd; // number of significant quotient digits we have
1683 int trial=0; // trial quotient digit
1684 int minadj; // minimum adjustment
1685 boolean trialgood; // Flag to indicate a good trail digit
1686 int md=0; // most sig digit in result
1687 int excp; // exceptions
1688
1689 // make sure we don't mix number with different precision
1690 if (field.getRadixDigits() != divisor.field.getRadixDigits()) {
1691 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1692 final Dfp result = newInstance(getZero());
1693 result.nans = QNAN;
1694 return dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
1695 }
1696
1697 Dfp result = newInstance(getZero());
1698
1699 /* handle special cases */
1700 if (nans != FINITE || divisor.nans != FINITE) {
1701 if (isNaN()) {
1702 return this;
1703 }
1704
1705 if (divisor.isNaN()) {
1706 return divisor;
1707 }
1708
1709 if (nans == INFINITE && divisor.nans == FINITE) {
1710 result = newInstance(this);
1711 result.sign = (byte) (sign * divisor.sign);
1712 return result;
1713 }
1714
1715 if (divisor.nans == INFINITE && nans == FINITE) {
1716 result = newInstance(getZero());
1717 result.sign = (byte) (sign * divisor.sign);
1718 return result;
1719 }
1720
1721 if (divisor.nans == INFINITE && nans == INFINITE) {
1722 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1723 result = newInstance(getZero());
1724 result.nans = QNAN;
1725 result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
1726 return result;
1727 }
1728 }
1729
1730 /* Test for divide by zero */
1731 if (divisor.mant[mant.length-1] == 0) {
1732 field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
1733 result = newInstance(getZero());
1734 result.sign = (byte) (sign * divisor.sign);
1735 result.nans = INFINITE;
1736 result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, divisor, result);
1737 return result;
1738 }
1739
1740 dividend = new int[mant.length+1]; // one extra digit needed
1741 quotient = new int[mant.length+2]; // two extra digits needed 1 for overflow, 1 for rounding
1742 remainder = new int[mant.length+1]; // one extra digit needed
1743
1744 /* Initialize our most significant digits to zero */
1745
1746 dividend[mant.length] = 0;
1747 quotient[mant.length] = 0;
1748 quotient[mant.length+1] = 0;
1749 remainder[mant.length] = 0;
1750
1751 /* copy our mantissa into the dividend, initialize the
1752 quotient while we are at it */
1753
1754 for (int i = 0; i < mant.length; i++) {
1755 dividend[i] = mant[i];
1756 quotient[i] = 0;
1757 remainder[i] = 0;
1758 }
1759
1760 /* outer loop. Once per quotient digit */
1761 nsqd = 0;
1762 for (qd = mant.length+1; qd >= 0; qd--) {
1763 /* Determine outer limits of our quotient digit */
1764
1765 // r = most sig 2 digits of dividend
1766 final int divMsb = dividend[mant.length]*RADIX+dividend[mant.length-1];
1767 int min = divMsb / (divisor.mant[mant.length-1]+1);
1768 int max = (divMsb + 1) / divisor.mant[mant.length-1];
1769
1770 trialgood = false;
1771 while (!trialgood) {
1772 // try the mean
1773 trial = (min+max)/2;
1774
1775 /* Multiply by divisor and store as remainder */
1776 int rh = 0;
1777 for (int i = 0; i < mant.length + 1; i++) {
1778 int dm = (i<mant.length)?divisor.mant[i]:0;
1779 final int r = (dm * trial) + rh;
1780 rh = r / RADIX;
1781 remainder[i] = r - rh * RADIX;
1782 }
1783
1784 /* subtract the remainder from the dividend */
1785 rh = 1; // carry in to aid the subtraction
1786 for (int i = 0; i < mant.length + 1; i++) {
1787 final int r = ((RADIX-1) - remainder[i]) + dividend[i] + rh;
1788 rh = r / RADIX;
1789 remainder[i] = r - rh * RADIX;
1790 }
1791
1792 /* Lets analyze what we have here */
1793 if (rh == 0) {
1794 // trial is too big -- negative remainder
1795 max = trial-1;
1796 continue;
1797 }
1798
1799 /* find out how far off the remainder is telling us we are */
1800 minadj = (remainder[mant.length] * RADIX)+remainder[mant.length-1];
1801 minadj = minadj / (divisor.mant[mant.length-1]+1);
1802
1803 if (minadj >= 2) {
1804 min = trial+minadj; // update the minimum
1805 continue;
1806 }
1807
1808 /* May have a good one here, check more thoroughly. Basically
1809 its a good one if it is less than the divisor */
1810 trialgood = false; // assume false
1811 for (int i = mant.length - 1; i >= 0; i--) {
1812 if (divisor.mant[i] > remainder[i]) {
1813 trialgood = true;
1814 }
1815 if (divisor.mant[i] < remainder[i]) {
1816 break;
1817 }
1818 }
1819
1820 if (remainder[mant.length] != 0) {
1821 trialgood = false;
1822 }
1823
1824 if (trialgood == false) {
1825 min = trial+1;
1826 }
1827 }
1828
1829 /* Great we have a digit! */
1830 quotient[qd] = trial;
1831 if (trial != 0 || nsqd != 0) {
1832 nsqd++;
1833 }
1834
1835 if (field.getRoundingMode() == DfpField.RoundingMode.ROUND_DOWN && nsqd == mant.length) {
1836 // We have enough for this mode
1837 break;
1838 }
1839
1840 if (nsqd > mant.length) {
1841 // We have enough digits
1842 break;
1843 }
1844
1845 /* move the remainder into the dividend while left shifting */
1846 dividend[0] = 0;
1847 for (int i = 0; i < mant.length; i++) {
1848 dividend[i + 1] = remainder[i];
1849 }
1850 }
1851
1852 /* Find the most sig digit */
1853 md = mant.length; // default
1854 for (int i = mant.length + 1; i >= 0; i--) {
1855 if (quotient[i] != 0) {
1856 md = i;
1857 break;
1858 }
1859 }
1860
1861 /* Copy the digits into the result */
1862 for (int i=0; i<mant.length; i++) {
1863 result.mant[mant.length-i-1] = quotient[md-i];
1864 }
1865
1866 /* Fixup the exponent. */
1867 result.exp = exp - divisor.exp + md - mant.length;
1868 result.sign = (byte) ((sign == divisor.sign) ? 1 : -1);
1869
1870 if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
1871 result.exp = 0;
1872 }
1873
1874 if (md > (mant.length-1)) {
1875 excp = result.round(quotient[md-mant.length]);
1876 } else {
1877 excp = result.round(0);
1878 }
1879
1880 if (excp != 0) {
1881 result = dotrap(excp, DIVIDE_TRAP, divisor, result);
1882 }
1883
1884 return result;
1885 }
1886
1887 /** Divide by a single digit less than radix.
1888 * Special case, so there are speed advantages. 0 <= divisor < radix
1889 * @param divisor divisor
1890 * @return quotient of this by divisor
1891 */
1892 public Dfp divide(int divisor) {
1893
1894 // Handle special cases
1895 if (nans != FINITE) {
1896 if (isNaN()) {
1897 return this;
1898 }
1899
1900 if (nans == INFINITE) {
1901 return newInstance(this);
1902 }
1903 }
1904
1905 // Test for divide by zero
1906 if (divisor == 0) {
1907 field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
1908 Dfp result = newInstance(getZero());
1909 result.sign = sign;
1910 result.nans = INFINITE;
1911 result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, getZero(), result);
1912 return result;
1913 }
1914
1915 // range check divisor
1916 if (divisor < 0 || divisor >= RADIX) {
1917 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1918 Dfp result = newInstance(getZero());
1919 result.nans = QNAN;
1920 result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, result, result);
1921 return result;
1922 }
1923
1924 Dfp result = newInstance(this);
1925
1926 int rl = 0;
1927 for (int i = mant.length-1; i >= 0; i--) {
1928 final int r = rl*RADIX + result.mant[i];
1929 final int rh = r / divisor;
1930 rl = r - rh * divisor;
1931 result.mant[i] = rh;
1932 }
1933
1934 if (result.mant[mant.length-1] == 0) {
1935 // normalize
1936 result.shiftLeft();
1937 final int r = rl * RADIX; // compute the next digit and put it in
1938 final int rh = r / divisor;
1939 rl = r - rh * divisor;
1940 result.mant[0] = rh;
1941 }
1942
1943 final int excp = result.round(rl * RADIX / divisor); // do the rounding
1944 if (excp != 0) {
1945 result = dotrap(excp, DIVIDE_TRAP, result, result);
1946 }
1947
1948 return result;
1949
1950 }
1951
1952 /** {@inheritDoc} */
1953 public Dfp reciprocal() {
1954 return field.getOne().divide(this);
1955 }
1956
1957 /** Compute the square root.
1958 * @return square root of the instance
1959 */
1960 public Dfp sqrt() {
1961
1962 // check for unusual cases
1963 if (nans == FINITE && mant[mant.length-1] == 0) {
1964 // if zero
1965 return newInstance(this);
1966 }
1967
1968 if (nans != FINITE) {
1969 if (nans == INFINITE && sign == 1) {
1970 // if positive infinity
1971 return newInstance(this);
1972 }
1973
1974 if (nans == QNAN) {
1975 return newInstance(this);
1976 }
1977
1978 if (nans == SNAN) {
1979 Dfp result;
1980
1981 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1982 result = newInstance(this);
1983 result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
1984 return result;
1985 }
1986 }
1987
1988 if (sign == -1) {
1989 // if negative
1990 Dfp result;
1991
1992 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1993 result = newInstance(this);
1994 result.nans = QNAN;
1995 result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
1996 return result;
1997 }
1998
1999 Dfp x = newInstance(this);
2000
2001 /* Lets make a reasonable guess as to the size of the square root */
2002 if (x.exp < -1 || x.exp > 1) {
2003 x.exp = this.exp / 2;
2004 }
2005
2006 /* Coarsely estimate the mantissa */
2007 switch (x.mant[mant.length-1] / 2000) {
2008 case 0:
2009 x.mant[mant.length-1] = x.mant[mant.length-1]/2+1;
2010 break;
2011 case 2:
2012 x.mant[mant.length-1] = 1500;
2013 break;
2014 case 3:
2015 x.mant[mant.length-1] = 2200;
2016 break;
2017 default:
2018 x.mant[mant.length-1] = 3000;
2019 }
2020
2021 Dfp dx = newInstance(x);
2022
2023 /* Now that we have the first pass estimate, compute the rest
2024 by the formula dx = (y - x*x) / (2x); */
2025
2026 Dfp px = getZero();
2027 Dfp ppx = getZero();
2028 while (x.unequal(px)) {
2029 dx = newInstance(x);
2030 dx.sign = -1;
2031 dx = dx.add(this.divide(x));
2032 dx = dx.divide(2);
2033 ppx = px;
2034 px = x;
2035 x = x.add(dx);
2036
2037 if (x.equals(ppx)) {
2038 // alternating between two values
2039 break;
2040 }
2041
2042 // if dx is zero, break. Note testing the most sig digit
2043 // is a sufficient test since dx is normalized
2044 if (dx.mant[mant.length-1] == 0) {
2045 break;
2046 }
2047 }
2048
2049 return x;
2050
2051 }
2052
2053 /** Get a string representation of the instance.
2054 * @return string representation of the instance
2055 */
2056 @Override
2057 public String toString() {
2058 if (nans != FINITE) {
2059 // if non-finite exceptional cases
2060 if (nans == INFINITE) {
2061 return (sign < 0) ? NEG_INFINITY_STRING : POS_INFINITY_STRING;
2062 } else {
2063 return NAN_STRING;
2064 }
2065 }
2066
2067 if (exp > mant.length || exp < -1) {
2068 return dfp2sci();
2069 }
2070
2071 return dfp2string();
2072
2073 }
2074
2075 /** Convert an instance to a string using scientific notation.
2076 * @return string representation of the instance in scientific notation
2077 */
2078 protected String dfp2sci() {
2079 char rawdigits[] = new char[mant.length * 4];
2080 char outputbuffer[] = new char[mant.length * 4 + 20];
2081 int p;
2082 int q;
2083 int e;
2084 int ae;
2085 int shf;
2086
2087 // Get all the digits
2088 p = 0;
2089 for (int i = mant.length - 1; i >= 0; i--) {
2090 rawdigits[p++] = (char) ((mant[i] / 1000) + '0');
2091 rawdigits[p++] = (char) (((mant[i] / 100) %10) + '0');
2092 rawdigits[p++] = (char) (((mant[i] / 10) % 10) + '0');
2093 rawdigits[p++] = (char) (((mant[i]) % 10) + '0');
2094 }
2095
2096 // Find the first non-zero one
2097 for (p = 0; p < rawdigits.length; p++) {
2098 if (rawdigits[p] != '0') {
2099 break;
2100 }
2101 }
2102 shf = p;
2103
2104 // Now do the conversion
2105 q = 0;
2106 if (sign == -1) {
2107 outputbuffer[q++] = '-';
2108 }
2109
2110 if (p != rawdigits.length) {
2111 // there are non zero digits...
2112 outputbuffer[q++] = rawdigits[p++];
2113 outputbuffer[q++] = '.';
2114
2115 while (p<rawdigits.length) {
2116 outputbuffer[q++] = rawdigits[p++];
2117 }
2118 } else {
2119 outputbuffer[q++] = '0';
2120 outputbuffer[q++] = '.';
2121 outputbuffer[q++] = '0';
2122 outputbuffer[q++] = 'e';
2123 outputbuffer[q++] = '0';
2124 return new String(outputbuffer, 0, 5);
2125 }
2126
2127 outputbuffer[q++] = 'e';
2128
2129 // Find the msd of the exponent
2130
2131 e = exp * 4 - shf - 1;
2132 ae = e;
2133 if (e < 0) {
2134 ae = -e;
2135 }
2136
2137 // Find the largest p such that p < e
2138 for (p = 1000000000; p > ae; p /= 10) {
2139 // nothing to do
2140 }
2141
2142 if (e < 0) {
2143 outputbuffer[q++] = '-';
2144 }
2145
2146 while (p > 0) {
2147 outputbuffer[q++] = (char)(ae / p + '0');
2148 ae = ae % p;
2149 p = p / 10;
2150 }
2151
2152 return new String(outputbuffer, 0, q);
2153
2154 }
2155
2156 /** Convert an instance to a string using normal notation.
2157 * @return string representation of the instance in normal notation
2158 */
2159 protected String dfp2string() {
2160 char buffer[] = new char[mant.length*4 + 20];
2161 int p = 1;
2162 int q;
2163 int e = exp;
2164 boolean pointInserted = false;
2165
2166 buffer[0] = ' ';
2167
2168 if (e <= 0) {
2169 buffer[p++] = '0';
2170 buffer[p++] = '.';
2171 pointInserted = true;
2172 }
2173
2174 while (e < 0) {
2175 buffer[p++] = '0';
2176 buffer[p++] = '0';
2177 buffer[p++] = '0';
2178 buffer[p++] = '0';
2179 e++;
2180 }
2181
2182 for (int i = mant.length - 1; i >= 0; i--) {
2183 buffer[p++] = (char) ((mant[i] / 1000) + '0');
2184 buffer[p++] = (char) (((mant[i] / 100) % 10) + '0');
2185 buffer[p++] = (char) (((mant[i] / 10) % 10) + '0');
2186 buffer[p++] = (char) (((mant[i]) % 10) + '0');
2187 if (--e == 0) {
2188 buffer[p++] = '.';
2189 pointInserted = true;
2190 }
2191 }
2192
2193 while (e > 0) {
2194 buffer[p++] = '0';
2195 buffer[p++] = '0';
2196 buffer[p++] = '0';
2197 buffer[p++] = '0';
2198 e--;
2199 }
2200
2201 if (!pointInserted) {
2202 // Ensure we have a radix point!
2203 buffer[p++] = '.';
2204 }
2205
2206 // Suppress leading zeros
2207 q = 1;
2208 while (buffer[q] == '0') {
2209 q++;
2210 }
2211 if (buffer[q] == '.') {
2212 q--;
2213 }
2214
2215 // Suppress trailing zeros
2216 while (buffer[p-1] == '0') {
2217 p--;
2218 }
2219
2220 // Insert sign
2221 if (sign < 0) {
2222 buffer[--q] = '-';
2223 }
2224
2225 return new String(buffer, q, p - q);
2226
2227 }
2228
2229 /** Raises a trap. This does not set the corresponding flag however.
2230 * @param type the trap type
2231 * @param what - name of routine trap occurred in
2232 * @param oper - input operator to function
2233 * @param result - the result computed prior to the trap
2234 * @return The suggested return value from the trap handler
2235 */
2236 public Dfp dotrap(int type, String what, Dfp oper, Dfp result) {
2237 Dfp def = result;
2238
2239 switch (type) {
2240 case DfpField.FLAG_INVALID:
2241 def = newInstance(getZero());
2242 def.sign = result.sign;
2243 def.nans = QNAN;
2244 break;
2245
2246 case DfpField.FLAG_DIV_ZERO:
2247 if (nans == FINITE && mant[mant.length-1] != 0) {
2248 // normal case, we are finite, non-zero
2249 def = newInstance(getZero());
2250 def.sign = (byte)(sign*oper.sign);
2251 def.nans = INFINITE;
2252 }
2253
2254 if (nans == FINITE && mant[mant.length-1] == 0) {
2255 // 0/0
2256 def = newInstance(getZero());
2257 def.nans = QNAN;
2258 }
2259
2260 if (nans == INFINITE || nans == QNAN) {
2261 def = newInstance(getZero());
2262 def.nans = QNAN;
2263 }
2264
2265 if (nans == INFINITE || nans == SNAN) {
2266 def = newInstance(getZero());
2267 def.nans = QNAN;
2268 }
2269 break;
2270
2271 case DfpField.FLAG_UNDERFLOW:
2272 if ( (result.exp+mant.length) < MIN_EXP) {
2273 def = newInstance(getZero());
2274 def.sign = result.sign;
2275 } else {
2276 def = newInstance(result); // gradual underflow
2277 }
2278 result.exp = result.exp + ERR_SCALE;
2279 break;
2280
2281 case DfpField.FLAG_OVERFLOW:
2282 result.exp = result.exp - ERR_SCALE;
2283 def = newInstance(getZero());
2284 def.sign = result.sign;
2285 def.nans = INFINITE;
2286 break;
2287
2288 default: def = result; break;
2289 }
2290
2291 return trap(type, what, oper, def, result);
2292
2293 }
2294
2295 /** Trap handler. Subclasses may override this to provide trap
2296 * functionality per IEEE 854-1987.
2297 *
2298 * @param type The exception type - e.g. FLAG_OVERFLOW
2299 * @param what The name of the routine we were in e.g. divide()
2300 * @param oper An operand to this function if any
2301 * @param def The default return value if trap not enabled
2302 * @param result The result that is specified to be delivered per
2303 * IEEE 854, if any
2304 * @return the value that should be return by the operation triggering the trap
2305 */
2306 protected Dfp trap(int type, String what, Dfp oper, Dfp def, Dfp result) {
2307 return def;
2308 }
2309
2310 /** Returns the type - one of FINITE, INFINITE, SNAN, QNAN.
2311 * @return type of the number
2312 */
2313 public int classify() {
2314 return nans;
2315 }
2316
2317 /** Creates an instance that is the same as x except that it has the sign of y.
2318 * abs(x) = dfp.copysign(x, dfp.one)
2319 * @param x number to get the value from
2320 * @param y number to get the sign from
2321 * @return a number with the value of x and the sign of y
2322 */
2323 public static Dfp copysign(final Dfp x, final Dfp y) {
2324 Dfp result = x.newInstance(x);
2325 result.sign = y.sign;
2326 return result;
2327 }
2328
2329 /** Returns the next number greater than this one in the direction of x.
2330 * If this==x then simply returns this.
2331 * @param x direction where to look at
2332 * @return closest number next to instance in the direction of x
2333 */
2334 public Dfp nextAfter(final Dfp x) {
2335
2336 // make sure we don't mix number with different precision
2337 if (field.getRadixDigits() != x.field.getRadixDigits()) {
2338 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
2339 final Dfp result = newInstance(getZero());
2340 result.nans = QNAN;
2341 return dotrap(DfpField.FLAG_INVALID, NEXT_AFTER_TRAP, x, result);
2342 }
2343
2344 // if this is greater than x
2345 boolean up = false;
2346 if (this.lessThan(x)) {
2347 up = true;
2348 }
2349
2350 if (compare(this, x) == 0) {
2351 return newInstance(x);
2352 }
2353
2354 if (lessThan(getZero())) {
2355 up = !up;
2356 }
2357
2358 final Dfp inc;
2359 Dfp result;
2360 if (up) {
2361 inc = newInstance(getOne());
2362 inc.exp = this.exp-mant.length+1;
2363 inc.sign = this.sign;
2364
2365 if (this.equals(getZero())) {
2366 inc.exp = MIN_EXP-mant.length;
2367 }
2368
2369 result = add(inc);
2370 } else {
2371 inc = newInstance(getOne());
2372 inc.exp = this.exp;
2373 inc.sign = this.sign;
2374
2375 if (this.equals(inc)) {
2376 inc.exp = this.exp-mant.length;
2377 } else {
2378 inc.exp = this.exp-mant.length+1;
2379 }
2380
2381 if (this.equals(getZero())) {
2382 inc.exp = MIN_EXP-mant.length;
2383 }
2384
2385 result = this.subtract(inc);
2386 }
2387
2388 if (result.classify() == INFINITE && this.classify() != INFINITE) {
2389 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
2390 result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
2391 }
2392
2393 if (result.equals(getZero()) && this.equals(getZero()) == false) {
2394 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
2395 result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
2396 }
2397
2398 return result;
2399
2400 }
2401
2402 /** Convert the instance into a double.
2403 * @return a double approximating the instance
2404 * @see #toSplitDouble()
2405 */
2406 public double toDouble() {
2407
2408 if (isInfinite()) {
2409 if (lessThan(getZero())) {
2410 return Double.NEGATIVE_INFINITY;
2411 } else {
2412 return Double.POSITIVE_INFINITY;
2413 }
2414 }
2415
2416 if (isNaN()) {
2417 return Double.NaN;
2418 }
2419
2420 Dfp y = this;
2421 boolean negate = false;
2422 int cmp0 = compare(this, getZero());
2423 if (cmp0 == 0) {
2424 return sign < 0 ? -0.0 : +0.0;
2425 } else if (cmp0 < 0) {
2426 y = negate();
2427 negate = true;
2428 }
2429
2430 /* Find the exponent, first estimate by integer log10, then adjust.
2431 Should be faster than doing a natural logarithm. */
2432 int exponent = (int)(y.log10() * 3.32);
2433 if (exponent < 0) {
2434 exponent--;
2435 }
2436
2437 Dfp tempDfp = DfpMath.pow(getTwo(), exponent);
2438 while (tempDfp.lessThan(y) || tempDfp.equals(y)) {
2439 tempDfp = tempDfp.multiply(2);
2440 exponent++;
2441 }
2442 exponent--;
2443
2444 /* We have the exponent, now work on the mantissa */
2445
2446 y = y.divide(DfpMath.pow(getTwo(), exponent));
2447 if (exponent > -1023) {
2448 y = y.subtract(getOne());
2449 }
2450
2451 if (exponent < -1074) {
2452 return 0;
2453 }
2454
2455 if (exponent > 1023) {
2456 return negate ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
2457 }
2458
2459
2460 y = y.multiply(newInstance(4503599627370496l)).rint();
2461 String str = y.toString();
2462 str = str.substring(0, str.length()-1);
2463 long mantissa = Long.parseLong(str);
2464
2465 if (mantissa == 4503599627370496L) {
2466 // Handle special case where we round up to next power of two
2467 mantissa = 0;
2468 exponent++;
2469 }
2470
2471 /* Its going to be subnormal, so make adjustments */
2472 if (exponent <= -1023) {
2473 exponent--;
2474 }
2475
2476 while (exponent < -1023) {
2477 exponent++;
2478 mantissa >>>= 1;
2479 }
2480
2481 long bits = mantissa | ((exponent + 1023L) << 52);
2482 double x = Double.longBitsToDouble(bits);
2483
2484 if (negate) {
2485 x = -x;
2486 }
2487
2488 return x;
2489
2490 }
2491
2492 /** Convert the instance into a split double.
2493 * @return an array of two doubles which sum represent the instance
2494 * @see #toDouble()
2495 */
2496 public double[] toSplitDouble() {
2497 double split[] = new double[2];
2498 long mask = 0xffffffffc0000000L;
2499
2500 split[0] = Double.longBitsToDouble(Double.doubleToLongBits(toDouble()) & mask);
2501 split[1] = subtract(newInstance(split[0])).toDouble();
2502
2503 return split;
2504 }
2505
2506 }