(28/48/50) Complete Elliptic Integrals

[Albert Chan]

(03-04-2024 06:39 PM)Albert Chan Wrote:  Lua code had been updated, instead of starting s = a², it start at 0 (see post #7)
The reason is to produce as accurate E(m) - K(m) as possible.

I just realized I never posted lua update. Here it is

Code:

function E.agm(a, b)
    local s, t = 0, 1
    repeat
        local b0, g = b, (a-b)/2
        a, b = b+g, sqrt(a*b)
        t = t + t
        s = s - t*g*g
    until b == b0 or not finite(b)
    return b, s
end

Code:

function E.K(m,c)               -- K(m), E(m)-K(m), E(m)
    if not c or m+c~=1 then c=1-m end   -- c only a hint
    if c==0 then return inf, -inf, 1 end
    local x, y = E.agm(1, sqrt(c))
    x = pi/(4*x)
    return x+x, x*(y-m), x+x*(y+c)
end

lua> eps = 1e-8
lua> E.K(1-eps, eps)
10.596634757087662      -9.596634706604487      1.000000050483174

Compare this to asymptote result for K(1-ε), E(1-ε)

lua> k0 = log(16/eps)/2
lua> k0 + eps/4*(k0-1)      -- ≈ K(1-eps)
10.59663475708766
lua> 1 + eps/4*(2*k0-1)    -- ≈ E(1-eps)
1.0000000504831736