Weird behavior of ArbNumerics returning nan

roi.holtzman · April 5, 2023, 2:35pm

I want to use elliptic integrals with complex arguments, and for that, I am using ArbNumerics.jl (as suggested in this post).

The issue is that I get nan for a reason I do not understand.
This code works great:

ArbNumerics.elliptic_e(ArbNumerics.ArbComplex(1.57079632679489661923132169164 - 3.783672704329450982374478608948im), ArbNumerics.ArbComplex(1/2))
# 15.57243365199836315804428747305 - 0.5034307962536973664952181708151im

I chose this weird number since it is the value of arcsin(22).
To get this, I defined the following function

function asin_mathematica(z)
    - 1im * log(1im * z + sqrt(1-z^2))
end

Which is the definition of arcsine by Mathematica (see Eq. (1) here)

So arcsin(22) is this value:

asin_mathematica(ArbNumerics.ArbComplex(22.0))
#  1.57079632679489661923132169164 - 3.783672704329450982374478608948im

Now comes the problem. When I run

ArbNumerics.elliptic_e(asin_mathematica(ArbNumerics.ArbComplex(22.0)), ArbNumerics.ArbComplex(1/2))
# nan + nanim

I get this nan.

I want to add some more peculiar information on this issue. If you take a value of the argument of the arcsine to be smaller than 1 then everything works fine. For example:

ArbNumerics.elliptic_e(
    asin_mathematica(ArbNumerics.ArbComplex(0.8)), 
    ArbNumerics.ArbComplex(1/2))
# 0.8681394584167536706972243469339 + 0e-38im

ArbNumerics.elliptic_e(
    asin_mathematica(ArbNumerics.ArbComplex(1.1)), 
    ArbNumerics.ArbComplex(1/2))
# nan + nanim

How can I fix this?

jd-foster · April 6, 2023, 1:44pm

This change of argument from below 1 to above 1 appears to cause the result to cross the “standard strip”, where the real component lies between -\pi/2 and \pi/2: see acb_elliptic.h – elliptic integrals and functions of complex variables — Arb 2.23.0 documentation
Perhaps consider handling the quasi-periodicity directly?

roi.holtzman · April 6, 2023, 2:52pm

Thanks @jd-foster !

Indeed, that causes the problem!

Thanks to your link, I tried to implement an auxiliary function to avoid this issue. But I am seeing more weird behavior!

Consider these two functions. The only difference is tin the 1st line where one is taking the angle phi as Float64, and the other as ArbReal.

function elliptic_e_fix1(phi, m)
    n, phi_rem = ArbNumerics.divrem(Float64(phi), π) # this line is the only difference
    phi_rem_complex = ArbNumerics.ArbComplex(phi_rem, imag(phi))
    return 2n * ArbNumerics.elliptic_e(m) + ArbNumerics.elliptic_e(phi_rem_complex, m)
end

function elliptic_e_fix2(phi, m)
    n, phi_rem = ArbNumerics.divrem(ArbNumerics.ArbReal(phi), π)  # this line is the only difference
    phi_rem_complex = ArbNumerics.ArbComplex(phi_rem, imag(phi))
    return 2n * ArbNumerics.elliptic_e(m) + ArbNumerics.elliptic_e(phi_rem_complex, m)
end

The function elliptic_e_fix2 is better, since it does not lose accuracy performing the divrem. However, the weird thing is this

elliptic_e_fix1(
    asin_mathematica(ArbNumerics.ArbComplex(1.1)),
    ArbNumerics.ArbComplex(1/2))
# 1.350643881047675464036245811974 - 0.3025914046608560182585117085402im

elliptic_e_fix2(
    asin_mathematica(ArbNumerics.ArbComplex(1.1)),
    ArbNumerics.ArbComplex(1/2))
# nan + nanim

whereas if putting the value of asin(1.1) both give the correct result:

asin_mathematica(ArbNumerics.ArbComplex(1.1))
# 1.57079632679489661923132169164 - 0.4435682543851153829493319664516im

elliptic_e_fix1(
    ArbNumerics.ArbComplex(1.57079632679489661923132169164 - 0.4435682543851153829493319664516im),
    ArbNumerics.ArbComplex(1/2))
# 1.350643881047675464036245811974 - 0.3025914046608560158890098484932im 

elliptic_e_fix2(
    ArbNumerics.ArbComplex(1.57079632679489661923132169164 - 0.4435682543851153829493319664516im),
    ArbNumerics.ArbComplex(1/2))
# 1.350643881047675464036245811974 - 0.3025914046608560158890098484932im

So it seems that something in the type ArbComplex in the remainder operation is causing the issue. But I do not know what it is. Any thoughts?

jd-foster · April 7, 2023, 3:52am

Consider also:

function elliptic_e_fix3(phi, m; bits=24)
           n, phi_rem = ArbNumerics.divrem(ArbFloat(phi, bits=bits), π) # this line is the only difference
           phi_rem_complex = ArbNumerics.ArbComplex(phi_rem, imag(phi))
           return 2n * ArbNumerics.elliptic_e(m) + ArbNumerics.elliptic_e(phi_rem_complex, m)
end

Then:

julia> for b in [24 + 8*i for i=0:20]
       print("Bits=$(b) : ")
       elliptic_e_fix3(
           asin_mathematica(ArbNumerics.ArbComplex(1.1)),
           ArbNumerics.ArbComplex(1/2); bits=b) |> println
       end
Bits=24 : 1.35064388104767211476735958579 - 0.302591404660856018258511708546im
Bits=32 : 1.35064388104767549892446337683 - 0.3025914046608560182585117085402im
Bits=40 : 1.350643881047675502467797973261 - 0.3025914046608560182585117085402im
Bits=48 : 1.350643881047675502519968464615 - 0.3025914046608560182585117085402im
Bits=56 : 1.350643881047675502520174335096 - 0.3025914046608560182585117085402im
Bits=64 : 1.350643881047675502520174733125 - 0.3025914046608560182585117085402im
Bits=72 : 1.350643881047675502520174735331 - 0.3025914046608560182585117085402im
Bits=80 : 1.350643881047675502520174735339 - 0.3025914046608560182585117085402im
Bits=88 : 1.350643881047675502520174735339 - 0.3025914046608560182585117085402im
Bits=96 : 1.350643881047675502520174735339 - 0.3025914046608560182585117085402im
Bits=104 : 1.350643881047675502520174735339 - 0e+0im
Bits=112 : nan + nanim
Bits=120 : nan + nanim
Bits=128 : nan + nanim
Bits=136 : nan + nanim
Bits=144 : nan + nanim
Bits=152 : nan + nanim
Bits=160 : nan + nanim
Bits=168 : nan + nanim
Bits=176 : nan + nanim
Bits=184 : nan + nanim

roi.holtzman · April 7, 2023, 4:06am

Interesting stuff!

I am still not sure how to solve the issue. Is it a bug?

If you have an idea for a work around I would love to hear!

jd-foster · April 7, 2023, 4:24am

Are you able to do calculations at a certain precision that avoid the nan, say, ArbFloat(x, bits=64)? Of course, this might be hard to determine in advance.
I’m not knowledgeable about the package to tell you if there is bug or why this occurs.

joeldahne · April 7, 2023, 7:13am

I believe the issue is that the function doesn’t support evaluation for phi with real part pi / 2 and non-zero imaginary part. I give the examples using https://github.com/kalmarek/Arblib.jl/ since that is what I’m more familiar with, it still uses Arb in the background.

julia> using Arblib

julia> Arblib.elliptic_e_inc!(Acb(), Acb(Arb(π) / 2, Arb(0)), Acb(1 // 2), 0)
[1.350643881047675502520174735338725841349522366924354545323253708857877890836 +/- 9.30e-76]

julia> Arblib.elliptic_e_inc!(Acb(), Acb(Arb(π) / 2, Arb(1)), Acb(1 // 2), 0)
nan + nanim

julia> Arblib.elliptic_e_inc!(Acb(), Acb(Arb(π) / 2, Arb(2.5)), Acb(1 // 2), 0)
nan + nanim

Note that it works fine if you use an approximation of pi

julia> Arblib.elliptic_e_inc!(Acb(), Acb(Float64(π) / 2, Arb(2.5)), Acb(1 // 2), 0)
[4.39449028010113575435763353672460317371090719423397674015371470868274922829 +/- 2.75e-75] + [0.50343079625369667427124870931999960962915777480350895883793450986892139466 +/- 6.29e-75]im

It is only an issue if the input overlaps pi / 2 (or is extremely close, depending on the precision). Note that this is an inherent property of the function, it is in general not continuous at pi / 2 and there is hence not always a well defined value for an input overlapping pi / 2.

Also, do you need to compute it in higher precision? Otherwise you can use the Float64 wrappers in Arb, they can be used from Arblib as

Arblib.fpwrap_elliptic_e_inc(0.5 + 0.6im, 1 / 2 + 0im, 0)

Note that it is still a bad idea to evaluate it with a real part near pi / 2, it will give different results depending on if you are slightly to the left or the right of it.

julia> Arblib.fpwrap_elliptic_e_inc(π / 2 + 3im, 1 / 2 + 0im, 0)
7.154175278771655 + 0.5034307962536968im

julia> Arblib.fpwrap_elliptic_e_inc(π / 2 + eps() + 3im, 1 / 2 + 0im, 0)
-4.452887516676304 + 0.5034307962536976im

roi.holtzman · April 7, 2023, 4:49pm

@joeldahne thanks for the very nice answer, and for introducing Arblib.jl.

That issue you pointed out with the instability of elliptic_e around real(π/2) is the key to what helped to find a solution!

y = \arcsin(x) cannot have Re(y) > π/2. But it seems that my function

function asin_mathematica(z)
    - 1im * log(1im * z + sqrt(1-z^2))
end

was returning Re(y) which are slightly larger than π/2. That is what caused the elliptic_e to return nan.

To fix it, I defined a function that makes sure the real part is not greater than π/2

function asin_mathematica(z)
    asinz = - 1im * log(1im * z + sqrt(1-z^2))
    # this makes sure that the real part is not greater than π/2
    re_asinz = min(π/2, real(asinz))
    return re_asinz + 1im*imag(asinz)
end

Using this function I don’t get nan anymore

ArbNumerics.elliptic_e(asin_mathematica(ArbNumerics.ArbComplex(1.1)), ArbNumerics.ArbComplex(1/2))
# 1.350643881047675464036245811974 - 0.3025914046608560182585117085402im

joeldahne · April 7, 2023, 6:33pm

Good that you found a working solution! It can be worth to note that this only works since Float64(pi) happens to be less than \pi, if you for example work at a different precision this might no longer be the case.

roi.holtzman · April 7, 2023, 6:44pm

That’s a good point!

roi.holtzman · April 7, 2023, 11:34pm

Do you have any idea on how to make sure that in any precision it works fine?

joeldahne · April 8, 2023, 12:30pm

One simple way would be to use min(prevfloat(oftype(asinz, π) / 2), asinz). Maybe a bit ugly, but should work in practice. If you want to evaluate at exactly \pi / 2 you can use the pi flag for Arblib.elliptic_e_inc!.