commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r992696 [2/4] - in /commons/proper/math/branches/MATH_2_X/src: main/java/org/apache/commons/math/dfp/ site/xdoc/ test/java/org/apache/commons/math/dfp/
Date Sat, 04 Sep 2010 22:57:31 GMT
Added: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpDec.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpDec.java?rev=992696&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpDec.java (added)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpDec.java Sat Sep  4 22:57:31 2010
@@ -0,0 +1,359 @@
+/*
+ * 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.math.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.
+ * @version $Revision$ $Date$
+ * @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} */
+    public Dfp newInstance() {
+        return new DfpDec(getField());
+    }
+
+    /** {@inheritDoc} */
+    public Dfp newInstance(final byte x) {
+        return new DfpDec(getField(), x);
+    }
+
+    /** {@inheritDoc} */
+    public Dfp newInstance(final int x) {
+        return new DfpDec(getField(), x);
+    }
+
+    /** {@inheritDoc} */
+    public Dfp newInstance(final long x) {
+        return new DfpDec(getField(), x);
+    }
+
+    /** {@inheritDoc} */
+    public Dfp newInstance(final double x) {
+        return new DfpDec(getField(), x);
+    }
+
+    /** {@inheritDoc} */
+    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} */
+    public Dfp newInstance(final String s) {
+        return new DfpDec(getField(), s);
+    }
+
+    /** {@inheritDoc} */
+    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} */
+    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} */
+    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(log10() - 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(log10());
+            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;
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpDec.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpDec.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java?rev=992696&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java (added)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java Sat Sep  4 22:57:31 2010
@@ -0,0 +1,742 @@
+/*
+ * 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.math.dfp;
+
+import org.apache.commons.math.Field;
+
+/** Field for Decimal floating point instances.
+ * @version $Revision$ $Date$
+ * @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 &radic;2. */
+    private static String sqr2String;
+
+    /** High precision string representation of &radic;2 / 2. */
+    private static String sqr2ReciprocalString;
+
+    /** High precision string representation of &radic;3. */
+    private static String sqr3String;
+
+    /** High precision string representation of &radic;3 / 3. */
+    private static String sqr3ReciprocalString;
+
+    /** High precision string representation of &pi;. */
+    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 &radic;2. */
+    private final Dfp sqr2;
+
+    /** A two elements {@link Dfp} array with value &radic;2 split in two pieces. */
+    private final Dfp[] sqr2Split;
+
+    /** A {@link Dfp} with value &radic;2 / 2. */
+    private final Dfp sqr2Reciprocal;
+
+    /** A {@link Dfp} with value &radic;3. */
+    private final Dfp sqr3;
+
+    /** A {@link Dfp} with value &radic;3 / 3. */
+    private final Dfp sqr3Reciprocal;
+
+    /** A {@link Dfp} with value &pi;. */
+    private final Dfp pi;
+
+    /** A two elements {@link Dfp} array with value &pi; 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)
+     */
+    public DfpField(final int decimalDigits, final boolean computeConstants) {
+
+        this.radixDigits = (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;
+    }
+
+    /** Get the constant 2.
+     * @return a {@link Dfp} with value 2
+     */
+    public Dfp getTwo() {
+        return two;
+    }
+
+    /** Get the constant &radic;2.
+     * @return a {@link Dfp} with value &radic;2
+     */
+    public Dfp getSqr2() {
+        return sqr2;
+    }
+
+    /** Get the constant &radic;2 split in two pieces.
+     * @return a {@link Dfp} with value &radic;2 split in two pieces
+     */
+    public Dfp[] getSqr2Split() {
+        return sqr2Split.clone();
+    }
+
+    /** Get the constant &radic;2 / 2.
+     * @return a {@link Dfp} with value &radic;2 / 2
+     */
+    public Dfp getSqr2Reciprocal() {
+        return sqr2Reciprocal;
+    }
+
+    /** Get the constant &radic;3.
+     * @return a {@link Dfp} with value &radic;3
+     */
+    public Dfp getSqr3() {
+        return sqr3;
+    }
+
+    /** Get the constant &radic;3 / 3.
+     * @return a {@link Dfp} with value &radic;3 / 3
+     */
+    public Dfp getSqr3Reciprocal() {
+        return sqr3Reciprocal;
+    }
+
+    /** Get the constant &pi;.
+     * @return a {@link Dfp} with value &pi;
+     */
+    public Dfp getPi() {
+        return pi;
+    }
+
+    /** Get the constant &pi; split in two pieces.
+     * @return a {@link Dfp} with value &pi; 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 &pi; by atan(1/&radic;(3)) = &pi;/6.
+     * @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 &pi;
+     */
+    private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) {
+
+        Dfp x = three;
+        x = x.sqrt();
+        x = one.divide(x);
+
+        Dfp denom = one;
+
+        Dfp py = new Dfp(x);
+        Dfp y  = new Dfp(x);
+
+        for (int i = 1; i < 10000; i++) {
+            x = x.divide(three);
+            denom = denom.add(two);
+            if ((i&1) != 0) {
+                y = y.subtract(x.divide(denom));
+            } else {
+                y = y.add(x.divide(denom));
+            }
+            if (y.equals(py)) {
+                break;
+            }
+            py = new Dfp(y);
+        }
+
+        return y.multiply(new Dfp(one.getField(), 6));
+
+    }
+
+    /** 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 = den + 2;
+            Dfp t = num.divide(den);
+            y = y.add(t);
+            if (y.equals(py)) {
+                break;
+            }
+            py = new Dfp(y);
+        }
+
+        return y.multiply(two);
+
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpField.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpMath.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpMath.java?rev=992696&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpMath.java (added)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpMath.java Sat Sep  4 22:57:31 2010
@@ -0,0 +1,969 @@
+/*
+ * 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.math.dfp;
+
+/** Mathematical routines for use with {@link Dfp}.
+ * The constants are defined in {@link DfpField}
+ * @version $Revision$ $Date$
+ * @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 &times; 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 = 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 = trial * 2;
+            } while (a>trial);
+
+            r = prevr;
+            trial = prevtrial;
+
+            a = 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> &times; 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, (byte) 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.equals(a) == false)) {
+            // negative, zero or NaN
+            a.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
+            return a.dotrap(DfpField.FLAG_INVALID, "ln", a, a.newInstance((byte)1, (byte) 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 StringBuffer().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 StringBuffer().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 = 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.equals(x)) {
+                // Test for NaNs
+                x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
+                return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x);
+            }
+            return x;
+        }
+
+        if (!x.equals(x) || !y.equals(y)) {
+            // Test for NaNs
+            x.getField().setIEEEFlagsBits(DfpField.FLAG_INVALID);
+            return x.dotrap(DfpField.FLAG_INVALID, POW_TRAP, x, x.newInstance((byte)1, (byte) 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, (byte)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, (byte)Dfp.INFINITE));
+                    }
+                } else {
+                    // Y is not odd integer
+                    if (y.greaterThan(zero)) {
+                        return x.newInstance(zero);
+                    } else {
+                        return x.newInstance(x.newInstance((byte)1, (byte)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, (byte) 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, (byte)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, (byte)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, (byte) 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) {
+            // if y is odd integer
+            if (y.rint().equals(y) && !y.remainder(two).equals(zero)) {
+                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);
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpMath.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/DfpMath.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Added: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/package.html
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/package.html?rev=992696&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/package.html (added)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/package.html Sat Sep  4 22:57:31 2010
@@ -0,0 +1,88 @@
+<html>
+<!--
+   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.
+  -->
+    <!-- $Revision$ $Date$ -->
+    <body>
+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 &times; mant &times; (radix)<sup>exp</sup>;</p>
+</pre>
+where sign is &plusmn;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>
+    </body>
+</html>

Propchange: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/package.html
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/dfp/package.html
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Modified: commons/proper/math/branches/MATH_2_X/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/site/xdoc/changes.xml?rev=992696&r1=992695&r2=992696&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/branches/MATH_2_X/src/site/xdoc/changes.xml Sat Sep  4 22:57:31 2010
@@ -52,6 +52,13 @@ The <action> type attribute can be add,u
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="2.2" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-412" due-to="Bill Rossi">
+        Added the dfp library providing arbitrary precision floating point computation in the spirit of
+        IEEE 854-1987 (not exactly as it uses base 1000 instead of base 10). In addition to finite numbers,
+        infinities and NaNs are available (but there are no subnormals). All IEEE 854-1987 rounding modes and
+        signaling flags are supported. The available operations are +, -, *, / and the available functions
+        are sqrt, sin, cos, tan, asin, acos, atan, exp, log.
+      </action>
       <action dev="luc" type="fix" issue="MATH-375" due-to="Bill Rossi">
         Added faster and more accurate version of traditional mathematical functions in a FastMath
         class intended to be a drop-in replacement for java.util.Math at source-level. Some functions

Added: commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/dfp/Decimal10.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/dfp/Decimal10.java?rev=992696&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/dfp/Decimal10.java (added)
+++ commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/dfp/Decimal10.java Sat Sep  4 22:57:31 2010
@@ -0,0 +1,90 @@
+/*
+ * 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.math.dfp;
+
+public class Decimal10 extends DfpDec {
+
+    Decimal10(final DfpField factory) {
+        super(factory);
+    }
+
+    Decimal10(final DfpField factory, final byte x) {
+        super(factory, x);
+    }
+
+    Decimal10(final DfpField factory, final int x) {
+        super(factory, x);
+    }
+
+    Decimal10(final DfpField factory, final long x) {
+        super(factory, x);
+    }
+
+    Decimal10(final DfpField factory, final double x) {
+        super(factory, x);
+    }
+
+    public Decimal10(final Dfp d) {
+        super(d);
+    }
+
+    public Decimal10(final DfpField factory, final String s) {
+        super(factory, s);
+    }
+
+    protected Decimal10(final DfpField factory, final byte sign, final byte nans) {
+        super(factory, sign, nans);
+    }
+
+    public Dfp newInstance() {
+        return new Decimal10(getField());
+    }
+
+    public Dfp newInstance(final byte x) {
+        return new Decimal10(getField(), x);
+    }
+
+    public Dfp newInstance(final int x) {
+        return new Decimal10(getField(), x);
+    }
+
+    public Dfp newInstance(final long x) {
+        return new Decimal10(getField(), x);
+    }
+
+    public Dfp newInstance(final double x) {
+        return new Decimal10(getField(), x);
+    }
+
+    public Dfp newInstance(final Dfp d) {
+        return new Decimal10(d);
+    }
+
+    public Dfp newInstance(final String s) {
+        return new Decimal10(getField(), s);
+    }
+
+    public Dfp newInstance(final byte sign, final byte nans) {
+        return new Decimal10(getField(), sign, nans);
+    }
+
+    protected int getDecimalDigits() {
+        return 10;
+    }
+
+}

Propchange: commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/dfp/Decimal10.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/branches/MATH_2_X/src/test/java/org/apache/commons/math/dfp/Decimal10.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision



Mime
View raw message