`BigFloat` incorrectly rounded to `Float16` (subnormals)

[nsajko]

julia_rounding_experiment.jl

distance(x, y) = abs(x - y)

function rounding_experiment(bf::BigFloat, r::AbstractFloat)
  p = prevfloat(r)
  n = nextfloat(r)

  r_dist = distance(bf, r)
  p_dist = distance(bf, p)
  n_dist = distance(bf, n)

  if !((r_dist ≤ p_dist) & (r_dist ≤ n_dist))
    print("rounding $r not correct: ")
    if p_dist ≤ n_dist
      println("$p is the correct rounding")
    else
      println("$n is the correct rounding")
    end
  end

  nothing
end

f16_rounding_experiment(bf::BigFloat) = rounding_experiment(bf, Float16(bf))

setprecision(BigFloat, 500)

f16_rounding_experiment(big"1.4901162082026128889687591176485489397376143775948511e-07")

julia> include("/home/nsajko/julia_rounding_experiment.jl")
rounding 1.0e-7 not correct: 2.0e-7 is the correct rounding

julia> versioninfo()
Julia Version 1.11.0-DEV.139
Commit bf00ff4a110 (2023-07-21 19:34 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD Ryzen 3 5300U with Radeon Graphics
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
  Threads: 1 on 8 virtual cores

[oscardssmith]

I think this is a dup of something. the problem is we convert to float64 first. I've meant to fix this for a while but haven't bothered yet

[nsajko]

Yeah, if the BigFloat is converted to Float64 before converting it to Float16, the result is the same.

Is there an easy workaround, even if not performant? I'm trying to revisit the exact rational conversion PR, and the many broken tests are bothering me.

[nsajko]

Yeah, if the BigFloat is converted to Float64 before converting it to Float16, the result is the same.

Wait a minute, this actually means the bug is not just caused by the conversion to Float64?

[nsajko]

I mean, f16_rounding_experiment(big(1.490116208202613e-7)) gives the same error report.

[oscardssmith]

The bug is caused by double rounding. The workaround is relatively easy and relatively performant. You just need to read the bits out of the Bigfloat mantisa and do the rounding yourself. It's conceptually pretty easy but somewhat annoying to implement.

[nsajko]

Currently the conversion from BigFloat to Float16 is implemented by first converting to Float32, not Float64. I suppose the bug would manifest less often if we first converted to Float64, without using Float32 at all.

Would that be OK as a quick fix?

[oscardssmith]

I'm pretty sure I made that PR and then reverted it because it was causing problems (but I forget exactly why).

[giordano]

#40315?