Improve accuracy of rem with Normed types (e.g. ::Float32 % N0f32)

[kimikage]

This fixes #150.

Float32 and Float64 versions can be unified, but the unified method may be difficult to read. So, I implemented them into two methods.
As I mentioned in #150 (comment), there are still numerical errors.

[kimikage]

Subsequent commits: master...kimikage:ctor_fixed

I want to modify some constructor-style conversion methods of Normed (e.g. issue #140), but that is not a bug fix.
So, what about releasing v0.8.0 after merging ctor_fixed?
In the following v0.9.0, let's change the arithmetic functions, i.e. let's start #152.

[timholy]

I'm not certain I understand the strategy, so some comments would be useful. Does this essentially correspond to an affine rescaling? If so does it fix values near 1.0 at the expense of others? (The answers to these questions don't block merging, but let's document this a bit more.)

[kimikage]

Does this essentially correspond to an affine rescaling?

No, this is based on the "pure" linear transformation, so this is effective for almost all values other than the exception mentioned above.

First, the root of this problem is that while the width of UInt32 is 32-bit, the mantissa of Float32 is only 23+1-bit. This means that rawone of the type with f > 24 cannot be exactly represented in Float32.

julia> rawone(Normed{UInt32, 24}) == Float32(0xFFFFFF)
true

julia> rawone(Normed{UInt32, 25}) == Float32(0x1FFFFFF)
false

Therefore, if f <= 24, we use the usual method.

    f <= 24 && return reinterpret(N, _unsafe_trunc(UInt32, round(rawone(N) * x)))

Well, the rawone is 2^f - 1, i.e.

x * rawone == x * 2^f  -  x
           == x * 2^(f + k - k)  -  x
           == x * 2^(f + k) * 2^-k  -  x * 2^(f + k) * 2^-(f + k)

where k is an arbitrary integer.
We calculate this with an integer type with the same bit width as rawtype T (i.e. 32-bit) in order to use registers efficiently. Since the bit width is finite, the choice of k has a trade-off between the accuracy and the overflow tolerance. I prioritized the accuracy and chose f + k == 24, i.e. the number of bits of mantissa.
Thus,

x * rawone == x * 2^24 * 2^(f - 24)  -  x * 2^24 * 2^-24
           == r * 2^(f - 24) - r * 2^-24

where r == x * 2^24. (0x18 == 24)

    r = _unsafe_trunc(UInt32, round(x * @f32(0x1p24)))
    reinterpret(N, r << UInt8(f - 24) - unsigned(signed(r) >> 0x18))

Do I make sense? And what and how much should I comment?

[kimikage]

BTW, we may want to rename T to N. A month ago I seemed to follow the existing rems.
Edit: Done.

[timholy]

That's a great explanation. I'd say just insert a comment linking to that post and call it good.

[codecov-io]

Codecov Report

Merging #166 into master will increase coverage by 0.19%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master     #166      +/-   ##
==========================================
+ Coverage   87.93%   88.12%   +0.19%     
==========================================
  Files           5        5              
  Lines         373      379       +6     
==========================================
+ Hits          328      334       +6     
  Misses         45       45

Impacted Files	Coverage Δ
src/normed.jl	89.83% <100%> (+0.35%)	⬆️

---

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 930f19a...5fa5e8e. Read the comment docs.