Skip to content

Commit

Permalink
Implement rem(float, float, RoundNearest) in Julia
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
ararslan committed Sep 28, 2021
1 parent f5d944c commit bc13ba2
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 5 deletions.
41 changes: 36 additions & 5 deletions base/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -880,11 +880,42 @@ 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)
sx = sign(x)
if p <= floatmax(T) / 2 # Ensure p+p won't overflow
x = rem(x, p + p) # now x < 2p
end
x = abs(x)
p = abs(p)
if p < 2 * floatmin(T) # Check whether dividing p by 2 will underflow
if x + x > p
x -= p
if x + x >= p
x -= p
end
end
else
p_half = T(0.5) * p
if x > p_half
x -= p
if x >= p_half
x -= p
end
end
end
return flipsign(x, sx)
end


"""
Expand Down
12 changes: 12 additions & 0 deletions test/numbers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2523,6 +2523,18 @@ 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))
if !(T == BigFloat && mode == RoundUp) # FIXME: Erroneously returns -Inf
@test rem(T(4), floatmin(T) * 2, mode) == 0
end
end
@test rem(nextfloat(typemin(T)), T(2), RoundToZero) == -0.0
@test rem(nextfloat(typemin(T)), T(2), RoundNearest) == -0.0
@test rem(nextfloat(typemin(T)), T(2), RoundDown) == 0.0
@test rem(nextfloat(typemin(T)), T(2), RoundUp) == 0.0
end

@testset "rem for $T RoundNearest" for T in (Int8, Int16, Int32, Int64, Int128)
Expand Down

0 comments on commit bc13ba2

Please sign in to comment.