View Javadoc
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;
18  
19  import org.hipparchus.util.FastMath;
20  
21  /** Duplication algorithm for Carlson symmetric forms.
22   * <p>
23   * The algorithms are described in B. C. Carlson 1995 paper
24   * "Numerical computation of real or complex elliptic integrals", with
25   * improvements described in the appendix of B. C. Carlson and James FitzSimons
26   * 2000 paper "Reduction theorems for elliptic integrands with the square root
27   * of two quadratic factors". They are also described in
28   * <a href="https://dlmf.nist.gov/19.36#i">section 19.36(i)</a>
29   * of Digital Library of Mathematical Functions.
30   * </p>
31   * @since 2.0
32   */
33  abstract class RealDuplication {
34  
35      /** Max number of iterations. */
36      private static final int M_MAX = 16;
37  
38      /** Symmetric variables of the integral, plus mean point. */
39      private final double[] initialVA;
40  
41      /** Convergence criterion. */
42      private final double q;
43  
44      /** Constructor.
45       * @param v symmetric variables of the integral
46       */
47      RealDuplication(final double... v) {
48  
49          final int n = v.length;
50          initialVA = new double[n + 1];
51          System.arraycopy(v, 0, initialVA, 0, n);
52          initialMeanPoint(initialVA);
53  
54          double max = 0;
55          final double a0 = initialVA[n];
56          for (final double vi : v) {
57              max = FastMath.max(max, FastMath.abs(a0 - vi));
58          }
59          this.q = convergenceCriterion(FastMath.ulp(1.0), max);
60  
61      }
62  
63      /** Get the i<sup>th</sup> symmetric variable.
64       * @param i index of the variable
65       * @return i<sup>th</sup> symmetric variable
66       */
67      protected double getVi(final int i) {
68          return initialVA[i];
69      }
70  
71      /** Compute initial mean point.
72       * <p>
73       * The initial mean point is put as the last array element
74       * </>
75       * @param va symmetric variables of the integral (plus placeholder for initial mean point)
76       */
77      protected abstract void initialMeanPoint(double[] va);
78  
79      /** Compute convergence criterion.
80       * @param r relative tolerance
81       * @param max max(|a0-v[i]|)
82       * @return convergence criterion
83       */
84      protected abstract double convergenceCriterion(double r, double max);
85  
86      /** Update reduced variables in place.
87       * <ul>
88       *  <li>vₘ₊₁|i] ← (vₘ[i] + λₘ) / 4</li>
89       *  <li>aₘ₊₁ ← (aₘ + λₘ) / 4</li>
90       * </ul>
91       * @param m iteration index
92       * @param vaM reduced variables and mean point (updated in place)
93       * @param sqrtM square roots of reduced variables
94       * @param fourM 4<sup>m</sup>
95       */
96      protected abstract void update(int m, double[] vaM, double[] sqrtM, double fourM);
97  
98      /** Evaluate integral.
99       * @param va0 initial symmetric variables and mean point of the integral
100      * @param aM reduced mean point
101      * @param fourM 4<sup>m</sup>
102      * @return integral value
103      */
104     protected abstract double evaluate(double[] va0, double aM, double fourM);
105 
106     /** Compute Carlson elliptic integral.
107      * @return Carlson elliptic integral
108      */
109     public double integral() {
110 
111         // duplication iterations
112         final int       n    = initialVA.length - 1;
113         final double[] vaM   = initialVA.clone();
114         final double[] sqrtM = new double[n];
115         double         fourM = 1.0;
116         for (int m = 0; m < M_MAX; ++m) {
117 
118             if (m > 0 && q < fourM * FastMath.abs(vaM[n])) {
119                 // convergence reached
120                 break;
121             }
122 
123             // apply duplication once more
124             // (we know that {Field}Complex.sqrt() returns the root with nonnegative real part)
125             for (int i = 0; i < n; ++i) {
126                 sqrtM[i] = FastMath.sqrt(vaM[i]);
127             }
128             update(m, vaM, sqrtM, fourM);
129 
130             fourM *= 4;
131 
132         }
133 
134         return evaluate(initialVA, vaM[n], fourM);
135 
136     }
137 
138 }