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.  *      https://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. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.dfp;

  25. import org.hipparchus.Field;

  27. /** Field for Decimal floating point instances.
  28.  */
  29. public class DfpField implements Field<Dfp> {

  31.     /** Enumerate for rounding modes. */
  32.     public enum RoundingMode {

  34.         /** Rounds toward zero (truncation). */
  35.         ROUND_DOWN,

  37.         /** Rounds away from zero if discarded digit is non-zero. */
  38.         ROUND_UP,

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

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

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

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

  54.         /** Rounds towards positive infinity. */
  55.         ROUND_CEIL,

  57.         /** Rounds towards negative infinity. */
  58.         ROUND_FLOOR

  60.     }

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

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

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

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

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

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

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

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

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

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

  91.     /** High precision string representation of &pi;. */
  92.     private static String piString;

  94.     /** High precision string representation of e. */
  95.     private static String eString;

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

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

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

  106.     /** The number of radix digits.
  107.      * Note these depend on the radix which is 10000 digits,
  108.      * so each one is equivalent to 4 decimal digits.
  109.      */
  110.     private final int radixDigits;

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

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

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

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

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

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

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

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

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

  139.     /** A {@link Dfp} for converting degrees to radians. */
  140.     private final Dfp degToRad;

  142.     /** A {@link Dfp} for converting radians to degrees. */
  143.     private final Dfp radToDeg;

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

  148.     /** A {@link Dfp} with value e. */
  149.     private final Dfp e;

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

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

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

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

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

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

  169.     /** Current rounding mode. */
  170.     private RoundingMode rMode;

  172.     /** IEEE 854-1987 signals. */
  173.     private int ieeeFlags;

  175.     /** Create a factory for the specified number of radix digits.
  176.      * <p>
  177.      * Note that since the {@link Dfp} class uses 10000 as its radix, each radix
  178.      * digit is equivalent to 4 decimal digits. This implies that asking for
  179.      * 13, 14, 15 or 16 decimal digits will really lead to a 4 radix 10000 digits in
  180.      * all cases.
  181.      * </p>
  182.      * @param decimalDigits minimal number of decimal digits.
  183.      */
  184.     public DfpField(final int decimalDigits) {
  185.         this(decimalDigits, true);
  186.     }

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

  201.         this.radixDigits = (decimalDigits < 13) ? 4 : (decimalDigits + 3) / 4;
  202.         this.rMode       = RoundingMode.ROUND_HALF_EVEN;
  203.         this.ieeeFlags   = 0;
  204.         this.zero        = new Dfp(this, 0);
  205.         this.one         = new Dfp(this, 1);
  206.         this.two         = new Dfp(this, 2);

  208.         if (computeConstants) {
  209.             // set up transcendental constants
  210.             synchronized (DfpField.class) {

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

  218.                 // set up the constants at current field accuracy
  219.                 sqr2           = new Dfp(this, sqr2String);
  220.                 sqr2Split      = split(sqr2String);
  221.                 sqr2Reciprocal = new Dfp(this, sqr2ReciprocalString);
  222.                 sqr3           = new Dfp(this, sqr3String);
  223.                 sqr3Reciprocal = new Dfp(this, sqr3ReciprocalString);
  224.                 pi             = new Dfp(this, piString);
  225.                 degToRad       = pi.divide(180.0);
  226.                 radToDeg       = pi.divide(180.0).reciprocal();
  227.                 piSplit        = split(piString);
  228.                 e              = new Dfp(this, eString);
  229.                 eSplit         = split(eString);
  230.                 ln2            = new Dfp(this, ln2String);
  231.                 ln2Split       = split(ln2String);
  232.                 ln5            = new Dfp(this, ln5String);
  233.                 ln5Split       = split(ln5String);
  234.                 ln10           = new Dfp(this, ln10String);

  236.             }
  237.         } else {
  238.             // dummy settings for unused constants
  239.             sqr2           = null;
  240.             sqr2Split      = null;
  241.             sqr2Reciprocal = null;
  242.             sqr3           = null;
  243.             sqr3Reciprocal = null;
  244.             pi             = null;
  245.             degToRad       = null;
  246.             radToDeg       = null;
  247.             piSplit        = null;
  248.             e              = null;
  249.             eSplit         = null;
  250.             ln2            = null;
  251.             ln2Split       = null;
  252.             ln5            = null;
  253.             ln5Split       = null;
  254.             ln10           = null;
  255.         }

  257.     }

  259.     /** Get the number of radix digits of the {@link Dfp} instances built by this factory.
  260.      * @return number of radix digits
  261.      */
  262.     public int getRadixDigits() {
  263.         return radixDigits;
  264.     }

  266.     /** Set the rounding mode.
  267.      *  If not set, the default value is {@link RoundingMode#ROUND_HALF_EVEN}.
  268.      * @param mode desired rounding mode
  269.      * Note that the rounding mode is common to all {@link Dfp} instances
  270.      * belonging to the current {@link DfpField} in the system and will
  271.      * affect all future calculations.
  272.      */
  273.     public void setRoundingMode(final RoundingMode mode) {
  274.         rMode = mode;
  275.     }

  277.     /** Get the current rounding mode.
  278.      * @return current rounding mode
  279.      */
  280.     public RoundingMode getRoundingMode() {
  281.         return rMode;
  282.     }

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

  299.     /** Clears the IEEE 854 status flags.
  300.      * @see #getIEEEFlags()
  301.      * @see #setIEEEFlags(int)
  302.      * @see #setIEEEFlagsBits(int)
  303.      * @see #FLAG_INVALID
  304.      * @see #FLAG_DIV_ZERO
  305.      * @see #FLAG_OVERFLOW
  306.      * @see #FLAG_UNDERFLOW
  307.      * @see #FLAG_INEXACT
  308.      */
  309.     public void clearIEEEFlags() {
  310.         ieeeFlags = 0;
  311.     }

  313.     /** Sets the IEEE 854 status flags.
  314.      * @param flags desired value for the flags
  315.      * @see #getIEEEFlags()
  316.      * @see #clearIEEEFlags()
  317.      * @see #setIEEEFlagsBits(int)
  318.      * @see #FLAG_INVALID
  319.      * @see #FLAG_DIV_ZERO
  320.      * @see #FLAG_OVERFLOW
  321.      * @see #FLAG_UNDERFLOW
  322.      * @see #FLAG_INEXACT
  323.      */
  324.     public void setIEEEFlags(final int flags) {
  325.         ieeeFlags = flags & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT);
  326.     }

  328.     /** Sets some bits in the IEEE 854 status flags, without changing the already set bits.
  329.      * <p>
  330.      * Calling this method is equivalent to call {@code setIEEEFlags(getIEEEFlags() | bits)}
  331.      * </p>
  332.      * @param bits bits to set
  333.      * @see #getIEEEFlags()
  334.      * @see #clearIEEEFlags()
  335.      * @see #setIEEEFlags(int)
  336.      * @see #FLAG_INVALID
  337.      * @see #FLAG_DIV_ZERO
  338.      * @see #FLAG_OVERFLOW
  339.      * @see #FLAG_UNDERFLOW
  340.      * @see #FLAG_INEXACT
  341.      */
  342.     public void setIEEEFlagsBits(final int bits) {
  343.         ieeeFlags |= bits & (FLAG_INVALID | FLAG_DIV_ZERO | FLAG_OVERFLOW | FLAG_UNDERFLOW | FLAG_INEXACT);
  344.     }

  346.     /** Makes a {@link Dfp} with a value of 0.
  347.      * @return a new {@link Dfp} with a value of 0
  348.      */
  349.     public Dfp newDfp() {
  350.         return new Dfp(this);
  351.     }

  353.     /** Create an instance from a byte 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 byte x) {
  358.         return new Dfp(this, x);
  359.     }

  361.     /** Create an instance from an int 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 int x) {
  366.         return new Dfp(this, x);
  367.     }

  369.     /** Create an instance from a long value.
  370.      * @param x value to convert to an instance
  371.      * @return a new {@link Dfp} with the same value as x
  372.      */
  373.     public Dfp newDfp(final long x) {
  374.         return new Dfp(this, x);
  375.     }

  377.     /** Create an instance from a double value.
  378.      * @param x value to convert to an instance
  379.      * @return a new {@link Dfp} with the same value as x
  380.      */
  381.     public Dfp newDfp(final double x) {
  382.         return new Dfp(this, x);
  383.     }

  385.     /** Copy constructor.
  386.      * @param d instance to copy
  387.      * @return a new {@link Dfp} with the same value as d
  388.      */
  389.     public Dfp newDfp(Dfp d) {
  390.         return new Dfp(d);
  391.     }

  393.     /** Create a {@link Dfp} given a String representation.
  394.      * @param s string representation of the instance
  395.      * @return a new {@link Dfp} parsed from specified string
  396.      */
  397.     public Dfp newDfp(final String s) {
  398.         return new Dfp(this, s);
  399.     }

  401.     /** Creates a {@link Dfp} with a non-finite value.
  402.      * @param sign sign of the Dfp to create
  403.      * @param nans code of the value, must be one of {@link Dfp#INFINITE},
  404.      * {@link Dfp#SNAN},  {@link Dfp#QNAN}
  405.      * @return a new {@link Dfp} with a non-finite value
  406.      */
  407.     public Dfp newDfp(final byte sign, final byte nans) {
  408.         return new Dfp(this, sign, nans);
  409.     }

  411.     /** Get the constant 0.
  412.      * @return a {@link Dfp} with value 0
  413.      */
  414.     @Override
  415.     public Dfp getZero() {
  416.         return zero;
  417.     }

  419.     /** Get the constant 1.
  420.      * @return a {@link Dfp} with value 1
  421.      */
  422.     @Override
  423.     public Dfp getOne() {
  424.         return one;
  425.     }

  427.     /** {@inheritDoc} */
  428.     @Override
  429.     public Class<Dfp> getRuntimeClass() {
  430.         return Dfp.class;
  431.     }

  433.     /** Get the constant 2.
  434.      * @return a {@link Dfp} with value 2
  435.      */
  436.     public Dfp getTwo() {
  437.         return two;
  438.     }

  440.     /** Get the constant &radic;2.
  441.      * @return a {@link Dfp} with value &radic;2
  442.      */
  443.     public Dfp getSqr2() {
  444.         return sqr2;
  445.     }

  447.     /** Get the constant &radic;2 split in two pieces.
  448.      * @return a {@link Dfp} with value &radic;2 split in two pieces
  449.      */
  450.     public Dfp[] getSqr2Split() {
  451.         return sqr2Split.clone();
  452.     }

  454.     /** Get the constant &radic;2 / 2.
  455.      * @return a {@link Dfp} with value &radic;2 / 2
  456.      */
  457.     public Dfp getSqr2Reciprocal() {
  458.         return sqr2Reciprocal;
  459.     }

  461.     /** Get the constant &radic;3.
  462.      * @return a {@link Dfp} with value &radic;3
  463.      */
  464.     public Dfp getSqr3() {
  465.         return sqr3;
  466.     }

  468.     /** Get the constant &radic;3 / 3.
  469.      * @return a {@link Dfp} with value &radic;3 / 3
  470.      */
  471.     public Dfp getSqr3Reciprocal() {
  472.         return sqr3Reciprocal;
  473.     }

  475.     /** Get the constant &pi;.
  476.      * @return a {@link Dfp} with value &pi;
  477.      */
  478.     public Dfp getPi() {
  479.         return pi;
  480.     }

  482.     /** Get the degrees to radians conversion factor.
  483.      * @return a {@link Dfp} for degrees to radians conversion factor
  484.      */
  485.     public Dfp getDegToRad() {
  486.         return degToRad;
  487.     }

  489.     /** Get the radians to degrees conversion factor.
  490.      * @return a {@link Dfp} for radians to degrees conversion factor
  491.      */
  492.     public Dfp getRadToDeg() {
  493.         return radToDeg;
  494.     }

  496.     /** Get the constant &pi; split in two pieces.
  497.      * @return a {@link Dfp} with value &pi; split in two pieces
  498.      */
  499.     public Dfp[] getPiSplit() {
  500.         return piSplit.clone();
  501.     }

  503.     /** Get the constant e.
  504.      * @return a {@link Dfp} with value e
  505.      */
  506.     public Dfp getE() {
  507.         return e;
  508.     }

  510.     /** Get the constant e split in two pieces.
  511.      * @return a {@link Dfp} with value e split in two pieces
  512.      */
  513.     public Dfp[] getESplit() {
  514.         return eSplit.clone();
  515.     }

  517.     /** Get the constant ln(2).
  518.      * @return a {@link Dfp} with value ln(2)
  519.      */
  520.     public Dfp getLn2() {
  521.         return ln2;
  522.     }

  524.     /** Get the constant ln(2) split in two pieces.
  525.      * @return a {@link Dfp} with value ln(2) split in two pieces
  526.      */
  527.     public Dfp[] getLn2Split() {
  528.         return ln2Split.clone();
  529.     }

  531.     /** Get the constant ln(5).
  532.      * @return a {@link Dfp} with value ln(5)
  533.      */
  534.     public Dfp getLn5() {
  535.         return ln5;
  536.     }

  538.     /** Get the constant ln(5) split in two pieces.
  539.      * @return a {@link Dfp} with value ln(5) split in two pieces
  540.      */
  541.     public Dfp[] getLn5Split() {
  542.         return ln5Split.clone();
  543.     }

  545.     /** Get the constant ln(10).
  546.      * @return a {@link Dfp} with value ln(10)
  547.      */
  548.     public Dfp getLn10() {
  549.         return ln10;
  550.     }

  552.     /** {@inheritDoc}
  553.      * <p>
  554.      * Two fields are considered equals if they have the same number
  555.      * of radix digits and the same rounding mode.
  556.      * </p>
  557.      */
  558.     @Override
  559.     public boolean equals(final Object other) {
  560.         if (this == other) {
  561.             return true;
  562.         } else if (other instanceof DfpField) {
  563.             DfpField rhs = (DfpField) other;
  564.             return getRadixDigits()  == rhs.getRadixDigits()  &&
  565.                    getRoundingMode() == rhs.getRoundingMode();
  566.         } else {
  567.             return false;
  568.         }
  569.     }

  571.     /** {@inheritDoc} */
  572.     @Override
  573.     public int hashCode() {
  574.         return 0xdf49a2ca ^ ((radixDigits << 16) & (rMode.ordinal() << 5) & ieeeFlags);
  575.     }

  577.     /** Breaks a string representation up into two {@link Dfp}'s.
  578.      * The split is such that the sum of them is equivalent to the input string,
  579.      * but has higher precision than using a single Dfp.
  580.      * @param a string representation of the number to split
  581.      * @return an array of two {@link Dfp Dfp} instances which sum equals a
  582.      */
  583.     private Dfp[] split(final String a) {
  584.       Dfp[] result = new Dfp[2];
  585.       boolean leading = true;
  586.       int sp = 0;
  587.       int sig = 0;

  589.       StringBuilder builder1 = new StringBuilder(a.length());

  591.       for (int i = 0; i < a.length(); i++) {
  592.         final char c = a.charAt(i);
  593.         builder1.append(c);

  595.         if (c >= '1' && c <= '9') {
  596.             leading = false;
  597.         }

  599.         if (c == '.') {
  600.           sig += (400 - sig) % 4;
  601.           leading = false;
  602.         }

  604.         if (sig == (radixDigits / 2) * 4) {
  605.           sp = i;
  606.           break;
  607.         }

  609.         if (c >= '0' &&c <= '9' && !leading) {
  610.             sig ++;
  611.         }
  612.       }

  614.       result[0] = new Dfp(this, builder1.substring(0, sp));

  616.       StringBuilder builder2 = new StringBuilder(a.length());
  617.       for (int i = 0; i < a.length(); i++) {
  618.           final char c = a.charAt(i);
  619.           if (c >= '0' && c <= '9' && i < sp) {
  620.               builder2.append('0');
  621.           } else {
  622.               builder2.append(c);
  623.         }
  624.       }

  626.       result[1] = new Dfp(this, builder2.toString());

  628.       return result;

  630.     }

  632.     /** Recompute the high precision string constants.
  633.      * @param highPrecisionDecimalDigits precision at which the string constants mus be computed
  634.      */
  635.     private static void computeStringConstants(final int highPrecisionDecimalDigits) {
  636.         synchronized (DfpField.class) {
  637.             if (sqr2String == null || sqr2String.length() < highPrecisionDecimalDigits - 3) {

  639.                 // recompute the string representation of the transcendental constants
  640.                 final DfpField highPrecisionField = new DfpField(highPrecisionDecimalDigits, false);
  641.                 final Dfp highPrecisionOne        = new Dfp(highPrecisionField, 1);
  642.                 final Dfp highPrecisionTwo        = new Dfp(highPrecisionField, 2);
  643.                 final Dfp highPrecisionThree      = new Dfp(highPrecisionField, 3);

  645.                 final Dfp highPrecisionSqr2 = highPrecisionTwo.sqrt();
  646.                 sqr2String           = highPrecisionSqr2.toString();
  647.                 sqr2ReciprocalString = highPrecisionOne.divide(highPrecisionSqr2).toString();

  649.                 final Dfp highPrecisionSqr3 = highPrecisionThree.sqrt();
  650.                 sqr3String           = highPrecisionSqr3.toString();
  651.                 sqr3ReciprocalString = highPrecisionOne.divide(highPrecisionSqr3).toString();

  653.                 piString   = computePi(highPrecisionOne, highPrecisionTwo, highPrecisionThree).toString();
  654.                 eString    = computeExp(highPrecisionOne, highPrecisionOne).toString();
  655.                 ln2String  = computeLn(highPrecisionTwo, highPrecisionOne, highPrecisionTwo).toString();
  656.                 ln5String  = computeLn(new Dfp(highPrecisionField, 5),  highPrecisionOne, highPrecisionTwo).toString();
  657.                 ln10String = computeLn(new Dfp(highPrecisionField, 10), highPrecisionOne, highPrecisionTwo).toString();

  659.             }
  660.         }
  661.     }

  663.     /** Compute &pi; using Jonathan and Peter Borwein quartic formula.
  664.      * @param one constant with value 1 at desired precision
  665.      * @param two constant with value 2 at desired precision
  666.      * @param three constant with value 3 at desired precision
  667.      * @return &pi;
  668.      */
  669.     private static Dfp computePi(final Dfp one, final Dfp two, final Dfp three) {

  671.         Dfp sqrt2   = two.sqrt();
  672.         Dfp yk      = sqrt2.subtract(one);
  673.         Dfp four    = two.add(two);
  674.         Dfp two2kp3 = two;
  675.         Dfp ak      = two.multiply(three.subtract(two.multiply(sqrt2)));

  677.         // The formula converges quartically. This means the number of correct
  678.         // digits is multiplied by 4 at each iteration! Five iterations are
  679.         // sufficient for about 160 digits, eight iterations give about
  680.         // 10000 digits (this has been checked) and 20 iterations more than
  681.         // 160 billions of digits (this has NOT been checked).
  682.         // So the limit here is considered sufficient for most purposes ...
  683.         for (int i = 1; i < 20; i++) {
  684.             final Dfp ykM1 = yk;

  686.             final Dfp y2         = yk.multiply(yk);
  687.             final Dfp oneMinusY4 = one.subtract(y2.multiply(y2));
  688.             final Dfp s          = oneMinusY4.sqrt().sqrt();
  689.             yk = one.subtract(s).divide(one.add(s));

  691.             two2kp3 = two2kp3.multiply(four);

  693.             final Dfp p = one.add(yk);
  694.             final Dfp p2 = p.multiply(p);
  695.             ak = ak.multiply(p2.multiply(p2)).subtract(two2kp3.multiply(yk).multiply(one.add(yk).add(yk.multiply(yk))));

  697.             if (yk.equals(ykM1)) {
  698.                 break;
  699.             }
  700.         }

  702.         return one.divide(ak);

  704.     }

  706.     /** Compute exp(a).
  707.      * @param a number for which we want the exponential
  708.      * @param one constant with value 1 at desired precision
  709.      * @return exp(a)
  710.      */
  711.     public static Dfp computeExp(final Dfp a, final Dfp one) {

  713.         Dfp y  = new Dfp(one);
  714.         Dfp py = new Dfp(one);
  715.         Dfp f  = new Dfp(one);
  716.         Dfp fi = new Dfp(one);
  717.         Dfp x  = new Dfp(one);

  719.         for (int i = 0; i < 10000; i++) {
  720.             x = x.multiply(a);
  721.             y = y.add(x.divide(f));
  722.             fi = fi.add(one);
  723.             f = f.multiply(fi);
  724.             if (y.equals(py)) {
  725.                 break;
  726.             }
  727.             py = new Dfp(y);
  728.         }

  730.         return y;

  732.     }


  735.     /** Compute ln(a).
  736.      *
  737.      *  Let f(x) = ln(x),
  738.      *
  739.      *  We know that f'(x) = 1/x, thus from Taylor's theorem we have:
  740.      *
  741.      *           -----          n+1         n
  742.      *  f(x) =   \           (-1)    (x - 1)
  743.      *           /          ----------------    for 1 &lt;= n &lt;= infinity
  744.      *           -----             n
  745.      *
  746.      *  or
  747.      *                       2        3       4
  748.      *                   (x-1)   (x-1)    (x-1)
  749.      *  ln(x) =  (x-1) - ----- + ------ - ------ + ...
  750.      *                     2       3        4
  751.      *
  752.      *  alternatively,
  753.      *
  754.      *                  2    3   4
  755.      *                 x    x   x
  756.      *  ln(x+1) =  x - -  + - - - + ...
  757.      *                 2    3   4
  758.      *
  759.      *  This series can be used to compute ln(x), but it converges too slowly.
  760.      *
  761.      *  If we substitute -x for x above, we get
  762.      *
  763.      *                   2    3    4
  764.      *                  x    x    x
  765.      *  ln(1-x) =  -x - -  - -  - - + ...
  766.      *                  2    3    4
  767.      *
  768.      *  Note that all terms are now negative.  Because the even powered ones
  769.      *  absorbed the sign.  Now, subtract the series above from the previous
  770.      *  one to get ln(x+1) - ln(1-x).  Note the even terms cancel out leaving
  771.      *  only the odd ones
  772.      *
  773.      *                             3     5      7
  774.      *                           2x    2x     2x
  775.      *  ln(x+1) - ln(x-1) = 2x + --- + --- + ---- + ...
  776.      *                            3     5      7
  777.      *
  778.      *  By the property of logarithms that ln(a) - ln(b) = ln (a/b) we have:
  779.      *
  780.      *                                3        5        7
  781.      *      x+1           /          x        x        x          \
  782.      *  ln ----- =   2 *  |  x  +   ----  +  ----  +  ---- + ...  |
  783.      *      x-1           \          3        5        7          /
  784.      *
  785.      *  But now we want to find ln(a), so we need to find the value of x
  786.      *  such that a = (x+1)/(x-1).   This is easily solved to find that
  787.      *  x = (a-1)/(a+1).
  788.      * @param a number for which we want the exponential
  789.      * @param one constant with value 1 at desired precision
  790.      * @param two constant with value 2 at desired precision
  791.      * @return ln(a)
  792.      */

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

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

  799.         Dfp y = new Dfp(x);
  800.         Dfp num = new Dfp(x);
  801.         Dfp py = new Dfp(y);
  802.         for (int i = 0; i < 10000; i++) {
  803.             num = num.multiply(x);
  804.             num = num.multiply(x);
  805.             den += 2;
  806.             Dfp t = num.divide(den);
  807.             y = y.add(t);
  808.             if (y.equals(py)) {
  809.                 break;
  810.             }
  811.             py = new Dfp(y);
  812.         }

  814.         return y.multiply(two);

  816.     }

  818.     /** Get extended field for accuracy conversion.
  819.      * @param digitsFactor multiplication factor for number of digits
  820.      * @param computeConstants if true, the transcendental constants for the given precision
  821.      * must be computed (setting this flag to false is RESERVED for the internal recursive call)
  822.      * @return field with extended precision
  823.      * @since 1.7
  824.      */
  825.     public DfpField getExtendedField(final int digitsFactor, final boolean computeConstants) {
  826.         final int oldDecimalDigits = getRadixDigits() * 4;
  827.         return new DfpField(oldDecimalDigits * digitsFactor, computeConstants);
  828.     }

  830. }
Created with JaCoCo 0.8.12.202403310830