Implement rem(float, float, RoundNearest) in Julia (JuliaLang#42380)

One of the few remaining uses of openlibm is implementing `rem` with
`RoundNearest` for `Float32` and `Float64`. This commit translates the
msun libm implementations of `__ieee754_remainder` and
`___ieee754_remainderf` to Julia to avoid relying on openlibm.

base/math.jl

@@ -880,11 +880,39 @@ function frexp(x::T) where T<:IEEEFloat
  return reinterpret(T, xu), k
 end
 
-rem(x::Float64, y::Float64, ::RoundingMode{:Nearest}) =
- ccall((:remainder, libm),Float64,(Float64,Float64),x,y)
-rem(x::Float32, y::Float32, ::RoundingMode{:Nearest}) =
- ccall((:remainderf, libm),Float32,(Float32,Float32),x,y)
-rem(x::Float16, y::Float16, r::RoundingMode{:Nearest}) = Float16(rem(Float32(x), Float32(y), r))
+# NOTE: This `rem` method is adapted from the msun `remainder` and `remainderf`
+# functions, which are under the following license:
+#
+# Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+#
+# Developed at SunSoft, a Sun Microsystems, Inc. business.
+# Permission to use, copy, modify, and distribute this
+# software is freely granted, provided that this notice
+# is preserved.
+function rem(x::T, p::T, ::RoundingMode{:Nearest}) where T<:IEEEFloat
+ (iszero(p) || !isfinite(x) || isnan(p)) && return T(NaN)
+ x == p && return copysign(zero(T), x)
+ oldx = x
+ x = abs(rem(x, 2p)) # 2p may overflow but that's okay
+ p = abs(p)
+ if p < 2 * floatmin(T) # Check whether dividing p by 2 will underflow
+ if 2x > p
+ x -= p
+ if 2x >= p
+ x -= p
+ end
+ end
+ else
+ p_half = p / 2
+ if x > p_half
+ x -= p
+ if x >= p_half
+ x -= p
+ end
+ end
+ end
+ return flipsign(x, oldx)
+end
 
 
 """

test/numbers.jl

@@ -2523,6 +2523,17 @@ end
  @test rem(T(-1.5), T(2), RoundNearest) == 0.5
  @test rem(T(-1.5), T(2), RoundDown) == 0.5
  @test rem(T(-1.5), T(2), RoundUp) == -1.5
+ for mode in [RoundToZero, RoundNearest, RoundDown, RoundUp]
+ @test isnan(rem(T(1), T(0), mode))
+ @test isnan(rem(T(Inf), T(2), mode))
+ @test isnan(rem(T(1), T(NaN), mode))
+ # FIXME: The broken case erroneously returns -Inf
+ @test rem(T(4), floatmin(T) * 2, mode) == 0 broken=(T == BigFloat && mode == RoundUp)
+ end
+ @test isequal(rem(nextfloat(typemin(T)), T(2), RoundToZero), -0.0)
+ @test isequal(rem(nextfloat(typemin(T)), T(2), RoundNearest), -0.0)
+ @test isequal(rem(nextfloat(typemin(T)), T(2), RoundDown), 0.0)
+ @test isequal(rem(nextfloat(typemin(T)), T(2), RoundUp), 0.0)
 end
 
 @testset "rem for $T RoundNearest" for T in (Int8, Int16, Int32, Int64, Int128)