Incomplete elliptic integral of the first kind errors on specific input

[kagalenko-m-b]

julia> phi=big"1.5707963267948966"
1.570796326794896600000000000000000000000000000000000000000000000000000000000003

julia> ellipticF(phi, 0.81)
ERROR: StackOverflowError:
Stacktrace:
 [1] /(x::BigFloat, y::BigFloat)
   @ Base.MPFR ./mpfr.jl:431
 [2] /(x::BigFloat, y::Irrational{:π})
   @ Base ./promotion.jl:391
 [3] ellipticF(phi::BigFloat, m::Float64)
   @ EllipticFunctions ~/.julia/packages/EllipticFunctions/PFtTG/src/EllipticFunctions.jl:762
 [4] ellipticF(phi::BigFloat, m::Float64) (repeats 20931 times)
   @ EllipticFunctions ~/.julia/packages/EllipticFunctions/PFtTG/src/EllipticFunctions.jl:765
 [5] top-level scope
   @ REPL[48]:1
julia> versioninfo()
Julia Version 1.8.5
Commit 17cfb8e65ea (2023-01-08 06:45 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × Intel(R) Core(TM) i3 CPU         540  @ 3.07GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, westmere)
  Threads: 4 on 4 virtual cores
Environment:
  JULIA_NUM_THREADS = 4

If, instead, the second argument is big"0.81", the computation hangs while CPU goes to 100%. If the last digit of phi=big".." is deleted, the computation returns normally.

[simonp0420]

@stla can comment further, but it is known that most of the functions in this package do not (yet) reliably support extended precision. See Issue #11 for status of extended precision support. I've corrected the checklist there.

[stla]

The epsilon in CarlsonRF is too small. Maybe removing the ^2 is fine.

[stla]

This is not due to epsilon. There's an infinite loop because of this strange thing:

> phi > pi/2
true
> phi - pi > -pi/2
false

[stla]

Ok I changed the code and it works. Now the code deals with the comparison of phi / pi with +/- 0.5.

[simonp0420]

phi - pi > -pi/2 You probably already realized this, but because pi is of type Irrational, it has indefinitely large precision, depending on how it's used. So in the expression phi - pi it will be evaluated accurately to the same precision as phi (256 bits or about 77 decimal digits). But in the right-hand side -pi/2 it will only be evaluated to Float64 precision (53 bits or about 15 decimal digits). Hence the error. There are lots of this type of thing in the code base that will have to be fixed to support extended precision.

[stla]

Ok. I already fixed in the elliptic integrals.

[kagalenko-m-b]

I have not looked at the code, but from the error message it seems to me that the iterative algorithm to compute the integral is implemented recursively, and there's no hard limit on the depth of recursion. Perhaps, rewriting it as a loop with conditional termination and fixed upper number of iterations would enhance the robustness of the software? It is one thing to not provide all correct significant digits, but failure to return any answer at all on certain inputs seems like a serious bug.

[stla]

It is fixed now.

[kagalenko-m-b]

Thanks ! Although there might be another multiprecision issue coming up :)