filters: Fletcher-32 as libhdf5 computes it

Our checksum reduced its sums with `% 65535`; libhdf5's
H5_checksum_fletcher32 folds them with `(s & 0xffff) + (s >> 16)`, which
leaves 0xffff where the modulo leaves 0. On about one chunk in 32768
libhdf5 refused the chunks we wrote and we refused the chunks it wrote.
Every release since v2.1.0 is affected.

clawhdf5_format::checksum::fletcher32 is a port of H5_checksum_fletcher32
and the filter's only implementation. Verification also accepts the
byte-swapped form libhdf5 accepts (1.6.2 and earlier) and the `% 65535`
form earlier releases wrote, so their files stay readable.

The new interop test compares the checksum with libhdf5's own function
(ctypes) on every 1- and 2-byte input and 40 000 random and fold-heavy
inputs, and moves fold-case chunks between h5py and FileBuilder/FileEditor
in both directions; with the old filters.rs the three file tests fail.

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

crates/clawhdf5-format/src/checksum.rs

@@ -14,6 +14,45 @@ pub fn jenkins_lookup3(data: &[u8]) -> u32 {
     hashlittle(data, 0)
 }
 
+/// HDF5's Fletcher-32 checksum, as the Fletcher-32 I/O filter (filter id 3)
+/// stores it after each chunk.
+///
+/// A line-for-line port of `H5_checksum_fletcher32` (H5checksum.c, libhdf5
+/// 1.8 through 1.14): big-endian 16-bit words summed in blocks of 360, each
+/// sum reduced after a block by the ones'-complement fold
+/// `(s & 0xffff) + (s >> 16)` rather than `% 65535`, an odd trailing byte
+/// taken as the high byte of a last word, and a final fold of both sums.
+/// The fold and `% 65535` differ whenever a sum is a non-zero multiple of
+/// 65535: the fold leaves 0xffff where the modulo gives 0, so the two
+/// disagree on about one chunk in 32768 and libhdf5 rejects the other's
+/// checksum. This must stay the only implementation.
+pub fn fletcher32(data: &[u8]) -> u32 {
+    let mut sum1: u32 = 0;
+    let mut sum2: u32 = 0;
+    // 360 words keep both sums inside 32 bits between folds (the bound
+    // libhdf5 uses: after a fold sum1 < 0x10200, so sum2 stays below
+    // 360 * 361 / 2 * 0xffff + 360 * 0x10200 + 0x1fffe < 2^32). The adds wrap
+    // like the C unsigned arithmetic all the same.
+    let (words, odd) = data.as_chunks::<2>();
+    for block in words.chunks(360) {
+        for w in block {
+            sum1 = sum1.wrapping_add((u32::from(w[0]) << 8) | u32::from(w[1]));
+            sum2 = sum2.wrapping_add(sum1);
+        }
+        sum1 = (sum1 & 0xffff) + (sum1 >> 16);
+        sum2 = (sum2 & 0xffff) + (sum2 >> 16);
+    }
+    if let [last] = odd {
+        sum1 = sum1.wrapping_add(u32::from(*last) << 8);
+        sum2 = sum2.wrapping_add(sum1);
+        sum1 = (sum1 & 0xffff) + (sum1 >> 16);
+        sum2 = (sum2 & 0xffff) + (sum2 >> 16);
+    }
+    sum1 = (sum1 & 0xffff) + (sum1 >> 16);
+    sum2 = (sum2 & 0xffff) + (sum2 >> 16);
+    (sum2 << 16) | sum1
+}
+
 /// Compute CRC32 (IEEE / ISO 3309) over data.
 ///
 /// When the `fast-checksum` feature is enabled, this uses hardware CRC32
@@ -207,6 +246,19 @@ fn hashlittle(data: &[u8], initval: u32) -> u32 {
 mod tests {
     use super::*;
 
+    /// Values of libhdf5's `H5_checksum_fletcher32` (h5py 3.x's bundled
+    /// libhdf5, called through ctypes). The first three are sums that are
+    /// multiples of 65535, where `% 65535` gave 0 instead of 0xffff.
+    #[test]
+    fn fletcher32_matches_libhdf5() {
+        assert_eq!(fletcher32(&[0x00, 0x01, 0xff, 0xfe]), 0x0001_ffff);
+        assert_eq!(fletcher32(&[0xff; 720]), 0xffff_ffff);
+        assert_eq!(fletcher32(&[0xff; 721]), 0xff00_ff00);
+        assert_eq!(fletcher32(&[0xff; 1441]), 0xff00_ff00);
+        assert_eq!(fletcher32(&[]), 0);
+        assert_eq!(fletcher32(&[7]), 0x0700_0700);
+    }
+
     #[test]
     fn empty_input() {
         // Empty input should return the initial state after no mixing

crates/clawhdf5-format/src/filters.rs

@@ -1678,56 +1678,6 @@ fn shuffle_compress_general(data: &[u8], n: usize, element_size: usize, result:
     }
 }
 
-/// Compute HDF5 Fletcher32 checksum over data.
-/// HDF5 uses a modified Fletcher32 that operates on 16-bit words.
-///
-/// Optimized with wider accumulators: processes blocks of 360 words before
-/// taking the modulo, reducing the number of expensive modulo operations.
-/// (360 is the maximum block size that avoids u32 overflow for sum2.)
-fn fletcher32_compute(data: &[u8]) -> u32 {
-    let mut sum1: u32 = 0;
-    let mut sum2: u32 = 0;
-
-    // Process in blocks of 360 16-bit words (720 bytes) to delay modulo.
-    // Max sum1 before mod: 360 * 65535 = 23_592_600 < u32::MAX
-    // Max sum2 before mod: 360 * 23_592_600 ~ 8.5B > u32::MAX, but actual
-    // sum2 accumulates incrementally, so worst case is 360*360*65535/2 which
-    // fits in u64. We use u32 with block size 360 which is safe.
-    const BLOCK_WORDS: usize = 360;
-    const BLOCK_BYTES: usize = BLOCK_WORDS * 2;
-
-    let mut offset = 0;
-    let len = data.len();
-
-    while offset + BLOCK_BYTES <= len {
-        let end = offset + BLOCK_BYTES;
-        let mut i = offset;
-        while i < end {
-            let val = ((data[i] as u32) << 8) | (data[i + 1] as u32);
-            sum1 += val;
-            sum2 += sum1;
-            i += 2;
-        }
-        sum1 %= 65535;
-        sum2 %= 65535;
-        offset = end;
-    }
-
-    // Handle remaining bytes
-    while offset < len {
-        let val = if offset + 1 < len {
-            ((data[offset] as u32) << 8) | (data[offset + 1] as u32)
-        } else {
-            (data[offset] as u32) << 8
-        };
-        sum1 = (sum1 + val) % 65535;
-        sum2 = (sum2 + sum1) % 65535;
-        offset += 2;
-    }
-
-    (sum2 << 16) | sum1
-}
-
 /// Verify Fletcher32 checksum and strip it from the data.
 /// The last 4 bytes are the stored checksum.
 fn fletcher32_verify(data: &[u8]) -> Result<Vec<u8>, FormatError> {
@@ -1749,8 +1699,18 @@ fn fletcher32_payload(data: &[u8]) -> Result<usize, FormatError> {
         data[data.len() - 2],
         data[data.len() - 1],
     ]);
-    let computed = fletcher32_compute(payload);
-    if stored != computed {
+    let computed = crate::checksum::fletcher32(payload);
+    // libhdf5 also accepts the checksum with the bytes of each 16-bit half
+    // swapped, which is how 1.6.2 and earlier stored it
+    // (H5Z__filter_fletcher32's `reversed_fletcher`).
+    let reversed = ((computed & 0x00ff_00ff) << 8) | ((computed >> 8) & 0x00ff_00ff);
+    // clawhdf5 v2.7.0 and earlier reduced the sums `% 65535`, which gives 0
+    // where libhdf5's fold gives 0xffff; accept that form too, so that files
+    // those releases wrote can still be read (and rewritten for libhdf5).
+    // It differs from `computed` only in a half that is 0xffff.
+    let half = |h: u32| if h == 0xffff { 0 } else { h };
+    let legacy = (half(computed >> 16) << 16) | half(computed & 0xffff);
+    if stored != computed && stored != reversed && stored != legacy {
         return Err(FormatError::Fletcher32Mismatch {
             expected: stored,
             computed,
@@ -1761,7 +1721,7 @@ fn fletcher32_payload(data: &[u8]) -> Result<usize, FormatError> {
 
 /// Append Fletcher32 checksum to data.
 fn fletcher32_append(data: &[u8]) -> Result<Vec<u8>, FormatError> {
-    let checksum = fletcher32_compute(data);
+    let checksum = crate::checksum::fletcher32(data);
     let mut result = data.to_vec();
     result.extend_from_slice(&checksum.to_le_bytes());
     Ok(result)