Improve accuracy of conversions from floating-point numbers (Fixes #102) #131

kimikage · 2019-10-28T15:40:41Z

This fixes a problem with the input range checking (#102).
Since this PR includes low-level operations, special attention is required. The additional tests are for range checking, not for accuracy or integrity.
This PR has little effect on the conversions to Normed{UInt8} and Normed{UInt16}.

kimikage · 2019-10-28T19:41:04Z

On the 32-bit versions, some pre-calculated constants do not agree with the dynamically calculated values.

julia> versioninfo()
Julia Version 1.4.0-DEV.378
Commit dbaa6ff660 (2019-10-28 12:55 UTC)
Platform Info:
  OS: Windows (i686-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8565U CPU @ 1.80GHz
  WORD_SIZE: 32
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)

julia> a(x) = bitstring(x*(typemax(UInt64)/(FixedPointNumbers.rawone(N11f53)))); # pre-calculated constant

julia> b(x) = bitstring((x*typemax(UInt64))/(FixedPointNumbers.rawone(N11f53))); # dynamically calculated

julia> a(1) # 1 is an identity element. On 64-bit systems, `a(1)` agrees with `b(1)`.
"0100000010100000000000000000000000000000000000000000000000000000"

julia> b(1)
"0100000010100000000000000000000000000000000000000000000000000001"

julia> Float64(BigFloat(typemax(UInt64))/BigFloat(FixedPointNumbers.rawone(N11f53))) # ideal result, but "usually", not obtained in Float64 calculations.
2048.0

julia> Float64(typemax(UInt64))/Float64(FixedPointNumbers.rawone(N11f53)) # usual result in Float64 calculations.
2048.0000000000005

julia> BigFloat(typemax(UInt64))/BigFloat(FixedPointNumbers.rawone(N11f53)) # this should be rounded to `2048.0`, so the ideal result is correct.
2048.000000000000227262653140769569055940417886493182912978948283864490173983107

Is this a problem that should be solved by modifying the tests?

kimikage · 2019-10-29T03:21:00Z

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)

kimikage · 2019-10-30T10:01:35Z

src/normed.jl

+    if T == UInt128 && f == 53
+        0 <= x <= Tf(3.777893186295717e22) || throw_converterror(U, x)
+    else
+        0 <= x <= Tf((typemax(T)-rawone(U))/rawone(U)+1) || throw_converterror(U, x)
+    end


(typemax(T)-rawone(U))/rawone(U)+1 may have slightly higher accuracy than typemax(T)/rawone(U) because the 1s in the dividend are reduced. However, this have no effect for UInt128. So, the specialization for Normed{UInt128, 53} is needed (f == 53 is one of singular points). FYI, using BigFloat hinders the constant pre-calculation and inline expansion.

Although this can be a temporary measure, I think the problem with pre-calculating constants is potentially troublesome.

codecov-io · 2019-10-30T10:29:05Z

Codecov Report

Merging #131 into master will increase coverage by 2.67%.
The diff coverage is 100%.

@@            Coverage Diff            @@
##           master    #131      +/-   ##
=========================================
+ Coverage   79.72%   82.4%   +2.67%     
=========================================
  Files           3       3              
  Lines         217     233      +16     
=========================================
+ Hits          173     192      +19     
+ Misses         44      41       -3

Impacted Files	Coverage Δ
src/FixedPointNumbers.jl	`76.92% <ø> (ø)`	⬆️
src/normed.jl	`86.23% <100%> (+5.8%)`	⬆️
src/fixed.jl	`82.6% <0%> (-0.38%)`	⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update da39318...d98d60c. Read the comment docs.

timholy

Wow, that's an impressive conversion. Does this have negative consequences for performance?

timholy · 2019-10-30T19:05:59Z

src/normed.jl

@@ -41,27 +41,61 @@ rawone(v) = reinterpret(one(v))
 function Normed{T,f}(x::Normed{T2}) where {T <: Unsigned,T2 <: Unsigned,f}
    U = Normed{T,f}
    y = round((rawone(U)/rawone(x))*reinterpret(x))
-    (0 <= y) & (y <= typemax(T)) || throw_converterror(U, x)
+    0 <= y <= typemax(T) || throw_converterror(U, x)


These are actually subtly different, and at least at one point there was a performance difference (doesn't seem to be one now, at least not on my current machine):

julia> function f1(x) 0.25 <= x <= 0.75 || throw(DomainError(x)) return x end f1 (generic function with 1 method) julia> function f2(x) ((0.25 <= x) & (x <= 0.75)) || throw(DomainError(x)) return x end f2 (generic function with 1 method) julia> @code_lowered f1(0.5) CodeInfo( 1 ─ %1 = 0.25 <= x └── goto #3 if not %1 2 ─ @_3 = x <= 0.75 └── goto #4 3 ─ @_3 = false 4 ┄ %6 = @_3 └── goto #6 if not %6 5 ─ goto #7 6 ─ %9 = Main.DomainError(x) └── Main.throw(%9) 7 ┄ return x ) julia> @code_lowered f2(0.5) CodeInfo( 1 ─ %1 = 0.25 <= x │ %2 = x <= 0.75 │ %3 = %1 & %2 └── goto #3 if not %3 2 ─ goto #4 3 ─ %6 = Main.DomainError(x) └── Main.throw(%6) 4 ┄ return x ) julia> using BenchmarkTools julia> @btime f1(x) setup=(x=rand()/10 + 0.3) 1.980 ns (0 allocations: 0 bytes) 0.34079130696459065 julia> @btime f2(x) setup=(x=rand()/10 + 0.3) 1.977 ns (0 allocations: 0 bytes) 0.30934306198706746

Yes. I guess, using & instead of && may reduce conditional branches. However, throw_converterror(U, x) is annotated with @noinline and throws the exception. This means that we cannot avoid conditional branches, and this is not SIMD-suitable. Depending on the type (and CPU), the comparison with zero can be specialized (this is the reason why I did not use zero()).

Of course, I follow the law of "if it works, don't touch it".:smile:

Yeah, you definitely need one branch, but one is usually better than two. For example when Matt & I were developing the array infrastructure we noticed it was better to check bounds on all indexes and combine with &, and then throw a BoundsError if necessary, rather than having one branch per dimension.

But this is definitely something that depends on the CPU and its ability to predict branches, so it can even be hard to benchmark. I would generally suspect that fewer branches is a good thing.

That's right. Nevertheless, I think it is another good practice to avoid something tricky and leave it to the compiler (LLVM) now, because this is usually not a pair of integer comparisons.

However, I do never intend to dispute about the comparison manner, which has no definitive answer. Instead I would care about the /, but it is another matter.

kimikage · 2019-10-31T02:38:31Z

Does this have negative consequences for performance?

Unfortunately, my software floating-point calculation is ~2.5x slower.
The following are the benchmark results based on the (inverse) Vec4 version of #129 (comment)

x86_64	before	after
`Float32`->`N0f8`	15.999 μs	14.099 μs
`Float32`->`N8f24`	20.000 μs	48.199 μs
`Float32`->`N0f32`	20.100 μs	49.100 μs
`Float64`->`N0f8`	15.999 μs	13.999 μs
`Float64`->`N0f32`	16.099 μs	14.299 μs
`Float64`->`N11f53`	41.899 μs	49.500 μs
`Float64`->`N0f64`	40.999 μs	68.900 μs

However, we cannot measure the conversion speed of what cannot be converted. As implied above, in applications where the speed is needed, rem is used and the conversions with range checking should not be used.

Edit:
As I mentioned in #125 (comment), * and / already have performance problems.

FixedPointNumbers.jl/src/normed.jl

Lines 93 to 94 in da39318

    
           *(x::T, y::T) where {T <: Normed} = convert(T,convert(floattype(T), x)*convert(floattype(T), y)) 
        
           /(x::T, y::T) where {T <: Normed} = convert(T,convert(floattype(T), x)/convert(floattype(T), y))

At least for *, I want to avoid the floating-point arithmetic. Although the discussion of overflow (#41) is undecided, since currently +, - and Fixed's * are the wraparound style, even breaking changes are reasonable.

The order of bug fixes is not so important, but this problem (#102) disturbs the checks or benchmarks for the modifications to mitigate the issue #129, and the modifications for #129 may be helpful for the PR #123 (issue #120).

timholy · 2019-10-31T09:38:08Z

Thanks for posting those numbers. Correctness is infinitely more important than speed, so I would support merging this.

kimikage · 2019-11-01T02:08:03Z

For the lines based on the current source code, I reverted the comparison manner. (cf. #131 (comment))

In the new conversion method, I used the a <= x <= b style. This decision is based on my policy and the rough benchmark. Honestly, the reason is that it is shorter in source codes. 😃
Although the inconsistency makes me somewhat itchy, I think it will also make more people notice this comparison issue.

timholy · 2019-11-01T09:09:57Z

I'm fine with whatever you choose here. Very grateful for this change!

kimikage · 2019-11-01T12:29:15Z

I am ready to rebase and merge this PR, but not yet ready for other issues. When do you bump the version number up?

timholy · 2019-11-01T14:38:49Z

Version numbers are cheap: any time it's potentially breaking, we should bump the minor version number. If the other fixes will follow quickly I'd say let's consolidate the churn and wait to tag a release until it's all done. If not, let's get a release out for the fixes we have.

As the community has moved to putting upper version bounds on packages, it's also more likely that breaking changes won't break things for users.

This fixes a problem with the input range checking.

kimikage · 2019-11-02T03:25:58Z

If this fix causes problems, consider using rem (% Nxfy) instead.

kimikage force-pushed the issue102 branch from f8b893c to a41e65c Compare October 29, 2019 00:45

kimikage force-pushed the issue102 branch from a41e65c to d98d60c Compare October 30, 2019 09:57

kimikage commented Oct 30, 2019

View reviewed changes

timholy reviewed Oct 30, 2019

View reviewed changes

timholy approved these changes Oct 31, 2019

View reviewed changes

kimikage force-pushed the issue102 branch from d98d60c to 118ad55 Compare November 1, 2019 01:12

kimikage added 2 commits November 2, 2019 11:49

Fix comments

2886070

Improve accuracy of conversions from floating-point numbers

547eb95

This fixes a problem with the input range checking.

kimikage force-pushed the issue102 branch from 118ad55 to 547eb95 Compare November 2, 2019 02:51

kimikage merged commit 70ae1d6 into JuliaMath:master Nov 2, 2019

kimikage deleted the issue102 branch November 2, 2019 03:22

kimikage mentioned this pull request Nov 2, 2019

can not convert float to the same or more number of bit fixed point #102

Closed

kimikage mentioned this pull request Nov 19, 2019

Optimizing Normed -> Normed conversions #140

Closed

kimikage mentioned this pull request Dec 15, 2019

Rounding errors in rem of Normed #150

Closed

kimikage mentioned this pull request Dec 26, 2019

[RFC] Acceptable input ranges #156

Closed

kimikage mentioned this pull request Mar 3, 2020

Improve HSV/HSL/HSI conversions (Fixes #378, #379) JuliaGraphics/Colors.jl#407

Merged

kimikage mentioned this pull request Apr 21, 2020

Small improvements to docstrings and error messages JuliaGraphics/ColorTypes.jl#183

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Improve accuracy of conversions from floating-point numbers (Fixes #102) #131

Improve accuracy of conversions from floating-point numbers (Fixes #102) #131

kimikage commented Oct 28, 2019 •

edited

Loading

kimikage commented Oct 28, 2019 •

edited

Loading

kimikage commented Oct 29, 2019

kimikage Oct 30, 2019

codecov-io commented Oct 30, 2019

timholy left a comment

timholy Oct 30, 2019

kimikage Oct 31, 2019 •

edited

Loading

timholy Oct 31, 2019

kimikage Oct 31, 2019

kimikage commented Oct 31, 2019 •

edited

Loading

timholy commented Oct 31, 2019

kimikage commented Nov 1, 2019

timholy commented Nov 1, 2019

kimikage commented Nov 1, 2019 •

edited

Loading

timholy commented Nov 1, 2019

kimikage commented Nov 2, 2019

Improve accuracy of conversions from floating-point numbers (Fixes #102) #131

Improve accuracy of conversions from floating-point numbers (Fixes #102) #131

Conversation

kimikage commented Oct 28, 2019 • edited Loading

kimikage commented Oct 28, 2019 • edited Loading

kimikage commented Oct 29, 2019

kimikage Oct 30, 2019

Choose a reason for hiding this comment

codecov-io commented Oct 30, 2019

Codecov Report

timholy left a comment

Choose a reason for hiding this comment

timholy Oct 30, 2019

Choose a reason for hiding this comment

kimikage Oct 31, 2019 • edited Loading

Choose a reason for hiding this comment

timholy Oct 31, 2019

Choose a reason for hiding this comment

kimikage Oct 31, 2019

Choose a reason for hiding this comment

kimikage commented Oct 31, 2019 • edited Loading

timholy commented Oct 31, 2019

kimikage commented Nov 1, 2019

timholy commented Nov 1, 2019

kimikage commented Nov 1, 2019 • edited Loading

timholy commented Nov 1, 2019

kimikage commented Nov 2, 2019

kimikage commented Oct 28, 2019 •

edited

Loading

kimikage commented Oct 28, 2019 •

edited

Loading

kimikage Oct 31, 2019 •

edited

Loading

kimikage commented Oct 31, 2019 •

edited

Loading

kimikage commented Nov 1, 2019 •

edited

Loading