// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-only
//
// Adversarial differential test: BSLZ4DecoderGPU vs the CPU reference (JFJochDecompressHperfPtr).
//
// The device decoder is a hand-written LZ4 parser, so "it works on our files" is not enough - it has
// to agree with LZ4_decompress_safe on every construct a conforming encoder may emit, and on the
// block sizes foreign writers (DECTRIS, the stock bitshuffle HDF5 filter) use rather than only the
// 16 kB our own compressor emits. The cases below therefore drive the compressor into the corners
// (incompressible data for long literal runs and 255-extension chains, long repeats for overlapping
// matches at offset 1, engineered offsets across the offset-vs-matchlen branch boundary, sizes that
// put every rem%8 value in the verbatim tail) and, where the compressor's heuristics will not
// produce a wanted sequence at all, hand-build valid LZ4 blocks and verify them with the reference
// decoder before feeding them to the GPU.
//
// The statistics printed by the coverage report are the point: they PROVE which paths were reached
// rather than assuming the generated data happened to reach them.

#include <catch2/catch_all.hpp>
#include "../common/CUDAWrapper.h"

#ifdef JFJOCH_USE_CUDA

#include <atomic>
#include <chrono>
#include <cstdio>
#include <cstring>
#include <map>
#include <mutex>
#include <random>
#include <sstream>
#include <thread>
#include <vector>

#include "../image_analysis/image_preprocessing/BSLZ4DecoderGPU.h"
#include "../compression/JFJochCompressor.h"
#include "../compression/JFJochDecompress.h"

namespace {

constexpr auto ALG = CompressionAlgorithm::BSHUF_LZ4;

size_t BlockElems(size_t elem_size) { return JFJochBitShuffleCompressor::BlockSize(ALG, elem_size); }
size_t DefaultBlockBytes() { return JFJochBitShuffleCompressor::DefaultBlockSizeBytes(ALG); }

CompressedImageMode ModeOf(size_t elem_size, bool is_signed) {
    switch (elem_size) {
        case 1: return is_signed ? CompressedImageMode::Int8  : CompressedImageMode::Uint8;
        case 2: return is_signed ? CompressedImageMode::Int16 : CompressedImageMode::Uint16;
        default: return is_signed ? CompressedImageMode::Int32 : CompressedImageMode::Uint32;
    }
}

void PutBE64(uint8_t *p, uint64_t v) { for (int i = 0; i < 8; i++) p[i] = (uint8_t)(v >> (8 * (7 - i))); }
void PutBE32(uint8_t *p, uint32_t v) { for (int i = 0; i < 4; i++) p[i] = (uint8_t)(v >> (8 * (3 - i))); }
uint64_t GetBE64(const uint8_t *p) { uint64_t v = 0; for (int i = 0; i < 8; i++) v = (v << 8) | p[i]; return v; }
uint32_t GetBE32(const uint8_t *p) { return ((uint32_t)p[0] << 24) | ((uint32_t)p[1] << 16) | ((uint32_t)p[2] << 8) | p[3]; }

// ---------------------------------------------------------------------------------------------
// LZ4 sequence statistics - so the test can PROVE which code paths it actually reached rather
// than assume the compressor produced them.
// ---------------------------------------------------------------------------------------------
struct Lz4Stats {
    size_t images = 0, blocks = 0, sequences = 0;
    size_t lit_ge15 = 0, lit_ext255 = 0, max_litlen = 0;
    size_t ml_ge19 = 0, ml_ext255 = 0, max_matchlen = 0;
    size_t overlap = 0, nonoverlap = 0;
    size_t off_eq_ml = 0, off_ml_plus1 = 0, off_ml_minus1 = 0;
    std::map<uint32_t, size_t> offsets_le64;
    uint32_t min_offset = 0xffffffffu;

    void Merge(const Lz4Stats &o) {
        images += o.images; blocks += o.blocks; sequences += o.sequences;
        lit_ge15 += o.lit_ge15; lit_ext255 += o.lit_ext255; max_litlen = std::max(max_litlen, o.max_litlen);
        ml_ge19 += o.ml_ge19; ml_ext255 += o.ml_ext255; max_matchlen = std::max(max_matchlen, o.max_matchlen);
        overlap += o.overlap; nonoverlap += o.nonoverlap;
        off_eq_ml += o.off_eq_ml; off_ml_plus1 += o.off_ml_plus1; off_ml_minus1 += o.off_ml_minus1;
        for (auto &kv : o.offsets_le64) offsets_le64[kv.first] += kv.second;
        min_offset = std::min(min_offset, o.min_offset);
    }
};

std::mutex g_stats_mtx;
Lz4Stats g_stats;

void AccumulateStats(const Lz4Stats &s) {
    std::lock_guard<std::mutex> lock(g_stats_mtx);
    g_stats.Merge(s);
}

// Walk one LZ4 block exactly as the kernel does and record what the sequences look like.
void ParseLz4Block(const uint8_t *p, size_t len, Lz4Stats &st) {
    const uint8_t *const end = p + len;
    while (p < end) {
        const uint32_t token = *p++;
        uint32_t litlen = token >> 4;
        if (litlen == 15) {
            st.lit_ge15++;
            uint32_t s, n_ext = 0;
            do { s = *p++; litlen += s; n_ext++; } while (s == 255 && p < end);
            if (n_ext > 1) st.lit_ext255++;
        }
        p += litlen;
        st.max_litlen = std::max<size_t>(st.max_litlen, litlen);
        if (p >= end) break;   // final sequence: literals only
        const uint32_t offset = (uint32_t)p[0] | ((uint32_t)p[1] << 8);
        p += 2;
        uint32_t ml = token & 0x0F;
        if (ml == 15) {
            st.ml_ge19++;
            uint32_t s, n_ext = 0;
            do { s = *p++; ml += s; n_ext++; } while (s == 255 && p < end);
            if (n_ext > 1) st.ml_ext255++;
        }
        ml += 4;
        st.sequences++;
        st.max_matchlen = std::max<size_t>(st.max_matchlen, ml);
        if (offset < ml) st.overlap++; else st.nonoverlap++;
        if (offset == ml) st.off_eq_ml++;
        if (offset == ml + 1) st.off_ml_plus1++;
        if (offset + 1 == ml) st.off_ml_minus1++;
        if (offset <= 64) st.offsets_le64[offset]++;
        st.min_offset = std::min(st.min_offset, offset);
    }
}

void ParseContainer(const uint8_t *c, size_t clen, size_t elem_size, Lz4Stats &st) {
    if (clen < 12) return;
    const size_t total = GetBE64(c);
    const size_t block_elems = GetBE32(c + 8) / elem_size;
    const size_t nelements = total / elem_size;
    const size_t nfull = nelements / block_elems;
    const size_t rem = nelements - nfull * block_elems;
    const size_t last = rem - rem % 8;
    size_t off = 12;
    st.images++;
    for (size_t i = 0; i < nfull + (last > 0 ? 1 : 0); i++) {
        if (off + 4 > clen) return;
        const uint32_t bl = GetBE32(c + off);
        off += 4;
        if (off + bl > clen) return;
        st.blocks++;
        ParseLz4Block(c + off, bl, st);
        off += bl;
    }
}

// ---------------------------------------------------------------------------------------------
// Byte-exact comparison with a diagnosable failure message.
// ---------------------------------------------------------------------------------------------
std::string DescribeMismatch(const uint8_t *expected, const uint8_t *got, size_t nbytes,
                             size_t block_bytes, size_t body_bytes, const std::string &label) {
    for (size_t i = 0; i < nbytes; i++) {
        if (expected[i] != got[i]) {
            size_t ndiff = 0;
            for (size_t j = i; j < nbytes; j++) if (expected[j] != got[j]) ndiff++;
            std::ostringstream os;
            os << label << ": FIRST DIFF at byte " << i
               << " expected 0x" << std::hex << (unsigned)expected[i]
               << " got 0x" << (unsigned)got[i] << std::dec
               << " | bitshuffle block " << (i / block_bytes)
               << " (offset within block " << (i % block_bytes) << " of " << block_bytes << ")"
               << (i >= body_bytes ? " [IN THE VERBATIM LEFTOVER TAIL]" : "")
               << " | total differing bytes " << ndiff << " of " << nbytes;
            return os.str();
        }
    }
    return {};
}

// ---------------------------------------------------------------------------------------------
// The core differential check: production compressor -> CPU decode and GPU decode -> must match.
// ---------------------------------------------------------------------------------------------
struct CaseResult {
    bool cpu_roundtrip_ok = false;
    bool gpu_matches_cpu = false;
    std::string message;
    size_t compressed_bytes = 0;
};

CaseResult RunCaseNoAssert(const std::string &label,
                           const std::vector<uint8_t> &image_bytes,
                           size_t elem_size, bool is_signed,
                           JFJochBitShuffleCompressor &compressor,
                           BSLZ4DecoderGPU &decoder,
                           CudaStream &stream,
                           bool collect_stats) {
    CaseResult r;
    const size_t nelements = image_bytes.size() / elem_size;
    const std::vector<uint8_t> compressed = compressor.Compress(image_bytes.data(), nelements, elem_size);
    r.compressed_bytes = compressed.size();

    if (collect_stats) {
        Lz4Stats st;
        ParseContainer(compressed.data(), compressed.size(), elem_size, st);
        AccumulateStats(st);
    }

    const CompressedImage image(compressed.data(), compressed.size(), nelements, 1,
                                ModeOf(elem_size, is_signed), ALG);

    // CPU reference - the exact call the host path makes.
    std::vector<uint8_t> cpu(image.GetUncompressedSize());
    image.GetUncompressed(cpu);
    r.cpu_roundtrip_ok = (cpu == image_bytes);

    CudaDevicePtr<uint8_t> gpu_out(image_bytes.size());
    // Poison the output so a byte the decoder never writes is caught rather than accidentally right.
    cudaMemset(gpu_out.get(), 0xA5, image_bytes.size());
    decoder.Decode(image, gpu_out.get());
    if (cudaStreamSynchronize(stream) != cudaSuccess) {
        r.message = label + ": cudaStreamSynchronize failed";
        return r;
    }
    std::vector<uint8_t> gpu(image_bytes.size());
    if (cudaMemcpy(gpu.data(), gpu_out.get(), image_bytes.size(), cudaMemcpyDeviceToHost) != cudaSuccess) {
        r.message = label + ": cudaMemcpy D2H failed";
        return r;
    }

    const size_t block_bytes = GetBE32(compressed.data() + 8);
    const size_t block_elems = block_bytes / elem_size;
    const size_t rem = nelements % block_elems;
    const size_t body_bytes = (nelements - rem % 8) * elem_size;

    const std::string msg = DescribeMismatch(cpu.data(), gpu.data(), image_bytes.size(),
                                             block_bytes, body_bytes, label);
    r.gpu_matches_cpu = msg.empty();
    if (!msg.empty()) r.message = msg;
    return r;
}

// Convenience wrapper that owns a fresh decoder/stream/compressor and asserts.
void RunCase(const std::string &label, const std::vector<uint8_t> &image_bytes,
             size_t elem_size, bool is_signed) {
    JFJochBitShuffleCompressor compressor(ALG);
    auto stream = std::make_shared<CudaStream>();
    BSLZ4DecoderGPU decoder(image_bytes.size(), stream);
    const auto r = RunCaseNoAssert(label, image_bytes, elem_size, is_signed,
                                   compressor, decoder, *stream, true);
    INFO(label << " (" << image_bytes.size() << " bytes -> " << r.compressed_bytes << " compressed)");
    if (!r.cpu_roundtrip_ok) FAIL(label << ": CPU reference did not round-trip (compressor/CPU bug, not GPU)");
    if (!r.gpu_matches_cpu) FAIL(r.message);
    SUCCEED();
}

// ---------------------------------------------------------------------------------------------
// Image generators (image domain).
// ---------------------------------------------------------------------------------------------
std::vector<uint8_t> GenIncompressible(size_t nbytes, uint32_t seed) {
    std::mt19937 rng(seed);
    std::vector<uint8_t> v(nbytes);
    for (auto &b : v) b = (uint8_t)(rng() & 0xFF);   // full 0..255 range
    return v;
}

std::vector<uint8_t> GenConstant(size_t nbytes, uint8_t val) { return std::vector<uint8_t>(nbytes, val); }

std::vector<uint8_t> GenRepeatedPattern(size_t nbytes, size_t period, uint32_t seed) {
    std::mt19937 rng(seed);
    std::vector<uint8_t> pat(period);
    for (auto &b : pat) b = (uint8_t)(rng() & 0xFF);
    std::vector<uint8_t> v(nbytes);
    for (size_t i = 0; i < nbytes; i++) v[i] = pat[i % period];
    return v;
}

// Incrementing element values - bit planes come out highly structured after bitshuffling.
std::vector<uint8_t> GenRamp(size_t nelements, size_t elem_size, uint64_t start, uint64_t step) {
    std::vector<uint8_t> v(nelements * elem_size);
    for (size_t i = 0; i < nelements; i++) {
        const uint64_t val = start + i * step;
        for (size_t j = 0; j < elem_size; j++) v[i * elem_size + j] = (uint8_t)(val >> (8 * j));
    }
    return v;
}

std::vector<uint8_t> GenDetectorLike(size_t nelements, size_t elem_size, uint32_t seed) {
    std::mt19937 rng(seed);
    std::vector<uint8_t> v(nelements * elem_size, 0);
    auto put = [&](size_t i, uint64_t val) {
        for (size_t j = 0; j < elem_size; j++) v[i * elem_size + j] = (uint8_t)(val >> (8 * j));
    };
    for (size_t i = nelements / 4; i < nelements / 2; i++) put(i, rng() % 7);
    for (size_t s = 0; s < 64 && nelements > 0; s++) {
        const size_t c = rng() % nelements;
        for (size_t d = 0; d < 9 && c + d < nelements; d++) put(c + d, 30000 + d);
    }
    for (size_t i = nelements * 3 / 4; i < nelements && i < nelements * 3 / 4 + 5000; i++) put(i, 42);
    return v;
}

// A grab-bag: random structure chosen per region, so a single case covers several regimes.
std::vector<uint8_t> GenMixed(size_t nelements, size_t elem_size, uint32_t seed) {
    std::mt19937 rng(seed);
    std::vector<uint8_t> v(nelements * elem_size, 0);
    size_t i = 0;
    while (i < nelements) {
        const size_t run = 1 + rng() % 4096;
        const size_t n = std::min(run, nelements - i);
        const int kind = rng() % 5;
        for (size_t k = 0; k < n; k++) {
            uint64_t val;
            switch (kind) {
                case 0: val = 0; break;
                case 1: val = rng(); break;                       // incompressible region
                case 2: val = 12345; break;                       // constant region
                case 3: val = k; break;                            // ramp region
                default: val = (rng() % 3); break;                 // low-entropy noise
            }
            for (size_t j = 0; j < elem_size; j++)
                v[(i + k) * elem_size + j] = (uint8_t)(val >> (8 * j));
        }
        i += n;
    }
    return v;
}

// ---------------------------------------------------------------------------------------------
// Shuffled-domain engineering: bitshuffle is a bit permutation, so to make LZ4 see a chosen byte
// pattern S we hand the compressor bitshuf_decode_block(S). encode(decode(S)) == S, which the
// helper verifies, so the LZ4 stage really does get the offsets we intend.
// ---------------------------------------------------------------------------------------------
std::vector<uint8_t> ImageFromShuffled(const std::vector<uint8_t> &shuffled, size_t nelements,
                                       size_t elem_size, bool verify) {
    std::vector<uint8_t> img(nelements * elem_size);
    std::vector<char> scratch(nelements * elem_size);
    REQUIRE(bitshuf_decode_block((char *)img.data(), (const char *)shuffled.data(),
                                 scratch.data(), nelements, elem_size) == 0);
    if (verify) {
        std::vector<uint8_t> re(nelements * elem_size);
        REQUIRE(bitshuf_encode_block((char *)re.data(), (const char *)img.data(),
                                     scratch.data(), nelements, elem_size) == 0);
        // Compare as a bool, not as vectors: Catch2 would otherwise stringify 16 kB of bytes.
        const bool encode_decode_is_identity = (re == shuffled);
        REQUIRE(encode_decode_is_identity);   // the construction is only meaningful if this holds
    }
    return img;
}

// A shuffled-domain buffer built from periodic regions of the requested periods, plus noise
// separators, so LZ4 emits matches at exactly those offsets.
std::vector<uint8_t> ShuffledWithPeriods(size_t nbytes, const std::vector<uint32_t> &periods,
                                         uint32_t seed) {
    std::mt19937 rng(seed);
    std::vector<uint8_t> s(nbytes);
    size_t i = 0;
    size_t pi = 0;
    while (i < nbytes) {
        // Noise separator so the periodic region starts fresh.
        const size_t nsep = std::min<size_t>(24 + rng() % 40, nbytes - i);
        for (size_t k = 0; k < nsep; k++) s[i + k] = (uint8_t)(rng() & 0xFF);
        i += nsep;
        if (i >= nbytes) break;
        const uint32_t P = periods[pi++ % periods.size()];
        // Region of period P; length varied so the resulting matchlen straddles the offset.
        const size_t reps = 1 + rng() % 12;
        const size_t len = std::min<size_t>((size_t)P * (1 + reps) + (rng() % 3), nbytes - i);
        for (size_t k = 0; k < len; k++)
            s[i + k] = (k < P) ? (uint8_t)(rng() & 0xFF) : s[i + k - P];
        i += len;
    }
    return s;
}

// Short NON-overlapping back-references at chosen distances: a chunk of L bytes repeated at
// distance P >= L, so offset >= matchlen.
std::vector<uint8_t> ShuffledShortMatches(size_t nbytes, const std::vector<uint32_t> &distances,
                                          uint32_t seed) {
    std::mt19937 rng(seed);
    std::vector<uint8_t> s(nbytes);
    for (auto &b : s) b = (uint8_t)(rng() & 0xFF);
    size_t i = 64;
    size_t di = 0;
    while (i + 128 < nbytes) {
        const uint32_t P = distances[di++ % distances.size()];
        if (P < 4) { i += 64; continue; }
        const uint32_t L = 4 + rng() % std::max<uint32_t>(1u, P - 3);   // 4 <= L <= P
        if (i < P) { i += 64; continue; }
        for (uint32_t k = 0; k < L && i + k < nbytes; k++) s[i + k] = s[i + k - P];
        i += L + 8 + rng() % 40;
    }
    return s;
}

// ---------------------------------------------------------------------------------------------
// Hand-built LZ4 blocks, so exact (offset, matchlen) pairs can be reached that the real
// compressor's heuristics never emit - in particular offset == matchlen +- 1.
// ---------------------------------------------------------------------------------------------
// Only the LAST sequence of an LZ4 block may be literals-only, so literals are staged and flushed
// by the match that follows them (or by Finish()).
class Lz4BlockBuilder {
    std::vector<uint8_t> pending;
public:
    std::vector<uint8_t> code;    // the encoded LZ4 block
    std::vector<uint8_t> plain;   // what it must decode to

    static void PutExt(std::vector<uint8_t> &v, uint32_t len) {
        uint32_t r = len - 15;
        while (r >= 255) { v.push_back(255); r -= 255; }
        v.push_back((uint8_t)r);
    }

    void AddLiterals(const std::vector<uint8_t> &l) { pending.insert(pending.end(), l.begin(), l.end()); }
    size_t OutputPos() const { return plain.size() + pending.size(); }
    size_t PendingLiterals() const { return pending.size(); }

    void Match(uint32_t offset, uint32_t matchlen) {
        const uint32_t litlen = (uint32_t)pending.size();
        const uint32_t mlcode = matchlen - 4;
        code.push_back((uint8_t)((std::min(litlen, 15u) << 4) | std::min(mlcode, 15u)));
        if (litlen >= 15) PutExt(code, litlen);
        code.insert(code.end(), pending.begin(), pending.end());
        plain.insert(plain.end(), pending.begin(), pending.end());
        pending.clear();
        code.push_back((uint8_t)(offset & 0xFF));
        code.push_back((uint8_t)(offset >> 8));
        if (mlcode >= 15) PutExt(code, mlcode);
        const size_t start = plain.size() - offset;
        for (uint32_t i = 0; i < matchlen; i++) {
            const uint8_t b = plain[start + i];   // may be a byte this very loop just wrote
            plain.push_back(b);
        }
    }

    // Every LZ4 block ends with a literals-only sequence; LZ4 also requires the last match to end
    // at least 12 bytes before the block end, so this tail must be >= 12 bytes.
    void Finish(const std::vector<uint8_t> &tail) {
        AddLiterals(tail);
        const uint32_t litlen = (uint32_t)pending.size();
        code.push_back((uint8_t)(std::min(litlen, 15u) << 4));
        if (litlen >= 15) PutExt(code, litlen);
        code.insert(code.end(), pending.begin(), pending.end());
        plain.insert(plain.end(), pending.begin(), pending.end());
        pending.clear();
    }
};

std::vector<uint8_t> MakeContainer(const std::vector<std::vector<uint8_t>> &block_codes,
                                   size_t total_bytes, size_t block_bytes) {
    std::vector<uint8_t> c(12);
    PutBE64(c.data(), total_bytes);
    PutBE32(c.data() + 8, (uint32_t)block_bytes);
    for (const auto &bc : block_codes) {
        uint8_t len[4];
        PutBE32(len, (uint32_t)bc.size());
        c.insert(c.end(), len, len + 4);
        c.insert(c.end(), bc.begin(), bc.end());
    }
    return c;
}

}  // namespace

// =============================================================================================
// (a) INCOMPRESSIBLE data - long pure-literal sequences, litlen >= 15 with 255-extension chains.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_Incompressible", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");
    size_t n_cases = 0;
    for (size_t es : {1u, 2u, 4u}) {
        for (bool sgn : {false, true}) {
            const size_t B = BlockElems(es);
            for (size_t nelem : {B / 2, B, B + 1, B + 7, B + 8, 3 * B, 3 * B + 13, 5 * B + 8 * 7 + 3}) {
                std::ostringstream lbl;
                lbl << "incompressible es=" << es << (sgn ? " signed" : " unsigned") << " nelem=" << nelem;
                RunCase(lbl.str(), GenIncompressible(nelem * es, (uint32_t)(1000 + n_cases)), es, sgn);
                n_cases++;
            }
        }
    }
    printf("[BSLZ4Fuzz] incompressible: %zu cases passed\n", n_cases);
}

// =============================================================================================
// (b) HIGHLY COMPRESSIBLE - very long matches, 255-extension match chains, offset==1 overlaps.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_HighlyCompressible", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");
    size_t n_cases = 0;
    for (size_t es : {1u, 2u, 4u}) {
        for (bool sgn : {false, true}) {
            const size_t B = BlockElems(es);
            for (size_t nelem : {B, B + 5, 4 * B, 4 * B + 8 * 11 + 6}) {
                std::ostringstream l1; l1 << "all-zero es=" << es << " n=" << nelem << (sgn ? " s" : " u");
                RunCase(l1.str(), GenConstant(nelem * es, 0x00), es, sgn);
                std::ostringstream l2; l2 << "all-0xFF es=" << es << " n=" << nelem << (sgn ? " s" : " u");
                RunCase(l2.str(), GenConstant(nelem * es, 0xFF), es, sgn);
                std::ostringstream l3; l3 << "all-0x5A es=" << es << " n=" << nelem << (sgn ? " s" : " u");
                RunCase(l3.str(), GenConstant(nelem * es, 0x5A), es, sgn);
                for (size_t period : {1u, 3u, 4u, 7u, 16u, 33u, 64u}) {
                    std::ostringstream l4;
                    l4 << "repeat-pattern p=" << period << " es=" << es << " n=" << nelem << (sgn ? " s" : " u");
                    RunCase(l4.str(), GenRepeatedPattern(nelem * es, period * es, (uint32_t)(2000 + n_cases)), es, sgn);
                    n_cases++;
                }
                n_cases += 3;
            }
        }
    }
    printf("[BSLZ4Fuzz] highly compressible: %zu cases passed\n", n_cases);
}

// =============================================================================================
// (c) MIXED / STEPPED, and engineered LZ4 offsets 1,2,3,4,5,15,16,17,31,32,33 in both the
//     offset >= matchlen and offset < matchlen branches.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_MixedAndEngineeredOffsets", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");
    const std::vector<uint32_t> targets = {1, 2, 3, 4, 5, 15, 16, 17, 31, 32, 33};
    size_t n_cases = 0;

    for (size_t es : {1u, 2u, 4u}) {
        for (bool sgn : {false, true}) {
            const size_t B = BlockElems(es);

            // Stepped data: after bitshuffling the bit planes are highly structured.
            for (uint64_t step : {1u, 2u, 3u, 256u, 65536u}) {
                std::ostringstream l; l << "ramp step=" << step << " es=" << es << (sgn ? " s" : " u");
                RunCase(l.str(), GenRamp(2 * B + 8 * 5 + 3, es, 0, step), es, sgn);
                n_cases++;
            }

            // Mixed regions.
            for (uint32_t seed = 0; seed < 4; seed++) {
                std::ostringstream l; l << "mixed es=" << es << " seed=" << seed << (sgn ? " s" : " u");
                RunCase(l.str(), GenMixed(3 * B + 8 * 3 + 5, es, 3000 + seed * 17 + (uint32_t)es), es, sgn);
                n_cases++;
            }

            // Detector-like.
            for (uint32_t seed = 0; seed < 3; seed++) {
                std::ostringstream l; l << "detector-like es=" << es << " seed=" << seed << (sgn ? " s" : " u");
                RunCase(l.str(), GenDetectorLike(4 * B + 8 * 9 + 7, es, 4000 + seed * 31 + (uint32_t)es), es, sgn);
                n_cases++;
            }

            // Engineered periodic offsets -> offset < matchlen (the overlapping branch).
            for (uint32_t seed = 0; seed < 4; seed++) {
                const size_t nelem = 3 * B;
                const auto shuffled = ShuffledWithPeriods(nelem * es, targets, 5000 + seed * 7 + (uint32_t)es);
                // Per-block construction, because bitshuffle works per block.
                std::vector<uint8_t> img;
                img.reserve(nelem * es);
                for (size_t b = 0; b < nelem / B; b++) {
                    std::vector<uint8_t> sblk(shuffled.begin() + b * B * es, shuffled.begin() + (b + 1) * B * es);
                    const auto part = ImageFromShuffled(sblk, B, es, b == 0 && seed == 0);
                    img.insert(img.end(), part.begin(), part.end());
                }
                std::ostringstream l; l << "engineered-periodic es=" << es << " seed=" << seed << (sgn ? " s" : " u");
                RunCase(l.str(), img, es, sgn);
                n_cases++;
            }

            // Engineered short back-references -> offset >= matchlen (the non-overlapping branch).
            for (uint32_t seed = 0; seed < 4; seed++) {
                const size_t nelem = 3 * B;
                const auto shuffled = ShuffledShortMatches(nelem * es, targets, 6000 + seed * 11 + (uint32_t)es);
                std::vector<uint8_t> img;
                img.reserve(nelem * es);
                for (size_t b = 0; b < nelem / B; b++) {
                    std::vector<uint8_t> sblk(shuffled.begin() + b * B * es, shuffled.begin() + (b + 1) * B * es);
                    const auto part = ImageFromShuffled(sblk, B, es, false);
                    img.insert(img.end(), part.begin(), part.end());
                }
                std::ostringstream l; l << "engineered-short-match es=" << es << " seed=" << seed << (sgn ? " s" : " u");
                RunCase(l.str(), img, es, sgn);
                n_cases++;
            }
        }
    }
    printf("[BSLZ4Fuzz] mixed/engineered: %zu cases passed\n", n_cases);
}

// =============================================================================================
// (d) SIZE SWEEP - every interesting last-block / leftover-tail configuration.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_SizeSweep", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");
    size_t n_cases = 0;
    for (size_t es : {1u, 2u, 4u}) {
        const size_t B = BlockElems(es);
        INFO("elem_size " << es << " block " << B << " elements = " << B * es << " bytes");
        std::vector<size_t> sizes;
        sizes.push_back(8);                       // the smallest legal block
        sizes.push_back(16);
        sizes.push_back(B);                       // exactly one full block, no remainder, no leftover
        sizes.push_back(2 * B);                   // exact multiple: no partial block, no leftover
        sizes.push_back(7 * B);
        for (size_t k = 1; k <= 8; k++) sizes.push_back(B + k);       // 1..8 leftover elements
        for (size_t r = 0; r < 8; r++) sizes.push_back(B + 128 + r);  // rem % 8 == r with a real last block
        for (size_t r = 0; r < 8; r++) sizes.push_back(3 * B + 8 * 37 + r);
        sizes.push_back(B - 8);                   // one partial block only, smaller than a full block
        sizes.push_back(B - 1);                   // partial block + leftover
        sizes.push_back(B / 2 + 3);

        for (bool sgn : {false, true}) {
            for (size_t nelem : sizes) {
                for (int gen = 0; gen < 3; gen++) {
                    std::vector<uint8_t> img;
                    if (gen == 0) img = GenIncompressible(nelem * es, (uint32_t)(7000 + n_cases));
                    else if (gen == 1) img = GenDetectorLike(nelem, es, (uint32_t)(8000 + n_cases));
                    else img = GenConstant(nelem * es, 0x77);
                    std::ostringstream l;
                    l << "size-sweep es=" << es << " nelem=" << nelem << " gen=" << gen << (sgn ? " s" : " u");
                    RunCase(l.str(), img, es, sgn);
                    n_cases++;
                }
            }
        }
    }
    printf("[BSLZ4Fuzz] size sweep: %zu cases passed\n", n_cases);
}

// Images with fewer than BSHUF_BLOCKED_MULT elements contain no bitshuffle block at all - only the
// verbatim tail. Documented behaviour probe, not an assertion about correctness.
TEST_CASE("BSLZ4Fuzz_SubBlockSizes", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");
    JFJochBitShuffleCompressor compressor(ALG);
    auto stream = std::make_shared<CudaStream>();
    BSLZ4DecoderGPU decoder(4096, stream);
    for (size_t nelem = 1; nelem <= 7; nelem++) {
        const auto img = GenIncompressible(nelem * 4, (uint32_t)(9000 + nelem));
        const auto compressed = compressor.Compress(img.data(), nelem, 4);
        const CompressedImage image(compressed.data(), compressed.size(), nelem, 1,
                                    CompressedImageMode::Uint32, ALG);
        std::vector<uint8_t> cpu;
        image.GetUncompressed(cpu);
        CHECK(cpu == img);   // CPU handles a block-less image
        // Fewer than 8 elements is entirely the verbatim tail - there is no bitshuffle block at all.
        // The device decoder must still produce the image rather than decline it, so that the caller
        // does not need a special case the host route would have handled.
        CudaDevicePtr<uint8_t> gpu_out(nelem * 4);
        cudaMemset(gpu_out.get(), 0xA5, nelem * 4);
        decoder.Decode(image, gpu_out.get());
        REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess);
        std::vector<uint8_t> gpu(nelem * 4);
        REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), nelem * 4, cudaMemcpyDeviceToHost) == cudaSuccess);
        CHECK(gpu == img);
    }
}

// =============================================================================================
// (f) LARGE REALISTIC FRAME + speed.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_LargeFrame", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");
    const size_t W = 4148, H = 4362, N = W * H, ES = 4;
    const auto img = GenDetectorLike(N, ES, 12345);

    JFJochBitShuffleCompressor compressor(ALG);
    const auto compressed = compressor.Compress(img.data(), N, ES);
    Lz4Stats st;
    ParseContainer(compressed.data(), compressed.size(), ES, st);
    AccumulateStats(st);

    const CompressedImage image(compressed.data(), compressed.size(), N, 1, CompressedImageMode::Uint32, ALG);
    std::vector<uint8_t> cpu;
    auto t0 = std::chrono::steady_clock::now();
    image.GetUncompressed(cpu);
    auto t1 = std::chrono::steady_clock::now();
    { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); }

    auto stream = std::make_shared<CudaStream>();
    BSLZ4DecoderGPU decoder(N * ES, stream);
    CudaDevicePtr<uint8_t> gpu_out(N * ES);
    decoder.Decode(image, gpu_out.get());
    REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess);
    std::vector<uint8_t> gpu(N * ES);
    REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), N * ES, cudaMemcpyDeviceToHost) == cudaSuccess);
    const auto msg = DescribeMismatch(cpu.data(), gpu.data(), N * ES, 16384, N * ES, "large frame");
    if (!msg.empty()) FAIL(msg);

    // Timed loop (upload + both kernels), after a warm-up.
    const int reps = 20;
    auto t2 = std::chrono::steady_clock::now();
    for (int i = 0; i < reps; i++) decoder.Decode(image, gpu_out.get());
    REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess);
    auto t3 = std::chrono::steady_clock::now();

    const double cpu_ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
    const double gpu_ms = std::chrono::duration<double, std::milli>(t3 - t2).count() / reps;
    printf("[BSLZ4Fuzz] large frame %zux%zu uint32 = %.1f MB, compressed %.1f MB (%.2fx), %zu blocks\n",
           W, H, N * ES / 1e6, compressed.size() / 1e6, (double)(N * ES) / compressed.size(), st.blocks);
    printf("[BSLZ4Fuzz] CPU decode (1 thread) %.2f ms   GPU decode %.3f ms   speedup %.1fx\n",
           cpu_ms, gpu_ms, cpu_ms / gpu_ms);
    SUCCEED();
}

// =============================================================================================
// (g) BUFFER REUSE - one decoder instance, many images of different sizes in sequence.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_DecoderReuse", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");
    const size_t big_elems = 2048 * 2048;          // 8 MB at uint32
    JFJochBitShuffleCompressor compressor(ALG);
    auto stream = std::make_shared<CudaStream>();
    BSLZ4DecoderGPU decoder(big_elems * 4, stream);

    struct Item { size_t nelem; size_t es; int gen; };
    const std::vector<Item> seq = {
        {big_elems, 4, 1},     // large
        {1000, 4, 0},          // tiny right after a large one
        {big_elems, 4, 0},     // large again
        {8, 4, 0},             // minimum
        {big_elems / 3 + 5, 2, 1},
        {big_elems, 4, 2},     // all-constant large
        {4096 + 3, 1, 0},
        {big_elems, 4, 1},
        {17, 2, 0},
        {big_elems * 2, 2, 1}, // same byte count, different element size
    };
    size_t idx = 0;
    for (const auto &it : seq) {
        std::vector<uint8_t> img;
        if (it.gen == 0) img = GenIncompressible(it.nelem * it.es, (uint32_t)(11000 + idx));
        else if (it.gen == 1) img = GenDetectorLike(it.nelem, it.es, (uint32_t)(12000 + idx));
        else img = GenConstant(it.nelem * it.es, 0x3C);
        std::ostringstream l; l << "reuse[" << idx << "] es=" << it.es << " nelem=" << it.nelem;
        const auto r = RunCaseNoAssert(l.str(), img, it.es, false, compressor, decoder, *stream, true);
        INFO(l.str());
        if (!r.cpu_roundtrip_ok) FAIL(l.str() << ": CPU reference did not round-trip");
        if (!r.gpu_matches_cpu) FAIL(r.message);
        idx++;
    }
    printf("[BSLZ4Fuzz] decoder reuse: %zu sequential images through ONE decoder passed\n", seq.size());
    SUCCEED();
}

// =============================================================================================
// (h) CONCURRENCY.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_ConcurrentSeparateDecoders", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");
    const int nthreads = 4;
    const int iters = 12;
    const size_t nelem = 1024 * 1024;

    std::vector<std::string> failures(nthreads);
    std::vector<int> passed(nthreads, 0);

    std::vector<std::thread> threads;
    for (int t = 0; t < nthreads; t++) {
        threads.emplace_back([&, t]() {
            try {
                JFJochBitShuffleCompressor compressor(ALG);
                auto stream = std::make_shared<CudaStream>();
                BSLZ4DecoderGPU decoder(nelem * 4, stream);
                for (int i = 0; i < iters; i++) {
                    const size_t n = nelem - (size_t)(t * 977 + i * 13);
                    std::vector<uint8_t> img = (i % 2)
                        ? GenDetectorLike(n, 4, (uint32_t)(20000 + t * 100 + i))
                        : GenIncompressible(n * 4, (uint32_t)(21000 + t * 100 + i));
                    std::ostringstream l; l << "thread " << t << " iter " << i << " nelem=" << n;
                    const auto r = RunCaseNoAssert(l.str(), img, 4, false, compressor, decoder, *stream, false);
                    if (!r.cpu_roundtrip_ok) { failures[t] = l.str() + ": CPU round-trip failed"; return; }
                    if (!r.gpu_matches_cpu) { failures[t] = r.message; return; }
                    passed[t]++;
                }
            } catch (const std::exception &e) {
                failures[t] = std::string("thread threw: ") + e.what();
            }
        });
    }
    for (auto &th : threads) th.join();

    int total = 0;
    for (int t = 0; t < nthreads; t++) {
        total += passed[t];
        if (!failures[t].empty()) FAIL(failures[t]);
    }
    printf("[BSLZ4Fuzz] concurrency (own decoder + own stream per thread): %d/%d decodes byte-exact\n",
           total, nthreads * iters);
    SUCCEED();
}

// =============================================================================================
// Hand-built LZ4 sequences: exact (offset, matchlen) pairs across the branch boundary.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_SyntheticOffsetMatrix", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");

    // A bitshuffle block is 16384 bytes, so an LZ4 offset can never exceed that here even though the
    // format allows 65535.
    const std::vector<uint32_t> offsets = {1, 2, 3, 4, 5, 6, 7, 8, 15, 16, 17, 30, 31, 32, 33, 34,
                                           63, 64, 65, 127, 128, 129, 255, 256, 257, 1024, 4095,
                                           8192, 16000, 16300};
    const std::vector<uint32_t> elem_sizes = {1, 2, 4};

    size_t n_cases = 0, n_seq = 0, n_skipped = 0;
    size_t n_overlap = 0, n_nonoverlap = 0, n_boundary = 0;
    for (uint32_t es : elem_sizes) {
        const size_t B = BlockElems(es);          // elements per block
        const size_t block_bytes = B * es;        // == 16384

        auto stream = std::make_shared<CudaStream>();
        BSLZ4DecoderGPU decoder(block_bytes, stream);
        CudaDevicePtr<uint8_t> gpu_out(block_bytes);

        for (uint32_t off : offsets) {
            // For each offset, a set of match lengths that straddles the offset >= matchlen boundary.
            std::vector<uint32_t> mls;
            for (int d : {-2, -1, 0, 1, 2}) {
                const long v = (long)off + d;
                if (v >= 4) mls.push_back((uint32_t)v);
            }
            mls.push_back(4);
            mls.push_back(18);                    // largest matchlen with no extension byte
            mls.push_back(19);                    // first matchlen needing an extension byte
            mls.push_back(4 + 15 + 255);          // exactly one full 255 extension byte
            mls.push_back(4 + 15 + 255 + 1);
            mls.push_back(4 + 15 + 255 * 2 + 7);  // a 255-extension CHAIN
            std::sort(mls.begin(), mls.end());
            mls.erase(std::unique(mls.begin(), mls.end()), mls.end());

            std::mt19937 rng(40000 + off * 31 + es);
            Lz4BlockBuilder bb;
            // Prime the output with enough literals that this offset is a legal back-reference.
            {
                std::vector<uint8_t> lits(std::max<size_t>(64, off + 64));
                for (auto &b : lits) b = (uint8_t)(rng() & 0xFF);
                bb.AddLiterals(lits);
            }
            if (bb.OutputPos() + 32 > block_bytes) { n_skipped++; continue; }
            size_t seq_in_block = 0;
            for (uint32_t ml : mls) {
                if (bb.OutputPos() + ml + 64 > block_bytes) break;
                if (off > bb.OutputPos()) continue;
                bb.Match(off, ml);
                n_seq++; seq_in_block++;
                // A few literals before the next sequence, so litlen varies too.
                const uint32_t litlen = (uint32_t)(rng() % 40);
                std::vector<uint8_t> lits(litlen);
                for (auto &b : lits) b = (uint8_t)(rng() & 0xFF);
                bb.AddLiterals(lits);
            }
            if (seq_in_block == 0) { n_skipped++; continue; }
            // LZ4 requires the block to end with a literal run and the last match to end at least
            // 12 bytes before the block end; pad out to exactly one full block.
            REQUIRE(bb.OutputPos() + 12 <= block_bytes);
            std::vector<uint8_t> tail(block_bytes - bb.OutputPos());
            for (auto &b : tail) b = (uint8_t)(rng() & 0xFF);
            bb.Finish(tail);
            REQUIRE(bb.plain.size() == block_bytes);

            // The hand-encoded stream must be a stream the real LZ4 decoder accepts, otherwise
            // this test would be measuring nonsense.
            std::vector<uint8_t> lz4_check(block_bytes);
            const int ret = LZ4_decompress_safe((const char *)bb.code.data(), (char *)lz4_check.data(),
                                                (int)bb.code.size(), (int)block_bytes);
            REQUIRE(ret == (int)block_bytes);
            { const bool lz4_ok = (lz4_check == bb.plain); REQUIRE(lz4_ok); }

            // Confirm the block really carries the offset under test, in both branches.
            Lz4Stats st;
            ParseLz4Block(bb.code.data(), bb.code.size(), st);
            REQUIRE(st.sequences == seq_in_block);
            REQUIRE(st.min_offset == off);
            n_overlap += st.overlap;
            n_nonoverlap += st.nonoverlap;
            n_boundary += st.off_eq_ml + st.off_ml_plus1 + st.off_ml_minus1;

            // bb.plain is the SHUFFLED image; the image itself is its bitshuffle inverse.
            const auto img = ImageFromShuffled(bb.plain, B, es, off == 1);

            const auto container = MakeContainer({bb.code}, block_bytes, block_bytes);
            const CompressedImage image(container.data(), container.size(), B, 1, ModeOf(es, false), ALG);
            REQUIRE(BSLZ4DecoderGPU::Supports(image));

            std::vector<uint8_t> cpu;
            image.GetUncompressed(cpu);
            { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); }

            cudaMemset(gpu_out.get(), 0xA5, block_bytes);
            decoder.Decode(image, gpu_out.get());
            REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess);
            std::vector<uint8_t> gpu(block_bytes);
            REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), block_bytes, cudaMemcpyDeviceToHost) == cudaSuccess);

            std::ostringstream l; l << "synthetic offset=" << off << " es=" << es;
            const auto msg = DescribeMismatch(cpu.data(), gpu.data(), block_bytes, block_bytes, block_bytes, l.str());
            if (!msg.empty()) FAIL(msg);
            n_cases++;
        }
    }
    printf("[BSLZ4Fuzz] synthetic hand-built LZ4: %zu blocks, %zu engineered (offset,matchlen) "
           "sequences passed (%zu offsets skipped as too large to prime in a 16 kB block)\n",
           n_cases, n_seq, n_skipped);
    printf("[BSLZ4Fuzz]   of those: %zu offset<matchlen (OVERLAP branch), %zu offset>=matchlen, "
           "%zu exactly on the offset==matchlen+-1 boundary\n", n_overlap, n_nonoverlap, n_boundary);
    SUCCEED();
}

// Very long literal runs with 255-extension chains, hand-built so the chain length is certain.
TEST_CASE("BSLZ4Fuzz_SyntheticLongLiterals", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");
    for (uint32_t es : {1u, 2u, 4u}) {
        const size_t B = BlockElems(es), block_bytes = B * es;
        auto stream = std::make_shared<CudaStream>();
        BSLZ4DecoderGPU decoder(block_bytes, stream);
        CudaDevicePtr<uint8_t> gpu_out(block_bytes);

        for (uint32_t litlen : {14u, 15u, 16u, 269u, 270u, 271u, 524u, 525u, 4000u, 8000u}) {
            std::mt19937 rng(50000 + litlen + es);
            Lz4BlockBuilder bb;
            // Prime and flush, so the long literal run below is a sequence of its own with exactly
            // the litlen under test and a legal back-reference after it.
            {
                std::vector<uint8_t> prime(512);
                for (auto &b : prime) b = (uint8_t)(rng() & 0xFF);
                bb.AddLiterals(prime);
                bb.Match(64, 64);
            }
            std::vector<uint8_t> lits(litlen);
            for (auto &b : lits) b = (uint8_t)(rng() & 0xFF);
            bb.AddLiterals(lits);
            REQUIRE(bb.PendingLiterals() == litlen);
            bb.Match(40, 300);           // this sequence carries exactly `litlen` literals
            REQUIRE(bb.OutputPos() + 12 <= block_bytes);
            std::vector<uint8_t> tail(block_bytes - bb.OutputPos());
            for (auto &b : tail) b = (uint8_t)(rng() & 0xFF);
            bb.Finish(tail);
            REQUIRE(bb.plain.size() == block_bytes);

            std::vector<uint8_t> chk(block_bytes);
            REQUIRE(LZ4_decompress_safe((const char *)bb.code.data(), (char *)chk.data(),
                                        (int)bb.code.size(), (int)block_bytes) == (int)block_bytes);
            { const bool lz4_ok = (chk == bb.plain); REQUIRE(lz4_ok); }

            // Confirm the encoder really emitted the literal run under test (and its 255-chain).
            Lz4Stats st;
            ParseLz4Block(bb.code.data(), bb.code.size(), st);
            REQUIRE(st.max_litlen >= litlen);

            const auto img = ImageFromShuffled(bb.plain, B, es, false);
            const auto container = MakeContainer({bb.code}, block_bytes, block_bytes);
            const CompressedImage image(container.data(), container.size(), B, 1, ModeOf(es, false), ALG);
            std::vector<uint8_t> cpu;
            image.GetUncompressed(cpu);
            { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); }

            cudaMemset(gpu_out.get(), 0xA5, block_bytes);
            decoder.Decode(image, gpu_out.get());
            REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess);
            std::vector<uint8_t> gpu(block_bytes);
            REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), block_bytes, cudaMemcpyDeviceToHost) == cudaSuccess);
            std::ostringstream l; l << "synthetic litlen=" << litlen << " es=" << es;
            const auto msg = DescribeMismatch(cpu.data(), gpu.data(), block_bytes, block_bytes, block_bytes, l.str());
            if (!msg.empty()) FAIL(msg);
        }
    }
    printf("[BSLZ4Fuzz] synthetic long-literal cases passed\n");
    SUCCEED();
}

// =============================================================================================
// FOREIGN BLOCK SIZES.
//
// Everything above goes through this repo's compressor, which always emits 16384-byte bitshuffle
// blocks. Real DECTRIS/EIGER files - and anything written by the stock bitshuffle HDF5 filter -
// use bshuf_default_block_size(), i.e. 8192 BYTES, half of that; the filter also lets a writer
// choose any multiple of 8. The GPU decoder reads the block size out of the container header, so
// it must cope with all of them, including block counts large enough to force the descriptor
// arrays to grow.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_ForeignBlockSizes", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");

    size_t n_ok = 0, n_refused = 0;
    for (uint32_t es : {1u, 2u, 4u}) {
        const size_t nelem = 300000;
        const auto img = GenDetectorLike(nelem, es, 70000 + es);

        std::vector<size_t> bsizes = {0, bshuf_default_block_size(es), 128, 512, 1024, 4096,
                                      8192 / es, 32768 / es, 49152 / es, 65536 / es};
        std::sort(bsizes.begin(), bsizes.end());
        bsizes.erase(std::unique(bsizes.begin(), bsizes.end()), bsizes.end());

        for (size_t bs : bsizes) {
            // bshuf_compress_lz4 emits only the blocks - the 12-byte container header (total size
            // BE64, block size in BYTES BE32) is the caller's, which is why JFJochCompressor writes
            // it by hand and every decoder starts reading at +12. Build it here the same way, or the
            // container is not the one the pipeline reads.
            const size_t block_elems = (bs != 0) ? bs : bshuf_default_block_size(es);
            const size_t block_bytes = block_elems * es;
            std::vector<uint8_t> out(12 + bshuf_compress_lz4_bound(nelem, es, block_elems));
            PutBE64(out.data(), (uint64_t)nelem * es);
            PutBE32(out.data() + 8, (uint32_t)block_bytes);
            const int64_t clen = bshuf_compress_lz4(img.data(), out.data() + 12, nelem, es, block_elems);
            REQUIRE(clen > 0);
            out.resize(12 + (size_t)clen);
            const size_t nblocks = (nelem + block_elems - 1) / block_elems;

            const CompressedImage image(out.data(), out.size(), nelem, 1, ModeOf(es, false), ALG);
            std::vector<uint8_t> cpu;
            image.GetUncompressed(cpu);
            { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); }

            auto stream = std::make_shared<CudaStream>();
            BSLZ4DecoderGPU decoder(nelem * es, stream);
            CudaDevicePtr<uint8_t> gpu_out(nelem * es);
            cudaMemset(gpu_out.get(), 0xA5, nelem * es);

            std::string what;
            try {
                decoder.Decode(image, gpu_out.get());
                const cudaError_t sync = cudaStreamSynchronize(*stream);
                if (sync != cudaSuccess) what = cudaGetErrorString(sync);
            } catch (const std::exception &e) { what = e.what(); }

            if (!what.empty()) {
                printf("[BSLZ4Fuzz] foreign block %6zu bytes, es=%u (%5zu blocks): GPU FAILED - %s\n",
                       block_bytes, es, nblocks, what.c_str());
                n_refused++;
                cudaGetLastError();   // clear, so later cases are not confused by a stale error
                continue;
            }

            std::vector<uint8_t> gpu(nelem * es);
            REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), nelem * es, cudaMemcpyDeviceToHost) == cudaSuccess);
            std::ostringstream l; l << "foreign block_bytes=" << block_bytes << " es=" << es;
            const auto msg = DescribeMismatch(cpu.data(), gpu.data(), nelem * es, block_bytes,
                                              nelem * es, l.str());
            if (!msg.empty()) FAIL(msg);
            printf("[BSLZ4Fuzz] foreign block %6zu bytes, es=%u (%5zu blocks): byte-exact\n",
                   block_bytes, es, nblocks);
            n_ok++;
        }
    }
    printf("[BSLZ4Fuzz] foreign bitshuffle block sizes: %zu byte-exact, %zu refused by the GPU\n",
           n_ok, n_refused);
    SUCCEED();
}

// =============================================================================================
// Targeted probe for a WARP-SYNCHRONY race.
//
// lz4_decode_blocks splits each copy across the 32 lanes with `for (i = lane; i < n; i += 32)`, so
// the trip count is lane-dependent and the lanes DIVERGE at the loop exit. There is no
// __syncwarp() before the next sequence's copy, and on Volta+ (this is sm_120) independent thread
// scheduling gives no implicit warp-lockstep guarantee. A lane that finished early can therefore
// run ahead and read a byte another lane has not written yet.
//
// The construction below makes exactly that overlap: a literal run of length L == 1 (mod 32), so
// the LAST byte is written by lane 0 on its SECOND trip while lanes 1..31 have already left the
// loop, followed immediately by a match at a tiny offset that reads that very byte.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_SyntheticWarpRaceProbe", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");

    const uint32_t es = 4;
    const size_t B = BlockElems(es), block_bytes = B * es;
    auto stream = std::make_shared<CudaStream>();
    BSLZ4DecoderGPU decoder(block_bytes, stream);
    CudaDevicePtr<uint8_t> gpu_out(block_bytes);

    size_t configs = 0, total_decodes = 0, bad_decodes = 0;

    // Literal-run lengths chosen around multiples of the 32-lane width, and matches whose source
    // is the tail the previous copy just produced.
    for (uint32_t L : {1u, 2u, 17u, 31u, 32u, 33u, 34u, 63u, 64u, 65u, 96u, 97u}) {
        for (uint32_t off : {1u, 2u, 3u, 4u, 5u, 8u, 15u, 16u, 17u, 31u, 32u, 33u}) {
            for (uint32_t ml : {4u, 5u, 33u, 64u, 65u}) {
                std::mt19937 rng(60000 + L * 977 + off * 31 + ml);
                Lz4BlockBuilder bb;
                {
                    std::vector<uint8_t> prime(256);
                    for (auto &b : prime) b = (uint8_t)(rng() & 0xFF);
                    bb.AddLiterals(prime);
                    bb.Match(128, 128);
                }
                // As many back-to-back [literals L][match off,ml] pairs as fit.
                while (bb.OutputPos() + L + ml + 32 < block_bytes) {
                    std::vector<uint8_t> lits(L);
                    for (auto &b : lits) b = (uint8_t)(rng() & 0xFF);
                    bb.AddLiterals(lits);
                    bb.Match(off, ml);
                }
                std::vector<uint8_t> tail(block_bytes - bb.OutputPos());
                for (auto &b : tail) b = (uint8_t)(rng() & 0xFF);
                bb.Finish(tail);
                REQUIRE(bb.plain.size() == block_bytes);

                std::vector<uint8_t> chk(block_bytes);
                REQUIRE(LZ4_decompress_safe((const char *)bb.code.data(), (char *)chk.data(),
                                            (int)bb.code.size(), (int)block_bytes) == (int)block_bytes);
                { const bool lz4_ok = (chk == bb.plain); REQUIRE(lz4_ok); }

                const auto img = ImageFromShuffled(bb.plain, B, es, false);
                const auto container = MakeContainer({bb.code}, block_bytes, block_bytes);
                const CompressedImage image(container.data(), container.size(), B, 1,
                                            CompressedImageMode::Uint32, ALG);
                std::vector<uint8_t> cpu;
                image.GetUncompressed(cpu);
                { const bool cpu_ok = (cpu == img); REQUIRE(cpu_ok); }

                // A race is probabilistic, so decode the same block many times.
                std::vector<uint8_t> gpu(block_bytes);
                for (int rep = 0; rep < 25; rep++) {
                    cudaMemset(gpu_out.get(), 0xA5, block_bytes);
                    decoder.Decode(image, gpu_out.get());
                    REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess);
                    REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), block_bytes,
                                       cudaMemcpyDeviceToHost) == cudaSuccess);
                    total_decodes++;
                    if (memcmp(gpu.data(), cpu.data(), block_bytes) != 0) {
                        bad_decodes++;
                        if (bad_decodes == 1) {
                            std::ostringstream l;
                            l << "warp-race probe L=" << L << " offset=" << off << " matchlen=" << ml
                              << " rep=" << rep;
                            FAIL(DescribeMismatch(cpu.data(), gpu.data(), block_bytes,
                                                  block_bytes, block_bytes, l.str()));
                        }
                    }
                }
                configs++;
            }
        }
    }
    printf("[BSLZ4Fuzz] warp-synchrony probe: %zu configurations x 25 repeats = %zu decodes, "
           "%zu corrupted\n", configs, total_decodes, bad_decodes);
    SUCCEED();
}

// =============================================================================================
// CORRUPT payloads. An HDF5 chunk carries no checksum unless Fletcher32 is on, so a bit flip on the
// wire or on disk reaches the decoder looking like a plausible container. The host decoder catches
// it - LZ4_decompress_safe insists the block decode to exactly its declared length - and the device
// decoder MUST do the same, because the buffers it decodes into are reused frame to frame: a block
// that stops early would otherwise leave the PREVIOUS image's bytes in place, and those bytes are
// the most significant byte-plane of the block, so they do not look like a missing corner. They
// look like real pixels several powers of two too bright.
//
// The invariant asserted here is the one that matters: whenever the CPU rejects a chunk, the GPU
// must reject it too. (The converse is not required - a flip can leave a stream that still decodes
// to the right length, and no LZ4 decoder can detect that.)
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_CorruptPayload", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");

    const size_t nelem = 200000, es = 4;
    const auto img = GenDetectorLike(nelem, es, 4242);
    JFJochBitShuffleCompressor compressor(ALG);
    const auto clean = compressor.Compress(img.data(), nelem, es);

    auto stream = std::make_shared<CudaStream>();
    BSLZ4DecoderGPU decoder(nelem * es, stream);
    CudaDevicePtr<uint8_t> gpu_out(nelem * es);

    std::mt19937 rng(777);
    size_t both_threw = 0, both_ok = 0, gpu_missed = 0, cpu_only_ok = 0;

    for (int trial = 0; trial < 400; trial++) {
        std::vector<uint8_t> bad = clean;
        // Corrupt inside the LZ4 payload area, never the 12-byte container header - a header flip is
        // a different (and already covered) failure mode.
        const size_t pos = 12 + rng() % (bad.size() - 12);
        bad[pos] ^= (uint8_t)(1u << (rng() % 8));

        const CompressedImage image(bad.data(), bad.size(), nelem, 1, CompressedImageMode::Uint32, ALG);

        bool cpu_threw = false;
        std::vector<uint8_t> cpu;
        try { image.GetUncompressed(cpu); } catch (const std::exception &) { cpu_threw = true; }

        // Prime the device buffers with a recognisable pattern, so a block that quietly stopped
        // early would show up as leftover poison rather than as plausible data.
        cudaMemset(gpu_out.get(), 0x5A, nelem * es);
        bool gpu_threw = false;
        try {
            decoder.Decode(image, gpu_out.get());
            REQUIRE(cudaStreamSynchronize(*stream) == cudaSuccess);
            decoder.ThrowIfDecodeFailed();
        } catch (const std::exception &) { gpu_threw = true; }

        // Whatever happened, it must not have been a device fault: a corrupt length must be clamped,
        // never walked off the end of a buffer.
        REQUIRE(cudaGetLastError() == cudaSuccess);

        if (cpu_threw && gpu_threw) both_threw++;
        else if (!cpu_threw && !gpu_threw) {
            both_ok++;
            std::vector<uint8_t> gpu(nelem * es);
            REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), nelem * es, cudaMemcpyDeviceToHost) == cudaSuccess);
            const auto msg = DescribeMismatch(cpu.data(), gpu.data(), nelem * es,
                                              DefaultBlockBytes(), nelem * es,
                                              "corrupt payload both-accepted");
            if (!msg.empty()) FAIL(msg);   // if both accept it, they must agree on the result
        }
        else if (cpu_threw && !gpu_threw) {
            gpu_missed++;
            // Whatever the GPU accepted, it must at least be a COMPLETE decode. The poison check is
            // the important one: a block that stopped early would leave 0x5A here, and in production
            // those bytes are the previous frame rather than poison.
            std::vector<uint8_t> gpu(nelem * es);
            REQUIRE(cudaMemcpy(gpu.data(), gpu_out.get(), nelem * es, cudaMemcpyDeviceToHost) == cudaSuccess);
            size_t poison = 0;
            for (size_t i = 0; i < nelem * es; i++) if (gpu[i] == 0x5A) poison++;
            INFO("corrupt flip at byte " << pos << " left " << poison << " undecoded bytes");
            REQUIRE(poison == 0);
        }
        else cpu_only_ok++;
    }

    printf("[BSLZ4Fuzz] corrupt payload, 400 single-bit flips: %zu both rejected, %zu both accepted "
           "(and agreed), %zu GPU MISSED what the CPU caught, %zu GPU stricter than CPU\n",
           both_threw, both_ok, gpu_missed, cpu_only_ok);

    // The device decoder reproduces the reference's length, bounds and parsing-restriction checks,
    // which accounts for the large majority of what LZ4_decompress_safe rejects. It does NOT
    // reproduce the reference's fast-loop/safe-loop selection exactly, so a small tail of corrupt
    // streams still decodes here that the host would refuse. Those are complete decodes - the block
    // is fully written, asserted above - differing from the true image in a handful of bytes, not
    // partial decodes leaking a previous frame, which is the failure mode that actually matters.
    // The bound is a regression guard: it must not silently get worse.
    CHECK(gpu_missed * 100 <= 400 * 2);   // <= 2% of single-bit flips
}

// Malformed CONTAINERS - the parts the host can reject before any work is queued. Each of these is
// rejected by the CPU decoder, so the device decoder must reject them too rather than decode
// something plausible-looking out of uninitialised memory.
TEST_CASE("BSLZ4Fuzz_MalformedContainer", "[BSLZ4Fuzz]") {
    if (get_gpu_count() == 0) SKIP("No CUDA GPU present");

    const size_t nelem = 100000, es = 4;
    const auto img = GenDetectorLike(nelem, es, 99);
    JFJochBitShuffleCompressor compressor(ALG);
    const auto clean = compressor.Compress(img.data(), nelem, es);

    auto stream = std::make_shared<CudaStream>();
    BSLZ4DecoderGPU decoder(nelem * es, stream);
    CudaDevicePtr<uint8_t> gpu_out(nelem * es);

    auto expect_throw = [&](const std::vector<uint8_t> &chunk, const char *what) {
        const CompressedImage image(chunk.data(), chunk.size(), nelem, 1, CompressedImageMode::Uint32, ALG);
        bool threw = false;
        try {
            decoder.Decode(image, gpu_out.get());
            cudaStreamSynchronize(*stream);
            decoder.ThrowIfDecodeFailed();
        } catch (const std::exception &) { threw = true; }
        cudaGetLastError();
        INFO(what);
        CHECK(threw);
    };

    { // Block size that is not a multiple of 8 elements: bitshuffle transposes 8 at a time, so the
      // un-transpose would leave part of every block unwritten.
        auto bad = clean; PutBE32(bad.data() + 8, 4 * es);
        expect_throw(bad, "block size not a multiple of 8 elements");
    }
    { // A one-element block size implies millions of blocks the chunk cannot possibly hold. Must be
      // rejected on the arithmetic, BEFORE it becomes a several-hundred-MB pinned allocation.
        auto bad = clean; PutBE32(bad.data() + 8, es);
        expect_throw(bad, "block size implies more blocks than the chunk holds");
    }
    { auto bad = clean; PutBE64(bad.data(), (uint64_t)nelem * es + 8);
      expect_throw(bad, "declared total size disagrees with the image"); }
    { auto bad = clean; PutBE32(bad.data() + 8, 0);
      expect_throw(bad, "zero block size"); }
    { auto bad = clean; bad.resize(11);
      expect_throw(bad, "shorter than the container header"); }
    { auto bad = clean; bad.resize(clean.size() / 2);
      expect_throw(bad, "truncated mid-stream"); }
    { auto bad = clean; bad.push_back(0); bad.push_back(0);
      expect_throw(bad, "trailing bytes after the last block"); }
    { // A block length field larger than what is left in the chunk.
        auto bad = clean; PutBE32(bad.data() + 12, (uint32_t)clean.size());
        expect_throw(bad, "block length runs past the chunk");
    }
    { auto bad = clean; PutBE32(bad.data() + 12, 0);
      expect_throw(bad, "zero-length block"); }
}

// =============================================================================================
// Coverage report - what the fuzz cases above actually made LZ4 emit.
// =============================================================================================
TEST_CASE("BSLZ4Fuzz_ZZ_CoverageReport", "[BSLZ4Fuzz]") {
    std::lock_guard<std::mutex> lock(g_stats_mtx);
    const auto &s = g_stats;
    printf("\n[BSLZ4Fuzz] ===== LZ4 CODE-PATH COVERAGE (production-compressor cases) =====\n");
    printf("  images parsed .................. %zu\n", s.images);
    printf("  bitshuffle blocks .............. %zu\n", s.blocks);
    printf("  LZ4 match sequences ............ %zu\n", s.sequences);
    printf("  litlen field == 15 (extended) .. %zu\n", s.lit_ge15);
    printf("  litlen 255-extension CHAIN ..... %zu\n", s.lit_ext255);
    printf("  max literal run ................ %zu bytes\n", s.max_litlen);
    printf("  matchlen field == 15 (extended)  %zu\n", s.ml_ge19);
    printf("  matchlen 255-extension CHAIN ... %zu\n", s.ml_ext255);
    printf("  max match length ............... %zu bytes\n", s.max_matchlen);
    printf("  offset >= matchlen (non-overlap) %zu\n", s.nonoverlap);
    printf("  offset <  matchlen (OVERLAP)  .. %zu\n", s.overlap);
    printf("  offset == matchlen ............. %zu\n", s.off_eq_ml);
    printf("  offset == matchlen + 1 ......... %zu\n", s.off_ml_plus1);
    printf("  offset == matchlen - 1 ......... %zu\n", s.off_ml_minus1);
    printf("  smallest offset seen ........... %u\n", s.min_offset);
    printf("  offsets <= 64 seen: ");
    for (auto &kv : s.offsets_le64) printf("%u(x%zu) ", kv.first, kv.second);
    printf("\n[BSLZ4Fuzz] =============================================================\n\n");
    SUCCEED();
}

#endif