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