diff --git a/crates/clawhdf5-format/Cargo.toml b/crates/clawhdf5-format/Cargo.toml index c1ec6df..21bc626 100644 --- a/crates/clawhdf5-format/Cargo.toml +++ b/crates/clawhdf5-format/Cargo.toml @@ -78,8 +78,11 @@ bzip2 = ["dep:bzip2", "std"] blosc = ["lz4_flex", "ruzstd", "snap", "deflate", "std"] # Blosc2 (32026), read-only: frames, B2ND arrays, and the Blosc codecs above. blosc2 = ["blosc"] +# ZFP (32013, H5Z-ZFP), read-only: every mode, for int32, int64, float and +# double fields of 1 to 4 dimensions. +zfp = [] # Every plugin filter above. -plugin-filters = ["lzf", "bitshuffle", "bzip2", "blosc", "blosc2"] +plugin-filters = ["lzf", "bitshuffle", "bzip2", "blosc", "blosc2", "zfp"] [[bench]] name = "parallel_decompress_bench" diff --git a/crates/clawhdf5-format/src/filter_pipeline.rs b/crates/clawhdf5-format/src/filter_pipeline.rs index 13ed69f..72f61bc 100644 --- a/crates/clawhdf5-format/src/filter_pipeline.rs +++ b/crates/clawhdf5-format/src/filter_pipeline.rs @@ -27,7 +27,8 @@ pub const FILTER_LZF: u16 = 32000; pub const FILTER_BLOSC: u16 = 32001; /// Bitshuffle, optionally with LZ4 or Zstandard (hdf5plugin's `Bitshuffle`). pub const FILTER_BITSHUFFLE: u16 = 32008; -/// ZFP lossy floating-point compression (hdf5plugin's `Zfp`). Not supported. +/// ZFP lossy (and lossless) compression of numeric arrays (H5Z-ZFP; +/// hdf5plugin's `Zfp`). Read-only, with the `zfp` feature. pub const FILTER_ZFP: u16 = 32013; /// Blosc 2 (hdf5plugin's `Blosc2`). pub const FILTER_BLOSC2: u16 = 32026; diff --git a/crates/clawhdf5-format/src/filter_registry.rs b/crates/clawhdf5-format/src/filter_registry.rs index 442e27c..7670a8b 100644 --- a/crates/clawhdf5-format/src/filter_registry.rs +++ b/crates/clawhdf5-format/src/filter_registry.rs @@ -6,7 +6,7 @@ //! build: the HDF5 standard filters (deflate, shuffle, Fletcher32, szip, //! N-Bit, scale-offset) and the plugin filters whose cargo features are //! enabled (LZ4, Zstandard, pcodec, LZF, bitshuffle, bzip2, blosc, -//! blosc2). +//! blosc2, zfp). //! [`builtin_filters`] lists them. //! * **Registered filters** (`std` only) — codecs the application supplies //! for any other ID with [`register_filter`] (a [`FilterCodec`], or just a @@ -162,7 +162,7 @@ pub fn known_filter(id: u16) -> Option<(&'static str, Option<&'static str>)> { 32001 => ("Blosc", Some("blosc")), 32004 => ("LZ4", Some("lz4")), 32008 => ("bitshuffle", Some("bitshuffle")), - 32013 => ("ZFP", None), + 32013 => ("ZFP", Some("zfp")), 32015 => ("Zstandard", Some("zstd")), 32019 => ("JPEG", None), 32022 => ("BitGroom", None), @@ -455,8 +455,10 @@ pub(crate) mod tests { let msg = FormatError::UnsupportedFilter(32026).to_string(); assert!(msg.contains("Blosc2") && msg.contains("`blosc2`"), "{msg}"); let msg = FormatError::UnsupportedFilter(32013).to_string(); + assert!(msg.contains("ZFP") && msg.contains("`zfp`"), "{msg}"); + let msg = FormatError::UnsupportedFilter(32019).to_string(); assert!( - msg.contains("ZFP") && msg.contains("not implemented"), + msg.contains("JPEG") && msg.contains("not implemented"), "{msg}" ); let msg = FormatError::UnsupportedFilter(32000).to_string(); diff --git a/crates/clawhdf5-format/src/filters.rs b/crates/clawhdf5-format/src/filters.rs index 40d6cd9..b261284 100644 --- a/crates/clawhdf5-format/src/filters.rs +++ b/crates/clawhdf5-format/src/filters.rs @@ -432,6 +432,13 @@ pub(crate) static BUILTIN_FILTERS: &[BuiltinFilter] = &[ decode: crate::filters_bitshuffle::bitshuffle_decode, encode: Some(crate::filters_bitshuffle::bitshuffle_encode), }, + #[cfg(feature = "zfp")] + BuiltinFilter { + id: crate::filter_pipeline::FILTER_ZFP, + name: "zfp", + decode: crate::filters_zfp::zfp_decode, + encode: None, + }, #[cfg(feature = "zstd")] BuiltinFilter { id: FILTER_ZSTD, diff --git a/crates/clawhdf5-format/src/filters_zfp.rs b/crates/clawhdf5-format/src/filters_zfp.rs new file mode 100644 index 0000000..7ac775c --- /dev/null +++ b/crates/clawhdf5-format/src/filters_zfp.rs @@ -0,0 +1,1083 @@ +//! ZFP (HDF5 filter 32013, `H5Z-ZFP`, hdf5plugin's `Zfp`) in pure Rust, +//! read-only: a port of the zfp 1.0.1 decoder (`src/zfp.c`, +//! `src/template/decode*.c`, `revdecode*.c`, `include/zfp/bitstream.inl`) +//! and of the decompression half of H5Z-ZFP 1.1.1's `H5Zzfp.c`. +//! +//! **The filter.** H5Z-ZFP keeps the zfp header in the filter's +//! `cd_values`, not in the chunks: `[0]` packs the zfp library version +//! (bits 16-31), the zfp codec version (bits 12-15; 0 before H5Z-ZFP 1.1.0, +//! when it is inferred from the library version) and the filter version +//! (bits 0-11); `[1..]` (at most 6 words) are a zfp bit stream holding the +//! full header written by `zfp_write_header`. Each chunk is the bare +//! compressed stream. The words are read as little-endian bytes; if the +//! magic is not there they are read byte-swapped (a big-endian writer), and +//! then the decoded values are byte-swapped too, as H5Z-ZFP does, because +//! the dataset's datatype is big-endian. +//! +//! **The header** (LSB-first bit stream; H5Z-ZFP requires zfp built with +//! 8-bit stream words, so the stream is plain bytes): 32 bits of magic +//! (`'z' 'f' 'p'` and codec version 5); 52 bits of field metadata (scalar +//! type: int32, int64, float, double; dimensionality 1-4; the sizes, 48 +//! bits shared between the dimensions, the fastest-varying first); and a +//! 12-bit mode, or 64 bits when the 12 are all ones: fixed rate (maxbits +//! per block), fixed precision (bit planes), fixed accuracy (the smallest +//! bit plane kept), reversible (lossless), or the four expert parameters +//! (minbits, maxbits, maxprec, minexp). +//! +//! **The stream.** The field is cut into blocks of 4^d values (partial at +//! the array's edges; the decoder decodes whole blocks and keeps the part +//! inside the array), decoded one after another. A floating-point block is +//! a "nonzero" bit, the block's common exponent (8 or 11 bits) and a block +//! of integers; an integer block of the array is that block of integers +//! directly. Integers are coded as negabinary bit planes, most significant +//! first, each plane as the bits of the values already significant plus a +//! unary run-length ("group test") code for the rest, within the budget of +//! maxbits and maxprec bits (the precision of a float block also depends on +//! its exponent and minexp). The coefficients are stored in order of +//! sequency, then run through the inverse decorrelating transform (a +//! lifting scheme along each dimension) and, for floats, scaled by the +//! block exponent. Blocks that used fewer than minbits bits are padded. +//! Reversible mode stores the precision in the block, uses an integer +//! Lorenzo transform, and stores float blocks either as integers scaled by +//! a common exponent or as their bit patterns. +//! +//! The decoder is deterministic and reproduces libzfp's output bit for bit +//! (libzfp built with its default `ZFP_ROUND_NEVER` rounding, as hdf5plugin +//! builds it): integer arithmetic wraps as libzfp's does, and floats are +//! scaled by an exact power of two. A stream that ends before the decoder +//! is done is an error (libzfp reads past the buffer). + +#[cfg(not(feature = "std"))] +use alloc::{format, vec, vec::Vec}; + +use crate::error::FormatError; +use crate::filter_registry::FilterContext; + +fn err(msg: &str) -> FormatError { + FormatError::DecompressionError(format!("zfp: {msg}")) +} + +const ZFP_MIN_BITS: u32 = 1; +const ZFP_MAX_BITS: u32 = 16658; +const ZFP_MAX_PREC: u32 = 64; +const ZFP_MIN_EXP: i32 = -1074; +const ZFP_META_BITS: u32 = 52; +const ZFP_MODE_SHORT_BITS: u32 = 12; +const ZFP_MODE_LONG_BITS: u32 = 64; +const ZFP_MODE_SHORT_MAX: u64 = (1 << ZFP_MODE_SHORT_BITS) - 2; +/// The zfp codec this decoder implements (zfp 0.5.0 to 1.0.1). +const ZFP_CODEC: u64 = 5; +/// `H5Z_ZFP_CD_NELMTS_MAX`: header words after the version word. +const H5Z_ZFP_CD_NELMTS_MAX: usize = 6; + +/// An LSB-first bit stream over bytes (zfp's `bitstream` with 8-bit words). +struct Bits<'a> { + data: &'a [u8], + pos: u64, +} + +fn truncated() -> FormatError { + err("compressed stream ends early") +} + +impl<'a> Bits<'a> { + fn new(data: &'a [u8]) -> Self { + Bits { data, pos: 0 } + } + + #[inline] + fn bit(&mut self) -> Result { + let byte = *self + .data + .get((self.pos >> 3) as usize) + .ok_or_else(truncated)?; + let v = (byte >> (self.pos & 7)) & 1; + self.pos += 1; + Ok(v as u32) + } + + /// The next `n` (0..=64) bits, the first in the least significant place. + fn bits(&mut self, n: u32) -> Result { + if n == 0 { + return Ok(0); + } + if self.pos + n as u64 > self.data.len() as u64 * 8 { + return Err(truncated()); + } + let mut v = 0u64; + let mut got = 0u32; + while got < n { + let byte = self.data[(self.pos >> 3) as usize] as u64; + let off = (self.pos & 7) as u32; + let take = (8 - off).min(n - got); + v |= ((byte >> off) & ((1u64 << take) - 1)) << got; + got += take; + self.pos += take as u64; + } + Ok(v) + } + + /// Skip `n` bits. Nothing is read, so skipping past the end is not an + /// error until a bit there is read. + fn skip(&mut self, n: u32) { + self.pos += n as u64; + } + + fn tell(&self) -> u64 { + self.pos + } +} + +/// A zfp stream's four parameters (`zfp_stream` minus the bit stream). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +struct Params { + minbits: u32, + maxbits: u32, + maxprec: u32, + minexp: i32, +} + +impl Params { + fn reversible(&self) -> bool { + self.minexp < ZFP_MIN_EXP + } + + /// `zfp_stream_set_params`. + fn new(minbits: u32, maxbits: u32, maxprec: u32, minexp: i32) -> Option { + if minbits > maxbits || !(0 < maxprec && maxprec <= 64) { + return None; + } + Some(Params { + minbits, + maxbits, + maxprec, + minexp, + }) + } + + /// `zfp_stream_set_mode`. + fn from_mode(mut mode: u64) -> Option { + let (minbits, maxbits, maxprec, minexp); + if mode <= ZFP_MODE_SHORT_MAX { + if mode < 2048 { + minbits = mode as u32 + 1; + maxbits = minbits; + maxprec = ZFP_MAX_PREC; + minexp = ZFP_MIN_EXP; + } else if mode < 2048 + 128 { + minbits = ZFP_MIN_BITS; + maxbits = ZFP_MAX_BITS; + maxprec = mode as u32 + 1 - 2048; + minexp = ZFP_MIN_EXP; + } else if mode == 2048 + 128 { + minbits = ZFP_MIN_BITS; + maxbits = ZFP_MAX_BITS; + maxprec = ZFP_MAX_PREC; + minexp = ZFP_MIN_EXP - 1; + } else { + minbits = ZFP_MIN_BITS; + maxbits = ZFP_MAX_BITS; + maxprec = ZFP_MAX_PREC; + minexp = mode as i32 + ZFP_MIN_EXP - (2048 + 128 + 1); + } + } else { + mode >>= 12; + minbits = (mode & 0x7fff) as u32 + 1; + mode >>= 15; + maxbits = (mode & 0x7fff) as u32 + 1; + mode >>= 15; + maxprec = (mode & 0x7f) as u32 + 1; + mode >>= 7; + minexp = (mode & 0x7fff) as i32 - 16495; + } + Params::new(minbits, maxbits, maxprec, minexp) + } + + /// `zfp_stream_compression_mode`: 1 expert, 2 fixed rate, 3 fixed + /// precision, 4 fixed accuracy, 5 reversible. + fn compression_mode(&self) -> u32 { + let p = self; + if p.minbits == ZFP_MIN_BITS + && p.maxbits == ZFP_MAX_BITS + && p.maxprec == ZFP_MAX_PREC + && p.minexp == ZFP_MIN_EXP + { + return 1; + } + if p.minbits == p.maxbits + && 1 <= p.maxbits + && p.maxbits <= ZFP_MAX_BITS + && p.maxprec >= ZFP_MAX_PREC + && p.minexp == ZFP_MIN_EXP + { + return 2; + } + if p.minbits <= ZFP_MIN_BITS + && p.maxbits >= ZFP_MAX_BITS + && p.maxprec >= 1 + && p.minexp == ZFP_MIN_EXP + { + return 3; + } + if p.minbits <= ZFP_MIN_BITS + && p.maxbits >= ZFP_MAX_BITS + && p.maxprec >= ZFP_MAX_PREC + && p.minexp >= ZFP_MIN_EXP + { + return 4; + } + if p.minbits <= ZFP_MIN_BITS + && p.maxbits >= ZFP_MAX_BITS + && p.maxprec >= ZFP_MAX_PREC + && p.minexp < ZFP_MIN_EXP + { + return 5; + } + 1 + } + + /// `zfp_stream_mode`: the mode word these parameters are written as. + fn mode(&self) -> u64 { + let p = self; + match p.compression_mode() { + 2 if p.maxbits <= 2048 => return (p.maxbits - 1) as u64, + 3 if p.maxprec <= 128 => return (p.maxprec - 1 + 2048) as u64, + 4 if p.minexp <= 843 => return (p.minexp - ZFP_MIN_EXP) as u64 + (2048 + 128 + 1), + 5 => return 2048 + 128, + _ => {} + } + let minbits = p.minbits.clamp(1, 0x8000) - 1; + let maxbits = p.maxbits.clamp(1, 0x8000) - 1; + let maxprec = p.maxprec.clamp(1, 0x80) - 1; + let minexp = (p.minexp.saturating_add(16495)).clamp(0, 0x7fff) as u64; + let mut mode = minexp; + mode = (mode << 7) + maxprec as u64; + mode = (mode << 15) + maxbits as u64; + mode = (mode << 15) + minbits as u64; + (mode << 12) + 0xfff + } +} + +/// Scalar types (`zfp_type`). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum ZType { + Int32, + Int64, + Float, + Double, +} + +impl ZType { + fn size(self) -> usize { + match self { + ZType::Int32 | ZType::Float => 4, + ZType::Int64 | ZType::Double => 8, + } + } +} + +/// A field's type and sizes (`zfp_field` without the data), from the +/// header's 52-bit metadata (`zfp_field_set_metadata`). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +struct Field { + ztype: ZType, + dims: usize, + /// nx (fastest), ny, nz, nw; 1 for the unused ones. + n: [usize; 4], +} + +impl Field { + fn from_meta(mut meta: u64) -> Field { + let ztype = match meta & 3 { + 0 => ZType::Int32, + 1 => ZType::Int64, + 2 => ZType::Float, + _ => ZType::Double, + }; + meta >>= 2; + let dims = (meta & 3) as usize + 1; + meta >>= 2; + let mut n = [1usize; 4]; + match dims { + // zfp limits 1-D sizes to 32 bits (of the 48 in the header). + 1 => n[0] = (meta & 0xffff_ffff) as usize + 1, + 2 => { + for v in n.iter_mut().take(2) { + *v = (meta & 0xff_ffff) as usize + 1; + meta >>= 24; + } + } + 3 => { + for v in n.iter_mut().take(3) { + *v = (meta & 0xffff) as usize + 1; + meta >>= 16; + } + } + _ => { + for v in n.iter_mut() { + *v = (meta & 0xfff) as usize + 1; + meta >>= 12; + } + } + } + Field { ztype, dims, n } + } + + fn values(&self) -> Option { + self.n.iter().try_fold(1usize, |a, &v| a.checked_mul(v)) + } + + fn blocks(&self) -> Option { + self.n + .iter() + .try_fold(1usize, |a, &v| a.checked_mul(v.div_ceil(4))) + } +} + +/// What the filter's `cd_values` say: the stream parameters, the field, +/// and whether the header was written big-endian. +#[derive(Debug)] +struct Header { + params: Params, + field: Field, + swap: bool, +} + +/// `zfp_read_header` with `ZFP_HEADER_MAGIC`. +fn read_magic(b: &mut Bits<'_>) -> Result { + Ok(b.bits(8)? == b'z' as u64 + && b.bits(8)? == b'f' as u64 + && b.bits(8)? == b'p' as u64 + && b.bits(8)? == ZFP_CODEC) +} + +/// `get_zfp_info_from_cd_values` and `zfp_codec_version_mismatch`. +fn parse_cd_values(cd: &[u32]) -> Result { + let (&version, words) = cd + .split_first() + .ok_or_else(|| err("no filter parameters"))?; + if words.len() > H5Z_ZFP_CD_NELMTS_MAX { + return Err(err("too many filter parameters")); + } + let bytes = |swap: bool| -> Vec { + words + .iter() + .flat_map(|w| { + if swap { + w.to_be_bytes() + } else { + w.to_le_bytes() + } + }) + .collect() + }; + let mut swap = false; + let mut hdr = bytes(false); + if !read_magic(&mut Bits::new(&hdr)).unwrap_or(false) { + swap = true; + hdr = bytes(true); + if !read_magic(&mut Bits::new(&hdr)).unwrap_or(false) { + return Err(err( + "no zfp header in the filter parameters, or a zfp codec other than 5", + )); + } + } + let mut b = Bits::new(&hdr); + let bad = |_| err("truncated zfp header in the filter parameters"); + b.skip(32); + let meta = b.bits(ZFP_META_BITS).map_err(bad)?; + let field = Field::from_meta(meta); + let mut mode = b.bits(ZFP_MODE_SHORT_BITS).map_err(bad)?; + if mode > ZFP_MODE_SHORT_MAX { + let rest = b + .bits(ZFP_MODE_LONG_BITS - ZFP_MODE_SHORT_BITS) + .map_err(bad)?; + mode += rest << ZFP_MODE_SHORT_BITS; + } + let params = Params::from_mode(mode).ok_or_else(|| err("invalid compression mode"))?; + // H5Z-ZFP hands the decoder `zfp_stream_mode()` of the header's + // parameters, which maps some expert settings to a standard mode's. + let params = Params::from_mode(params.mode()).ok_or_else(|| err("invalid compression mode"))?; + + // Data written by a newer codec than this decoder's. + let h5z_version = version & 0xfff; + let codec = (version >> 12) & 0xf; + let writer_codec = if h5z_version < 0x110 { + let v = (version >> 16) << 4; + if v < 0x0500 { + 4 + } else if v < 0x1000 { + (v & 0x0f00) >> 8 + } else { + 5 + } + } else { + codec + }; + if writer_codec as u64 > ZFP_CODEC { + return Err(err("data written by a newer zfp codec")); + } + Ok(Header { + params, + field, + swap, + }) +} + +/// Coefficient order by sequency (zfp's `perm_1` .. `perm_4`). +const PERM_1: [u8; 4] = [0, 1, 2, 3]; +const PERM_2: [u8; 16] = [0, 1, 4, 5, 2, 8, 6, 9, 3, 12, 10, 7, 13, 11, 14, 15]; +const PERM_3: [u8; 64] = [ + 0, 1, 4, 16, 20, 17, 5, 2, 8, 32, 21, 6, 18, 24, 9, 33, 36, 3, 12, 48, 22, 25, 37, 40, 34, 10, + 7, 19, 28, 13, 49, 52, 41, 38, 26, 23, 29, 53, 11, 35, 44, 14, 50, 56, 42, 27, 39, 45, 30, 54, + 57, 60, 51, 15, 43, 46, 58, 61, 55, 31, 62, 59, 47, 63, +]; +const PERM_4: [u8; 256] = [ + 0, 1, 4, 16, 64, 5, 80, 17, 68, 65, 20, 2, 8, 32, 128, 84, 81, 69, 21, 6, 18, 66, 24, 72, 9, + 96, 33, 36, 129, 132, 144, 3, 12, 48, 192, 85, 82, 70, 22, 73, 25, 88, 37, 100, 97, 148, 145, + 133, 10, 160, 34, 136, 130, 40, 7, 19, 67, 28, 76, 13, 112, 49, 52, 193, 196, 208, 86, 89, 101, + 149, 161, 137, 41, 134, 38, 164, 26, 152, 146, 104, 98, 74, 83, 71, 23, 77, 29, 92, 53, 116, + 113, 212, 209, 197, 11, 35, 131, 44, 140, 14, 176, 50, 56, 194, 200, 224, 90, 165, 102, 153, + 150, 105, 168, 162, 138, 42, 87, 93, 117, 213, 27, 75, 99, 39, 135, 147, 108, 45, 141, 156, 30, + 78, 177, 180, 54, 114, 120, 57, 198, 210, 216, 201, 225, 228, 15, 240, 51, 204, 195, 60, 169, + 166, 154, 106, 91, 103, 151, 109, 157, 94, 181, 118, 121, 214, 217, 229, 163, 139, 43, 142, 46, + 172, 58, 184, 178, 232, 226, 202, 241, 205, 61, 199, 55, 244, 31, 220, 211, 124, 115, 79, 170, + 167, 155, 107, 158, 110, 173, 122, 185, 182, 233, 230, 218, 95, 245, 119, 221, 215, 125, 242, + 206, 62, 203, 59, 248, 47, 236, 227, 188, 179, 143, 171, 174, 186, 234, 246, 222, 126, 219, + 123, 249, 111, 237, 231, 189, 183, 159, 252, 243, 207, 63, 175, 250, 187, 238, 235, 190, 253, + 247, 223, 127, 254, 251, 239, 191, 255, +]; + +fn perm(dims: usize) -> &'static [u8] { + match dims { + 1 => &PERM_1, + 2 => &PERM_2, + 3 => &PERM_3, + _ => &PERM_4, + } +} + +/// `with_maxbits`: whether the block's bit budget may run out before its +/// precision does. +fn with_maxbits(maxbits: u32, maxprec: u32, size: u32) -> bool { + ((maxprec + 1) * size).wrapping_sub(1) > maxbits +} + +/// Apply `lift` to every line of a block of `dims` dimensions, one axis at +/// a time from the slowest to the fastest, as zfp's `inv_xform` does. +fn xform(p: &mut [T], dims: usize, lift: impl Fn(&mut [T], usize, usize)) { + let size = 1usize << (2 * dims); + for axis in (0..dims).rev() { + let s = 1usize << (2 * axis); + for i in 0..size { + if (i / s).is_multiple_of(4) { + lift(p, i, s); + } + } + } +} + +/// The integer codec (`decode.c`, `revdecode.c`) for 32- or 64-bit +/// integers. +macro_rules! int_codec { + ($m:ident, $int:ty, $uint:ty, $nbmask:expr, $pbits:expr) => { + mod $m { + use super::{Bits, FormatError, perm, with_maxbits, xform}; + + pub(super) type Int = $int; + type UInt = $uint; + const INTPREC: u32 = <$uint>::BITS; + const NBMASK: UInt = $nbmask; + const PBITS: u32 = $pbits; + + fn uint2int(x: UInt) -> Int { + (x ^ NBMASK).wrapping_sub(NBMASK) as Int + } + + fn inv_lift(p: &mut [Int], at: usize, s: usize) { + let (mut x, mut y, mut z, mut w) = (p[at], p[at + s], p[at + 2 * s], p[at + 3 * s]); + y = y.wrapping_add(w >> 1); + w = w.wrapping_sub(y >> 1); + y = y.wrapping_add(w); + w = w.wrapping_shl(1); + w = w.wrapping_sub(y); + z = z.wrapping_add(x); + x = x.wrapping_shl(1); + x = x.wrapping_sub(z); + y = y.wrapping_add(z); + z = z.wrapping_shl(1); + z = z.wrapping_sub(y); + w = w.wrapping_add(x); + x = x.wrapping_shl(1); + x = x.wrapping_sub(w); + p[at] = x; + p[at + s] = y; + p[at + 2 * s] = z; + p[at + 3 * s] = w; + } + + fn rev_inv_lift(p: &mut [Int], at: usize, s: usize) { + let (x, mut y, mut z, mut w) = (p[at], p[at + s], p[at + 2 * s], p[at + 3 * s]); + w = w.wrapping_add(z); + z = z.wrapping_add(y); + w = w.wrapping_add(z); + y = y.wrapping_add(x); + z = z.wrapping_add(y); + w = w.wrapping_add(z); + p[at + s] = y; + p[at + 2 * s] = z; + p[at + 3 * s] = w; + } + + fn kmin(maxprec: u32) -> u32 { + INTPREC.saturating_sub(maxprec) + } + + /// `decode_few_ints` (size <= 64) and `decode_many_ints`: bit + /// planes within a budget of `maxbits` bits. + fn decode_ints_rate( + s: &mut Bits<'_>, + maxbits: u32, + maxprec: u32, + data: &mut [UInt], + ) -> Result { + let size = data.len(); + let kmin = kmin(maxprec); + let mut bits = maxbits; + data.fill(0); + let mut k = INTPREC; + let mut n = 0usize; + while bits != 0 && k > kmin { + k -= 1; + // step 1: the first n bits of bit plane k + let m = n.min(bits as usize); + bits -= m as u32; + if size <= 64 { + let mut x = s.bits(m as u32)?; + // step 2: unary run-length decode the rest + while bits != 0 && n < size { + bits -= 1; + if s.bit()? != 0 { + while bits != 0 && n < size - 1 { + bits -= 1; + if s.bit()? != 0 { + break; + } + n += 1; + } + x += 1u64 << n; + } else { + break; + } + n += 1; + } + // step 3: deposit the bit plane + let mut i = 0; + while x != 0 { + data[i] += ((x & 1) as UInt) << k; + x >>= 1; + i += 1; + } + } else { + for v in data.iter_mut().take(m) { + if s.bit()? != 0 { + *v += (1 as UInt) << k; + } + } + while bits != 0 && n < size { + bits -= 1; + if s.bit()? != 0 { + while bits != 0 && n < size - 1 { + bits -= 1; + if s.bit()? != 0 { + break; + } + n += 1; + } + data[n] += (1 as UInt) << k; + } else { + break; + } + n += 1; + } + } + } + Ok(maxbits - bits) + } + + /// `decode_few_ints_prec` and `decode_many_ints_prec`: whole bit + /// planes, no bit budget. + fn decode_ints_prec( + s: &mut Bits<'_>, + maxprec: u32, + data: &mut [UInt], + ) -> Result { + let size = data.len(); + let start = s.tell(); + let kmin = kmin(maxprec); + data.fill(0); + let mut k = INTPREC; + let mut n = 0usize; + while k > kmin { + k -= 1; + if size <= 64 { + let mut x = s.bits(n as u32)?; + while n < size && s.bit()? != 0 { + while n < size - 1 && s.bit()? == 0 { + n += 1; + } + x += 1u64 << n; + n += 1; + } + let mut i = 0; + while x != 0 { + data[i] += ((x & 1) as UInt) << k; + x >>= 1; + i += 1; + } + } else { + for v in data.iter_mut().take(n) { + if s.bit()? != 0 { + *v += (1 as UInt) << k; + } + } + while n < size && s.bit()? != 0 { + while n < size - 1 && s.bit()? == 0 { + n += 1; + } + data[n] += (1 as UInt) << k; + n += 1; + } + } + } + Ok((s.tell() - start) as u32) + } + + /// `decode_ints`. + fn decode_ints( + s: &mut Bits<'_>, + maxbits: u32, + maxprec: u32, + data: &mut [UInt], + ) -> Result { + if with_maxbits(maxbits, maxprec, data.len() as u32) { + decode_ints_rate(s, maxbits, maxprec, data) + } else { + decode_ints_prec(s, maxprec, data) + } + } + + /// Undo the coefficient order and the negabinary mapping. + fn inv_order(ublock: &[UInt], iblock: &mut [Int], dims: usize) { + for (&u, &p) in ublock.iter().zip(perm(dims)) { + iblock[p as usize] = uint2int(u); + } + } + + /// `decode_block` for integers: a block of `4^dims` values in + /// `iblock`. + pub(super) fn decode_block( + s: &mut Bits<'_>, + dims: usize, + minbits: u32, + maxbits: u32, + maxprec: u32, + iblock: &mut [Int], + ) -> Result { + let mut ublock = [0 as UInt; 256]; + let ublock = &mut ublock[..iblock.len()]; + let mut bits = decode_ints(s, maxbits, maxprec, ublock)?; + if bits < minbits { + s.skip(minbits - bits); + bits = minbits; + } + inv_order(ublock, iblock, dims); + xform(iblock, dims, inv_lift); + Ok(bits) + } + + /// `rev_decode_block` for integers. + pub(super) fn rev_decode_block( + s: &mut Bits<'_>, + dims: usize, + minbits: u32, + maxbits: u32, + iblock: &mut [Int], + ) -> Result { + let mut bits = PBITS; + let prec = s.bits(PBITS)? as u32 + 1; + let mut ublock = [0 as UInt; 256]; + let ublock = &mut ublock[..iblock.len()]; + bits += decode_ints(s, maxbits.wrapping_sub(bits), prec, ublock)?; + if bits < minbits { + s.skip(minbits - bits); + bits = minbits; + } + inv_order(ublock, iblock, dims); + xform(iblock, dims, rev_inv_lift); + Ok(bits) + } + } + }; +} + +int_codec!(int32, i32, u32, 0xaaaa_aaaa, 5); +int_codec!(int64, i64, u64, 0xaaaa_aaaa_aaaa_aaaa, 6); + +/// `ldexp(1, e)` for `f32`, exactly: 0 below the smallest subnormal (the +/// exact value, at most half of it, rounds to 0). +fn pow2_f32(e: i32) -> f32 { + if e >= -126 { + f32::from_bits(((e.min(128) + 127) as u32) << 23) + } else if e >= -149 { + f32::from_bits(1 << (e + 149)) + } else { + 0.0 + } +} + +/// `ldexp(1, e)` for `f64`, exactly, as [`pow2_f32`]. +fn pow2_f64(e: i32) -> f64 { + if e >= -1022 { + f64::from_bits(((e.min(1024) + 1023) as u64) << 52) + } else if e >= -1074 { + f64::from_bits(1 << (e + 1074)) + } else { + 0.0 + } +} + +/// The floating-point codec (`decodef.c`, `revdecodef.c`). +macro_rules! float_codec { + ($m:ident, $f:ty, $ints:ident, $ebits:expr, $tcmask:expr, $pow2:ident) => { + mod $m { + use super::{Bits, FormatError, Params, $ints, $pow2}; + use $ints::Int; + + const EBITS: u32 = $ebits; + const EBIAS: i32 = (1 << (EBITS - 1)) - 1; + const INTBITS: i32 = (core::mem::size_of::<$f>() * 8) as i32; + + /// `precision`: the bit planes to decode for a block whose + /// largest exponent is `maxexp`. + fn precision(maxexp: i32, maxprec: u32, minexp: i32, dims: usize) -> u32 { + let p = maxexp as i64 - minexp as i64 + 2 * dims as i64 + 2; + maxprec.min(p.max(0).min(u32::MAX as i64) as u32) + } + + /// `inv_cast`: scale the block's integers by its exponent. + fn inv_cast(iblock: &[Int], fblock: &mut [$f], emax: i32) { + let s = $pow2(emax - (INTBITS - 2)); + for (f, &i) in fblock.iter_mut().zip(iblock) { + *f = s * i as $f; + } + } + + fn decode_lossy( + s: &mut Bits<'_>, + p: &Params, + dims: usize, + fblock: &mut [$f], + ) -> Result<(), FormatError> { + let mut bits = 1u32; + if s.bit()? != 0 { + let mut iblock = [0 as Int; 256]; + let iblock = &mut iblock[..fblock.len()]; + bits += EBITS; + let emax = s.bits(EBITS)? as i32 - EBIAS; + let maxprec = precision(emax, p.maxprec, p.minexp, dims); + $ints::decode_block( + s, + dims, + p.minbits - bits.min(p.minbits), + p.maxbits.wrapping_sub(bits), + maxprec, + iblock, + )?; + inv_cast(iblock, fblock, emax); + } else { + fblock.fill(0.0); + if p.minbits > bits { + s.skip(p.minbits - bits); + } + } + Ok(()) + } + + fn decode_reversible( + s: &mut Bits<'_>, + p: &Params, + dims: usize, + fblock: &mut [$f], + ) -> Result<(), FormatError> { + let mut bits = 1u32; + if s.bit()? != 0 { + let mut iblock = [0 as Int; 256]; + let iblock = &mut iblock[..fblock.len()]; + bits += 1; + if s.bit()? != 0 { + // the values' bit patterns, as sign-magnitude + // integers stored in two's complement + $ints::rev_decode_block( + s, + dims, + p.minbits - bits.min(p.minbits), + p.maxbits.wrapping_sub(bits), + iblock, + )?; + for (f, &i) in fblock.iter_mut().zip(iblock.iter()) { + let i = if i < 0 { i ^ $tcmask } else { i }; + *f = <$f>::from_bits(i as _); + } + } else { + bits += EBITS; + let emax = s.bits(EBITS)? as i32 - EBIAS; + $ints::rev_decode_block( + s, + dims, + p.minbits - bits.min(p.minbits), + p.maxbits.wrapping_sub(bits), + iblock, + )?; + if emax != -EBIAS { + inv_cast(iblock, fblock, emax); + } else { + fblock.fill(0.0); + } + } + } else { + fblock.fill(0.0); + if p.minbits > bits { + s.skip(p.minbits - bits); + } + } + Ok(()) + } + + /// `zfp_decode_block` for floats: `4^dims` values. + pub(super) fn decode_block( + s: &mut Bits<'_>, + p: &Params, + dims: usize, + fblock: &mut [$f], + ) -> Result<(), FormatError> { + if p.reversible() { + decode_reversible(s, p, dims, fblock) + } else { + decode_lossy(s, p, dims, fblock) + } + } + } + }; +} + +float_codec!(float32, f32, int32, 8, 0x7fff_ffff, pow2_f32); +float_codec!(float64, f64, int64, 11, 0x7fff_ffff_ffff_ffff, pow2_f64); + +/// A decoded scalar, stored into the output in little- or big-endian order. +trait Scalar: Copy + Default { + const SIZE: usize; + fn store(self, out: &mut [u8], big_endian: bool); +} + +macro_rules! scalar { + ($t:ty) => { + impl Scalar for $t { + const SIZE: usize = core::mem::size_of::<$t>(); + fn store(self, out: &mut [u8], big_endian: bool) { + out.copy_from_slice(&if big_endian { + self.to_be_bytes() + } else { + self.to_le_bytes() + }); + } + } + }; +} +scalar!(i32); +scalar!(i64); +scalar!(f32); +scalar!(f64); + +/// `zfp_decompress`: decode the field's blocks in order (x fastest) and +/// keep the part of each inside the array. +fn decompress( + s: &mut Bits<'_>, + field: &Field, + out: &mut [u8], + big_endian: bool, + mut decode: impl FnMut(&mut Bits<'_>, &mut [T]) -> Result<(), FormatError>, +) -> Result<(), FormatError> { + let [nx, ny, nz, nw] = field.n; + let bsize = 1usize << (2 * field.dims); + let mut block = [T::default(); 256]; + let block = &mut block[..bsize]; + let ext = |n: usize, at: usize, d: usize| { + if d < field.dims { (n - at).min(4) } else { 1 } + }; + for w in (0..nw).step_by(4) { + for z in (0..nz).step_by(4) { + for y in (0..ny).step_by(4) { + for x in (0..nx).step_by(4) { + decode(s, block)?; + let (bx, by, bz, bw) = + (ext(nx, x, 0), ext(ny, y, 1), ext(nz, z, 2), ext(nw, w, 3)); + for j in 0..bw { + for k in 0..bz { + for i in 0..by { + let row = (((w + j) * nz + z + k) * ny + y + i) * nx + x; + let q = 64 * j + 16 * k + 4 * i; + for (h, &v) in block[q..q + bx].iter().enumerate() { + let at = (row + h) * T::SIZE; + v.store(&mut out[at..at + T::SIZE], big_endian); + } + } + } + } + } + } + } + } + Ok(()) +} + +/// Decode one chunk: the H5Z-ZFP filter's decompression direction. +pub(crate) fn zfp_decode(input: &[u8], ctx: &FilterContext<'_>) -> Result, FormatError> { + zfp_decompress(input, ctx.client_data(), ctx.max_output) +} + +/// Decode `input`, one chunk compressed by H5Z-ZFP with the filter +/// parameters `cd`. `max_output` is the chunk's size in bytes (0 when +/// unknown, in which case the ceiling is 256 MiB); a field of another size +/// is an error. +pub fn zfp_decompress(input: &[u8], cd: &[u32], max_output: usize) -> Result, FormatError> { + let h = parse_cd_values(cd)?; + let field = &h.field; + let size = field + .values() + .and_then(|n| n.checked_mul(field.ztype.size())) + .ok_or_else(|| err("field too large"))?; + if max_output != 0 && size != max_output { + return Err(err(&format!( + "the filter parameters describe {size} bytes, the chunk holds {max_output}" + ))); + } + if size > crate::filters::MAX_DECOMPRESS_SIZE.max(max_output) { + return Err(err("field too large")); + } + // Every block takes at least one bit of the stream. + let blocks = field.blocks().ok_or_else(|| err("field too large"))?; + if blocks as u64 > input.len() as u64 * 8 { + return Err(truncated()); + } + let mut out = vec![0u8; size]; + let mut s = Bits::new(input); + let p = h.params; + let dims = field.dims; + let be = h.swap; + match field.ztype { + ZType::Float => decompress::(&mut s, field, &mut out, be, |s, b| { + float32::decode_block(s, &p, dims, b) + })?, + ZType::Double => decompress::(&mut s, field, &mut out, be, |s, b| { + float64::decode_block(s, &p, dims, b) + })?, + ZType::Int32 => decompress::(&mut s, field, &mut out, be, |s, b| { + if p.reversible() { + int32::rev_decode_block(s, dims, p.minbits, p.maxbits, b) + } else { + int32::decode_block(s, dims, p.minbits, p.maxbits, p.maxprec, b) + } + .map(|_| ()) + })?, + ZType::Int64 => decompress::(&mut s, field, &mut out, be, |s, b| { + if p.reversible() { + int64::rev_decode_block(s, dims, p.minbits, p.maxbits, b) + } else { + int64::decode_block(s, dims, p.minbits, p.maxbits, p.maxprec, b) + } + .map(|_| ()) + })?, + } + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn modes_round_trip() { + // Fixed rate 16 bits/value in 2-D: 256 bits per block. + let p = Params::from_mode(255).unwrap(); + assert_eq!((p.minbits, p.maxbits), (256, 256)); + assert_eq!(p.mode(), 255); + // Reversible. + let p = Params::from_mode(2176).unwrap(); + assert!(p.reversible()); + assert_eq!(p.mode(), 2176); + // Fixed accuracy 2^-10. + let p = Params::from_mode((-10 - ZFP_MIN_EXP) as u64 + 2177).unwrap(); + assert_eq!(p.minexp, -10); + // An expert maxbits above ZFP_MAX_BITS is fixed precision to + // H5Z-ZFP, which hands zfp the standard mode's parameters. + let long = Params::new(1, 20000, 20, ZFP_MIN_EXP).unwrap().mode(); + assert_eq!(long, 2048 + 19); + let p = Params { + minbits: 3, + maxbits: 20000, + maxprec: 20, + minexp: -3, + }; + assert_eq!(Params::from_mode(p.mode()), Some(p)); + } + + #[test] + fn pow2_is_exact() { + // 2^e by repeated doubling or halving, each step exact down to the + // smallest subnormal (the next halving rounds to 0, as ldexp does). + let pow2 = |e: i32| { + let mut v = 1f64; + for _ in 0..e.unsigned_abs() { + v = if e < 0 { v / 2.0 } else { v * 2.0 }; + } + v + }; + for e in -160..128 { + let want = if e < -149 { 0.0 } else { pow2(e) }; + assert_eq!(pow2_f32(e) as f64, want, "{e}"); + } + for e in -1080..1000 { + assert_eq!(pow2_f64(e), pow2(e), "{e}"); + } + } + + /// The filter parameters of the conformance corpus's `h5ex_d_zfp.h5` + /// (H5Z-ZFP 1.0.1, zfp 0.5.5; no codec version in `cd_values[0]`): a + /// 2-D float field of 4x8 chunks. + #[test] + fn corpus_header_parses() { + let cd = [ + 5570817u32, 91252346, 805306486, 4293918720, 3767009280, 493487, + ]; + let h = parse_cd_values(&cd).unwrap(); + assert!(!h.swap); + assert_eq!(h.field.ztype, ZType::Float); + assert_eq!(h.field.dims, 2); + assert_eq!(h.field.n[..2], [8, 4]); + } + + #[test] + fn hostile_headers_are_errors() { + assert!(parse_cd_values(&[]).is_err()); + assert!(parse_cd_values(&[0]).is_err()); + assert!(parse_cd_values(&[0, 0x0570_667a]).is_err()); + assert!(parse_cd_values(&[0; 8]).is_err()); + // A newer codec. + let mut cd = [ + 5570817u32, 91252346, 805306486, 4293918720, 3767009280, 493487, + ]; + cd[0] = (0x1100 << 16) | (6 << 12) | 0x111; + assert!(parse_cd_values(&cd).is_err()); + } +} diff --git a/crates/clawhdf5-format/src/lib.rs b/crates/clawhdf5-format/src/lib.rs index a197c9f..7042a13 100644 --- a/crates/clawhdf5-format/src/lib.rs +++ b/crates/clawhdf5-format/src/lib.rs @@ -94,6 +94,8 @@ mod filters_bzip2; #[cfg(feature = "lzf")] pub mod filters_lzf; mod filters_szip; +#[cfg(feature = "zfp")] +pub mod filters_zfp; pub mod fixed_array; pub mod float16; pub mod fractal_heap; diff --git a/crates/clawhdf5-format/tests/zfp_alloc_bounds.rs b/crates/clawhdf5-format/tests/zfp_alloc_bounds.rs new file mode 100644 index 0000000..90119a6 --- /dev/null +++ b/crates/clawhdf5-format/tests/zfp_alloc_bounds.rs @@ -0,0 +1,274 @@ +//! Crafted ZFP filter parameters and streams cannot make the decoder panic +//! or allocate out of proportion to the chunk it decodes. +//! +//! The field size comes from the filter's `cd_values` (up to 2^48 values), +//! not from the chunk: the decoder allocates the output only when it matches +//! the chunk's size (or, when that is unknown, is within the 256 MiB +//! ceiling), and only when the stream is long enough to hold a bit per +//! block. Peak heap use is measured with a counting global allocator; the +//! tests share it, so each holds `SERIAL` for its whole run. +#![cfg(feature = "zfp")] + +use std::alloc::{GlobalAlloc, Layout, System}; +use std::sync::Mutex; +use std::sync::atomic::{AtomicUsize, Ordering}; + +use clawhdf5_format::filters_zfp::zfp_decompress; + +struct Counting; + +static CURRENT: AtomicUsize = AtomicUsize::new(0); +static PEAK: AtomicUsize = AtomicUsize::new(0); +static SERIAL: Mutex<()> = Mutex::new(()); + +unsafe impl GlobalAlloc for Counting { + unsafe fn alloc(&self, layout: Layout) -> *mut u8 { + let p = unsafe { System.alloc(layout) }; + if !p.is_null() { + let now = CURRENT.fetch_add(layout.size(), Ordering::Relaxed) + layout.size(); + PEAK.fetch_max(now, Ordering::Relaxed); + } + p + } + + unsafe fn alloc_zeroed(&self, layout: Layout) -> *mut u8 { + let p = unsafe { System.alloc_zeroed(layout) }; + if !p.is_null() { + let now = CURRENT.fetch_add(layout.size(), Ordering::Relaxed) + layout.size(); + PEAK.fetch_max(now, Ordering::Relaxed); + } + p + } + + unsafe fn dealloc(&self, ptr: *mut u8, layout: Layout) { + unsafe { System.dealloc(ptr, layout) }; + CURRENT.fetch_sub(layout.size(), Ordering::Relaxed); + } +} + +#[global_allocator] +static ALLOC: Counting = Counting; + +/// Bytes allocated at the peak of `f`, above what was live when it started. +fn peak_during(f: impl FnOnce() -> T) -> (T, usize) { + let base = CURRENT.load(Ordering::Relaxed); + PEAK.store(base, Ordering::Relaxed); + let out = f(); + (out, PEAK.load(Ordering::Relaxed).saturating_sub(base)) +} + +fn lock() -> std::sync::MutexGuard<'static, ()> { + SERIAL.lock().unwrap_or_else(|e| e.into_inner()) +} + +/// What decoding may hold: the output (at most the limit, and at most a +/// 4-D block of doubles, 2 KiB, per bit of input), and a little more. +fn bound(max_output: usize, input: &[u8]) -> usize { + let limit = if max_output == 0 { + 256 << 20 + } else { + max_output + }; + limit.min(input.len() * 8 * 2048) + 4096 +} + +/// An LSB-first bit writer. +#[derive(Default)] +struct Bits { + v: Vec, + n: usize, +} + +impl Bits { + fn put(&mut self, x: u64, bits: usize) { + for i in 0..bits { + if self.n.is_multiple_of(8) { + self.v.push(0); + } + if (x >> i) & 1 == 1 { + *self.v.last_mut().unwrap() |= 1 << (self.n % 8); + } + self.n += 1; + } + } +} + +/// H5Z-ZFP `cd_values`: a version word and a zfp header for a field of +/// `ztype` (0 int32, 1 int64, 2 float, 3 double) and sizes `n` (fastest +/// first), with `mode` (12 bits, or 64 when `long`). +fn cd_values(ztype: u64, n: &[u64], mode: u64, long: bool) -> Vec { + let mut b = Bits::default(); + for c in b"zfp" { + b.put(*c as u64, 8); + } + b.put(5, 8); + let dims = n.len(); + let mut meta = 0u64; + let bits = [48, 24, 16, 12][dims - 1]; + for &v in n.iter().rev() { + meta = (meta << bits) + v - 1; + } + meta = (meta << 2) + dims as u64 - 1; + meta = (meta << 2) + ztype; + b.put(meta, 52); + b.put(mode, if long { 64 } else { 12 }); + b.v.resize(b.v.len().div_ceil(4) * 4, 0); + let mut cd = vec![0x1001_1111u32]; + cd.extend( + b.v.chunks(4) + .map(|w| u32::from_le_bytes(w.try_into().unwrap())), + ); + cd +} + +fn elem(ztype: u64) -> usize { + if ztype & 1 == 0 { 4 } else { 8 } +} + +/// A 1-D field of 2^32 doubles (32 GiB) in a 1-byte chunk: refused for +/// its size, with or without the chunk size known, before anything is +/// allocated. +#[test] +fn huge_fields_are_refused_without_allocating() { + let _g = lock(); + for (n, ztype) in [ + (vec![1u64 << 32], 3), + (vec![1 << 24, 1 << 24], 3), + (vec![4096; 4], 1), + ] { + let cd = cd_values(ztype, &n, 2176, false); + for max_output in [0usize, 1 << 20] { + let (r, peak) = peak_during(|| zfp_decompress(&[0xff], &cd, max_output)); + assert!(r.is_err(), "{n:?}: decoded {:?} bytes", r.map(|v| v.len())); + assert!(peak < 4096, "{n:?}: peak {peak} bytes"); + } + } +} + +/// A field of the chunk's size whose stream is too short for its blocks +/// is refused before the output is allocated. +#[test] +fn short_streams_are_refused_before_allocating() { + let _g = lock(); + let n = [1u64 << 18]; + let cd = cd_values(2, &n, 2176, false); + let size = (1 << 18) * 4; + let (r, peak) = peak_during(|| zfp_decompress(&[0u8; 100], &cd, size)); + assert!(r.is_err()); + assert!(peak < 4096, "peak {peak} bytes"); + // A stream with a bit per block: all-zero blocks, which decode. + let input = vec![0u8; (1 << 16) / 8]; + let out = zfp_decompress(&input, &cd, size).unwrap(); + assert_eq!(out, vec![0u8; size]); +} + +/// xorshift64*: deterministic, so a failure reproduces. +struct Rng(u64); + +impl Rng { + fn next(&mut self) -> u64 { + let mut x = self.0; + x ^= x >> 12; + x ^= x << 25; + x ^= x >> 27; + self.0 = x; + x.wrapping_mul(0x2545_F491_4F6C_DD1D) + } + + fn below(&mut self, n: u64) -> u64 { + self.next() % n.max(1) + } +} + +/// A mode word: one of the four short forms, or the 64-bit expert form +/// with parameters at and past their edges (maxbits below a float block's +/// exponent, minbits past maxbits, precision 0, minexp at the reversible +/// boundary). +fn mode(rng: &mut Rng) -> (u64, bool) { + match rng.below(6) { + 0 => (rng.below(2048), false), + 1 => (2048 + rng.below(128), false), + 2 => (2176, false), + 3 => (2177 + rng.below(4094 - 2177 + 1), false), + _ => { + fn pick(rng: &mut Rng, v: [u64; 6]) -> u64 { + v[rng.below(6) as usize] + } + let r = [rng.below(2000), rng.below(0x8000), rng.below(40)]; + let minbits = pick(rng, [0, 1, 11, r[0], r[1], 1]); + let r = [rng.below(0x8000), rng.below(40)]; + let maxbits = pick(rng, [minbits, minbits + r[1], 7, 11, r[0], 0x7fff]); + let maxprec = rng.below(0x80); + let r = [rng.below(0x8000), rng.below(400)]; + let minexp = pick( + rng, + [ + r[0], + 16495 - 1074, + 16495 - 1075, + 16495 + r[1] - 200, + 16495 - 1074, + r[0], + ], + ); + let m = ((((minexp << 7) + maxprec) << 15) + maxbits) << 15; + ((m + minbits) << 12 | 0xfff, true) + } + } +} + +/// Random fields, modes and streams (random bytes, runs of ones, mostly +/// zeros; of every length), decoded with the chunk size known and not, +/// and with header words that are random too. +#[test] +fn fuzzed_headers_and_streams_stay_within_the_allocation_bound() { + let _g = lock(); + let mut rng = Rng(0x2f9); + let mut decoded = 0; + for i in 0..20_000 { + let ztype = rng.below(4); + let dims = 1 + rng.below(4) as usize; + let max = [300, 40, 14, 8][dims - 1]; + let n: Vec = (0..dims).map(|_| 1 + rng.below(max)).collect(); + let (m, long) = mode(&mut rng); + let mut cd = cd_values(ztype, &n, m, long); + if rng.below(10) == 0 { + let at = rng.below(cd.len() as u64) as usize; + cd[at] ^= 1 << rng.below(32); + } + if rng.below(20) == 0 { + cd.truncate(rng.below(cd.len() as u64 + 1) as usize); + } + let len = match rng.below(4) { + 0 => rng.below(8), + 1 => rng.below(300), + _ => rng.below(20_000), + } as usize; + let input: Vec = match rng.below(3) { + 0 => (0..len).map(|_| rng.next() as u8).collect(), + 1 => (0..len) + .map(|_| [0, 0xff, rng.next() as u8][rng.below(3) as usize]) + .collect(), + _ => (0..len) + .map(|_| [0, 0, 0, 1, 0x80, rng.next() as u8][rng.below(6) as usize]) + .collect(), + }; + let size = n.iter().product::() as usize * elem(ztype); + let max_output = if rng.below(4) == 0 { 0 } else { size }; + let (r, peak) = peak_during(|| zfp_decompress(&input, &cd, max_output)); + if let Ok(out) = &r { + decoded += 1; + if max_output != 0 { + assert_eq!(out.len(), max_output, "iteration {i}"); + } + } + assert!( + peak <= bound(max_output, &input), + "iteration {i}: peak {peak} bytes for {n:?} from {} bytes ({:?})", + input.len(), + r.map(|v| v.len()) + ); + } + // Most inputs are streams zfp decodes without running out. + assert!(decoded > 5_000, "only {decoded} decoded"); +}