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.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 }