CarlsonEllipticIntegral.java

  1. /*
  2.  * Licensed to the Hipparchus project 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 Hipparchus project 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.  */
  17. package org.hipparchus.special.elliptic.carlson;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.complex.Complex;
  21. import org.hipparchus.complex.FieldComplex;
  22. import org.hipparchus.util.FastMath;

  24. /** Elliptic integrals in Carlson symmetric form.
  25.  * <p>
  26.  * This utility class computes the various symmetric elliptic
  27.  * integrals defined as:
  28.  * \[
  29.  *   \left\{\begin{align}
  30.  *   R_F(x,y,z)   &amp;= \frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{s(t)}\\
  31.  *   R_J(x,y,z,p) &amp;= \frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{s(t)(t+p)}\\
  32.  *   R_G(x,y,z)   &amp;= \frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
  33.                      \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t\\
  34.  *   R_D(x,y,z)   &amp;= R_J(x,y,z,z)\\
  35.  *   R_C(x,y)     &amp;= R_F(x,y,y)
  36.  *   \end{align}\right.
  37.  * \]
  38.  * </p>
  39.  * <p>
  40.  * where
  41.  * \[
  42.  *   s(t) = \sqrt{t+x}\sqrt{t+y}\sqrt{t+z}
  43.  * \]
  44.  * </p>
  45.  * <p>
  46.  * The algorithms used are based on the duplication method as described in
  47.  * B. C. Carlson 1995 paper "Numerical computation of real or complex
  48.  * elliptic integrals", with the improvements described in the appendix
  49.  * of B. C. Carlson and James FitzSimons 2000 paper "Reduction theorems
  50.  * for elliptic integrands with the square root of two quadratic factors".
  51.  * They are also described in <a href="https://dlmf.nist.gov/19.36#i">section 19.36(i)</a>
  52.  * of Digital Library of Mathematical Functions.
  53.  * </p>
  54.  * <p>
  55.  * <em>
  56.  * Beware that when computing elliptic integrals in the complex plane,
  57.  * many issues arise due to branch cuts. See the
  58.  * <a href="https://www.hipparchus.org/hipparchus-core/special.html#Elliptic_functions_and_integrals">user guide</a>
  59.  * for a thorough explanation.
  60.  * </em>
  61.  * </p>
  62.  * @since 2.0
  63.  */
  64. public class CarlsonEllipticIntegral {

  66.     /** Private constructor for a utility class.
  67.      */
  68.     private CarlsonEllipticIntegral() {
  69.     }

  71.     /** Compute Carlson elliptic integral R<sub>C</sub>.
  72.      * <p>
  73.      * The Carlson elliptic integral R<sub>C</sub>is defined as
  74.      * \[
  75.      *   R_C(x,y,z)=R_F(x,y,y)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}(t+y)}
  76.      * \]
  77.      * </p>
  78.      * @param x first symmetric variable of the integral
  79.      * @param y second symmetric variable of the integral
  80.      * @return Carlson elliptic integral R<sub>C</sub>
  81.      */
  82.     public static double rC(final double x, final double y) {
  83.         if (y < 0) {
  84.             // y is on the branch cut, we must use a transformation to get the Cauchy principal value
  85.             // see equation 2.14 in Carlson[1995]
  86.             final double xMy = x - y;
  87.             return FastMath.sqrt(x / xMy) * new RcRealDuplication(xMy, -y).integral();
  88.         } else {
  89.             return new RcRealDuplication(x, y).integral();
  90.         }
  91.     }

  93.     /** Compute Carlson elliptic integral R<sub>C</sub>.
  94.      * <p>
  95.      * The Carlson elliptic integral R<sub>C</sub>is defined as
  96.      * \[
  97.      *   R_C(x,y,z)=R_F(x,y,y)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}(t+y)}
  98.      * \]
  99.      * </p>
  100.      * @param x first symmetric variable of the integral
  101.      * @param y second symmetric variable of the integral
  102.      * @param <T> type of the field elements
  103.      * @return Carlson elliptic integral R<sub>C</sub>
  104.      */
  105.     public static <T extends CalculusFieldElement<T>> T rC(final T x, final T y) {
  106.         if (y.getReal() < 0) {
  107.             // y is on the branch cut, we must use a transformation to get the Cauchy principal value
  108.             // see equation 2.14 in Carlson[1995]
  109.             final T xMy = x.subtract(y);
  110.             return FastMath.sqrt(x.divide(xMy)).multiply(new RcFieldDuplication<>(xMy, y.negate()).integral());
  111.         } else {
  112.             return new RcFieldDuplication<>(x, y).integral();
  113.         }
  114.     }

  116.     /** Compute Carlson elliptic integral R<sub>C</sub>.
  117.      * <p>
  118.      * The Carlson elliptic integral R<sub>C</sub>is defined as
  119.      * \[
  120.      *   R_C(x,y,z)=R_F(x,y,y)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}(t+y)}
  121.      * \]
  122.      * </p>
  123.      * @param x first symmetric variable of the integral
  124.      * @param y second symmetric variable of the integral
  125.      * @return Carlson elliptic integral R<sub>C</sub>
  126.      */
  127.     public static Complex rC(final Complex x, final Complex y) {
  128.         if (y.getImaginaryPart() == 0 && y.getRealPart() < 0) {
  129.             // y is on the branch cut, we must use a transformation to get the Cauchy principal value
  130.             // see equation 2.14 in Carlson[1995]
  131.             final Complex xMy = x.subtract(y);
  132.             return FastMath.sqrt(x.divide(xMy)).multiply(new RcFieldDuplication<>(xMy, y.negate()).integral());
  133.         } else {
  134.             return new RcFieldDuplication<>(x, y).integral();
  135.         }
  136.     }

  138.     /** Compute Carlson elliptic integral R<sub>C</sub>.
  139.      * <p>
  140.      * The Carlson elliptic integral R<sub>C</sub>is defined as
  141.      * \[
  142.      *   R_C(x,y,z)=R_F(x,y,y)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}(t+y)}
  143.      * \]
  144.      * </p>
  145.      * @param x first symmetric variable of the integral
  146.      * @param y second symmetric variable of the integral
  147.      * @param <T> type of the field elements
  148.      * @return Carlson elliptic integral R<sub>C</sub>
  149.      */
  150.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> rC(final FieldComplex<T> x, final FieldComplex<T> y) {
  151.         if (y.getImaginaryPart().isZero() && y.getRealPart().getReal() < 0) {
  152.             // y is on the branch cut, we must use a transformation to get the Cauchy principal value
  153.             // see equation 2.14 in Carlson[1995]
  154.             final FieldComplex<T> xMy = x.subtract(y);
  155.             return FastMath.sqrt(x.divide(xMy)).multiply(new RcFieldDuplication<>(xMy, y.negate()).integral());
  156.         } else {
  157.             return new RcFieldDuplication<>(x, y).integral();
  158.         }
  159.     }

  161.     /** Compute Carlson elliptic integral R<sub>F</sub>.
  162.      * <p>
  163.      * The Carlson elliptic integral R<sub>F</sub> is defined as
  164.      * \[
  165.      *   R_F(x,y,z)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}}
  166.      * \]
  167.      * </p>
  168.      * @param x first symmetric variable of the integral
  169.      * @param y second symmetric variable of the integral
  170.      * @param z third symmetric variable of the integral
  171.      * @return Carlson elliptic integral R<sub>F</sub>
  172.      */
  173.     public static double rF(final double x, final double y, final double z) {
  174.         return new RfRealDuplication(x, y, z).integral();
  175.     }

  177.     /** Compute Carlson elliptic integral R<sub>F</sub>.
  178.      * <p>
  179.      * The Carlson elliptic integral R<sub>F</sub> is defined as
  180.      * \[
  181.      *   R_F(x,y,z)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}}
  182.      * \]
  183.      * </p>
  184.      * @param x first symmetric variable of the integral
  185.      * @param y second symmetric variable of the integral
  186.      * @param z third symmetric variable of the integral
  187.      * @param <T> type of the field elements
  188.      * @return Carlson elliptic integral R<sub>F</sub>
  189.      */
  190.     public static <T extends CalculusFieldElement<T>> T rF(final T x, final T y, final T z) {
  191.         return new RfFieldDuplication<>(x, y, z).integral();
  192.     }

  194.     /** Compute Carlson elliptic integral R<sub>F</sub>.
  195.      * <p>
  196.      * The Carlson elliptic integral R<sub>F</sub> is defined as
  197.      * \[
  198.      *   R_F(x,y,z)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}}
  199.      * \]
  200.      * </p>
  201.      * @param x first symmetric variable of the integral
  202.      * @param y second symmetric variable of the integral
  203.      * @param z third symmetric variable of the integral
  204.      * @return Carlson elliptic integral R<sub>F</sub>
  205.      */
  206.     public static Complex rF(final Complex x, final Complex y, final Complex z) {
  207.         return new RfFieldDuplication<>(x, y, z).integral();
  208.     }

  210.     /** Compute Carlson elliptic integral R<sub>F</sub>.
  211.      * <p>
  212.      * The Carlson elliptic integral R<sub>F</sub> is defined as
  213.      * \[
  214.      *   R_F(x,y,z)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}}
  215.      * \]
  216.      * </p>
  217.      * @param x first symmetric variable of the integral
  218.      * @param y second symmetric variable of the integral
  219.      * @param z third symmetric variable of the integral
  220.      * @param <T> type of the field elements
  221.      * @return Carlson elliptic integral R<sub>F</sub>
  222.      */
  223.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> rF(final FieldComplex<T> x, final FieldComplex<T> y, final FieldComplex<T> z) {
  224.         return new RfFieldDuplication<>(x, y, z).integral();
  225.     }

  227.     /** Compute Carlson elliptic integral R<sub>J</sub>.
  228.      * <p>
  229.      * The Carlson elliptic integral R<sub>J</sub> is defined as
  230.      * \[
  231.      *   R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
  232.      * \]
  233.      * </p>
  234.      * @param x first symmetric variable of the integral
  235.      * @param y second symmetric variable of the integral
  236.      * @param z third symmetric variable of the integral
  237.      * @param p fourth <em>not</em> symmetric variable of the integral
  238.      * @return Carlson elliptic integral R<sub>J</sub>
  239.      */
  240.     public static double rJ(final double x, final double y, final double z, final double p) {
  241.         final double delta = (p - x) * (p - y) * (p - z);
  242.         return rJ(x, y, z, p, delta);
  243.     }

  245.     /** Compute Carlson elliptic integral R<sub>J</sub>.
  246.      * <p>
  247.      * The Carlson elliptic integral R<sub>J</sub> is defined as
  248.      * \[
  249.      *   R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
  250.      * \]
  251.      * </p>
  252.      * @param x first symmetric variable of the integral
  253.      * @param y second symmetric variable of the integral
  254.      * @param z third symmetric variable of the integral
  255.      * @param p fourth <em>not</em> symmetric variable of the integral
  256.      * @param delta precomputed value of (p-x)(p-y)(p-z)
  257.      * @return Carlson elliptic integral R<sub>J</sub>
  258.      */
  259.     public static double rJ(final double x, final double y, final double z, final double p, final double delta) {
  260.         return new RjRealDuplication(x, y, z, p, delta).integral();
  261.     }

  263.     /** Compute Carlson elliptic integral R<sub>J</sub>.
  264.      * <p>
  265.      * The Carlson elliptic integral R<sub>J</sub> is defined as
  266.      * \[
  267.      *   R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
  268.      * \]
  269.      * </p>
  270.      * @param x first symmetric variable of the integral
  271.      * @param y second symmetric variable of the integral
  272.      * @param z third symmetric variable of the integral
  273.      * @param p fourth <em>not</em> symmetric variable of the integral
  274.      * @param <T> type of the field elements
  275.      * @return Carlson elliptic integral R<sub>J</sub>
  276.      */
  277.     public static <T extends CalculusFieldElement<T>> T rJ(final T x, final T y, final T z, final T p) {
  278.         final T delta = p.subtract(x).multiply(p.subtract(y)).multiply(p.subtract(z));
  279.         return new RjFieldDuplication<>(x, y, z, p, delta). integral();
  280.     }

  282.     /** Compute Carlson elliptic integral R<sub>J</sub>.
  283.      * <p>
  284.      * The Carlson elliptic integral R<sub>J</sub> is defined as
  285.      * \[
  286.      *   R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
  287.      * \]
  288.      * </p>
  289.      * @param x first symmetric variable of the integral
  290.      * @param y second symmetric variable of the integral
  291.      * @param z third symmetric variable of the integral
  292.      * @param p fourth <em>not</em> symmetric variable of the integral
  293.      * @param delta precomputed value of (p-x)(p-y)(p-z)
  294.      * @param <T> type of the field elements
  295.      * @return Carlson elliptic integral R<sub>J</sub>
  296.      */
  297.     public static <T extends CalculusFieldElement<T>> T rJ(final T x, final T y, final T z, final T p, final T delta) {
  298.         return new RjFieldDuplication<>(x, y, z, p, delta).integral();
  299.     }

  301.     /** Compute Carlson elliptic integral R<sub>J</sub>.
  302.      * <p>
  303.      * The Carlson elliptic integral R<sub>J</sub> is defined as
  304.      * \[
  305.      *   R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
  306.      * \]
  307.      * </p>
  308.      * @param x first symmetric variable of the integral
  309.      * @param y second symmetric variable of the integral
  310.      * @param z third symmetric variable of the integral
  311.      * @param p fourth <em>not</em> symmetric variable of the integral
  312.      * @return Carlson elliptic integral R<sub>J</sub>
  313.      */
  314.     public static Complex rJ(final Complex x, final Complex y, final Complex z, final Complex p) {
  315.         final Complex delta = p.subtract(x).multiply(p.subtract(y)).multiply(p.subtract(z));
  316.         return new RjFieldDuplication<>(x, y, z, p, delta).integral();
  317.     }

  319.     /** Compute Carlson elliptic integral R<sub>J</sub>.
  320.      * <p>
  321.      * The Carlson elliptic integral R<sub>J</sub> is defined as
  322.      * \[
  323.      *   R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
  324.      * \]
  325.      * </p>
  326.      * @param x first symmetric variable of the integral
  327.      * @param y second symmetric variable of the integral
  328.      * @param z third symmetric variable of the integral
  329.      * @param p fourth <em>not</em> symmetric variable of the integral
  330.      * @param delta precomputed value of (p-x)(p-y)(p-z)
  331.      * @return Carlson elliptic integral R<sub>J</sub>
  332.      */
  333.     public static Complex rJ(final Complex x, final Complex y, final Complex z, final Complex p, final Complex delta) {
  334.         return new RjFieldDuplication<>(x, y, z, p, delta).integral();
  335.     }

  337.     /** Compute Carlson elliptic integral R<sub>J</sub>.
  338.      * <p>
  339.      * The Carlson elliptic integral R<sub>J</sub> is defined as
  340.      * \[
  341.      *   R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
  342.      * \]
  343.      * </p>
  344.      * @param x first symmetric variable of the integral
  345.      * @param y second symmetric variable of the integral
  346.      * @param z third symmetric variable of the integral
  347.      * @param p fourth <em>not</em> symmetric variable of the integral
  348.      * @param <T> type of the field elements
  349.      * @return Carlson elliptic integral R<sub>J</sub>
  350.      */
  351.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> rJ(final FieldComplex<T> x, final FieldComplex<T> y,
  352.                                                                          final FieldComplex<T> z, final FieldComplex<T> p) {
  353.         final FieldComplex<T> delta = p.subtract(x).multiply(p.subtract(y)).multiply(p.subtract(z));
  354.         return new RjFieldDuplication<>(x, y, z, p, delta).integral();
  355.     }

  357.     /** Compute Carlson elliptic integral R<sub>J</sub>.
  358.      * <p>
  359.      * The Carlson elliptic integral R<sub>J</sub> is defined as
  360.      * \[
  361.      *   R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
  362.      * \]
  363.      * </p>
  364.      * @param x first symmetric variable of the integral
  365.      * @param y second symmetric variable of the integral
  366.      * @param z third symmetric variable of the integral
  367.      * @param p fourth <em>not</em> symmetric variable of the integral
  368.      * @param delta precomputed value of (p-x)(p-y)(p-z)
  369.      * @param <T> type of the field elements
  370.      * @return Carlson elliptic integral R<sub>J</sub>
  371.      */
  372.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> rJ(final FieldComplex<T> x, final FieldComplex<T> y,
  373.                                                                          final FieldComplex<T> z, final FieldComplex<T> p,
  374.                                                                          final FieldComplex<T> delta) {
  375.         return new RjFieldDuplication<>(x, y, z, p, delta).integral();
  376.     }

  378.     /** Compute Carlson elliptic integral R<sub>D</sub>.
  379.      * <p>
  380.      * The Carlson elliptic integral R<sub>D</sub> is defined as
  381.      * \[
  382.      *   R_D(x,y,z)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+z)}
  383.      * \]
  384.      * </p>
  385.      * @param x first symmetric variable of the integral
  386.      * @param y second symmetric variable of the integral
  387.      * @param z third symmetric variable of the integral
  388.      * @return Carlson elliptic integral R<sub>D</sub>
  389.      */
  390.     public static double rD(final double x, final double y, final double z) {
  391.         return new RdRealDuplication(x, y, z).integral();
  392.     }

  394.     /** Compute Carlson elliptic integral R<sub>D</sub>.
  395.      * <p>
  396.      * The Carlson elliptic integral R<sub>D</sub> is defined as
  397.      * \[
  398.      *   R_D(x,y,z)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+z)}
  399.      * \]
  400.      * </p>
  401.      * @param x first symmetric variable of the integral
  402.      * @param y second symmetric variable of the integral
  403.      * @param z third symmetric variable of the integral
  404.      * @param <T> type of the field elements
  405.      * @return Carlson elliptic integral R<sub>D</sub>
  406.      */
  407.     public static <T extends CalculusFieldElement<T>> T rD(final T x, final T y, final T z) {
  408.         return new RdFieldDuplication<>(x, y, z).integral();
  409.     }

  411.     /** Compute Carlson elliptic integral R<sub>D</sub>.
  412.      * <p>
  413.      * The Carlson elliptic integral R<sub>D</sub> is defined as
  414.      * \[
  415.      *   R_D(x,y,z)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+z)}
  416.      * \]
  417.      * </p>
  418.      * @param x first symmetric variable of the integral
  419.      * @param y second symmetric variable of the integral
  420.      * @param z third symmetric variable of the integral
  421.      * @return Carlson elliptic integral R<sub>D</sub>
  422.      */
  423.     public static Complex rD(final Complex x, final Complex y, final Complex z) {
  424.         return new RdFieldDuplication<>(x, y, z).integral();
  425.     }

  427.     /** Compute Carlson elliptic integral R<sub>D</sub>.
  428.      * <p>
  429.      * The Carlson elliptic integral R<sub>D</sub> is defined as
  430.      * \[
  431.      *   R_D(x,y,z)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+z)}
  432.      * \]
  433.      * </p>
  434.      * @param x first symmetric variable of the integral
  435.      * @param y second symmetric variable of the integral
  436.      * @param z third symmetric variable of the integral
  437.      * @param <T> type of the field elements
  438.      * @return Carlson elliptic integral R<sub>D</sub>
  439.      */
  440.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> rD(final FieldComplex<T> x, final FieldComplex<T> y,
  441.                                                                          final FieldComplex<T> z) {
  442.         return new RdFieldDuplication<>(x, y, z).integral();
  443.     }

  445.     /** Compute Carlson elliptic integral R<sub>G</sub>.
  446.      * <p>
  447.      * The Carlson elliptic integral R<sub>G</sub>is defined as
  448.      * \[
  449.      *   R_{G}(x,y,z)=\frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
  450.      *                \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t
  451.      * \]
  452.      * </p>
  453.      * @param x first symmetric variable of the integral
  454.      * @param y second symmetric variable of the integral
  455.      * @param z second symmetric variable of the integral
  456.      * @return Carlson elliptic integral R<sub>G</sub>
  457.      */
  458.     public static double rG(final double x, final double y, final double z) {
  459.         return generalComputeRg(x, y, z);
  460.     }

  462.     /** Compute Carlson elliptic integral R<sub>G</sub>.
  463.      * <p>
  464.      * The Carlson elliptic integral R<sub>G</sub>is defined as
  465.      * \[
  466.      *   R_{G}(x,y,z)=\frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
  467.      *                \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t
  468.      * \]
  469.      * </p>
  470.      * @param x first symmetric variable of the integral
  471.      * @param y second symmetric variable of the integral
  472.      * @param z second symmetric variable of the integral
  473.      * @param <T> type of the field elements
  474.      * @return Carlson elliptic integral R<sub>G</sub>
  475.      */
  476.     public static <T extends CalculusFieldElement<T>> T rG(final T x, final T y, final T z) {
  477.         return generalComputeRg(x, y, z);
  478.     }

  480.     /** Compute Carlson elliptic integral R<sub>G</sub>.
  481.      * <p>
  482.      * The Carlson elliptic integral R<sub>G</sub>is defined as
  483.      * \[
  484.      *   R_{G}(x,y,z)=\frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
  485.      *                \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t
  486.      * \]
  487.      * </p>
  488.      * @param x first symmetric variable of the integral
  489.      * @param y second symmetric variable of the integral
  490.      * @param z second symmetric variable of the integral
  491.      * @return Carlson elliptic integral R<sub>G</sub>
  492.      */
  493.     public static Complex rG(final Complex x, final Complex y, final Complex z) {
  494.         return generalComputeRg(x, y, z);
  495.     }

  497.     /** Compute Carlson elliptic integral R<sub>G</sub>.
  498.      * <p>
  499.      * The Carlson elliptic integral R<sub>G</sub>is defined as
  500.      * \[
  501.      *   R_{G}(x,y,z)=\frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
  502.      *                \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t
  503.      * \]
  504.      * </p>
  505.      * @param x first symmetric variable of the integral
  506.      * @param y second symmetric variable of the integral
  507.      * @param z second symmetric variable of the integral
  508.      * @param <T> type of the field elements
  509.      * @return Carlson elliptic integral R<sub>G</sub>
  510.      */
  511.     public static <T extends CalculusFieldElement<T>> FieldComplex<T> rG(final FieldComplex<T> x,
  512.                                                                          final FieldComplex<T> y,
  513.                                                                          final FieldComplex<T> z) {
  514.         return generalComputeRg(x, y, z);
  515.     }

  517.     /** Compute Carlson elliptic integral R<sub>G</sub> in the general case.
  518.      * @param x first symmetric variable of the integral
  519.      * @param y second symmetric variable of the integral
  520.      * @param z third symmetric variable of the integral
  521.      * @return Carlson elliptic integral R<sub>G</sub>
  522.      */
  523.     private static double generalComputeRg(final double x, final double y, final double z) {
  524.         // permute parameters if needed to avoid cancellations
  525.         if (x <= y) {
  526.             if (y <= z) {
  527.                 // x ≤ y ≤ z
  528.                 return permutedComputeRg(x, z, y);
  529.             } else if (x <= z) {
  530.                 // x ≤ z < y
  531.                 return permutedComputeRg(x, y, z);
  532.             } else {
  533.                 // z < x ≤ y
  534.                 return permutedComputeRg(z, y, x);
  535.             }
  536.         } else if (x <= z) {
  537.             // y < x ≤ z
  538.             return permutedComputeRg(y, z, x);
  539.         } else if (y <= z) {
  540.             // y ≤ z < x
  541.             return permutedComputeRg(y, x, z);
  542.         } else {
  543.             // z < y < x
  544.             return permutedComputeRg(z, x, y);
  545.         }
  546.     }

  548.     /** Compute Carlson elliptic integral R<sub>G</sub> in the general case.
  549.      * @param x first symmetric variable of the integral
  550.      * @param y second symmetric variable of the integral
  551.      * @param z third symmetric variable of the integral
  552.      * @param <T> type of the field elements (really {@link Complex} or {@link FieldComplex})
  553.      * @return Carlson elliptic integral R<sub>G</sub>
  554.      */
  555.     private static <T extends CalculusFieldElement<T>> T generalComputeRg(final T x, final T y, final T z) {
  556.         // permute parameters if needed to avoid cancellations
  557.         final double xR = x.getReal();
  558.         final double yR = y.getReal();
  559.         final double zR = z.getReal();
  560.         if (xR <= yR) {
  561.             if (yR <= zR) {
  562.                 // x ≤ y ≤ z
  563.                 return permutedComputeRg(x, z, y);
  564.             } else if (xR <= zR) {
  565.                 // x ≤ z < y
  566.                 return permutedComputeRg(x, y, z);
  567.             } else {
  568.                 // z < x ≤ y
  569.                 return permutedComputeRg(z, y, x);
  570.             }
  571.         } else if (xR <= zR) {
  572.             // y < x ≤ z
  573.             return permutedComputeRg(y, z, x);
  574.         } else if (yR <= zR) {
  575.             // y ≤ z < x
  576.             return permutedComputeRg(y, x, z);
  577.         } else {
  578.             // z < y < x
  579.             return permutedComputeRg(z, x, y);
  580.         }
  581.     }

  583.     /** Compute Carlson elliptic integral R<sub>G</sub> with already permuted variables to avoid cancellations.
  584.      * @param x first symmetric variable of the integral
  585.      * @param y second symmetric variable of the integral
  586.      * @param z third symmetric variable of the integral
  587.      * @return Carlson elliptic integral R<sub>G</sub>
  588.      */
  589.     private static double permutedComputeRg(final double x, final double y, final double z) {
  590.         // permute parameters if needed to avoid divisions by zero
  591.         if (z == 0) {
  592.             return x == 0 ? safeComputeRg(z, x, y) : safeComputeRg(y, z, x);
  593.         } else {
  594.             return safeComputeRg(x, y, z);
  595.         }
  596.     }

  598.     /** Compute Carlson elliptic integral R<sub>G</sub> with already permuted variables to avoid cancellations.
  599.      * @param x first symmetric variable of the integral
  600.      * @param y second symmetric variable of the integral
  601.      * @param z third symmetric variable of the integral
  602.      * @param <T> type of the field elements (really {@link Complex} or {@link FieldComplex})
  603.      * @return Carlson elliptic integral R<sub>G</sub>
  604.      */
  605.     private static <T extends CalculusFieldElement<T>> T permutedComputeRg(final T x, final T y, final T z) {
  606.         // permute parameters if needed to avoid divisions by zero
  607.         if (z.isZero()) {
  608.             return x.isZero() ? safeComputeRg(z, x, y) : safeComputeRg(y, z, x);
  609.         } else {
  610.             return safeComputeRg(x, y, z);
  611.         }
  612.     }

  614.     /** Compute Carlson elliptic integral R<sub>G</sub> with non-zero third variable.
  615.      * @param x first symmetric variable of the integral
  616.      * @param y second symmetric variable of the integral
  617.      * @param z third symmetric variable of the integral
  618.      * @see <a href="https://dlmf.nist.gov/19.21#E10">Digital Library of Mathematical Functions, equation 19.21.10</a>
  619.      * @return Carlson elliptic integral R<sub>G</sub>
  620.      */
  621.     private static double safeComputeRg(final double x, final double y, final double z) {

  623.         // contribution of the R_F integral
  624.         final double termF = new RfRealDuplication(x, y, z).integral() * z;

  626.         // contribution of the R_D integral
  627.         final double termD = (x - z) * (y - z) * new RdRealDuplication(x, y, z).integral() / 3;

  629.         // contribution of the square roots
  630.         final double termS = FastMath.sqrt(x * y / z);

  632.         // equation 19.21.10
  633.         return (termF - termD + termS) * 0.5;

  635.     }

  637.     /** Compute Carlson elliptic integral R<sub>G</sub> with non-zero third variable.
  638.      * @param x first symmetric variable of the integral
  639.      * @param y second symmetric variable of the integral
  640.      * @param z third symmetric variable of the integral
  641.      * @param <T> type of the field elements (really {@link Complex} or {@link FieldComplex})
  642.      * @see <a href="https://dlmf.nist.gov/19.21#E10">Digital Library of Mathematical Functions, equation 19.21.10</a>
  643.      * @return Carlson elliptic integral R<sub>G</sub>
  644.      */
  645.     private static <T extends CalculusFieldElement<T>> T safeComputeRg(final T x, final T y, final T z) {

  647.         // contribution of the R_F integral
  648.         final T termF = new RfFieldDuplication<>(x, y, z).integral().multiply(z);

  650.         // contribution of the R_D integral
  651.         final T termD = x.subtract(z).multiply(y.subtract(z)).multiply(new RdFieldDuplication<>(x, y, z).integral()).divide(3);

  653.         // contribution of the square roots
  654.         // BEWARE: this term MUST be computed as √x√y/√z with all square roots selected with positive real part
  655.         // and NOT as √(xy/z), otherwise sign errors may occur
  656.         final T termS = x.sqrt().multiply(y.sqrt()).divide(z.sqrt());

  658.         // equation 19.21.10
  659.         return termF.subtract(termD).add(termS).multiply(0.5);

  661.     }

  663. }
Created with JaCoCo 0.8.12.202403310830