Fix for K=64

src/gfloat/round.py

@@ -62,16 +62,15 @@ def round_float(
         # Extract exponent
         expval = int(np.floor(np.log2(vpos)))
 
-        assert expval > -1024 + p  # not yet tested for float64 near-subnormals
-
         # Effective precision, accounting for right shift for subnormal values
         if fi.has_subnormals:
             expval = max(expval, 1 - bias)
 
         # Lift to "integer * 2^e"
         expval = expval - p + 1
 
-        fsignificand = vpos * 2.0**-expval
+        # use ldexp instead of vpos*2**-expval to avoid overflow
+        fsignificand = math.ldexp(vpos, -expval)
 
         # Round
         isignificand = math.floor(fsignificand)