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.CalculusFieldElement;
20 import org.hipparchus.complex.Complex;
21 import org.hipparchus.complex.FieldComplex;
22 import org.hipparchus.util.FastMath;
23
24 /** Elliptic integrals in Carlson symmetric form.
25 * <p>
26 * This utility class computes the various symmetric elliptic
27 * integrals defined as:
28 * \[
29 * \left\{\begin{align}
30 * R_F(x,y,z) &= \frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{s(t)}\\
31 * R_J(x,y,z,p) &= \frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{s(t)(t+p)}\\
32 * R_G(x,y,z) &= \frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
33 \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t\\
34 * R_D(x,y,z) &= R_J(x,y,z,z)\\
35 * R_C(x,y) &= R_F(x,y,y)
36 * \end{align}\right.
37 * \]
38 * </p>
39 * <p>
40 * where
41 * \[
42 * s(t) = \sqrt{t+x}\sqrt{t+y}\sqrt{t+z}
43 * \]
44 * </p>
45 * <p>
46 * The algorithms used are based on the duplication method as described in
47 * B. C. Carlson 1995 paper "Numerical computation of real or complex
48 * elliptic integrals", with the improvements described in the appendix
49 * of B. C. Carlson and James FitzSimons 2000 paper "Reduction theorems
50 * for elliptic integrands with the square root of two quadratic factors".
51 * They are also described in <a href="https://dlmf.nist.gov/19.36#i">section 19.36(i)</a>
52 * of Digital Library of Mathematical Functions.
53 * </p>
54 * <p>
55 * <em>
56 * Beware that when computing elliptic integrals in the complex plane,
57 * many issues arise due to branch cuts. See the
58 * <a href="https://www.hipparchus.org/hipparchus-core/special.html#Elliptic_functions_and_integrals">user guide</a>
59 * for a thorough explanation.
60 * </em>
61 * </p>
62 * @since 2.0
63 */
64 public class CarlsonEllipticIntegral {
65
66 /** Private constructor for a utility class.
67 */
68 private CarlsonEllipticIntegral() {
69 }
70
71 /** Compute Carlson elliptic integral R<sub>C</sub>.
72 * <p>
73 * The Carlson elliptic integral R<sub>C</sub>is defined as
74 * \[
75 * R_C(x,y,z)=R_F(x,y,y)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}(t+y)}
76 * \]
77 * </p>
78 * @param x first symmetric variable of the integral
79 * @param y second symmetric variable of the integral
80 * @return Carlson elliptic integral R<sub>C</sub>
81 */
82 public static double rC(final double x, final double y) {
83 if (y < 0) {
84 // y is on the branch cut, we must use a transformation to get the Cauchy principal value
85 // see equation 2.14 in Carlson[1995]
86 final double xMy = x - y;
87 return FastMath.sqrt(x / xMy) * new RcRealDuplication(xMy, -y).integral();
88 } else {
89 return new RcRealDuplication(x, y).integral();
90 }
91 }
92
93 /** Compute Carlson elliptic integral R<sub>C</sub>.
94 * <p>
95 * The Carlson elliptic integral R<sub>C</sub>is defined as
96 * \[
97 * R_C(x,y,z)=R_F(x,y,y)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}(t+y)}
98 * \]
99 * </p>
100 * @param x first symmetric variable of the integral
101 * @param y second symmetric variable of the integral
102 * @param <T> type of the field elements
103 * @return Carlson elliptic integral R<sub>C</sub>
104 */
105 public static <T extends CalculusFieldElement<T>> T rC(final T x, final T y) {
106 if (y.getReal() < 0) {
107 // y is on the branch cut, we must use a transformation to get the Cauchy principal value
108 // see equation 2.14 in Carlson[1995]
109 final T xMy = x.subtract(y);
110 return FastMath.sqrt(x.divide(xMy)).multiply(new RcFieldDuplication<>(xMy, y.negate()).integral());
111 } else {
112 return new RcFieldDuplication<>(x, y).integral();
113 }
114 }
115
116 /** Compute Carlson elliptic integral R<sub>C</sub>.
117 * <p>
118 * The Carlson elliptic integral R<sub>C</sub>is defined as
119 * \[
120 * R_C(x,y,z)=R_F(x,y,y)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}(t+y)}
121 * \]
122 * </p>
123 * @param x first symmetric variable of the integral
124 * @param y second symmetric variable of the integral
125 * @return Carlson elliptic integral R<sub>C</sub>
126 */
127 public static Complex rC(final Complex x, final Complex y) {
128 if (y.getImaginaryPart() == 0 && y.getRealPart() < 0) {
129 // y is on the branch cut, we must use a transformation to get the Cauchy principal value
130 // see equation 2.14 in Carlson[1995]
131 final Complex xMy = x.subtract(y);
132 return FastMath.sqrt(x.divide(xMy)).multiply(new RcFieldDuplication<>(xMy, y.negate()).integral());
133 } else {
134 return new RcFieldDuplication<>(x, y).integral();
135 }
136 }
137
138 /** Compute Carlson elliptic integral R<sub>C</sub>.
139 * <p>
140 * The Carlson elliptic integral R<sub>C</sub>is defined as
141 * \[
142 * R_C(x,y,z)=R_F(x,y,y)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}(t+y)}
143 * \]
144 * </p>
145 * @param x first symmetric variable of the integral
146 * @param y second symmetric variable of the integral
147 * @param <T> type of the field elements
148 * @return Carlson elliptic integral R<sub>C</sub>
149 */
150 public static <T extends CalculusFieldElement<T>> FieldComplex<T> rC(final FieldComplex<T> x, final FieldComplex<T> y) {
151 if (y.getImaginaryPart().isZero() && y.getRealPart().getReal() < 0) {
152 // y is on the branch cut, we must use a transformation to get the Cauchy principal value
153 // see equation 2.14 in Carlson[1995]
154 final FieldComplex<T> xMy = x.subtract(y);
155 return FastMath.sqrt(x.divide(xMy)).multiply(new RcFieldDuplication<>(xMy, y.negate()).integral());
156 } else {
157 return new RcFieldDuplication<>(x, y).integral();
158 }
159 }
160
161 /** Compute Carlson elliptic integral R<sub>F</sub>.
162 * <p>
163 * The Carlson elliptic integral R<sub>F</sub> is defined as
164 * \[
165 * R_F(x,y,z)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}}
166 * \]
167 * </p>
168 * @param x first symmetric variable of the integral
169 * @param y second symmetric variable of the integral
170 * @param z third symmetric variable of the integral
171 * @return Carlson elliptic integral R<sub>F</sub>
172 */
173 public static double rF(final double x, final double y, final double z) {
174 return new RfRealDuplication(x, y, z).integral();
175 }
176
177 /** Compute Carlson elliptic integral R<sub>F</sub>.
178 * <p>
179 * The Carlson elliptic integral R<sub>F</sub> is defined as
180 * \[
181 * R_F(x,y,z)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}}
182 * \]
183 * </p>
184 * @param x first symmetric variable of the integral
185 * @param y second symmetric variable of the integral
186 * @param z third symmetric variable of the integral
187 * @param <T> type of the field elements
188 * @return Carlson elliptic integral R<sub>F</sub>
189 */
190 public static <T extends CalculusFieldElement<T>> T rF(final T x, final T y, final T z) {
191 return new RfFieldDuplication<>(x, y, z).integral();
192 }
193
194 /** Compute Carlson elliptic integral R<sub>F</sub>.
195 * <p>
196 * The Carlson elliptic integral R<sub>F</sub> is defined as
197 * \[
198 * R_F(x,y,z)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}}
199 * \]
200 * </p>
201 * @param x first symmetric variable of the integral
202 * @param y second symmetric variable of the integral
203 * @param z third symmetric variable of the integral
204 * @return Carlson elliptic integral R<sub>F</sub>
205 */
206 public static Complex rF(final Complex x, final Complex y, final Complex z) {
207 return new RfFieldDuplication<>(x, y, z).integral();
208 }
209
210 /** Compute Carlson elliptic integral R<sub>F</sub>.
211 * <p>
212 * The Carlson elliptic integral R<sub>F</sub> is defined as
213 * \[
214 * R_F(x,y,z)=\frac{1}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}}
215 * \]
216 * </p>
217 * @param x first symmetric variable of the integral
218 * @param y second symmetric variable of the integral
219 * @param z third symmetric variable of the integral
220 * @param <T> type of the field elements
221 * @return Carlson elliptic integral R<sub>F</sub>
222 */
223 public static <T extends CalculusFieldElement<T>> FieldComplex<T> rF(final FieldComplex<T> x, final FieldComplex<T> y, final FieldComplex<T> z) {
224 return new RfFieldDuplication<>(x, y, z).integral();
225 }
226
227 /** Compute Carlson elliptic integral R<sub>J</sub>.
228 * <p>
229 * The Carlson elliptic integral R<sub>J</sub> is defined as
230 * \[
231 * R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
232 * \]
233 * </p>
234 * @param x first symmetric variable of the integral
235 * @param y second symmetric variable of the integral
236 * @param z third symmetric variable of the integral
237 * @param p fourth <em>not</em> symmetric variable of the integral
238 * @return Carlson elliptic integral R<sub>J</sub>
239 */
240 public static double rJ(final double x, final double y, final double z, final double p) {
241 final double delta = (p - x) * (p - y) * (p - z);
242 return rJ(x, y, z, p, delta);
243 }
244
245 /** Compute Carlson elliptic integral R<sub>J</sub>.
246 * <p>
247 * The Carlson elliptic integral R<sub>J</sub> is defined as
248 * \[
249 * R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
250 * \]
251 * </p>
252 * @param x first symmetric variable of the integral
253 * @param y second symmetric variable of the integral
254 * @param z third symmetric variable of the integral
255 * @param p fourth <em>not</em> symmetric variable of the integral
256 * @param delta precomputed value of (p-x)(p-y)(p-z)
257 * @return Carlson elliptic integral R<sub>J</sub>
258 */
259 public static double rJ(final double x, final double y, final double z, final double p, final double delta) {
260 return new RjRealDuplication(x, y, z, p, delta).integral();
261 }
262
263 /** Compute Carlson elliptic integral R<sub>J</sub>.
264 * <p>
265 * The Carlson elliptic integral R<sub>J</sub> is defined as
266 * \[
267 * R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
268 * \]
269 * </p>
270 * @param x first symmetric variable of the integral
271 * @param y second symmetric variable of the integral
272 * @param z third symmetric variable of the integral
273 * @param p fourth <em>not</em> symmetric variable of the integral
274 * @param <T> type of the field elements
275 * @return Carlson elliptic integral R<sub>J</sub>
276 */
277 public static <T extends CalculusFieldElement<T>> T rJ(final T x, final T y, final T z, final T p) {
278 final T delta = p.subtract(x).multiply(p.subtract(y)).multiply(p.subtract(z));
279 return new RjFieldDuplication<>(x, y, z, p, delta). integral();
280 }
281
282 /** Compute Carlson elliptic integral R<sub>J</sub>.
283 * <p>
284 * The Carlson elliptic integral R<sub>J</sub> is defined as
285 * \[
286 * R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
287 * \]
288 * </p>
289 * @param x first symmetric variable of the integral
290 * @param y second symmetric variable of the integral
291 * @param z third symmetric variable of the integral
292 * @param p fourth <em>not</em> symmetric variable of the integral
293 * @param delta precomputed value of (p-x)(p-y)(p-z)
294 * @param <T> type of the field elements
295 * @return Carlson elliptic integral R<sub>J</sub>
296 */
297 public static <T extends CalculusFieldElement<T>> T rJ(final T x, final T y, final T z, final T p, final T delta) {
298 return new RjFieldDuplication<>(x, y, z, p, delta).integral();
299 }
300
301 /** Compute Carlson elliptic integral R<sub>J</sub>.
302 * <p>
303 * The Carlson elliptic integral R<sub>J</sub> is defined as
304 * \[
305 * R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
306 * \]
307 * </p>
308 * @param x first symmetric variable of the integral
309 * @param y second symmetric variable of the integral
310 * @param z third symmetric variable of the integral
311 * @param p fourth <em>not</em> symmetric variable of the integral
312 * @return Carlson elliptic integral R<sub>J</sub>
313 */
314 public static Complex rJ(final Complex x, final Complex y, final Complex z, final Complex p) {
315 final Complex delta = p.subtract(x).multiply(p.subtract(y)).multiply(p.subtract(z));
316 return new RjFieldDuplication<>(x, y, z, p, delta).integral();
317 }
318
319 /** Compute Carlson elliptic integral R<sub>J</sub>.
320 * <p>
321 * The Carlson elliptic integral R<sub>J</sub> is defined as
322 * \[
323 * R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
324 * \]
325 * </p>
326 * @param x first symmetric variable of the integral
327 * @param y second symmetric variable of the integral
328 * @param z third symmetric variable of the integral
329 * @param p fourth <em>not</em> symmetric variable of the integral
330 * @param delta precomputed value of (p-x)(p-y)(p-z)
331 * @return Carlson elliptic integral R<sub>J</sub>
332 */
333 public static Complex rJ(final Complex x, final Complex y, final Complex z, final Complex p, final Complex delta) {
334 return new RjFieldDuplication<>(x, y, z, p, delta).integral();
335 }
336
337 /** Compute Carlson elliptic integral R<sub>J</sub>.
338 * <p>
339 * The Carlson elliptic integral R<sub>J</sub> is defined as
340 * \[
341 * R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
342 * \]
343 * </p>
344 * @param x first symmetric variable of the integral
345 * @param y second symmetric variable of the integral
346 * @param z third symmetric variable of the integral
347 * @param p fourth <em>not</em> symmetric variable of the integral
348 * @param <T> type of the field elements
349 * @return Carlson elliptic integral R<sub>J</sub>
350 */
351 public static <T extends CalculusFieldElement<T>> FieldComplex<T> rJ(final FieldComplex<T> x, final FieldComplex<T> y,
352 final FieldComplex<T> z, final FieldComplex<T> p) {
353 final FieldComplex<T> delta = p.subtract(x).multiply(p.subtract(y)).multiply(p.subtract(z));
354 return new RjFieldDuplication<>(x, y, z, p, delta).integral();
355 }
356
357 /** Compute Carlson elliptic integral R<sub>J</sub>.
358 * <p>
359 * The Carlson elliptic integral R<sub>J</sub> is defined as
360 * \[
361 * R_J(x,y,z,p)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+p)}
362 * \]
363 * </p>
364 * @param x first symmetric variable of the integral
365 * @param y second symmetric variable of the integral
366 * @param z third symmetric variable of the integral
367 * @param p fourth <em>not</em> symmetric variable of the integral
368 * @param delta precomputed value of (p-x)(p-y)(p-z)
369 * @param <T> type of the field elements
370 * @return Carlson elliptic integral R<sub>J</sub>
371 */
372 public static <T extends CalculusFieldElement<T>> FieldComplex<T> rJ(final FieldComplex<T> x, final FieldComplex<T> y,
373 final FieldComplex<T> z, final FieldComplex<T> p,
374 final FieldComplex<T> delta) {
375 return new RjFieldDuplication<>(x, y, z, p, delta).integral();
376 }
377
378 /** Compute Carlson elliptic integral R<sub>D</sub>.
379 * <p>
380 * The Carlson elliptic integral R<sub>D</sub> is defined as
381 * \[
382 * R_D(x,y,z)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+z)}
383 * \]
384 * </p>
385 * @param x first symmetric variable of the integral
386 * @param y second symmetric variable of the integral
387 * @param z third symmetric variable of the integral
388 * @return Carlson elliptic integral R<sub>D</sub>
389 */
390 public static double rD(final double x, final double y, final double z) {
391 return new RdRealDuplication(x, y, z).integral();
392 }
393
394 /** Compute Carlson elliptic integral R<sub>D</sub>.
395 * <p>
396 * The Carlson elliptic integral R<sub>D</sub> is defined as
397 * \[
398 * R_D(x,y,z)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+z)}
399 * \]
400 * </p>
401 * @param x first symmetric variable of the integral
402 * @param y second symmetric variable of the integral
403 * @param z third symmetric variable of the integral
404 * @param <T> type of the field elements
405 * @return Carlson elliptic integral R<sub>D</sub>
406 */
407 public static <T extends CalculusFieldElement<T>> T rD(final T x, final T y, final T z) {
408 return new RdFieldDuplication<>(x, y, z).integral();
409 }
410
411 /** Compute Carlson elliptic integral R<sub>D</sub>.
412 * <p>
413 * The Carlson elliptic integral R<sub>D</sub> is defined as
414 * \[
415 * R_D(x,y,z)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+z)}
416 * \]
417 * </p>
418 * @param x first symmetric variable of the integral
419 * @param y second symmetric variable of the integral
420 * @param z third symmetric variable of the integral
421 * @return Carlson elliptic integral R<sub>D</sub>
422 */
423 public static Complex rD(final Complex x, final Complex y, final Complex z) {
424 return new RdFieldDuplication<>(x, y, z).integral();
425 }
426
427 /** Compute Carlson elliptic integral R<sub>D</sub>.
428 * <p>
429 * The Carlson elliptic integral R<sub>D</sub> is defined as
430 * \[
431 * R_D(x,y,z)=\frac{3}{2}\int_{0}^{\infty}\frac{\mathrm{d}t}{\sqrt{t+x}\sqrt{t+y}\sqrt{t+z}(t+z)}
432 * \]
433 * </p>
434 * @param x first symmetric variable of the integral
435 * @param y second symmetric variable of the integral
436 * @param z third symmetric variable of the integral
437 * @param <T> type of the field elements
438 * @return Carlson elliptic integral R<sub>D</sub>
439 */
440 public static <T extends CalculusFieldElement<T>> FieldComplex<T> rD(final FieldComplex<T> x, final FieldComplex<T> y,
441 final FieldComplex<T> z) {
442 return new RdFieldDuplication<>(x, y, z).integral();
443 }
444
445 /** Compute Carlson elliptic integral R<sub>G</sub>.
446 * <p>
447 * The Carlson elliptic integral R<sub>G</sub>is defined as
448 * \[
449 * R_{G}(x,y,z)=\frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
450 * \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t
451 * \]
452 * </p>
453 * @param x first symmetric variable of the integral
454 * @param y second symmetric variable of the integral
455 * @param z second symmetric variable of the integral
456 * @return Carlson elliptic integral R<sub>G</sub>
457 */
458 public static double rG(final double x, final double y, final double z) {
459 return generalComputeRg(x, y, z);
460 }
461
462 /** Compute Carlson elliptic integral R<sub>G</sub>.
463 * <p>
464 * The Carlson elliptic integral R<sub>G</sub>is defined as
465 * \[
466 * R_{G}(x,y,z)=\frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
467 * \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t
468 * \]
469 * </p>
470 * @param x first symmetric variable of the integral
471 * @param y second symmetric variable of the integral
472 * @param z second symmetric variable of the integral
473 * @param <T> type of the field elements
474 * @return Carlson elliptic integral R<sub>G</sub>
475 */
476 public static <T extends CalculusFieldElement<T>> T rG(final T x, final T y, final T z) {
477 return generalComputeRg(x, y, z);
478 }
479
480 /** Compute Carlson elliptic integral R<sub>G</sub>.
481 * <p>
482 * The Carlson elliptic integral R<sub>G</sub>is defined as
483 * \[
484 * R_{G}(x,y,z)=\frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
485 * \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t
486 * \]
487 * </p>
488 * @param x first symmetric variable of the integral
489 * @param y second symmetric variable of the integral
490 * @param z second symmetric variable of the integral
491 * @return Carlson elliptic integral R<sub>G</sub>
492 */
493 public static Complex rG(final Complex x, final Complex y, final Complex z) {
494 return generalComputeRg(x, y, z);
495 }
496
497 /** Compute Carlson elliptic integral R<sub>G</sub>.
498 * <p>
499 * The Carlson elliptic integral R<sub>G</sub>is defined as
500 * \[
501 * R_{G}(x,y,z)=\frac{1}{4}\int_{0}^{\infty}\frac{1}{s(t)}
502 * \left(\frac{x}{t+x}+\frac{y}{t+y}+\frac{z}{t+z}\right)t\mathrm{d}t
503 * \]
504 * </p>
505 * @param x first symmetric variable of the integral
506 * @param y second symmetric variable of the integral
507 * @param z second symmetric variable of the integral
508 * @param <T> type of the field elements
509 * @return Carlson elliptic integral R<sub>G</sub>
510 */
511 public static <T extends CalculusFieldElement<T>> FieldComplex<T> rG(final FieldComplex<T> x,
512 final FieldComplex<T> y,
513 final FieldComplex<T> z) {
514 return generalComputeRg(x, y, z);
515 }
516
517 /** Compute Carlson elliptic integral R<sub>G</sub> in the general case.
518 * @param x first symmetric variable of the integral
519 * @param y second symmetric variable of the integral
520 * @param z third symmetric variable of the integral
521 * @return Carlson elliptic integral R<sub>G</sub>
522 */
523 private static double generalComputeRg(final double x, final double y, final double z) {
524 // permute parameters if needed to avoid cancellations
525 if (x <= y) {
526 if (y <= z) {
527 // x ≤ y ≤ z
528 return permutedComputeRg(x, z, y);
529 } else if (x <= z) {
530 // x ≤ z < y
531 return permutedComputeRg(x, y, z);
532 } else {
533 // z < x ≤ y
534 return permutedComputeRg(z, y, x);
535 }
536 } else if (x <= z) {
537 // y < x ≤ z
538 return permutedComputeRg(y, z, x);
539 } else if (y <= z) {
540 // y ≤ z < x
541 return permutedComputeRg(y, x, z);
542 } else {
543 // z < y < x
544 return permutedComputeRg(z, x, y);
545 }
546 }
547
548 /** Compute Carlson elliptic integral R<sub>G</sub> in the general case.
549 * @param x first symmetric variable of the integral
550 * @param y second symmetric variable of the integral
551 * @param z third symmetric variable of the integral
552 * @param <T> type of the field elements (really {@link Complex} or {@link FieldComplex})
553 * @return Carlson elliptic integral R<sub>G</sub>
554 */
555 private static <T extends CalculusFieldElement<T>> T generalComputeRg(final T x, final T y, final T z) {
556 // permute parameters if needed to avoid cancellations
557 final double xR = x.getReal();
558 final double yR = y.getReal();
559 final double zR = z.getReal();
560 if (xR <= yR) {
561 if (yR <= zR) {
562 // x ≤ y ≤ z
563 return permutedComputeRg(x, z, y);
564 } else if (xR <= zR) {
565 // x ≤ z < y
566 return permutedComputeRg(x, y, z);
567 } else {
568 // z < x ≤ y
569 return permutedComputeRg(z, y, x);
570 }
571 } else if (xR <= zR) {
572 // y < x ≤ z
573 return permutedComputeRg(y, z, x);
574 } else if (yR <= zR) {
575 // y ≤ z < x
576 return permutedComputeRg(y, x, z);
577 } else {
578 // z < y < x
579 return permutedComputeRg(z, x, y);
580 }
581 }
582
583 /** Compute Carlson elliptic integral R<sub>G</sub> with already permuted variables to avoid cancellations.
584 * @param x first symmetric variable of the integral
585 * @param y second symmetric variable of the integral
586 * @param z third symmetric variable of the integral
587 * @return Carlson elliptic integral R<sub>G</sub>
588 */
589 private static double permutedComputeRg(final double x, final double y, final double z) {
590 // permute parameters if needed to avoid divisions by zero
591 if (z == 0) {
592 return x == 0 ? safeComputeRg(z, x, y) : safeComputeRg(y, z, x);
593 } else {
594 return safeComputeRg(x, y, z);
595 }
596 }
597
598 /** Compute Carlson elliptic integral R<sub>G</sub> with already permuted variables to avoid cancellations.
599 * @param x first symmetric variable of the integral
600 * @param y second symmetric variable of the integral
601 * @param z third symmetric variable of the integral
602 * @param <T> type of the field elements (really {@link Complex} or {@link FieldComplex})
603 * @return Carlson elliptic integral R<sub>G</sub>
604 */
605 private static <T extends CalculusFieldElement<T>> T permutedComputeRg(final T x, final T y, final T z) {
606 // permute parameters if needed to avoid divisions by zero
607 if (z.isZero()) {
608 return x.isZero() ? safeComputeRg(z, x, y) : safeComputeRg(y, z, x);
609 } else {
610 return safeComputeRg(x, y, z);
611 }
612 }
613
614 /** Compute Carlson elliptic integral R<sub>G</sub> with non-zero third variable.
615 * @param x first symmetric variable of the integral
616 * @param y second symmetric variable of the integral
617 * @param z third symmetric variable of the integral
618 * @see <a href="https://dlmf.nist.gov/19.21#E10">Digital Library of Mathematical Functions, equation 19.21.10</a>
619 * @return Carlson elliptic integral R<sub>G</sub>
620 */
621 private static double safeComputeRg(final double x, final double y, final double z) {
622
623 // contribution of the R_F integral
624 final double termF = new RfRealDuplication(x, y, z).integral() * z;
625
626 // contribution of the R_D integral
627 final double termD = (x - z) * (y - z) * new RdRealDuplication(x, y, z).integral() / 3;
628
629 // contribution of the square roots
630 final double termS = FastMath.sqrt(x * y / z);
631
632 // equation 19.21.10
633 return (termF - termD + termS) * 0.5;
634
635 }
636
637 /** Compute Carlson elliptic integral R<sub>G</sub> with non-zero third variable.
638 * @param x first symmetric variable of the integral
639 * @param y second symmetric variable of the integral
640 * @param z third symmetric variable of the integral
641 * @param <T> type of the field elements (really {@link Complex} or {@link FieldComplex})
642 * @see <a href="https://dlmf.nist.gov/19.21#E10">Digital Library of Mathematical Functions, equation 19.21.10</a>
643 * @return Carlson elliptic integral R<sub>G</sub>
644 */
645 private static <T extends CalculusFieldElement<T>> T safeComputeRg(final T x, final T y, final T z) {
646
647 // contribution of the R_F integral
648 final T termF = new RfFieldDuplication<>(x, y, z).integral().multiply(z);
649
650 // contribution of the R_D integral
651 final T termD = x.subtract(z).multiply(y.subtract(z)).multiply(new RdFieldDuplication<>(x, y, z).integral()).divide(3);
652
653 // contribution of the square roots
654 // BEWARE: this term MUST be computed as √x√y/√z with all square roots selected with positive real part
655 // and NOT as √(xy/z), otherwise sign errors may occur
656 final T termS = x.sqrt().multiply(y.sqrt()).divide(z.sqrt());
657
658 // equation 19.21.10
659 return termF.subtract(termD).add(termS).multiply(0.5);
660
661 }
662
663 }