Computation of elliptic_pi

Principal value of elliptic_pi for $n > 1$

This note covers how Maxima should evaluate the incomplete elliptic integral of the third kind,

$$ \Pi(n;\varphi\thinspace|\thinspace m) = \int_0^{\varphi} \frac{d\theta}{\left(1 - n\sin^2\theta\right)\sqrt{1 - m\sin^2\theta}}, $$

for real $n > 1$, where the integrand has poles on the path of integration. It explains:

• what the existing code returns;
• what Mathematica, Wolfram|Alpha and mpmath return;
• how the Sokhotski–Plemelj formula connects the two;
• how the result was implemented and tested.

Contents

1. The problem
2. What the Carlson code computes
3. The Sokhotski–Plemelj formula
4. Applying it to Π
5. The same result as a contour indentation
6. Regions of the real parameter space
7. Implementation
8. Tests
9. Known limitations
10. References

1. The problem

When $n > 1$, the factor $1 - n\sin^2\theta$ vanishes at

$$ \theta^{\ast} = \arcsin\frac{1}{\sqrt{n}} $$

and at $\pi - \theta^{\ast}$ and their translates by $k\pi$. Once $\varphi > \theta^{\ast}$ the integral does not exist in the ordinary sense. There are two common conventions for assigning it a value:

• The Cauchy principal value (CPV). This is real.
• The analytic continuation from one side of the branch cut $n \in (1, \infty)$. This is complex.

Mathematica, Wolfram|Alpha and mpmath use the analytic continuation. For example,

$$ \Pi(1.5\thinspace|\thinspace 0.5) = -0.45672031345290990\ldots \ -\ 2.72069904635132678\ldots\thinspace i . $$

Before this change, Maxima returned only the real part, $-0.45672\ldots$, which is the CPV.

2. What the Carlson code computes

Maxima evaluates $\Pi$ through Carlson's symmetric integrals:

$$ \Pi(n;\varphi\thinspace|\thinspace m) = \sin\varphi\ R_F\negthinspace\left(c^2,\thinspace 1 - m s^2,\thinspace 1\right) + \frac{n}{3}\thinspace s^3\thinspace R_J\negthinspace\left(c^2,\thinspace 1 - m s^2,\thinspace 1,\thinspace 1 - n s^2\right), $$

where $s = \sin\varphi$ and $c = \cos\varphi$. When $n s^2 > 1$, the fourth argument $p = 1 - n s^2$ of $R_J$ is negative.

bf-rj handles real $x, y, z \ge 0$ with $p < 0$ through Carlson (1995), eq. (33), which is the same as DLMF 19.20.14:

$$ (y + q)\thinspace R_J(x, y, z, -q) = (p - y)\thinspace R_J(x, y, z, p) - 3 R_F(x, y, z) + 3\sqrt{\frac{xyz}{xz + pq}}\ R_C(xz + pq,\thinspace pq), $$

with

$$ p = y + \frac{(z - y)(y - x)}{y + q}, $$

and $x, y, z$ permuted so that $(z - y)(y - x) > 0$. Carlson states that this gives the Cauchy principal value, which is real. So the existing code returns the CPV of $\Pi$ whenever this branch is taken.

The real part was verified independently: a brute-force symmetric cut-out around the pole converges to the same value. For $\Pi(1.5;\thinspace 1.2\thinspace|\thinspace 0.5)$, both give $0.75009395785731911\ldots$.

Carlson (1995) does not say how to get the imaginary part, and it does not treat the complete case $x = 0$. The next sections derive the imaginary part directly from the integral.

3. The Sokhotski–Plemelj formula

Let $a < 0 < b$, and let $f$ be differentiable near $0$ (Hölder continuity is enough). Then

$$ \lim_{\varepsilon \to 0^{+}} \int_a^b \frac{f(x)}{x \pm i\varepsilon}\thinspace dx = \mathrm{PV}\negthinspace\int_a^b \frac{f(x)}{x}\thinspace dx \ \mp\ i\pi f(0). $$

The real part of the limit is the principal value. The imaginary part is $\pi$ times the regular factor of the integrand, evaluated at the singularity.

To see why, split the kernel into real and imaginary parts:

$$ \frac{1}{x \pm i\varepsilon} = \frac{x}{x^2 + \varepsilon^2} \mp i\thinspace\frac{\varepsilon}{x^2 + \varepsilon^2} . $$

• The real part equals $1/x$ except within about $\varepsilon$ of $0$, where it stays bounded and is odd. Integrated against $f$, it tends to the principal value.
• The imaginary part is a peak of width about $\varepsilon$ at $0$, with total area $\int_{-\infty}^{\infty} \varepsilon/(x^2 + \varepsilon^2)\thinspace dx = \pi$. As $\varepsilon \to 0$ it only sees $f(0)$, so it contributes $\mp i\pi f(0)$.

4. Applying it to Π

Write the integrand as $h/g$, where

$$ g(\theta) = 1 - n\sin^2\theta, \qquad h(\theta) = \frac{1}{\sqrt{1 - m\sin^2\theta}} . $$

For $n > 1$, $g$ has simple zeros at $\theta^{\ast}$, $\pi - \theta^{\ast}$, and their translates by $k\pi$. For $m < n$, $h$ is smooth near each of them.

Step 1: move $n$ off the real axis

Replace $n$ by $n - i\varepsilon$:

$$ 1 - (n - i\varepsilon)\sin^2\theta = g(\theta) + i\varepsilon\sin^2\theta . $$

Near the pole, $1 - n\sin^2\theta \approx 0$, so $\sin^2\theta \approx 1/n > 0$. The denominator is therefore $g(\theta) + i\varepsilon'$ with $\varepsilon' = \varepsilon\sin^2\theta \approx \varepsilon/n > 0$, which has the $+i\varepsilon$ form of §3. Away from the poles $g \ne 0$, and the extra term vanishes in the limit.

Step 2: apply the formula at each pole

Take a small interval $(\theta_k - \delta,\ \theta_k + \delta)$ around a pole $\theta_k$, and substitute $u = g(\theta)$, so that $d\theta = du/g'(\theta)$. At $\theta^{\ast}$, $g' < 0$, so $u$ decreases from $b = g(\theta^{\ast} - \delta) > 0$ to $a = g(\theta^{\ast} + \delta) < 0$ as $\theta$ increases. Reversing the limits absorbs the sign:

$$ \int_{\theta^{\ast} - \delta}^{\theta^{\ast} + \delta} \frac{h(\theta)}{g(\theta) + i\varepsilon'}\thinspace d\theta = \int_{a}^{b} \frac{h(\theta(u))}{\lvert g'(\theta(u))\rvert}\cdot\frac{du}{u + i\varepsilon'} . $$

At poles where $g' > 0$, such as $\pi - \theta^{\ast}$, $u$ increases and the same expression results directly. This is now the integral of §3, with $f(u) = h(\theta(u))/\lvert g'(\theta(u))\rvert$ and $f(0) = h(\theta_k)/\lvert g'(\theta_k)\rvert$. The parts of the path away from the poles have no singularity and need nothing special. Adding everything up gives

$$ \Pi(n - i0;\varphi\thinspace|\thinspace m) = \mathrm{PV}\negthinspace\int_0^{\varphi} \frac{h(\theta)}{g(\theta)}\thinspace d\theta \ -\ i\pi \sum_{\theta_k \in (0,\varphi)} \frac{h(\theta_k)}{\lvert g'(\theta_k)\rvert} . $$

Step 3: evaluate at $\theta^{\ast}$

From $\sin\theta^{\ast} = 1/\sqrt{n}$ and $\cos\theta^{\ast} = \sqrt{(n-1)/n}$:

$$ \lvert g'(\theta^{\ast})\rvert = 2n\sin\theta^{\ast}\cos\theta^{\ast} = 2\sqrt{n - 1}, \qquad h(\theta^{\ast}) = \frac{1}{\sqrt{1 - m/n}} = \sqrt{\frac{n}{n - m}} . $$

Each pole therefore contributes

$$ -\thinspace i\pi\thinspace\frac{\sqrt{n/(n - m)}}{2\sqrt{n - 1}} = -\thinspace i\thinspace\frac{\pi}{2}\sqrt{\frac{n}{(n - 1)(n - m)}} . $$

Step 4: count the poles

At every zero of $g$, the values of $\sin^2\theta$ and $\lvert\sin\theta\cos\theta\rvert$ are the same, so every pole contributes the same amount. With

$$ \sigma = \frac{\pi}{2}\sqrt{\frac{n}{(n - 1)(n - m)}}, $$

the imaginary part is $-\sigma$ times the number of poles crossed between $0$ and $\varphi$, counted with the sign of $\varphi$. Integrating over negative $\varphi$ runs the path backwards, which is why $\Pi$ is odd in $\varphi$.

Check. For $n = 1.5$ and $m = 0.5$:

$$ \sigma = \frac{\pi}{2}\sqrt{\frac{1.5}{0.5 \cdot 1.0}} = \frac{\pi}{2}\sqrt{3} = 2.72069904635132678\ldots $$

This matches the imaginary part of $\Pi(1.5\thinspace|\thinspace 0.5)$ in §1.

5. The same result as a contour indentation

With $n \to n - i\varepsilon$, the pole moves to

$$ \theta \approx \theta^{\ast} + i\varepsilon\thinspace\frac{\sin^2\theta^{\ast}}{\lvert g'(\theta^{\ast})\rvert}, $$

just above the real axis. The path along the real axis therefore passes below the pole. That is equivalent to indenting the path with a small semicircle below the pole. The semicircle goes counterclockwise around the pole, so it contributes $+i\pi$ times the residue:

$$ \mathop{\mathrm{Res}}_{\theta = \theta^{\ast}} \frac{h}{g} = \frac{h(\theta^{\ast})}{g'(\theta^{\ast})} = -\frac{1}{2}\sqrt{\frac{n}{(n - 1)(n - m)}}, \qquad i\pi \cdot \mathop{\mathrm{Res}} = -\thinspace i\sigma . $$

The result is the same as in §4.

Taking $n + i0$ instead puts the pole below the path and gives the complex conjugate. The CPV is the average of the two one-sided limits. Mathematica, Wolfram|Alpha and mpmath all use the $n - i0$ side.

6. Regions of the real parameter space

The Sokhotski–Plemelj term is needed exactly where bf-rj takes its CPV branch, that is, where $p = 1 - n\sin^2\theta < 0$ while $y = 1 - m\sin^2\theta \ge 0$. Everywhere else, the general duplication algorithm bf-rj1 handles the continuation and already returns the $n - i0$ value.

Let $\theta_b = \arcsin(1/\sqrt{m})$ be the branch point of $\sqrt{1 - m\sin^2\theta}$, defined when $m > 1$.

Region	Behaviour	Action
$n < 1$	No pole on the path	Nothing to add
$n = 1$	Double zero at $\pi/2$, so the integral diverges for $\varphi \ge \pi/2$	None; no principal value exists
$n > 1,\ m \le 1$	Poles at $\pm\theta^{\ast} + k\pi$	Add $-i\sigma$ for each pole crossed
$n > 1,\ 1 < m < n$, and $\theta^{\ast} < \lvert\varphi_r\rvert \le \theta_b$	Pole crossed before the branch point	Add $-i\sigma\thinspace\mathrm{sgn}(\varphi_r)$
$n > 1,\ 1 < m < n$, and $\lvert\varphi_r\rvert > \theta_b$	Past the branch point	Nothing to add; bf-rj1 already gives the $n - i0$ value
$n > 1,\ m > n$	Branch point comes before the pole, where $h$ is imaginary	Nothing to add; bf-rj1 already gives the $n - i0$ value
$m = n$	Pole and branch point coincide, giving $g^{-3/2}$, which is not integrable	None; only a Hadamard finite part exists
$n\sin^2\varphi = 1$	Endpoint on the pole, giving a logarithmic divergence	None

Here $\varphi_r$ is $\varphi$ reduced to $[-\pi/2, \pi/2]$, defined in §7.

The cycles of $\pi$ need separate treatment. For $m \le 1$, the complete integral is evaluated through the CPV branch, so each cycle needs $2 \times (-i\sigma)$. For $m > 1$, $1 - m < 0$, so the complete integral goes through bf-rj1 and already includes the term.

These rules were found by comparing the uncorrected code with mpmath over a grid of $n$, $m$ and $\varphi$. The table follows the rule above.

7. Implementation

bf-elliptic-pi reduces $\varphi$ by multiples of $\pi$,

$$ \Pi(n;\thinspace\varphi_r + k\pi\thinspace|\thinspace m) = 2k\thinspace\Pi(n\thinspace|\thinspace m) + \Pi(n;\thinspace\varphi_r\thinspace|\thinspace m), $$

and sums the CPV pieces. The new function bf-elliptic-pi-residue, in the bigfloat package in src/ellipt.lisp, computes the imaginary correction:

(defun bf-elliptic-pi-residue (n phi m)
  (if (and (realp n) (realp phi) (realp m) (> n 1))
      (let* ((bf-pi (%pi phi))
             (cycles (round phi bf-pi))
             (rem (- phi (* cycles bf-pi)))
             (s2 (expt (sin rem) 2))
             (crossings (+ (if (<= m 1)
                               (* 2 cycles)
                               0)
                           (if (and (> (* n s2) 1) (<= (* m s2) 1))
                               (signum rem)
                               0))))
        (if (zerop crossings)
            0
            (* (complex 0 -1) crossings (/ bf-pi 2)
               (sqrt (/ n (* (- n 1) (- n m)))))))
      0))

The elliptic_pi simplifier adds this correction in its real-float and real-bigfloat branches. The complex branches are unchanged.

((float-numerical-eval-p n phi m)
 (add (elliptic-pi ($float n) ($float phi) ($float m))
      (to (bigfloat::bf-elliptic-pi-residue ($float n) ($float phi) ($float m)))))
...
((bigfloat-numerical-eval-p n phi m)
 (to (bigfloat:+ (bigfloat::bf-elliptic-pi (bigfloat:to n)
                                           (bigfloat:to phi)
                                           (bigfloat:to m))
                 (bigfloat::bf-elliptic-pi-residue (bigfloat:to n)
                                                   (bigfloat:to phi)
                                                   (bigfloat:to m)))))

8. Tests

New tests in tests/rtest_elliptic.mac check each case in both double float and bigfloat at fpprec: 32, against 50-digit references from mpmath's ellippi.

• All inputs are dyadic rationals, so the float and bigfloat arguments are exact.
• Each tolerance is the error observed in testing, rounded up to five significant digits.

$(n,\ \varphi,\ m)$	What it covers	Reference value
$(1.5,\ \pi/2,\ 0.5)$	Complete integral, $0 < m < 1$	$-0.456720313452909897 - 2.720699046351326776\thinspace i$
$(1.5,\ 1.25,\ 0.5)$	Pole crossed	$0.545248053746195541 - 2.720699046351326776\thinspace i$
$(1.5,\ 0.625,\ 0.5)$	Pole not reached; imaginary part is exactly $0$	$0.827387372389447886$
$(1.5,\ -1.25,\ 0.5)$	Odd in $\varphi$	$-0.545248053746195541 + 2.720699046351326776\thinspace i$
$(1.5,\ 2.5,\ 0.5)$	Reduction by $\pi$	$-1.779327655628395173 - 5.441398092702653552\thinspace i$
$(1.5,\ 4.0,\ 0.5)$	Reduction by $\pi$	$0.871635485678527766 - 5.441398092702653552\thinspace i$
$(2,\ 1,\ -0.75)$	$m < 0$	$0.779446017152295514 - 1.339579608485734033\thinspace i$
$(1.5,\ \pi/2,\ -0.75)$	Complete integral, $m < 0$	$0.233606481142978194 - 1.813799364234217851\thinspace i$
$(0.5,\ -1.25,\ 2)$	$n < 1$, negative $\varphi$, complex result	$-1.532738349029281103 + 1.460480317267677154\thinspace i$
$(4,\ 0.625,\ 2)$	$1 < m < n$, pole crossed before the branch point	$0.555008279584908678 - 1.282549830161864096\thinspace i$
$(4,\ 1,\ 2)$	$1 < m < n$, past the branch point	$-0.221709571883221198 - 0.751955515548111328\thinspace i$
$(1.5,\ 1.25,\ 3)$	$m > n$	$-0.780810492104561147 - 1.131673674679173545\thinspace i$

The last tests check the sign convention without relying on reference values:

• $\Pi(n - i\varepsilon;\thinspace\varphi\thinspace|\thinspace m)$, computed through the complex path, must equal the new real-$n$ result.
• $\Pi(n + i\varepsilon;\thinspace\varphi\thinspace|\thinspace m)$ must equal its complex conjugate.

They use $\varepsilon = 10^{-14}$ in float and $10^{-30}$ in bigfloat.

With the patch, all 288 tests in the file pass. Without the source changes, 11 of the new tests fail. The four that still pass cover regions where the rule adds nothing: the no-pole case, its imaginary-part check, past the branch point, and $m > n$.

These results come from an SBCL build. Tolerances may need adjusting for other Lisps.

A wider comparison of 336 cases also matched mpmath:

• $n \in \lbrace 1.1,\ 1.5,\ 4\rbrace$
• $m \in \lbrace -2,\ 0,\ 0.5,\ 0.99,\ 1.2,\ 2,\ 3,\ 6\rbrace$
• 14 values of $\varphi$, including negative values and several cycles of $\pi$

The worst relative errors were:

	$m \le 1$	$m > 1$
Bigfloat at fpprec: 30	$6 \times 10^{-27}$	$5 \times 10^{-26}$
Double float	$8 \times 10^{-16}$	$4 \times 10^{-15}$

9. Known limitations

• $m = 1$. bf-elliptic-pi always evaluates the complete integral, which is infinite at $m = 1$, so even $\varphi < \pi/2$ fails with 0.0/0.0. This also happens on master.
• Complex arguments. The correction applies only to real $n$, $\varphi$ and $m$. Complex inputs go through bf-rj1 unchanged. Continuity with the real-axis value from below is tested, but that is the only check of the complex path here.

10. References

• B. C. Carlson, "Numerical computation of real or complex elliptic integrals", Numerical Algorithms 10 (1995) 13–26. arXiv:math/9409227. Eq. (33) gives the Cauchy principal value of $R_J$ for $p < 0$.
• NIST Digital Library of Mathematical Functions, §19.20(iii) for $R_J$ with $p < 0$ (19.20.14), and §19.25 for $\Pi$ in terms of Carlson integrals.
• N. I. Muskhelishvili, Singular Integral Equations, and F. D. Gakhov, Boundary Value Problems. Both are standard treatments of the Sokhotski–Plemelj formulas. Section numbers have not been verified.
• F. Johansson et al., mpmath, ellippi, used for the reference values.