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) <[email protected]>

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!(
+        if !is_native_le_float(datatype, FloatFormat::Double) {
-            datatype,
-            Datatype::FloatingPoint {
-                size: 8,
-                byte_order: DatatypeByteOrder::LittleEndian,
-                ..
-            }
-        ) {
             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!(
+        if !is_native_le_float(datatype, FloatFormat::Single) {
-            datatype,
-            Datatype::FloatingPoint {
-                size: 4,
-                byte_order: DatatypeByteOrder::LittleEndian,
-                ..
-            }
-        ) {
             return None;
         }
         if !raw.len().is_multiple_of(4) {
@@ -919,20 +905,19 @@ pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result<Vec<f64>, FormatEr
 
     // Fast path: native-endian f64 — single bulk memcpy
     #[cfg(target_endian = "little")]
-    if matches!(
+    if is_native_le_float(datatype, FloatFormat::Double) {
-        datatype,
-        Datatype::FloatingPoint {
-            size: 8,
-            byte_order: DatatypeByteOrder::LittleEndian,
-            ..
-        }
-    ) {
         return Ok(native_le_to_vec::<f64>(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<f64, FormatError> {
     match dt {
-        Datatype::FloatingPoint { size, .. } => match size {
+        Datatype::FloatingPoint { .. } => Ok(FloatFormat::of(dt)?.decode(bytes, order)),
-            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::FixedPoint {
             size,
             signed,
@@ -1139,25 +1113,11 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result<Vec<f32>, FormatEr
 
     // Fast path: native-endian f32 — single bulk memcpy
     #[cfg(target_endian = "little")]
-    if matches!(
+    if is_native_le_float(datatype, FloatFormat::Single) {
-        datatype,
-        Datatype::FloatingPoint {
-            size: 4,
-            byte_order: DatatypeByteOrder::LittleEndian,
-            ..
-        }
-    ) {
         return Ok(native_le_to_vec::<f32>(raw, count));
     }
-    // Little-endian half precision (numpy float16): widen directly.
+    // Little-endian IEEE half precision (numpy float16): widen directly.
-    if matches!(
+    if is_native_le_float(datatype, FloatFormat::Half) {
-        datatype,
-        Datatype::FloatingPoint {
-            size: 2,
-            byte_order: DatatypeByteOrder::LittleEndian,
-            ..
-        }
-    ) {
         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<Vec<f32>, 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<FloatFormat, FormatError> {
+        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<u8> = [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.

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], "<u2")
+# 8-bit patterns, all 256 of them
+fp8 = np.arange(256, dtype="u1")
+types = {{
+    "bf16_le": (custom(h5py.h5t.IEEE_F32LE, (15, 7, 8, 0, 7), 127, 2), bf16),
+    "bf16_be": (custom(h5py.h5t.IEEE_F32BE, (15, 7, 8, 0, 7), 127, 2), bf16.byteswap()),
+    "e4m3": (custom(h5py.h5t.IEEE_F32LE, (7, 3, 4, 0, 3), 7, 1), fp8),
+    "e5m2": (custom(h5py.h5t.IEEE_F32LE, (7, 2, 5, 0, 2), 15, 1), fp8),
+}}
+with h5py.File("{path}", "w") as f:
+    for name, (t, raw) in types.items():
+        write(f, name, t, raw)
+    vals = np.array([1.5, -2.25, np.inf, -np.inf, 0.0, -0.0, 6e-8, 65504, 1e-40, 1e300])
+    f.create_dataset("f2_be", data=vals.astype(">f2"))
+    f.create_dataset("f4_be", data=vals.astype(">f4"))
+    f.create_dataset("f8_le", data=vals.astype("<f8"))
+with h5py.File("{path}", "r") as f:
+    for name in list(types) + ["f2_be", "f4_be", "f8_le"]:
+        v = f[name].astype("<f8")[()]
+        # NaN payloads differ between converters; compare NaN as one value.
+        v[np.isnan(v)] = np.nan
+        emit(name, v.view("<u8"))
+"#,
+        path = path.display()
+    );
+    let expected = run_python(&script);
+
+    let canonical = |v: f64| if v.is_nan() { f64::NAN } else { v };
+    let file = File::open(&path).unwrap();
+    for name in [
+        "bf16_le", "bf16_be", "e4m3", "e5m2", "f2_be", "f4_be", "f8_le",
+    ] {
+        let ds = file.dataset(name).unwrap();
+        let want: Vec<u64> = parse(&expected[name]);
+        let got: Vec<u64> = 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<u32> = ds
+            .read_f32()
+            .unwrap()
+            .into_iter()
+            .map(|v| if v.is_nan() { f32::NAN } else { v }.to_bits())
+            .collect();
+        let want32: Vec<u32> = 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");
+    }
+}