Skip to content

Commit

Permalink
Add native julia fmod (#47501)
Browse files Browse the repository at this point in the history 
* Add native julia rem

Co-authored-by: Alex Arslan <ararslan@comcast.net>
(cherry picked from commit cf5ae03)
  • Loading branch information
@gbaraldi @KristofferC
gbaraldi authored and KristofferC committed Dec 8, 2022
1 parent dba443d commit ce7a372
Show file tree
Hide file tree
Showing 2 changed files with 227 additions and 2 deletions.
106 changes: 104 additions & 2 deletions base/float.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,8 @@ exponent_one(::Type{Float16}) = 0x3c00
exponent_half(::Type{Float16}) = 0x3800
significand_mask(::Type{Float16}) = 0x03ff

mantissa(x::T) where {T} = reinterpret(Unsigned, x) & significand_mask(T)

for T in (Float16, Float32, Float64)
@eval significand_bits(::Type{$T}) = $(trailing_ones(significand_mask(T)))
@eval exponent_bits(::Type{$T}) = $(sizeof(T)*8 - significand_bits(T) - 1)
Expand Down Expand Up @@ -414,9 +416,109 @@ muladd(x::T, y::T, z::T) where {T<:IEEEFloat} = muladd_float(x, y, z)
# TODO: faster floating point fld?
# TODO: faster floating point mod?

rem(x::T, y::T) where {T<:IEEEFloat} = rem_float(x, y)
function unbiased_exponent(x::T) where {T<:IEEEFloat}
return (reinterpret(Unsigned, x) & exponent_mask(T)) >> significand_bits(T)
end

function explicit_mantissa_noinfnan(x::T) where {T<:IEEEFloat}
m = mantissa(x)
issubnormal(x) || (m |= significand_mask(T) + uinttype(T)(1))
return m
end

function _to_float(number::U, ep) where {U<:Unsigned}
F = floattype(U)
S = signed(U)
epint = unsafe_trunc(S,ep)
lz::signed(U) = unsafe_trunc(S, Core.Intrinsics.ctlz_int(number) - U(exponent_bits(F)))
number <<= lz
epint -= lz
bits = U(0)
if epint >= 0
bits = number & significand_mask(F)
bits |= ((epint + S(1)) << significand_bits(F)) & exponent_mask(F)
else
bits = (number >> -epint) & significand_mask(F)
end
return reinterpret(F, bits)
end

@assume_effects :terminates_locally :nothrow function rem_internal(x::T, y::T) where {T<:IEEEFloat}
xuint = reinterpret(Unsigned, x)
yuint = reinterpret(Unsigned, y)
if xuint <= yuint
if xuint < yuint
return x
end
return zero(T)
end

e_x = unbiased_exponent(x)
e_y = unbiased_exponent(y)
# Most common case where |y| is "very normal" and |x/y| < 2^EXPONENT_WIDTH
if e_y > (significand_bits(T)) && (e_x - e_y) <= (exponent_bits(T))
m_x = explicit_mantissa_noinfnan(x)
m_y = explicit_mantissa_noinfnan(y)
d = urem_int((m_x << (e_x - e_y)), m_y)
iszero(d) && return zero(T)
return _to_float(d, e_y - uinttype(T)(1))
end
# Both are subnormals
if e_x == 0 && e_y == 0
return reinterpret(T, urem_int(xuint, yuint) & significand_mask(T))
end

m_x = explicit_mantissa_noinfnan(x)
e_x -= uinttype(T)(1)
m_y = explicit_mantissa_noinfnan(y)
lz_m_y = uinttype(T)(exponent_bits(T))
if e_y > 0
e_y -= uinttype(T)(1)
else
m_y = mantissa(y)
lz_m_y = Core.Intrinsics.ctlz_int(m_y)
end

tz_m_y = Core.Intrinsics.cttz_int(m_y)
sides_zeroes_cnt = lz_m_y + tz_m_y

# n>0
exp_diff = e_x - e_y
# Shift hy right until the end or n = 0
right_shift = min(exp_diff, tz_m_y)
m_y >>= right_shift
exp_diff -= right_shift
e_y += right_shift
# Shift hx left until the end or n = 0
left_shift = min(exp_diff, uinttype(T)(exponent_bits(T)))
m_x <<= left_shift
exp_diff -= left_shift

m_x = urem_int(m_x, m_y)
iszero(m_x) && return zero(T)
iszero(exp_diff) && return _to_float(m_x, e_y)

while exp_diff > sides_zeroes_cnt
exp_diff -= sides_zeroes_cnt
m_x <<= sides_zeroes_cnt
m_x = urem_int(m_x, m_y)
end
m_x <<= exp_diff
m_x = urem_int(m_x, m_y)
return _to_float(m_x, e_y)
end

function rem(x::T, y::T) where {T<:IEEEFloat}
if isfinite(x) && !iszero(x) && isfinite(y) && !iszero(y)
return copysign(rem_internal(abs(x), abs(y)), x)
elseif isinf(x) || isnan(y) || iszero(y) # y can still be Inf
return T(NaN)
else
return x
end
end

function mod(x::T, y::T) where T<:AbstractFloat
function mod(x::T, y::T) where {T<:AbstractFloat}
r = rem(x,y)
if r == 0
copysign(r,y)
Expand Down
123 changes: 123 additions & 0 deletions test/numbers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2929,3 +2929,126 @@ end
@test false == ceil(Bool, -0.7)
end
end

@testset "modf" begin
@testset "remd" begin
denorm_min = nextfloat(0.0)
minfloat = floatmin(Float64)
maxfloat = floatmax(Float64)
values = [3.0,denorm_min,-denorm_min, minfloat,
-minfloat, maxfloat, -maxfloat]
# rem (0, y) == 0 for y != 0.
for val in values
@test isequal(rem(0.0, val), 0.0)
end
# rem (-0, y) == -0 for y != 0.
for val in values
@test isequal(rem(-0.0, val), -0.0)
end
# rem (+Inf, y) == NaN
values2 = [3.0,-1.1,0.0,-0.0,denorm_min,minfloat,
maxfloat,Inf,-Inf]
for val in values2
@test isequal(rem(Inf, val), NaN)
end
# rem (-Inf, y) == NaN
for val in values2
@test isequal(rem(-Inf, val), NaN)
end
# rem (x, +0) == NaN
values3 = values2[begin:end-2]
for val in values3
@test isequal(rem(val, 0.0), NaN)
end
# rem (x, -0) == NaN
for val in values3
@test isequal(rem(val, -0.0), NaN)
end
# rem (x, +Inf) == x for x not infinite.
@test isequal(rem(0.0, Inf), 0.0)
@test isequal(rem(-0.0, Inf), -0.0)
@test isequal(rem(denorm_min, Inf), denorm_min)
@test isequal(rem(minfloat, Inf), minfloat)
@test isequal(rem(maxfloat, Inf), maxfloat)
@test isequal(rem(3.0, Inf), 3.0)
# rem (x, -Inf) == x for x not infinite.
@test isequal(rem(0.0, -Inf), 0.0)
@test isequal(rem(-0.0, -Inf), -0.0)
@test isequal(rem(denorm_min, -Inf), denorm_min)
@test isequal(rem(minfloat, -Inf), minfloat)
@test isequal(rem(maxfloat, -Inf), maxfloat)
@test isequal(rem(3.0, -Inf), 3.0)
#NaN tests
@test isequal(rem(0.0, NaN), NaN)
@test isequal(rem(1.0, NaN), NaN)
@test isequal(rem(Inf, NaN), NaN)
@test isequal(rem(NaN, 0.0), NaN)
@test isequal(rem(NaN, 1.0), NaN)
@test isequal(rem(NaN, Inf), NaN)
@test isequal(rem(NaN, NaN), NaN)
#Sign tests
@test isequal(rem(6.5, 2.25), 2.0)
@test isequal(rem(-6.5, 2.25), -2.0)
@test isequal(rem(6.5, -2.25), 2.0)
@test isequal(rem(-6.5, -2.25), -2.0)
values4 = [maxfloat,-maxfloat,minfloat,-minfloat,
denorm_min, -denorm_min]
for val in values4
@test isequal(rem(maxfloat,val), 0.0)
end
for val in values4
@test isequal(rem(-maxfloat,val), -0.0)
end
@test isequal(rem(minfloat, maxfloat), minfloat)
@test isequal(rem(minfloat, -maxfloat), minfloat)
values5 = values4[begin+2:end]
for val in values5
@test isequal(rem(minfloat,val), 0.0)
end
@test isequal(rem(-minfloat, maxfloat), -minfloat)
@test isequal(rem(-minfloat, -maxfloat), -minfloat)
for val in values5
@test isequal(rem(-minfloat,val), -0.0)
end
values6 = values4[begin:end-2]
for val in values6
@test isequal(rem(denorm_min,val), denorm_min)
end
@test isequal(rem(denorm_min, denorm_min), 0.0)
@test isequal(rem(denorm_min, -denorm_min), 0.0)
for val in values6
@test isequal(rem(-denorm_min,val), -denorm_min)
end
@test isequal(rem(-denorm_min, denorm_min), -0.0)
@test isequal(rem(-denorm_min, -denorm_min), -0.0)
#Max value tests
values7 = [0x3p-1074,-0x3p-1074,0x3p-1073,-0x3p-1073]
for val in values7
@test isequal(rem(0x1p1023,val), 0x1p-1073)
end
@test isequal(rem(0x1p1023, 0x3p-1022), 0x1p-1021)
@test isequal(rem(0x1p1023, -0x3p-1022), 0x1p-1021)
for val in values7
@test isequal(rem(-0x1p1023,val), -0x1p-1073)
end
@test isequal(rem(-0x1p1023, 0x3p-1022), -0x1p-1021)
@test isequal(rem(-0x1p1023, -0x3p-1022), -0x1p-1021)

end

@testset "remf" begin
@test isequal(rem(Float32(0x1p127), Float32(0x3p-149)), Float32(0x1p-149))
@test isequal(rem(Float32(0x1p127), -Float32(0x3p-149)), Float32(0x1p-149))
@test isequal(rem(Float32(0x1p127), Float32(0x3p-148)), Float32(0x1p-147))
@test isequal(rem(Float32(0x1p127), -Float32(0x3p-148)), Float32(0x1p-147))
@test isequal(rem(Float32(0x1p127), Float32(0x3p-126)), Float32(0x1p-125))
@test isequal(rem(Float32(0x1p127), -Float32(0x3p-126)), Float32(0x1p-125))
@test isequal(rem(-Float32(0x1p127), Float32(0x3p-149)), -Float32(0x1p-149))
@test isequal(rem(-Float32(0x1p127), -Float32(0x3p-149)), -Float32(0x1p-149))
@test isequal(rem(-Float32(0x1p127), Float32(0x3p-148)), -Float32(0x1p-147))
@test isequal(rem(-Float32(0x1p127), -Float32(0x3p-148)), -Float32(0x1p-147))
@test isequal(rem(-Float32(0x1p127), Float32(0x3p-126)), -Float32(0x1p-125))
@test isequal(rem(-Float32(0x1p127), -Float32(0x3p-126)), -Float32(0x1p-125))
end

end

0 comments on commit ce7a372

Please sign in to comment.