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; - } -}
