Special Functions (Elliptic.jl) and arbitrary precision

I do not understand how to work in arbitrary precision in Julia.

I encountered the problem when using a function from the Elliptic.jl package

         K(1-10^-23)

returns
Inf

while for example Mathematica accurately computes a value of approx- 30.168, using 200 digits of working precision.

If my understanding is right, the BigFloat type shoud allow arbitrary precision computations in Julia.
I tried

setprecision(500)
a = BigFloat(1-10^-23)
K(a)

but still obtain Inf as a result. As a matter of fact

 BigFloat(1-10^-23) == 1

returns True.

Could somebody please explain what is that I am doing wrong? Thanks a lot

You identify the source of the problem correctly: precision is lost before converting to BigFloat.

You should convert something to BigFloat before performing that calculation, eg make

a = 1 - BigFloat(10)^-23

or similar.

Thanks for your reply, but I did try that one too

      a  = 1 -BigFloat(10)^-23

and while I now noticed that indeed, a =! 1

still

        K(a)

diverges.
At least now I can define the argument correctly. Could it be something to do with the function itself?
Thanks

This is probably a problem with the package? I see here that K converts everything to Float64.

3 Likes

As an alternative, the DoubleFloats package exports the ellipk function, but it seems to give a different answer than Mathematica:

julia> using DoubleFloats; ellipk(1-Double64(10)^(-23))
27.866022930551416
2 Likes

Another alternative: https://github.com/JuliaMath/SpecialFunctions.jl/blob/37e5999c5fe0f778f2de2e41b6beeff339a21cb8/src/elliptic.jl

This precision loss appears to be a known issue with the package, but it has not been addressed in two-and-a-half years.

1 Like
julia> using ArbNumerics
julia> setextrabits(64); setprecision(ArbFloat, bits=256);
julia> x = 1 - ArbFloat(10)^(-23)
0.99999999999999999999999
julia> elliptic_k(x)
27.866022930551415985041433136843867902103629991741788896439319076571988943551
3 Likes

Nice! This should be the one then.

Also, @StarecZosima could you kindly check whether your Mathematica is doing something strange? Mine gives the same answer of ~27.866.

Thanks to all very much, I got it working now, good to learn about these additional packages.
@EHBaozi, I inadvertendly typed 10^-25 and not 10^-23 when checking. Indeed in Mathemica

 EllipticK [1-10^-23] = 27.866022930551415985

so matching what you got and Julia’s ArbNumerics.

Thanks once more