Skip to content
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

Jump to bottom

Add native julia fmod #47501

Merged
merged 21 commits into from
Dec 6, 2022
Merged

Add native julia fmod #47501

Show file tree
Hide file tree
Changes from 10 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
124 changes: 122 additions & 2 deletions base/float.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,18 +88,26 @@ exponent_mask(::Type{Float64}) = 0x7ff0_0000_0000_0000
exponent_one(::Type{Float64}) = 0x3ff0_0000_0000_0000
exponent_half(::Type{Float64}) = 0x3fe0_0000_0000_0000
significand_mask(::Type{Float64}) = 0x000f_ffff_ffff_ffff
mantissa_width(::Type{Float64}) = UInt64(52)
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
exponent_width(::Type{Float64}) = UInt64(11)

sign_mask(::Type{Float32}) = 0x8000_0000
exponent_mask(::Type{Float32}) = 0x7f80_0000
exponent_one(::Type{Float32}) = 0x3f80_0000
exponent_half(::Type{Float32}) = 0x3f00_0000
significand_mask(::Type{Float32}) = 0x007f_ffff
mantissa_width(::Type{Float32}) = UInt32(23)
exponent_width(::Type{Float32}) = UInt32(8)

sign_mask(::Type{Float16}) = 0x8000
exponent_mask(::Type{Float16}) = 0x7c00
exponent_one(::Type{Float16}) = 0x3c00
exponent_half(::Type{Float16}) = 0x3800
significand_mask(::Type{Float16}) = 0x03ff
mantissa_width(::Type{Float16}) = UInt16(10)
exponent_width(::Type{Float16}) = UInt16(5)

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)))
Expand Down Expand Up @@ -414,9 +422,121 @@ 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}
exp = (reinterpret(Unsigned, x) & exponent_mask(T)) >> mantissa_width(T)
return exp
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
end

function explicit_mantissa_noinfnan(x::T) where {T<:IEEEFloat}
subnormal = !issubnormal(x) ? significand_mask(T) + uinttype(T)(1) : uinttype(T)(0)
return subnormal | (significand_mask(T) & reinterpret(Unsigned, x))
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
end

function make_value(number::T, ep) where {T<:Unsigned}
ararslan marked this conversation as resolved.
Show resolved Hide resolved
Tfloat = floattype(T)
Tint = signed(T)
epint = unsafe_trunc(Tint,ep)
lz::signed(T) = Core.Intrinsics.ctlz_int(number) - exponent_width(Tfloat)
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
number <<= lz
epint -= lz
bits = T(0)
if epint >= T(0)
mant = number & significand_mask(Tfloat) #this might not be necessary
bits |= mant
epint = ((epint + Tint(1)) << mantissa_width(Tfloat)) & exponent_mask(Tfloat)
bits |= epint
else
number = number >> -epint
mant = number & significand_mask(Tfloat)
bits |= mant
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
end
return reinterpret(Tfloat, bits)
end

function fmod_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 > mantissa_width(T) && (e_x - e_y) <= exponent_width(T)
m_x = explicit_mantissa_noinfnan(x)
m_y = explicit_mantissa_noinfnan(y)
d = (e_x == e_y) ? (m_x - m_y) : (m_x << (e_x - e_y)) % m_y
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
iszero(d) && return zero(T)
return make_value(d, e_y - uinttype(T)(1))
end
# Both are subnormals
if e_x == uinttype(T)(0) && e_y == uinttype(T)(0)
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
bits = uinttype(T)(0)
number = xuint % yuint
mant = number & significand_mask(T) #this might not be necessary
bits |= mant
return reinterpret(T, bits)
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
end

m_x = explicit_mantissa_noinfnan(x)
e_x -= uinttype(T)(1)
m_y = explicit_mantissa_noinfnan(y)
lz_m_y = exponent_width(T)
if e_y > uinttype(T)(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::uinttype(T) = e_x - e_y
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
# Shift hy right until the end or n = 0
right_shift = exp_diff < tz_m_y ? exp_diff : tz_m_y
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
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, exponent_width(T))
m_x <<= left_shift
exp_diff -= left_shift

m_x %= m_y
iszero(m_x) && return zero(T)
iszero(exp_diff) && return make_value(m_x, e_y)

while exp_diff > sides_zeroes_cnt
exp_diff -= sides_zeroes_cnt
m_x <<= sides_zeroes_cnt
m_x %= m_y
end
m_x <<= exp_diff
m_x %= m_y
return make_value(m_x, e_y)
end

function fmod(x::T, y::T) where {T<:IEEEFloat}
if isfinite(x) && !iszero(x) && isfinite(y) && !iszero(y)
xabs = abs(x)
yabs = abs(y)
res = fmod_internal(xabs, yabs)
return copysign(res,x)
end
(isnan(x) || isnan(y)) && return T(NaN)
(isinf(x) || iszero(y)) && return T(NaN)
return x
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
end

rem(x::T, y::T) where {T<:IEEEFloat} = fmod(x, y)
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved

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)
gbaraldi marked this conversation as resolved.
Show resolved Hide resolved
@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