Rounding errors in rem of Normed

[kimikage]

As I mentioned in #102, PR #131 did not change the rem, so this is a reminder.

To avoid impact on the performance, this PR does not change the rem.

FixedPointNumbers.jl/src/normed.jl

Lines 67 to 70 in da39318

	rem(x::T, ::Type{T}) where {T <: Normed} = x
	rem(x::Normed, ::Type{T}) where {T <: Normed} = reinterpret(T, _unsafe_trunc(rawtype(T), round((rawone(T)/rawone(x))*reinterpret(x))))
	rem(x::Real, ::Type{T}) where {T <: Normed} = reinterpret(T, _unsafe_trunc(rawtype(T), round(rawone(T)*x)))
	rem(x::Float16, ::Type{T}) where {T <: Normed} = rem(Float32(x), T) # avoid overflow

Therefore, this does not fix the problems such as:

julia> using Colors

julia> white = parse(RGB{Float32}, "white")
RGB{Float32}(1.0f0,1.0f0,1.0f0)

julia> convert(RGB{N0f32}, white)
RGB{N0f32}(0.0,0.0,0.0)

Originally posted by @kimikage in #131 (comment)

julia> using FixedPointNumbers

julia> for f = 1:32
           F = Normed{UInt32,f}
           rem_one = 1.0f0 % F
           rem_one == one(F) || @show rem_one
       end
rem_one = 1.00000003N7f25
rem_one = 1.00000001N6f26
rem_one = 1.000000007N5f27
rem_one = 1.000000004N4f28
rem_one = 1.000000002N3f29
rem_one = 1.0000000009N2f30
rem_one = 1.0000000005N1f31
rem_one = 0.0N0f32

[kimikage]

Perhaps the following is more accurate than the current implementation, and less accurate than PR #131.

function Base.rem(x::Float32, ::Type{T}) where {f, T <: Normed{UInt32,f}}
    f <= 24 && return reinterpret(T, unsafe_trunc(UInt32, round(rawone(T)*x)))
    r = unsafe_trunc(UInt32, round(Int32, x * Float32(0x1p24)))
    reinterpret(T, r << UInt8(f - 24) - r >> 0x18)
end

[kimikage]

BTW, I reported this issue because I noticed that the documents of Colors.jl include the example with N0f32.
JuliaGraphics/Colors.jl#324
http://juliagraphics.github.io/Colors.jl/latest/namedcolors/

I think it would be easier to fix the bug than add a note to the document.

[timholy]

Definitely worth fixing. Though I can't imagine why the docs use N0f32, from a perceptual standpoint that's serious overkill. I posted a query about the reasons for that choice in JuliaGraphics/Colors.jl#324

[kimikage]

minor adjustment:

function rem(x::Float32, ::Type{T}) where {f, T <: Normed{UInt32,f}}
    f <= 24 && return reinterpret(T, _unsafe_trunc(UInt32, round(rawone(T) * x)))
    r = _unsafe_trunc(UInt32, round(x * @f32(0x1p24)))
    reinterpret(T, r << UInt8(f - 24) - unsigned(signed(r) >> 0x18))
end
function rem(x::Float64, ::Type{T}) where {f, T <: Normed{UInt64,f}}
    f <= 53 && return reinterpret(T, _unsafe_trunc(UInt64, round(rawone(T) * x)))
    r = _unsafe_trunc(UInt64, round(x * 0x1p53))
    reinterpret(T, r << UInt8(f - 53) - unsigned(signed(r) >> 0x35))
end

However, similar problem remains due to the overflow instead of rounding error.

julia> 64.0f0 % N0f32 # OK
0.9999999853N0f32

julia> 127.99999f0 % N0f32 # OK
0.9999923413N0f32

julia> 128.0f0 % N0f32 # internal overflow
2.98e-8N0f32

julia> 128.1f0 % N0f32 # (somewhat) low accuracy
0.1000061333N0f32

Since % T is often used in constructors (cf. ColorTypes.jl), the following benchmark is not so practical.

julia> @btime x .% N8f24 setup=(x=rand(Float32,64,64));
  3.962 μs (1 allocation: 16.13 KiB)

julia> @btime x .% N0f32 setup=(x=rand(Float32,64,64));
  4.050 μs (1 allocation: 16.13 KiB)

julia> @btime x .% N11f53 setup=(x=rand(Float64,64,64));
  6.080 μs (2 allocations: 32.08 KiB)

julia> @btime x .% N0f64 setup=(x=rand(Float64,64,64));
  6.860 μs (2 allocations: 32.08 KiB)