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>C</sub> elliptic integral.
23 * @since 2.0
24 */
25 class RcRealDuplication extends RealDuplication {
26
27 /** Constant term in R<sub>C</sub> polynomial. */
28 static final double S0 = 80080;
29
30 /** Coefficient of s² in R<sub>C</sub> polynomial. */
31 static final double S2 = 24024;
32
33 /** Coefficient of s³ in R<sub>C</sub> polynomial. */
34 static final double S3 = 11440;
35
36 /** Coefficient of s⁴ in R<sub>C</sub> polynomial. */
37 static final double S4 = 30030;
38
39 /** Coefficient of s⁵ in R<sub>C</sub> polynomial. */
40 static final double S5 = 32760;
41
42 /** Coefficient of s⁶ in R<sub>C</sub> polynomial. */
43 static final double S6 = 61215;
44
45 /** Coefficient of s⁷ in R<sub>C</sub> polynomial. */
46 static final double S7 = 90090;
47
48 /** Denominator in R<sub>C</sub> polynomial. */
49 static final double DENOMINATOR = 80080;
50
51 /** Simple constructor.
52 * @param x first symmetric variable of the integral
53 * @param y second symmetric variable of the integral
54 */
55 RcRealDuplication(final double x, final double y) {
56 super(x, y);
57 }
58
59 /** {@inheritDoc} */
60 @Override
61 protected void initialMeanPoint(final double[] va) {
62 va[2] = (va[0] + va[1] * 2) / 3.0;
63 }
64
65 /** {@inheritDoc} */
66 @Override
67 protected double convergenceCriterion(final double r, final double max) {
68 return max / FastMath.sqrt(FastMath.sqrt(FastMath.sqrt(r * 3.0)));
69 }
70
71 /** {@inheritDoc} */
72 @Override
73 protected void update(final int m, final double[] vaM, final double[] sqrtM, final double fourM) {
74 final double lambdaA = sqrtM[0] * sqrtM[1] * 2;
75 final double lambdaB = vaM[1];
76 vaM[0] = MathArrays.linearCombination(0.25, vaM[0], 0.25, lambdaA, 0.25, lambdaB); // xₘ
77 vaM[1] = MathArrays.linearCombination(0.25, vaM[1], 0.25, lambdaA, 0.25, lambdaB); // yₘ
78 vaM[2] = MathArrays.linearCombination(0.25, vaM[2], 0.25, lambdaA, 0.25, lambdaB); // aₘ
79 }
80
81 /** {@inheritDoc} */
82 @Override
83 protected double evaluate(final double[] va0, final double aM, final double fourM) {
84
85 // compute the single polynomial independent variable
86 final double s = (va0[1] - va0[2]) / (aM * fourM);
87
88 // evaluate integral using equation 2.13 in Carlson[1995]
89 final double poly = ((((((S7 * s + S6) * s + S5) * s + S4) * s + S3) * s + S2) * s * s + S0) / DENOMINATOR;
90 return poly / FastMath.sqrt(aM);
91
92 }
93
94 }