1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.hipparchus.special.elliptic.carlson;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.random.RandomGenerator;
21 import org.hipparchus.random.Well19937a;
22 import org.hipparchus.random.Well19937c;
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 public abstract class CarlsonEllipticIntegralAbstractComplexTest<T extends CalculusFieldElement<T>> {
30
31 protected abstract T buildComplex(double realPart);
32 protected abstract T buildComplex(double realPart, double imaginaryPart);
33 protected abstract T rF(T x, T y, T z);
34 protected abstract T rC(T x, T y);
35 protected abstract T rJ(T x, T y, T z, T p);
36 protected abstract T rD(T x, T y, T z);
37 protected abstract T rG(T x, T y, T z);
38
39 private void check(double expectedReal, double expectedImaginary, T result, double tol) {
40 assertEquals(0, buildComplex(expectedReal, expectedImaginary).subtract(result).norm(), tol);
41 }
42
43 @Test
44 public void testNoConvergenceRf() {
45 assertTrue(rF(buildComplex(1), buildComplex(2), buildComplex(Double.NaN)).isNaN());
46 }
47
48 @Test
49 public void testDlmfRf() {
50 T rf = rF(buildComplex(1), buildComplex(2), buildComplex(4));
51 check(0.6850858166, 0.0, rf, 1.0e-10);
52 }
53
54 @Test
55 public void testCarlson1995rF() {
56
57 T rf1 = rF(buildComplex(1), buildComplex(2), buildComplex(0));
58 check( 1.3110287771461, 0.0, rf1, 1.0e-13);
59
60 T rf2 = rF(buildComplex(0.5), buildComplex(1), buildComplex(0));
61 check( 1.8540746773014, 0.0, rf2, 1.0e-13);
62
63 T rf3 = rF(buildComplex(-1, 1), buildComplex(0, 1), buildComplex(0));
64 check( 0.79612586584234, -1.2138566698365, rf3, 1.0e-13);
65
66 T rf4 = rF(buildComplex(2), buildComplex(3), buildComplex(4));
67 check( 0.58408284167715, 0.0, rf4, 1.0e-13);
68
69 T rf5 = rF(buildComplex(0, 1), buildComplex(0, -1), buildComplex(2));
70 check( 1.0441445654064, 0.0, rf5, 1.0e-13);
71
72 T rf6 = rF(buildComplex(-1, 1), buildComplex(0, 1), buildComplex(1, -1));
73 check( 0.93912050218619, -0.53296252018635, rf6, 1.0e-13);
74
75 }
76
77 @Test
78 public void testRfAlongImaginaryAxis() {
79 final T x = buildComplex(0, 1);
80 final T yN = buildComplex(0, -1 - 1.0e-13);
81 final T y0 = buildComplex(0, -1);
82 final T yP = buildComplex(0, -1 + 1.0e-13);
83 final T z = buildComplex(0);
84 check(1.8540746773013255, +2.12e-14, rF(x, yN, z), 1.0e-15);
85 check(1.8540746773013719, 0.0, rF(x, y0, z), 1.0e-15);
86 check(1.8540746773014183, -2.11e-14, rF(x, yP, z), 1.0e-15);
87 }
88
89 @Test
90 public void testCarlson1995ConsistencyRf() {
91 RandomGenerator random = new Well19937c(0x57f2689b3f4028b4l);
92 for (int i = 0; i < 10000; ++i) {
93 T x = buildComplex(random.nextDouble() * 3);
94 T y = buildComplex(random.nextDouble() * 3);
95 T lambda = buildComplex(random.nextDouble() * 6 - 3, random.nextDouble() * 3);
96 T mu = x.multiply(y).divide(lambda);
97 T rfL = rF(x.add(lambda), y.add(lambda), lambda);
98 T rfM = rF(x.add(mu), y.add(mu), mu);
99 T rf0 = rF(x, y, buildComplex(0));
100 assertEquals(0.0, rfL.add(rfM).subtract(rf0).norm(), 2.0e-14);
101 }
102 }
103
104 @Test
105 public void testNoConvergenceRc() {
106 assertTrue(rC(buildComplex(1), buildComplex(Double.NaN)).isNaN());
107 }
108
109 @Test
110 public void testCarlson1995rC() {
111
112 T rc1 = rC(buildComplex(0), buildComplex(0.25));
113 check(FastMath.PI, 0.0, rc1, 1.0e-15);
114
115 T rc2 = rC(buildComplex(2.25), buildComplex(2));
116 check(FastMath.log(2), 0.0, rc2, 1.0e-15);
117
118 T rc3 = rC(buildComplex(0), buildComplex(0, 1));
119 check( 1.1107207345396, -1.1107207345396, rc3, 1.0e-13);
120
121 T rc4 = rC(buildComplex(0, -1), buildComplex(0, 1));
122 check( 1.2260849569072, -0.34471136988768, rc4, 1.0e-13);
123
124 T rc5 = rC(buildComplex(0.25), buildComplex(-2));
125 check(FastMath.log(2) / 3.0, 0.0, rc5, 1.0e-15);
126
127 T rc6 = rC(buildComplex(0, 1), buildComplex(-1));
128 check( 0.77778596920447, 0.19832484993429, rc6, 1.0e-13);
129
130 }
131
132 @Test
133 public void testCarlson1995ConsistencyRc() {
134 RandomGenerator random = new Well19937c(0xf1170b6fc1a199cal);
135 for (int i = 0; i < 10000; ++i) {
136 T x = buildComplex(random.nextDouble() * 3);
137 T lambda = buildComplex(random.nextDouble() * 6 - 3, random.nextDouble() * 3);
138 T mu = x.square().divide(lambda);
139 T rcL = rC(lambda, x.add(lambda));
140 T rcM = rC(mu, x.add(mu));
141 T rc0 = rC(buildComplex(0), x);
142 assertEquals(0.0, rcL.add(rcM).subtract(rc0).norm(), 3.0e-14);
143 }
144 }
145
146 @Test
147 public void testRfRc() {
148 RandomGenerator random = new Well19937a(0x7e8041334a8c20edl);
149 for (int i = 0; i < 10000; ++i) {
150 final T x = buildComplex(6 * random.nextDouble() - 3,
151 6 * random.nextDouble() - 3);
152 final T y = buildComplex(6 * random.nextDouble() - 3,
153 6 * random.nextDouble() - 3);
154 final T rf = rF(x, y, y);
155 final T rc = rC(x, y);
156 assertEquals(0.0, rf.subtract(rc).norm(), 4.0e-15);
157 }
158 }
159
160 @Test
161 public void testNoConvergenceRj() {
162 assertTrue(rJ(buildComplex(1), buildComplex(1), buildComplex(1), buildComplex(Double.NaN)).isNaN());
163 }
164
165 @Test
166 public void testCarlson1995rJ() {
167
168 T rj01 = rJ(buildComplex(0), buildComplex(1), buildComplex(2), buildComplex(3));
169 check(0.77688623778582, 0.0, rj01, 1.0e-13);
170
171 T rj02 = rJ(buildComplex(2), buildComplex(3), buildComplex(4), buildComplex(5));
172 check( 0.14297579667157, 0.0, rj02, 1.0e-13);
173
174 T rj03 = rJ(buildComplex(2), buildComplex(3), buildComplex(4), buildComplex(-1, 1));
175 check( 0.13613945827771, -0.38207561624427, rj03, 1.0e-13);
176
177 T rj04 = rJ(buildComplex(0, 1), buildComplex(0, -1), buildComplex(0), buildComplex(2));
178 check( 1.6490011662711, 0.0, rj04, 1.0e-13);
179
180 T rj05 = rJ(buildComplex(-1, 1), buildComplex(-1, -1), buildComplex(1), buildComplex(2));
181 check( 0.94148358841220, 0.0, rj05, 1.0e-13);
182
183 T rj06 = rJ(buildComplex(0, 1), buildComplex(0, -1), buildComplex(0), buildComplex(1, -1));
184 check( 1.8260115229009, 1.2290661908643, rj06, 1.0e-13);
185
186 T rj07 = rJ(buildComplex(-1, 1), buildComplex(-1, -1), buildComplex(1), buildComplex(-3, 1));
187 check(-0.61127970812028, -1.0684038390007, rj07, 1.0e-13);
188
189 T rj08 = rJ(buildComplex(-1, 1), buildComplex(-2, -1), buildComplex(0, -1), buildComplex(-1, 1));
190 check( 1.8249027393704, -1.2218475784827, rj08, 1.0e-13);
191
192 T rj09 = rJ(buildComplex(2), buildComplex(3), buildComplex(4), buildComplex(-0.5));
193 check( 0.24723819703052, -0.7509842836891, rj09, 1.0e-13);
194
195 T rj10 = rJ(buildComplex(2), buildComplex(3), buildComplex(4), buildComplex(-5));
196 check(-0.12711230042964, -0.2099064885453, rj10, 1.0e-13);
197
198 }
199
200 @Test
201 public void testCarlson1995ConsistencyRj() {
202 RandomGenerator random = new Well19937c(0x4af7bb722712e64el);
203 for (int i = 0; i < 10000; ++i) {
204 T x = buildComplex(random.nextDouble() * 3);
205 T y = buildComplex(random.nextDouble() * 3);
206 T p = buildComplex(random.nextDouble() * 3);
207 T lambda = buildComplex(random.nextDouble() * 6 - 3, random.nextDouble() * 3);
208 T mu = x.multiply(y).divide(lambda);
209 T a = p.multiply(p).multiply(lambda.add(mu).add(x).add(y));
210 T b = p.multiply(p.add(lambda)).multiply(p.add(mu));
211 T rjL = rJ(x.add(lambda), y.add(lambda), lambda, p.add(lambda));
212 T rjM = rJ(x.add(mu), y.add(mu), mu, p.add(mu));
213 T rj0 = rJ(x, y, buildComplex(0), p);
214 T rc = rC(a, b);
215 assertEquals(0.0, rjL.add(rjM).subtract(rj0.subtract(rc.multiply(3))).norm(), 3.0e-13);
216 }
217 }
218
219 @Test
220 public void testNoConvergenceRd() {
221 assertTrue(rD(buildComplex(1), buildComplex(1), buildComplex(Double.NaN)).isNaN());
222 }
223
224 @Test
225 public void testCarlson1995rD() {
226
227 T rd1 = rD(buildComplex(0), buildComplex(2), buildComplex(1));
228 check(1.7972103521034, 0.0, rd1, 1.0e-13);
229
230 T rd2 = rD(buildComplex(2), buildComplex(3), buildComplex(4));
231 check( 0.16510527294261, 0.0, rd2, 1.0e-13);
232
233 T rd3 = rD(buildComplex(0, 1), buildComplex(0, -1), buildComplex(2));
234 check( 0.65933854154220, 0.0, rd3, 1.0e-13);
235
236 T rd4 = rD(buildComplex(0), buildComplex(0, 1), buildComplex(0, -1));
237 check( 1.2708196271910, 2.7811120159521, rd4, 1.0e-13);
238
239 T rd5 = rD(buildComplex(0), buildComplex(-1, 1), buildComplex(0, 1));
240 check(-1.8577235439239, -0.96193450888830, rd5, 1.0e-13);
241
242 T rd6 = rD(buildComplex(-2, -1), buildComplex(0, -1), buildComplex(-1, 1));
243 check( 1.8249027393704, -1.2218475784827, rd6, 1.0e-13);
244
245 }
246
247 @Test
248 public void testCarlson1995ConsistencyRd() {
249 RandomGenerator random = new Well19937c(0x17dea97eeb78206al);
250 for (int i = 0; i < 10000; ++i) {
251 T x = buildComplex(random.nextDouble() * 3);
252 T y = buildComplex(random.nextDouble() * 3);
253 T lambda = buildComplex(random.nextDouble() * 6 - 3, random.nextDouble() * 3);
254 T mu = x.multiply(y).divide(lambda);
255 T rdL = rD(lambda, x.add(lambda), y.add(lambda));
256 T rdM = rD(mu, x.add(mu), y.add(mu));
257 T rd0 = rD(buildComplex(0), x, y);
258 T frac = y.multiply(x.add(y).add(lambda).add(mu).sqrt()).reciprocal().multiply(3);
259 assertEquals(0.0, rdL.add(rdM).subtract(rd0.subtract(frac)).norm(), 9.0e-12);
260 }
261 }
262
263 @Test
264 public void testRdNonSymmetry1() {
265 RandomGenerator random = new Well19937c(0x66db170b5ee1afc2l);
266 int countWrongRoot = 0;
267 for (int i = 0; i < 10000; ++i) {
268 T x = buildComplex(random.nextDouble() * 2 - 1, random.nextDouble() * 2 - 1);
269 T y = buildComplex(random.nextDouble() * 2 - 1, random.nextDouble() * 2 - 1);
270 T z = buildComplex(random.nextDouble() * 2 - 1, random.nextDouble() * 2 - 1);
271 if (x.isZero() || y.isZero()) {
272 continue;
273 }
274
275
276
277
278 T lhs = x.subtract(y).multiply(rD(y, z, x)).add(z.subtract(y).multiply(rD(x, y, z)));
279 T rootGlobal = y.divide(x.multiply(z)).sqrt();
280 T rootSeparated = y.sqrt().divide(x.sqrt().multiply(z.sqrt()));
281 T rhsGlobal = rF(x, y, z).subtract(rootGlobal).multiply(3);
282 T rhsSeparated = rF(x, y, z).subtract(rootSeparated).multiply(3);
283 if (lhs.subtract(rhsGlobal).norm() > 1.0e-3) {
284 ++countWrongRoot;
285
286 assertTrue(lhs.subtract(rhsGlobal).norm() > 0.1);
287 }
288 assertEquals(0.0, lhs.subtract(rhsSeparated).norm(), 1.0e-10);
289 }
290 assertTrue(countWrongRoot > 3300);
291 }
292
293 @Test
294 public void testRdNonSymmetry2() {
295 RandomGenerator random = new Well19937c(0x1a8994acc807438dl);
296 int countWrongRoot = 0;
297 for (int i = 0; i < 10000; ++i) {
298 T x = buildComplex(random.nextDouble() * 2 - 1, random.nextDouble() * 2 - 1);
299 T y = buildComplex(random.nextDouble() * 2 - 1, random.nextDouble() * 2 - 1);
300 T z = buildComplex(random.nextDouble() * 2 - 1, random.nextDouble() * 2 - 1);
301 if (x.isZero() || y.isZero() || z.isZero()) {
302 continue;
303 }
304
305
306
307
308 T lhs = rD(y, z, x).add(rD(z, x, y)).add(rD(x, y, z));
309 T rootGlobal = x.multiply(y.multiply(z)).sqrt();
310 T rootSeparated = x.sqrt().multiply(y.sqrt().multiply(z.sqrt()));
311 T rhsGlobal = rootGlobal.reciprocal().multiply(3);
312 T rhsSeparated = rootSeparated.reciprocal().multiply(3);
313 if (lhs.subtract(rhsGlobal).norm() > 1.0e-3) {
314 ++countWrongRoot;
315
316 assertTrue(lhs.subtract(rhsGlobal).norm() > 3.0);
317 }
318 assertEquals(0.0, lhs.subtract(rhsSeparated).norm(), 2.0e-11);
319 }
320 assertTrue(countWrongRoot > 3300);
321 }
322
323 @Test
324 public void testCarlson1995rG() {
325
326 T rg1 = rG(buildComplex(0), buildComplex(16), buildComplex(16));
327 check(FastMath.PI, 0.0, rg1, 1.0e-13);
328
329 T rg2 = rG(buildComplex(2), buildComplex(3), buildComplex(4));
330 check(1.7255030280692, 0.0, rg2, 1.0e-13);
331
332 T rg3 = rG(buildComplex(0), buildComplex(0, 1), buildComplex(0, -1));
333 check( 0.42360654239699, 0.0, rg3, 1.0e-13);
334
335 T rg4 = rG(buildComplex(-1, 1), buildComplex(0, 1), buildComplex(0));
336 check(0.44660591677018, 0.70768352357515, rg4, 1.0e-13);
337
338 T rg5 = rG(buildComplex(0, -1), buildComplex(-1, 1), buildComplex(0, 1));
339 check(0.36023392184473, 0.40348623401722, rg5, 1.0e-13);
340
341 T rg6 = rG(buildComplex(0), buildComplex(0.0796), buildComplex(4));
342 check( 1.0284758090288, 0.0, rg6, 1.0e-13);
343
344 }
345
346 @Test
347 public void testAlternateRG() {
348 RandomGenerator random = new Well19937c(0xa2946e4a55d133a6l);
349 for (int i = 0; i < 10000; ++i) {
350 T x = buildComplex(random.nextDouble() * 3);
351 T y = buildComplex(random.nextDouble() * 3);
352 T z = buildComplex(random.nextDouble() * 3);
353 assertEquals(0.0, rG(x, y, z).subtract(rgAlternateImplementation(x, y, z)).norm(), 2.0e-15);
354 }
355 }
356
357 @Test
358 public void testRgBuggySquareRoot() {
359
360
361 T x = buildComplex(FastMath.scalb(7745000, -24), -0.5625);
362 T y = buildComplex(-0.3125, -0.6875);
363 T z = buildComplex( 0.9375, 0.25);
364
365
366 assertEquals(0.0, rG(x, y, z). subtract(rgAlternateImplementation(x, y, z)).norm(), 2.1e-16);
367 assertEquals(0.0, buggyRG(x, y, z).subtract(rgAlternateImplementation(x, y, z)).norm(), 2.0e-16);
368
369
370
371
372 x = buildComplex(FastMath.scalb(7744999, -24), -0.5625);
373 assertEquals(0.0, rG(x, y, z). subtract(rgAlternateImplementation(x, y, z)).norm(), 2.5e-16);
374 assertEquals(0.75258, buggyRG(x, y, z).subtract(rgAlternateImplementation(x, y, z)).norm(), 1.0e-5);
375
376 }
377
378 private T buggyRG(final T x, final T y, final T z) {
379 final T termF = new RfFieldDuplication<>(x, y, z).integral().multiply(z);
380 final T termD = x.subtract(z).multiply(y.subtract(z)).multiply(new RdFieldDuplication<>(x, y, z).integral()).divide(3);
381 final T termS = x.multiply(y).divide(z).sqrt();
382 return termF.subtract(termD).add(termS).multiply(0.5);
383 }
384
385 private T rgAlternateImplementation(final T x, final T y, final T z) {
386
387 return d(x, y, z).add(d(y, z, x)).add(d(z, x, y)).divide(6);
388 }
389
390 private T d(final T u, final T v, final T w) {
391 return u.isZero() ? u : u.multiply(v.add(w)).multiply(new RdFieldDuplication<>(v, w, u).integral());
392 }
393
394 }