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