Perimeter of Ellipse

[Albert Chan]

We can get E(e),K(e) by recursively going for E(h),K(h), like this.
Note: code use parameter m, not modulus k, m = k^2
Note: code quit when m is small enough for taylor linear term estimate.

E(m=ε) ≈ pi/2 - (pi/8)*ε ,     K(m=ε) ≈ pi/2 + (pi/8)*ε

Code:

from cmath import sqrt, pi
    
def EK(m, verbal=False, c=pi/2):
    e = abs(m)
    if 1+e*e == 1: m *= 0.25*c; return c-m, c+m
    b = sqrt(1-m)
    if not b: return 1, 99  # K(1) assumed 99 (should be inf)
    h = m / ((b+2)*b+1)     # = (1-b)/(1+b)
    e, k = EK(h*h, verbal)
    hk = h*k
    e = e-k + hk + b*e      # = 2/(1+h)*e - (1-h)*k
    k = hk + k              # = (1+h)*k
    if verbal: print 'm = %s\n\tE(m) = %s\n\tK(m) = %s' % (m, e, k)
    return e, k

>>> e, k = EK(0.96, verbal=True)
m = (2.89332317725e-05+0j)
      E(m) = (1.57078496468+0j)
      K(m) = (1.57080768903+0j)
m = (0.0212862362522+0j)
      E(m) = (1.56240357945+0j)
      K(m) = (1.57925700384+0j)
m = (0.444444444444+0j)
      E(m) = (1.3781039379+0j)
      K(m) = (1.80966749549+0j)
m = 0.96
      E(m) = (1.05050222698+0j)
      K(m) = (3.01611249248+0j)
>>> 4 * 50 * abs(e)                       # ellipse_perimeter(10,50), error = -3 ULP
210.10044539689011

>>> e, k = EK(2, verbal=True)        # see https://www.hpmuseum.org/forum/thread-15...#pid132745
m = (5.57959210499e-05-0j)
      E(m) = (1.57077441556+0j)
      K(m) = (1.57081823849+0j)
m = (0.0294372515229-0j)
      E(m) = (1.55917174457+0j)
      K(m) = (1.58255172722+0j)
m = (-1-0j)
      E(m) = (1.91009889451+0j)
      K(m) = (1.31102877715+0j)
m = 2
      E(m) = (0.599070117368+0.599070117368j)
      K(m) = (1.31102877715-1.31102877715j)