# Weird behavior of ArbNumerics returning nan

**URL:** <https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114>\
**Category:** Numerics\
**Tags:** question, numerics, complex-numbers\
**Created:** [April 5, 2023, 2:35pm UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114 "2023-04-05T14:35:55Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [April 5, 2023, 2:35pm UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/1 "2023-04-05T14:35:55Z")

</div>

I want to use elliptic integrals with complex arguments, and for that, I am using `ArbNumerics.jl` (as suggested in [this post](https://discourse.julialang.org/t/elliptic-functions-of-complex-argument/51876/4)).

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

```julia
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

```julia
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](https://mathworld.wolfram.com/InverseSine.html))

So `arcsin(22)` is this value:

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

```

Now comes the problem. When I run

```julia
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:

```julia
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?

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [April 6, 2023, 1:44pm UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/2 "2023-04-06T13:44:28Z")

</div>

> [@roi.holtzman](#):
>
> `asin_mathematica(ArbNumerics.ArbComplex(1.1))`

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](https://arblib.org/acb_elliptic.html?highlight=elliptic_e#c.acb_elliptic_e_inc)  
Perhaps consider handling the quasi-periodicity directly?

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [April 6, 2023, 2:52pm UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/3 "2023-04-06T14:52:59Z")

</div>

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`.

```julia
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**

```julia
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:

```julia
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?

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [April 7, 2023, 3:52am UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/4 "2023-04-07T03:52:07Z")

</div>

Consider also:

```julia
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
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

```

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [April 7, 2023, 4:06am UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/5 "2023-04-07T04:06:19Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [April 7, 2023, 4:24am UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/6 "2023-04-07T04:24:36Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![joeldahne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joeldahne/32/45861_2.png) [@joeldahne](https://discourse.julialang.org/u/joeldahne)\
**Post date:** [April 7, 2023, 7:13am UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/7 "2023-04-07T07:13:37Z")

</div>

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
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
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

```julia
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
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

```

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [April 7, 2023, 4:49pm UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/8 "2023-04-07T16:49:10Z")

</div>

@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

```julia
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

```julia
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

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

```

---

<div class="post-metadata">

**Author:** ![joeldahne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joeldahne/32/45861_2.png) [@joeldahne](https://discourse.julialang.org/u/joeldahne)\
**Post date:** [April 7, 2023, 6:33pm UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/9 "2023-04-07T18:33:33Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [April 7, 2023, 6:44pm UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/10 "2023-04-07T18:44:00Z")

</div>

That’s a good point!

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [April 7, 2023, 11:34pm UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/11 "2023-04-07T23:34:11Z")

</div>

> [@joeldahne](#):
>
> 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π\pi, if you for example work at a different precision this might no longer be the case.

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

---

<div class="post-metadata">

**Author:** ![joeldahne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joeldahne/32/45861_2.png) [@joeldahne](https://discourse.julialang.org/u/joeldahne)\
**Post date:** [April 8, 2023, 12:30pm UTC](https://discourse.julialang.org/t/weird-behavior-of-arbnumerics-returning-nan/97114/12 "2023-04-08T12:30:41Z")

</div>

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!`.