DfpField.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. package org.apache.commons.math4.legacy.core.dfp;

  20. import org.apache.commons.math4.legacy.core.Field;
  21. import org.apache.commons.math4.legacy.core.FieldElement;

  23. /** Field for Decimal floating point instances.
  24.  * @since 2.2
  25.  */
  26. public class DfpField implements Field<Dfp> {

  28.     /** Enumerate for rounding modes. */
  29.     public enum RoundingMode {

  31.         /** Rounds toward zero (truncation). */
  32.         ROUND_DOWN,

  34.         /** Rounds away from zero if discarded digit is non-zero. */
  35.         ROUND_UP,

  37.         /** Rounds towards nearest unless both are equidistant in which case it rounds away from zero. */
  38.         ROUND_HALF_UP,

  40.         /** Rounds towards nearest unless both are equidistant in which case it rounds toward zero. */
  41.         ROUND_HALF_DOWN,

  43.         /** Rounds towards nearest unless both are equidistant in which case it rounds toward the even neighbor.
  44.          * This is the default as  specified by IEEE 854-1987
  45.          */
  46.         ROUND_HALF_EVEN,

  48.         /** Rounds towards nearest unless both are equidistant in which case it rounds toward the odd neighbor.  */
  49.         ROUND_HALF_ODD,

  51.         /** Rounds towards positive infinity. */
  52.         ROUND_CEIL,

  54.         /** Rounds towards negative infinity. */
  55.         ROUND_FLOOR;
  56.     }

  58.     /** IEEE 854-1987 flag for invalid operation. */
  59.     public static final int FLAG_INVALID   =  1;

  61.     /** IEEE 854-1987 flag for division by zero. */
  62.     public static final int FLAG_DIV_ZERO  =  2;

  64.     /** IEEE 854-1987 flag for overflow. */
  65.     public static final int FLAG_OVERFLOW  =  4;

  67.     /** IEEE 854-1987 flag for underflow. */
  68.     public static final int FLAG_UNDERFLOW =  8;

  70.     /** IEEE 854-1987 flag for inexact result. */
  71.     public static final int FLAG_INEXACT   = 16;

  73.     /** High precision string representation of &radic;2. */
  74.     private static String sqr2String;

  76.     // Note: the static strings are set up (once) by the ctor and @GuardedBy("DfpField.class")

  78.     /** High precision string representation of &radic;2 / 2. */
  79.     private static String sqr2ReciprocalString;

  81.     /** High precision string representation of &radic;3. */
  82.     private static String sqr3String;

  84.     /** High precision string representation of &radic;3 / 3. */
  85.     private static String sqr3ReciprocalString;

  87.     /** High precision string representation of &pi;. */
  88.     private static String piString;

  90.     /** High precision string representation of e. */
  91.     private static String eString;

  93.     /** High precision string representation of ln(2). */
  94.     private static String ln2String;

  96.     /** High precision string representation of ln(5). */
  97.     private static String ln5String;

  99.     /** High precision string representation of ln(10). */
  100.     private static String ln10String;

  102.     /** The number of radix digits.
  103.      * Note these depend on the radix which is 10000 digits,
  104.      * so each one is equivalent to 4 decimal digits.
  105.      */
  106.     private final int radixDigits;

  108.     /** A {@link Dfp} with value 0. */
  109.     private final Dfp zero;

  111.     /** A {@link Dfp} with value 1. */
  112.     private final Dfp one;

  114.     /** A {@link Dfp} with value 2. */
  115.     private final Dfp two;

  117.     /** A {@link Dfp} with value &radic;2. */
  118.     private final Dfp sqr2;

  120.     /** A two elements {@link Dfp} array with value &radic;2 split in two pieces. */
  121.     private final Dfp[] sqr2Split;

  123.     /** A {@link Dfp} with value &radic;2 / 2. */
  124.     private final Dfp sqr2Reciprocal;

  126.     /** A {@link Dfp} with value &radic;3. */
  127.     private final Dfp sqr3;

  129.     /** A {@link Dfp} with value &radic;3 / 3. */
  130.     private final Dfp sqr3Reciprocal;

  132.     /** A {@link Dfp} with value &pi;. */
  133.     private final Dfp pi;

  135.     /** A two elements {@link Dfp} array with value &pi; split in two pieces. */
  136.     private final Dfp[] piSplit;

  138.     /** A {@link Dfp} with value e. */
  139.     private final Dfp e;

  141.     /** A two elements {@link Dfp} array with value e split in two pieces. */
  142.     private final Dfp[] eSplit;

  144.     /** A {@link Dfp} with value ln(2). */
  145.     private final Dfp ln2;

  147.     /** A two elements {@link Dfp} array with value ln(2) split in two pieces. */
  148.     private final Dfp[] ln2Split;

  150.     /** A {@link Dfp} with value ln(5). */
  151.     private final Dfp ln5;

  153.     /** A two elements {@link Dfp} array with value ln(5) split in two pieces. */
  154.     private final Dfp[] ln5Split;

  156.     /** A {@link Dfp} with value ln(10). */
  157.     private final Dfp ln10;

  159.     /** Current rounding mode. */
  160.     private RoundingMode rMode;

  162.     /** IEEE 854-1987 signals. */
  163.     private int ieeeFlags;

  165.     /** Create a factory for the specified number of radix digits.
  166.      * <p>
  167.      * Note that since the {@link Dfp} class uses 10000 as its radix, each radix
  168.      * digit is equivalent to 4 decimal digits. This implies that asking for
  169.      * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in
  170.      * all cases.
  171.      * </p>
  172.      * @param decimalDigits minimal number of decimal digits.
  173.      */
  174.     public DfpField(final int decimalDigits) {
  175.         this(decimalDigits, true);
  176.     }

  178.     /** Create a factory for the specified number of radix digits.
  179.      * <p>
  180.      * Note that since the {@link Dfp} class uses 10000 as its radix, each radix
  181.      * digit is equivalent to 4 decimal digits. This implies that asking for
  182.      * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in
  183.      * all cases.
  184.      * </p>
  185.      * @param decimalDigits minimal number of decimal digits
  186.      * @param computeConstants if true, the transcendental constants for the given precision
  187.      * must be computed (setting this flag to false is RESERVED for the internal recursive call)
  188.      */
  189.     private DfpField(final int decimalDigits, final boolean computeConstants) {

  191.         this.radixDigits = (decimalDigits < 13) ? 4 : (decimalDigits + 3) / 4;
  192.         this.rMode       = RoundingMode.ROUND_HALF_EVEN;
  193.         this.ieeeFlags   = 0;
  194.         this.zero        = new Dfp(this, 0);
  195.         this.one         = new Dfp(this, 1);
  196.         this.two         = new Dfp(this, 2);

  198.         if (computeConstants) {
  199.             // set up transcendental constants
  200.             synchronized (DfpField.class) {

  202.                 // as a heuristic to circumvent Table-Maker's Dilemma, we set the string
  203.                 // representation of the constants to be at least 3 times larger than the
  204.                 // number of decimal digits, also as an attempt to really compute these
  205.                 // constants only once, we set a minimum number of digits
  206.                 computeStringConstants((decimalDigits < 67) ? 200 : (3 * decimalDigits));

  208.                 // set up the constants at current field accuracy
  209.                 sqr2           = new Dfp(this, sqr2String);
  210.                 sqr2Split      = split(sqr2String);
  211.                 sqr2Reciprocal = new Dfp(this, sqr2ReciprocalString);
  212.                 sqr3           = new Dfp(this, sqr3String);
  213.                 sqr3Reciprocal = new Dfp(this, sqr3ReciprocalString);
  214.                 pi             = new Dfp(this, piString);
  215.                 piSplit        = split(piString);
  216.                 e              = new Dfp(this, eString);
  217.                 eSplit         = split(eString);
  218.                 ln2            = new Dfp(this, ln2String);
  219.                 ln2Split       = split(ln2String);
  220.                 ln5            = new Dfp(this, ln5String);
  221.                 ln5Split       = split(ln5String);
  222.                 ln10           = new Dfp(this, ln10String);
  223.             }
  224.         } else {
  225.             // dummy settings for unused constants
  226.             sqr2           = null;
  227.             sqr2Split      = null;
  228.             sqr2Reciprocal = null;
  229.             sqr3           = null;
  230.             sqr3Reciprocal = null;
  231.             pi             = null;
  232.             piSplit        = null;
  233.             e              = null;
  234.             eSplit         = null;
  235.             ln2            = null;
  236.             ln2Split       = null;
  237.             ln5            = null;
  238.             ln5Split       = null;
  239.             ln10           = null;
  240.         }
  241.     }

  243.     /** Get the number of radix digits of the {@link Dfp} instances built by this factory.
  244.      * @return number of radix digits
  245.      */
  246.     public int getRadixDigits() {
  247.         return radixDigits;
  248.     }

  250.     /** Set the rounding mode.
  251.      *  If not set, the default value is {@link RoundingMode#ROUND_HALF_EVEN}.
  252.      * @param mode desired rounding mode
  253.      * Note that the rounding mode is common to all {@link Dfp} instances
  254.      * belonging to the current {@link DfpField} in the system and will
  255.      * affect all future calculations.
  256.      */
  257.     public void setRoundingMode(final RoundingMode mode) {
  258.         rMode = mode;
  259.     }

  261.     /** Get the current rounding mode.
  262.      * @return current rounding mode
  263.      */
  264.     public RoundingMode getRoundingMode() {
  265.         return rMode;
  266.     }

  268.     /** Get the IEEE 854 status flags.
  269.      * @return IEEE 854 status flags
  270.      * @see #clearIEEEFlags()
  271.      * @see #setIEEEFlags(int)
  272.      * @see #setIEEEFlagsBits(int)
  273.      * @see #FLAG_INVALID
  274.      * @see #FLAG_DIV_ZERO
  275.      * @see #FLAG_OVERFLOW
  276.      * @see #FLAG_UNDERFLOW
  277.      * @see #FLAG_INEXACT
  278.      */
  279.     public int getIEEEFlags() {
  280.         return ieeeFlags;
  281.     }

  283.     /** Clears the IEEE 854 status flags.
  284.      * @see #getIEEEFlags()
  285.      * @see #setIEEEFlags(int)
  286.      * @see #setIEEEFlagsBits(int)
  287.      * @see #FLAG_INVALID
  288.      * @see #FLAG_DIV_ZERO
  289.      * @see #FLAG_OVERFLOW
  290.      * @see #FLAG_UNDERFLOW
  291.      * @see #FLAG_INEXACT
  292.      */
  293.     public void clearIEEEFlags() {
  294.         ieeeFlags = 0;
  295.     }

  297.     /** Sets the IEEE 854 status flags.
  298.      * @param flags desired value for the flags
  299.      * @see #getIEEEFlags()
  300.      * @see #clearIEEEFlags()
  301.      * @see #setIEEEFlagsBits(int)
  302.      * @see #FLAG_INVALID
  303.      * @see #FLAG_DIV_ZERO
  304.      * @see #FLAG_OVERFLOW
  305.      * @see #FLAG_UNDERFLOW
  306.      * @see #FLAG_INEXACT
  307.      */
  308.     public void setIEEEFlags(final int flags) {
  309.         ieeeFlags = flags & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT);
  310.     }

  312.     /** Sets some bits in the IEEE 854 status flags, without changing the already set bits.
  313.      * <p>
  314.      * Calling this method is equivalent to call {@code setIEEEFlags(getIEEEFlags() | bits)}
  315.      * </p>
  316.      * @param bits bits to set
  317.      * @see #getIEEEFlags()
  318.      * @see #clearIEEEFlags()
  319.      * @see #setIEEEFlags(int)
  320.      * @see #FLAG_INVALID
  321.      * @see #FLAG_DIV_ZERO
  322.      * @see #FLAG_OVERFLOW
  323.      * @see #FLAG_UNDERFLOW
  324.      * @see #FLAG_INEXACT
  325.      */
  326.     public void setIEEEFlagsBits(final int bits) {
  327.         ieeeFlags |= bits & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT);
  328.     }

  330.     /** Makes a {@link Dfp} with a value of 0.
  331.      * @return a new {@link Dfp} with a value of 0
  332.      */
  333.     public Dfp newDfp() {
  334.         return new Dfp(this);
  335.     }

  337.     /** Create an instance from a byte value.
  338.      * @param x value to convert to an instance
  339.      * @return a new {@link Dfp} with the same value as x
  340.      */
  341.     public Dfp newDfp(final byte x) {
  342.         return new Dfp(this, x);
  343.     }

  345.     /** Create an instance from an int value.
  346.      * @param x value to convert to an instance
  347.      * @return a new {@link Dfp} with the same value as x
  348.      */
  349.     public Dfp newDfp(final int x) {
  350.         return new Dfp(this, x);
  351.     }

  353.     /** Create an instance from a long value.
  354.      * @param x value to convert to an instance
  355.      * @return a new {@link Dfp} with the same value as x
  356.      */
  357.     public Dfp newDfp(final long x) {
  358.         return new Dfp(this, x);
  359.     }

  361.     /** Create an instance from a double value.
  362.      * @param x value to convert to an instance
  363.      * @return a new {@link Dfp} with the same value as x
  364.      */
  365.     public Dfp newDfp(final double x) {
  366.         return new Dfp(this, x);
  367.     }

  369.     /** Copy constructor.
  370.      * @param d instance to copy
  371.      * @return a new {@link Dfp} with the same value as d
  372.      */
  373.     public Dfp newDfp(Dfp d) {
  374.         return new Dfp(d);
  375.     }

  377.     /** Create a {@link Dfp} given a String representation.
  378.      * @param s string representation of the instance
  379.      * @return a new {@link Dfp} parsed from specified string
  380.      */
  381.     public Dfp newDfp(final String s) {
  382.         return new Dfp(this, s);
  383.     }

  385.     /** Creates a {@link Dfp} with a non-finite value.
  386.      * @param sign sign of the Dfp to create
  387.      * @param nans code of the value, must be one of {@link Dfp#INFINITE},
  388.      * {@link Dfp#SNAN},  {@link Dfp#QNAN}
  389.      * @return a new {@link Dfp} with a non-finite value
  390.      */
  391.     public Dfp newDfp(final byte sign, final byte nans) {
  392.         return new Dfp(this, sign, nans);
  393.     }

  395.     /** Get the constant 0.
  396.      * @return a {@link Dfp} with value 0
  397.      */
  398.     @Override
  399.     public Dfp getZero() {
  400.         return zero;
  401.     }

  403.     /** Get the constant 1.
  404.      * @return a {@link Dfp} with value 1
  405.      */
  406.     @Override
  407.     public Dfp getOne() {
  408.         return one;
  409.     }

  411.     /** {@inheritDoc} */
  412.     @Override
  413.     public Class<? extends FieldElement<Dfp>> getRuntimeClass() {
  414.         return Dfp.class;
  415.     }

  417.     /** Get the constant 2.
  418.      * @return a {@link Dfp} with value 2
  419.      */
  420.     public Dfp getTwo() {
  421.         return two;
  422.     }

  424.     /** Get the constant &radic;2.
  425.      * @return a {@link Dfp} with value &radic;2
  426.      */
  427.     public Dfp getSqr2() {
  428.         return sqr2;
  429.     }

  431.     /** Get the constant &radic;2 split in two pieces.
  432.      * @return a {@link Dfp} with value &radic;2 split in two pieces
  433.      */
  434.     public Dfp[] getSqr2Split() {
  435.         return sqr2Split.clone();
  436.     }

  438.     /** Get the constant &radic;2 / 2.
  439.      * @return a {@link Dfp} with value &radic;2 / 2
  440.      */
  441.     public Dfp getSqr2Reciprocal() {
  442.         return sqr2Reciprocal;
  443.     }

  445.     /** Get the constant &radic;3.
  446.      * @return a {@link Dfp} with value &radic;3
  447.      */
  448.     public Dfp getSqr3() {
  449.         return sqr3;
  450.     }

  452.     /** Get the constant &radic;3 / 3.
  453.      * @return a {@link Dfp} with value &radic;3 / 3
  454.      */
  455.     public Dfp getSqr3Reciprocal() {
  456.         return sqr3Reciprocal;
  457.     }

  459.     /** Get the constant &pi;.
  460.      * @return a {@link Dfp} with value &pi;
  461.      */
  462.     public Dfp getPi() {
  463.         return pi;
  464.     }

  466.     /** Get the constant &pi; split in two pieces.
  467.      * @return a {@link Dfp} with value &pi; split in two pieces
  468.      */
  469.     public Dfp[] getPiSplit() {
  470.         return piSplit.clone();
  471.     }

  473.     /** Get the constant e.
  474.      * @return a {@link Dfp} with value e
  475.      */
  476.     public Dfp getE() {
  477.         return e;
  478.     }

  480.     /** Get the constant e split in two pieces.
  481.      * @return a {@link Dfp} with value e split in two pieces
  482.      */
  483.     public Dfp[] getESplit() {
  484.         return eSplit.clone();
  485.     }

  487.     /** Get the constant ln(2).
  488.      * @return a {@link Dfp} with value ln(2)
  489.      */
  490.     public Dfp getLn2() {
  491.         return ln2;
  492.     }

  494.     /** Get the constant ln(2) split in two pieces.
  495.      * @return a {@link Dfp} with value ln(2) split in two pieces
  496.      */
  497.     public Dfp[] getLn2Split() {
  498.         return ln2Split.clone();
  499.     }

  501.     /** Get the constant ln(5).
  502.      * @return a {@link Dfp} with value ln(5)
  503.      */
  504.     public Dfp getLn5() {
  505.         return ln5;
  506.     }

  508.     /** Get the constant ln(5) split in two pieces.
  509.      * @return a {@link Dfp} with value ln(5) split in two pieces
  510.      */
  511.     public Dfp[] getLn5Split() {
  512.         return ln5Split.clone();
  513.     }

  515.     /** Get the constant ln(10).
  516.      * @return a {@link Dfp} with value ln(10)
  517.      */
  518.     public Dfp getLn10() {
  519.         return ln10;
  520.     }

  522.     /** Breaks a string representation up into two {@link Dfp}'s.
  523.      * The split is such that the sum of them is equivalent to the input string,
  524.      * but has higher precision than using a single Dfp.
  525.      * @param a string representation of the number to split
  526.      * @return an array of two {@link Dfp Dfp} instances which sum equals a
  527.      */
  528.     private Dfp[] split(final String a) {
  529.         Dfp[] result = new Dfp[2];
  530.         boolean leading = true;
  531.         int sp = 0;
  532.         int sig = 0;

  534.         char[] buf = new char[a.length()];

  536.         for (int i = 0; i < buf.length; i++) {
  537.             buf[i] = a.charAt(i);

  539.             if (buf[i] >= '1' && buf[i] <= '9') {
  540.                 leading = false;
  541.             }

  543.             if (buf[i] == '.') {
  544.                 sig += (400 - sig) % 4;
  545.                 leading = false;
  546.             }

  548.             if (sig == (radixDigits / 2) * 4) {
  549.                 sp = i;
  550.                 break;
  551.             }

  553.             if (buf[i] >= '0' && buf[i] <= '9' && !leading) {
  554.                 sig++;
  555.             }
  556.         }

  558.         result[0] = new Dfp(this, String.valueOf(buf, 0, sp));

  560.         for (int i = 0; i < buf.length; i++) {
  561.             buf[i] = a.charAt(i);
  562.             if (buf[i] >= '0' && buf[i] <= '9' && i < sp) {
  563.                 buf[i] = '0';
  564.             }
  565.         }

  567.         result[1] = new Dfp(this, String.valueOf(buf));

  569.         return result;
  570.     }

  572.     /** Recompute the high precision string constants.
  573.      * @param highPrecisionDecimalDigits precision at which the string constants mus be computed
  574.      */
  575.     private static void computeStringConstants(final int highPrecisionDecimalDigits) {
  576.         if (sqr2String == null || sqr2String.length() < highPrecisionDecimalDigits - 3) {

  578.             // recompute the string representation of the transcendental constants
  579.             final DfpField highPrecisionField = new DfpField(highPrecisionDecimalDigits, false);
  580.             final Dfp highPrecisionOne        = new Dfp(highPrecisionField, 1);
  581.             final Dfp highPrecisionTwo        = new Dfp(highPrecisionField, 2);
  582.             final Dfp highPrecisionThree      = new Dfp(highPrecisionField, 3);

  584.             final Dfp highPrecisionSqr2 = highPrecisionTwo.sqrt();
  585.             sqr2String           = highPrecisionSqr2.toString();
  586.             sqr2ReciprocalString = highPrecisionOne.divide(highPrecisionSqr2).toString();

  588.             final Dfp highPrecisionSqr3 = highPrecisionThree.sqrt();
  589.             sqr3String           = highPrecisionSqr3.toString();
  590.             sqr3ReciprocalString = highPrecisionOne.divide(highPrecisionSqr3).toString();

  592.             piString   = computePi(highPrecisionOne, highPrecisionTwo, highPrecisionThree).toString();
  593.             eString    = computeExp(highPrecisionOne, highPrecisionOne).toString();
  594.             ln2String  = computeLn(highPrecisionTwo, highPrecisionOne, highPrecisionTwo).toString();
  595.             ln5String  = computeLn(new Dfp(highPrecisionField, 5),  highPrecisionOne, highPrecisionTwo).toString();
  596.             ln10String = computeLn(new Dfp(highPrecisionField, 10), highPrecisionOne, highPrecisionTwo).toString();
  597.         }
  598.     }

  600.     /** Compute &pi; using Jonathan and Peter Borwein quartic formula.
  601.      * @param one constant with value 1 at desired precision
  602.      * @param two constant with value 2 at desired precision
  603.      * @param three constant with value 3 at desired precision
  604.      * @return &pi;
  605.      */
  606.     private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) {

  608.         Dfp sqrt2   = two.sqrt();
  609.         Dfp yk      = sqrt2.subtract(one);
  610.         Dfp four    = two.add(two);
  611.         Dfp two2kp3 = two;
  612.         Dfp ak      = two.multiply(three.subtract(two.multiply(sqrt2)));

  614.         // The formula converges quartically. This means the number of correct
  615.         // digits is multiplied by 4 at each iteration! Five iterations are
  616.         // sufficient for about 160 digits, eight iterations give about
  617.         // 10000 digits (this has been checked) and 20 iterations more than
  618.         // 160 billions of digits (this has NOT been checked).
  619.         // So the limit here is considered sufficient for most purposes ...
  620.         for (int i = 1; i < 20; i++) {
  621.             final Dfp ykM1 = yk;

  623.             final Dfp y2         = yk.multiply(yk);
  624.             final Dfp oneMinusY4 = one.subtract(y2.multiply(y2));
  625.             final Dfp s          = oneMinusY4.sqrt().sqrt();
  626.             yk = one.subtract(s).divide(one.add(s));

  628.             two2kp3 = two2kp3.multiply(four);

  630.             final Dfp p = one.add(yk);
  631.             final Dfp p2 = p.multiply(p);
  632.             ak = ak.multiply(p2.multiply(p2)).subtract(two2kp3.multiply(yk).multiply(one.add(yk).add(yk.multiply(yk))));

  634.             if (yk.equals(ykM1)) {
  635.                 break;
  636.             }
  637.         }

  639.         return one.divide(ak);
  640.     }

  642.     /** Compute exp(a).
  643.      * @param a number for which we want the exponential
  644.      * @param one constant with value 1 at desired precision
  645.      * @return exp(a)
  646.      */
  647.     public static Dfp computeExp(final Dfp a, final Dfp one) {

  649.         Dfp y  = new Dfp(one);
  650.         Dfp py = new Dfp(one);
  651.         Dfp f  = new Dfp(one);
  652.         Dfp fi = new Dfp(one);
  653.         Dfp x  = new Dfp(one);

  655.         for (int i = 0; i < 10000; i++) {
  656.             x = x.multiply(a);
  657.             y = y.add(x.divide(f));
  658.             fi = fi.add(one);
  659.             f = f.multiply(fi);
  660.             if (y.equals(py)) {
  661.                 break;
  662.             }
  663.             py = new Dfp(y);
  664.         }

  666.         return y;
  667.     }


  670.     /** Compute ln(a).
  671.      *
  672.      *  <pre>{@code
  673.      *  Let f(x) = ln(x),
  674.      *
  675.      *  We know that f'(x) = 1/x, thus from Taylor's theorem we have:
  676.      *
  677.      *           -----          n+1         n
  678.      *  f(x) =   \           (-1)    (x - 1)
  679.      *           /          ----------------    for 1 <= n <= infinity
  680.      *           -----             n
  681.      *
  682.      *  or
  683.      *                       2        3       4
  684.      *                   (x-1)   (x-1)    (x-1)
  685.      *  ln(x) =  (x-1) - ----- + ------ - ------ + ...
  686.      *                     2       3        4
  687.      *
  688.      *  alternatively,
  689.      *
  690.      *                  2    3   4
  691.      *                 x    x   x
  692.      *  ln(x+1) =  x - -  + - - - + ...
  693.      *                 2    3   4
  694.      *
  695.      *  This series can be used to compute ln(x), but it converges too slowly.
  696.      *
  697.      *  If we substitute -x for x above, we get
  698.      *
  699.      *                   2    3    4
  700.      *                  x    x    x
  701.      *  ln(1-x) =  -x - -  - -  - - + ...
  702.      *                  2    3    4
  703.      *
  704.      *  Note that all terms are now negative.  Because the even powered ones
  705.      *  absorbed the sign.  Now, subtract the series above from the previous
  706.      *  one to get ln(x+1) - ln(1-x).  Note the even terms cancel out leaving
  707.      *  only the odd ones
  708.      *
  709.      *                             3     5      7
  710.      *                           2x    2x     2x
  711.      *  ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
  712.      *                            3     5      7
  713.      *
  714.      *  By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
  715.      *
  716.      *                                3        5        7
  717.      *      x+1           /          x        x        x          \
  718.      *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
  719.      *      x-1           \          3        5        7          /
  720.      *
  721.      *  But now we want to find ln(a), so we need to find the value of x
  722.      *  such that a = (x+1)/(x-1).   This is easily solved to find that
  723.      *  x = (a-1)/(a+1).
  724.      * }</pre>
  725.      * @param a number for which we want the exponential
  726.      * @param one constant with value 1 at desired precision
  727.      * @param two constant with value 2 at desired precision
  728.      * @return ln(a)
  729.      */

  731.     public static Dfp computeLn(final Dfp a, final Dfp one, final Dfp two) {

  733.         int den = 1;
  734.         Dfp x = a.add(new Dfp(a.getField(), -1)).divide(a.add(one));

  736.         Dfp y = new Dfp(x);
  737.         Dfp num = new Dfp(x);
  738.         Dfp py = new Dfp(y);
  739.         for (int i = 0; i < 10000; i++) {
  740.             num = num.multiply(x);
  741.             num = num.multiply(x);
  742.             den += 2;
  743.             Dfp t = num.divide(den);
  744.             y = y.add(t);
  745.             if (y.equals(py)) {
  746.                 break;
  747.             }
  748.             py = new Dfp(y);
  749.         }

  751.         return y.multiply(two);
  752.     }
  753. }
Created with JaCoCo 0.8.10.202304240956