Improve accuracy and speed of Float16 conversions from Normed{UInt16} and Normed{UInt8}

src/normed.jl

@@ -142,10 +142,24 @@ function (::Type{T})(x::Normed) where {T <: AbstractFloat}
     convert(T, y)  # needed for types like Float16 which promote arithmetic to Float32
 end
 
-function Base.Float16(x::Normed{Ti,f}) where {Ti <: Union{UInt8, UInt16, UInt32}, f}
-    f == 1 ? Float16(x.i) : Float16(Float32(x))
-end
-function Base.Float16(x::Normed{Ti,f}) where {Ti <: Union{UInt64, UInt128}, f}
+function Base.Float16(x::Normed{T,f}) where {T <: Union{UInt8, UInt32}, f}
+    f == 1 && return Float16(x.i)
+    T === UInt8 && return Float16(Float32(x.i) * Float32(inv_rawone(x)))
+    Float16(Float32(x))
+end
+function Base.Float16(x::Normed{UInt16,f}) where f
+    f == 1 && return Float16(x.i)
+    f16(int::Integer, sf32::Float32) = Float16(Float32(int) * sf32)
+    f ==  3 && return Float16(Float32(x))
+    f ==  4 && return Float16(Float32(x))
+    f == 12 && return f16(signed(widemul(x.i, 0x1001) >> 0x4 | (x.i > 0x7)), @f32(0x1p-20))
+    f == 13 && return f16(signed(widemul(x.i, 0x2001) >> 0x5 | (x.i > 0x7)), @f32(0x1p-21))
+    f == 14 && return f16(signed(widemul(x.i, 0x4001) >> 0x6 | (x.i > 0x7)), @f32(0x1p-22))
+    f == 15 && return f16(x.i - (x.i in (0x7ff7, 0xffee)) + (x.i == 0x800f), 1 / 32767f0)
+    f == 16 && return f16(x.i - (x.i in (0xffaf, 0xffef)), 1 / 65535f0)
+    f16(x.i, Float32(inv_rawone(x)))
+end
+function Base.Float16(x::Normed{T,f}) where {T <: Union{UInt64, UInt128}, f}
     f == 1 ? Float16(x.i) : Float16(Float64(x))
 end
 

test/common.jl

@@ -4,6 +4,35 @@ using Base.Checked
 
 SP = VERSION >= v"1.6.0-DEV.771" ? " " : "" # JuliaLang/julia #37085
 
+# issue #246
+from_bf(::Type{Tf}, x::BigFloat) where Tf= Tf(x)
+function from_bf(::Type{Float16}, x::BigFloat) # `x` is assumed to be finite.
+    sbit = signbit(x) ? 0x8000 : 0x0000
+    iszero(x) && return reinterpret(Float16, sbit)
+    num, pow, _ = Base.decompose(x) # x == num * 2^pow
+    n = abs(num)
+    e = ndigits(n, base=2) - 1 + pow # x == 1.xyz... * 2^e
+    e >= 16 && return reinterpret(Float16, sbit | 0x7c00) # inf
+    e < -25 && return reinterpret(Float16, sbit) # less than half of the min subnormal
+    q = max(e - 10, -24) # exponent of the least significant bit
+    s = pow - q
+    if s >= 0
+        r = n << s
+    else
+        sh = -s
+        r = n >> sh
+        rem = n & ((big(1) << sh) - 1)
+        half = big(1) << (sh - 1)
+        if rem > half || (rem == half && isodd(r))
+            r += 1
+        end
+    end
+    # a carry in `r` naturally propagates into the exponent field
+    bits = ((q + 0x18) << 0xa) + Int(r)
+    bits >= 0x7c00 && return reinterpret(Float16, sbit | 0x7c00) # inf
+    return reinterpret(Float16, sbit | UInt16(bits))
+end
+
 """
     target(X::Type, Ss...; ex = :default)
 

test/normed.jl

@@ -258,7 +258,7 @@ end
             float_err = zero(Tf)
             for i = typemin(T):typemax(T)
                 b_expected = i / BigFloat(FixedPointNumbers.rawone(N))
-                f_expected = Tf(promote_type(Tf, Float32)(b_expected)) # workaround for issue #246
+                f_expected = from_bf(Tf, b_expected)
                 isinf(f_expected) && break # for Float16(::Normed{UInt16,1})
                 f_actual = Tf(reinterpret(N, i))
                 float_err += abs(f_actual - f_expected)
@@ -270,7 +270,7 @@ end
             error_count = 0
             for i in vcat(T(0x00):T(0xFF), (typemax(T)-0xFF):typemax(T))
                 b_expected = i / BigFloat(FixedPointNumbers.rawone(N))
-                f_expected = Tf(promote_type(Tf, Float32)(b_expected)) # workaround for issue #246
+                f_expected = from_bf(Tf, b_expected)
                 isinf(f_expected) && break # for Float16() and Float32()
                 f_actual = Tf(reinterpret(N, i))
                 f_actual == f_expected && continue