Add native julia fmod

base/float.jl

@@ -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)
+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)))
@@ -414,9 +422,122 @@ 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
+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))
+end
+
+function make_value(number::T, ep) where {T<:Unsigned}
+    Tfloat = floattype(T)
+    Tint = signed(T)
+    epint = unsafe_trunc(Tint,ep)
+    lz::signed(T) = Core.Intrinsics.ctlz_int(number) - exponent_width(Tfloat)
+    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
+    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
+        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)
+        bits = uinttype(T)(0)
+        number = xuint % yuint
+        mant = number & significand_mask(T) #this might not be necessary
+        bits |= mant
+        return reinterpret(T, bits)
+    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
+    # Shift hy right until the end or n = 0
+    right_shift = exp_diff < tz_m_y ? 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 = exp_diff < exponent_width(T) ? 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 (!isnan(x) && isfinite(x) && !iszero(x) && !isnan(y) && 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)
+    isinf(y) && return x
+    return x
+end
+
+rem(x::T, y::T) where {T<:IEEEFloat} = fmod(x, y)
 
-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)

test/numbers.jl

@@ -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