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.random.RandomGenerator;
20  import org.hipparchus.random.Well19937a;
21  import org.hipparchus.random.Well19937c;
22  import org.hipparchus.util.FastMath;
23  import org.junit.jupiter.api.Test;
24  
25  import static org.junit.jupiter.api.Assertions.assertEquals;
26  import static org.junit.jupiter.api.Assertions.assertTrue;
27  
28  class CarlsonEllipticIntegralRealTest {
29  
30      @Test
31      void testNoConvergenceRf() {
32          assertTrue(Double.isNaN(CarlsonEllipticIntegral.rF(1, 2, Double.NaN)));
33      }
34  
35      @Test
36      void testDlmfRf() {
37          double rf = CarlsonEllipticIntegral.rF(1, 2, 4);
38          assertEquals(0.6850858166, rf, 1.0e-10);
39      }
40  
41      @Test
42      void testCarlson1995rF() {
43  
44          double rf1 = CarlsonEllipticIntegral.rF(1, 2, 0);
45          assertEquals( 1.3110287771461, rf1, 1.0e-13);
46  
47          double rf2 = CarlsonEllipticIntegral.rF(0.5, 1, 0);
48          assertEquals( 1.8540746773014, rf2, 1.0e-13);
49  
50          double rf4 = CarlsonEllipticIntegral.rF(2, 3, 4);
51          assertEquals( 0.58408284167715, rf4, 1.0e-13);
52  
53      }
54  
55      @Test
56      void testCarlson1995ConsistencyRf() {
57          RandomGenerator random = new Well19937c(0x57f2689b3f4028b4l);
58          for (int i = 0; i < 10000; ++i) {
59              double x      = random.nextDouble() * 3;
60              double y      = random.nextDouble() * 3;
61              double lambda = random.nextDouble() * 3;
62              double mu     = x * y / lambda;
63              double rfL    = CarlsonEllipticIntegral.rF(x + lambda, y + lambda, lambda);
64              double rfM    = CarlsonEllipticIntegral.rF(x + mu,     y + mu,     mu);
65              double rf0    = CarlsonEllipticIntegral.rF(x,             y,             0);
66              assertEquals(0.0, FastMath.abs(rfL + rfM - rf0), 2.0e-14);
67          }
68      }
69  
70      @Test
71      void testNoConvergenceRc() {
72          assertTrue(Double.isNaN(CarlsonEllipticIntegral.rC(1, Double.NaN)));
73      }
74  
75      @Test
76      void testCarlson1995rC() {
77  
78          double rc1 = CarlsonEllipticIntegral.rC(0, 0.25);
79          assertEquals(FastMath.PI, rc1, 1.0e-15);
80  
81          double rc2 = CarlsonEllipticIntegral.rC(2.25, 2);
82          assertEquals(FastMath.log(2), rc2, 1.0e-15);
83  
84          double rc5 = CarlsonEllipticIntegral.rC(0.25, -2);
85          assertEquals(FastMath.log(2) / 3.0, rc5, 1.0e-15);
86  
87      }
88  
89      @Test
90      void testCarlson1995ConsistencyRc() {
91          RandomGenerator random = new Well19937c(0xf1170b6fc1a199cal);
92          for (int i = 0; i < 10000; ++i) {
93              double x      = random.nextDouble() * 3;
94              double lambda = random.nextDouble() * 3;
95              double mu     = x * x / lambda;
96              double rcL    = CarlsonEllipticIntegral.rC(lambda,          x + lambda);
97              double rcM    = CarlsonEllipticIntegral.rC(mu,              x + mu);
98              double rc0    = CarlsonEllipticIntegral.rC(0, x);
99              assertEquals(0.0, FastMath.abs(rcL + rcM - rc0), 3.0e-14);
100         }
101     }
102 
103     @Test
104     void testRfRc() {
105         RandomGenerator random = new Well19937a(0x7e8041334a8c20edl);
106         for (int i = 0; i < 10000; ++i) {
107             final double x = 3 * random.nextDouble();
108             final double y = 3 * random.nextDouble();
109             final double rf = CarlsonEllipticIntegral.rF(x, y, y);
110             final double rc = CarlsonEllipticIntegral.rC(x, y);
111             assertEquals(0.0, FastMath.abs(rf - rc), 4.0e-15);
112         }
113     }
114 
115     @Test
116     void testNoConvergenceRj() {
117         assertTrue(Double.isNaN(CarlsonEllipticIntegral.rJ(1, 1, 1, Double.NaN)));
118     }
119 
120     @Test
121     void testCarlson1995rJ() {
122 
123         double rj01 = CarlsonEllipticIntegral.rJ(0, 1, 2, 3);
124         assertEquals(0.77688623778582, rj01, 1.0e-13);
125 
126         double rj02 = CarlsonEllipticIntegral.rJ(2, 3, 4, 5);
127         assertEquals( 0.14297579667157, rj02, 1.0e-13);
128 
129     }
130 
131     @Test
132     void testCarlson1995ConsistencyRj() {
133         RandomGenerator random = new Well19937c(0x4af7bb722712e64el);
134         for (int i = 0; i < 10000; ++i) {
135             double x      = random.nextDouble() * 3;
136             double y      = random.nextDouble() * 3;
137             double p      = random.nextDouble() * 3;
138             double lambda = random.nextDouble() * 3;
139             double mu     = x * y / lambda;
140             double a      = p * p * (lambda + mu + x + y);
141             double b      = p * (p + lambda) * (p + mu);
142             double rjL    = CarlsonEllipticIntegral.rJ(x + lambda, y + lambda, lambda,  p + lambda);
143             double rjM    = CarlsonEllipticIntegral.rJ(x + mu,     y + mu,     mu,      p + mu);
144             double rj0    = CarlsonEllipticIntegral.rJ(x,          y,          0,       p);
145             double rc     = CarlsonEllipticIntegral.rC(a, b);
146             assertEquals(0.0, FastMath.abs(rjL + rjM - (rj0 - rc * 3)), 3.0e-13);
147         }
148     }
149 
150     @Test
151     void testNoConvergenceRd() {
152         assertTrue(Double.isNaN(CarlsonEllipticIntegral.rD(1, 1, Double.NaN)));
153     }
154 
155     @Test
156     void testCarlson1995rD() {
157 
158         double rd1 = CarlsonEllipticIntegral.rD(0, 2, 1);
159         assertEquals(1.7972103521034, rd1, 1.0e-13);
160 
161         double rd2 = CarlsonEllipticIntegral.rD(2, 3, 4);
162         assertEquals( 0.16510527294261, rd2, 1.0e-13);
163 
164     }
165 
166     @Test
167     void testCarlson1995ConsistencyRd() {
168         RandomGenerator random = new Well19937c(0x17dea97eeb78206al);
169         for (int i = 0; i < 10000; ++i) {
170             double x      = random.nextDouble() * 3;
171             double y      = random.nextDouble() * 3;
172             double lambda = random.nextDouble() * 3;
173             double mu     = x * y / lambda;
174             double rdL    = CarlsonEllipticIntegral.rD(lambda,          x + lambda, y + lambda);
175             double rdM    = CarlsonEllipticIntegral.rD(mu,              x + mu,     y + mu);
176             double rd0    = CarlsonEllipticIntegral.rD(0,               x,          y);
177             double frac   = 3 / (y * FastMath.sqrt(x + y + lambda + mu));
178             assertEquals(0.0, FastMath.abs(rdL + rdM - rd0 + frac), 9.0e-12);
179         }
180     }
181 
182     @Test
183     void testRdNonSymmetry1() {
184         RandomGenerator random = new Well19937c(0x66db170b5ee1afc2l);
185         for (int i = 0; i < 10000; ++i) {
186             double x = random.nextDouble();
187             double y = random.nextDouble();
188             double z = random.nextDouble();
189             if (x == 0 || y == 0) {
190                 continue;
191             }
192             // this is DLMF equation 19.21.7
193             double lhs = (x - y) * CarlsonEllipticIntegral.rD(y, z, x) + (z - y) * CarlsonEllipticIntegral.rD(x, y, z);
194             double rhs = (CarlsonEllipticIntegral.rF(x, y, z) - FastMath.sqrt(y / (x * z))) * 3;
195             assertEquals(0.0, FastMath.abs(lhs - rhs), 1.0e-10);
196         }
197     }
198 
199     @Test
200     void testRdNonSymmetry2() {
201         RandomGenerator random = new Well19937c(0x1a8994acc807438dl);
202         for (int i = 0; i < 10000; ++i) {
203             double x = random.nextDouble();
204             double y = random.nextDouble();
205             double z = random.nextDouble();
206             if (x == 0 || y == 0 || z == 0) {
207                 continue;
208             }
209             // this is DLMF equation 19.21.8
210             double lhs = CarlsonEllipticIntegral.rD(y, z, x) + CarlsonEllipticIntegral.rD(z, x, y) + CarlsonEllipticIntegral.rD(x, y, z);
211             double rhs = 3 / FastMath.sqrt(x * y * z);
212             assertEquals(0.0, FastMath.abs(lhs - rhs), 2.0e-11);
213         }
214     }
215 
216     @Test
217     void testCarlson1995rG() {
218 
219         double rg1 = CarlsonEllipticIntegral.rG(0, 16, 16);
220         assertEquals(FastMath.PI, rg1, 1.0e-13);
221 
222         double rg2 = CarlsonEllipticIntegral.rG(2, 3, 4);
223         assertEquals(1.7255030280692, rg2, 1.0e-13);
224 
225         double rg6 = CarlsonEllipticIntegral.rG(0, 0.0796, 4);
226         assertEquals( 1.0284758090288, rg6, 1.0e-13);
227 
228     }
229 
230     @Test
231     void testAlternateRG() {
232         RandomGenerator random = new Well19937c(0xa2946e4a55d133a6l);
233         for (int i = 0; i < 10000; ++i) {
234             double x = random.nextDouble() * 3;
235             double y = random.nextDouble() * 3;
236             double z = random.nextDouble() * 3;
237             assertEquals(0.0, FastMath.abs(CarlsonEllipticIntegral.rG(x, y, z) - rgAlternateImplementation(x, y, z)), 2.0e-15);
238         }
239     }
240 
241     private double rgAlternateImplementation(final double x, final double y, final double z) {
242         // this implementation uses DLFM equation 19.21.11
243         return (d(x, y, z) + d(y, z, x) + d(z, x, y)) / 6;
244     }
245 
246     private double d(final double u, final double v, final double w) {
247         return u == 0 ? u : u * (v + w) * (new RdRealDuplication(v, w, u).integral());
248     }
249 
250 }