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:05 UTC
[35/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/DfpDec.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/dfp/DfpDec.java b/src/main/java/org/apache/commons/math3/dfp/DfpDec.java
deleted file mode 100644
index 20875c0..0000000
--- a/src/main/java/org/apache/commons/math3/dfp/DfpDec.java
+++ /dev/null
@@ -1,368 +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;
-
-/** Subclass of {@link Dfp} which hides the radix-10000 artifacts of the superclass.
- * This should give outward appearances of being a decimal number with DIGITS*4-3
- * decimal digits. This class can be subclassed to appear to be an arbitrary number
- * of decimal digits less than DIGITS*4-3.
- * @since 2.2
- */
-public class DfpDec extends Dfp {
-
- /** Makes an instance with a value of zero.
- * @param factory factory linked to this instance
- */
- protected DfpDec(final DfpField factory) {
- super(factory);
- }
-
- /** Create an instance from a byte value.
- * @param factory factory linked to this instance
- * @param x value to convert to an instance
- */
- protected DfpDec(final DfpField factory, byte x) {
- super(factory, x);
- }
-
- /** Create an instance from an int value.
- * @param factory factory linked to this instance
- * @param x value to convert to an instance
- */
- protected DfpDec(final DfpField factory, int x) {
- super(factory, x);
- }
-
- /** Create an instance from a long value.
- * @param factory factory linked to this instance
- * @param x value to convert to an instance
- */
- protected DfpDec(final DfpField factory, long x) {
- super(factory, x);
- }
-
- /** Create an instance from a double value.
- * @param factory factory linked to this instance
- * @param x value to convert to an instance
- */
- protected DfpDec(final DfpField factory, double x) {
- super(factory, x);
- round(0);
- }
-
- /** Copy constructor.
- * @param d instance to copy
- */
- public DfpDec(final Dfp d) {
- super(d);
- round(0);
- }
-
- /** Create an instance from a String representation.
- * @param factory factory linked to this instance
- * @param s string representation of the instance
- */
- protected DfpDec(final DfpField factory, final String s) {
- super(factory, s);
- round(0);
- }
-
- /** Creates an instance with a non-finite value.
- * @param factory factory linked to this instance
- * @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 DfpDec(final DfpField factory, final byte sign, final byte nans) {
- super(factory, sign, nans);
- }
-
- /** {@inheritDoc} */
- @Override
- public Dfp newInstance() {
- return new DfpDec(getField());
- }
-
- /** {@inheritDoc} */
- @Override
- public Dfp newInstance(final byte x) {
- return new DfpDec(getField(), x);
- }
-
- /** {@inheritDoc} */
- @Override
- public Dfp newInstance(final int x) {
- return new DfpDec(getField(), x);
- }
-
- /** {@inheritDoc} */
- @Override
- public Dfp newInstance(final long x) {
- return new DfpDec(getField(), x);
- }
-
- /** {@inheritDoc} */
- @Override
- public Dfp newInstance(final double x) {
- return new DfpDec(getField(), x);
- }
-
- /** {@inheritDoc} */
- @Override
- public Dfp newInstance(final Dfp d) {
-
- // make sure we don't mix number with different precision
- if (getField().getRadixDigits() != d.getField().getRadixDigits()) {
- getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = newInstance(getZero());
- result.nans = QNAN;
- return dotrap(DfpField.FLAG_INVALID, "newInstance", d, result);
- }
-
- return new DfpDec(d);
-
- }
-
- /** {@inheritDoc} */
- @Override
- public Dfp newInstance(final String s) {
- return new DfpDec(getField(), s);
- }
-
- /** {@inheritDoc} */
- @Override
- public Dfp newInstance(final byte sign, final byte nans) {
- return new DfpDec(getField(), sign, nans);
- }
-
- /** Get the number of decimal digits this class is going to represent.
- * Default implementation returns {@link #getRadixDigits()}*4-3. Subclasses can
- * override this to return something less.
- * @return number of decimal digits this class is going to represent
- */
- protected int getDecimalDigits() {
- return getRadixDigits() * 4 - 3;
- }
-
- /** {@inheritDoc} */
- @Override
- protected int round(int in) {
-
- int msb = mant[mant.length-1];
- if (msb == 0) {
- // special case -- this == zero
- return 0;
- }
-
- int cmaxdigits = mant.length * 4;
- int lsbthreshold = 1000;
- while (lsbthreshold > msb) {
- lsbthreshold /= 10;
- cmaxdigits --;
- }
-
-
- final int digits = getDecimalDigits();
- final int lsbshift = cmaxdigits - digits;
- final int lsd = lsbshift / 4;
-
- lsbthreshold = 1;
- for (int i = 0; i < lsbshift % 4; i++) {
- lsbthreshold *= 10;
- }
-
- final int lsb = mant[lsd];
-
- if (lsbthreshold <= 1 && digits == 4 * mant.length - 3) {
- return super.round(in);
- }
-
- int discarded = in; // not looking at this after this point
- final int n;
- if (lsbthreshold == 1) {
- // look to the next digit for rounding
- n = (mant[lsd-1] / 1000) % 10;
- mant[lsd-1] %= 1000;
- discarded |= mant[lsd-1];
- } else {
- n = (lsb * 10 / lsbthreshold) % 10;
- discarded |= lsb % (lsbthreshold/10);
- }
-
- for (int i = 0; i < lsd; i++) {
- discarded |= mant[i]; // need to know if there are any discarded bits
- mant[i] = 0;
- }
-
- mant[lsd] = lsb / lsbthreshold * lsbthreshold;
-
- final boolean inc;
- switch (getField().getRoundingMode()) {
- case ROUND_DOWN:
- inc = false;
- break;
-
- case ROUND_UP:
- inc = (n != 0) || (discarded != 0); // round up if n!=0
- break;
-
- case ROUND_HALF_UP:
- inc = n >= 5; // round half up
- break;
-
- case ROUND_HALF_DOWN:
- inc = n > 5; // round half down
- break;
-
- case ROUND_HALF_EVEN:
- inc = (n > 5) ||
- (n == 5 && discarded != 0) ||
- (n == 5 && discarded == 0 && ((lsb / lsbthreshold) & 1) == 1); // round half-even
- break;
-
- case ROUND_HALF_ODD:
- inc = (n > 5) ||
- (n == 5 && discarded != 0) ||
- (n == 5 && discarded == 0 && ((lsb / lsbthreshold) & 1) == 0); // round half-odd
- break;
-
- case ROUND_CEIL:
- inc = (sign == 1) && (n != 0 || discarded != 0); // round ceil
- break;
-
- case ROUND_FLOOR:
- default:
- inc = (sign == -1) && (n != 0 || discarded != 0); // round floor
- break;
- }
-
- if (inc) {
- // increment if necessary
- int rh = lsbthreshold;
- for (int i = lsd; i < mant.length; i++) {
- final int r = mant[i] + rh;
- rh = r / RADIX;
- mant[i] = r % 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
- getField().setIEEEFlagsBits(DfpField.FLAG_UNDERFLOW);
- return DfpField.FLAG_UNDERFLOW;
- }
-
- if (exp > MAX_EXP) {
- // Overflow
- getField().setIEEEFlagsBits(DfpField.FLAG_OVERFLOW);
- return DfpField.FLAG_OVERFLOW;
- }
-
- if (n != 0 || discarded != 0) {
- // Inexact
- getField().setIEEEFlagsBits(DfpField.FLAG_INEXACT);
- return DfpField.FLAG_INEXACT;
- }
- return 0;
- }
-
- /** {@inheritDoc} */
- @Override
- public Dfp nextAfter(Dfp x) {
-
- final String trapName = "nextAfter";
-
- // make sure we don't mix number with different precision
- if (getField().getRadixDigits() != x.getField().getRadixDigits()) {
- getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = newInstance(getZero());
- result.nans = QNAN;
- return dotrap(DfpField.FLAG_INVALID, trapName, x, result);
- }
-
- boolean up = false;
- Dfp result;
- Dfp inc;
-
- // if this is greater than x
- if (this.lessThan(x)) {
- up = true;
- }
-
- if (equals(x)) {
- return newInstance(x);
- }
-
- if (lessThan(getZero())) {
- up = !up;
- }
-
- if (up) {
- inc = power10(intLog10() - getDecimalDigits() + 1);
- inc = copysign(inc, this);
-
- if (this.equals(getZero())) {
- inc = power10K(MIN_EXP-mant.length-1);
- }
-
- if (inc.equals(getZero())) {
- result = copysign(newInstance(getZero()), this);
- } else {
- result = add(inc);
- }
- } else {
- inc = power10(intLog10());
- inc = copysign(inc, this);
-
- if (this.equals(inc)) {
- inc = inc.divide(power10(getDecimalDigits()));
- } else {
- inc = inc.divide(power10(getDecimalDigits() - 1));
- }
-
- if (this.equals(getZero())) {
- inc = power10K(MIN_EXP-mant.length-1);
- }
-
- if (inc.equals(getZero())) {
- result = copysign(newInstance(getZero()), this);
- } else {
- result = subtract(inc);
- }
- }
-
- if (result.classify() == INFINITE && this.classify() != INFINITE) {
- getField().setIEEEFlagsBits(DfpField.FLAG_INEXACT);
- result = dotrap(DfpField.FLAG_INEXACT, trapName, x, result);
- }
-
- if (result.equals(getZero()) && this.equals(getZero()) == false) {
- getField().setIEEEFlagsBits(DfpField.FLAG_INEXACT);
- result = dotrap(DfpField.FLAG_INEXACT, trapName, x, result);
- }
-
- return result;
- }
-
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/dfp/DfpField.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/dfp/DfpField.java b/src/main/java/org/apache/commons/math3/dfp/DfpField.java
deleted file mode 100644
index fcdec82..0000000
--- a/src/main/java/org/apache/commons/math3/dfp/DfpField.java
+++ /dev/null
@@ -1,757 +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 org.apache.commons.math3.Field;
-import org.apache.commons.math3.FieldElement;
-
-/** Field for Decimal floating point instances.
- * @since 2.2
- */
-public class DfpField implements Field<Dfp> {
-
- /** Enumerate for rounding modes. */
- public enum RoundingMode {
-
- /** Rounds toward zero (truncation). */
- ROUND_DOWN,
-
- /** Rounds away from zero if discarded digit is non-zero. */
- ROUND_UP,
-
- /** Rounds towards nearest unless both are equidistant in which case it rounds away from zero. */
- ROUND_HALF_UP,
-
- /** Rounds towards nearest unless both are equidistant in which case it rounds toward zero. */
- ROUND_HALF_DOWN,
-
- /** Rounds towards nearest unless both are equidistant in which case it rounds toward the even neighbor.
- * This is the default as specified by IEEE 854-1987
- */
- ROUND_HALF_EVEN,
-
- /** Rounds towards nearest unless both are equidistant in which case it rounds toward the odd neighbor. */
- ROUND_HALF_ODD,
-
- /** Rounds towards positive infinity. */
- ROUND_CEIL,
-
- /** Rounds towards negative infinity. */
- ROUND_FLOOR;
-
- }
-
- /** IEEE 854-1987 flag for invalid operation. */
- public static final int FLAG_INVALID = 1;
-
- /** IEEE 854-1987 flag for division by zero. */
- public static final int FLAG_DIV_ZERO = 2;
-
- /** IEEE 854-1987 flag for overflow. */
- public static final int FLAG_OVERFLOW = 4;
-
- /** IEEE 854-1987 flag for underflow. */
- public static final int FLAG_UNDERFLOW = 8;
-
- /** IEEE 854-1987 flag for inexact result. */
- public static final int FLAG_INEXACT = 16;
-
- /** High precision string representation of √2. */
- private static String sqr2String;
-
- // Note: the static strings are set up (once) by the ctor and @GuardedBy("DfpField.class")
-
- /** High precision string representation of √2 / 2. */
- private static String sqr2ReciprocalString;
-
- /** High precision string representation of √3. */
- private static String sqr3String;
-
- /** High precision string representation of √3 / 3. */
- private static String sqr3ReciprocalString;
-
- /** High precision string representation of π. */
- private static String piString;
-
- /** High precision string representation of e. */
- private static String eString;
-
- /** High precision string representation of ln(2). */
- private static String ln2String;
-
- /** High precision string representation of ln(5). */
- private static String ln5String;
-
- /** High precision string representation of ln(10). */
- private static String ln10String;
-
- /** The number of radix digits.
- * Note these depend on the radix which is 10000 digits,
- * so each one is equivalent to 4 decimal digits.
- */
- private final int radixDigits;
-
- /** A {@link Dfp} with value 0. */
- private final Dfp zero;
-
- /** A {@link Dfp} with value 1. */
- private final Dfp one;
-
- /** A {@link Dfp} with value 2. */
- private final Dfp two;
-
- /** A {@link Dfp} with value √2. */
- private final Dfp sqr2;
-
- /** A two elements {@link Dfp} array with value √2 split in two pieces. */
- private final Dfp[] sqr2Split;
-
- /** A {@link Dfp} with value √2 / 2. */
- private final Dfp sqr2Reciprocal;
-
- /** A {@link Dfp} with value √3. */
- private final Dfp sqr3;
-
- /** A {@link Dfp} with value √3 / 3. */
- private final Dfp sqr3Reciprocal;
-
- /** A {@link Dfp} with value π. */
- private final Dfp pi;
-
- /** A two elements {@link Dfp} array with value π split in two pieces. */
- private final Dfp[] piSplit;
-
- /** A {@link Dfp} with value e. */
- private final Dfp e;
-
- /** A two elements {@link Dfp} array with value e split in two pieces. */
- private final Dfp[] eSplit;
-
- /** A {@link Dfp} with value ln(2). */
- private final Dfp ln2;
-
- /** A two elements {@link Dfp} array with value ln(2) split in two pieces. */
- private final Dfp[] ln2Split;
-
- /** A {@link Dfp} with value ln(5). */
- private final Dfp ln5;
-
- /** A two elements {@link Dfp} array with value ln(5) split in two pieces. */
- private final Dfp[] ln5Split;
-
- /** A {@link Dfp} with value ln(10). */
- private final Dfp ln10;
-
- /** Current rounding mode. */
- private RoundingMode rMode;
-
- /** IEEE 854-1987 signals. */
- private int ieeeFlags;
-
- /** Create a factory for the specified number of radix digits.
- * <p>
- * Note that since the {@link Dfp} class uses 10000 as its radix, each radix
- * digit is equivalent to 4 decimal digits. This implies that asking for
- * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in
- * all cases.
- * </p>
- * @param decimalDigits minimal number of decimal digits.
- */
- public DfpField(final int decimalDigits) {
- this(decimalDigits, true);
- }
-
- /** Create a factory for the specified number of radix digits.
- * <p>
- * Note that since the {@link Dfp} class uses 10000 as its radix, each radix
- * digit is equivalent to 4 decimal digits. This implies that asking for
- * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in
- * all cases.
- * </p>
- * @param decimalDigits minimal number of decimal digits
- * @param computeConstants if true, the transcendental constants for the given precision
- * must be computed (setting this flag to false is RESERVED for the internal recursive call)
- */
- private DfpField(final int decimalDigits, final boolean computeConstants) {
-
- this.radixDigits = (decimalDigits < 13) ? 4 : (decimalDigits + 3) / 4;
- this.rMode = RoundingMode.ROUND_HALF_EVEN;
- this.ieeeFlags = 0;
- this.zero = new Dfp(this, 0);
- this.one = new Dfp(this, 1);
- this.two = new Dfp(this, 2);
-
- if (computeConstants) {
- // set up transcendental constants
- synchronized (DfpField.class) {
-
- // as a heuristic to circumvent Table-Maker's Dilemma, we set the string
- // representation of the constants to be at least 3 times larger than the
- // number of decimal digits, also as an attempt to really compute these
- // constants only once, we set a minimum number of digits
- computeStringConstants((decimalDigits < 67) ? 200 : (3 * decimalDigits));
-
- // set up the constants at current field accuracy
- sqr2 = new Dfp(this, sqr2String);
- sqr2Split = split(sqr2String);
- sqr2Reciprocal = new Dfp(this, sqr2ReciprocalString);
- sqr3 = new Dfp(this, sqr3String);
- sqr3Reciprocal = new Dfp(this, sqr3ReciprocalString);
- pi = new Dfp(this, piString);
- piSplit = split(piString);
- e = new Dfp(this, eString);
- eSplit = split(eString);
- ln2 = new Dfp(this, ln2String);
- ln2Split = split(ln2String);
- ln5 = new Dfp(this, ln5String);
- ln5Split = split(ln5String);
- ln10 = new Dfp(this, ln10String);
-
- }
- } else {
- // dummy settings for unused constants
- sqr2 = null;
- sqr2Split = null;
- sqr2Reciprocal = null;
- sqr3 = null;
- sqr3Reciprocal = null;
- pi = null;
- piSplit = null;
- e = null;
- eSplit = null;
- ln2 = null;
- ln2Split = null;
- ln5 = null;
- ln5Split = null;
- ln10 = null;
- }
-
- }
-
- /** Get the number of radix digits of the {@link Dfp} instances built by this factory.
- * @return number of radix digits
- */
- public int getRadixDigits() {
- return radixDigits;
- }
-
- /** Set the rounding mode.
- * If not set, the default value is {@link RoundingMode#ROUND_HALF_EVEN}.
- * @param mode desired rounding mode
- * Note that the rounding mode is common to all {@link Dfp} instances
- * belonging to the current {@link DfpField} in the system and will
- * affect all future calculations.
- */
- public void setRoundingMode(final RoundingMode mode) {
- rMode = mode;
- }
-
- /** Get the current rounding mode.
- * @return current rounding mode
- */
- public RoundingMode getRoundingMode() {
- return rMode;
- }
-
- /** Get the IEEE 854 status flags.
- * @return IEEE 854 status flags
- * @see #clearIEEEFlags()
- * @see #setIEEEFlags(int)
- * @see #setIEEEFlagsBits(int)
- * @see #FLAG_INVALID
- * @see #FLAG_DIV_ZERO
- * @see #FLAG_OVERFLOW
- * @see #FLAG_UNDERFLOW
- * @see #FLAG_INEXACT
- */
- public int getIEEEFlags() {
- return ieeeFlags;
- }
-
- /** Clears the IEEE 854 status flags.
- * @see #getIEEEFlags()
- * @see #setIEEEFlags(int)
- * @see #setIEEEFlagsBits(int)
- * @see #FLAG_INVALID
- * @see #FLAG_DIV_ZERO
- * @see #FLAG_OVERFLOW
- * @see #FLAG_UNDERFLOW
- * @see #FLAG_INEXACT
- */
- public void clearIEEEFlags() {
- ieeeFlags = 0;
- }
-
- /** Sets the IEEE 854 status flags.
- * @param flags desired value for the flags
- * @see #getIEEEFlags()
- * @see #clearIEEEFlags()
- * @see #setIEEEFlagsBits(int)
- * @see #FLAG_INVALID
- * @see #FLAG_DIV_ZERO
- * @see #FLAG_OVERFLOW
- * @see #FLAG_UNDERFLOW
- * @see #FLAG_INEXACT
- */
- public void setIEEEFlags(final int flags) {
- ieeeFlags = flags & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT);
- }
-
- /** Sets some bits in the IEEE 854 status flags, without changing the already set bits.
- * <p>
- * Calling this method is equivalent to call {@code setIEEEFlags(getIEEEFlags() | bits)}
- * </p>
- * @param bits bits to set
- * @see #getIEEEFlags()
- * @see #clearIEEEFlags()
- * @see #setIEEEFlags(int)
- * @see #FLAG_INVALID
- * @see #FLAG_DIV_ZERO
- * @see #FLAG_OVERFLOW
- * @see #FLAG_UNDERFLOW
- * @see #FLAG_INEXACT
- */
- public void setIEEEFlagsBits(final int bits) {
- ieeeFlags |= bits & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT);
- }
-
- /** Makes a {@link Dfp} with a value of 0.
- * @return a new {@link Dfp} with a value of 0
- */
- public Dfp newDfp() {
- return new Dfp(this);
- }
-
- /** Create an instance from a byte value.
- * @param x value to convert to an instance
- * @return a new {@link Dfp} with the same value as x
- */
- public Dfp newDfp(final byte x) {
- return new Dfp(this, x);
- }
-
- /** Create an instance from an int value.
- * @param x value to convert to an instance
- * @return a new {@link Dfp} with the same value as x
- */
- public Dfp newDfp(final int x) {
- return new Dfp(this, x);
- }
-
- /** Create an instance from a long value.
- * @param x value to convert to an instance
- * @return a new {@link Dfp} with the same value as x
- */
- public Dfp newDfp(final long x) {
- return new Dfp(this, x);
- }
-
- /** Create an instance from a double value.
- * @param x value to convert to an instance
- * @return a new {@link Dfp} with the same value as x
- */
- public Dfp newDfp(final double x) {
- return new Dfp(this, x);
- }
-
- /** Copy constructor.
- * @param d instance to copy
- * @return a new {@link Dfp} with the same value as d
- */
- public Dfp newDfp(Dfp d) {
- return new Dfp(d);
- }
-
- /** Create a {@link Dfp} given a String representation.
- * @param s string representation of the instance
- * @return a new {@link Dfp} parsed from specified string
- */
- public Dfp newDfp(final String s) {
- return new Dfp(this, s);
- }
-
- /** Creates a {@link Dfp} with a non-finite value.
- * @param sign sign of the Dfp to create
- * @param nans code of the value, must be one of {@link Dfp#INFINITE},
- * {@link Dfp#SNAN}, {@link Dfp#QNAN}
- * @return a new {@link Dfp} with a non-finite value
- */
- public Dfp newDfp(final byte sign, final byte nans) {
- return new Dfp(this, sign, nans);
- }
-
- /** Get the constant 0.
- * @return a {@link Dfp} with value 0
- */
- public Dfp getZero() {
- return zero;
- }
-
- /** Get the constant 1.
- * @return a {@link Dfp} with value 1
- */
- public Dfp getOne() {
- return one;
- }
-
- /** {@inheritDoc} */
- public Class<? extends FieldElement<Dfp>> getRuntimeClass() {
- return Dfp.class;
- }
-
- /** Get the constant 2.
- * @return a {@link Dfp} with value 2
- */
- public Dfp getTwo() {
- return two;
- }
-
- /** Get the constant √2.
- * @return a {@link Dfp} with value √2
- */
- public Dfp getSqr2() {
- return sqr2;
- }
-
- /** Get the constant √2 split in two pieces.
- * @return a {@link Dfp} with value √2 split in two pieces
- */
- public Dfp[] getSqr2Split() {
- return sqr2Split.clone();
- }
-
- /** Get the constant √2 / 2.
- * @return a {@link Dfp} with value √2 / 2
- */
- public Dfp getSqr2Reciprocal() {
- return sqr2Reciprocal;
- }
-
- /** Get the constant √3.
- * @return a {@link Dfp} with value √3
- */
- public Dfp getSqr3() {
- return sqr3;
- }
-
- /** Get the constant √3 / 3.
- * @return a {@link Dfp} with value √3 / 3
- */
- public Dfp getSqr3Reciprocal() {
- return sqr3Reciprocal;
- }
-
- /** Get the constant π.
- * @return a {@link Dfp} with value π
- */
- public Dfp getPi() {
- return pi;
- }
-
- /** Get the constant π split in two pieces.
- * @return a {@link Dfp} with value π split in two pieces
- */
- public Dfp[] getPiSplit() {
- return piSplit.clone();
- }
-
- /** Get the constant e.
- * @return a {@link Dfp} with value e
- */
- public Dfp getE() {
- return e;
- }
-
- /** Get the constant e split in two pieces.
- * @return a {@link Dfp} with value e split in two pieces
- */
- public Dfp[] getESplit() {
- return eSplit.clone();
- }
-
- /** Get the constant ln(2).
- * @return a {@link Dfp} with value ln(2)
- */
- public Dfp getLn2() {
- return ln2;
- }
-
- /** Get the constant ln(2) split in two pieces.
- * @return a {@link Dfp} with value ln(2) split in two pieces
- */
- public Dfp[] getLn2Split() {
- return ln2Split.clone();
- }
-
- /** Get the constant ln(5).
- * @return a {@link Dfp} with value ln(5)
- */
- public Dfp getLn5() {
- return ln5;
- }
-
- /** Get the constant ln(5) split in two pieces.
- * @return a {@link Dfp} with value ln(5) split in two pieces
- */
- public Dfp[] getLn5Split() {
- return ln5Split.clone();
- }
-
- /** Get the constant ln(10).
- * @return a {@link Dfp} with value ln(10)
- */
- public Dfp getLn10() {
- return ln10;
- }
-
- /** Breaks a string representation up into two {@link Dfp}'s.
- * The split is such that the sum of them is equivalent to the input string,
- * but has higher precision than using a single Dfp.
- * @param a string representation of the number to split
- * @return an array of two {@link Dfp Dfp} instances which sum equals a
- */
- private Dfp[] split(final String a) {
- Dfp result[] = new Dfp[2];
- boolean leading = true;
- int sp = 0;
- int sig = 0;
-
- char[] buf = new char[a.length()];
-
- for (int i = 0; i < buf.length; i++) {
- buf[i] = a.charAt(i);
-
- if (buf[i] >= '1' && buf[i] <= '9') {
- leading = false;
- }
-
- if (buf[i] == '.') {
- sig += (400 - sig) % 4;
- leading = false;
- }
-
- if (sig == (radixDigits / 2) * 4) {
- sp = i;
- break;
- }
-
- if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
- sig ++;
- }
- }
-
- result[0] = new Dfp(this, new String(buf, 0, sp));
-
- for (int i = 0; i < buf.length; i++) {
- buf[i] = a.charAt(i);
- if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
- buf[i] = '0';
- }
- }
-
- result[1] = new Dfp(this, new String(buf));
-
- return result;
-
- }
-
- /** Recompute the high precision string constants.
- * @param highPrecisionDecimalDigits precision at which the string constants mus be computed
- */
- private static void computeStringConstants(final int highPrecisionDecimalDigits) {
- if (sqr2String == null || sqr2String.length() < highPrecisionDecimalDigits - 3) {
-
- // recompute the string representation of the transcendental constants
- final DfpField highPrecisionField = new DfpField(highPrecisionDecimalDigits, false);
- final Dfp highPrecisionOne = new Dfp(highPrecisionField, 1);
- final Dfp highPrecisionTwo = new Dfp(highPrecisionField, 2);
- final Dfp highPrecisionThree = new Dfp(highPrecisionField, 3);
-
- final Dfp highPrecisionSqr2 = highPrecisionTwo.sqrt();
- sqr2String = highPrecisionSqr2.toString();
- sqr2ReciprocalString = highPrecisionOne.divide(highPrecisionSqr2).toString();
-
- final Dfp highPrecisionSqr3 = highPrecisionThree.sqrt();
- sqr3String = highPrecisionSqr3.toString();
- sqr3ReciprocalString = highPrecisionOne.divide(highPrecisionSqr3).toString();
-
- piString = computePi(highPrecisionOne, highPrecisionTwo, highPrecisionThree).toString();
- eString = computeExp(highPrecisionOne, highPrecisionOne).toString();
- ln2String = computeLn(highPrecisionTwo, highPrecisionOne, highPrecisionTwo).toString();
- ln5String = computeLn(new Dfp(highPrecisionField, 5), highPrecisionOne, highPrecisionTwo).toString();
- ln10String = computeLn(new Dfp(highPrecisionField, 10), highPrecisionOne, highPrecisionTwo).toString();
-
- }
- }
-
- /** Compute π using Jonathan and Peter Borwein quartic formula.
- * @param one constant with value 1 at desired precision
- * @param two constant with value 2 at desired precision
- * @param three constant with value 3 at desired precision
- * @return π
- */
- private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) {
-
- Dfp sqrt2 = two.sqrt();
- Dfp yk = sqrt2.subtract(one);
- Dfp four = two.add(two);
- Dfp two2kp3 = two;
- Dfp ak = two.multiply(three.subtract(two.multiply(sqrt2)));
-
- // The formula converges quartically. This means the number of correct
- // digits is multiplied by 4 at each iteration! Five iterations are
- // sufficient for about 160 digits, eight iterations give about
- // 10000 digits (this has been checked) and 20 iterations more than
- // 160 billions of digits (this has NOT been checked).
- // So the limit here is considered sufficient for most purposes ...
- for (int i = 1; i < 20; i++) {
- final Dfp ykM1 = yk;
-
- final Dfp y2 = yk.multiply(yk);
- final Dfp oneMinusY4 = one.subtract(y2.multiply(y2));
- final Dfp s = oneMinusY4.sqrt().sqrt();
- yk = one.subtract(s).divide(one.add(s));
-
- two2kp3 = two2kp3.multiply(four);
-
- final Dfp p = one.add(yk);
- final Dfp p2 = p.multiply(p);
- ak = ak.multiply(p2.multiply(p2)).subtract(two2kp3.multiply(yk).multiply(one.add(yk).add(yk.multiply(yk))));
-
- if (yk.equals(ykM1)) {
- break;
- }
- }
-
- return one.divide(ak);
-
- }
-
- /** Compute exp(a).
- * @param a number for which we want the exponential
- * @param one constant with value 1 at desired precision
- * @return exp(a)
- */
- public static Dfp computeExp(final Dfp a, final Dfp one) {
-
- Dfp y = new Dfp(one);
- Dfp py = new Dfp(one);
- Dfp f = new Dfp(one);
- Dfp fi = new Dfp(one);
- Dfp x = new Dfp(one);
-
- for (int i = 0; i < 10000; i++) {
- x = x.multiply(a);
- y = y.add(x.divide(f));
- fi = fi.add(one);
- f = f.multiply(fi);
- if (y.equals(py)) {
- break;
- }
- py = new Dfp(y);
- }
-
- return y;
-
- }
-
-
- /** Compute ln(a).
- *
- * Let f(x) = ln(x),
- *
- * We know that f'(x) = 1/x, thus from Taylor's theorem we have:
- *
- * ----- n+1 n
- * f(x) = \ (-1) (x - 1)
- * / ---------------- for 1 <= n <= infinity
- * ----- n
- *
- * or
- * 2 3 4
- * (x-1) (x-1) (x-1)
- * ln(x) = (x-1) - ----- + ------ - ------ + ...
- * 2 3 4
- *
- * alternatively,
- *
- * 2 3 4
- * x x x
- * ln(x+1) = x - - + - - - + ...
- * 2 3 4
- *
- * This series can be used to compute ln(x), but it converges too slowly.
- *
- * If we substitute -x for x above, we get
- *
- * 2 3 4
- * x x x
- * ln(1-x) = -x - - - - - - + ...
- * 2 3 4
- *
- * Note that all terms are now negative. Because the even powered ones
- * absorbed the sign. Now, subtract the series above from the previous
- * one to get ln(x+1) - ln(1-x). Note the even terms cancel out leaving
- * only the odd ones
- *
- * 3 5 7
- * 2x 2x 2x
- * ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
- * 3 5 7
- *
- * By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
- *
- * 3 5 7
- * x+1 / x x x \
- * ln ----- = 2 * | x + ---- + ---- + ---- + ... |
- * x-1 \ 3 5 7 /
- *
- * But now we want to find ln(a), so we need to find the value of x
- * such that a = (x+1)/(x-1). This is easily solved to find that
- * x = (a-1)/(a+1).
- * @param a number for which we want the exponential
- * @param one constant with value 1 at desired precision
- * @param two constant with value 2 at desired precision
- * @return ln(a)
- */
-
- public static Dfp computeLn(final Dfp a, final Dfp one, final Dfp two) {
-
- int den = 1;
- Dfp x = a.add(new Dfp(a.getField(), -1)).divide(a.add(one));
-
- Dfp y = new Dfp(x);
- Dfp num = new Dfp(x);
- Dfp py = new Dfp(y);
- for (int i = 0; i < 10000; i++) {
- num = num.multiply(x);
- num = num.multiply(x);
- den += 2;
- Dfp t = num.divide(den);
- y = y.add(t);
- if (y.equals(py)) {
- break;
- }
- py = new Dfp(y);
- }
-
- return y.multiply(two);
-
- }
-
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/dfp/DfpMath.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/dfp/DfpMath.java b/src/main/java/org/apache/commons/math3/dfp/DfpMath.java
deleted file mode 100644
index 3b19cb6..0000000
--- a/src/main/java/org/apache/commons/math3/dfp/DfpMath.java
+++ /dev/null
@@ -1,967 +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;
-
-/** Mathematical routines for use with {@link Dfp}.
- * The constants are defined in {@link DfpField}
- * @since 2.2
- */
-public class DfpMath {
-
- /** Name for traps triggered by pow. */
- private static final String POW_TRAP = "pow";
-
- /**
- * Private Constructor.
- */
- private DfpMath() {
- }
-
- /** Breaks a string representation up into two dfp's.
- * <p>The two dfp are such that the sum of them is equivalent
- * to the input string, but has higher precision than using a
- * single dfp. This is useful for improving accuracy of
- * exponentiation and critical multiplies.
- * @param field field to which the Dfp must belong
- * @param a string representation to split
- * @return an array of two {@link Dfp} which sum is a
- */
- protected static Dfp[] split(final DfpField field, final String a) {
- Dfp result[] = new Dfp[2];
- char[] buf;
- boolean leading = true;
- int sp = 0;
- int sig = 0;
-
- buf = new char[a.length()];
-
- for (int i = 0; i < buf.length; i++) {
- buf[i] = a.charAt(i);
-
- if (buf[i] >= '1' && buf[i] <= '9') {
- leading = false;
- }
-
- if (buf[i] == '.') {
- sig += (400 - sig) % 4;
- leading = false;
- }
-
- if (sig == (field.getRadixDigits() / 2) * 4) {
- sp = i;
- break;
- }
-
- if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
- sig ++;
- }
- }
-
- result[0] = field.newDfp(new String(buf, 0, sp));
-
- for (int i = 0; i < buf.length; i++) {
- buf[i] = a.charAt(i);
- if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
- buf[i] = '0';
- }
- }
-
- result[1] = field.newDfp(new String(buf));
-
- return result;
- }
-
- /** Splits a {@link Dfp} into 2 {@link Dfp}'s such that their sum is equal to the input {@link Dfp}.
- * @param a number to split
- * @return two elements array containing the split number
- */
- protected static Dfp[] split(final Dfp a) {
- final Dfp[] result = new Dfp[2];
- final Dfp shift = a.multiply(a.power10K(a.getRadixDigits() / 2));
- result[0] = a.add(shift).subtract(shift);
- result[1] = a.subtract(result[0]);
- return result;
- }
-
- /** Multiply two numbers that are split in to two pieces that are
- * meant to be added together.
- * Use binomial multiplication so ab = a0 b0 + a0 b1 + a1 b0 + a1 b1
- * Store the first term in result0, the rest in result1
- * @param a first factor of the multiplication, in split form
- * @param b second factor of the multiplication, in split form
- * @return a × b, in split form
- */
- protected static Dfp[] splitMult(final Dfp[] a, final Dfp[] b) {
- final Dfp[] result = new Dfp[2];
-
- result[1] = a[0].getZero();
- result[0] = a[0].multiply(b[0]);
-
- /* If result[0] is infinite or zero, don't compute result[1].
- * Attempting to do so may produce NaNs.
- */
-
- if (result[0].classify() == Dfp.INFINITE || result[0].equals(result[1])) {
- return result;
- }
-
- result[1] = a[0].multiply(b[1]).add(a[1].multiply(b[0])).add(a[1].multiply(b[1]));
-
- return result;
- }
-
- /** Divide two numbers that are split in to two pieces that are meant to be added together.
- * Inverse of split multiply above:
- * (a+b) / (c+d) = (a/c) + ( (bc-ad)/(c**2+cd) )
- * @param a dividend, in split form
- * @param b divisor, in split form
- * @return a / b, in split form
- */
- protected static Dfp[] splitDiv(final Dfp[] a, final Dfp[] b) {
- final Dfp[] result;
-
- result = new Dfp[2];
-
- result[0] = a[0].divide(b[0]);
- result[1] = a[1].multiply(b[0]).subtract(a[0].multiply(b[1]));
- result[1] = result[1].divide(b[0].multiply(b[0]).add(b[0].multiply(b[1])));
-
- return result;
- }
-
- /** Raise a split base to the a power.
- * @param base number to raise
- * @param a power
- * @return base<sup>a</sup>
- */
- protected static Dfp splitPow(final Dfp[] base, int a) {
- boolean invert = false;
-
- Dfp[] r = new Dfp[2];
-
- Dfp[] result = new Dfp[2];
- result[0] = base[0].getOne();
- result[1] = base[0].getZero();
-
- if (a == 0) {
- // Special case a = 0
- return result[0].add(result[1]);
- }
-
- if (a < 0) {
- // If a is less than zero
- invert = true;
- a = -a;
- }
-
- // Exponentiate by successive squaring
- do {
- r[0] = new Dfp(base[0]);
- r[1] = new Dfp(base[1]);
- int trial = 1;
-
- int prevtrial;
- while (true) {
- prevtrial = trial;
- trial *= 2;
- if (trial > a) {
- break;
- }
- r = splitMult(r, r);
- }
-
- trial = prevtrial;
-
- a -= trial;
- result = splitMult(result, r);
-
- } while (a >= 1);
-
- result[0] = result[0].add(result[1]);
-
- if (invert) {
- result[0] = base[0].getOne().divide(result[0]);
- }
-
- return result[0];
-
- }
-
- /** Raises base to the power a by successive squaring.
- * @param base number to raise
- * @param a power
- * @return base<sup>a</sup>
- */
- public static Dfp pow(Dfp base, int a)
- {
- boolean invert = false;
-
- Dfp result = base.getOne();
-
- if (a == 0) {
- // Special case
- return result;
- }
-
- if (a < 0) {
- invert = true;
- a = -a;
- }
-
- // Exponentiate by successive squaring
- do {
- Dfp r = new Dfp(base);
- Dfp prevr;
- int trial = 1;
- int prevtrial;
-
- do {
- prevr = new Dfp(r);
- prevtrial = trial;
- r = r.multiply(r);
- trial *= 2;
- } while (a>trial);
-
- r = prevr;
- trial = prevtrial;
-
- a -= trial;
- result = result.multiply(r);
-
- } while (a >= 1);
-
- if (invert) {
- result = base.getOne().divide(result);
- }
-
- return base.newInstance(result);
-
- }
-
- /** Computes e to the given power.
- * a is broken into two parts, such that a = n+m where n is an integer.
- * We use pow() to compute e<sup>n</sup> and a Taylor series to compute
- * e<sup>m</sup>. We return e*<sup>n</sup> × e<sup>m</sup>
- * @param a power at which e should be raised
- * @return e<sup>a</sup>
- */
- public static Dfp exp(final Dfp a) {
-
- final Dfp inta = a.rint();
- final Dfp fraca = a.subtract(inta);
-
- final int ia = inta.intValue();
- if (ia > 2147483646) {
- // return +Infinity
- return a.newInstance((byte)1, Dfp.INFINITE);
- }
-
- if (ia < -2147483646) {
- // return 0;
- return a.newInstance();
- }
-
- final Dfp einta = splitPow(a.getField().getESplit(), ia);
- final Dfp efraca = expInternal(fraca);
-
- return einta.multiply(efraca);
- }
-
- /** Computes e to the given power.
- * Where -1 < a < 1. Use the classic Taylor series. 1 + x**2/2! + x**3/3! + x**4/4! ...
- * @param a power at which e should be raised
- * @return e<sup>a</sup>
- */
- protected static Dfp expInternal(final Dfp a) {
- Dfp y = a.getOne();
- Dfp x = a.getOne();
- Dfp fact = a.getOne();
- Dfp py = new Dfp(y);
-
- for (int i = 1; i < 90; i++) {
- x = x.multiply(a);
- fact = fact.divide(i);
- y = y.add(x.multiply(fact));
- if (y.equals(py)) {
- break;
- }
- py = new Dfp(y);
- }
-
- return y;
- }
-
- /** Returns the natural logarithm of a.
- * a is first split into three parts such that a = (10000^h)(2^j)k.
- * ln(a) is computed by ln(a) = ln(5)*h + ln(2)*(h+j) + ln(k)
- * k is in the range 2/3 < k <4/3 and is passed on to a series expansion.
- * @param a number from which logarithm is requested
- * @return log(a)
- */
- public static Dfp log(Dfp a) {
- int lr;
- Dfp x;
- int ix;
- int p2 = 0;
-
- // Check the arguments somewhat here
- if (a.equals(a.getZero()) || a.lessThan(a.getZero()) || a.isNaN()) {
- // negative, zero or NaN
- a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
- return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte)1, Dfp.QNAN));
- }
-
- if (a.classify() == Dfp.INFINITE) {
- return a;
- }
-
- x = new Dfp(a);
- lr = x.log10K();
-
- x = x.divide(pow(a.newInstance(10000), lr)); /* This puts x in the range 0-10000 */
- ix = x.floor().intValue();
-
- while (ix > 2) {
- ix >>= 1;
- p2++;
- }
-
-
- Dfp[] spx = split(x);
- Dfp[] spy = new Dfp[2];
- spy[0] = pow(a.getTwo(), p2); // use spy[0] temporarily as a divisor
- spx[0] = spx[0].divide(spy[0]);
- spx[1] = spx[1].divide(spy[0]);
-
- spy[0] = a.newInstance("1.33333"); // Use spy[0] for comparison
- while (spx[0].add(spx[1]).greaterThan(spy[0])) {
- spx[0] = spx[0].divide(2);
- spx[1] = spx[1].divide(2);
- p2++;
- }
-
- // X is now in the range of 2/3 < x < 4/3
- Dfp[] spz = logInternal(spx);
-
- spx[0] = a.newInstance(new StringBuilder().append(p2+4*lr).toString());
- spx[1] = a.getZero();
- spy = splitMult(a.getField().getLn2Split(), spx);
-
- spz[0] = spz[0].add(spy[0]);
- spz[1] = spz[1].add(spy[1]);
-
- spx[0] = a.newInstance(new StringBuilder().append(4*lr).toString());
- spx[1] = a.getZero();
- spy = splitMult(a.getField().getLn5Split(), spx);
-
- spz[0] = spz[0].add(spy[0]);
- spz[1] = spz[1].add(spy[1]);
-
- return a.newInstance(spz[0].add(spz[1]));
-
- }
-
- /** Computes the natural log of a number between 0 and 2.
- * Let f(x) = ln(x),
- *
- * We know that f'(x) = 1/x, thus from Taylor's theorum we have:
- *
- * ----- n+1 n
- * f(x) = \ (-1) (x - 1)
- * / ---------------- for 1 <= n <= infinity
- * ----- n
- *
- * or
- * 2 3 4
- * (x-1) (x-1) (x-1)
- * ln(x) = (x-1) - ----- + ------ - ------ + ...
- * 2 3 4
- *
- * alternatively,
- *
- * 2 3 4
- * x x x
- * ln(x+1) = x - - + - - - + ...
- * 2 3 4
- *
- * This series can be used to compute ln(x), but it converges too slowly.
- *
- * If we substitute -x for x above, we get
- *
- * 2 3 4
- * x x x
- * ln(1-x) = -x - - - - - - + ...
- * 2 3 4
- *
- * Note that all terms are now negative. Because the even powered ones
- * absorbed the sign. Now, subtract the series above from the previous
- * one to get ln(x+1) - ln(1-x). Note the even terms cancel out leaving
- * only the odd ones
- *
- * 3 5 7
- * 2x 2x 2x
- * ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
- * 3 5 7
- *
- * By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
- *
- * 3 5 7
- * x+1 / x x x \
- * ln ----- = 2 * | x + ---- + ---- + ---- + ... |
- * x-1 \ 3 5 7 /
- *
- * But now we want to find ln(a), so we need to find the value of x
- * such that a = (x+1)/(x-1). This is easily solved to find that
- * x = (a-1)/(a+1).
- * @param a number from which logarithm is requested, in split form
- * @return log(a)
- */
- protected static Dfp[] logInternal(final Dfp a[]) {
-
- /* Now we want to compute x = (a-1)/(a+1) but this is prone to
- * loss of precision. So instead, compute x = (a/4 - 1/4) / (a/4 + 1/4)
- */
- Dfp t = a[0].divide(4).add(a[1].divide(4));
- Dfp x = t.add(a[0].newInstance("-0.25")).divide(t.add(a[0].newInstance("0.25")));
-
- Dfp y = new Dfp(x);
- Dfp num = new Dfp(x);
- Dfp py = new Dfp(y);
- int den = 1;
- for (int i = 0; i < 10000; i++) {
- num = num.multiply(x);
- num = num.multiply(x);
- den += 2;
- t = num.divide(den);
- y = y.add(t);
- if (y.equals(py)) {
- break;
- }
- py = new Dfp(y);
- }
-
- y = y.multiply(a[0].getTwo());
-
- return split(y);
-
- }
-
- /** Computes x to the y power.<p>
- *
- * Uses the following method:<p>
- *
- * <ol>
- * <li> Set u = rint(y), v = y-u
- * <li> Compute a = v * ln(x)
- * <li> Compute b = rint( a/ln(2) )
- * <li> Compute c = a - b*ln(2)
- * <li> x<sup>y</sup> = x<sup>u</sup> * 2<sup>b</sup> * e<sup>c</sup>
- * </ol>
- * if |y| > 1e8, then we compute by exp(y*ln(x)) <p>
- *
- * <b>Special Cases</b><p>
- * <ul>
- * <li> if y is 0.0 or -0.0 then result is 1.0
- * <li> if y is 1.0 then result is x
- * <li> if y is NaN then result is NaN
- * <li> if x is NaN and y is not zero then result is NaN
- * <li> if |x| > 1.0 and y is +Infinity then result is +Infinity
- * <li> if |x| < 1.0 and y is -Infinity then result is +Infinity
- * <li> if |x| > 1.0 and y is -Infinity then result is +0
- * <li> if |x| < 1.0 and y is +Infinity then result is +0
- * <li> if |x| = 1.0 and y is +/-Infinity then result is NaN
- * <li> if x = +0 and y > 0 then result is +0
- * <li> if x = +Inf and y < 0 then result is +0
- * <li> if x = +0 and y < 0 then result is +Inf
- * <li> if x = +Inf and y > 0 then result is +Inf
- * <li> if x = -0 and y > 0, finite, not odd integer then result is +0
- * <li> if x = -0 and y < 0, finite, and odd integer then result is -Inf
- * <li> if x = -Inf and y > 0, finite, and odd integer then result is -Inf
- * <li> if x = -0 and y < 0, not finite odd integer then result is +Inf
- * <li> if x = -Inf and y > 0, not finite odd integer then result is +Inf
- * <li> if x < 0 and y > 0, finite, and odd integer then result is -(|x|<sup>y</sup>)
- * <li> if x < 0 and y > 0, finite, and not integer then result is NaN
- * </ul>
- * @param x base to be raised
- * @param y power to which base should be raised
- * @return x<sup>y</sup>
- */
- public static Dfp pow(Dfp x, final Dfp y) {
-
- // make sure we don't mix number with different precision
- if (x.getField().getRadixDigits() != y.getField().getRadixDigits()) {
- x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
- final Dfp result = x.newInstance(x.getZero());
- result.nans = Dfp.QNAN;
- return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, result);
- }
-
- final Dfp zero = x.getZero();
- final Dfp one = x.getOne();
- final Dfp two = x.getTwo();
- boolean invert = false;
- int ui;
-
- /* Check for special cases */
- if (y.equals(zero)) {
- return x.newInstance(one);
- }
-
- if (y.equals(one)) {
- if (x.isNaN()) {
- // Test for NaNs
- x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
- return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x);
- }
- return x;
- }
-
- if (x.isNaN() || y.isNaN()) {
- // Test for NaNs
- x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
- return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
- }
-
- // X == 0
- if (x.equals(zero)) {
- if (Dfp.copysign(one, x).greaterThan(zero)) {
- // X == +0
- if (y.greaterThan(zero)) {
- return x.newInstance(zero);
- } else {
- return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
- }
- } else {
- // X == -0
- if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
- // If y is odd integer
- if (y.greaterThan(zero)) {
- return x.newInstance(zero.negate());
- } else {
- return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
- }
- } else {
- // Y is not odd integer
- if (y.greaterThan(zero)) {
- return x.newInstance(zero);
- } else {
- return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
- }
- }
- }
- }
-
- if (x.lessThan(zero)) {
- // Make x positive, but keep track of it
- x = x.negate();
- invert = true;
- }
-
- if (x.greaterThan(one) && y.classify() == Dfp.INFINITE) {
- if (y.greaterThan(zero)) {
- return y;
- } else {
- return x.newInstance(zero);
- }
- }
-
- if (x.lessThan(one) && y.classify() == Dfp.INFINITE) {
- if (y.greaterThan(zero)) {
- return x.newInstance(zero);
- } else {
- return x.newInstance(Dfp.copysign(y, one));
- }
- }
-
- if (x.equals(one) && y.classify() == Dfp.INFINITE) {
- x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
- return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
- }
-
- if (x.classify() == Dfp.INFINITE) {
- // x = +/- inf
- if (invert) {
- // negative infinity
- if (y.classify() == Dfp.FINITE && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
- // If y is odd integer
- if (y.greaterThan(zero)) {
- return x.newInstance(x.newInstance((byte)-1, Dfp.INFINITE));
- } else {
- return x.newInstance(zero.negate());
- }
- } else {
- // Y is not odd integer
- if (y.greaterThan(zero)) {
- return x.newInstance(x.newInstance((byte)1, Dfp.INFINITE));
- } else {
- return x.newInstance(zero);
- }
- }
- } else {
- // positive infinity
- if (y.greaterThan(zero)) {
- return x;
- } else {
- return x.newInstance(zero);
- }
- }
- }
-
- if (invert && !y.rint().equals(y)) {
- x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
- return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, Dfp.QNAN));
- }
-
- // End special cases
-
- Dfp r;
- if (y.lessThan(x.newInstance(100000000)) && y.greaterThan(x.newInstance(-100000000))) {
- final Dfp u = y.rint();
- ui = u.intValue();
-
- final Dfp v = y.subtract(u);
-
- if (v.unequal(zero)) {
- final Dfp a = v.multiply(log(x));
- final Dfp b = a.divide(x.getField().getLn2()).rint();
-
- final Dfp c = a.subtract(b.multiply(x.getField().getLn2()));
- r = splitPow(split(x), ui);
- r = r.multiply(pow(two, b.intValue()));
- r = r.multiply(exp(c));
- } else {
- r = splitPow(split(x), ui);
- }
- } else {
- // very large exponent. |y| > 1e8
- r = exp(log(x).multiply(y));
- }
-
- if (invert && y.rint().equals(y) && !y.remainder(two).equals(zero)) {
- // if y is odd integer
- r = r.negate();
- }
-
- return x.newInstance(r);
-
- }
-
- /** Computes sin(a) Used when 0 < a < pi/4.
- * Uses the classic Taylor series. x - x**3/3! + x**5/5! ...
- * @param a number from which sine is desired, in split form
- * @return sin(a)
- */
- protected static Dfp sinInternal(Dfp a[]) {
-
- Dfp c = a[0].add(a[1]);
- Dfp y = c;
- c = c.multiply(c);
- Dfp x = y;
- Dfp fact = a[0].getOne();
- Dfp py = new Dfp(y);
-
- for (int i = 3; i < 90; i += 2) {
- x = x.multiply(c);
- x = x.negate();
-
- fact = fact.divide((i-1)*i); // 1 over fact
- y = y.add(x.multiply(fact));
- if (y.equals(py)) {
- break;
- }
- py = new Dfp(y);
- }
-
- return y;
-
- }
-
- /** Computes cos(a) Used when 0 < a < pi/4.
- * Uses the classic Taylor series for cosine. 1 - x**2/2! + x**4/4! ...
- * @param a number from which cosine is desired, in split form
- * @return cos(a)
- */
- protected static Dfp cosInternal(Dfp a[]) {
- final Dfp one = a[0].getOne();
-
-
- Dfp x = one;
- Dfp y = one;
- Dfp c = a[0].add(a[1]);
- c = c.multiply(c);
-
- Dfp fact = one;
- Dfp py = new Dfp(y);
-
- for (int i = 2; i < 90; i += 2) {
- x = x.multiply(c);
- x = x.negate();
-
- fact = fact.divide((i - 1) * i); // 1 over fact
-
- y = y.add(x.multiply(fact));
- if (y.equals(py)) {
- break;
- }
- py = new Dfp(y);
- }
-
- return y;
-
- }
-
- /** computes the sine of the argument.
- * @param a number from which sine is desired
- * @return sin(a)
- */
- public static Dfp sin(final Dfp a) {
- final Dfp pi = a.getField().getPi();
- final Dfp zero = a.getField().getZero();
- boolean neg = false;
-
- /* First reduce the argument to the range of +/- PI */
- Dfp x = a.remainder(pi.multiply(2));
-
- /* if x < 0 then apply identity sin(-x) = -sin(x) */
- /* This puts x in the range 0 < x < PI */
- if (x.lessThan(zero)) {
- x = x.negate();
- neg = true;
- }
-
- /* Since sine(x) = sine(pi - x) we can reduce the range to
- * 0 < x < pi/2
- */
-
- if (x.greaterThan(pi.divide(2))) {
- x = pi.subtract(x);
- }
-
- Dfp y;
- if (x.lessThan(pi.divide(4))) {
- Dfp c[] = new Dfp[2];
- c[0] = x;
- c[1] = zero;
-
- //y = sinInternal(c);
- y = sinInternal(split(x));
- } else {
- final Dfp c[] = new Dfp[2];
- final Dfp[] piSplit = a.getField().getPiSplit();
- c[0] = piSplit[0].divide(2).subtract(x);
- c[1] = piSplit[1].divide(2);
- y = cosInternal(c);
- }
-
- if (neg) {
- y = y.negate();
- }
-
- return a.newInstance(y);
-
- }
-
- /** computes the cosine of the argument.
- * @param a number from which cosine is desired
- * @return cos(a)
- */
- public static Dfp cos(Dfp a) {
- final Dfp pi = a.getField().getPi();
- final Dfp zero = a.getField().getZero();
- boolean neg = false;
-
- /* First reduce the argument to the range of +/- PI */
- Dfp x = a.remainder(pi.multiply(2));
-
- /* if x < 0 then apply identity cos(-x) = cos(x) */
- /* This puts x in the range 0 < x < PI */
- if (x.lessThan(zero)) {
- x = x.negate();
- }
-
- /* Since cos(x) = -cos(pi - x) we can reduce the range to
- * 0 < x < pi/2
- */
-
- if (x.greaterThan(pi.divide(2))) {
- x = pi.subtract(x);
- neg = true;
- }
-
- Dfp y;
- if (x.lessThan(pi.divide(4))) {
- Dfp c[] = new Dfp[2];
- c[0] = x;
- c[1] = zero;
-
- y = cosInternal(c);
- } else {
- final Dfp c[] = new Dfp[2];
- final Dfp[] piSplit = a.getField().getPiSplit();
- c[0] = piSplit[0].divide(2).subtract(x);
- c[1] = piSplit[1].divide(2);
- y = sinInternal(c);
- }
-
- if (neg) {
- y = y.negate();
- }
-
- return a.newInstance(y);
-
- }
-
- /** computes the tangent of the argument.
- * @param a number from which tangent is desired
- * @return tan(a)
- */
- public static Dfp tan(final Dfp a) {
- return sin(a).divide(cos(a));
- }
-
- /** computes the arc-tangent of the argument.
- * @param a number from which arc-tangent is desired
- * @return atan(a)
- */
- protected static Dfp atanInternal(final Dfp a) {
-
- Dfp y = new Dfp(a);
- Dfp x = new Dfp(y);
- Dfp py = new Dfp(y);
-
- for (int i = 3; i < 90; i += 2) {
- x = x.multiply(a);
- x = x.multiply(a);
- x = x.negate();
- y = y.add(x.divide(i));
- if (y.equals(py)) {
- break;
- }
- py = new Dfp(y);
- }
-
- return y;
-
- }
-
- /** computes the arc tangent of the argument
- *
- * Uses the typical taylor series
- *
- * but may reduce arguments using the following identity
- * tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
- *
- * since tan(PI/8) = sqrt(2)-1,
- *
- * atan(x) = atan( (x - sqrt(2) + 1) / (1+x*sqrt(2) - x) + PI/8.0
- * @param a number from which arc-tangent is desired
- * @return atan(a)
- */
- public static Dfp atan(final Dfp a) {
- final Dfp zero = a.getField().getZero();
- final Dfp one = a.getField().getOne();
- final Dfp[] sqr2Split = a.getField().getSqr2Split();
- final Dfp[] piSplit = a.getField().getPiSplit();
- boolean recp = false;
- boolean neg = false;
- boolean sub = false;
-
- final Dfp ty = sqr2Split[0].subtract(one).add(sqr2Split[1]);
-
- Dfp x = new Dfp(a);
- if (x.lessThan(zero)) {
- neg = true;
- x = x.negate();
- }
-
- if (x.greaterThan(one)) {
- recp = true;
- x = one.divide(x);
- }
-
- if (x.greaterThan(ty)) {
- Dfp sty[] = new Dfp[2];
- sub = true;
-
- sty[0] = sqr2Split[0].subtract(one);
- sty[1] = sqr2Split[1];
-
- Dfp[] xs = split(x);
-
- Dfp[] ds = splitMult(xs, sty);
- ds[0] = ds[0].add(one);
-
- xs[0] = xs[0].subtract(sty[0]);
- xs[1] = xs[1].subtract(sty[1]);
-
- xs = splitDiv(xs, ds);
- x = xs[0].add(xs[1]);
-
- //x = x.subtract(ty).divide(dfp.one.add(x.multiply(ty)));
- }
-
- Dfp y = atanInternal(x);
-
- if (sub) {
- y = y.add(piSplit[0].divide(8)).add(piSplit[1].divide(8));
- }
-
- if (recp) {
- y = piSplit[0].divide(2).subtract(y).add(piSplit[1].divide(2));
- }
-
- if (neg) {
- y = y.negate();
- }
-
- return a.newInstance(y);
-
- }
-
- /** computes the arc-sine of the argument.
- * @param a number from which arc-sine is desired
- * @return asin(a)
- */
- public static Dfp asin(final Dfp a) {
- return atan(a.divide(a.getOne().subtract(a.multiply(a)).sqrt()));
- }
-
- /** computes the arc-cosine of the argument.
- * @param a number from which arc-cosine is desired
- * @return acos(a)
- */
- public static Dfp acos(Dfp a) {
- Dfp result;
- boolean negative = false;
-
- if (a.lessThan(a.getZero())) {
- negative = true;
- }
-
- a = Dfp.copysign(a, a.getOne()); // absolute value
-
- result = atan(a.getOne().subtract(a.multiply(a)).sqrt().divide(a));
-
- if (negative) {
- result = a.getField().getPi().subtract(result);
- }
-
- return a.newInstance(result);
- }
-
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/dfp/UnivariateDfpFunction.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/dfp/UnivariateDfpFunction.java b/src/main/java/org/apache/commons/math3/dfp/UnivariateDfpFunction.java
deleted file mode 100644
index b627a32..0000000
--- a/src/main/java/org/apache/commons/math3/dfp/UnivariateDfpFunction.java
+++ /dev/null
@@ -1,39 +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;
-
-/**
- * An interface representing a univariate {@link Dfp} function.
- *
- */
-public interface UnivariateDfpFunction {
-
- /**
- * Compute the value of the function.
- *
- * @param x Point at which the function value should be computed.
- * @return the value.
- * @throws IllegalArgumentException when the activated method itself can
- * ascertain that preconditions, specified in the API expressed at the
- * level of the activated method, have been violated. In the vast
- * majority of cases where Commons-Math throws IllegalArgumentException,
- * it is the result of argument checking of actual parameters immediately
- * passed to a method.
- */
- Dfp value(Dfp x);
-
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/dfp/package-info.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/dfp/package-info.java b/src/main/java/org/apache/commons/math3/dfp/package-info.java
deleted file mode 100644
index 42a4b48..0000000
--- a/src/main/java/org/apache/commons/math3/dfp/package-info.java
+++ /dev/null
@@ -1,88 +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.
- */
-/**
- *
- * 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 keep 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 add 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>
- *
- */
-package org.apache.commons.math3.dfp;
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/AbstractIntegerDistribution.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/distribution/AbstractIntegerDistribution.java b/src/main/java/org/apache/commons/math3/distribution/AbstractIntegerDistribution.java
deleted file mode 100644
index 82a96c5..0000000
--- a/src/main/java/org/apache/commons/math3/distribution/AbstractIntegerDistribution.java
+++ /dev/null
@@ -1,253 +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.distribution;
-
-import java.io.Serializable;
-
-import org.apache.commons.math3.exception.MathInternalError;
-import org.apache.commons.math3.exception.NotStrictlyPositiveException;
-import org.apache.commons.math3.exception.NumberIsTooLargeException;
-import org.apache.commons.math3.exception.OutOfRangeException;
-import org.apache.commons.math3.exception.util.LocalizedFormats;
-import org.apache.commons.math3.random.RandomGenerator;
-import org.apache.commons.math3.random.RandomDataImpl;
-import org.apache.commons.math3.util.FastMath;
-
-/**
- * Base class for integer-valued discrete distributions. Default
- * implementations are provided for some of the methods that do not vary
- * from distribution to distribution.
- *
- */
-public abstract class AbstractIntegerDistribution implements IntegerDistribution, Serializable {
-
- /** Serializable version identifier */
- private static final long serialVersionUID = -1146319659338487221L;
-
- /**
- * RandomData instance used to generate samples from the distribution.
- * @deprecated As of 3.1, to be removed in 4.0. Please use the
- * {@link #random} instance variable instead.
- */
- @Deprecated
- protected final RandomDataImpl randomData = new RandomDataImpl();
-
- /**
- * RNG instance used to generate samples from the distribution.
- * @since 3.1
- */
- protected final RandomGenerator random;
-
- /**
- * @deprecated As of 3.1, to be removed in 4.0. Please use
- * {@link #AbstractIntegerDistribution(RandomGenerator)} instead.
- */
- @Deprecated
- protected AbstractIntegerDistribution() {
- // Legacy users are only allowed to access the deprecated "randomData".
- // New users are forbidden to use this constructor.
- random = null;
- }
-
- /**
- * @param rng Random number generator.
- * @since 3.1
- */
- protected AbstractIntegerDistribution(RandomGenerator rng) {
- random = rng;
- }
-
- /**
- * {@inheritDoc}
- *
- * The default implementation uses the identity
- * <p>{@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}</p>
- */
- public double cumulativeProbability(int x0, int x1) throws NumberIsTooLargeException {
- if (x1 < x0) {
- throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
- x0, x1, true);
- }
- return cumulativeProbability(x1) - cumulativeProbability(x0);
- }
-
- /**
- * {@inheritDoc}
- *
- * The default implementation returns
- * <ul>
- * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
- * <li>{@link #getSupportUpperBound()} for {@code p = 1}, and</li>
- * <li>{@link #solveInverseCumulativeProbability(double, int, int)} for
- * {@code 0 < p < 1}.</li>
- * </ul>
- */
- public int inverseCumulativeProbability(final double p) throws OutOfRangeException {
- if (p < 0.0 || p > 1.0) {
- throw new OutOfRangeException(p, 0, 1);
- }
-
- int lower = getSupportLowerBound();
- if (p == 0.0) {
- return lower;
- }
- if (lower == Integer.MIN_VALUE) {
- if (checkedCumulativeProbability(lower) >= p) {
- return lower;
- }
- } else {
- lower -= 1; // this ensures cumulativeProbability(lower) < p, which
- // is important for the solving step
- }
-
- int upper = getSupportUpperBound();
- if (p == 1.0) {
- return upper;
- }
-
- // use the one-sided Chebyshev inequality to narrow the bracket
- // cf. AbstractRealDistribution.inverseCumulativeProbability(double)
- final double mu = getNumericalMean();
- final double sigma = FastMath.sqrt(getNumericalVariance());
- final boolean chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
- Double.isInfinite(sigma) || Double.isNaN(sigma) || sigma == 0.0);
- if (chebyshevApplies) {
- double k = FastMath.sqrt((1.0 - p) / p);
- double tmp = mu - k * sigma;
- if (tmp > lower) {
- lower = ((int) FastMath.ceil(tmp)) - 1;
- }
- k = 1.0 / k;
- tmp = mu + k * sigma;
- if (tmp < upper) {
- upper = ((int) FastMath.ceil(tmp)) - 1;
- }
- }
-
- return solveInverseCumulativeProbability(p, lower, upper);
- }
-
- /**
- * This is a utility function used by {@link
- * #inverseCumulativeProbability(double)}. It assumes {@code 0 < p < 1} and
- * that the inverse cumulative probability lies in the bracket {@code
- * (lower, upper]}. The implementation does simple bisection to find the
- * smallest {@code p}-quantile <code>inf{x in Z | P(X<=x) >= p}</code>.
- *
- * @param p the cumulative probability
- * @param lower a value satisfying {@code cumulativeProbability(lower) < p}
- * @param upper a value satisfying {@code p <= cumulativeProbability(upper)}
- * @return the smallest {@code p}-quantile of this distribution
- */
- protected int solveInverseCumulativeProbability(final double p, int lower, int upper) {
- while (lower + 1 < upper) {
- int xm = (lower + upper) / 2;
- if (xm < lower || xm > upper) {
- /*
- * Overflow.
- * There will never be an overflow in both calculation methods
- * for xm at the same time
- */
- xm = lower + (upper - lower) / 2;
- }
-
- double pm = checkedCumulativeProbability(xm);
- if (pm >= p) {
- upper = xm;
- } else {
- lower = xm;
- }
- }
- return upper;
- }
-
- /** {@inheritDoc} */
- public void reseedRandomGenerator(long seed) {
- random.setSeed(seed);
- randomData.reSeed(seed);
- }
-
- /**
- * {@inheritDoc}
- *
- * The default implementation uses the
- * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
- * inversion method</a>.
- */
- public int sample() {
- return inverseCumulativeProbability(random.nextDouble());
- }
-
- /**
- * {@inheritDoc}
- *
- * The default implementation generates the sample by calling
- * {@link #sample()} in a loop.
- */
- public int[] sample(int sampleSize) {
- if (sampleSize <= 0) {
- throw new NotStrictlyPositiveException(
- LocalizedFormats.NUMBER_OF_SAMPLES, sampleSize);
- }
- int[] out = new int[sampleSize];
- for (int i = 0; i < sampleSize; i++) {
- out[i] = sample();
- }
- return out;
- }
-
- /**
- * Computes the cumulative probability function and checks for {@code NaN}
- * values returned. Throws {@code MathInternalError} if the value is
- * {@code NaN}. Rethrows any exception encountered evaluating the cumulative
- * probability function. Throws {@code MathInternalError} if the cumulative
- * probability function returns {@code NaN}.
- *
- * @param argument input value
- * @return the cumulative probability
- * @throws MathInternalError if the cumulative probability is {@code NaN}
- */
- private double checkedCumulativeProbability(int argument)
- throws MathInternalError {
- double result = Double.NaN;
- result = cumulativeProbability(argument);
- if (Double.isNaN(result)) {
- throw new MathInternalError(LocalizedFormats
- .DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN, argument);
- }
- return result;
- }
-
- /**
- * For a random variable {@code X} whose values are distributed according to
- * this distribution, this method returns {@code log(P(X = x))}, where
- * {@code log} is the natural logarithm. In other words, this method
- * represents the logarithm of the probability mass function (PMF) for the
- * distribution. Note that due to the floating point precision and
- * under/overflow issues, this method will for some distributions be more
- * precise and faster than computing the logarithm of
- * {@link #probability(int)}.
- * <p>
- * The default implementation simply computes the logarithm of {@code probability(x)}.</p>
- *
- * @param x the point at which the PMF is evaluated
- * @return the logarithm of the value of the probability mass function at {@code x}
- */
- public double logProbability(int x) {
- return FastMath.log(probability(x));
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/AbstractMultivariateRealDistribution.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/distribution/AbstractMultivariateRealDistribution.java b/src/main/java/org/apache/commons/math3/distribution/AbstractMultivariateRealDistribution.java
deleted file mode 100644
index a1dfd64..0000000
--- a/src/main/java/org/apache/commons/math3/distribution/AbstractMultivariateRealDistribution.java
+++ /dev/null
@@ -1,70 +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.distribution;
-
-import org.apache.commons.math3.exception.NotStrictlyPositiveException;
-import org.apache.commons.math3.exception.util.LocalizedFormats;
-import org.apache.commons.math3.random.RandomGenerator;
-
-/**
- * Base class for multivariate probability distributions.
- *
- * @since 3.1
- */
-public abstract class AbstractMultivariateRealDistribution
- implements MultivariateRealDistribution {
- /** RNG instance used to generate samples from the distribution. */
- protected final RandomGenerator random;
- /** The number of dimensions or columns in the multivariate distribution. */
- private final int dimension;
-
- /**
- * @param rng Random number generator.
- * @param n Number of dimensions.
- */
- protected AbstractMultivariateRealDistribution(RandomGenerator rng,
- int n) {
- random = rng;
- dimension = n;
- }
-
- /** {@inheritDoc} */
- public void reseedRandomGenerator(long seed) {
- random.setSeed(seed);
- }
-
- /** {@inheritDoc} */
- public int getDimension() {
- return dimension;
- }
-
- /** {@inheritDoc} */
- public abstract double[] sample();
-
- /** {@inheritDoc} */
- public double[][] sample(final int sampleSize) {
- if (sampleSize <= 0) {
- throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES,
- sampleSize);
- }
- final double[][] out = new double[sampleSize][dimension];
- for (int i = 0; i < sampleSize; i++) {
- out[i] = sample();
- }
- return out;
- }
-}