From 417c9516ca69e9908c10d914b5adc4844a8ff299 Mon Sep 17 00:00:00 2001 From: osobh Date: Fri, 25 Sep 2026 21:06:19 -0500 Subject: [PATCH] fix(format): decode floats by their datatype fields, not their size Every 2-byte float was decoded as IEEE half, so bfloat16 (HDF5 2.0's H5T_FLOAT_BFLOAT16*, or any custom 8-bit-exponent type) read wrong: 1.5 as 1.9375, +inf as NaN. 1-byte FP8 floats were refused. Read the exponent/mantissa location and size and the bias from the datatype message: IEEE half/single/double keep their existing paths (half still through clawhdf5_format::float16), any other IEEE-style layout up to 64 bits whose values fit f64 (bfloat16, FP8 E4M3/E5M2, ...) is decoded generically, and the bulk-copy and zero-copy fast paths now require the IEEE layout rather than just the size. Datatypes with fields that describe no float still fall back to IEEE by size; layouts that cannot be represented in f64 (x87 80-bit, binary128) remain an error. Co-Authored-By: Claude Opus 5.5 (1M context) --- crates/clawhdf5-format/src/data_read.rs | 298 ++++++++++++++---- .../tests/numeric_conversion_interop.rs | 85 +++++ 2 files changed, 320 insertions(+), 63 deletions(-) diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index 01060a6..1309080 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -773,14 +773,7 @@ pub fn read_as_f64_zerocopy<'a>(raw: &'a [u8], datatype: &Datatype) -> Option<&' // Only native LE f64 is eligible #[cfg(target_endian = "little")] { - if !matches!( - datatype, - Datatype::FloatingPoint { - size: 8, - byte_order: DatatypeByteOrder::LittleEndian, - .. - } - ) { + if !is_native_le_float(datatype, FloatFormat::Double) { return None; } if !raw.len().is_multiple_of(8) { @@ -809,14 +802,7 @@ pub fn read_as_f64_zerocopy<'a>(raw: &'a [u8], datatype: &Datatype) -> Option<&' pub fn read_as_f32_zerocopy<'a>(raw: &'a [u8], datatype: &Datatype) -> Option<&'a [f32]> { #[cfg(target_endian = "little")] { - if !matches!( - datatype, - Datatype::FloatingPoint { - size: 4, - byte_order: DatatypeByteOrder::LittleEndian, - .. - } - ) { + if !is_native_le_float(datatype, FloatFormat::Single) { return None; } if !raw.len().is_multiple_of(4) { @@ -919,20 +905,19 @@ pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result, FormatEr // Fast path: native-endian f64 — single bulk memcpy #[cfg(target_endian = "little")] - if matches!( - datatype, - Datatype::FloatingPoint { - size: 8, - byte_order: DatatypeByteOrder::LittleEndian, - .. - } - ) { + if is_native_le_float(datatype, FloatFormat::Double) { return Ok(native_le_to_vec::(raw, count)); } let order = get_byte_order(datatype); let mut result = Vec::with_capacity(count); - + if let Datatype::FloatingPoint { .. } = datatype { + let format = FloatFormat::of(datatype)?; + for chunk in raw.chunks_exact(elem_size) { + result.push(format.decode(chunk, &order)); + } + return Ok(result); + } for i in 0..count { let chunk = &raw[i * elem_size..(i + 1) * elem_size]; let val = convert_to_f64(chunk, datatype, &order)?; @@ -947,18 +932,7 @@ fn convert_to_f64( order: &DatatypeByteOrder, ) -> Result { match dt { - Datatype::FloatingPoint { size, .. } => match size { - 4 => { - let v = read_f32_bytes(bytes, order); - Ok(v as f64) - } - 8 => Ok(read_f64_bytes(bytes, order)), - 2 => Ok(read_f16_bytes(bytes, order) as f64), - _ => Err(FormatError::DataSizeMismatch { - expected: 8, - actual: *size as usize, - }), - }, + Datatype::FloatingPoint { .. } => Ok(FloatFormat::of(dt)?.decode(bytes, order)), Datatype::FixedPoint { size, signed, @@ -1139,25 +1113,11 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result, FormatEr // Fast path: native-endian f32 — single bulk memcpy #[cfg(target_endian = "little")] - if matches!( - datatype, - Datatype::FloatingPoint { - size: 4, - byte_order: DatatypeByteOrder::LittleEndian, - .. - } - ) { + if is_native_le_float(datatype, FloatFormat::Single) { return Ok(native_le_to_vec::(raw, count)); } - // Little-endian half precision (numpy float16): widen directly. - if matches!( - datatype, - Datatype::FloatingPoint { - size: 2, - byte_order: DatatypeByteOrder::LittleEndian, - .. - } - ) { + // Little-endian IEEE half precision (numpy float16): widen directly. + if is_native_le_float(datatype, FloatFormat::Half) { let (halves, _) = raw[..count * 2].as_chunks::<2>(); return Ok(halves .iter() @@ -1167,18 +1127,22 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result, FormatEr let order = get_byte_order(datatype); let mut result = Vec::with_capacity(count); + if let Datatype::FloatingPoint { .. } = datatype { + let format = FloatFormat::of(datatype)?; + for chunk in raw.chunks_exact(elem_size) { + result.push(match format { + FloatFormat::Single => read_f32_bytes(chunk, &order), + FloatFormat::Half => read_f16_bytes(chunk, &order), + // Double rounds; every other supported layout (bfloat16, FP8) + // is exact in f32. + _ => format.decode(chunk, &order) as f32, + }); + } + return Ok(result); + } for i in 0..count { let chunk = &raw[i * elem_size..(i + 1) * elem_size]; match datatype { - Datatype::FloatingPoint { size: 4, .. } => { - result.push(read_f32_bytes(chunk, &order)); - } - Datatype::FloatingPoint { size: 8, .. } => { - result.push(read_f64_bytes(chunk, &order) as f32); - } - Datatype::FloatingPoint { size: 2, .. } => { - result.push(read_f16_bytes(chunk, &order)); - } Datatype::FixedPoint { signed: true, size, @@ -1693,6 +1657,174 @@ fn reorder_bytes(bytes: &[u8], order: &DatatypeByteOrder) -> [u8; 8] { buf } +/// How the bits of a floating-point datatype are laid out, read from the +/// datatype message's fields rather than assumed from its size (a 2-byte +/// float may be IEEE half or bfloat16). +#[derive(Debug, Clone, Copy, PartialEq)] +enum FloatFormat { + /// IEEE-754 binary16. + Half, + /// IEEE-754 binary32. + Single, + /// IEEE-754 binary64. + Double, + /// Any other IEEE-style layout (implied leading mantissa bit, all-ones + /// exponent for infinity/NaN) whose values are all exact in `f64`: + /// bfloat16, the FP8 formats, and similar. + Other(FloatLayout), +} + +#[derive(Debug, Clone, Copy, PartialEq)] +struct FloatLayout { + exponent_location: u32, + exponent_size: u32, + mantissa_location: u32, + mantissa_size: u32, + exponent_bias: u32, +} + +impl FloatFormat { + fn of(dt: &Datatype) -> Result { + let Datatype::FloatingPoint { + size, + exponent_location, + exponent_size, + mantissa_location, + mantissa_size, + exponent_bias, + .. + } = dt + else { + return Err(FormatError::TypeMismatch { + expected: "FloatingPoint", + actual: datatype_name(dt), + }); + }; + let layout = FloatLayout { + exponent_location: u32::from(*exponent_location), + exponent_size: u32::from(*exponent_size), + mantissa_location: u32::from(*mantissa_location), + mantissa_size: u32::from(*mantissa_size), + exponent_bias: *exponent_bias, + }; + let fields = ( + layout.exponent_location, + layout.exponent_size, + layout.mantissa_location, + layout.mantissa_size, + layout.exponent_bias, + ); + let bits = size.saturating_mul(8); + // The sign bit is not kept in `Datatype`; every standard layout has it + // directly above the exponent, with the mantissa below. + let well_formed = layout.exponent_size > 0 + && layout.mantissa_size > 0 + && layout.mantissa_location + layout.mantissa_size <= layout.exponent_location + && layout.exponent_location + layout.exponent_size < bits; + match (size, fields) { + (2, (10, 5, 0, 10, 15)) => Ok(FloatFormat::Half), + (4, (23, 8, 0, 23, 127)) => Ok(FloatFormat::Single), + (8, (52, 11, 0, 52, 1023)) => Ok(FloatFormat::Double), + _ if well_formed + && *size <= 8 + && layout.exponent_size <= 11 + && layout.mantissa_size <= 52 => + { + Ok(FloatFormat::Other(layout)) + } + // Fields that cannot describe any float (e.g. left zeroed by a + // hand-built datatype): fall back to the IEEE type of that size. + (2, _) if !well_formed => Ok(FloatFormat::Half), + (4, _) if !well_formed => Ok(FloatFormat::Single), + (8, _) if !well_formed => Ok(FloatFormat::Double), + // x87 80-bit extended, binary128, ...: not representable in f64. + _ => Err(FormatError::TypeMismatch { + expected: "floating point of at most 64 bits (IEEE-style layout)", + actual: "FloatingPoint", + }), + } + } + + fn decode(self, bytes: &[u8], order: &DatatypeByteOrder) -> f64 { + match self { + FloatFormat::Half => f64::from(read_f16_bytes(bytes, order)), + FloatFormat::Single => f64::from(read_f32_bytes(bytes, order)), + FloatFormat::Double => read_f64_bytes(bytes, order), + FloatFormat::Other(layout) => { + layout.decode(read_unsigned_int(bytes, bytes.len(), order)) + } + } + } +} + +impl FloatLayout { + /// Decode the value held in the low `size * 8` bits of `bits`. + fn decode(self, bits: u64) -> f64 { + let field = |location: u32, size: u32| (bits >> location) & ((1u64 << size) - 1); + let exponent = field(self.exponent_location, self.exponent_size); + let mantissa = field(self.mantissa_location, self.mantissa_size); + let negative = field(self.exponent_location + self.exponent_size, 1) == 1; + let max_exponent = (1u64 << self.exponent_size) - 1; + let magnitude = if exponent == max_exponent { + if mantissa == 0 { + f64::INFINITY + } else { + f64::NAN + } + } else { + let bias = i64::from(self.exponent_bias); + let msize = i64::from(self.mantissa_size); + // value = significand * 2^power, with an implied leading 1 unless + // the number is subnormal (exponent field 0). + let (significand, power) = if exponent == 0 { + (mantissa, 1 - bias - msize) + } else { + ( + mantissa | (1u64 << self.mantissa_size), + exponent as i64 - bias - msize, + ) + }; + scale_by_pow2(significand as f64, power) + }; + if negative { -magnitude } else { magnitude } + } +} + +/// `x * 2^power` without `std` (no `powi`/`libm`). `x` is a non-negative +/// integer below 2^53, so it is exact. +fn scale_by_pow2(x: f64, power: i64) -> f64 { + if x == 0.0 || power < -1200 { + return 0.0; + } + if power > 1100 { + return f64::INFINITY; + } + let pow2 = |p: i64| f64::from_bits(((p + 1023) as u64) << 52); + let mut x = x; + let mut power = power; + while power > 1023 { + x *= pow2(1023); + power -= 1023; + } + while power < -1022 { + x *= pow2(-1022); + power += 1022; + } + x * pow2(power) +} + +/// Whether `datatype` is the little-endian IEEE float `format`, whose bytes +/// can be copied straight into native values on a little-endian target. +fn is_native_le_float(datatype: &Datatype, format: FloatFormat) -> bool { + matches!( + datatype, + Datatype::FloatingPoint { + byte_order: DatatypeByteOrder::LittleEndian, + .. + } + ) && FloatFormat::of(datatype).is_ok_and(|f| f == format) +} + fn read_f64_bytes(bytes: &[u8], order: &DatatypeByteOrder) -> f64 { let buf = reorder_bytes(bytes, order); f64::from_le_bytes(buf) @@ -1976,6 +2108,46 @@ mod tests { ); } + #[test] + fn bfloat16_and_fp8_decode_by_fields() { + // bfloat16 is a 2-byte float that is not IEEE half. + let bf16 = Datatype::FloatingPoint { + size: 2, + byte_order: DatatypeByteOrder::LittleEndian, + bit_offset: 0, + bit_precision: 16, + exponent_location: 7, + exponent_size: 8, + mantissa_location: 0, + mantissa_size: 7, + exponent_bias: 127, + }; + let raw: Vec = [0x3FC0u16, 0xC010, 0x7F80, 0x0001] + .iter() + .flat_map(|v| v.to_le_bytes()) + .collect(); + let got = read_as_f64(&raw, &bf16).unwrap(); + assert_eq!(&got[..3], &[1.5, -2.25, f64::INFINITY]); + assert_eq!(got[3], 2f64.powi(-133)); // smallest subnormal + assert_eq!(read_as_f32(&raw, &bf16).unwrap()[..2], [1.5, -2.25]); + + // FP8 E4M3: 1, -1, 2, 0, NaN (IEEE-style, as libhdf5 treats it). + let e4m3 = Datatype::FloatingPoint { + size: 1, + byte_order: DatatypeByteOrder::LittleEndian, + bit_offset: 0, + bit_precision: 8, + exponent_location: 3, + exponent_size: 4, + mantissa_location: 0, + mantissa_size: 3, + exponent_bias: 7, + }; + let got = read_as_f64(&[0x38, 0xB8, 0x40, 0x00, 0x7E], &e4m3).unwrap(); + assert_eq!(&got[..4], &[1.0, -1.0, 2.0, 0.0]); + assert!(got[4].is_nan()); + } + #[test] fn full_width_signed_unchanged() { // Regression: full-width 32-bit signed must be unaffected. diff --git a/crates/clawhdf5/tests/numeric_conversion_interop.rs b/crates/clawhdf5/tests/numeric_conversion_interop.rs index e75700a..556f2a9 100644 --- a/crates/clawhdf5/tests/numeric_conversion_interop.rs +++ b/crates/clawhdf5/tests/numeric_conversion_interop.rs @@ -188,3 +188,88 @@ with h5py.File("{path}", "r") as f: let i8be = file.dataset("i8be").unwrap(); assert_eq!(i8be.read_u64().unwrap(), vec![(1 << 40) + 5, 0, 7, 0]); } + +#[test] +fn floats_decode_by_their_datatype_fields() { + // Every 2-byte float used to decode as IEEE half, so bfloat16 1.5 read as + // 1.9375 and +inf as NaN; 1-byte FP8 floats were refused. + skip_if_no_python!(); + let dir = tempfile::tempdir().unwrap(); + let path = dir.path().join("float_layouts.h5"); + let script = format!( + r#"{PRELUDE} +def custom(base, fields, bias, size): + # fields: (sign pos, exponent pos, exponent size, mantissa pos, mantissa size) + t = base.copy() + t.set_fields(*fields) + t.set_ebias(bias) + t.set_precision(size * 8) + t.set_size(size) + return t + +def write(f, name, ftype, raw): + raw = np.ascontiguousarray(raw) + space = h5py.h5s.create_simple(raw.shape) + ds = h5py.h5d.create(f.id, name.encode(), ftype, space) + ds.write(h5py.h5s.ALL, h5py.h5s.ALL, raw, mtype=ftype) + +# bfloat16: 1.5, -2.25, +inf, 0, 3.140625, 1, -0, smallest subnormal, +# largest finite, NaN +bf16 = np.array([0x3FC0, 0xC010, 0x7F80, 0x0000, 0x4049, 0x3F80, 0x8000, 0x0001, + 0x7F7F, 0x7FC1], "f2")) + f.create_dataset("f4_be", data=vals.astype(">f4")) + f.create_dataset("f8_le", data=vals.astype(" = parse(&expected[name]); + let got: Vec = ds + .read_f64() + .unwrap() + .into_iter() + .map(|v| canonical(v).to_bits()) + .collect(); + assert_eq!(got, want, "{name} as f64"); + // f32 reads agree too (every value here is exact in f32 except the + // f64 dataset, which rounds like `as f32`). + let got32: Vec = ds + .read_f32() + .unwrap() + .into_iter() + .map(|v| if v.is_nan() { f32::NAN } else { v }.to_bits()) + .collect(); + let want32: Vec = want + .iter() + .map(|&b| canonical(f64::from_bits(b)) as f32) + .map(|v| if v.is_nan() { f32::NAN } else { v }.to_bits()) + .collect(); + assert_eq!(got32, want32, "{name} as f32"); + } +}