Implement modf in Julia (JuliaLang#42392)

Currently `modf` calls into libm, and is one of the few remaining
holdouts that use libm. This commit reimplements it in Julia using a
simple, naive implementation that matches the IEEE behavior and has
performance on par with a port of the musl implementation and with the
existing implementation.

Co-Authored-By: Simon Byrne <[email protected]>
Co-Authored-By: woclass <[email protected]>

base/math.jl

@@ -932,16 +932,11 @@ julia> modf(-3.5)
 """
 modf(x) = isinf(x) ? (flipsign(zero(x), x), x) : (rem(x, one(x)), trunc(x))
 
-function modf(x::Float32)
- temp = Ref{Float32}()
- f = ccall((:modff, libm), Float32, (Float32, Ptr{Float32}), x, temp)
- f, temp[]
-end
-
-function modf(x::Float64)
- temp = Ref{Float64}()
- f = ccall((:modf, libm), Float64, (Float64, Ptr{Float64}), x, temp)
- f, temp[]
+function modf(x::T) where T<:IEEEFloat
+ isinf(x) && return (copysign(zero(T), x), x)
+ ix = trunc(x)
+ rx = copysign(x - ix, x)
+ return (rx, ix)
 end
 
 @inline function ^(x::Float64, y::Float64)

test/math.jl

@@ -663,14 +663,14 @@ end
 end
 
 @testset "modf" begin
- @testset "$elty" for elty in (Float16, Float32, Float64)
- @test modf( convert(elty,1.2) )[1] ≈ convert(elty,0.2)
- @test modf( convert(elty,1.2) )[2] ≈ convert(elty,1.0)
- @test modf( convert(elty,1.0) )[1] ≈ convert(elty,0.0)
- @test modf( convert(elty,1.0) )[2] ≈ convert(elty,1.0)
- @test isequal(modf( convert(elty,-Inf) ), (-0.0, -Inf))
- @test isequal(modf( convert(elty,Inf) ), (0.0, Inf))
- @test isequal(modf( convert(elty,NaN) ), (NaN, NaN))
+ @testset "$T" for T in (Float16, Float32, Float64)
+ @test modf(T(1.25)) === (T(0.25), T(1.0))
+ @test modf(T(1.0)) === (T(0.0), T(1.0))
+ @test modf(T(-Inf)) === (T(-0.0), T(-Inf))
+ @test modf(T(Inf)) === (T(0.0), T(Inf))
+ @test modf(T(NaN)) === (T(NaN), T(NaN))
+ @test modf(T(-0.0)) === (T(-0.0), T(-0.0))
+ @test modf(T(-1.0)) === (T(-0.0), T(-1.0))
  end
 end