/*
    * Licensed to the Hipparchus project under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      https://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.hipparchus.special.elliptic.legendre;
  
  import java.util.function.DoubleFunction;
  import java.util.function.Function;
  
  import org.hipparchus.CalculusFieldElement;
  import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
  import org.hipparchus.complex.Complex;
  import org.hipparchus.complex.ComplexUnivariateIntegrator;
  import org.hipparchus.complex.FieldComplex;
  import org.hipparchus.complex.FieldComplexUnivariateIntegrator;
  import org.hipparchus.special.elliptic.carlson.CarlsonEllipticIntegral;
  import org.hipparchus.util.FastMath;
  import org.hipparchus.util.MathUtils;
  
  /** Complete and incomplete elliptic integrals in Legendre form.
   * <p>
   * The elliptic integrals are related to Jacobi elliptic functions.
   * </p>
   * <p>
   * <em>
   * Beware that when computing elliptic integrals in the complex plane,
   * many issues arise due to branch cuts. See the
   * <a href="https://www.hipparchus.org/hipparchus-core/special.html#Elliptic_functions_and_integrals">user guide</a>
   * for a thorough explanation.
   * </em>
   * </p>
   * <p>
   * There are different conventions to interpret the arguments of
   * Legendre elliptic integrals. In mathematical texts, these conventions show
   * up using the separator between arguments. So for example for the incomplete
   * integral of the first kind F we have:
   * </p>
   * <ul>
   *   <li>F(φ, k): the first argument φ is an angle and the second argument k
   *       is the elliptic modulus: this is the trigonometric form of the integral</li>
   *   <li>F(φ; m): the first argument φ is an angle and the second argument m=k²
   *       is the parameter: this is also a trigonometric form of the integral</li>
   *   <li>F(x|m): the first argument x=sin(φ) is not an angle anymore and the
   *       second argument m=k² is the parameter: this is the Legendre form</li>
   *   <li>F(φ\α): the first argument φ is an angle and the second argument α is the
   *       modular angle</li>
   * </ul>
   * <p>
   * As we have no separator in a method call, we have to adopt one convention
   * and stick to it. In Hipparchus, we adopted the Legendre form (i.e. F(x|m),
   * with x=sin(φ) and m=k². These conventions are consistent with Wolfram Alpha
   * functions EllipticF, EllipticE, ElliptiPI…
   * </p>
   * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
   * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
   * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
   * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
   * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
   * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
   * @since 2.0
   */
  public class LegendreEllipticIntegral {
  
      /** Private constructor for a utility class.
       */
      private LegendreEllipticIntegral() {
          // nothing to do
      }
  
      /** Get the nome q.
       * @param m parameter (m=k² where k is the elliptic modulus)
       * @return nome q
       */
      public static double nome(final double m) {
          if (m < 1.0e-16) {
              // first terms of infinite series in Abramowitz and Stegun 17.3.21
              final double m16 = m * 0.0625;
              return m16 * (1 + 8 * m16);
          } else {
              return FastMath.exp(-FastMath.PI * bigKPrime(m) / bigK(m));
          }
      }
  
      /** Get the nome q.
       * @param m parameter (m=k² where k is the elliptic modulus)
       * @param <T> the type of the field elements
       * @return nome q
      */
     public static <T extends CalculusFieldElement<T>> T nome(final T m) {
         final T one = m.getField().getOne();
         if (m.norm() < 100 * one.ulp().getReal()) {
             // first terms of infinite series in Abramowitz and Stegun 17.3.21
             final T m16 = m.multiply(0.0625);
             return m16.multiply(m16.multiply(8).add(1));
         } else {
             return FastMath.exp(bigKPrime(m).divide(bigK(m)).multiply(one.getPi().negate()));
         }
     }
 
     /** Get the complete elliptic integral of the first kind K(m).
      * <p>
      * The complete elliptic integral of the first kind K(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * it corresponds to the real quarter-period of Jacobi elliptic functions
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral of the first kind K(m)
      * @see #bigKPrime(double)
      * @see #bigF(double, double)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static double bigK(final double m) {
         if (m < 1.0e-8) {
             // first terms of infinite series in Abramowitz and Stegun 17.3.11
             return (1 + 0.25 * m) * MathUtils.SEMI_PI;
         } else {
             return CarlsonEllipticIntegral.rF(0, 1.0 - m, 1);
         }
     }
 
     /** Get the complete elliptic integral of the first kind K(m).
      * <p>
      * The complete elliptic integral of the first kind K(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * it corresponds to the real quarter-period of Jacobi elliptic functions
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral of the first kind K(m)
      * @see #bigKPrime(CalculusFieldElement)
      * @see #bigF(CalculusFieldElement, CalculusFieldElement)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> T bigK(final T m) {
         final T zero = m.getField().getZero();
         final T one  = m.getField().getOne();
         if (m.norm() < 1.0e7 * one.ulp().getReal()) {
 
             // first terms of infinite series in Abramowitz and Stegun 17.3.11
             return one.add(m.multiply(0.25)).multiply(zero.getPi().multiply(0.5));
 
         } else {
             return CarlsonEllipticIntegral.rF(zero, one.subtract(m), one);
         }
     }
 
     /** Get the complete elliptic integral of the first kind K(m).
      * <p>
      * The complete elliptic integral of the first kind K(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * it corresponds to the real quarter-period of Jacobi elliptic functions
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral of the first kind K(m)
      * @see #bigKPrime(Complex)
      * @see #bigF(Complex, Complex)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigK(final Complex m) {
         if (m.norm() < 1.0e-8) {
             // first terms of infinite series in Abramowitz and Stegun 17.3.11
             return Complex.ONE.add(m.multiply(0.25)).multiply(MathUtils.SEMI_PI);
         } else {
             return CarlsonEllipticIntegral.rF(Complex.ZERO, Complex.ONE.subtract(m), Complex.ONE);
         }
     }
 
     /** Get the complete elliptic integral of the first kind K(m).
      * <p>
      * The complete elliptic integral of the first kind K(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * it corresponds to the real quarter-period of Jacobi elliptic functions
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral of the first kind K(m)
      * @see #bigKPrime(FieldComplex)
      * @see #bigF(FieldComplex, FieldComplex)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigK(final FieldComplex<T> m) {
         final FieldComplex<T> zero = m.getField().getZero();
         final FieldComplex<T> one  = m.getField().getOne();
         if (m.norm() < 1.0e7 * one.ulp().getReal()) {
 
             // first terms of infinite series in Abramowitz and Stegun 17.3.11
             return one.add(m.multiply(0.25)).multiply(zero.getPi().multiply(0.5));
 
         } else {
             return CarlsonEllipticIntegral.rF(zero, one.subtract(m), one);
         }
     }
 
     /** Get the complete elliptic integral of the first kind K'(m).
      * <p>
      * The complete elliptic integral of the first kind K'(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-(1-m) \sin^2\theta}}
      * \]
      * it corresponds to the imaginary quarter-period of Jacobi elliptic functions
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral of the first kind K'(m)
      * @see #bigK(double)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static double bigKPrime(final double m) {
         return CarlsonEllipticIntegral.rF(0, m, 1);
     }
 
     /** Get the complete elliptic integral of the first kind K'(m).
      * <p>
      * The complete elliptic integral of the first kind K'(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-(1-m) \sin^2\theta}}
      * \]
      * it corresponds to the imaginary quarter-period of Jacobi elliptic functions
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral of the first kind K'(m)
      * @see #bigK(CalculusFieldElement)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> T bigKPrime(final T m) {
         final T zero = m.getField().getZero();
         final T one  = m.getField().getOne();
         return CarlsonEllipticIntegral.rF(zero, m, one);
     }
 
     /** Get the complete elliptic integral of the first kind K'(m).
      * <p>
      * The complete elliptic integral of the first kind K'(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-(1-m) \sin^2\theta}}
      * \]
      * it corresponds to the imaginary quarter-period of Jacobi elliptic functions
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral of the first kind K'(m)
      * @see #bigK(Complex)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigKPrime(final Complex m) {
         return CarlsonEllipticIntegral.rF(Complex.ZERO, m, Complex.ONE);
     }
 
     /** Get the complete elliptic integral of the first kind K'(m).
      * <p>
      * The complete elliptic integral of the first kind K'(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-(1-m) \sin^2\theta}}
      * \]
      * it corresponds to the imaginary quarter-period of Jacobi elliptic functions
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral of the first kind K'(m)
      * @see #bigK(FieldComplex)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html">Complete Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigKPrime(final FieldComplex<T> m) {
         final FieldComplex<T> zero = m.getField().getZero();
         final FieldComplex<T> one  = m.getField().getOne();
         return CarlsonEllipticIntegral.rF(zero, m, one);
     }
 
     /** Get the complete elliptic integral of the second kind E(m).
      * <p>
      * The complete elliptic integral of the second kind E(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral of the second kind E(m)
      * @see #bigE(double, double)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static double bigE(final double m) {
         return CarlsonEllipticIntegral.rG(0, 1 - m, 1) * 2;
     }
 
     /** Get the complete elliptic integral of the second kind E(m).
      * <p>
      * The complete elliptic integral of the second kind E(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral of the second kind E(m)
      * @see #bigE(CalculusFieldElement, CalculusFieldElement)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> T bigE(final T m) {
         final T zero = m.getField().getZero();
         final T one  = m.getField().getOne();
         return CarlsonEllipticIntegral.rG(zero, one.subtract(m), one).multiply(2);
     }
 
     /** Get the complete elliptic integral of the second kind E(m).
      * <p>
      * The complete elliptic integral of the second kind E(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral of the second kind E(m)
      * @see #bigE(Complex, Complex)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigE(final Complex m) {
         return CarlsonEllipticIntegral.rG(Complex.ZERO,
                                           Complex.ONE.subtract(m),
                                           Complex.ONE).multiply(2);
     }
 
     /** Get the complete elliptic integral of the second kind E(m).
      * <p>
      * The complete elliptic integral of the second kind E(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral of the second kind E(m)
      * @see #bigE(FieldComplex, FieldComplex)
      * @see <a href="https://mathworld.wolfram.com/CompleteEllipticIntegraloftheSecondKind.html">Complete Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigE(final FieldComplex<T> m) {
         final FieldComplex<T> zero = m.getField().getZero();
         final FieldComplex<T> one  = m.getField().getOne();
         return CarlsonEllipticIntegral.rG(zero, one.subtract(m), one).multiply(2);
     }
 
     /** Get the complete elliptic integral D(m) = [K(m) - E(m)]/m.
      * <p>
      * The complete elliptic integral D(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral D(m)
      * @see #bigD(double, double)
      */
     public static double bigD(final double m) {
         return CarlsonEllipticIntegral.rD(0, 1 - m, 1) / 3;
     }
 
     /** Get the complete elliptic integral D(m) = [K(m) - E(m)]/m.
      * <p>
      * The complete elliptic integral D(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral D(m)
      * @see #bigD(CalculusFieldElement, CalculusFieldElement)
      */
     public static <T extends CalculusFieldElement<T>> T bigD(final T m) {
         final T zero = m.getField().getZero();
         final T one  = m.getField().getOne();
         return CarlsonEllipticIntegral.rD(zero, one.subtract(m), one).divide(3);
     }
 
     /** Get the complete elliptic integral D(m) = [K(m) - E(m)]/m.
      * <p>
      * The complete elliptic integral D(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral D(m)
      * @see #bigD(Complex, Complex)
      */
     public static Complex bigD(final Complex m) {
         return CarlsonEllipticIntegral.rD(Complex.ZERO, Complex.ONE.subtract(m), Complex.ONE).divide(3);
     }
 
     /** Get the complete elliptic integral D(m) = [K(m) - E(m)]/m.
      * <p>
      * The complete elliptic integral D(m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral D(m)
      * @see #bigD(FieldComplex, FieldComplex)
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigD(final FieldComplex<T> m) {
         final FieldComplex<T> zero = m.getField().getZero();
         final FieldComplex<T> one  = m.getField().getOne();
         return CarlsonEllipticIntegral.rD(zero, one.subtract(m), one).divide(3);
     }
 
     /** Get the complete elliptic integral of the third kind Π(n, m).
      * <p>
      * The complete elliptic integral of the third kind Π(n, m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param n elliptic characteristic
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral of the third kind Π(n, m)
      * @see #bigPi(double, double, double)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static double bigPi(final double n, final double m) {
         final double kPrime2 = 1 - m;
         final double delta   = n * (m - n) * (n - 1);
         return CarlsonEllipticIntegral.rF(0, kPrime2, 1) +
                CarlsonEllipticIntegral.rJ(0, kPrime2, 1, 1 - n, delta) * n / 3;
     }
 
     /** Get the complete elliptic integral of the third kind Π(n, m).
      * <p>
      * The complete elliptic integral of the third kind Π(n, m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param n elliptic characteristic
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral of the third kind Π(n, m)
      * @see #bigPi(CalculusFieldElement, CalculusFieldElement, CalculusFieldElement)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> T bigPi(final T n, final T m) {
         final T zero    = m.getField().getZero();
         final T one     = m.getField().getOne();
         final T kPrime2 = one.subtract(m);
         final T delta   = n.multiply(m.subtract(n)).multiply(n.subtract(1));
         return CarlsonEllipticIntegral.rF(zero, kPrime2, one).
                add(CarlsonEllipticIntegral.rJ(zero, kPrime2, one, one.subtract(n), delta).multiply(n).divide(3));
     }
 
     /** Get the complete elliptic integral of the third kind Π(n, m).
      * <p>
      * The complete elliptic integral of the third kind Π(n, m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param n elliptic characteristic
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return complete elliptic integral of the third kind Π(n, m)
      * @see #bigPi(Complex, Complex, Complex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigPi(final Complex n, final Complex m) {
         final Complex kPrime2 = Complex.ONE.subtract(m);
         final Complex delta   = n.multiply(m.subtract(n)).multiply(n.subtract(1));
         return CarlsonEllipticIntegral.rF(Complex.ZERO, kPrime2, Complex.ONE).
                add(CarlsonEllipticIntegral.rJ(Complex.ZERO, kPrime2, Complex.ONE, Complex.ONE.subtract(n), delta).multiply(n).divide(3));
     }
 
     /** Get the complete elliptic integral of the third kind Π(n, m).
      * <p>
      * The complete elliptic integral of the third kind Π(n, m) is
      * \[
      *    \int_0^{\frac{\pi}{2}} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param n elliptic characteristic
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return complete elliptic integral of the third kind Π(n, m)
      * @see #bigPi(FieldComplex, FieldComplex, FieldComplex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigPi(final FieldComplex<T> n, final FieldComplex<T> m) {
         final FieldComplex<T> zero    = m.getField().getZero();
         final FieldComplex<T> one     = m.getField().getOne();
         final FieldComplex<T> kPrime2 = one.subtract(m);
         final FieldComplex<T> delta   = n.multiply(m.subtract(n)).multiply(n.subtract(1));
         return CarlsonEllipticIntegral.rF(zero, kPrime2, one).
                add(CarlsonEllipticIntegral.rJ(zero, kPrime2, one, one.subtract(n), delta).multiply(n).divide(3));
     }
 
     /** Get the incomplete elliptic integral of the first kind F(φ, m).
      * <p>
      * The incomplete elliptic integral of the first kind F(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return incomplete elliptic integral of the first kind F(φ, m)
      * @see #bigK(double)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static double bigF(final double phi, final double m) {
 
         // argument reduction
         final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, n -> bigK(n));
 
         // integrate part between 0 and π/2
         final double cM1 = ar.csc2 - 1.0;
         final double cMm = ar.csc2 - m;
         final double incomplete =  CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2);
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;
 
     }
 
     /** Get the incomplete elliptic integral of the first kind F(φ, m).
      * <p>
      * The incomplete elliptic integral of the first kind F(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral of the first kind F(φ, m)
      * @see #bigK(CalculusFieldElement)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> T bigF(final T phi, final T m) {
 
         // argument reduction
         final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, n -> bigK(n));
 
         // integrate part between 0 and π/2
         final T cM1        = ar.csc2.subtract(1);
         final T cMm        = ar.csc2.subtract(m);
         final T incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2);
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the first kind F(φ, m).
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the first kind F(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return incomplete elliptic integral of the first kind F(φ, m)
      * @see #bigK(Complex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigF(final Complex phi, final Complex m) {
 
         // argument reduction
         final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, n -> bigK(n));
 
         // integrate part between 0 and π/2
         final Complex cM1        = ar.csc2.subtract(1);
         final Complex cMm        = ar.csc2.subtract(m);
         final Complex incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2);
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the first kind F(φ, m) using numerical integration.
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the first kind F(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on numerical integration.
      * If integration path comes too close to a pole of the integrand, then integration will fail
      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
      * even for very large {@code maxEval}. This is normal behavior.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param integrator integrator to use
      * @param maxEval maximum number of evaluations (real and imaginary
      * parts are evaluated separately, so up to twice this number may be used)
      * @return incomplete elliptic integral of the first kind F(φ, m)
      * @see #bigK(Complex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigF(final Complex phi, final Complex m,
                                final ComplexUnivariateIntegrator integrator, final int maxEval) {
         return integrator.integrate(maxEval, new First<>(m), phi.getField().getZero(), phi);
     }
 
     /** Get the incomplete elliptic integral of the first kind F(φ, m).
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the first kind F(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral of the first kind F(φ, m)
      * @see #bigK(CalculusFieldElement)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigF(final FieldComplex<T> phi, final FieldComplex<T> m) {
 
         // argument reduction
         final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, n -> bigK(n));
 
         // integrate part between 0 and π/2
         final FieldComplex<T> cM1        = ar.csc2.subtract(1);
         final FieldComplex<T> cMm        = ar.csc2.subtract(m);
         final FieldComplex<T> incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2);
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the first kind F(φ, m).
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the first kind F(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on numerical integration.
      * If integration path comes too close to a pole of the integrand, then integration will fail
      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
      * even for very large {@code maxEval}. This is normal behavior.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param integrator integrator to use
      * @param maxEval maximum number of evaluations (real and imaginary
      * parts are evaluated separately, so up to twice this number may be used)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral of the first kind F(φ, m)
      * @see #bigK(CalculusFieldElement)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html">Elliptic Integrals of the First Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigF(final FieldComplex<T> phi, final FieldComplex<T> m,
                                                                            final FieldComplexUnivariateIntegrator<T> integrator,
                                                                            final int maxEval) {
         return integrator.integrate(maxEval, new First<>(m), phi.getField().getZero(), phi);
     }
 
     /** Get the incomplete elliptic integral of the second kind E(φ, m).
      * <p>
      * The incomplete elliptic integral of the second kind E(φ, m) is
      * \[
      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return incomplete elliptic integral of the second kind E(φ, m)
      * @see #bigE(double)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static double bigE(final double phi, final double m) {
 
         // argument reduction
         final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, n -> bigE(n));
 
         // integrate part between 0 and π/2
         final double cM1        = ar.csc2 - 1.0;
         final double cMm        = ar.csc2 - m;
         final double incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2) -
                                   CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2) * (m / 3);
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;
 
     }
 
     /** Get the incomplete elliptic integral of the second kind E(φ, m).
      * <p>
      * The incomplete elliptic integral of the second kind E(φ, m) is
      * \[
      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral of the second kind E(φ, m)
      * @see #bigE(CalculusFieldElement)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> T bigE(final T phi, final T m) {
 
         // argument reduction
         final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, n -> bigE(n));
 
         // integrate part between 0 and π/2
         final T cM1        = ar.csc2.subtract(1);
         final T cMm        = ar.csc2.subtract(m);
         final T incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
                              subtract(CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).multiply(m.divide(3)));
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the second kind E(φ, m).
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the second kind E(φ, m) is
      * \[
      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return incomplete elliptic integral of the second kind E(φ, m)
      * @see #bigE(Complex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigE(final Complex phi, final Complex m) {
 
         // argument reduction
         final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, n -> bigE(n));
 
         // integrate part between 0 and π/2
         final Complex cM1        = ar.csc2.subtract(1);
         final Complex cMm        = ar.csc2.subtract(m);
         final Complex incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
                                    subtract(CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).multiply(m.divide(3)));
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the second kind E(φ, m) using numerical integration.
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the second kind E(φ, m) is
      * \[
      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on numerical integration.
      * If integration path comes too close to a pole of the integrand, then integration will fail
      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
      * even for very large {@code maxEval}. This is normal behavior.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param integrator integrator to use
      * @param maxEval maximum number of evaluations (real and imaginary
      * parts are evaluated separately, so up to twice this number may be used)
      * @return incomplete elliptic integral of the second kind E(φ, m)
      * @see #bigE(Complex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigE(final Complex phi, final Complex m,
                                final ComplexUnivariateIntegrator integrator, final int maxEval) {
         return integrator.integrate(maxEval, new Second<>(m), phi.getField().getZero(), phi);
     }
 
     /** Get the incomplete elliptic integral of the second kind E(φ, m).
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the second kind E(φ, m) is
      * \[
      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral of the second kind E(φ, m)
      * @see #bigE(FieldComplex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigE(final FieldComplex<T> phi, final FieldComplex<T> m) {
 
         // argument reduction
         final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, n -> bigE(n));
 
         // integrate part between 0 and π/2
         final FieldComplex<T> cM1        = ar.csc2.subtract(1);
         final FieldComplex<T> cMm        = ar.csc2.subtract(m);
         final FieldComplex<T> incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
                                            subtract(CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).multiply(m.divide(3)));
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the second kind E(φ, m).
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the second kind E(φ, m) is
      * \[
      *    \int_0^{\phi} \sqrt{1-m \sin^2\theta} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on numerical integration.
      * If integration path comes too close to a pole of the integrand, then integration will fail
      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
      * even for very large {@code maxEval}. This is normal behavior.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param integrator integrator to use
      * @param maxEval maximum number of evaluations (real and imaginary
      * parts are evaluated separately, so up to twice this number may be used)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral of the second kind E(φ, m)
      * @see #bigE(FieldComplex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheSecondKind.html">Elliptic Integrals of the Second Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigE(final FieldComplex<T> phi, final FieldComplex<T> m,
                                                                            final FieldComplexUnivariateIntegrator<T> integrator,
                                                                            final int maxEval) {
         return integrator.integrate(maxEval, new Second<>(m), phi.getField().getZero(), phi);
     }
 
     /** Get the incomplete elliptic integral D(φ, m) = [F(φ, m) - E(φ, m)]/m.
      * <p>
      * The incomplete elliptic integral D(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return incomplete elliptic integral D(φ, m)
      * @see #bigD(double)
      */
     public static double bigD(final double phi, final double m) {
 
         // argument reduction
         final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, n -> bigD(n));
 
         // integrate part between 0 and π/2
         final double cM1        = ar.csc2 - 1.0;
         final double cMm        = ar.csc2 - m;
         final double incomplete = CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2) / 3;
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;
 
     }
 
     /** Get the incomplete elliptic integral D(φ, m) = [F(φ, m) - E(φ, m)]/m.
      * <p>
      * The incomplete elliptic integral D(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral D(φ, m)
      * @see #bigD(CalculusFieldElement)
      */
     public static <T extends CalculusFieldElement<T>> T bigD(final T phi, final T m) {
 
         // argument reduction
         final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, n -> bigD(n));
 
         // integrate part between 0 and π/2
         final T cM1        = ar.csc2.subtract(1);
         final T cMm        = ar.csc2.subtract(m);
         final T incomplete = CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).divide(3);
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral D(φ, m) = [F(φ, m) - E(φ, m)]/m.
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral D(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return incomplete elliptic integral D(φ, m)
      * @see #bigD(Complex)
      */
     public static Complex bigD(final Complex phi, final Complex m) {
 
         // argument reduction
         final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, n -> bigD(n));
 
         // integrate part between 0 and π/2
         final Complex cM1        = ar.csc2.subtract(1);
         final Complex cMm        = ar.csc2.subtract(m);
         final Complex incomplete = CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).divide(3);
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral D(φ, m) = [F(φ, m) - E(φ, m)]/m.
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral D(φ, m) is
      * \[
      *    \int_0^{\phi} \frac{\sin^2\theta}{\sqrt{1-m \sin^2\theta}} d\theta
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral D(φ, m)
      * @see #bigD(CalculusFieldElement)
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigD(final FieldComplex<T> phi, final FieldComplex<T> m) {
 
         // argument reduction
         final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, n -> bigD(n));
 
         // integrate part between 0 and π/2
         final FieldComplex<T> cM1        = ar.csc2.subtract(1);
         final FieldComplex<T> cMm        = ar.csc2.subtract(m);
         final FieldComplex<T> incomplete = CarlsonEllipticIntegral.rD(cM1, cMm, ar.csc2).divide(3);
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
      * <p>
      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param n elliptic characteristic
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
      * @see #bigPi(double, double)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static double bigPi(final double n, final double phi, final double m) {
 
         // argument reduction
         final DoubleArgumentReduction ar = new DoubleArgumentReduction(phi, m, parameter -> bigPi(n, parameter));
 
         // integrate part between 0 and π/2
         final double cM1        = ar.csc2 - 1.0;
         final double cMm        = ar.csc2 - m;
         final double cMn        = ar.csc2 - n;
         final double delta      = n * (m - n) * (n - 1);
         final double incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2) +
                                   CarlsonEllipticIntegral.rJ(cM1, cMm, ar.csc2, cMn, delta) * n / 3;
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete - incomplete : ar.complete + incomplete;
 
     }
 
     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
      * <p>
      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param n elliptic characteristic
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
      * @see #bigPi(CalculusFieldElement, CalculusFieldElement)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> T bigPi(final T n, final T phi, final T m) {
 
         // argument reduction
         final FieldArgumentReduction<T> ar = new FieldArgumentReduction<>(phi, m, parameter -> bigPi(n, parameter));
 
         // integrate part between 0 and π/2
         final T cM1        = ar.csc2.subtract(1);
         final T cMm        = ar.csc2.subtract(m);
         final T cMn        = ar.csc2.subtract(n);
         final T delta      = n.multiply(m.subtract(n)).multiply(n.subtract(1));
         final T incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
                              add(CarlsonEllipticIntegral.rJ(cM1, cMm, ar.csc2, cMn, delta).multiply(n).divide(3));
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param n elliptic characteristic
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
      * @see #bigPi(Complex, Complex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigPi(final Complex n, final Complex phi, final Complex m) {
 
         // argument reduction
         final FieldArgumentReduction<Complex> ar = new FieldArgumentReduction<>(phi, m, parameter -> bigPi(n, parameter));
 
         // integrate part between 0 and π/2
         final Complex cM1        = ar.csc2.subtract(1);
         final Complex cMm        = ar.csc2.subtract(m);
         final Complex cMn        = ar.csc2.subtract(n);
         final Complex delta      = n.multiply(m.subtract(n)).multiply(n.subtract(1));
         final Complex incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
                                    add(CarlsonEllipticIntegral.rJ(cM1, cMm, ar.csc2, cMn, delta).multiply(n).divide(3));
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m) using numerical integration.
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on numerical integration.
      * If integration path comes too close to a pole of the integrand, then integration will fail
      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
      * even for very large {@code maxEval}. This is normal behavior.
      * </p>
      * @param n elliptic characteristic
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param integrator integrator to use
      * @param maxEval maximum number of evaluations (real and imaginary
      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
      * @see #bigPi(Complex, Complex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static Complex bigPi(final Complex n, final Complex phi, final Complex m,
                                 final ComplexUnivariateIntegrator integrator, final int maxEval) {
          return integrator.integrate(maxEval, new Third<>(n, m), phi.getField().getZero(), phi);
     }
 
     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on {@link CarlsonEllipticIntegral
      * Carlson elliptic integrals}.
      * </p>
      * @param n elliptic characteristic
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
      * @see #bigPi(FieldComplex, FieldComplex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigPi(final FieldComplex<T> n,
                                                                             final FieldComplex<T> phi,
                                                                             final FieldComplex<T> m) {
 
         // argument reduction
         final FieldArgumentReduction<FieldComplex<T>> ar = new FieldArgumentReduction<>(phi, m, parameter -> bigPi(n, parameter));
 
         // integrate part between 0 and π/2
         final FieldComplex<T> cM1        = ar.csc2.subtract(1);
         final FieldComplex<T> cMm        = ar.csc2.subtract(m);
         final FieldComplex<T> cMn        = ar.csc2.subtract(n);
         final FieldComplex<T> delta      = n.multiply(m.subtract(n)).multiply(n.subtract(1));
         final FieldComplex<T> incomplete = CarlsonEllipticIntegral.rF(cM1, cMm, ar.csc2).
                                            add(CarlsonEllipticIntegral.rJ(cM1, cMm, ar.csc2, cMn, delta).multiply(n).divide(3));
 
         // combine complete and incomplete parts
         return ar.negate ? ar.complete.subtract(incomplete) : ar.complete.add(incomplete);
 
     }
 
     /** Get the incomplete elliptic integral of the third kind Π(n, φ, m).
      * <p>
      * <em>
      * BEWARE! Elliptic integrals for complex numbers in the incomplete case
      * are considered experimental for now, they have known issues.
      * </em>
      * </p>
      * <p>
      * The incomplete elliptic integral of the third kind Π(n, φ, m) is
      * \[
      *    \int_0^{\phi} \frac{d\theta}{\sqrt{1-m \sin^2\theta}(1-n \sin^2\theta)}
      * \]
      * </p>
      * <p>
      * The algorithm for evaluating the functions is based on numerical integration.
      * If integration path comes too close to a pole of the integrand, then integration will fail
      * with a {@link org.hipparchus.exception.MathIllegalStateException MathIllegalStateException}
      * even for very large {@code maxEval}. This is normal behavior.
      * </p>
      * @param n elliptic characteristic
      * @param phi amplitude (i.e. upper bound of the integral)
      * @param m parameter (m=k² where k is the elliptic modulus)
      * @param integrator integrator to use
      * @param maxEval maximum number of evaluations (real and imaginary
      * parts are evaluated separately, so up to twice this number may be used)
      * @param <T> the type of the field elements
      * @return incomplete elliptic integral of the third kind Π(n, φ, m)
      * @see #bigPi(FieldComplex, FieldComplex)
      * @see <a href="https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html">Elliptic Integrals of the Third Kind (MathWorld)</a>
      * @see <a href="https://en.wikipedia.org/wiki/Elliptic_integral">Elliptic Integrals (Wikipedia)</a>
      */
     public static <T extends CalculusFieldElement<T>> FieldComplex<T> bigPi(final FieldComplex<T> n,
                                                                             final FieldComplex<T> phi,
                                                                             final FieldComplex<T> m,
                                                                             final FieldComplexUnivariateIntegrator<T> integrator,
                                                                             final int maxEval) {
         return integrator.integrate(maxEval, new Third<>(n, m), phi.getField().getZero(), phi);
     }
 
     /** Argument reduction for an incomplete integral. */
     private static class DoubleArgumentReduction {
 
         /** Complete part. */
         private final double complete;
 
         /** Squared cosecant of the Jacobi amplitude. */
         private final double csc2;
 
         /** Indicator for negated Jacobi amplitude. */
         private boolean negate;
 
         /** Simple constructor.
          * @param phi amplitude (i.e. upper bound of the integral)
          * @param m parameter (m=k² where k is the elliptic modulus)
          * @param integral provider for complete integral
          */
         DoubleArgumentReduction(final double phi, final double m, final DoubleFunction<Double> integral) {
             final double sin = FastMath.sin(phi);
             final int    p   = (int) FastMath.rint(phi / FastMath.PI);
             complete         = p == 0 ? 0 : integral.apply(m) * 2 * p;
             negate           = sin < 0 ^ (p & 0x1) == 1;
             csc2             = 1.0 / (sin * sin);
         }
 
     }
 
     /** Argument reduction for an incomplete integral.
      * @param <T> type fo the field elements
      */
     private static class FieldArgumentReduction<T extends CalculusFieldElement<T>> {
 
         /** Complete part. */
         private final T complete;
 
         /** Squared cosecant of the Jacobi amplitude. */
         private final T csc2;
 
         /** Indicator for negated Jacobi amplitude. */
         private boolean negate;
 
         /** Simple constructor.
          * @param phi amplitude (i.e. upper bound of the integral)
          * @param m parameter (m=k² where k is the elliptic modulus)
          * @param integral provider for complete integral
          */
         FieldArgumentReduction(final T phi, final T m, final Function<T, T> integral) {
             final T   sin = FastMath.sin(phi);
             final int p   = (int) FastMath.rint(phi.getReal() / FastMath.PI);
             complete      = p == 0 ? phi.getField().getZero() : integral.apply(m).multiply(2 * p);
             negate        = sin.getReal() < 0 ^ (p & 0x1) == 1;
             csc2          = sin.multiply(sin).reciprocal();
         }
 
     }
 
     /** Integrand for elliptic integrals of the first kind.
      * @param <T> type of the field elements
      */
     private static class First<T extends CalculusFieldElement<T>> implements CalculusFieldUnivariateFunction<T> {
 
         /** Parameter. */
         private final T m;
 
         /** Simple constructor.
          * @param m parameter (m=k² where k is the elliptic modulus)
          */
         First(final T m) {
             this.m = m;
         }
 
         /** {@inheritDoc} */
         @Override
         public T value(final T theta) {
             final T sin  = theta.sin();
             final T sin2 = sin.multiply(sin);
             return sin2.multiply(m).negate().add(1).sqrt().reciprocal();
         }
 
     }
 
     /** Integrand for elliptic integrals of the second kind.
      * @param <T> type of the field elements
      */
     private static class Second<T extends CalculusFieldElement<T>> implements CalculusFieldUnivariateFunction<T> {
 
         /** Parameter. */
         private final T m;
 
         /** Simple constructor.
          * @param m parameter (m=k² where k is the elliptic modulus)
          */
         Second(final T m) {
             this.m = m;
         }
 
         /** {@inheritDoc} */
         @Override
         public T value(final T theta) {
             final T sin = theta.sin();
             final T sin2 = sin.multiply(sin);
             return sin2.multiply(m).negate().add(1).sqrt();
         }
 
     }
 
     /** Integrand for elliptic integrals of the third kind.
      * @param <T> type of the field elements
      */
     private static class Third<T extends CalculusFieldElement<T>> implements CalculusFieldUnivariateFunction<T> {
 
         /** Elliptic characteristic. */
         private final T n;
 
         /** Parameter. */
         private final T m;
 
         /** Simple constructor.
          * @param n elliptic characteristic
          * @param m parameter (m=k² where k is the elliptic modulus)
          */
         Third(final T n, final T m) {
             this.n = n;
             this.m = m;
         }
 
         /** {@inheritDoc} */
         @Override
         public T value(final T theta) {
             final T sin  = theta.sin();
             final T sin2 = sin.multiply(sin);
             final T d1   = sin2.multiply(m).negate().add(1).sqrt();
             final T da   = sin2.multiply(n).negate().add(1);
             return d1.multiply(da).reciprocal();
         }
 
     }
 }