Skip to content

Commit cf5ae03

Browse files
gbaraldiararslan
gbaraldi
and
ararslan
authored
Add native julia fmod (#47501)
* Add native julia rem Co-authored-by: Alex Arslan <[email protected]>
1 parent c3bf180 commit cf5ae03

File tree

2 files changed

+227
-2
lines changed

2 files changed

+227
-2
lines changed

‎base/float.jl‎

Lines changed: 104 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -101,6 +101,8 @@ exponent_one(::Type{Float16}) = 0x3c00
101101
exponent_half(::Type{Float16}) = 0x3800
102102
significand_mask(::Type{Float16}) = 0x03ff
103103

104+
mantissa(x::T) where {T} = reinterpret(Unsigned, x) & significand_mask(T)
105+
104106
for T in (Float16, Float32, Float64)
105107
@eval significand_bits(::Type{$T}) = $(trailing_ones(significand_mask(T)))
106108
@eval exponent_bits(::Type{$T}) = $(sizeof(T)*8 - significand_bits(T) - 1)
@@ -414,9 +416,109 @@ muladd(x::T, y::T, z::T) where {T<:IEEEFloat} = muladd_float(x, y, z)
414416
# TODO: faster floating point fld?
415417
# TODO: faster floating point mod?
416418

417-
rem(x::T, y::T) where {T<:IEEEFloat} = rem_float(x, y)
419+
function unbiased_exponent(x::T) where {T<:IEEEFloat}
420+
return (reinterpret(Unsigned, x) & exponent_mask(T)) >> significand_bits(T)
421+
end
422+
423+
function explicit_mantissa_noinfnan(x::T) where {T<:IEEEFloat}
424+
m = mantissa(x)
425+
issubnormal(x) || (m |= significand_mask(T) + uinttype(T)(1))
426+
return m
427+
end
428+
429+
function _to_float(number::U, ep) where {U<:Unsigned}
430+
F = floattype(U)
431+
S = signed(U)
432+
epint = unsafe_trunc(S,ep)
433+
lz::signed(U) = unsafe_trunc(S, Core.Intrinsics.ctlz_int(number) - U(exponent_bits(F)))
434+
number <<= lz
435+
epint -= lz
436+
bits = U(0)
437+
if epint >= 0
438+
bits = number & significand_mask(F)
439+
bits |= ((epint + S(1)) << significand_bits(F)) & exponent_mask(F)
440+
else
441+
bits = (number >> -epint) & significand_mask(F)
442+
end
443+
return reinterpret(F, bits)
444+
end
445+
446+
@assume_effects :terminates_locally :nothrow function rem_internal(x::T, y::T) where {T<:IEEEFloat}
447+
xuint = reinterpret(Unsigned, x)
448+
yuint = reinterpret(Unsigned, y)
449+
if xuint <= yuint
450+
if xuint < yuint
451+
return x
452+
end
453+
return zero(T)
454+
end
455+
456+
e_x = unbiased_exponent(x)
457+
e_y = unbiased_exponent(y)
458+
# Most common case where |y| is "very normal" and |x/y| < 2^EXPONENT_WIDTH
459+
if e_y > (significand_bits(T)) && (e_x - e_y) <= (exponent_bits(T))
460+
m_x = explicit_mantissa_noinfnan(x)
461+
m_y = explicit_mantissa_noinfnan(y)
462+
d = urem_int((m_x << (e_x - e_y)), m_y)
463+
iszero(d) && return zero(T)
464+
return _to_float(d, e_y - uinttype(T)(1))
465+
end
466+
# Both are subnormals
467+
if e_x == 0 && e_y == 0
468+
return reinterpret(T, urem_int(xuint, yuint) & significand_mask(T))
469+
end
470+
471+
m_x = explicit_mantissa_noinfnan(x)
472+
e_x -= uinttype(T)(1)
473+
m_y = explicit_mantissa_noinfnan(y)
474+
lz_m_y = uinttype(T)(exponent_bits(T))
475+
if e_y > 0
476+
e_y -= uinttype(T)(1)
477+
else
478+
m_y = mantissa(y)
479+
lz_m_y = Core.Intrinsics.ctlz_int(m_y)
480+
end
481+
482+
tz_m_y = Core.Intrinsics.cttz_int(m_y)
483+
sides_zeroes_cnt = lz_m_y + tz_m_y
484+
485+
# n>0
486+
exp_diff = e_x - e_y
487+
# Shift hy right until the end or n = 0
488+
right_shift = min(exp_diff, tz_m_y)
489+
m_y >>= right_shift
490+
exp_diff -= right_shift
491+
e_y += right_shift
492+
# Shift hx left until the end or n = 0
493+
left_shift = min(exp_diff, uinttype(T)(exponent_bits(T)))
494+
m_x <<= left_shift
495+
exp_diff -= left_shift
496+
497+
m_x = urem_int(m_x, m_y)
498+
iszero(m_x) && return zero(T)
499+
iszero(exp_diff) && return _to_float(m_x, e_y)
500+
501+
while exp_diff > sides_zeroes_cnt
502+
exp_diff -= sides_zeroes_cnt
503+
m_x <<= sides_zeroes_cnt
504+
m_x = urem_int(m_x, m_y)
505+
end
506+
m_x <<= exp_diff
507+
m_x = urem_int(m_x, m_y)
508+
return _to_float(m_x, e_y)
509+
end
510+
511+
function rem(x::T, y::T) where {T<:IEEEFloat}
512+
if isfinite(x) && !iszero(x) && isfinite(y) && !iszero(y)
513+
return copysign(rem_internal(abs(x), abs(y)), x)
514+
elseif isinf(x) || isnan(y) || iszero(y) # y can still be Inf
515+
return T(NaN)
516+
else
517+
return x
518+
end
519+
end
418520

419-
function mod(x::T, y::T) where T<:AbstractFloat
521+
function mod(x::T, y::T) where {T<:AbstractFloat}
420522
r = rem(x,y)
421523
if r == 0
422524
copysign(r,y)

‎test/numbers.jl‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2929,3 +2929,126 @@ end
29292929
@test false == ceil(Bool, -0.7)
29302930
end
29312931
end
2932+
2933+
@testset "modf" begin
2934+
@testset "remd" begin
2935+
denorm_min = nextfloat(0.0)
2936+
minfloat = floatmin(Float64)
2937+
maxfloat = floatmax(Float64)
2938+
values = [3.0,denorm_min,-denorm_min, minfloat,
2939+
-minfloat, maxfloat, -maxfloat]
2940+
# rem (0, y) == 0 for y != 0.
2941+
for val in values
2942+
@test isequal(rem(0.0, val), 0.0)
2943+
end
2944+
# rem (-0, y) == -0 for y != 0.
2945+
for val in values
2946+
@test isequal(rem(-0.0, val), -0.0)
2947+
end
2948+
# rem (+Inf, y) == NaN
2949+
values2 = [3.0,-1.1,0.0,-0.0,denorm_min,minfloat,
2950+
maxfloat,Inf,-Inf]
2951+
for val in values2
2952+
@test isequal(rem(Inf, val), NaN)
2953+
end
2954+
# rem (-Inf, y) == NaN
2955+
for val in values2
2956+
@test isequal(rem(-Inf, val), NaN)
2957+
end
2958+
# rem (x, +0) == NaN
2959+
values3 = values2[begin:end-2]
2960+
for val in values3
2961+
@test isequal(rem(val, 0.0), NaN)
2962+
end
2963+
# rem (x, -0) == NaN
2964+
for val in values3
2965+
@test isequal(rem(val, -0.0), NaN)
2966+
end
2967+
# rem (x, +Inf) == x for x not infinite.
2968+
@test isequal(rem(0.0, Inf), 0.0)
2969+
@test isequal(rem(-0.0, Inf), -0.0)
2970+
@test isequal(rem(denorm_min, Inf), denorm_min)
2971+
@test isequal(rem(minfloat, Inf), minfloat)
2972+
@test isequal(rem(maxfloat, Inf), maxfloat)
2973+
@test isequal(rem(3.0, Inf), 3.0)
2974+
# rem (x, -Inf) == x for x not infinite.
2975+
@test isequal(rem(0.0, -Inf), 0.0)
2976+
@test isequal(rem(-0.0, -Inf), -0.0)
2977+
@test isequal(rem(denorm_min, -Inf), denorm_min)
2978+
@test isequal(rem(minfloat, -Inf), minfloat)
2979+
@test isequal(rem(maxfloat, -Inf), maxfloat)
2980+
@test isequal(rem(3.0, -Inf), 3.0)
2981+
#NaN tests
2982+
@test isequal(rem(0.0, NaN), NaN)
2983+
@test isequal(rem(1.0, NaN), NaN)
2984+
@test isequal(rem(Inf, NaN), NaN)
2985+
@test isequal(rem(NaN, 0.0), NaN)
2986+
@test isequal(rem(NaN, 1.0), NaN)
2987+
@test isequal(rem(NaN, Inf), NaN)
2988+
@test isequal(rem(NaN, NaN), NaN)
2989+
#Sign tests
2990+
@test isequal(rem(6.5, 2.25), 2.0)
2991+
@test isequal(rem(-6.5, 2.25), -2.0)
2992+
@test isequal(rem(6.5, -2.25), 2.0)
2993+
@test isequal(rem(-6.5, -2.25), -2.0)
2994+
values4 = [maxfloat,-maxfloat,minfloat,-minfloat,
2995+
denorm_min, -denorm_min]
2996+
for val in values4
2997+
@test isequal(rem(maxfloat,val), 0.0)
2998+
end
2999+
for val in values4
3000+
@test isequal(rem(-maxfloat,val), -0.0)
3001+
end
3002+
@test isequal(rem(minfloat, maxfloat), minfloat)
3003+
@test isequal(rem(minfloat, -maxfloat), minfloat)
3004+
values5 = values4[begin+2:end]
3005+
for val in values5
3006+
@test isequal(rem(minfloat,val), 0.0)
3007+
end
3008+
@test isequal(rem(-minfloat, maxfloat), -minfloat)
3009+
@test isequal(rem(-minfloat, -maxfloat), -minfloat)
3010+
for val in values5
3011+
@test isequal(rem(-minfloat,val), -0.0)
3012+
end
3013+
values6 = values4[begin:end-2]
3014+
for val in values6
3015+
@test isequal(rem(denorm_min,val), denorm_min)
3016+
end
3017+
@test isequal(rem(denorm_min, denorm_min), 0.0)
3018+
@test isequal(rem(denorm_min, -denorm_min), 0.0)
3019+
for val in values6
3020+
@test isequal(rem(-denorm_min,val), -denorm_min)
3021+
end
3022+
@test isequal(rem(-denorm_min, denorm_min), -0.0)
3023+
@test isequal(rem(-denorm_min, -denorm_min), -0.0)
3024+
#Max value tests
3025+
values7 = [0x3p-1074,-0x3p-1074,0x3p-1073,-0x3p-1073]
3026+
for val in values7
3027+
@test isequal(rem(0x1p1023,val), 0x1p-1073)
3028+
end
3029+
@test isequal(rem(0x1p1023, 0x3p-1022), 0x1p-1021)
3030+
@test isequal(rem(0x1p1023, -0x3p-1022), 0x1p-1021)
3031+
for val in values7
3032+
@test isequal(rem(-0x1p1023,val), -0x1p-1073)
3033+
end
3034+
@test isequal(rem(-0x1p1023, 0x3p-1022), -0x1p-1021)
3035+
@test isequal(rem(-0x1p1023, -0x3p-1022), -0x1p-1021)
3036+
3037+
end
3038+
3039+
@testset "remf" begin
3040+
@test isequal(rem(Float32(0x1p127), Float32(0x3p-149)), Float32(0x1p-149))
3041+
@test isequal(rem(Float32(0x1p127), -Float32(0x3p-149)), Float32(0x1p-149))
3042+
@test isequal(rem(Float32(0x1p127), Float32(0x3p-148)), Float32(0x1p-147))
3043+
@test isequal(rem(Float32(0x1p127), -Float32(0x3p-148)), Float32(0x1p-147))
3044+
@test isequal(rem(Float32(0x1p127), Float32(0x3p-126)), Float32(0x1p-125))
3045+
@test isequal(rem(Float32(0x1p127), -Float32(0x3p-126)), Float32(0x1p-125))
3046+
@test isequal(rem(-Float32(0x1p127), Float32(0x3p-149)), -Float32(0x1p-149))
3047+
@test isequal(rem(-Float32(0x1p127), -Float32(0x3p-149)), -Float32(0x1p-149))
3048+
@test isequal(rem(-Float32(0x1p127), Float32(0x3p-148)), -Float32(0x1p-147))
3049+
@test isequal(rem(-Float32(0x1p127), -Float32(0x3p-148)), -Float32(0x1p-147))
3050+
@test isequal(rem(-Float32(0x1p127), Float32(0x3p-126)), -Float32(0x1p-125))
3051+
@test isequal(rem(-Float32(0x1p127), -Float32(0x3p-126)), -Float32(0x1p-125))
3052+
end
3053+
3054+
end

0 commit comments

Comments
 (0)