Julia: sinpi and cospi return incorrect results when using `--math-mode=fast`

[AboAmmar]

Julia 1.0:

I observed that when starting up Julia with the --math-mode=fast option, I get incorrect results for sinpi and cospi functions.

With the default math mode:

julia> a = rand(5)
5-element Array{Float64,1}:
 0.5897641978923696
 0.2459464145866952
 0.06373272357797632
 0.26244049804714087
 0.7254898186800567

julia> cospi.(a)
5-element Array{Float64,1}:
 -0.27827964958957724
  0.7160540045153001
  0.9800223981473272
  0.6789380012385138
 -0.650617399829044

julia> sinpi.(a)
5-element Array{Float64,1}:
 0.9605000971495538
 0.6980448858186719
 0.19888715174581242
 0.7341956077737403
 0.7594056880480248

while with --math-mode=fast option:

julia> cospi.(a)
5-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
 1.0

julia> sinpi.(a)
5-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0

[vtjnash]

Yes, that’s what it means to run in fast math mode. If you wanted correct answers, you should pass the ieee flag.

[StefanKarpinski]

Fast math mode should be off by no more than a few ulps; this not even close to correct:

# normal: `julia`

julia> sin(0.25pi)
0.7071067811865475

julia> sinpi(0.25)
0.7071067811865476

# fast math: `julia --math-mode=fast`

julia> sin(0.25pi)
0.7071067811865475

julia> sinpi(0.25)
0.0

[StefanKarpinski]

I suspect that this is a consequence of fast math mode performing some optimization which the implementation of sinpi and cospi relies on not happening—i.e. they are relying on correct IEEE behavior. So technically, this is what you're asking for but it seems like the better behavior here would be to replace sinpi(x) with sin(pi*x) in fast math mode instead of doing range reduction.

[tomaklutfu]

If you look at llvm ir of cospi(...::Float64) in fast mode it just returns 1.00 and it goes away if --inline=no flag is used together.

[yuyichao]

https://github.com/JuliaLang/julia/blob/6b2d211d22288ef7e66fc0f5bd21d8db952f28ba/base/special/trig.jl#L826 is 0 under fast math. That's pretty much the expected result...

[Keno]

... In which we discover why global fast math flags are a horrible idea. I guess we could have a @nonfastmath macro to say "no really, I need IEEE behavior here".

[vtjnash]

Fast math mode should be off by no more than a few ulps;

That's a cute idea, but not what the flag means.

[StefanKarpinski]

... In which we discover why global fast math flags are a horrible idea.

💯%. However, we have it so should make this less broken.

This patch makes @fastmath sinpi(x) call sin(pi*x) instead of calling the sinpi implementation. However, it doesn't fix the behavior with the global --math-mode=fast flag on. Does that flag have a different effect than wrapping everything in @fastmath?

[StefanKarpinski]

Stop closing the issue, @vtjnash.

[vtjnash]

dup https://github.com/JuliaLang/julia/issues/24568
dup https://github.com/JuliaLang/julia/issues/21375
ref https://github.com/JuliaLang/julia/issues/25028

[vtjnash]

Does that flag have a different effect than

dup https://github.com/JuliaLang/julia/issues/26828

[vtjnash]

Why can't we close duplicate issues with extensive prior discussion?

[StefanKarpinski]

Forgot to include the patch:

diff --git a/base/fastmath.jl b/base/fastmath.jl
index 0dccf8ba4b..482050da81 100644
--- a/base/fastmath.jl
+++ b/base/fastmath.jl
@@ -61,6 +61,7 @@ const fast_op =
          :cbrt => :cbrt_fast,
          :cis => :cis_fast,
          :cos => :cos_fast,
+         :cospi => :cospi_fast,
          :cosh => :cosh_fast,
          :exp10 => :exp10_fast,
          :exp2 => :exp2_fast,
@@ -75,6 +76,7 @@ const fast_op =
          :min => :min_fast,
          :minmax => :minmax_fast,
          :sin => :sin_fast,
+         :sinpi => :sinpi_fast,
          :sincos => :sincos_fast,
          :sinh => :sinh_fast,
          :sqrt => :sqrt_fast,
@@ -310,6 +312,9 @@ sincos_fast(v) = (sin_fast(v), cos_fast(v))
     min_fast(x::T, y::T) where {T<:FloatTypes} = ifelse(y > x, x, y)
     minmax_fast(x::T, y::T) where {T<:FloatTypes} = ifelse(y > x, (x,y), (y,x))

+    cospi_fast(x::Union{Real,Complex}) = cos(pi*x)
+    sinpi_fast(x::Union{Real,Complex}) = sin(pi*x)
+
     # complex numbers

     function cis_fast(x::T) where {T<:FloatTypes}

With this patch applied we have this which is amusing:

julia> sinpi(0.25)
0.0

julia> @fastmath sinpi(0.25)
0.7071067811865475

It seems like one approach to improving this situation would be automatically applying the @fastmath transformation when in --math-mode=fast since then at least any specific behaviors that @fastmath applies will happen in the global mode too, which seems better.

[StefanKarpinski]

Why can't we close duplicate issues with extensive prior discussion?

Because I reopened it for discussion and asked you to stop closing it.

[vtjnash]

Please stop re-opening it then, and do discussion on discourse

[StefanKarpinski]

First off, do not tell me how to use GitHub on the JuliaLang project—it is rude and disrespectful.

The problem here is not the general brokenness of --math-mode=fast, which I am extremely well aware of. The specific problem here is that sinpi and cospi are more broken under --math-mode=fast than necessary. I have proposed a specific fix to this problem: apply the patch I've posted here and change --math-mode=fast to apply @fastmath universally. That seems like it would fix this issue.

[StefanKarpinski]

Is there some reason we could not apply @fastmath everywhere when --math-mode=fast is requested? Is that a terrible idea for some reason I'm not seeing? It seems potentially surprising that it isn't done that way already. I understand that --math-mode=fast is asking LLVM for non-IEEE optimizations while @fastmath is doing Julia-level transformations instead but it seems like asking for the former would imply that you want the latter done implicitly for all user code.

[vtjnash]

It would hinder our inference in more cases (we shouldn't run inference on floating point numbers when fast-math might get specified), wouldn't alter the case where the user calls sinpi directly, and leads to weird scope issues (@fastmath is a syntactic, not semantic, transform—doc issue #26828—so applying it would not preserve behavior). It also can't really fix the correctness issues with that flag, just handle a different subset of cases until the compiler becomes smarter (or does something dumber)...

What does it mean to "apply a macro everywhere"?

[StefanKarpinski]

What does it mean to "apply a macro everywhere"?

I guess it would mean apply the macro to the input of eval before evaling.