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 }