6
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ The AGM-based method is really only useful for moduli in the traditional range $-1<k<1$; for everything outside that range, one uses the reciprocal or imaginary modulus transformations, as you've noted in your answer. $\endgroup$
    – J. M.
    Dec 11, 2016 at 9:49

1 Answer 1

6
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.