Add elliptic integrals to ArbFloats.jl?

[PerezHz]

Hi Jeffrey! I'm looking for a way to evaluate the complete elliptic integral of the first kind, K(k), in Julia; there is Elliptic.jl but that package only allows double precision output. @simonbyrne and @ChrisRackauckas pointed to me that in the Arb library there is support for elliptic integrals so I wanted to ask you if this could be added to ArbFloats.jl? Many thanks in advance!

[JeffreySarnoff]

The elliptic functions that Arb implements are over the domain of Arb's version of complex variables. ArbFloats has no support for ArbComplex today. Nonetheless, the lack of support for Elliptic functions has bothered me for some time. I will take another look at it and let you know. Thanks for the inquiry, Jeffrey Sarnoff

…

On Wed, Oct 18, 2017 at 12:48 PM, Jorge Pérez ***@***.***> wrote: Hi Jeffrey! I'm looking for a way to evaluate the complete elliptic integral of the first kind, K(k), in Julia; there is Elliptic.jl but that package only allows double precision output. @simonbyrne <https://github.com/simonbyrne> and @ChrisRackauckas <https://github.com/chrisrackauckas> pointed to me that in the Arb library there is support for elliptic functions <http://fredrikj.net/blog/2017/02/arb-2-10-0-released/> so I wanted to ask you if this could be added to ArbFloats.jl? Many thanks in advance! — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#20>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABmqxnlKANF7I-5mdnFU4wTfS4vvMQT_ks5stivwgaJpZM4P-BQk> .

[JeffreySarnoff]

meanwhile, for x in 0..1

using ArbFloats

arb_precision = 750 # bits, select: 100..3200
setprecision(ArbFloat, arb_precision)

function ellK(x::Float64)
    @assert 0 <= x <= 1.0
    x == 1.0 && return ArbFloat(Inf)

    arb_one    = one(ArbFloat)
    arb_halfpi = ArbFloat(pi)/2
        
    k = sqrt( abs(ArbFloat(x) - arb_one) )
    k = agm(arb_one, k)
    k = arb_halfpi / k

    return k
end

[PerezHz]

Thanks a lot Jeffrey!! Works like a charm!!

[JeffreySarnoff]

I am happy to know that.

[JeffreySarnoff]

notes to myself this version is better when called with an ArbFloat

using ArbFloats

function ellipticK(x::T) where T
    !(zero(T) <= x <= one(T)) && throw(DomainError())
    x == 1.0 && return ArbFloat(Inf)

    arb_one = one(ArbFloat)
    arb_halfpi = ArbFloat(pi)/2
    
    k = T<:ArbFloat ? x : ArbFloat(x)

    k = arb_one - k
    k = sqrt(k)
    k = agm(arb_one, k)
    k = arb_halfpi / k

    return k
end

setprecision(ArbFloat, 256) # bits
x = ArbFloat(0.5)

k = ellipticK(x)
string("K(", x, ") = ", smartstring(K_onehalf))
"1.854074677301371918433850347195260046217598823521766905585928045056021777₋"