Elliptic integrals

The MeanFieldHomogenization.Elliptic submodule provides type-generic Legendre and Carlson integrals that flow end-to-end through automatic differentiation (ForwardDiff.Dual), arbitrary precision (BigFloat) and symbolic scalars (SymPy.Sym, Symbolics.Num).

Legendre elliptic integrals

Throughout the submodule, m denotes the parameter (not the modulus), with the convention   and   .

Complete integrals

Incomplete integrals

Both reduce to the complete integrals at  .

Arithmetic–geometric mean (AGM)

For any Number subtype other than Float64, the complete integrals are evaluated using the classical AGM recursion (Abramowitz & Stegun 17.6, NIST DLMF 19.8). Starting from

the sequence

converges quadratically, and

Eight to twelve iterations are typically enough to reach Float64 precision; BigFloat needs only a few more.

Carlson symmetric forms

Incomplete integrals are delegated to Carlson's symmetric integrals (Carlson 1995, Numerical computation of real or complex elliptic integrals):

With these two primitives one recovers the Legendre integrals via

R_F and R_D are implemented by the duplication theorem (Carlson 1995, §2): each iteration halves the relative spread of until their ratio approaches unity, at which point a fifth-order Taylor series in the Carlson invariants (and for ) is used. The recursion is arithmetic-only — no branch cuts, no transcendentals — so it extends to any Number subtype.

Dispatch table

Scalar type	Backend
Float64	Elliptic.jl (GSL C binding, fastest)
ForwardDiff.Dual	pure-Julia AGM / Carlson (derivatives work)
BigFloat, generic	AGM / Carlson
SymPy.Sym	sympy.elliptic_{k,e,f} via the SymPy weak ext
Symbolics.Num	AGM / Carlson (verbose; use simplify if needed)

Why a SymPy weak extension?

The AGM unrolls ~60 nested sqrt(a*b) operations. On a SymPy.Sym input this builds a deeply nested symbolic tree that overflows SymPy's pretty-printer. The weak extension MeanFieldHomogenizationSymPyExt (loaded automatically whenever SymPy is loaded alongside MeanFieldHomogenization) routes ell_K, ell_E, ell_F on Sym arguments directly to sympy.elliptic_{k,e,f}, returning the native closed form instead.

Special cases

•  :   ;   .

•  :     (logarithmic divergence);  . The Float64 fast path through Elliptic.jl throws at  .

•  : incomplete integrals coincide with complete ones.

References

• M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, §17.6, Dover 1972.

• NIST Digital Library of Mathematical Functions, §19.

• B.C. Carlson, Numerical computation of real or complex elliptic integrals, Numerical Algorithms 10 (1995) 13–26.

• B.C. Carlson and E.M. Notis, Algorithms for incomplete elliptic integrals, ACM Transactions on Mathematical Software 7 (1981) 398–403 — public-domain SLATEC routines DRF / DRD.