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