Ensure thread-safety for trigonometric functions

[OlivierHnt]

Closes #620

This PR makes sure we no longer rely on Julia default code to compare with π due to JuliaLang/julia#52862.
This should ensure that trigonometric functions are thread-safe.

The PR also reworks the inner function _quadrant.

Now I consistently get the same precision before and after calling sin on multiple threads (see also #612):

julia> x = interval(BigFloat, 1)
[1.0, 1.0]₂₅₆_com

julia> y = sin(x)
[0.84147, 0.841472]₂₅₆_com

julia> Threads.@threads for _ in 1:1000
           z = sin(x)
           @assert isequal_interval(y, z)
       end

julia> precision(BigFloat)
256

The implementation is based on @Joel-Dahne, @lbenet and @dpsanders suggestions.

Hopefully everything is in order! 🙂

[codecov-commenter]

Codecov Report

Attention: 8 lines in your changes are missing coverage. Please review.

Comparison is base (caebf92) 83.64% compared to head (4c9f38d) 83.45%.
Report is 2 commits behind head on master.

Files	Patch %	Lines
src/intervals/construction.jl	0.00%	8 Missing ⚠️

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #621      +/-   ##
==========================================
- Coverage   83.64%   83.45%   -0.19%     
==========================================
  Files          25       25              
  Lines        2140     2152      +12     
==========================================
+ Hits         1790     1796       +6     
- Misses        350      356       +6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[Joel-Dahne]

I don't think rem2pi can be assumed to give correct, or even faithful rounding.

For BigFloat it is defined by

rem2pi(x::BigFloat, r::RoundingMode) = rem(x, 2*BigFloat(pi), r)

which I don't think would guarantee the correct rounding since $2\pi$ is not computed to any higher precision and also not rounded in any specific direction. There is even a TODO-comment

# TODO: use a higher-precision BigFloat(pi) here?

So the one who implemented most likely didn't prove that it would be correct. I'm not sure how easy it would be to find a counter example though.

For Float64 it has a quite complicated implementation and it probably gives very good results in general. But it was also easy to find a counterexample to the correctness. If we define

g = (rng, N) -> begin
    for _ in 1:N
        x = reinterpret(Float64, rand(rng, UInt))
        rem2pi(x, RoundDown) == Float64(rem2pi(big(x), RoundDown)) || error("x = $x")
    end
end

then it was enough to run

g(MersenneTwister(42), 1000)

to get the counterexample x = 8.197839233774979e66. For this x we have that the correctly rounded result should be 1.9664659792957133, but rem2pi(x) gives us 1.9664659792957135 (so not even faithfully rounded).

I don't think this should be regarded as a bug in Julia though, most numerical methods don't give correct rounding.

[Joel-Dahne]

This issue was of course already present in the previous implementation. It is maybe reasonable to merge this PR, which fixes some issues, and leave this issue for another day.

[Joel-Dahne]

Of course it is probably extremely rare that any rounding errors in rem2pi actually affects the returned quadrant. It wouldn't surprise me if it actually gives correct results for all Float64 values.

[OlivierHnt]

Thx for peeking at the code of rem2pi. Yes perhaps this can be addressed in a separate PR. Then perhaps I should go back to using rem2pi(..., RoundNearest) since the version in this PR is flawed too, and less performant.

[Joel-Dahne]

To get a correct implementation I believe the easiest, but certainly not the fastest, is to more or less compute the remainder fully in interval arithmetic. Then you get a proper enclosure of it at least, if that enclosure is not enough to determine the quadrant you can fall back to some more pessimistic version. This would guarantee faithful, but not necessarily exactly correct determination of the quadrant.

For Float64 there is an implementation in CRlibm here. This is likely very fast and would also guarantee a correct result. I'm not sure if it is part of the public interface though...

[OlivierHnt]

Ah great thx for the reference.

I'll merge this PR (which now only focuses on thread-safety) and open an issue for _quadrant not being reliable.

[dpsanders]

We could ask the CRlibm people if we can transcribe their code into Julia. (I think the license is problematic for compatibility with our MIT license.)

I believe this is a standard algorithm, so alternatively somebody could implement it based on the relevant paper.