@@ -298,6 +298,30 @@ fn converter(mut pn PrepNumber) u64 {
298298 s0 = q0
299299 }
300300 }
301+
302+ // Handle subnormal (denormalized) numbers - very small numbers near zero
303+ //
304+ // Normal floats have an implicit leading 1 bit in their mantissa (like 1.xxxxx).
305+ // When numbers get too small (binexp < -1022), we can't represent them normally.
306+ // Instead, we use subnormals: set exponent to 0 and shift the mantissa right,
307+ // losing precision gradually. This prevents abrupt underflow to zero.
308+ //
309+ // Example: 1.23e-308 is smaller than the minimum normal float, so we:
310+ // 1. Keep the normalized mantissa from s2 and s1
311+ // 2. Shift it right to "denormalize" it (the leading 1 moves into the mantissa)
312+ // 3. Round correctly using the bits that were shifted out
313+ // 4. Return with exponent = 0 (subnormal marker)
314+ if binexp < - 1022 && (s2 | s1 ) != 0 {
315+ shift := - 1022 - binexp
316+ if shift > 60 {
317+ return if pn.negative { double_minus_zero } else { double_plus_zero }
318+ }
319+ shifted := ((u64 (s2 ) << 32 ) | u64 (s1 )) >> u32 (shift)
320+ q := (shifted >> 8 ) +
321+ u64 ((shifted >> 7 ) & 1 != 0 && ((shifted & 0x7F ) != 0 || (shifted >> 8 ) & 1 != 0 ))
322+ return (q & 0x000FFFFFFFFFFFFF ) | (u64 (pn.negative) << 63 )
323+ }
324+
301325 // rounding if needed
302326 /*
303327 * "round half to even" algorithm
@@ -374,6 +398,8 @@ fn converter(mut pn PrepNumber) u64 {
374398 result = double_plus_infinity
375399 }
376400 } else if binexp < 1 {
401+ // Should not reach here for subnormals anymore (handled earlier)
402+ // This is now only for true zeros
377403 if pn.negative {
378404 result = double_minus_zero
379405 } else {
0 commit comments