[whahnsr]

Summary

Use four independent accumulators to construct bit planes in
encode_few_ints, reducing the serial accumulation dependency.

Each input bit occupies a distinct output position, preserving
the encoded bitstream.

Performance

Measured on an Intel Core i7-4770S (GCC 13.3, four OpenMP threads),
compressing a 256^3 float array at 8 bits/value.

Each measurement comprises 20 compressions; results are averages
of three measurements.

Metric	Original	Optimized	Reduction
Elapsed time	5.107 s	4.003 s	21.6%
CPU user time	17.39 s	13.14 s	24.4%

Profiling identified bit-plane extraction as the dominant hotspot,
accounting for 68% of sampled CPU cycles before optimization.

Validation

• Both standard regression tests pass.
• 22 bit-for-bit comparisons against the original implementation pass.
• Tested 1D through 4D float and double arrays.
• Tested serial and four-thread OpenMP compression.
• Tested multiple fixed compression rates.

AI assistance

The optimization was developed with assistance from ChatGPT.
The author reviewed the change and performed the builds,
performance measurements, and correctness validation.

[lindstro]

@whahnsr Thank you for this contribution. Although I can envision some vectorization being enabled by this partial unrolling of the loop, I found that it actually consistently slowed down compression on a Mac M2 using Apple Clang 17.0.0. Running the speed test (see the examples directory), throughput dropped from 1971 MB/s to 1736 MB/s. Because the speed test uses just a single block of data, which could skew things like branch prediction, I also ran testzfp large float 3d, which gave similar results; throughput in MB/s dropped from {619.7, 446.5, 359.8} to {578.2, 416.4, 324.5} for rates {2, 8, 32} bits/value. These values are the median throughput of five runs.

It's clearly possible that we'd observe different performance behavior using other compilers and ISAs, but I'm quite reluctant to adopt a change that complicates the code without a clear and consistent performance benefit.

I will say that we are quite aware that bit extraction is a bottleneck and have devised some fairly sophisticated solutions, including one based on fast divide-and-conquer binary matrix transposition (think of the bit planes as columns of a binary matrix that is stored as rows). As currently written, there are up to $4^d p$ UInt data accesses, where $p \leq 64$ is the precision and $d$ the number of dimensions. Each such access extracts a single bit, which is very wasteful. Using fast matrix transposition, you can reduce the number of accesses to $4^d + d p$. You can also perform this transposition on demand using lazy evaluation, where bit planes are processed as soon as they are available. Particularly, you don't want to perform the full transposition when rate or precision is low as you won't need the last several rows (bit planes) of the transposed matrix.

The code is ugly, though, and looks something like this for 3D blocks of 64-bit integers:

#define swap(x, y, l) \
  do { \
    const uint64 m[] = { \
      0x5555555555555555ul, \
      0x3333333333333333ul, \
      0x0f0f0f0f0f0f0f0ful, \
      0x00ff00ff00ff00fful, \
      0x0000ffff0000fffful, \
      0x00000000fffffffful, \
    }; \
    uint s = 1u << (l); \
    uint64 t = ((x) ^ ((y) >> s)) & m[(l)]; \
    (x) ^= t; \
    (y) ^= t << s; \
  } while (0)

static void _t2(fwd_transpose, UInt, 3)(uint64* bp, const UInt* data)
{
  /* working space for 64 bit planes */
  uint64 a00, a01, a02, a03, a04, a05, a06, a07;
  uint64 a08, a09, a0a, a0b, a0c, a0d, a0e, a0f;
  uint64 a10, a11, a12, a13, a14, a15, a16, a17;
  uint64 a18, a19, a1a, a1b, a1c, a1d, a1e, a1f;
  uint64 a20, a21, a22, a23, a24, a25, a26, a27;
  uint64 a28, a29, a2a, a2b, a2c, a2d, a2e, a2f;
  uint64 a30, a31, a32, a33, a34, a35, a36, a37;
  uint64 a38, a39, a3a, a3b, a3c, a3d, a3e, a3f;

  /* copy 64x64 matrix from input */
  a00 = data[0x3f-0x00]; a01 = data[0x3f-0x01]; a02 = data[0x3f-0x02]; a03 = data[0x3f-0x03];
  a04 = data[0x3f-0x04]; a05 = data[0x3f-0x05]; a06 = data[0x3f-0x06]; a07 = data[0x3f-0x07];
  a08 = data[0x3f-0x08]; a09 = data[0x3f-0x09]; a0a = data[0x3f-0x0a]; a0b = data[0x3f-0x0b];
  a0c = data[0x3f-0x0c]; a0d = data[0x3f-0x0d]; a0e = data[0x3f-0x0e]; a0f = data[0x3f-0x0f];
  a10 = data[0x3f-0x10]; a11 = data[0x3f-0x11]; a12 = data[0x3f-0x12]; a13 = data[0x3f-0x13];
  a14 = data[0x3f-0x14]; a15 = data[0x3f-0x15]; a16 = data[0x3f-0x16]; a17 = data[0x3f-0x17];
  a18 = data[0x3f-0x18]; a19 = data[0x3f-0x19]; a1a = data[0x3f-0x1a]; a1b = data[0x3f-0x1b];
  a1c = data[0x3f-0x1c]; a1d = data[0x3f-0x1d]; a1e = data[0x3f-0x1e]; a1f = data[0x3f-0x1f];
  a20 = data[0x3f-0x20]; a21 = data[0x3f-0x21]; a22 = data[0x3f-0x22]; a23 = data[0x3f-0x23];
  a24 = data[0x3f-0x24]; a25 = data[0x3f-0x25]; a26 = data[0x3f-0x26]; a27 = data[0x3f-0x27];
  a28 = data[0x3f-0x28]; a29 = data[0x3f-0x29]; a2a = data[0x3f-0x2a]; a2b = data[0x3f-0x2b];
  a2c = data[0x3f-0x2c]; a2d = data[0x3f-0x2d]; a2e = data[0x3f-0x2e]; a2f = data[0x3f-0x2f];
  a30 = data[0x3f-0x30]; a31 = data[0x3f-0x31]; a32 = data[0x3f-0x32]; a33 = data[0x3f-0x33];
  a34 = data[0x3f-0x34]; a35 = data[0x3f-0x35]; a36 = data[0x3f-0x36]; a37 = data[0x3f-0x37];
  a38 = data[0x3f-0x38]; a39 = data[0x3f-0x39]; a3a = data[0x3f-0x3a]; a3b = data[0x3f-0x3b];
  a3c = data[0x3f-0x3c]; a3d = data[0x3f-0x3d]; a3e = data[0x3f-0x3e]; a3f = data[0x3f-0x3f];

  swap(a00, a20, 5);
  swap(a01, a21, 5);
  swap(a02, a22, 5);
  swap(a03, a23, 5);
  swap(a04, a24, 5);
  swap(a05, a25, 5);
  swap(a06, a26, 5);
  swap(a07, a27, 5);
  swap(a08, a28, 5);
  swap(a09, a29, 5);
  swap(a0a, a2a, 5);
  swap(a0b, a2b, 5);
  swap(a0c, a2c, 5);
  swap(a0d, a2d, 5);
  swap(a0e, a2e, 5);
  swap(a0f, a2f, 5);
  swap(a10, a30, 5);
  swap(a11, a31, 5);
  swap(a12, a32, 5);
  swap(a13, a33, 5);
  swap(a14, a34, 5);
  swap(a15, a35, 5);
  swap(a16, a36, 5);
  swap(a17, a37, 5);
  swap(a18, a38, 5);
  swap(a19, a39, 5);
  swap(a1a, a3a, 5);
  swap(a1b, a3b, 5);
  swap(a1c, a3c, 5);
  swap(a1d, a3d, 5);
  swap(a1e, a3e, 5);
  swap(a1f, a3f, 5);

  swap(a00, a10, 4);
  swap(a01, a11, 4);
  swap(a02, a12, 4);
  swap(a03, a13, 4);
  swap(a04, a14, 4);
  swap(a05, a15, 4);
  swap(a06, a16, 4);
  swap(a07, a17, 4);
  swap(a08, a18, 4);
  swap(a09, a19, 4);
  swap(a0a, a1a, 4);
  swap(a0b, a1b, 4);
  swap(a0c, a1c, 4);
  swap(a0d, a1d, 4);
  swap(a0e, a1e, 4);
  swap(a0f, a1f, 4);

  swap(a00, a08, 3);
  swap(a01, a09, 3);
  swap(a02, a0a, 3);
  swap(a03, a0b, 3);
  swap(a04, a0c, 3);
  swap(a05, a0d, 3);
  swap(a06, a0e, 3);
  swap(a07, a0f, 3);

  swap(a00, a04, 2);
  swap(a01, a05, 2);
  swap(a02, a06, 2);
  swap(a03, a07, 2);

  swap(a00, a02, 1);
  swap(a01, a03, 1);

  swap(a00, a01, 0);

  bp[0x00] = a00;
  bp[0x01] = a01;

  swap(a02, a03, 0);

  bp[0x02] = a02;
  bp[0x03] = a03;

  swap(a04, a06, 1);
  swap(a05, a07, 1);

  swap(a04, a05, 0);

  bp[0x04] = a04;
  bp[0x05] = a05;

  swap(a06, a07, 0);

  bp[0x06] = a06;
  bp[0x07] = a07;

  swap(a08, a0c, 2);
  swap(a09, a0d, 2);
  swap(a0a, a0e, 2);
  swap(a0b, a0f, 2);

  swap(a08, a0a, 1);
  swap(a09, a0b, 1);

  swap(a08, a09, 0);

  bp[0x08] = a08;
  bp[0x09] = a09;

  swap(a0a, a0b, 0);

  bp[0x0a] = a0a;
  bp[0x0b] = a0b;

  swap(a0c, a0e, 1);
  swap(a0d, a0f, 1);

  swap(a0c, a0d, 0);

  bp[0x0c] = a0c;
  bp[0x0d] = a0d;

  swap(a0e, a0f, 0);

  bp[0x0e] = a0e;
  bp[0x0f] = a0f;

  swap(a10, a18, 3);
  swap(a11, a19, 3);
  swap(a12, a1a, 3);
  swap(a13, a1b, 3);
  swap(a14, a1c, 3);
  swap(a15, a1d, 3);
  swap(a16, a1e, 3);
  swap(a17, a1f, 3);

  swap(a10, a14, 2);
  swap(a11, a15, 2);
  swap(a12, a16, 2);
  swap(a13, a17, 2);

  swap(a10, a12, 1);
  swap(a11, a13, 1);

  swap(a10, a11, 0);

  bp[0x10] = a10;
  bp[0x11] = a11;

  swap(a12, a13, 0);

  bp[0x12] = a12;
  bp[0x13] = a13;

  swap(a14, a16, 1);
  swap(a15, a17, 1);

  swap(a14, a15, 0);

  bp[0x14] = a14;
  bp[0x15] = a15;

  swap(a16, a17, 0);

  bp[0x16] = a16;
  bp[0x17] = a17;

  swap(a18, a1c, 2);
  swap(a19, a1d, 2);
  swap(a1a, a1e, 2);
  swap(a1b, a1f, 2);

  swap(a18, a1a, 1);
  swap(a19, a1b, 1);

  swap(a18, a19, 0);

  bp[0x18] = a18;
  bp[0x19] = a19;

  swap(a1a, a1b, 0);

  bp[0x1a] = a1a;
  bp[0x1b] = a1b;

  swap(a1c, a1e, 1);
  swap(a1d, a1f, 1);

  swap(a1c, a1d, 0);

  bp[0x1c] = a1c;
  bp[0x1d] = a1d;

  swap(a1e, a1f, 0);

  bp[0x1e] = a1e;
  bp[0x1f] = a1f;

  swap(a20, a30, 4);
  swap(a21, a31, 4);
  swap(a22, a32, 4);
  swap(a23, a33, 4);
  swap(a24, a34, 4);
  swap(a25, a35, 4);
  swap(a26, a36, 4);
  swap(a27, a37, 4);
  swap(a28, a38, 4);
  swap(a29, a39, 4);
  swap(a2a, a3a, 4);
  swap(a2b, a3b, 4);
  swap(a2c, a3c, 4);
  swap(a2d, a3d, 4);
  swap(a2e, a3e, 4);
  swap(a2f, a3f, 4);

  swap(a20, a28, 3);
  swap(a21, a29, 3);
  swap(a22, a2a, 3);
  swap(a23, a2b, 3);
  swap(a24, a2c, 3);
  swap(a25, a2d, 3);
  swap(a26, a2e, 3);
  swap(a27, a2f, 3);

  swap(a20, a24, 2);
  swap(a21, a25, 2);
  swap(a22, a26, 2);
  swap(a23, a27, 2);

  swap(a20, a22, 1);
  swap(a21, a23, 1);

  swap(a20, a21, 0);

  bp[0x20] = a20;
  bp[0x21] = a21;

  swap(a22, a23, 0);

  bp[0x22] = a22;
  bp[0x23] = a23;

  swap(a24, a26, 1);
  swap(a25, a27, 1);

  swap(a24, a25, 0);

  bp[0x24] = a24;
  bp[0x25] = a25;

  swap(a26, a27, 0);

  bp[0x26] = a26;
  bp[0x27] = a27;

  swap(a28, a2c, 2);
  swap(a29, a2d, 2);
  swap(a2a, a2e, 2);
  swap(a2b, a2f, 2);

  swap(a28, a2a, 1);
  swap(a29, a2b, 1);

  swap(a28, a29, 0);

  bp[0x28] = a28;
  bp[0x29] = a29;

  swap(a2a, a2b, 0);

  bp[0x2a] = a2a;
  bp[0x2b] = a2b;

  swap(a2c, a2e, 1);
  swap(a2d, a2f, 1);

  swap(a2c, a2d, 0);

  bp[0x2c] = a2c;
  bp[0x2d] = a2d;

  swap(a2e, a2f, 0);

  bp[0x2e] = a2e;
  bp[0x2f] = a2f;

  swap(a30, a38, 3);
  swap(a31, a39, 3);
  swap(a32, a3a, 3);
  swap(a33, a3b, 3);
  swap(a34, a3c, 3);
  swap(a35, a3d, 3);
  swap(a36, a3e, 3);
  swap(a37, a3f, 3);

  swap(a30, a34, 2);
  swap(a31, a35, 2);
  swap(a32, a36, 2);
  swap(a33, a37, 2);

  swap(a30, a32, 1);
  swap(a31, a33, 1);

  swap(a30, a31, 0);

  bp[0x30] = a30;
  bp[0x31] = a31;

  swap(a32, a33, 0);

  bp[0x32] = a32;
  bp[0x33] = a33;

  swap(a34, a36, 1);
  swap(a35, a37, 1);

  swap(a34, a35, 0);

  bp[0x34] = a34;
  bp[0x35] = a35;

  swap(a36, a37, 0);

  bp[0x36] = a36;
  bp[0x37] = a37;

  swap(a38, a3c, 2);
  swap(a39, a3d, 2);
  swap(a3a, a3e, 2);
  swap(a3b, a3f, 2);

  swap(a38, a3a, 1);
  swap(a39, a3b, 1);

  swap(a38, a39, 0);

  bp[0x38] = a38;
  bp[0x39] = a39;

  swap(a3a, a3b, 0);

  bp[0x3a] = a3a;
  bp[0x3b] = a3b;

  swap(a3c, a3e, 1);
  swap(a3d, a3f, 1);

  swap(a3c, a3d, 0);

  bp[0x3c] = a3c;
  bp[0x3d] = a3d;

  swap(a3e, a3f, 0);

  bp[0x3e] = a3e;
  bp[0x3f] = a3f;
}

Instead of storing the result in bp (for bit plane), one could call a function or macro to encode that bit plane, but that also requires keeping around enough state--essentially all local variables in the encode_few_ints() function--if calling another function, or executing the per bit plane loop body inline each time there is an assignment to bp. That becomes many hundreds of lines of code that no human could possibly understand or maintain. Though we've thought of providing this as an optional implementation when speed really matters.

[whahnsr]

Thanks, Peter, for the careful benchmarking and detailed technical explanation.

The regression on Apple M2, despite the improvement measured on my Intel/GCC system, makes a good case against adopting this change as-is, particularly given the added complexity.

I appreciate the explanation of the binary matrix transposition approach. It's an interesting direction for future investigation.

I'll close this PR for now. Thanks again for taking the time to evaluate it.