You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2015/02/16 23:40:06 UTC
[36/82] [partial] [math] Update for next development iteration:
commons-math4
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/dfp/Dfp.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/dfp/Dfp.java b/src/main/java/org/apache/commons/math3/dfp/Dfp.java
deleted file mode 100644
index 19d79a8..0000000
--- a/src/main/java/org/apache/commons/math3/dfp/Dfp.java
+++ /dev/null
@@ -1,2882 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-package org.apache.commons.math3.dfp;
-
-import java.util.Arrays;
-
-import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.exception.DimensionMismatchException;
-import org.apache.commons.math3.util.FastMath;
-
-/**
- * Decimal floating point library for Java
- *
- * <p>Another floating point class. This one is built using radix 10000
- * which is 10<sup>4</sup>, so its almost decimal.</p>
- *
- * <p>The design goals here are:
- * <ol>
- * <li>Decimal math, or close to it</li>
- * <li>Settable precision (but no mix between numbers using different settings)</li>
- * <li>Portability. Code should be kept as portable as possible.</li>
- * <li>Performance</li>
- * <li>Accuracy - Results should always be +/- 1 ULP for basic
- * algebraic operation</li>
- * <li>Comply with IEEE 854-1987 as much as possible.
- * (See IEEE 854-1987 notes below)</li>
- * </ol></p>
- *
- * <p>Trade offs:
- * <ol>
- * <li>Memory foot print. I'm using more memory than necessary to
- * represent numbers to get better performance.</li>
- * <li>Digits are bigger, so rounding is a greater loss. So, if you
- * really need 12 decimal digits, better use 4 base 10000 digits
- * there can be one partially filled.</li>
- * </ol></p>
- *
- * <p>Numbers are represented in the following form:
- * <pre>
- * n = sign × mant × (radix)<sup>exp</sup>;</p>
- * </pre>
- * where sign is ±1, mantissa represents a fractional number between
- * zero and one. mant[0] is the least significant digit.
- * exp is in the range of -32767 to 32768</p>
- *
- * <p>IEEE 854-1987 Notes and differences</p>
- *
- * <p>IEEE 854 requires the radix to be either 2 or 10. The radix here is
- * 10000, so that requirement is not met, but it is possible that a
- * subclassed can be made to make it behave as a radix 10
- * number. It is my opinion that if it looks and behaves as a radix
- * 10 number then it is one and that requirement would be met.</p>
- *
- * <p>The radix of 10000 was chosen because it should be faster to operate
- * on 4 decimal digits at once instead of one at a time. Radix 10 behavior
- * can be realized by adding an additional rounding step to ensure that
- * the number of decimal digits represented is constant.</p>
- *
- * <p>The IEEE standard specifically leaves out internal data encoding,
- * so it is reasonable to conclude that such a subclass of this radix
- * 10000 system is merely an encoding of a radix 10 system.</p>
- *
- * <p>IEEE 854 also specifies the existence of "sub-normal" numbers. This
- * class does not contain any such entities. The most significant radix
- * 10000 digit is always non-zero. Instead, we support "gradual underflow"
- * by raising the underflow flag for numbers less with exponent less than
- * expMin, but don't flush to zero until the exponent reaches MIN_EXP-digits.
- * Thus the smallest number we can represent would be:
- * 1E(-(MIN_EXP-digits-1)*4), eg, for digits=5, MIN_EXP=-32767, that would
- * be 1e-131092.</p>
- *
- * <p>IEEE 854 defines that the implied radix point lies just to the right
- * of the most significant digit and to the left of the remaining digits.
- * This implementation puts the implied radix point to the left of all
- * digits including the most significant one. The most significant digit
- * here is the one just to the right of the radix point. This is a fine
- * detail and is really only a matter of definition. Any side effects of
- * this can be rendered invisible by a subclass.</p>
- * @see DfpField
- * @since 2.2
- */
-public class Dfp implements RealFieldElement<Dfp> {
-
- /** The radix, or base of this system. Set to 10000 */
- public static final int RADIX = 10000;
-
- /** The minimum exponent before underflow is signaled. Flush to zero
- * occurs at minExp-DIGITS */
- public static final int MIN_EXP = -32767;
-
- /** The maximum exponent before overflow is signaled and results flushed
- * to infinity */
- public static final int MAX_EXP = 32768;
-
- /** The amount under/overflows are scaled by before going to trap handler */
- public static final int ERR_SCALE = 32760;
-
- /** Indicator value for normal finite numbers. */
- public static final byte FINITE = 0;
-
- /** Indicator value for Infinity. */
- public static final byte INFINITE = 1;
-
- /** Indicator value for signaling NaN. */
- public static final byte SNAN = 2;
-
- /** Indicator value for quiet NaN. */
- public static final byte QNAN = 3;
-
- /** String for NaN representation. */
- private static final String NAN_STRING = "NaN";
-
- /** String for positive infinity representation. */
- private static final String POS_INFINITY_STRING = "Infinity";
-
- /** String for negative infinity representation. */
- private static final String NEG_INFINITY_STRING = "-Infinity";
-
- /** Name for traps triggered by addition. */
- private static final String ADD_TRAP = "add";
-
- /** Name for traps triggered by multiplication. */
- private static final String MULTIPLY_TRAP = "multiply";
-
- /** Name for traps triggered by division. */
- private static final String DIVIDE_TRAP = "divide";
-
- /** Name for traps triggered by square root. */
- private static final String SQRT_TRAP = "sqrt";
-
- /** Name for traps triggered by alignment. */
- private static final String ALIGN_TRAP = "align";
-
- /** Name for traps triggered by truncation. */
- private static final String TRUNC_TRAP = "trunc";
-
- /** Name for traps triggered by nextAfter. */
- private static final String NEXT_AFTER_TRAP = "nextAfter";
-
- /** Name for traps triggered by lessThan. */
- private static final String LESS_THAN_TRAP = "lessThan";
-
- /** Name for traps triggered by greaterThan. */
- private static final String GREATER_THAN_TRAP = "greaterThan";
-
- /** Name for traps triggered by newInstance. */
- private static final String NEW_INSTANCE_TRAP = "newInstance";
-
- /** Mantissa. */
- protected int[] mant;
-
- /** Sign bit: 1 for positive, -1 for negative. */
- protected byte sign;
-
- /** Exponent. */
- protected int exp;
-
- /** Indicator for non-finite / non-number values. */
- protected byte nans;
-
- /** Factory building similar Dfp's. */
- private final DfpField field;
-
- /** Makes an instance with a value of zero.
- * @param field field to which this instance belongs
- */
- protected Dfp(final DfpField field) {
- mant = new int[field.getRadixDigits()];
- sign = 1;
- exp = 0;
- nans = FINITE;
- this.field = field;
- }
-
- /** Create an instance from a byte value.
- * @param field field to which this instance belongs
- * @param x value to convert to an instance
- */
- protected Dfp(final DfpField field, byte x) {
- this(field, (long) x);
- }
-
- /** Create an instance from an int value.
- * @param field field to which this instance belongs
- * @param x value to convert to an instance
- */
- protected Dfp(final DfpField field, int x) {
- this(field, (long) x);
- }
-
- /** Create an instance from a long value.
- * @param field field to which this instance belongs
- * @param x value to convert to an instance
- */
- protected Dfp(final DfpField field, long x) {
-
- // initialize as if 0
- mant = new int[field.getRadixDigits()];
- nans = FINITE;
- this.field = field;
-
- boolean isLongMin = false;
- if (x == Long.MIN_VALUE) {
- // special case for Long.MIN_VALUE (-9223372036854775808)
- // we must shift it before taking its absolute value
- isLongMin = true;
- ++x;
- }
-
- // set the sign
- if (x < 0) {
- sign = -1;
- x = -x;
- } else {
- sign = 1;
- }
-
- exp = 0;
- while (x != 0) {
- System.arraycopy(mant, mant.length - exp, mant, mant.length - 1 - exp, exp);
- mant[mant.length - 1] = (int) (x % RADIX);
- x /= RADIX;
- exp++;
- }
-
- if (isLongMin) {
- // remove the shift added for Long.MIN_VALUE
- // we know in this case that fixing the last digit is sufficient
- for (int i = 0; i < mant.length - 1; i++) {
- if (mant[i] != 0) {
- mant[i]++;
- break;
- }
- }
- }
- }
-
- /** Create an instance from a double value.
- * @param field field to which this instance belongs
- * @param x value to convert to an instance
- */
- protected Dfp(final DfpField field, double x) {
-
- // initialize as if 0
- mant = new int[field.getRadixDigits()];
- sign = 1;
- exp = 0;
- nans = FINITE;
- this.field = field;
-
- long bits = Double.doubleToLongBits(x);
- long mantissa = bits & 0x000fffffffffffffL;
- int exponent = (int) ((bits & 0x7ff0000000000000L) >> 52) - 1023;
-
- if (exponent == -1023) {
- // Zero or sub-normal
- if (x == 0) {
- // make sure 0 has the right sign
- if ((bits & 0x8000000000000000L) != 0) {
- sign = -1;
- }
- return;
- }
-
- exponent++;
-
- // Normalize the subnormal number
- while ( (mantissa & 0x0010000000000000L) == 0) {
- exponent--;
- mantissa <<= 1;
- }
- mantissa &= 0x000fffffffffffffL;
- }
-
- if (exponent == 1024) {
- // infinity or NAN
- if (x != x) {
- sign = (byte) 1;
- nans = QNAN;
- } else if (x < 0) {
- sign = (byte) -1;
- nans = INFINITE;
- } else {
- sign = (byte) 1;
- nans = INFINITE;
- }
- return;
- }
-
- Dfp xdfp = new Dfp(field, mantissa);
- xdfp = xdfp.divide(new Dfp(field, 4503599627370496l)).add(field.getOne()); // Divide by 2^52, then add one
- xdfp = xdfp.multiply(DfpMath.pow(field.getTwo(), exponent));
-
- if ((bits & 0x8000000000000000L) != 0) {
- xdfp = xdfp.negate();
- }
-
- System.arraycopy(xdfp.mant, 0, mant, 0, mant.length);
- sign = xdfp.sign;
- exp = xdfp.exp;
- nans = xdfp.nans;
-
- }
-
- /** Copy constructor.
- * @param d instance to copy
- */
- public Dfp(final Dfp d) {
- mant = d.mant.clone();
- sign = d.sign;
- exp = d.exp;
- nans = d.nans;
- field = d.field;
- }
-
- /** Create an instance from a String representation.
- * @param field field to which this instance belongs
- * @param s string representation of the instance
- */
- protected Dfp(final DfpField field, final String s) {
-
- // initialize as if 0
- mant = new int[field.getRadixDigits()];
- sign = 1;
- exp = 0;
- nans = FINITE;
- this.field = field;
-
- boolean decimalFound = false;
- final int rsize = 4; // size of radix in decimal digits
- final int offset = 4; // Starting offset into Striped
- final char[] striped = new char[getRadixDigits() * rsize + offset * 2];
-
- // Check some special cases
- if (s.equals(POS_INFINITY_STRING)) {
- sign = (byte) 1;
- nans = INFINITE;
- return;
- }
-
- if (s.equals(NEG_INFINITY_STRING)) {
- sign = (byte) -1;
- nans = INFINITE;
- return;
- }
-
- if (s.equals(NAN_STRING)) {
- sign = (byte) 1;
- nans = QNAN;
- return;
- }
-
- // Check for scientific notation
- int p = s.indexOf("e");
- if (p == -1) { // try upper case?
- p = s.indexOf("E");
- }
-
- final String fpdecimal;
- int sciexp = 0;
- if (p != -1) {
- // scientific notation
- fpdecimal = s.substring(0, p);
- String fpexp = s.substring(p+1);
- boolean negative = false;
-
- for (int i=0; i<fpexp.length(); i++)
- {
- if (fpexp.charAt(i) == '-')
- {
- negative = true;
- continue;
- }
- if (fpexp.charAt(i) >= '0' && fpexp.charAt(i) <= '9') {
- sciexp = sciexp * 10 + fpexp.charAt(i) - '0';
- }
- }
-
- if (negative) {
- sciexp = -sciexp;
- }
- } else {
- // normal case
- fpdecimal = s;
- }
-
- // If there is a minus sign in the number then it is negative
- if (fpdecimal.indexOf("-") != -1) {
- sign = -1;
- }
-
- // First off, find all of the leading zeros, trailing zeros, and significant digits
- p = 0;
-
- // Move p to first significant digit
- int decimalPos = 0;
- for (;;) {
- if (fpdecimal.charAt(p) >= '1' && fpdecimal.charAt(p) <= '9') {
- break;
- }
-
- if (decimalFound && fpdecimal.charAt(p) == '0') {
- decimalPos--;
- }
-
- if (fpdecimal.charAt(p) == '.') {
- decimalFound = true;
- }
-
- p++;
-
- if (p == fpdecimal.length()) {
- break;
- }
- }
-
- // Copy the string onto Stripped
- int q = offset;
- striped[0] = '0';
- striped[1] = '0';
- striped[2] = '0';
- striped[3] = '0';
- int significantDigits=0;
- for(;;) {
- if (p == (fpdecimal.length())) {
- break;
- }
-
- // Don't want to run pass the end of the array
- if (q == mant.length*rsize+offset+1) {
- break;
- }
-
- if (fpdecimal.charAt(p) == '.') {
- decimalFound = true;
- decimalPos = significantDigits;
- p++;
- continue;
- }
-
- if (fpdecimal.charAt(p) < '0' || fpdecimal.charAt(p) > '9') {
- p++;
- continue;
- }
-
- striped[q] = fpdecimal.charAt(p);
- q++;
- p++;
- significantDigits++;
- }
-
-
- // If the decimal point has been found then get rid of trailing zeros.
- if (decimalFound && q != offset) {
- for (;;) {
- q--;
- if (q == offset) {
- break;
- }
- if (striped[q] == '0') {
- significantDigits--;
- } else {
- break;
- }
- }
- }
-
- // special case of numbers like "0.00000"
- if (decimalFound && significantDigits == 0) {
- decimalPos = 0;
- }
-
- // Implicit decimal point at end of number if not present
- if (!decimalFound) {
- decimalPos = q-offset;
- }
-
- // Find the number of significant trailing zeros
- q = offset; // set q to point to first sig digit
- p = significantDigits-1+offset;
-
- while (p > q) {
- if (striped[p] != '0') {
- break;
- }
- p--;
- }
-
- // Make sure the decimal is on a mod 10000 boundary
- int i = ((rsize * 100) - decimalPos - sciexp % rsize) % rsize;
- q -= i;
- decimalPos += i;
-
- // Make the mantissa length right by adding zeros at the end if necessary
- while ((p - q) < (mant.length * rsize)) {
- for (i = 0; i < rsize; i++) {
- striped[++p] = '0';
- }
- }
-
- // Ok, now we know how many trailing zeros there are,
- // and where the least significant digit is
- for (i = mant.length - 1; i >= 0; i--) {
- mant[i] = (striped[q] - '0') * 1000 +
- (striped[q+1] - '0') * 100 +
- (striped[q+2] - '0') * 10 +
- (striped[q+3] - '0');
- q += 4;
- }
-
-
- exp = (decimalPos+sciexp) / rsize;
-
- if (q < striped.length) {
- // Is there possible another digit?
- round((striped[q] - '0')*1000);
- }
-
- }
-
- /** Creates an instance with a non-finite value.
- * @param field field to which this instance belongs
- * @param sign sign of the Dfp to create
- * @param nans code of the value, must be one of {@link #INFINITE},
- * {@link #SNAN}, {@link #QNAN}
- */
- protected Dfp(final DfpField field, final byte sign, final byte nans) {
- this.field = field;
- this.mant = new int[field.getRadixDigits()];
- this.sign = sign;
- this.exp = 0;
- this.nans = nans;
- }
-
- /** Create an instance with a value of 0.
- * Use this internally in preference to constructors to facilitate subclasses
- * @return a new instance with a value of 0
- */
- public Dfp newInstance() {
- return new Dfp(getField());
- }
-
- /** Create an instance from a byte value.
- * @param x value to convert to an instance
- * @return a new instance with value x
- */
- public Dfp newInstance(final byte x) {
- return new Dfp(getField(), x);
- }
-
- /** Create an instance from an int value.
- * @param x value to convert to an instance
- * @return a new instance with value x
- */
- public Dfp newInstance(final int x) {
- return new Dfp(getField(), x);
- }
-
- /** Create an instance from a long value.
- * @param x value to convert to an instance
- * @return a new instance with value x
- */
- public Dfp newInstance(final long x) {
- return new Dfp(getField(), x);
- }
-
- /** Create an instance from a double value.
- * @param x value to convert to an instance
- * @return a new instance with value x
- */
- public Dfp newInstance(final double x) {
- return new Dfp(getField(), x);
- }
-
- /** Create an instance by copying an existing one.
- * Use this internally in preference to constructors to facilitate subclasses.
- * @param d instance to copy
- * @return a new instance with the same value as d
- */
- public Dfp newInstance(final Dfp d) {
-
- // make sure we don't mix number with different precision
- if (field.getRadixDigits() != d.field.getRadixDigits()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = newInstance(getZero());
- result.nans = QNAN;
- return dotrap(DfpField.FLAG_INVALID, NEW_INSTANCE_TRAP, d, result);
- }
-
- return new Dfp(d);
-
- }
-
- /** Create an instance from a String representation.
- * Use this internally in preference to constructors to facilitate subclasses.
- * @param s string representation of the instance
- * @return a new instance parsed from specified string
- */
- public Dfp newInstance(final String s) {
- return new Dfp(field, s);
- }
-
- /** Creates an instance with a non-finite value.
- * @param sig sign of the Dfp to create
- * @param code code of the value, must be one of {@link #INFINITE},
- * {@link #SNAN}, {@link #QNAN}
- * @return a new instance with a non-finite value
- */
- public Dfp newInstance(final byte sig, final byte code) {
- return field.newDfp(sig, code);
- }
-
- /** Get the {@link org.apache.commons.math3.Field Field} (really a {@link DfpField}) to which the instance belongs.
- * <p>
- * The field is linked to the number of digits and acts as a factory
- * for {@link Dfp} instances.
- * </p>
- * @return {@link org.apache.commons.math3.Field Field} (really a {@link DfpField}) to which the instance belongs
- */
- public DfpField getField() {
- return field;
- }
-
- /** Get the number of radix digits of the instance.
- * @return number of radix digits
- */
- public int getRadixDigits() {
- return field.getRadixDigits();
- }
-
- /** Get the constant 0.
- * @return a Dfp with value zero
- */
- public Dfp getZero() {
- return field.getZero();
- }
-
- /** Get the constant 1.
- * @return a Dfp with value one
- */
- public Dfp getOne() {
- return field.getOne();
- }
-
- /** Get the constant 2.
- * @return a Dfp with value two
- */
- public Dfp getTwo() {
- return field.getTwo();
- }
-
- /** Shift the mantissa left, and adjust the exponent to compensate.
- */
- protected void shiftLeft() {
- for (int i = mant.length - 1; i > 0; i--) {
- mant[i] = mant[i-1];
- }
- mant[0] = 0;
- exp--;
- }
-
- /* Note that shiftRight() does not call round() as that round() itself
- uses shiftRight() */
- /** Shift the mantissa right, and adjust the exponent to compensate.
- */
- protected void shiftRight() {
- for (int i = 0; i < mant.length - 1; i++) {
- mant[i] = mant[i+1];
- }
- mant[mant.length - 1] = 0;
- exp++;
- }
-
- /** Make our exp equal to the supplied one, this may cause rounding.
- * Also causes de-normalized numbers. These numbers are generally
- * dangerous because most routines assume normalized numbers.
- * Align doesn't round, so it will return the last digit destroyed
- * by shifting right.
- * @param e desired exponent
- * @return last digit destroyed by shifting right
- */
- protected int align(int e) {
- int lostdigit = 0;
- boolean inexact = false;
-
- int diff = exp - e;
-
- int adiff = diff;
- if (adiff < 0) {
- adiff = -adiff;
- }
-
- if (diff == 0) {
- return 0;
- }
-
- if (adiff > (mant.length + 1)) {
- // Special case
- Arrays.fill(mant, 0);
- exp = e;
-
- field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
- dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
-
- return 0;
- }
-
- for (int i = 0; i < adiff; i++) {
- if (diff < 0) {
- /* Keep track of loss -- only signal inexact after losing 2 digits.
- * the first lost digit is returned to add() and may be incorporated
- * into the result.
- */
- if (lostdigit != 0) {
- inexact = true;
- }
-
- lostdigit = mant[0];
-
- shiftRight();
- } else {
- shiftLeft();
- }
- }
-
- if (inexact) {
- field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
- dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
- }
-
- return lostdigit;
-
- }
-
- /** Check if instance is less than x.
- * @param x number to check instance against
- * @return true if instance is less than x and neither are NaN, false otherwise
- */
- public boolean lessThan(final Dfp x) {
-
- // make sure we don't mix number with different precision
- if (field.getRadixDigits() != x.field.getRadixDigits()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = newInstance(getZero());
- result.nans = QNAN;
- dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, result);
- return false;
- }
-
- /* if a nan is involved, signal invalid and return false */
- if (isNaN() || x.isNaN()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, newInstance(getZero()));
- return false;
- }
-
- return compare(this, x) < 0;
- }
-
- /** Check if instance is greater than x.
- * @param x number to check instance against
- * @return true if instance is greater than x and neither are NaN, false otherwise
- */
- public boolean greaterThan(final Dfp x) {
-
- // make sure we don't mix number with different precision
- if (field.getRadixDigits() != x.field.getRadixDigits()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = newInstance(getZero());
- result.nans = QNAN;
- dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, result);
- return false;
- }
-
- /* if a nan is involved, signal invalid and return false */
- if (isNaN() || x.isNaN()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, newInstance(getZero()));
- return false;
- }
-
- return compare(this, x) > 0;
- }
-
- /** Check if instance is less than or equal to 0.
- * @return true if instance is not NaN and less than or equal to 0, false otherwise
- */
- public boolean negativeOrNull() {
-
- if (isNaN()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
- return false;
- }
-
- return (sign < 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
-
- }
-
- /** Check if instance is strictly less than 0.
- * @return true if instance is not NaN and less than or equal to 0, false otherwise
- */
- public boolean strictlyNegative() {
-
- if (isNaN()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
- return false;
- }
-
- return (sign < 0) && ((mant[mant.length - 1] != 0) || isInfinite());
-
- }
-
- /** Check if instance is greater than or equal to 0.
- * @return true if instance is not NaN and greater than or equal to 0, false otherwise
- */
- public boolean positiveOrNull() {
-
- if (isNaN()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
- return false;
- }
-
- return (sign > 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
-
- }
-
- /** Check if instance is strictly greater than 0.
- * @return true if instance is not NaN and greater than or equal to 0, false otherwise
- */
- public boolean strictlyPositive() {
-
- if (isNaN()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
- return false;
- }
-
- return (sign > 0) && ((mant[mant.length - 1] != 0) || isInfinite());
-
- }
-
- /** Get the absolute value of instance.
- * @return absolute value of instance
- * @since 3.2
- */
- public Dfp abs() {
- Dfp result = newInstance(this);
- result.sign = 1;
- return result;
- }
-
- /** Check if instance is infinite.
- * @return true if instance is infinite
- */
- public boolean isInfinite() {
- return nans == INFINITE;
- }
-
- /** Check if instance is not a number.
- * @return true if instance is not a number
- */
- public boolean isNaN() {
- return (nans == QNAN) || (nans == SNAN);
- }
-
- /** Check if instance is equal to zero.
- * @return true if instance is equal to zero
- */
- public boolean isZero() {
-
- if (isNaN()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
- return false;
- }
-
- return (mant[mant.length - 1] == 0) && !isInfinite();
-
- }
-
- /** Check if instance is equal to x.
- * @param other object to check instance against
- * @return true if instance is equal to x and neither are NaN, false otherwise
- */
- @Override
- public boolean equals(final Object other) {
-
- if (other instanceof Dfp) {
- final Dfp x = (Dfp) other;
- if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
- return false;
- }
-
- return compare(this, x) == 0;
- }
-
- return false;
-
- }
-
- /**
- * Gets a hashCode for the instance.
- * @return a hash code value for this object
- */
- @Override
- public int hashCode() {
- return 17 + (sign << 8) + (nans << 16) + exp + Arrays.hashCode(mant);
- }
-
- /** Check if instance is not equal to x.
- * @param x number to check instance against
- * @return true if instance is not equal to x and neither are NaN, false otherwise
- */
- public boolean unequal(final Dfp x) {
- if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
- return false;
- }
-
- return greaterThan(x) || lessThan(x);
- }
-
- /** Compare two instances.
- * @param a first instance in comparison
- * @param b second instance in comparison
- * @return -1 if a<b, 1 if a>b and 0 if a==b
- * Note this method does not properly handle NaNs or numbers with different precision.
- */
- private static int compare(final Dfp a, final Dfp b) {
- // Ignore the sign of zero
- if (a.mant[a.mant.length - 1] == 0 && b.mant[b.mant.length - 1] == 0 &&
- a.nans == FINITE && b.nans == FINITE) {
- return 0;
- }
-
- if (a.sign != b.sign) {
- if (a.sign == -1) {
- return -1;
- } else {
- return 1;
- }
- }
-
- // deal with the infinities
- if (a.nans == INFINITE && b.nans == FINITE) {
- return a.sign;
- }
-
- if (a.nans == FINITE && b.nans == INFINITE) {
- return -b.sign;
- }
-
- if (a.nans == INFINITE && b.nans == INFINITE) {
- return 0;
- }
-
- // Handle special case when a or b is zero, by ignoring the exponents
- if (b.mant[b.mant.length-1] != 0 && a.mant[b.mant.length-1] != 0) {
- if (a.exp < b.exp) {
- return -a.sign;
- }
-
- if (a.exp > b.exp) {
- return a.sign;
- }
- }
-
- // compare the mantissas
- for (int i = a.mant.length - 1; i >= 0; i--) {
- if (a.mant[i] > b.mant[i]) {
- return a.sign;
- }
-
- if (a.mant[i] < b.mant[i]) {
- return -a.sign;
- }
- }
-
- return 0;
-
- }
-
- /** Round to nearest integer using the round-half-even method.
- * That is round to nearest integer unless both are equidistant.
- * In which case round to the even one.
- * @return rounded value
- * @since 3.2
- */
- public Dfp rint() {
- return trunc(DfpField.RoundingMode.ROUND_HALF_EVEN);
- }
-
- /** Round to an integer using the round floor mode.
- * That is, round toward -Infinity
- * @return rounded value
- * @since 3.2
- */
- public Dfp floor() {
- return trunc(DfpField.RoundingMode.ROUND_FLOOR);
- }
-
- /** Round to an integer using the round ceil mode.
- * That is, round toward +Infinity
- * @return rounded value
- * @since 3.2
- */
- public Dfp ceil() {
- return trunc(DfpField.RoundingMode.ROUND_CEIL);
- }
-
- /** Returns the IEEE remainder.
- * @param d divisor
- * @return this less n × d, where n is the integer closest to this/d
- * @since 3.2
- */
- public Dfp remainder(final Dfp d) {
-
- final Dfp result = this.subtract(this.divide(d).rint().multiply(d));
-
- // IEEE 854-1987 says that if the result is zero, then it carries the sign of this
- if (result.mant[mant.length-1] == 0) {
- result.sign = sign;
- }
-
- return result;
-
- }
-
- /** Does the integer conversions with the specified rounding.
- * @param rmode rounding mode to use
- * @return truncated value
- */
- protected Dfp trunc(final DfpField.RoundingMode rmode) {
- boolean changed = false;
-
- if (isNaN()) {
- return newInstance(this);
- }
-
- if (nans == INFINITE) {
- return newInstance(this);
- }
-
- if (mant[mant.length-1] == 0) {
- // a is zero
- return newInstance(this);
- }
-
- /* If the exponent is less than zero then we can certainly
- * return zero */
- if (exp < 0) {
- field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
- Dfp result = newInstance(getZero());
- result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
- return result;
- }
-
- /* If the exponent is greater than or equal to digits, then it
- * must already be an integer since there is no precision left
- * for any fractional part */
-
- if (exp >= mant.length) {
- return newInstance(this);
- }
-
- /* General case: create another dfp, result, that contains the
- * a with the fractional part lopped off. */
-
- Dfp result = newInstance(this);
- for (int i = 0; i < mant.length-result.exp; i++) {
- changed |= result.mant[i] != 0;
- result.mant[i] = 0;
- }
-
- if (changed) {
- switch (rmode) {
- case ROUND_FLOOR:
- if (result.sign == -1) {
- // then we must increment the mantissa by one
- result = result.add(newInstance(-1));
- }
- break;
-
- case ROUND_CEIL:
- if (result.sign == 1) {
- // then we must increment the mantissa by one
- result = result.add(getOne());
- }
- break;
-
- case ROUND_HALF_EVEN:
- default:
- final Dfp half = newInstance("0.5");
- Dfp a = subtract(result); // difference between this and result
- a.sign = 1; // force positive (take abs)
- if (a.greaterThan(half)) {
- a = newInstance(getOne());
- a.sign = sign;
- result = result.add(a);
- }
-
- /** If exactly equal to 1/2 and odd then increment */
- if (a.equals(half) && result.exp > 0 && (result.mant[mant.length-result.exp]&1) != 0) {
- a = newInstance(getOne());
- a.sign = sign;
- result = result.add(a);
- }
- break;
- }
-
- field.setIEEEFlagsBits(DfpField.FLAG_INEXACT); // signal inexact
- result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
- return result;
- }
-
- return result;
- }
-
- /** Convert this to an integer.
- * If greater than 2147483647, it returns 2147483647. If less than -2147483648 it returns -2147483648.
- * @return converted number
- */
- public int intValue() {
- Dfp rounded;
- int result = 0;
-
- rounded = rint();
-
- if (rounded.greaterThan(newInstance(2147483647))) {
- return 2147483647;
- }
-
- if (rounded.lessThan(newInstance(-2147483648))) {
- return -2147483648;
- }
-
- for (int i = mant.length - 1; i >= mant.length - rounded.exp; i--) {
- result = result * RADIX + rounded.mant[i];
- }
-
- if (rounded.sign == -1) {
- result = -result;
- }
-
- return result;
- }
-
- /** Get the exponent of the greatest power of 10000 that is
- * less than or equal to the absolute value of this. I.E. if
- * this is 10<sup>6</sup> then log10K would return 1.
- * @return integer base 10000 logarithm
- */
- public int log10K() {
- return exp - 1;
- }
-
- /** Get the specified power of 10000.
- * @param e desired power
- * @return 10000<sup>e</sup>
- */
- public Dfp power10K(final int e) {
- Dfp d = newInstance(getOne());
- d.exp = e + 1;
- return d;
- }
-
- /** Get the exponent of the greatest power of 10 that is less than or equal to abs(this).
- * @return integer base 10 logarithm
- * @since 3.2
- */
- public int intLog10() {
- if (mant[mant.length-1] > 1000) {
- return exp * 4 - 1;
- }
- if (mant[mant.length-1] > 100) {
- return exp * 4 - 2;
- }
- if (mant[mant.length-1] > 10) {
- return exp * 4 - 3;
- }
- return exp * 4 - 4;
- }
-
- /** Return the specified power of 10.
- * @param e desired power
- * @return 10<sup>e</sup>
- */
- public Dfp power10(final int e) {
- Dfp d = newInstance(getOne());
-
- if (e >= 0) {
- d.exp = e / 4 + 1;
- } else {
- d.exp = (e + 1) / 4;
- }
-
- switch ((e % 4 + 4) % 4) {
- case 0:
- break;
- case 1:
- d = d.multiply(10);
- break;
- case 2:
- d = d.multiply(100);
- break;
- default:
- d = d.multiply(1000);
- }
-
- return d;
- }
-
- /** Negate the mantissa of this by computing the complement.
- * Leaves the sign bit unchanged, used internally by add.
- * Denormalized numbers are handled properly here.
- * @param extra ???
- * @return ???
- */
- protected int complement(int extra) {
-
- extra = RADIX-extra;
- for (int i = 0; i < mant.length; i++) {
- mant[i] = RADIX-mant[i]-1;
- }
-
- int rh = extra / RADIX;
- extra -= rh * RADIX;
- for (int i = 0; i < mant.length; i++) {
- final int r = mant[i] + rh;
- rh = r / RADIX;
- mant[i] = r - rh * RADIX;
- }
-
- return extra;
- }
-
- /** Add x to this.
- * @param x number to add
- * @return sum of this and x
- */
- public Dfp add(final Dfp x) {
-
- // make sure we don't mix number with different precision
- if (field.getRadixDigits() != x.field.getRadixDigits()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = newInstance(getZero());
- result.nans = QNAN;
- return dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
- }
-
- /* handle special cases */
- if (nans != FINITE || x.nans != FINITE) {
- if (isNaN()) {
- return this;
- }
-
- if (x.isNaN()) {
- return x;
- }
-
- if (nans == INFINITE && x.nans == FINITE) {
- return this;
- }
-
- if (x.nans == INFINITE && nans == FINITE) {
- return x;
- }
-
- if (x.nans == INFINITE && nans == INFINITE && sign == x.sign) {
- return x;
- }
-
- if (x.nans == INFINITE && nans == INFINITE && sign != x.sign) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- Dfp result = newInstance(getZero());
- result.nans = QNAN;
- result = dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
- return result;
- }
- }
-
- /* copy this and the arg */
- Dfp a = newInstance(this);
- Dfp b = newInstance(x);
-
- /* initialize the result object */
- Dfp result = newInstance(getZero());
-
- /* Make all numbers positive, but remember their sign */
- final byte asign = a.sign;
- final byte bsign = b.sign;
-
- a.sign = 1;
- b.sign = 1;
-
- /* The result will be signed like the arg with greatest magnitude */
- byte rsign = bsign;
- if (compare(a, b) > 0) {
- rsign = asign;
- }
-
- /* Handle special case when a or b is zero, by setting the exponent
- of the zero number equal to the other one. This avoids an alignment
- which would cause catastropic loss of precision */
- if (b.mant[mant.length-1] == 0) {
- b.exp = a.exp;
- }
-
- if (a.mant[mant.length-1] == 0) {
- a.exp = b.exp;
- }
-
- /* align number with the smaller exponent */
- int aextradigit = 0;
- int bextradigit = 0;
- if (a.exp < b.exp) {
- aextradigit = a.align(b.exp);
- } else {
- bextradigit = b.align(a.exp);
- }
-
- /* complement the smaller of the two if the signs are different */
- if (asign != bsign) {
- if (asign == rsign) {
- bextradigit = b.complement(bextradigit);
- } else {
- aextradigit = a.complement(aextradigit);
- }
- }
-
- /* add the mantissas */
- int rh = 0; /* acts as a carry */
- for (int i = 0; i < mant.length; i++) {
- final int r = a.mant[i]+b.mant[i]+rh;
- rh = r / RADIX;
- result.mant[i] = r - rh * RADIX;
- }
- result.exp = a.exp;
- result.sign = rsign;
-
- /* handle overflow -- note, when asign!=bsign an overflow is
- * normal and should be ignored. */
-
- if (rh != 0 && (asign == bsign)) {
- final int lostdigit = result.mant[0];
- result.shiftRight();
- result.mant[mant.length-1] = rh;
- final int excp = result.round(lostdigit);
- if (excp != 0) {
- result = dotrap(excp, ADD_TRAP, x, result);
- }
- }
-
- /* normalize the result */
- for (int i = 0; i < mant.length; i++) {
- if (result.mant[mant.length-1] != 0) {
- break;
- }
- result.shiftLeft();
- if (i == 0) {
- result.mant[0] = aextradigit+bextradigit;
- aextradigit = 0;
- bextradigit = 0;
- }
- }
-
- /* result is zero if after normalization the most sig. digit is zero */
- if (result.mant[mant.length-1] == 0) {
- result.exp = 0;
-
- if (asign != bsign) {
- // Unless adding 2 negative zeros, sign is positive
- result.sign = 1; // Per IEEE 854-1987 Section 6.3
- }
- }
-
- /* Call round to test for over/under flows */
- final int excp = result.round(aextradigit + bextradigit);
- if (excp != 0) {
- result = dotrap(excp, ADD_TRAP, x, result);
- }
-
- return result;
- }
-
- /** Returns a number that is this number with the sign bit reversed.
- * @return the opposite of this
- */
- public Dfp negate() {
- Dfp result = newInstance(this);
- result.sign = (byte) - result.sign;
- return result;
- }
-
- /** Subtract x from this.
- * @param x number to subtract
- * @return difference of this and a
- */
- public Dfp subtract(final Dfp x) {
- return add(x.negate());
- }
-
- /** Round this given the next digit n using the current rounding mode.
- * @param n ???
- * @return the IEEE flag if an exception occurred
- */
- protected int round(int n) {
- boolean inc = false;
- switch (field.getRoundingMode()) {
- case ROUND_DOWN:
- inc = false;
- break;
-
- case ROUND_UP:
- inc = n != 0; // round up if n!=0
- break;
-
- case ROUND_HALF_UP:
- inc = n >= 5000; // round half up
- break;
-
- case ROUND_HALF_DOWN:
- inc = n > 5000; // round half down
- break;
-
- case ROUND_HALF_EVEN:
- inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 1); // round half-even
- break;
-
- case ROUND_HALF_ODD:
- inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 0); // round half-odd
- break;
-
- case ROUND_CEIL:
- inc = sign == 1 && n != 0; // round ceil
- break;
-
- case ROUND_FLOOR:
- default:
- inc = sign == -1 && n != 0; // round floor
- break;
- }
-
- if (inc) {
- // increment if necessary
- int rh = 1;
- for (int i = 0; i < mant.length; i++) {
- final int r = mant[i] + rh;
- rh = r / RADIX;
- mant[i] = r - rh * RADIX;
- }
-
- if (rh != 0) {
- shiftRight();
- mant[mant.length-1] = rh;
- }
- }
-
- // check for exceptional cases and raise signals if necessary
- if (exp < MIN_EXP) {
- // Gradual Underflow
- field.setIEEEFlagsBits(DfpField.FLAG_UNDERFLOW);
- return DfpField.FLAG_UNDERFLOW;
- }
-
- if (exp > MAX_EXP) {
- // Overflow
- field.setIEEEFlagsBits(DfpField.FLAG_OVERFLOW);
- return DfpField.FLAG_OVERFLOW;
- }
-
- if (n != 0) {
- // Inexact
- field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
- return DfpField.FLAG_INEXACT;
- }
-
- return 0;
-
- }
-
- /** Multiply this by x.
- * @param x multiplicand
- * @return product of this and x
- */
- public Dfp multiply(final Dfp x) {
-
- // make sure we don't mix number with different precision
- if (field.getRadixDigits() != x.field.getRadixDigits()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = newInstance(getZero());
- result.nans = QNAN;
- return dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
- }
-
- Dfp result = newInstance(getZero());
-
- /* handle special cases */
- if (nans != FINITE || x.nans != FINITE) {
- if (isNaN()) {
- return this;
- }
-
- if (x.isNaN()) {
- return x;
- }
-
- if (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] != 0) {
- result = newInstance(this);
- result.sign = (byte) (sign * x.sign);
- return result;
- }
-
- if (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] != 0) {
- result = newInstance(x);
- result.sign = (byte) (sign * x.sign);
- return result;
- }
-
- if (x.nans == INFINITE && nans == INFINITE) {
- result = newInstance(this);
- result.sign = (byte) (sign * x.sign);
- return result;
- }
-
- if ( (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] == 0) ||
- (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] == 0) ) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- result = newInstance(getZero());
- result.nans = QNAN;
- result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
- return result;
- }
- }
-
- int[] product = new int[mant.length*2]; // Big enough to hold even the largest result
-
- for (int i = 0; i < mant.length; i++) {
- int rh = 0; // acts as a carry
- for (int j=0; j<mant.length; j++) {
- int r = mant[i] * x.mant[j]; // multiply the 2 digits
- r += product[i+j] + rh; // add to the product digit with carry in
-
- rh = r / RADIX;
- product[i+j] = r - rh * RADIX;
- }
- product[i+mant.length] = rh;
- }
-
- // Find the most sig digit
- int md = mant.length * 2 - 1; // default, in case result is zero
- for (int i = mant.length * 2 - 1; i >= 0; i--) {
- if (product[i] != 0) {
- md = i;
- break;
- }
- }
-
- // Copy the digits into the result
- for (int i = 0; i < mant.length; i++) {
- result.mant[mant.length - i - 1] = product[md - i];
- }
-
- // Fixup the exponent.
- result.exp = exp + x.exp + md - 2 * mant.length + 1;
- result.sign = (byte)((sign == x.sign)?1:-1);
-
- if (result.mant[mant.length-1] == 0) {
- // if result is zero, set exp to zero
- result.exp = 0;
- }
-
- final int excp;
- if (md > (mant.length-1)) {
- excp = result.round(product[md-mant.length]);
- } else {
- excp = result.round(0); // has no effect except to check status
- }
-
- if (excp != 0) {
- result = dotrap(excp, MULTIPLY_TRAP, x, result);
- }
-
- return result;
-
- }
-
- /** Multiply this by a single digit x.
- * @param x multiplicand
- * @return product of this and x
- */
- public Dfp multiply(final int x) {
- if (x >= 0 && x < RADIX) {
- return multiplyFast(x);
- } else {
- return multiply(newInstance(x));
- }
- }
-
- /** Multiply this by a single digit 0<=x<radix.
- * There are speed advantages in this special case.
- * @param x multiplicand
- * @return product of this and x
- */
- private Dfp multiplyFast(final int x) {
- Dfp result = newInstance(this);
-
- /* handle special cases */
- if (nans != FINITE) {
- if (isNaN()) {
- return this;
- }
-
- if (nans == INFINITE && x != 0) {
- result = newInstance(this);
- return result;
- }
-
- if (nans == INFINITE && x == 0) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- result = newInstance(getZero());
- result.nans = QNAN;
- result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, newInstance(getZero()), result);
- return result;
- }
- }
-
- /* range check x */
- if (x < 0 || x >= RADIX) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- result = newInstance(getZero());
- result.nans = QNAN;
- result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, result, result);
- return result;
- }
-
- int rh = 0;
- for (int i = 0; i < mant.length; i++) {
- final int r = mant[i] * x + rh;
- rh = r / RADIX;
- result.mant[i] = r - rh * RADIX;
- }
-
- int lostdigit = 0;
- if (rh != 0) {
- lostdigit = result.mant[0];
- result.shiftRight();
- result.mant[mant.length-1] = rh;
- }
-
- if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
- result.exp = 0;
- }
-
- final int excp = result.round(lostdigit);
- if (excp != 0) {
- result = dotrap(excp, MULTIPLY_TRAP, result, result);
- }
-
- return result;
- }
-
- /** Divide this by divisor.
- * @param divisor divisor
- * @return quotient of this by divisor
- */
- public Dfp divide(Dfp divisor) {
- int dividend[]; // current status of the dividend
- int quotient[]; // quotient
- int remainder[];// remainder
- int qd; // current quotient digit we're working with
- int nsqd; // number of significant quotient digits we have
- int trial=0; // trial quotient digit
- int minadj; // minimum adjustment
- boolean trialgood; // Flag to indicate a good trail digit
- int md=0; // most sig digit in result
- int excp; // exceptions
-
- // make sure we don't mix number with different precision
- if (field.getRadixDigits() != divisor.field.getRadixDigits()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = newInstance(getZero());
- result.nans = QNAN;
- return dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
- }
-
- Dfp result = newInstance(getZero());
-
- /* handle special cases */
- if (nans != FINITE || divisor.nans != FINITE) {
- if (isNaN()) {
- return this;
- }
-
- if (divisor.isNaN()) {
- return divisor;
- }
-
- if (nans == INFINITE && divisor.nans == FINITE) {
- result = newInstance(this);
- result.sign = (byte) (sign * divisor.sign);
- return result;
- }
-
- if (divisor.nans == INFINITE && nans == FINITE) {
- result = newInstance(getZero());
- result.sign = (byte) (sign * divisor.sign);
- return result;
- }
-
- if (divisor.nans == INFINITE && nans == INFINITE) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- result = newInstance(getZero());
- result.nans = QNAN;
- result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
- return result;
- }
- }
-
- /* Test for divide by zero */
- if (divisor.mant[mant.length-1] == 0) {
- field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
- result = newInstance(getZero());
- result.sign = (byte) (sign * divisor.sign);
- result.nans = INFINITE;
- result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, divisor, result);
- return result;
- }
-
- dividend = new int[mant.length+1]; // one extra digit needed
- quotient = new int[mant.length+2]; // two extra digits needed 1 for overflow, 1 for rounding
- remainder = new int[mant.length+1]; // one extra digit needed
-
- /* Initialize our most significant digits to zero */
-
- dividend[mant.length] = 0;
- quotient[mant.length] = 0;
- quotient[mant.length+1] = 0;
- remainder[mant.length] = 0;
-
- /* copy our mantissa into the dividend, initialize the
- quotient while we are at it */
-
- for (int i = 0; i < mant.length; i++) {
- dividend[i] = mant[i];
- quotient[i] = 0;
- remainder[i] = 0;
- }
-
- /* outer loop. Once per quotient digit */
- nsqd = 0;
- for (qd = mant.length+1; qd >= 0; qd--) {
- /* Determine outer limits of our quotient digit */
-
- // r = most sig 2 digits of dividend
- final int divMsb = dividend[mant.length]*RADIX+dividend[mant.length-1];
- int min = divMsb / (divisor.mant[mant.length-1]+1);
- int max = (divMsb + 1) / divisor.mant[mant.length-1];
-
- trialgood = false;
- while (!trialgood) {
- // try the mean
- trial = (min+max)/2;
-
- /* Multiply by divisor and store as remainder */
- int rh = 0;
- for (int i = 0; i < mant.length + 1; i++) {
- int dm = (i<mant.length)?divisor.mant[i]:0;
- final int r = (dm * trial) + rh;
- rh = r / RADIX;
- remainder[i] = r - rh * RADIX;
- }
-
- /* subtract the remainder from the dividend */
- rh = 1; // carry in to aid the subtraction
- for (int i = 0; i < mant.length + 1; i++) {
- final int r = ((RADIX-1) - remainder[i]) + dividend[i] + rh;
- rh = r / RADIX;
- remainder[i] = r - rh * RADIX;
- }
-
- /* Lets analyze what we have here */
- if (rh == 0) {
- // trial is too big -- negative remainder
- max = trial-1;
- continue;
- }
-
- /* find out how far off the remainder is telling us we are */
- minadj = (remainder[mant.length] * RADIX)+remainder[mant.length-1];
- minadj /= divisor.mant[mant.length-1] + 1;
-
- if (minadj >= 2) {
- min = trial+minadj; // update the minimum
- continue;
- }
-
- /* May have a good one here, check more thoroughly. Basically
- its a good one if it is less than the divisor */
- trialgood = false; // assume false
- for (int i = mant.length - 1; i >= 0; i--) {
- if (divisor.mant[i] > remainder[i]) {
- trialgood = true;
- }
- if (divisor.mant[i] < remainder[i]) {
- break;
- }
- }
-
- if (remainder[mant.length] != 0) {
- trialgood = false;
- }
-
- if (trialgood == false) {
- min = trial+1;
- }
- }
-
- /* Great we have a digit! */
- quotient[qd] = trial;
- if (trial != 0 || nsqd != 0) {
- nsqd++;
- }
-
- if (field.getRoundingMode() == DfpField.RoundingMode.ROUND_DOWN && nsqd == mant.length) {
- // We have enough for this mode
- break;
- }
-
- if (nsqd > mant.length) {
- // We have enough digits
- break;
- }
-
- /* move the remainder into the dividend while left shifting */
- dividend[0] = 0;
- for (int i = 0; i < mant.length; i++) {
- dividend[i + 1] = remainder[i];
- }
- }
-
- /* Find the most sig digit */
- md = mant.length; // default
- for (int i = mant.length + 1; i >= 0; i--) {
- if (quotient[i] != 0) {
- md = i;
- break;
- }
- }
-
- /* Copy the digits into the result */
- for (int i=0; i<mant.length; i++) {
- result.mant[mant.length-i-1] = quotient[md-i];
- }
-
- /* Fixup the exponent. */
- result.exp = exp - divisor.exp + md - mant.length;
- result.sign = (byte) ((sign == divisor.sign) ? 1 : -1);
-
- if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
- result.exp = 0;
- }
-
- if (md > (mant.length-1)) {
- excp = result.round(quotient[md-mant.length]);
- } else {
- excp = result.round(0);
- }
-
- if (excp != 0) {
- result = dotrap(excp, DIVIDE_TRAP, divisor, result);
- }
-
- return result;
- }
-
- /** Divide by a single digit less than radix.
- * Special case, so there are speed advantages. 0 <= divisor < radix
- * @param divisor divisor
- * @return quotient of this by divisor
- */
- public Dfp divide(int divisor) {
-
- // Handle special cases
- if (nans != FINITE) {
- if (isNaN()) {
- return this;
- }
-
- if (nans == INFINITE) {
- return newInstance(this);
- }
- }
-
- // Test for divide by zero
- if (divisor == 0) {
- field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
- Dfp result = newInstance(getZero());
- result.sign = sign;
- result.nans = INFINITE;
- result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, getZero(), result);
- return result;
- }
-
- // range check divisor
- if (divisor < 0 || divisor >= RADIX) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- Dfp result = newInstance(getZero());
- result.nans = QNAN;
- result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, result, result);
- return result;
- }
-
- Dfp result = newInstance(this);
-
- int rl = 0;
- for (int i = mant.length-1; i >= 0; i--) {
- final int r = rl*RADIX + result.mant[i];
- final int rh = r / divisor;
- rl = r - rh * divisor;
- result.mant[i] = rh;
- }
-
- if (result.mant[mant.length-1] == 0) {
- // normalize
- result.shiftLeft();
- final int r = rl * RADIX; // compute the next digit and put it in
- final int rh = r / divisor;
- rl = r - rh * divisor;
- result.mant[0] = rh;
- }
-
- final int excp = result.round(rl * RADIX / divisor); // do the rounding
- if (excp != 0) {
- result = dotrap(excp, DIVIDE_TRAP, result, result);
- }
-
- return result;
-
- }
-
- /** {@inheritDoc} */
- public Dfp reciprocal() {
- return field.getOne().divide(this);
- }
-
- /** Compute the square root.
- * @return square root of the instance
- * @since 3.2
- */
- public Dfp sqrt() {
-
- // check for unusual cases
- if (nans == FINITE && mant[mant.length-1] == 0) {
- // if zero
- return newInstance(this);
- }
-
- if (nans != FINITE) {
- if (nans == INFINITE && sign == 1) {
- // if positive infinity
- return newInstance(this);
- }
-
- if (nans == QNAN) {
- return newInstance(this);
- }
-
- if (nans == SNAN) {
- Dfp result;
-
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- result = newInstance(this);
- result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
- return result;
- }
- }
-
- if (sign == -1) {
- // if negative
- Dfp result;
-
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- result = newInstance(this);
- result.nans = QNAN;
- result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
- return result;
- }
-
- Dfp x = newInstance(this);
-
- /* Lets make a reasonable guess as to the size of the square root */
- if (x.exp < -1 || x.exp > 1) {
- x.exp = this.exp / 2;
- }
-
- /* Coarsely estimate the mantissa */
- switch (x.mant[mant.length-1] / 2000) {
- case 0:
- x.mant[mant.length-1] = x.mant[mant.length-1]/2+1;
- break;
- case 2:
- x.mant[mant.length-1] = 1500;
- break;
- case 3:
- x.mant[mant.length-1] = 2200;
- break;
- default:
- x.mant[mant.length-1] = 3000;
- }
-
- Dfp dx = newInstance(x);
-
- /* Now that we have the first pass estimate, compute the rest
- by the formula dx = (y - x*x) / (2x); */
-
- Dfp px = getZero();
- Dfp ppx = getZero();
- while (x.unequal(px)) {
- dx = newInstance(x);
- dx.sign = -1;
- dx = dx.add(this.divide(x));
- dx = dx.divide(2);
- ppx = px;
- px = x;
- x = x.add(dx);
-
- if (x.equals(ppx)) {
- // alternating between two values
- break;
- }
-
- // if dx is zero, break. Note testing the most sig digit
- // is a sufficient test since dx is normalized
- if (dx.mant[mant.length-1] == 0) {
- break;
- }
- }
-
- return x;
-
- }
-
- /** Get a string representation of the instance.
- * @return string representation of the instance
- */
- @Override
- public String toString() {
- if (nans != FINITE) {
- // if non-finite exceptional cases
- if (nans == INFINITE) {
- return (sign < 0) ? NEG_INFINITY_STRING : POS_INFINITY_STRING;
- } else {
- return NAN_STRING;
- }
- }
-
- if (exp > mant.length || exp < -1) {
- return dfp2sci();
- }
-
- return dfp2string();
-
- }
-
- /** Convert an instance to a string using scientific notation.
- * @return string representation of the instance in scientific notation
- */
- protected String dfp2sci() {
- char rawdigits[] = new char[mant.length * 4];
- char outputbuffer[] = new char[mant.length * 4 + 20];
- int p;
- int q;
- int e;
- int ae;
- int shf;
-
- // Get all the digits
- p = 0;
- for (int i = mant.length - 1; i >= 0; i--) {
- rawdigits[p++] = (char) ((mant[i] / 1000) + '0');
- rawdigits[p++] = (char) (((mant[i] / 100) %10) + '0');
- rawdigits[p++] = (char) (((mant[i] / 10) % 10) + '0');
- rawdigits[p++] = (char) (((mant[i]) % 10) + '0');
- }
-
- // Find the first non-zero one
- for (p = 0; p < rawdigits.length; p++) {
- if (rawdigits[p] != '0') {
- break;
- }
- }
- shf = p;
-
- // Now do the conversion
- q = 0;
- if (sign == -1) {
- outputbuffer[q++] = '-';
- }
-
- if (p != rawdigits.length) {
- // there are non zero digits...
- outputbuffer[q++] = rawdigits[p++];
- outputbuffer[q++] = '.';
-
- while (p<rawdigits.length) {
- outputbuffer[q++] = rawdigits[p++];
- }
- } else {
- outputbuffer[q++] = '0';
- outputbuffer[q++] = '.';
- outputbuffer[q++] = '0';
- outputbuffer[q++] = 'e';
- outputbuffer[q++] = '0';
- return new String(outputbuffer, 0, 5);
- }
-
- outputbuffer[q++] = 'e';
-
- // Find the msd of the exponent
-
- e = exp * 4 - shf - 1;
- ae = e;
- if (e < 0) {
- ae = -e;
- }
-
- // Find the largest p such that p < e
- for (p = 1000000000; p > ae; p /= 10) {
- // nothing to do
- }
-
- if (e < 0) {
- outputbuffer[q++] = '-';
- }
-
- while (p > 0) {
- outputbuffer[q++] = (char)(ae / p + '0');
- ae %= p;
- p /= 10;
- }
-
- return new String(outputbuffer, 0, q);
-
- }
-
- /** Convert an instance to a string using normal notation.
- * @return string representation of the instance in normal notation
- */
- protected String dfp2string() {
- char buffer[] = new char[mant.length*4 + 20];
- int p = 1;
- int q;
- int e = exp;
- boolean pointInserted = false;
-
- buffer[0] = ' ';
-
- if (e <= 0) {
- buffer[p++] = '0';
- buffer[p++] = '.';
- pointInserted = true;
- }
-
- while (e < 0) {
- buffer[p++] = '0';
- buffer[p++] = '0';
- buffer[p++] = '0';
- buffer[p++] = '0';
- e++;
- }
-
- for (int i = mant.length - 1; i >= 0; i--) {
- buffer[p++] = (char) ((mant[i] / 1000) + '0');
- buffer[p++] = (char) (((mant[i] / 100) % 10) + '0');
- buffer[p++] = (char) (((mant[i] / 10) % 10) + '0');
- buffer[p++] = (char) (((mant[i]) % 10) + '0');
- if (--e == 0) {
- buffer[p++] = '.';
- pointInserted = true;
- }
- }
-
- while (e > 0) {
- buffer[p++] = '0';
- buffer[p++] = '0';
- buffer[p++] = '0';
- buffer[p++] = '0';
- e--;
- }
-
- if (!pointInserted) {
- // Ensure we have a radix point!
- buffer[p++] = '.';
- }
-
- // Suppress leading zeros
- q = 1;
- while (buffer[q] == '0') {
- q++;
- }
- if (buffer[q] == '.') {
- q--;
- }
-
- // Suppress trailing zeros
- while (buffer[p-1] == '0') {
- p--;
- }
-
- // Insert sign
- if (sign < 0) {
- buffer[--q] = '-';
- }
-
- return new String(buffer, q, p - q);
-
- }
-
- /** Raises a trap. This does not set the corresponding flag however.
- * @param type the trap type
- * @param what - name of routine trap occurred in
- * @param oper - input operator to function
- * @param result - the result computed prior to the trap
- * @return The suggested return value from the trap handler
- */
- public Dfp dotrap(int type, String what, Dfp oper, Dfp result) {
- Dfp def = result;
-
- switch (type) {
- case DfpField.FLAG_INVALID:
- def = newInstance(getZero());
- def.sign = result.sign;
- def.nans = QNAN;
- break;
-
- case DfpField.FLAG_DIV_ZERO:
- if (nans == FINITE && mant[mant.length-1] != 0) {
- // normal case, we are finite, non-zero
- def = newInstance(getZero());
- def.sign = (byte)(sign*oper.sign);
- def.nans = INFINITE;
- }
-
- if (nans == FINITE && mant[mant.length-1] == 0) {
- // 0/0
- def = newInstance(getZero());
- def.nans = QNAN;
- }
-
- if (nans == INFINITE || nans == QNAN) {
- def = newInstance(getZero());
- def.nans = QNAN;
- }
-
- if (nans == INFINITE || nans == SNAN) {
- def = newInstance(getZero());
- def.nans = QNAN;
- }
- break;
-
- case DfpField.FLAG_UNDERFLOW:
- if ( (result.exp+mant.length) < MIN_EXP) {
- def = newInstance(getZero());
- def.sign = result.sign;
- } else {
- def = newInstance(result); // gradual underflow
- }
- result.exp += ERR_SCALE;
- break;
-
- case DfpField.FLAG_OVERFLOW:
- result.exp -= ERR_SCALE;
- def = newInstance(getZero());
- def.sign = result.sign;
- def.nans = INFINITE;
- break;
-
- default: def = result; break;
- }
-
- return trap(type, what, oper, def, result);
-
- }
-
- /** Trap handler. Subclasses may override this to provide trap
- * functionality per IEEE 854-1987.
- *
- * @param type The exception type - e.g. FLAG_OVERFLOW
- * @param what The name of the routine we were in e.g. divide()
- * @param oper An operand to this function if any
- * @param def The default return value if trap not enabled
- * @param result The result that is specified to be delivered per
- * IEEE 854, if any
- * @return the value that should be return by the operation triggering the trap
- */
- protected Dfp trap(int type, String what, Dfp oper, Dfp def, Dfp result) {
- return def;
- }
-
- /** Returns the type - one of FINITE, INFINITE, SNAN, QNAN.
- * @return type of the number
- */
- public int classify() {
- return nans;
- }
-
- /** Creates an instance that is the same as x except that it has the sign of y.
- * abs(x) = dfp.copysign(x, dfp.one)
- * @param x number to get the value from
- * @param y number to get the sign from
- * @return a number with the value of x and the sign of y
- */
- public static Dfp copysign(final Dfp x, final Dfp y) {
- Dfp result = x.newInstance(x);
- result.sign = y.sign;
- return result;
- }
-
- /** Returns the next number greater than this one in the direction of x.
- * If this==x then simply returns this.
- * @param x direction where to look at
- * @return closest number next to instance in the direction of x
- */
- public Dfp nextAfter(final Dfp x) {
-
- // make sure we don't mix number with different precision
- if (field.getRadixDigits() != x.field.getRadixDigits()) {
- field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = newInstance(getZero());
- result.nans = QNAN;
- return dotrap(DfpField.FLAG_INVALID, NEXT_AFTER_TRAP, x, result);
- }
-
- // if this is greater than x
- boolean up = false;
- if (this.lessThan(x)) {
- up = true;
- }
-
- if (compare(this, x) == 0) {
- return newInstance(x);
- }
-
- if (lessThan(getZero())) {
- up = !up;
- }
-
- final Dfp inc;
- Dfp result;
- if (up) {
- inc = newInstance(getOne());
- inc.exp = this.exp-mant.length+1;
- inc.sign = this.sign;
-
- if (this.equals(getZero())) {
- inc.exp = MIN_EXP-mant.length;
- }
-
- result = add(inc);
- } else {
- inc = newInstance(getOne());
- inc.exp = this.exp;
- inc.sign = this.sign;
-
- if (this.equals(inc)) {
- inc.exp = this.exp-mant.length;
- } else {
- inc.exp = this.exp-mant.length+1;
- }
-
- if (this.equals(getZero())) {
- inc.exp = MIN_EXP-mant.length;
- }
-
- result = this.subtract(inc);
- }
-
- if (result.classify() == INFINITE && this.classify() != INFINITE) {
- field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
- result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
- }
-
- if (result.equals(getZero()) && this.equals(getZero()) == false) {
- field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
- result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
- }
-
- return result;
-
- }
-
- /** Convert the instance into a double.
- * @return a double approximating the instance
- * @see #toSplitDouble()
- */
- public double toDouble() {
-
- if (isInfinite()) {
- if (lessThan(getZero())) {
- return Double.NEGATIVE_INFINITY;
- } else {
- return Double.POSITIVE_INFINITY;
- }
- }
-
- if (isNaN()) {
- return Double.NaN;
- }
-
- Dfp y = this;
- boolean negate = false;
- int cmp0 = compare(this, getZero());
- if (cmp0 == 0) {
- return sign < 0 ? -0.0 : +0.0;
- } else if (cmp0 < 0) {
- y = negate();
- negate = true;
- }
-
- /* Find the exponent, first estimate by integer log10, then adjust.
- Should be faster than doing a natural logarithm. */
- int exponent = (int)(y.intLog10() * 3.32);
- if (exponent < 0) {
- exponent--;
- }
-
- Dfp tempDfp = DfpMath.pow(getTwo(), exponent);
- while (tempDfp.lessThan(y) || tempDfp.equals(y)) {
- tempDfp = tempDfp.multiply(2);
- exponent++;
- }
- exponent--;
-
- /* We have the exponent, now work on the mantissa */
-
- y = y.divide(DfpMath.pow(getTwo(), exponent));
- if (exponent > -1023) {
- y = y.subtract(getOne());
- }
-
- if (exponent < -1074) {
- return 0;
- }
-
- if (exponent > 1023) {
- return negate ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
- }
-
-
- y = y.multiply(newInstance(4503599627370496l)).rint();
- String str = y.toString();
- str = str.substring(0, str.length()-1);
- long mantissa = Long.parseLong(str);
-
- if (mantissa == 4503599627370496L) {
- // Handle special case where we round up to next power of two
- mantissa = 0;
- exponent++;
- }
-
- /* Its going to be subnormal, so make adjustments */
- if (exponent <= -1023) {
- exponent--;
- }
-
- while (exponent < -1023) {
- exponent++;
- mantissa >>>= 1;
- }
-
- long bits = mantissa | ((exponent + 1023L) << 52);
- double x = Double.longBitsToDouble(bits);
-
- if (negate) {
- x = -x;
- }
-
- return x;
-
- }
-
- /** Convert the instance into a split double.
- * @return an array of two doubles which sum represent the instance
- * @see #toDouble()
- */
- public double[] toSplitDouble() {
- double split[] = new double[2];
- long mask = 0xffffffffc0000000L;
-
- split[0] = Double.longBitsToDouble(Double.doubleToLongBits(toDouble()) & mask);
- split[1] = subtract(newInstance(split[0])).toDouble();
-
- return split;
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public double getReal() {
- return toDouble();
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp add(final double a) {
- return add(newInstance(a));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp subtract(final double a) {
- return subtract(newInstance(a));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp multiply(final double a) {
- return multiply(newInstance(a));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp divide(final double a) {
- return divide(newInstance(a));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp remainder(final double a) {
- return remainder(newInstance(a));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public long round() {
- return FastMath.round(toDouble());
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp signum() {
- if (isNaN() || isZero()) {
- return this;
- } else {
- return newInstance(sign > 0 ? +1 : -1);
- }
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp copySign(final Dfp s) {
- if ((sign >= 0 && s.sign >= 0) || (sign < 0 && s.sign < 0)) { // Sign is currently OK
- return this;
- }
- return negate(); // flip sign
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp copySign(final double s) {
- long sb = Double.doubleToLongBits(s);
- if ((sign >= 0 && sb >= 0) || (sign < 0 && sb < 0)) { // Sign is currently OK
- return this;
- }
- return negate(); // flip sign
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp scalb(final int n) {
- return multiply(DfpMath.pow(getTwo(), n));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp hypot(final Dfp y) {
- return multiply(this).add(y.multiply(y)).sqrt();
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp cbrt() {
- return rootN(3);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp rootN(final int n) {
- return (sign >= 0) ?
- DfpMath.pow(this, getOne().divide(n)) :
- DfpMath.pow(negate(), getOne().divide(n)).negate();
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp pow(final double p) {
- return DfpMath.pow(this, newInstance(p));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp pow(final int n) {
- return DfpMath.pow(this, n);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp pow(final Dfp e) {
- return DfpMath.pow(this, e);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp exp() {
- return DfpMath.exp(this);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp expm1() {
- return DfpMath.exp(this).subtract(getOne());
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp log() {
- return DfpMath.log(this);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp log1p() {
- return DfpMath.log(this.add(getOne()));
- }
-
-// TODO: deactivate this implementation (and return type) in 4.0
- /** Get the exponent of the greatest power of 10 that is less than or equal to abs(this).
- * @return integer base 10 logarithm
- * @deprecated as of 3.2, replaced by {@link #intLog10()}, in 4.0 the return type
- * will be changed to Dfp
- */
- @Deprecated
- public int log10() {
- return intLog10();
- }
-
-// TODO: activate this implementation (and return type) in 4.0
-// /** {@inheritDoc}
-// * @since 3.2
-// */
-// public Dfp log10() {
-// return DfpMath.log(this).divide(DfpMath.log(newInstance(10)));
-// }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp cos() {
- return DfpMath.cos(this);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp sin() {
- return DfpMath.sin(this);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp tan() {
- return DfpMath.tan(this);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp acos() {
- return DfpMath.acos(this);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp asin() {
- return DfpMath.asin(this);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp atan() {
- return DfpMath.atan(this);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp atan2(final Dfp x)
- throws DimensionMismatchException {
-
- // compute r = sqrt(x^2+y^2)
- final Dfp r = x.multiply(x).add(multiply(this)).sqrt();
-
- if (x.sign >= 0) {
-
- // compute atan2(y, x) = 2 atan(y / (r + x))
- return getTwo().multiply(divide(r.add(x)).atan());
-
- } else {
-
- // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
- final Dfp tmp = getTwo().multiply(divide(r.subtract(x)).atan());
- final Dfp pmPi = newInstance((tmp.sign <= 0) ? -FastMath.PI : FastMath.PI);
- return pmPi.subtract(tmp);
-
- }
-
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp cosh() {
- return DfpMath.exp(this).add(DfpMath.exp(negate())).divide(2);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp sinh() {
- return DfpMath.exp(this).subtract(DfpMath.exp(negate())).divide(2);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp tanh() {
- final Dfp ePlus = DfpMath.exp(this);
- final Dfp eMinus = DfpMath.exp(negate());
- return ePlus.subtract(eMinus).divide(ePlus.add(eMinus));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp acosh() {
- return multiply(this).subtract(getOne()).sqrt().add(this).log();
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp asinh() {
- return multiply(this).add(getOne()).sqrt().add(this).log();
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp atanh() {
- return getOne().add(this).divide(getOne().subtract(this)).log().divide(2);
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp linearCombination(final Dfp[] a, final Dfp[] b)
- throws DimensionMismatchException {
- if (a.length != b.length) {
- throw new DimensionMismatchException(a.length, b.length);
- }
- Dfp r = getZero();
- for (int i = 0; i < a.length; ++i) {
- r = r.add(a[i].multiply(b[i]));
- }
- return r;
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp linearCombination(final double[] a, final Dfp[] b)
- throws DimensionMismatchException {
- if (a.length != b.length) {
- throw new DimensionMismatchException(a.length, b.length);
- }
- Dfp r = getZero();
- for (int i = 0; i < a.length; ++i) {
- r = r.add(b[i].multiply(a[i]));
- }
- return r;
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2) {
- return a1.multiply(b1).add(a2.multiply(b2));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2) {
- return b1.multiply(a1).add(b2.multiply(a2));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp linearCombination(final Dfp a1, final Dfp b1,
- final Dfp a2, final Dfp b2,
- final Dfp a3, final Dfp b3) {
- return a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp linearCombination(final double a1, final Dfp b1,
- final double a2, final Dfp b2,
- final double a3, final Dfp b3) {
- return b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2,
- final Dfp a3, final Dfp b3, final Dfp a4, final Dfp b4) {
- return a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3)).add(a4.multiply(b4));
- }
-
- /** {@inheritDoc}
- * @since 3.2
- */
- public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2,
- final double a3, final Dfp b3, final double a4, final Dfp b4) {
- return b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3)).add(b4.multiply(a4));
- }
-
-}