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 import org.hipparchus.util.MathArrays;
21
22 /** Duplication algorithm for Carlson R<sub>F</sub> elliptic integral.
23 * @since 2.0
24 */
25 class RfRealDuplication extends RealDuplication {
26
27 /** Max number of iterations in the AGM scale. */
28 static final int AGM_MAX = 32;
29
30 /** Constant term in R<sub>F</sub> polynomial. */
31 static final double CONSTANT = 240240;
32
33 /** Coefficient of E₂ in R<sub>F</sub> polynomial. */
34 static final double E2 = -24024;
35
36 /** Coefficient of E₃ in R<sub>F</sub> polynomial. */
37 static final double E3 = 17160;
38
39 /** Coefficient of E₂² in R<sub>F</sub> polynomial. */
40 static final double E2_E2 = 10010;
41
42 /** Coefficient of E₂E₃ in R<sub>F</sub> polynomial. */
43 static final double E2_E3 = -16380;
44
45 /** Coefficient of E₃² in R<sub>F</sub> polynomial. */
46 static final double E3_E3 = 6930;
47
48 /** Coefficient of E₂³ in R<sub>F</sub> polynomial. */
49 static final double E2_E2_E2 = -5775;
50
51 /** Denominator in R<sub>F</sub> polynomial. */
52 static final double DENOMINATOR = 240240;
53
54 /** Simple constructor.
55 * @param x first symmetric variable of the integral
56 * @param y second symmetric variable of the integral
57 * @param z third symmetric variable of the integral
58 */
59 RfRealDuplication(final double x, final double y, final double z) {
60 super(x, y, z);
61 }
62
63 /** {@inheritDoc} */
64 @Override
65 protected void initialMeanPoint(final double[] va) {
66 va[3] = (va[0] + va[1] + va[2]) / 3.0;
67 }
68
69 /** {@inheritDoc} */
70 @Override
71 protected double convergenceCriterion(final double r, final double max) {
72 return max / FastMath.sqrt(FastMath.sqrt(FastMath.sqrt(r * 3.0)));
73 }
74
75 /** {@inheritDoc} */
76 @Override
77 protected void update(final int m, final double[] vaM, final double[] sqrtM, final double fourM) {
78
79 // equation 2.3 in Carlson[1995]
80 final double lambdaA = sqrtM[0] * sqrtM[1];
81 final double lambdaB = sqrtM[0] * sqrtM[2];
82 final double lambdaC = sqrtM[1] * sqrtM[2];
83
84 // equations 2.3 and 2.4 in Carlson[1995]
85 vaM[0] = MathArrays.linearCombination(0.25, vaM[0], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // xₘ
86 vaM[1] = MathArrays.linearCombination(0.25, vaM[1], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // yₘ
87 vaM[2] = MathArrays.linearCombination(0.25, vaM[2], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // zₘ
88 vaM[3] = MathArrays.linearCombination(0.25, vaM[3], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // aₘ
89
90 }
91
92 /** {@inheritDoc} */
93 @Override
94 protected double evaluate(final double[] va0, final double aM, final double fourM) {
95
96 // compute symmetric differences
97 final double inv = 1.0 / (aM * fourM);
98 final double bigX = (va0[3] - va0[0]) * inv;
99 final double bigY = (va0[3] - va0[1]) * inv;
100 final double bigZ = -(bigX + bigY);
101
102 // compute elementary symmetric functions (we already know e1 = 0 by construction)
103 final double e2 = bigX * bigY - bigZ * bigZ;
104 final double e3 = bigX * bigY * bigZ;
105
106 final double e2e2 = e2 * e2;
107 final double e2e3 = e2 * e3;
108 final double e3e3 = e3 * e3;
109 final double e2e2e2 = e2e2 * e2;
110
111 // evaluate integral using equation 19.36.1 in DLMF
112 // (which add more terms than equation 2.7 in Carlson[1995])
113 final double poly = (e2e2e2 * E2_E2_E2 +
114 e3e3 * E3_E3 +
115 e2e3 * E2_E3 +
116 e2e2 * E2_E2 +
117 e3 * E3 +
118 e2 * E2 +
119 CONSTANT) /
120 DENOMINATOR;
121 return poly / FastMath.sqrt(aM);
122
123 }
124
125 /** {@inheritDoc} */
126 @Override
127 public double integral() {
128 final double x = getVi(0);
129 final double y = getVi(1);
130 final double z = getVi(2);
131 if (x == 0) {
132 return completeIntegral(y, z);
133 } else if (y == 0) {
134 return completeIntegral(x, z);
135 } else if (z == 0) {
136 return completeIntegral(x, y);
137 } else {
138 return super.integral();
139 }
140 }
141
142 /** Compute Carlson complete elliptic integral R<sub>F</sub>(u, v, 0).
143 * @param x first symmetric variable of the integral
144 * @param y second symmetric variable of the integral
145 * @return Carlson complete elliptic integral R<sub>F</sub>(u, v, 0)
146 */
147 private double completeIntegral(final double x, final double y) {
148
149 double xM = FastMath.sqrt(x);
150 double yM = FastMath.sqrt(y);
151
152 // iterate down
153 for (int i = 1; i < AGM_MAX; ++i) {
154
155 final double xM1 = xM;
156 final double yM1 = yM;
157
158 // arithmetic mean
159 xM = (xM1 + yM1) * 0.5;
160
161 // geometric mean
162 yM = FastMath.sqrt(xM1 * yM1);
163
164 // convergence (by the inequality of arithmetic and geometric means, this is non-negative)
165 if (FastMath.abs(xM - yM) <= 4 * FastMath.ulp(xM)) {
166 // convergence has been reached
167 break;
168 }
169
170 }
171
172 return FastMath.PI / (xM + yM);
173
174 }
175
176 }