View Javadoc
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>D</sub> elliptic integral.
25   * @param <T> type of the field elements (really {@link Complex} or {@link FieldComplex})
26   * @since 2.0
27   */
28  class RdFieldDuplication<T extends CalculusFieldElement<T>> extends FieldDuplication<T> {
29  
30      /** Partial sum. */
31      private T sum;
32  
33      /** Simple constructor.
34       * @param x first symmetric variable of the integral
35       * @param y second symmetric variable of the integral
36       * @param z third symmetric variable of the integral
37       */
38      RdFieldDuplication(final T x, final T y, final T z) {
39          super(x, y, z);
40          sum = x.getField().getZero();
41      }
42  
43      /** {@inheritDoc} */
44      @Override
45      protected void initialMeanPoint(final T[] va) {
46          va[3] = va[0].add(va[1]).add(va[2].multiply(3)).divide(5.0);
47      }
48  
49      /** {@inheritDoc} */
50      @Override
51      protected T convergenceCriterion(final T r, final T max) {
52          return max.divide(FastMath.sqrt(FastMath.sqrt(FastMath.sqrt(r.multiply(0.25)))));
53      }
54  
55      /** {@inheritDoc} */
56      @Override
57      protected void update(final int m, final T[] vaM, final T[] sqrtM, final  double fourM) {
58  
59          // equation 2.29 in Carlson[1995]
60          final T lambdaA = sqrtM[0].multiply(sqrtM[1]);
61          final T lambdaB = sqrtM[0].multiply(sqrtM[2]);
62          final T lambdaC = sqrtM[1].multiply(sqrtM[2]);
63  
64          // running sum in equation 2.34 in Carlson[1995]
65          final T lambda = lambdaA.add(lambdaB).add(lambdaC);
66          sum = sum.add(vaM[2].add(lambda).multiply(sqrtM[2]).multiply(fourM).reciprocal());
67  
68          // equations 2.29 and 2.30 in Carlson[1995]
69          vaM[0] = vaM[0].linearCombination(0.25, vaM[0], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // xₘ
70          vaM[1] = vaM[1].linearCombination(0.25, vaM[1], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // yₘ
71          vaM[2] = vaM[2].linearCombination(0.25, vaM[2], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // zₘ
72          vaM[3] = vaM[3].linearCombination(0.25, vaM[3], 0.25, lambdaA, 0.25, lambdaB, 0.25, lambdaC); // aₘ
73  
74      }
75  
76      /** {@inheritDoc} */
77      @Override
78      protected T evaluate(final T[] va0, final T aM, final  double fourM) {
79  
80          // compute symmetric differences
81          final T inv   = aM.multiply(fourM).reciprocal();
82          final T bigX  = va0[3].subtract(va0[0]).multiply(inv);
83          final T bigY  = va0[3].subtract(va0[1]).multiply(inv);
84          final T bigZ  = bigX.add(bigY).divide(-3);
85          final T bigXY = bigX.multiply(bigY);
86          final T bigZ2 = bigZ.multiply(bigZ);
87  
88          // compute elementary symmetric functions (we already know e1 = 0 by construction)
89          final T e2  = bigXY.subtract(bigZ2.multiply(6));
90          final T e3  = bigXY.multiply(3).subtract(bigZ2.multiply(8)).multiply(bigZ);
91          final T e4  = bigXY.subtract(bigZ2).multiply(3).multiply(bigZ2);
92          final T e5  = bigXY.multiply(bigZ2).multiply(bigZ);
93  
94          final T e2e2   = e2.multiply(e2);
95          final T e2e3   = e2.multiply(e3);
96          final T e2e4   = e2.multiply(e4);
97          final T e2e5   = e2.multiply(e5);
98          final T e3e3   = e3.multiply(e3);
99          final T e3e4   = e3.multiply(e4);
100         final T e2e2e2 = e2e2.multiply(e2);
101         final T e2e2e3 = e2e2.multiply(e3);
102 
103         // evaluate integral using equation 19.36.1 in DLMF
104         // (which add more terms than equation 2.7 in Carlson[1995])
105         final T poly = e3e4.add(e2e5).multiply(RdRealDuplication.E3_E4_P_E2_E5).
106                        add(e2e2e3.multiply(RdRealDuplication.E2_E2_E3)).
107                        add(e2e4.multiply(RdRealDuplication.E2_E4)).
108                        add(e3e3.multiply(RdRealDuplication.E3_E3)).
109                        add(e2e2e2.multiply(RdRealDuplication.E2_E2_E2)).
110                        add(e5.multiply(RdRealDuplication.E5)).
111                        add(e2e3.multiply(RdRealDuplication.E2_E3)).
112                        add(e4.multiply(RdRealDuplication.E4)).
113                        add(e2e2.multiply(RdRealDuplication.E2_E2)).
114                        add(e3.multiply(RdRealDuplication.E3)).
115                        add(e2.multiply(RdRealDuplication.E2)).
116                        add(RdRealDuplication.CONSTANT).
117                        divide(RdRealDuplication.DENOMINATOR);
118         final T polyTerm = poly.divide(aM.multiply(FastMath.sqrt(aM)).multiply(fourM));
119 
120         return polyTerm.add(sum.multiply(3));
121 
122     }
123 
124 }