Float8_4 subnormal bug?

[bjarthur]

the README says there is one. is that still the case?

[milankl]

I was wondering that too. Maybe there are some tests that are skipped? Seems I haven’t been good at documenting this. Sorry. But I doubt it has been resolved.

[milankl]

So there are only 8 subnormals in Float8_4, 0000_0001 to 0000_0111

julia> for i in 0b0000_0001:0b0000_0111
           f = Float8_4(i)
           println(bitstring(f)," ",f)
       end
00000001 Float8_4(0.001953125)
00000010 Float8_4(0.00390625)
00000011 Float8_4(0.005859375)
00000100 Float8_4(0.0078125)
00000101 Float8_4(0.009765625)
00000110 Float8_4(0.01171875)
00000111 Float8_4(0.013671875)

Just found this

julia> x = rand(Float32,100)*floatmin(Float8_4);
julia> sort!(x)

julia> cat(x,Float32.(Float8_4.(x)),bitstring.(Float8_4.(x)),dims=2)
100×3 Matrix{Any}:
# true value, rounded to Float8_4, as bitstring
 0.000104143  0.0         "00000000"
 0.000186678  0.0         "00000000"
 0.000538887  0.00390625  "00000010"    # should be 0x00
 0.000785412  0.00390625  "00000010"    # should be 0x00
 0.00083075   0.00585938  "00000011"    # should be 0x00
 0.00103349   0.00390625  "00000010"    # should be 0x01
 0.00112769   0.0117188   "00000110"    # should be 0x01
 0.00133011   0.00195312  "00000001"
 0.00156096   0.00195312  "00000001"
 0.00160453   0.0078125   "00000100"    # should be 0x01
 0.00234231   0.00195312  "00000001"
 0.00237345   0.00195312  "00000001"
 0.00257766   0.00195312  "00000001"
 0.00260994   0.00195312  "00000001"
 0.00262435   0.00195312  "00000001"
 0.00278183   0.00195312  "00000001"
 0.00281697   0.00195312  "00000001"
 0.00293266   0.00390625  "00000010"
 0.00296257   0.00390625  "00000010"
 0.00329266   0.00390625  "00000010"
 0.0034002    0.00390625  "00000010"
 0.00341041   0.00390625  "00000010"
 0.00370326   0.00390625  "00000010"
 0.00383257   0.00390625  "00000010"
 0.00384757   0.00390625  "00000010"
 0.00418744   0.00390625  "00000010"
 0.0042767    0.00390625  "00000010"
 0.00457281   0.00390625  "00000010"
 0.00461949   0.00390625  "00000010"
 0.00474257   0.00390625  "00000010"
 0.00552679   0.00585938  "00000011"
 0.00565881   0.00585938  "00000011"
 0.00572785   0.00585938  "00000011"
 0.00575996   0.00585938  "00000011"
 0.00585618   0.00585938  "00000011"
 0.00596875   0.00585938  "00000011"
 0.00601179   0.00585938  "00000011"
 0.00604342   0.00585938  "00000011"
 0.00604443   0.00585938  "00000011"
 0.00604547   0.00585938  "00000011"
 0.0060702    0.00585938  "00000011"
 0.00613414   0.00585938  "00000011"
 0.00614859   0.00585938  "00000011"
 0.00630236   0.00585938  "00000011"
 0.00631475   0.00585938  "00000011"
 0.0064428    0.00585938  "00000011"
 0.00684184   0.0078125   "00000100"
 0.00713566   0.0078125   "00000100"
 0.00775759   0.0078125   "00000100"
 0.00847213   0.0078125   "00000100"
 0.00860756   0.0078125   "00000100"
 0.00868505   0.0078125   "00000100"
 0.00870331   0.0078125   "00000100"
 0.00883958   0.00976562  "00000101"
 0.00908779   0.00976562  "00000101"
 0.00909926   0.00976562  "00000101"
 0.00919047   0.00976562  "00000101"
 0.00943508   0.00976562  "00000101"
 0.00948936   0.00976562  "00000101"
 0.00963979   0.00976562  "00000101"
 0.00989018   0.00976562  "00000101"
 0.0101677    0.00976562  "00000101"
 0.0102963    0.00976562  "00000101"
 0.0103968    0.00976562  "00000101"
 0.0104353    0.00976562  "00000101"
 0.0104995    0.00976562  "00000101"
 0.01055      0.00976562  "00000101"
 0.0105936    0.00976562  "00000101"
 0.0107531    0.0117188   "00000110"
 0.0108297    0.0117188   "00000110"
 0.0108384    0.0117188   "00000110"
 0.0108987    0.0117188   "00000110"
 0.0109998    0.0117188   "00000110"
 0.0111872    0.0117188   "00000110"
 0.0112731    0.0117188   "00000110"
 0.0113971    0.0117188   "00000110"
 0.0114834    0.0117188   "00000110"
 0.0120236    0.0117188   "00000110"
 0.0120299    0.0117188   "00000110"
 0.0120968    0.0117188   "00000110"
 0.0121795    0.0117188   "00000110"
 0.0121866    0.0117188   "00000110"
 0.0128721    0.0136719   "00000111"
 0.0130633    0.0136719   "00000111"
 0.0132557    0.0136719   "00000111"
 0.013483     0.0136719   "00000111"
 0.0136537    0.0136719   "00000111"
 0.0136634    0.0136719   "00000111"
 0.0137801    0.0136719   "00000111"
 0.0138396    0.0136719   "00000111"
 0.0138615    0.0136719   "00000111"
 0.0139297    0.0136719   "00000111"
 0.0139765    0.0136719   "00000111"
 0.0140286    0.0136719   "00000111"
 0.0141437    0.0136719   "00000111"
 0.0146549    0.015625    "00001000"
 0.0150794    0.015625    "00001000"
 0.0152748    0.015625    "00001000"
 0.0153345    0.015625    "00001000"
 0.0156243    0.015625    "00001000"

the rest looks fine to me.

Update: Just repeated this for Float8 which looks good though!

[bjarthur]

IIUC it's the conversion of Float32 to Float8_4 that is the problem. this method lives in a file ending in "old.jl". perhaps it's time to update Float8_4(::Float32) to use @JeffreySarnoff 's lookup table algorithm instead? fix the bug and improve performance in one fell swoop.

[AntonOresten]

The subnormals bug seems to come from using a small constant shift in the "very small -> zero" block, which leaves residual mantissa bits and can round up to a spurious subnormal even when |x| << min_subnormal.

I tried replacing the fixed shift in the e < e_limits[1] branch with a shift that depends on how far below the subnormal threshold the Float32 exponent is:

# before
if e < e_limits[1]                      # Very small numbers map to +- zero
    basetable[i|0x000+1] = zero(T)
    basetable[i|0x100+1] = -zero(T)
    shifttable[i|0x000+1] = n_significant_bits(T)+1
    shifttable[i|0x100+1] = n_significant_bits(T)+1

# after
if e < e_limits[1]                      # Very small numbers map round to zero or smallest subnormal
    basetable[i|0x000+1] = zero(T)
    basetable[i|0x100+1] = -zero(T)
    shifttable[i|0x000+1] = -e+e_shift_subnorm
    shifttable[i|0x100+1] = -e+e_shift_subnorm

EDIT: actually I just realized this is the same as the e < e_limits[2] branch.

Comparison before and after

julia> Random.seed!(0); x = sort!(rand(Float32,100)*floatmin(Float8_4));

julia> a = hcat(x, Float32.(Float8_4.(x)), bitstring.(Float8_4.(x)))
100×3 Matrix{Any}:
 0.000400645  0.00390625  "00000010"
 0.000463245  0.00390625  "00000010"
 0.000821631  0.0078125   "00000100"
 0.000915817  0.0136719   "00000111"
 0.00135854   0.00585938  "00000011"
 0.00151019   0.015625    "00001000"
 0.00155841   0.0         "00000000"
 0.00176078   0.00195312  "00000001"
 0.00179476   0.0         "00000000"
 0.00183332   0.0         "00000000"
 0.00183718   0.00390625  "00000010"
 0.00203049   0.00195312  "00000001"
 0.00274331   0.00195312  "00000001"
 0.00278716   0.00195312  "00000001"
 0.00291116   0.00195312  "00000001"
 0.00307447   0.00390625  "00000010"
 0.00332867   0.00390625  "00000010"
 0.00339464   0.00390625  "00000010"
 ⋮                        
 0.0130657    0.0136719   "00000111"
 0.0135685    0.0136719   "00000111"
 0.0135907    0.0136719   "00000111"
 0.0136306    0.0136719   "00000111"
 0.0136485    0.0136719   "00000111"
 0.0140136    0.0136719   "00000111"
 0.0140955    0.0136719   "00000111"
 0.0141938    0.0136719   "00000111"
 0.0142334    0.0136719   "00000111"
 0.0144193    0.0136719   "00000111"
 0.0144471    0.0136719   "00000111"
 0.0146861    0.015625    "00001000"
 0.0146966    0.015625    "00001000"
 0.0147704    0.015625    "00001000"
 0.0148779    0.015625    "00001000"
 0.0152605    0.015625    "00001000"
 0.0153877    0.015625    "00001000"

julia> Random.seed!(0); x = sort!(rand(Float32,100)*floatmin(Float8_4));

julia> b = hcat(x, Float32.(Float8_4.(x)), bitstring.(Float8_4.(x)))
100×3 Matrix{Any}:
 0.000400645  0.0         "00000000"
 0.000463245  0.0         "00000000"
 0.000821631  0.0         "00000000"
 0.000915817  0.0         "00000000"
 0.00135854   0.00195312  "00000001"
 0.00151019   0.00195312  "00000001"
 0.00155841   0.00195312  "00000001"
 0.00176078   0.00195312  "00000001"
 0.00179476   0.00195312  "00000001"
 0.00183332   0.00195312  "00000001"
 0.00183718   0.00195312  "00000001"
 0.00203049   0.00195312  "00000001"
 0.00274331   0.00195312  "00000001"
 0.00278716   0.00195312  "00000001"
 0.00291116   0.00195312  "00000001"
 0.00307447   0.00390625  "00000010"
 0.00332867   0.00390625  "00000010"
 0.00339464   0.00390625  "00000010"
 ⋮                        
 0.0130657    0.0136719   "00000111"
 0.0135685    0.0136719   "00000111"
 0.0135907    0.0136719   "00000111"
 0.0136306    0.0136719   "00000111"
 0.0136485    0.0136719   "00000111"
 0.0140136    0.0136719   "00000111"
 0.0140955    0.0136719   "00000111"
 0.0141938    0.0136719   "00000111"
 0.0142334    0.0136719   "00000111"
 0.0144193    0.0136719   "00000111"
 0.0144471    0.0136719   "00000111"
 0.0146861    0.015625    "00001000"
 0.0146966    0.015625    "00001000"
 0.0147704    0.015625    "00001000"
 0.0148779    0.015625    "00001000"
 0.0152605    0.015625    "00001000"
 0.0153877    0.015625    "00001000"

julia> count(a[:,3] .!== b[:,3])
10

10 mismatches because 0.00176078 0.00195312 "00000001" was correct by chance.

The new behavior appears to be rounding correctly.

I can open a PR if this is satisfactory. I should also note that I've been tinkering with Microfloats.jl.