CHANGELOG.md

@@ -299,6 +299,12 @@
   flag; Fletcher32 ahead of deflate (NetCDF-4's order). Unknown-message flags
   follow libhdf5 (`tbogus.h5`): "fail if unknown" is refused, "fail if unknown
   and writing" is ignored by a reader.
+- `clawhdf5-format` reader — dense groups and attributes (links or
+  attributes kept in a fractal heap indexed by a v2 B-tree):
+  - A link heap larger than the root indirect block's direct rows (512 KiB
+    with libhdf5's defaults: a few thousand long link names, or ~20 000 short
+    ones) could not be listed: child indirect blocks were given the wrong
+    number of rows, so every link stored in one was unreachable.
 - `clawhdf5-format` writer — **files libhdf5 rejects or reads wrong:**
   - Extensible Array (one unlimited dimension): chunks from index 244 on were
     written but never indexed and read as 0, by libhdf5 and by us.

crates/clawhdf5-format/src/fractal_heap.rs

@@ -418,35 +418,35 @@ impl FractalHeapHeader {
         }
 
         // If we have indirect block rows
+        // A child indirect block in row r spans exactly that row's block size
+        // of heap space, so it has as many rows as a table of that total size
+        // needs (not `row - start_indirect + 1`, which undercounts and makes
+        // every object past the root's direct rows unreachable).
         for row in start_indirect..nrows_usize {
-            let _block_size = self.block_size_for_row(row);
-            let child_nrows = row - start_indirect + 1;
+            let child_space = self.block_size_for_row(row);
+            let child_nrows = self.rows_for_size(child_space);
 
             for _col in 0..tw {
                 let child_addr = read_offset(file_data, pos, offset_size)?;
                 pos += offset_size as usize;
 
-                if !is_undefined(child_addr, offset_size) {
-                    // Calculate total heap space covered by this indirect block child
-                    let total_child_space = self.indirect_block_heap_size(child_nrows);
-                    let block_end = current_heap_offset + total_child_space;
-                    if target_offset >= current_heap_offset && target_offset < block_end {
-                        return self.read_from_indirect_block(
-                            file_data,
-                            child_addr as usize,
-                            child_nrows as u16,
-                            current_heap_offset,
-                            target_offset,
-                            length,
-                            offset_size,
-                            depth_remaining - 1,
-                        );
-                    }
-                    current_heap_offset += total_child_space;
-                } else {
-                    let total_child_space = self.indirect_block_heap_size(child_nrows);
-                    current_heap_offset += total_child_space;
+                let block_end = current_heap_offset.saturating_add(child_space);
+                if !is_undefined(child_addr, offset_size)
+                    && target_offset >= current_heap_offset
+                    && target_offset < block_end
+                {
+                    return self.read_from_indirect_block(
+                        file_data,
+                        child_addr as usize,
+                        child_nrows,
+                        current_heap_offset,
+                        target_offset,
+                        length,
+                        offset_size,
+                        depth_remaining - 1,
+                    );
                 }
+                current_heap_offset = block_end;
             }
         }
 
@@ -475,25 +475,24 @@ impl FractalHeapHeader {
         log2 + 2
     }
 
+    /// Rows an indirect block needs to span `size` bytes of heap space:
+    /// `log2(size) - log2(starting_block_size * table_width) + 1`, as
+    /// libhdf5's `H5HF__dtable_size_to_rows`.
+    fn rows_for_size(&self, size: u64) -> u16 {
+        let log2 = |v: u64| 63u32.saturating_sub(v.max(1).leading_zeros());
+        let first_row_bits = log2(self.starting_block_size) + log2(u64::from(self.table_width));
+        (log2(size).saturating_sub(first_row_bits) + 1) as u16
+    }
+
     /// Get block size for a given row in the doubling table.
     fn block_size_for_row(&self, row: usize) -> u64 {
         let sbs = self.starting_block_size;
         if row <= 1 {
             sbs
         } else {
-            sbs * (1u64 << (row - 1))
+            sbs.saturating_mul(1u64.checked_shl((row - 1) as u32).unwrap_or(u64::MAX))
         }
     }
-
-    /// Total heap space covered by an indirect block with the given number of rows.
-    fn indirect_block_heap_size(&self, nrows: usize) -> u64 {
-        let tw = self.table_width as u64;
-        let mut total = 0u64;
-        for row in 0..nrows {
-            total += self.block_size_for_row(row) * tw;
-        }
-        total
-    }
 }
 
 #[cfg(test)]

crates/clawhdf5/tests/dense_storage_interop.rs

@@ -0,0 +1,143 @@
+//! Dense ("new-style") link and attribute storage written by libhdf5 (via
+//! h5py): groups whose links live in a fractal heap indexed by a v2 B-tree,
+//! and objects whose attributes do. Every listing and value is compared with
+//! what h5py itself reports for the same file.
+//!
+//! Skipped when python3 with h5py is unavailable, unless
+//! `CLAWHDF5_REQUIRE_INTEROP=1`.
+
+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"])
+        .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;
+        }
+    };
+}
+
+fn run_python(script: &str) -> String {
+    let output = Command::new(python())
+        .args(["-c", script])
+        .output()
+        .expect("failed to run python");
+    if !output.status.success() {
+        panic!(
+            "Python script failed:\nSTDOUT: {}\nSTDERR: {}",
+            String::from_utf8_lossy(&output.stdout),
+            String::from_utf8_lossy(&output.stderr)
+        );
+    }
+    String::from_utf8_lossy(&output.stdout).trim().to_string()
+}
+
+/// Have h5py write a file with `body` (which sees `f`, `h5py` and `np`),
+/// returning the temp dir holding it and its path.
+fn h5py_file(body: &str) -> (tempfile::TempDir, String) {
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("dense.h5").display().to_string();
+    let script = format!(
+        "import h5py, numpy as np\n\
+         with h5py.File(r'{path}', 'w', libver='latest') as f:\n{}",
+        indent(body)
+    );
+    run_python(&script);
+    (dir, path)
+}
+
+fn indent(body: &str) -> String {
+    body.lines()
+        .map(|l| format!("    {l}\n"))
+        .collect::<String>()
+}
+
+/// The datasets and groups h5py lists in `group`, sorted: links h5py can
+/// resolve (hard and soft), without dangling soft links or external links.
+fn h5py_listing(path: &str, group: &str) -> (Vec<String>, Vec<String>) {
+    let out = run_python(&format!(
+        "import h5py\n\
+         ds, gs = [], []\n\
+         with h5py.File(r'{path}', 'r') as f:\n\
+         \x20   g = f[{group:?}]\n\
+         \x20   for k in g.keys():\n\
+         \x20       if isinstance(g.get(k, getlink=True), h5py.ExternalLink):\n\
+         \x20           continue\n\
+         \x20       try:\n\
+         \x20           o = g[k]\n\
+         \x20       except Exception:\n\
+         \x20           continue\n\
+         \x20       (ds if isinstance(o, h5py.Dataset) else gs).append(k)\n\
+         print('\\x1f'.join(sorted(ds)))\n\
+         print('\\x1f'.join(sorted(gs)))\n"
+    ));
+    let mut lines = out.lines();
+    let split = |l: Option<&str>| -> Vec<String> {
+        l.unwrap_or("")
+            .split('\x1f')
+            .filter(|s| !s.is_empty())
+            .map(str::to_string)
+            .collect()
+    };
+    let ds = split(lines.next());
+    let gs = split(lines.next());
+    (ds, gs)
+}
+
+fn our_listing(path: &str, group: &str) -> (Vec<String>, Vec<String>) {
+    let f = File::open(path).unwrap();
+    let g = f.group(group).unwrap();
+    let mut ds = g.datasets().unwrap();
+    let mut gs = g.groups().unwrap();
+    ds.sort();
+    gs.sort();
+    (ds, gs)
+}
+
+fn assert_same_listing(path: &str, group: &str) {
+    let ours = our_listing(path, group);
+    let theirs = h5py_listing(path, group);
+    assert_eq!(ours.0.len(), theirs.0.len(), "dataset count in {group}");
+    assert_eq!(ours, theirs, "listing of {group}");
+}
+
+#[test]
+fn dense_group_whose_heap_outgrows_the_root_direct_rows() {
+    skip_if_no_python!();
+    // Long link names make the link heap larger than the root indirect
+    // block's direct rows can hold (512 KiB with h5py's defaults), so links
+    // live in child indirect blocks. Those were sized from the wrong row
+    // count, and every link past the direct rows was unreachable.
+    let (_dir, path) = h5py_file(
+        "t = f.create_dataset('t', data=[1.0])\n\
+         g = f.create_group('g')\n\
+         for i in range(2500):\n\
+         \x20   g['n%05d_' % i + 'x' * 240] = t\n",
+    );
+    assert_same_listing(&path, "g");
+    let f = File::open(&path).unwrap();
+    let last = format!("g/n02499_{}", "x".repeat(240));
+    assert_eq!(f.dataset(&last).unwrap().read_f64().unwrap(), vec![1.0]);
+}

docs/known-issues.md

@@ -77,6 +77,8 @@ the VDS item, which is marked.
 - **Groups and links:**
   - Groups with a user-defined link type (e.g. 187) cannot be listed.
   - Dense groups with more than about 22 000 links cannot be listed.
+    *Partly fixed 2026-09-25:* a link heap past its root block's direct rows
+    (child indirect blocks) is now read.
   - Soft links are left out of `datasets()`.
 - **Dense attributes:** a large attribute stored as a fractal-heap "huge"
   object makes every attribute on the object fail. This affects real NetCDF