From 95dcb0445452dcb936de616945866bacaebbf9f0 Mon Sep 17 00:00:00 2001
From: osobh <om.sobh@gmail.com>
Date: Fri, 25 Sep 2026 21:15:49 -0500
Subject: [PATCH] fix(format): scale-offset float decode with libhdf5's
 arithmetic

D-scale floats were rebuilt as `minval + code / 10^D` in f64 and then
rounded to f32 once, but libhdf5 (H5Z_scaleoffset_modify_3/4 with
`float`/`powf`) computes `(float)(int)code / powf(10, D) + min` in single
precision. The two differ by 1 ULP for some values: le_data.h5
/Scale_offset_float_data_{le,be} gave 1.6663332 (0x3fd54a69) where
libhdf5 gives 1.6663333 (0x3fd54a6a). Use f32 arithmetic for 4-byte
floats and `(double)(long)code / pow(10, D) + min` for 8-byte ones.

Test: scaleoffset_float_dscale_matches_libhdf5_bits (le_data.h5 float
LE/BE and double chunks, bit-exact against h5py); failed before.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
---
 crates/clawhdf5-format/src/filters.rs | 54 ++++++++++++++++++++++++++-
 1 file changed, 53 insertions(+), 1 deletion(-)

diff --git a/crates/clawhdf5-format/src/filters.rs b/crates/clawhdf5-format/src/filters.rs
index 54e56d4..94fa15f 100644
--- a/crates/clawhdf5-format/src/filters.rs
+++ b/crates/clawhdf5-format/src/filters.rs
@@ -246,8 +246,20 @@ fn scaleoffset_decompress(
                     fill_value
                 } else if is_escale {
                     minval + code as f64 * powi_f64(2.0, scale_factor)
+                } else if elem_size == 4 {
+                    // H5Z_scaleoffset_modify_3/4 for `float`: the code is
+                    // read as an `int` and everything is single precision,
+                    // `(float)code / powf(10, D) + min`. Doing it in f64 and
+                    // rounding once at the end is off by 1 ULP at times.
+                    let d = if scale_factor >= 0 {
+                        powi_f64(10.0, scale_factor) as f32
+                    } else {
+                        1.0 / powi_f64(10.0, -scale_factor) as f32
+                    };
+                    ((code as u32 as i32) as f32 / d + minval as f32) as f64
                 } else {
-                    minval + code as f64 / powi_f64(10.0, scale_factor)
+                    // ... and for `double`: `(double)(long)code / pow(10, D) + min`.
+                    (code as i64) as f64 / powi_f64(10.0, scale_factor) + minval
                 }
             })
             .collect();
@@ -1630,6 +1642,46 @@ mod tests {
         }
     }
 
+    /// `le_data.h5` scale-offset (D-scale, D = 3, fill -2.2) chunks, decoded
+    /// bit for bit as libhdf5 does: single-precision arithmetic for `float`
+    /// (we computed in f64 and rounded once, which was 1 ULP off for e.g.
+    /// 1.6663333: `694ad53f` instead of `6a4ad53f`), double for `double`.
+    #[test]
+    fn scaleoffset_float_dscale_matches_libhdf5_bits() {
+        let file: &[u8] = include_bytes!("../tests/fixtures/filters/le_data.h5");
+        let cd = |size: u32, order: u32, fill_lo: u32, fill_hi: u32| {
+            let mut cd = vec![0, 3, 12, 1, size, 0, order, 1, fill_lo, fill_hi];
+            cd.resize(20, 0);
+            cd
+        };
+        let f32_le = cd(4, 0, 0xC00C_CCCD, 0);
+        let f32_be = cd(4, 1, 0xC00C_CCCD, 0);
+        let f64_le = cd(8, 0, 2576980378, 3221330329);
+        #[rustfmt::skip]
+        let cases: [(usize, usize, &[u32], &str); 6] = [
+            (2816, 38, &f32_le, "abaaaa3ed2942a3fec0a803fd2942a3fec0a803fabaaaa3fec0a803fabaaaa3f694ad53fabaaaa3f694ad53f76050040"),
+            (2854, 38, &f32_le, "abaaaa3f6a4ad53f760500406a4ad53f7605004056551540760500405655154034a52a405655154034a52a4076054040"),
+            (712, 38, &f32_be, "3eaaaaab3f2a94d23f800aec3f2a94d23f800aec3faaaaab3f800aec3faaaaab3fd54a693faaaaab3fd54a6940000576"),
+            (750, 38, &f32_be, "3faaaaab3fd54a6a400005763fd54a6a40000576401555564000057640155556402aa53440155556402aa53440400576"),
+            (2050, 38, &f64_le, concat!(
+                "555555555555d53fb9d75c489a52e53fce3e7c865d01f03fb9d75c489a52e53fce3e7c865d01f03f555555555555f53f",
+                "ce3e7c865d01f03f555555555555f53fdc6b2e244da9fa3f555555555555f53fdc6b2e244da9fa3f671f3ec3ae000040")),
+            (2088, 38, &f64_le, concat!(
+                "555555555555f53fdc6b2e244da9fa3f671f3ec3ae000040dc6b2e244da9fa3f671f3ec3ae000040aaaaaaaaaaaa0240",
+                "671f3ec3ae000040aaaaaaaaaaaa0240ee351792a6540540aaaaaaaaaaaa0240ee351792a6540540671f3ec3ae000840")),
+        ];
+        for (off, len, cd, want) in cases {
+            let pipeline = FilterPipeline {
+                version: 2,
+                filters: vec![one_filter(FILTER_SCALEOFFSET, cd.to_vec())],
+            };
+            let want = unhex(want);
+            let got =
+                decompress_chunk(&file[off..off + len], &pipeline, want.len(), cd[4]).unwrap();
+            assert_eq!(got, want, "chunk at {off}");
+        }
+    }
+
     /// `le_data.h5` `/Nbit_float_data_{le,be}` chunk (0,0): a 20-bit float
     /// (offset 7) packed by N-Bit. The filter must reproduce libhdf5's
     /// decoded bytes in the *file* datatype (h5py `DatasetID.read` with the