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