| /* |
| * 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: 995987 $ $Date: 2010-09-10 23:24:15 +0200 (ven. 10 sept. 2010) $ |
| * @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; |
| |
| /** 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; |
| } |
| |
| /** 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 = den + 2; |
| Dfp t = num.divide(den); |
| y = y.add(t); |
| if (y.equals(py)) { |
| break; |
| } |
| py = new Dfp(y); |
| } |
| |
| return y.multiply(two); |
| |
| } |
| |
| } |
