Numerical computation of the complex elliptic integral $E(k)$ for medium $|k|$

Question

I have implemented Carlson's algorithm for $E(k)$ from Numerical computation of real or complex elliptic integrals (available from ArXiv eprint, see also DLMF). It is essentially his formula (46) which simultaneously computes $K(k), E(k)$ for complex $k$ using the arithmetic-geometric sequence $a_m, b_m$ to get

$$K(k) = \frac{\pi}{a_n+b_n}, \quad E(k) = \left( \left(\frac{a_0+b_0}{2} \right)^2 -\sum\limits_{m=1}^{n} 2^{m-2}(a_m - b_m)^{2}\right) K(k)$$

Unfortunately with IEEE 64 or 80 bit floating point arithmetic the expression for $E(k)$ is accurate only up to $|k| \approx 1000,\,$ for larger arguments severe rounding/cancellation problems occur (e.g. for $k=10^3-10i$ the relative error in the real part is $17\epsilon,\;$for $k=10^4-10i$ it is $150\epsilon,\;$ and for $k=10^6-10i$ it is $24000\epsilon$). For very large $|k|$ I could use the logarithmic form, but is seems relative useless if there is no reliable algorithm for medium $|k|$.

Is there a known stable algorithm for this range, ie. which can compute $E(k)$ for $|k| > 100$ without using multi-precision?

Answer 1

After some searching I found two formulas to solve the accuracy problem. The two modifications are

In range $|k| \ge 40$ the (asymptotic) series from http://functions.wolfram.com/08.01.06.0010.02 is used (it is valid for $|k|>1,\,$ but convergence is suboptimal for small $k$): \begin{multline*} E(z)=\sqrt{-z} +\frac{\ln(-z)}{\sqrt{-z}}\sum_{n=0}^\infty \frac{\left(-\frac{1}{2}\right)_{n+1}^2 z^{-n}}{n!(n+1)!}\\ +\frac{1}{\sqrt{-z}}\sum_{n=0}^\infty \frac{\left(-\frac{1}{2}\right)_{n+1}^2 z^{-n}}{n!(n+1)!} \Big(2\psi(n+1)-2\psi(n+\tfrac{1}{2})+\frac{1}{n+1}\Big) \end{multline*} Note Wolfram's different definition of $E$ so we actually have $z = k^2$. To avoid overflows in computing $-k^2$ for very large $k,$ $\sqrt{-z}\,$ is computed in-line using three cases and $\ln(-z)\,$ is based on $\ln(k).$ Up to six terms are used in each sum, but only one for $|k| > \epsilon^{-1}$.

For certain values on the branch cut, i.e. $k \in \mathbb{R}, 4 < |k| \le 40\,$ the Legendre relation http://dlmf.nist.gov/19.7.E3

$$E\left(1/k\right)=(1/k)\left(E\left(k\right)\pm iE\left(k^{\prime}\right)-{k^{\prime}}^{2} K\left(k\right)\mp ik^{2} K\left(k^{\prime}\right)\right)$$ is applied, where special care is taken to avoid cancellation for the sums in the real and imaginary parts.

After these two changes the relative errors are bounded by about $5 \epsilon$ for IEEE double.