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.complex.Complex;
21 import org.hipparchus.complex.FieldComplex;
22 import org.hipparchus.util.FastMath;
23
24 /** Duplication algorithm for Carlson R<sub>C</sub> elliptic integral.
25 * @param <T> type of the field elements (really {@link Complex} or {@link FieldComplex})
26 * @since 2.0
27 */
28 class RcFieldDuplication<T extends CalculusFieldElement<T>> extends FieldDuplication<T> {
29
30 /** Simple constructor.
31 * @param x first symmetric variable of the integral
32 * @param y second symmetric variable of the integral
33 */
34 RcFieldDuplication(final T x, final T y) {
35 super(x, y);
36 }
37
38 /** {@inheritDoc} */
39 @Override
40 protected void initialMeanPoint(final T[] va) {
41 va[2] = va[0].add(va[1].multiply(2)).divide(3.0);
42 }
43
44 /** {@inheritDoc} */
45 @Override
46 protected T convergenceCriterion(final T r, final T max) {
47 return max.divide(FastMath.sqrt(FastMath.sqrt(FastMath.sqrt(r.multiply(3.0)))));
48 }
49
50 /** {@inheritDoc} */
51 @Override
52 protected void update(final int m, final T[] vaM, final T[] sqrtM, final double fourM) {
53 final T lambdaA = sqrtM[0].multiply(sqrtM[1]).multiply(2);
54 final T lambdaB = vaM[1];
55 vaM[0] = vaM[0].linearCombination(0.25, vaM[0], 0.25, lambdaA, 0.25, lambdaB); // xₘ
56 vaM[1] = vaM[1].linearCombination(0.25, vaM[1], 0.25, lambdaA, 0.25, lambdaB); // yₘ
57 vaM[2] = vaM[2].linearCombination(0.25, vaM[2], 0.25, lambdaA, 0.25, lambdaB); // aₘ
58 }
59
60 /** {@inheritDoc} */
61 @Override
62 protected T evaluate(final T[] va0, final T aM, final double fourM) {
63
64 // compute the single polynomial independent variable
65 final T s = va0[1].subtract(va0[2]).divide(aM.multiply(fourM));
66
67 // evaluate integral using equation 2.13 in Carlson[1995]
68 final T poly = s.multiply(RcRealDuplication.S7).
69 add(RcRealDuplication.S6).multiply(s).
70 add(RcRealDuplication.S5).multiply(s).
71 add(RcRealDuplication.S4).multiply(s).
72 add(RcRealDuplication.S3).multiply(s).
73 add(RcRealDuplication.S2).multiply(s).
74 multiply(s).
75 add(RcRealDuplication.S0).
76 divide(RcRealDuplication.DENOMINATOR);
77 return poly.divide(FastMath.sqrt(aM));
78
79 }
80
81 }