source: src/main/java/agents/anac/y2019/harddealer/math3/dfp/Dfp.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

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