fix(format): convert float data read as integers instead of returning bit patterns

read_i32/read_i64/read_u64 on a floating-point dataset reinterpreted the
IEEE bits (1.5 read as i64 was 4609434218613702656). Convert like
libhdf5's hard conversions instead: truncate toward zero and saturate at
the target range; NaN reads as 0.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>

crates/clawhdf5-format/src/data_read.rs

@@ -982,7 +982,74 @@ fn convert_to_f64(
     }
 }
 
+/// One numeric element as stored, before conversion to the caller's type.
+#[derive(Debug, Clone, Copy, PartialEq)]
+enum Scalar {
+    Signed(i64),
+    Unsigned(u64),
+    Float(f64),
+}
+
+impl Scalar {
+    /// Float to integer conversions follow libhdf5's hard conversions:
+    /// truncate toward zero, and saturate a value outside the target range to
+    /// its minimum or maximum. NaN converts to 0 (libhdf5 leaves that case to
+    /// the C cast, whose result is platform-dependent).
+    fn to_i64(self) -> i64 {
+        match self {
+            Scalar::Signed(v) => v,
+            Scalar::Unsigned(v) => v as i64,
+            Scalar::Float(v) => v as i64,
+        }
+    }
+
+    fn to_u64(self) -> u64 {
+        match self {
+            Scalar::Signed(v) => v as u64,
+            Scalar::Unsigned(v) => v,
+            Scalar::Float(v) => v as u64,
+        }
+    }
+
+    fn to_i32(self) -> i32 {
+        match self {
+            Scalar::Signed(v) => v as i32,
+            Scalar::Unsigned(v) => v as i32,
+            Scalar::Float(v) => v as i32,
+        }
+    }
+}
+
+/// Decode one element of a numeric datatype.
+fn decode_scalar(
+    bytes: &[u8],
+    dt: &Datatype,
+    order: &DatatypeByteOrder,
+) -> Result<Scalar, FormatError> {
+    match dt {
+        Datatype::FixedPoint {
+            size,
+            signed,
+            bit_offset,
+            bit_precision,
+            ..
+        } => {
+            let full = read_unsigned_int(bytes, *size as usize, order);
+            let (off, prec) = effective_bits(*size as usize, *bit_offset, *bit_precision);
+            Ok(if *signed {
+                Scalar::Signed(extract_signed(full, off, prec))
+            } else {
+                Scalar::Unsigned(extract_unsigned(full, off, prec))
+            })
+        }
+        _ => convert_to_f64(bytes, dt, order).map(Scalar::Float),
+    }
+}
+
 /// Convert raw bytes to `i64` values.
+///
+/// Floating-point data is converted the way libhdf5 converts it: truncated
+/// toward zero, saturating at the target type's range, with NaN read as 0.
 pub fn read_as_i64(raw: &[u8], datatype: &Datatype) -> Result<Vec<i64>, FormatError> {
     if let Datatype::Array { base_type, .. } = datatype {
         return read_as_i64(raw, base_type);
@@ -1014,17 +1081,18 @@ pub fn read_as_i64(raw: &[u8], datatype: &Datatype) -> Result<Vec<i64>, FormatEr
     }
 
     let order = get_byte_order(datatype);
-    let (off, prec) = fixed_bits(datatype);
     let mut result = Vec::with_capacity(count);
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
-        let full = read_unsigned_int(chunk, elem_size, &order);
+        result.push(decode_scalar(chunk, datatype, &order)?.to_i64());
-        result.push(extract_signed(full, off, prec));
     }
     Ok(result)
 }
 
 /// Convert raw bytes to `u64` values.
+///
+/// Floating-point data is converted the way libhdf5 converts it: truncated
+/// toward zero, saturating at the target type's range, with NaN read as 0.
 pub fn read_as_u64(raw: &[u8], datatype: &Datatype) -> Result<Vec<u64>, FormatError> {
     if let Datatype::Array { base_type, .. } = datatype {
         return read_as_u64(raw, base_type);
@@ -1039,12 +1107,10 @@ pub fn read_as_u64(raw: &[u8], datatype: &Datatype) -> Result<Vec<u64>, FormatEr
     }
     let count = raw.len() / elem_size;
     let order = get_byte_order(datatype);
-    let (off, prec) = fixed_bits(datatype);
     let mut result = Vec::with_capacity(count);
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
-        let full = read_unsigned_int(chunk, elem_size, &order);
+        result.push(decode_scalar(chunk, datatype, &order)?.to_u64());
-        result.push(extract_unsigned(full, off, prec));
     }
     Ok(result)
 }
@@ -1140,6 +1206,9 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result<Vec<f32>, FormatEr
 }
 
 /// Convert raw bytes to `i32` values.
+///
+/// Floating-point data is converted the way libhdf5 converts it: truncated
+/// toward zero, saturating at the target type's range, with NaN read as 0.
 pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result<Vec<i32>, FormatError> {
     if let Datatype::Array { base_type, .. } = datatype {
         return read_as_i32(raw, base_type);
@@ -1170,12 +1239,10 @@ pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result<Vec<i32>, FormatEr
     }
 
     let order = get_byte_order(datatype);
-    let (off, prec) = fixed_bits(datatype);
     let mut result = Vec::with_capacity(count);
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
-        let full = read_unsigned_int(chunk, elem_size, &order);
+        result.push(decode_scalar(chunk, datatype, &order)?.to_i32());
-        result.push(extract_signed(full, off, prec) as i32);
     }
     Ok(result)
 }
@@ -1666,20 +1733,6 @@ fn effective_bits(size: usize, bit_offset: u16, bit_precision: u16) -> (u32, u32
     (bit_offset as u32, prec)
 }
 
-/// `(bit_offset, bit_precision)` for a fixed-point datatype, full width for
-/// other types.
-fn fixed_bits(datatype: &Datatype) -> (u32, u32) {
-    match datatype {
-        Datatype::FixedPoint {
-            size,
-            bit_offset,
-            bit_precision,
-            ..
-        } => effective_bits(*size as usize, *bit_offset, *bit_precision),
-        _ => (0, 0),
-    }
-}
-
 /// Whether a datatype occupies its full storage width (bit offset 0, precision
 /// == size·8), in which case the bulk-copy fast read paths apply. Non
 /// fixed-point types are treated as full width.
@@ -1892,6 +1945,27 @@ mod tests {
         assert_eq!(read_as_u64(&raw, &dt).unwrap(), vec![4095, 1, 2048]);
     }
 
+    #[test]
+    fn float_to_int_truncates_and_saturates() {
+        // Values libhdf5 hands to an undefined C cast: NaN reads as 0 and
+        // exactly 2^63 saturates instead of wrapping to i64::MIN.
+        let dt = make_f64_le_type();
+        let vals = [f64::NAN, 2f64.powi(63), -2.5, 2.0f64.powi(64)];
+        let raw: Vec<u8> = vals.iter().flat_map(|v| v.to_le_bytes()).collect();
+        assert_eq!(
+            read_as_i64(&raw, &dt).unwrap(),
+            vec![0, i64::MAX, -2, i64::MAX]
+        );
+        assert_eq!(
+            read_as_u64(&raw, &dt).unwrap(),
+            vec![0, 1 << 63, 0, u64::MAX]
+        );
+        assert_eq!(
+            read_as_i32(&raw, &dt).unwrap(),
+            vec![0, i32::MAX, -2, i32::MAX]
+        );
+    }
+
     #[test]
     fn full_width_signed_unchanged() {
         // Regression: full-width 32-bit signed must be unaffected.

crates/clawhdf5/tests/numeric_conversion_interop.rs

@@ -0,0 +1,129 @@
+//! Numeric conversions on read, checked against h5py/libhdf5.
+//!
+//! h5py writes each file and prints what libhdf5 converts the data to
+//! (`Dataset.astype`); the typed readers must return the same values.
+//! Skipped when python3 with h5py is unavailable, unless
+//! `CLAWHDF5_REQUIRE_INTEROP=1`.
+
+use std::collections::HashMap;
+use std::process::Command;
+
+use clawhdf5::File;
+
+fn python() -> String {
+    std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string())
+}
+
+fn interop_required() -> bool {
+    std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1")
+}
+
+fn python_available() -> bool {
+    Command::new(python())
+        .args(["-c", "import h5py, numpy"])
+        .output()
+        .map(|o| o.status.success())
+        .unwrap_or(false)
+}
+
+macro_rules! skip_if_no_python {
+    () => {
+        if !python_available() {
+            assert!(
+                !interop_required(),
+                "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with h5py is not available"
+            );
+            eprintln!("SKIP: python3 with h5py not available");
+            return;
+        }
+    };
+}
+
+/// Run `script` (which writes the file at `path`) and return its stdout as
+/// `key -> values`, one `key v1 v2 ...` line per key.
+fn run_python(script: &str) -> HashMap<String, Vec<String>> {
+    let output = Command::new(python())
+        .args(["-c", script])
+        .output()
+        .expect("failed to run python");
+    assert!(
+        output.status.success(),
+        "python failed:\n{}",
+        String::from_utf8_lossy(&output.stderr)
+    );
+    String::from_utf8_lossy(&output.stdout)
+        .lines()
+        .filter_map(|line| {
+            let mut words = line.split_whitespace().map(str::to_string);
+            Some((words.next()?, words.collect()))
+        })
+        .collect()
+}
+
+fn parse<T: std::str::FromStr>(values: &[String]) -> Vec<T>
+where
+    T::Err: std::fmt::Debug,
+{
+    values.iter().map(|v| v.parse().unwrap()).collect()
+}
+
+/// Python prelude: `emit(key, array)` prints one line of integers.
+const PRELUDE: &str = r#"
+import h5py, numpy as np
+def emit(key, arr):
+    print(key, *[int(v) for v in np.asarray(arr).ravel()])
+"#;
+
+#[test]
+fn float_dataset_read_as_integers_converts_like_libhdf5() {
+    // read_i32/read_i64/read_u64 on a float dataset used to return the raw
+    // IEEE bit patterns (1.5 read as i64 was 4609434218613702656).
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("float_to_int.h5");
+    let script = format!(
+        r#"{PRELUDE}
+vals = [1.5, -2.75, 3e9, 1e300, -1e300, -0.5, 0.0, 7.99, np.inf, -np.inf, 1e19, -1e19]
+with h5py.File("{path}", "w") as f:
+    for name, dt in (("f8", "<f8"), ("f8be", ">f8"), ("f4", "<f4")):
+        f.create_dataset(name, data=np.array(vals).astype(dt))
+    # libhdf5's half conversions are not saturating (an infinite half becomes
+    # INT_MIN whatever its sign, a negative one wraps as u64), so the half
+    # case stays finite and its u64 read is checked separately below.
+    f.create_dataset("f2", data=np.array([1.5, -2.75, -0.5, 0.0, 7.99, 65504, -65504], "<f2"))
+with h5py.File("{path}", "r") as f:
+    for name in ("f8", "f8be", "f4", "f2"):
+        d = f[name]
+        emit(name + ":i32", d.astype("<i4")[()])
+        emit(name + ":i64", d.astype("<i8")[()])
+        emit(name + ":u64", d.astype("<u8")[()])
+"#,
+        path = path.display()
+    );
+    let expected = run_python(&script);
+
+    let file = File::open(&path).unwrap();
+    for name in ["f8", "f8be", "f4", "f2"] {
+        let ds = file.dataset(name).unwrap();
+        assert_eq!(
+            ds.read_i32().unwrap(),
+            parse::<i32>(&expected[&format!("{name}:i32")]),
+            "{name} as i32"
+        );
+        assert_eq!(
+            ds.read_i64().unwrap(),
+            parse::<i64>(&expected[&format!("{name}:i64")]),
+            "{name} as i64"
+        );
+        if name != "f2" {
+            assert_eq!(
+                ds.read_u64().unwrap(),
+                parse::<u64>(&expected[&format!("{name}:u64")]),
+                "{name} as u64"
+            );
+        }
+    }
+    // Negative values saturate at 0 rather than wrapping.
+    let f2 = file.dataset("f2").unwrap();
+    assert_eq!(f2.read_u64().unwrap(), vec![1, 0, 0, 0, 7, 65504, 0]);
+}