Fix rem2pi for NaN

[mgautam98]

Fix for #32888

This will also fix the mod2pi.

There was a discussion that Domain Error should be thrown or NaN should be returned. To be consistent with rem I returned NaN

[mgautam98]

@simeonschaub can you please review it?

[kimikage]

I guess this can be optimized:

Suggested change

	isnan(x) && return NaN
	abs(x) < pi && return x
	abs(x) >= pi || return x

[mgautam98]

@kimikage Yes. but won't it make it hard to understand and debug later?

[kimikage]

That's right. However, for such low-level functions, I prefer the speed. Is it not enough to just add comments into the source?

[StefanKarpinski]

Worth checking if LLVM figures that out and eliminates the isnan check.

[simeonschaub]

It seems like LLVM with optimize=true can't eliminate the fcmp ord double %0, 0.000000e+00 instruction and extra branch from isnan here. So I think it makes sense to use >= here to replace the isnan check and just add an explanatory comment about the NaN case.

[oscardssmith]

what's the performance impact of this pr?

[giordano]

About 5% according to

julia> using BenchmarkTools

julia> @benchmark x .= mod2pi.(x) setup=(x=randn(10_000))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  64.060 μs … 177.371 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     66.712 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   68.426 μs ±   6.567 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅█▇▅▇▇▄▄▄▂        ▁▁▁▁▂▁                 ▁                  ▂
  ▇██████████▆▅▅▁▅▆▇█████████▇▄▅▄▄▄▅▃▅▅██▇███▇▄▄▃▃▁▄▃▄▄▄▅▅▇▇██ █
  64.1 μs       Histogram: log(frequency) by time        99 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @eval Base.Math function rem2pi(x::Float64, ::RoundingMode{:Down})
           isnan(x) && return NaN

           if x < pi4o2_h
               if x >= 0
                   return x
               elseif x > -pi4o2_h
                   return add22condh(x,0.0,pi4o2_h,pi4o2_l)
               end
           end

           n,y = rem_pio2_kernel(x)

           if iseven(n)
               if n & 2 == 2 # n % 4 == 2: add pi
                   return add22condh(y.hi,y.lo,pi2o2_h,pi2o2_l)
               else          # n % 4 == 0: add 0 or 2pi
                   if y.hi > 0
                       return y.hi+y.lo
                   else      # negative: add 2pi
                       return add22condh(y.hi,y.lo,pi4o2_h,pi4o2_l)
                   end
               end
           else
               if n & 2 == 2 # n % 4 == 3: add 3pi/2
                   return add22condh(y.hi,y.lo,pi3o2_h,pi3o2_l)
               else          # n % 4 == 1: add pi/2
                   return add22condh(y.hi,y.lo,pi1o2_h,pi1o2_l)
               end
           end
       end
rem2pi (generic function with 9 methods)

julia> @benchmark x .= mod2pi.(x) setup=(x=randn(10_000))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  67.314 μs … 171.758 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     69.521 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   72.705 μs ±   8.752 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▂██▅▇▅▄▅▂    ▁ ▂▁▁▁▁                                         ▂
  ██████████▇▇████████▇▇▇▇█▇▇▇████▇▆▆▅▇▅▆▆▆▆▇▇▇▆▅▆▅▆▆▆▆▄▅▄▅▄▅▇ █
  67.3 μs       Histogram: log(frequency) by time       113 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

Unless you have smarter ways to deal with non-finite values, this seems a decent trade-off.

@JeffreySarnoff just reported on Slack that also Inf is wrong

julia> mod2pi.((-Inf, Inf, NaN))
(1.1963318081441687, 5.086853499035418, 1.3470949413735402)

[oscardssmith]

removing back-port labels since #46163 is a more complete fix.

[fredrikekre]

Pretty sure you need both to make the other patch apply cleanly (and avoid manual work).