Skip to content

Add native julia fmod #47501

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 21 commits into from
Dec 6, 2022
Merged

Add native julia fmod #47501

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
e29c20a
Add native julia rem
gbaraldi Nov 8, 2022
64751b5
Add Type checks and tests
gbaraldi Nov 8, 2022
765b0cf
Apply fixes and simplify tests.
gbaraldi Nov 9, 2022
cba6228
FIx a couple details
gbaraldi Nov 9, 2022
018c102
Apply fixes and simplify tests.
gbaraldi Nov 9, 2022
6301997
Also check if x is equal to 0
gbaraldi Nov 9, 2022
bdf07e7
Use reinterpret directly
gbaraldi Nov 10, 2022
4c01d97
Update base/float.jl
gbaraldi Nov 10, 2022
f95d60c
Update base/float.jl
gbaraldi Nov 10, 2022
4fa3c6c
Update base/float.jl
gbaraldi Nov 12, 2022
3dab2d9
Update base/float.jl
gbaraldi Nov 13, 2022
46a0312
Update base/float.jl
gbaraldi Nov 13, 2022
1334782
Update base/float.jl
gbaraldi Nov 13, 2022
7a9d9b0
Update base/float.jl
gbaraldi Nov 13, 2022
39b7673
Simplify and apply suggestions
gbaraldi Nov 13, 2022
11fbbc4
Make tests pass again
gbaraldi Nov 13, 2022
7f09718
Small changes that make it a bit faster
gbaraldi Nov 14, 2022
650bbf8
Add effects assumption
gbaraldi Nov 14, 2022
3358b80
Small changes
gbaraldi Nov 15, 2022
226c4a5
Merge branch 'master' into fast-rem
gbaraldi Nov 25, 2022
bff1ae5
Merge branch 'master' into fast-rem
ViralBShah Dec 6, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 104 additions & 2 deletions base/float.jl
View file
Open in desktop
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
Open in desktop
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