diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml new file mode 100644 index 0000000..70bcb30 --- /dev/null +++ b/.github/workflows/format.yml @@ -0,0 +1,21 @@ +name: C/C++ Formatting + +on: [push, pull_request] + +permissions: + contents: read + +jobs: + clang-format: + runs-on: ubuntu-24.04 + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: '3.12' + - name: Install pinned formatter + run: python -m pip install clang-format==18.1.3 + - name: Check tracked owned C/C++ formatting + env: + CLANG_FORMAT: clang-format + run: bash scripts/format.sh --check diff --git a/CMakeLists.txt b/CMakeLists.txt index e385c25..ab5c53f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,6 +22,7 @@ endif() option(GF256_BUILD_TESTS "Build unit tests" ON) option(GF256_BUILD_TOOLS "Build CLI tools" ON) +option(GF256_BUILD_MONTE_CARLO "Build the product Monte Carlo CLI (also in test profiles)" ON) option(GF256_BUILD_BENCHMARKS "Build benchmarks" ON) option(GF256_BUILD_REFERENCE_BENCHMARKS "Build external RS baselines" OFF) option(GF256_ENABLE_NATIVE_ISA "Compile owned sources for the build host" OFF) @@ -62,8 +63,10 @@ add_library(gf256_core STATIC src/reed_solomon/decoder/high_rate.cc src/reed_solomon/lch_encoder.cc src/reed_solomon/lch_decoder.cc + src/reed_solomon/strong_weak_rs_product_code.cc ) target_compile_features(gf256_core PUBLIC cxx_std_20) +set_target_properties(gf256_core PROPERTIES POSITION_INDEPENDENT_CODE ON) target_include_directories(gf256_core PUBLIC ${PROJECT_SOURCE_DIR}/include PRIVATE ${PROJECT_SOURCE_DIR}/src) @@ -85,6 +88,40 @@ if(GF256_ENABLE_CODEWORD_CANTOR_AFFINE_EXPERIMENT) GF256_ENABLE_CODEWORD_CANTOR_AFFINE_EXPERIMENT=1) endif() +if(GF256_BUILD_MONTE_CARLO OR GF256_BUILD_TOOLS OR GF256_BUILD_TESTS) + FetchContent_Declare(nlohmann_json + GIT_REPOSITORY https://github.com/nlohmann/json.git + GIT_TAG v3.12.0 + GIT_SHALLOW TRUE) + set(JSON_BuildTests OFF CACHE BOOL "" FORCE) + set(JSON_Install OFF CACHE BOOL "" FORCE) + FetchContent_MakeAvailable(nlohmann_json) +endif() + +if(GF256_BUILD_MONTE_CARLO OR GF256_BUILD_TOOLS OR GF256_BUILD_TESTS) + add_executable(rs-product-test tools/rs_product_test.cc) + target_link_libraries(rs-product-test PRIVATE gf256_core nlohmann_json::nlohmann_json) + include(GNUInstallDirs) + install(TARGETS rs-product-test + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT product-test) +endif() + +if(GF256_BUILD_MONTE_CARLO) + find_package(Threads REQUIRED) + add_library(product_monte_carlo_trials OBJECT tools/product_monte_carlo_native.cc) + target_link_libraries(product_monte_carlo_trials PRIVATE gf256_core) + set_target_properties(product_monte_carlo_trials PROPERTIES POSITION_INDEPENDENT_CODE ON) + if(GF256_ENABLE_NATIVE_ISA AND CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang") + target_compile_options(product_monte_carlo_trials PRIVATE -march=native) + endif() + add_executable(rs-product-monte-carlo tools/product_monte_carlo.cc) + target_link_libraries(rs-product-monte-carlo PRIVATE product_monte_carlo_trials + gf256_core Threads::Threads nlohmann_json::nlohmann_json) + include(GNUInstallDirs) + install(TARGETS rs-product-monte-carlo + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT product-monte-carlo) +endif() + if(GF256_BUILD_TOOLS) set(XXHASH_BUILD_XXHSUM OFF CACHE BOOL "" FORCE) set(INDICATORS_BUILD_TESTS OFF CACHE BOOL "" FORCE) @@ -113,7 +150,7 @@ if(GF256_BUILD_TOOLS) target_link_libraries(lch-rs PRIVATE gf256_core xxhash indicators::indicators Threads::Threads) include(GNUInstallDirs) - install(TARGETS lch-rs RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) + install(TARGETS lch-rs RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT lch-rs) endif() if(GF256_ENABLE_GFNI512_RADIX8_EXPERIMENT) @@ -137,12 +174,18 @@ endif() if(GF256_BUILD_TESTS) FetchContent_MakeAvailable(googletest) enable_testing() + find_package(Python3 3.9 REQUIRED COMPONENTS Interpreter) + add_test(NAME RSProductTestCli + COMMAND ${Python3_EXECUTABLE} ${PROJECT_SOURCE_DIR}/tests/rs_product_test.py + $) + set_tests_properties(RSProductTestCli PROPERTIES TIMEOUT 60) add_executable(gf_unittests tests/field_tests.cc) target_link_libraries(gf_unittests PRIVATE gf256_core gtest gtest_main) add_test(NAME Unittests COMMAND gf_unittests) - add_executable(lch_rs_unittests tests/lch_tests.cc tests/rs_tests.cc) + add_executable(lch_rs_unittests tests/lch_tests.cc tests/rs_tests.cc + tests/product_code_tests.cc) target_include_directories(lch_rs_unittests PRIVATE ${PROJECT_SOURCE_DIR}/src) target_link_libraries(lch_rs_unittests PRIVATE gf256_core gtest gtest_main) if(GF256_ENABLE_NATIVE_ISA AND CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang") @@ -159,6 +202,41 @@ if(GF256_BUILD_TESTS) GF256_ENABLE_CODEWORD_CANTOR_AFFINE_EXPERIMENT=1) endif() add_test(NAME LCHRSUnittests COMMAND lch_rs_unittests) + if(GF256_BUILD_MONTE_CARLO) + find_package(Python3 3.9 REQUIRED COMPONENTS Interpreter) + add_library(product_monte_carlo_native SHARED $) + set_target_properties(product_monte_carlo_native PROPERTIES PREFIX "") + target_link_libraries(product_monte_carlo_native PRIVATE gf256_core Threads::Threads) + configure_file(tests/reference/product_monte_carlo.py product_monte_carlo_reference.py COPYONLY) + add_executable(product_monte_carlo_data_tests tests/product_monte_carlo_data_tests.cc) + target_include_directories(product_monte_carlo_data_tests PRIVATE tools) + target_link_libraries(product_monte_carlo_data_tests PRIVATE nlohmann_json::nlohmann_json gtest_main) + target_compile_features(product_monte_carlo_data_tests PRIVATE cxx_std_20) + add_test(NAME ProductMonteCarloData COMMAND product_monte_carlo_data_tests) + add_executable(product_monte_carlo_fault_cli tools/product_monte_carlo.cc) + target_compile_definitions(product_monte_carlo_fault_cli PRIVATE GF256_MC_TEST_HOOKS=1) + target_link_libraries(product_monte_carlo_fault_cli PRIVATE product_monte_carlo_trials + gf256_core Threads::Threads nlohmann_json::nlohmann_json) + add_test(NAME ProductMonteCarloCli + COMMAND ${Python3_EXECUTABLE} ${PROJECT_SOURCE_DIR}/tests/product_monte_carlo_test.py + ${PROJECT_BINARY_DIR}/rs-product-monte-carlo) + set_tests_properties(ProductMonteCarloCli PROPERTIES TIMEOUT 180) + add_test(NAME ProductMonteCarloLegacy + COMMAND ${Python3_EXECUTABLE} ${PROJECT_SOURCE_DIR}/tests/product_monte_carlo_legacy_test.py + ${PROJECT_BINARY_DIR}/product_monte_carlo_reference.py) + set_tests_properties(ProductMonteCarloLegacy PROPERTIES TIMEOUT 180) + add_test(NAME ProductMonteCarloPlot + COMMAND ${Python3_EXECUTABLE} ${PROJECT_SOURCE_DIR}/tests/plot_product_monte_carlo_test.py) + if(CMAKE_CXX_FLAGS MATCHES "fsanitize=address") + # ctypes loads an instrumented DSO into an otherwise uninstrumented Python. + execute_process(COMMAND ${CMAKE_CXX_COMPILER} -print-file-name=libasan.so + OUTPUT_VARIABLE GF256_ASAN_RUNTIME OUTPUT_STRIP_TRAILING_WHITESPACE) + set_tests_properties(ProductMonteCarloLegacy PROPERTIES + ENVIRONMENT "LD_PRELOAD=${GF256_ASAN_RUNTIME}") + set_tests_properties(ProductMonteCarloCli PROPERTIES + ENVIRONMENT "MC_REFERENCE_LD_PRELOAD=${GF256_ASAN_RUNTIME}") + endif() + endif() if(GF256_BUILD_TOOLS) add_test(NAME LCHRSCli @@ -170,8 +248,10 @@ endif() if(GF256_BUILD_BENCHMARKS) FetchContent_MakeAvailable(googlebenchmark) - add_executable(benchmarks benchmarks/matrix_multiplication.cc) + add_executable(benchmarks benchmarks/matrix_multiplication.cc + benchmarks/strong_weak_rs_product_code_benchmarks.cc) target_link_libraries(benchmarks PRIVATE gf256_core benchmark benchmark_main) + target_include_directories(benchmarks PRIVATE ${PROJECT_SOURCE_DIR}/src) add_executable(rs_verbose_benchmarks benchmarks/lch_rs_benchmarks.cc diff --git a/CMakePresets.json b/CMakePresets.json index 485e146..5ff8a78 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -6,6 +6,22 @@ "patch": 0 }, "configurePresets": [ + { + "name": "experimental", + "displayName": "Experimental Tools", + "description": "Native Release build of the RS product Monte Carlo and fixture tools", + "binaryDir": "${sourceDir}/build/experimental-preset", + "cacheVariables": { + "CMAKE_BUILD_TYPE": "Release", + "CMAKE_EXPORT_COMPILE_COMMANDS": "ON", + "GF256_BUILD_MONTE_CARLO": "ON", + "GF256_BUILD_TOOLS": "OFF", + "GF256_BUILD_TESTS": "OFF", + "GF256_BUILD_BENCHMARKS": "OFF", + "GF256_BUILD_REFERENCE_BENCHMARKS": "OFF", + "GF256_ENABLE_NATIVE_ISA": "ON" + } + }, { "name": "release", "displayName": "Release", @@ -61,6 +77,15 @@ } ], "buildPresets": [ + { + "name": "experimental", + "displayName": "Build Experimental Tools", + "configurePreset": "experimental", + "targets": [ + "rs-product-monte-carlo", + "rs-product-test" + ] + }, { "name": "release", "displayName": "Build Release", diff --git a/benchmarks/reference/xdrs_benchmarks.cc b/benchmarks/reference/xdrs_benchmarks.cc index ceff514..8f5c05e 100644 --- a/benchmarks/reference/xdrs_benchmarks.cc +++ b/benchmarks/reference/xdrs_benchmarks.cc @@ -71,10 +71,11 @@ void InitializeXDRS(XDRSTables& tables, unsigned bytes, unsigned k) { ::function::init_dec(); } -bool VerifyXDRSRecovery( - unsigned k, unsigned bytes, bool low_rate, - const std::vector>& data, - const std::vector>& recovery) { +bool VerifyXDRSRecovery(unsigned k, + unsigned bytes, + bool low_rate, + const std::vector>& data, + const std::vector>& recovery) { const unsigned recovery_count = Size - k; auto data_pointers = ConstPointers(data); auto recovery_pointers = ConstPointers(recovery); diff --git a/benchmarks/strong_weak_rs_product_code_benchmarks.cc b/benchmarks/strong_weak_rs_product_code_benchmarks.cc new file mode 100644 index 0000000..9d66e9a --- /dev/null +++ b/benchmarks/strong_weak_rs_product_code_benchmarks.cc @@ -0,0 +1,303 @@ +#include +#include +#include +#include +#include +#include + +#include "benchmark/benchmark.h" +#include "reed_solomon/product_code_internal.h" +#include "reed_solomon/strong_weak_rs_product_code.h" + +namespace { + +void BenchmarkProductCorrectionBSC(benchmark::State& state, + int batch_passes = -1, + size_t weak_n = 256, + int optimizations = -1, + size_t weak_r = 2) { + using gf2p8::Element; + using gf2p8::rs::ProductCorrectionResult; + using gf2p8::rs::ProductTermination; + constexpr size_t kN = 256; + constexpr size_t kStrongK = 224; + const size_t kWeakK = weak_n - weak_r; + const size_t kInformationBytes = kStrongK * kWeakK; + const size_t kBlockBytes = kN * weak_n; + constexpr size_t kCorpusCount = 64; + constexpr size_t kPassLimit = 16; + constexpr uint32_t kSeed = 0x5b5c0224; + gf2p8::rs::StrongWeakRSProductCode code(kN, kStrongK, weak_n, kWeakK); + const auto correct = [&](std::vector& block) { + if (optimizations >= 0) { + return gf2p8::rs::detail::ProductCorrectionAccess::Experiment( + code, block, gf2p8::rs::ProductDecodeOptions{kPassLimit}, + optimizations); + } + return batch_passes < 0 + ? code.Correct(block, kPassLimit) + : gf2p8::rs::detail::ProductCorrectionAccess::Correct( + code, block, kPassLimit, batch_passes); + }; + if (!code.Valid() || code.BlockSize() != kBlockBytes) { + state.SkipWithError("invalid product-code dimensions"); + return; + } + + std::mt19937 messages(kSeed); + std::mt19937 channel(kSeed ^ 0x9e3779b9U); + std::vector> original(kCorpusCount, + std::vector(kBlockBytes)); + auto corrupted = original; + auto work = original; + std::vector results(kCorpusCount); + uint64_t channel_bits = 0, channel_symbols = 0; + for (size_t sample = 0; sample < kCorpusCount; ++sample) { + auto& block = original[sample]; + for (size_t row = 0; row < kStrongK; ++row) { + for (size_t col = 0; col < kWeakK; ++col) { + block[row * weak_n + col] = static_cast(messages()); + } + } + if (code.Encode(block) != gf2p8::lch::Status::ok) { + state.SkipWithError("product corpus encoding failed"); + return; + } + work[sample] = block; + const auto clean = code.Correct(work[sample], kPassLimit); + if (!clean.all_zero_syndromes || clean.changed_symbols != 0 || + work[sample] != block) { + state.SkipWithError("encoded product corpus is not valid"); + return; + } + corrupted[sample] = block; + for (auto& value : corrupted[sample]) { + const auto before = value; + for (unsigned bit = 0; bit < 8; ++bit) { + // Reject the incomplete residue range: exact Bernoulli(1/200), + // reproducible across standard libraries, with no fixed error count. + uint32_t draw; + do { + draw = static_cast(channel()); + } while (draw >= 4294967200U); + if (draw % 200 == 0) { + value ^= static_cast(1U << bit); + ++channel_bits; + } + } + channel_symbols += value != before; + } + } + + uint64_t data_residual_bits = 0, codeword_residual_bits = 0; + uint64_t message_failures = 0, block_failures = 0, validity_failures = 0; + uint64_t valid_wrong_blocks = 0, pass_limits = 0, passes = 0; + uint64_t strong_lines = 0, weak_lines = 0; + uint64_t accepted_bits = 0, accepted_symbols = 0; + uint64_t strong_bits = 0, strong_symbols = 0, weak_bits = 0, weak_symbols = 0; + size_t max_passes = 0; + // Validate this exact channel corpus against the retained single path outside + // timing, including scheduling outcomes, not just final message recovery. + for (size_t sample = 0; sample < kCorpusCount; ++sample) { + auto reference = corrupted[sample]; + const auto expected = gf2p8::rs::detail::ProductCorrectionAccess::Correct( + code, reference, kPassLimit, 0, false); + work[sample] = corrupted[sample]; + const auto actual = correct(work[sample]); + if (work[sample] != reference || + actual.termination != expected.termination || + actual.all_zero_syndromes != expected.all_zero_syndromes || + actual.directional_passes != expected.directional_passes || + actual.strong_lines_visited != expected.strong_lines_visited || + actual.weak_lines_visited != expected.weak_lines_visited || + actual.changed_symbols != expected.changed_symbols || + actual.changed_bits != expected.changed_bits || + actual.strong_changed_bits != expected.strong_changed_bits || + actual.weak_changed_bits != expected.weak_changed_bits || + actual.strong_changed_symbols != expected.strong_changed_symbols || + actual.weak_changed_symbols != expected.weak_changed_symbols) { + state.SkipWithError("batch/single corpus differential mismatch"); + return; + } + } + for (auto _ : state) { + state.PauseTiming(); + for (size_t sample = 0; sample < kCorpusCount; ++sample) { + std::copy(corrupted[sample].begin(), corrupted[sample].end(), + work[sample].begin()); + } + state.ResumeTiming(); + for (size_t sample = 0; sample < kCorpusCount; ++sample) { + results[sample] = correct(work[sample]); + benchmark::DoNotOptimize(results[sample]); + benchmark::ClobberMemory(); + } + state.PauseTiming(); + for (size_t sample = 0; sample < kCorpusCount; ++sample) { + const auto& result = results[sample]; + if (result.termination == ProductTermination::invalid_argument) { + state.SkipWithError("product correction rejected corpus input"); + break; + } + uint64_t data_bits = 0, block_bits = 0; + for (size_t pos = 0; pos < kBlockBytes; ++pos) { + const auto bits = std::popcount( + static_cast(work[sample][pos] ^ original[sample][pos])); + block_bits += bits; + if (pos / weak_n < kStrongK && pos % weak_n < kWeakK) { + data_bits += bits; + } + } + data_residual_bits += data_bits; + codeword_residual_bits += block_bits; + message_failures += data_bits != 0; + block_failures += block_bits != 0; + validity_failures += !result.all_zero_syndromes; + valid_wrong_blocks += result.all_zero_syndromes && block_bits != 0; + pass_limits += result.termination == ProductTermination::pass_limit; + passes += result.directional_passes; + max_passes = std::max(max_passes, result.directional_passes); + strong_lines += result.strong_lines_visited; + weak_lines += result.weak_lines_visited; + accepted_bits += result.changed_bits; + accepted_symbols += result.changed_symbols; + strong_bits += result.strong_changed_bits; + strong_symbols += result.strong_changed_symbols; + weak_bits += result.weak_changed_bits; + weak_symbols += result.weak_changed_symbols; + } + state.ResumeTiming(); + if (state.skipped()) { + break; + } + } + + const double sweeps = static_cast(state.iterations()); + if (sweeps == 0 || state.skipped()) { + return; + } + const double blocks = sweeps * kCorpusCount; + state.counters["corpus_blocks"] = kCorpusCount; + state.counters["timed_blocks"] = blocks; + state.counters["seed"] = kSeed; + state.counters["information_bytes_per_block"] = kInformationBytes; + state.counters["target_bit_probability"] = 0.005; + state.counters["channel_codeword_BER"] = + static_cast(channel_bits) / (kCorpusCount * kBlockBytes * 8); + state.counters["channel_flipped_bits"] = static_cast(channel_bits); + state.counters["channel_corrupted_bytes"] = + static_cast(channel_symbols); + state.counters["corpus_data_residual_bits"] = data_residual_bits / sweeps; + state.counters["corpus_codeword_residual_bits"] = + codeword_residual_bits / sweeps; + state.counters["data_residual_BER"] = + data_residual_bits / (blocks * kInformationBytes * 8); + state.counters["codeword_residual_BER"] = + codeword_residual_bits / (blocks * kBlockBytes * 8); + // Counts per unique corpus, not inflated by timed replays of the same noise. + state.counters["corpus_message_failures"] = message_failures / sweeps; + state.counters["corpus_block_failures"] = block_failures / sweeps; + state.counters["corpus_validity_failures"] = validity_failures / sweeps; + state.counters["corpus_valid_wrong_blocks"] = valid_wrong_blocks / sweeps; + state.counters["message_recovery_fraction"] = 1 - message_failures / blocks; + state.counters["block_recovery_fraction"] = 1 - block_failures / blocks; + state.counters["corpus_pass_limits"] = pass_limits / sweeps; + state.counters["mean_directional_passes"] = passes / blocks; + state.counters["max_directional_passes"] = static_cast(max_passes); + state.counters["mean_strong_lines"] = strong_lines / blocks; + state.counters["mean_weak_lines"] = weak_lines / blocks; + state.counters["mean_accepted_bit_changes"] = accepted_bits / blocks; + state.counters["mean_accepted_byte_changes"] = accepted_symbols / blocks; + state.counters["mean_strong_accepted_bit_changes"] = strong_bits / blocks; + state.counters["mean_strong_accepted_byte_changes"] = strong_symbols / blocks; + state.counters["mean_weak_accepted_bit_changes"] = weak_bits / blocks; + state.counters["mean_weak_accepted_byte_changes"] = weak_symbols / blocks; + state.counters["pass_limit"] = kPassLimit; + state.SetItemsProcessed(state.iterations() * kCorpusCount); + state.SetBytesProcessed(state.iterations() * kCorpusCount * + kInformationBytes); + state.SetLabel( + "64 blocks/iteration; information bytes=Kstrong*Kweak; " + "Correct includes exact final validity; tuned backend"); +} + +const auto* kProductCorrectionBSC = benchmark::RegisterBenchmark( + "LCH/Owned/StrongWeakRSProductCode/Correct/BSC005/" + "Nstrong:256/Kstrong:224/Nweak:256/Kweak:254", + [](benchmark::State& state) { BenchmarkProductCorrectionBSC(state); }); + +const auto kProductExperiments = [] { + for (int variant : {0, 1, 2, 3, 4, 6, 7}) { + for (size_t n : {256u, 175u}) { + const auto name = "LCH/Owned/StrongWeakRSProductCode/Experiment/" + + std::to_string(variant) + "/" + std::to_string(n); + benchmark::RegisterBenchmark(name.c_str(), [=](benchmark::State& state) { + BenchmarkProductCorrectionBSC(state, -1, n, variant); + }); + } + } + return true; +}(); + +const auto* kProductCorrectionSingle = benchmark::RegisterBenchmark( + "LCH/Owned/StrongWeakRSProductCode/Correct/BSC005/Single", + [](benchmark::State& state) { BenchmarkProductCorrectionBSC(state, 0); }); +const auto* kProductCorrectionR4 = benchmark::RegisterBenchmark( + "LCH/Owned/StrongWeakRSProductCode/R4/Direct256252", + [](benchmark::State& state) { + BenchmarkProductCorrectionBSC(state, -1, 256, -1, 4); + }); +const auto* kProductCorrectionGenericR4 = benchmark::RegisterBenchmark( + "LCH/Owned/StrongWeakRSProductCode/R4/Generic256252", + [](benchmark::State& state) { + BenchmarkProductCorrectionBSC(state, 2, 256, -1, 4); + }); +const auto* kProductCorrectionShortened = benchmark::RegisterBenchmark( + "LCH/Owned/StrongWeakRSProductCode/Correct/BSC005/" + "Nstrong:256/Kstrong:224/Nweak:175/Kweak:173", + [](benchmark::State& state) { + BenchmarkProductCorrectionBSC(state, -1, 175); + }); +const auto* kProductCorrectionStrongBatch = benchmark::RegisterBenchmark( + "LCH/Owned/StrongWeakRSProductCode/Correct/BSC005/StrongBatch", + [](benchmark::State& state) { BenchmarkProductCorrectionBSC(state, 1); }); +const auto* kProductCorrectionBothBatch = benchmark::RegisterBenchmark( + "LCH/Owned/StrongWeakRSProductCode/Correct/BSC005/BothBatch", + [](benchmark::State& state) { BenchmarkProductCorrectionBSC(state, 2); }); + +const auto* kProductCorrectionGenericShortened = benchmark::RegisterBenchmark( + "LCH/Owned/StrongWeakRSProductCode/Correct/BSC005/Generic175", + [](benchmark::State& state) { + BenchmarkProductCorrectionBSC(state, 2, 175); + }); + +void BenchmarkProductEncode(benchmark::State& state) { + const size_t n = state.range(0); + gf2p8::rs::StrongWeakRSProductCode code(256, 224, n, n - 2); + std::vector block(code.BlockSize()); + std::mt19937 random(0x5b5c0224); + for (auto& value : block) { + value = static_cast(random()); + } + if (code.Encode(block) != gf2p8::lch::Status::ok) { + state.SkipWithError("encoding failed"); + return; + } + const auto original = block; + for (auto _ : state) { + benchmark::DoNotOptimize(code.Encode(block)); + benchmark::ClobberMemory(); + } + if (block != original || !code.Correct(block).all_zero_syndromes) { + state.SkipWithError("encoding mismatch"); + } + state.SetBytesProcessed(state.iterations() * 224 * (n - 2)); +} + +const auto* kProductEncode = + benchmark::RegisterBenchmark("LCH/Owned/StrongWeakRSProductCode/Encode", + BenchmarkProductEncode) + ->Arg(256) + ->Arg(175); + +} // namespace diff --git a/docs/strong_weak_rs_acceleration_report.md b/docs/strong_weak_rs_acceleration_report.md new file mode 100644 index 0000000..afe4149 --- /dev/null +++ b/docs/strong_weak_rs_acceleration_report.md @@ -0,0 +1,956 @@ +# Strong-Weak Reed-Solomon Product Code Acceleration Report + +## Applied R=4 Weak Code (2026-09-10) + +**Implemented and retained:** `StrongWeakRSProductCode(256,224,256,252)` now +encodes and corrects with a direct two-error weak decoder. This section supersedes +the earlier statements that R=4 was deferred. Strong kernels and behavior are +unchanged. Default construction remains 256/224 x 256/254, and all existing R=2 +shortenings, including 175/173, retain their coordinate mapping and behavior. +R=4 support is deliberately restricted to full RS(256,252); smaller or shortened +R=4 codes, including the low-rate K=R boundary, are rejected rather than assumed +equivalent. Other already-supported strong dimensions can use this full weak code. + +### Exact Native Parity Checks + +All arithmetic below is in native Cantor GF(256); addition is XOR. Public weak +rows are `[252 data][4 parity]`. Native evaluation points for public data index +`j` are `j+4`, and parity index `252+i` has point `i`. Thus public index `p` +maps to native byte `(p+4) mod 256`. This is an index permutation, not field +addition by integer 4. + +The full native LCH code is evaluation of polynomials of degree at most 251. +Changing from its monic novel basis to ordinary monomials preserves that +polynomial space. For the full field, the vanishing polynomial is +\(P(X)=X^{256}+X\), with \(P'(X)=1\). Consequently all dual evaluation weights +are one. Equivalently, \(\sum_{x\in GF(256)}x^m=0\) for \(0\leq m\leq254\), +including \(m=0\), since 256 is zero in characteristic two. Therefore these +four ordinary power moments vanish on every encoded row: + +\[ + S_j=\sum_x c(x)x^j,\qquad j=0,1,2,3. +\] + +The four check rows have rank four: a nonzero polynomial of degree at most +three cannot vanish at all 256 distinct points. Their kernel thus has dimension +252 and is **exactly** the native code, not merely a necessary validity test. +These moments are not asserted to be the individual novel-basis syndrome +entries used internally by generic FDMA. They are an independent complete +parity-check system for the same code. Scalar LCH encoding and the untouched +generic `CorrectCodeword(LCHDecoder(252,4), ...)` supply independent checks. + +### Locator, Magnitudes, and Degeneracies + +For errors at distinct points \(x_1,x_2\) with nonzero magnitudes \(e_1,e_2\), +\(S_j=e_1x_1^j+e_2x_2^j\). Write the locator as +\(L(X)=X^2+aX+b\), where \(a=x_1+x_2\) and \(b=x_1x_2\). Its recurrence gives + +\[ + \begin{pmatrix}S_1&S_0\\S_2&S_1\end{pmatrix} + \begin{pmatrix}a\\b\end{pmatrix} + =\begin{pmatrix}S_2\\S_3\end{pmatrix},\qquad + D=S_1^2+S_0S_2=e_1e_2(x_1+x_2)^2. +\] + +Thus a genuine two-error pattern always has nonzero determinant, even when +equal magnitudes make \(S_0=0\). No division by \(S_0\) occurs on this branch: + +\[ + a=(S_1S_2+S_0S_3)/D,\qquad + b=(S_1S_3+S_2^2)/D. +\] + +For \(a\ne0\), substitute \(X=ay\) and solve +\(y^2+y=b/a^2\). The Artin-Schreier map has kernel \(\{0,1\}\) and its +128-element image consists exactly of trace-zero field elements. A private +512-byte immutable table stores a representative for each solvable value and +`-1` for insoluble values, keeping solvable zero distinct from failure. It is +generated once, thread-safely, using native `MultiplyCantor`; there is no global +mutable field setup or large multiplication table. The roots are \(ay\) and +\(ay+a\), and their magnitudes are + +\[ + e_1=(S_1+S_0x_2)/a,\qquad e_2=S_0+e_1. +\] + +The implementation handles every branch explicitly: + +- All four moments zero: already a codeword, not proof of original content. +- Nonzero syndrome with `D=0`: require `S0!=0`, propose `e=S0`, `x=S1/S0`, + and verify all four moments. Rank-one-looking but inconsistent moments fail. +- `D!=0` with `a=0`: reject the repeated-root locator; square roots cannot + produce two distinct error positions. +- Insoluble quadratic or zero proposed magnitude: reject. +- For either candidate size, XOR its contribution out of all four moments and + require zero before reporting success. No input bytes are modified here. + +All byte-valued roots belong to this full mother code; there are no omitted +positions. That fact is specific to the retained full-length R=4 scope. The +existing R=2 virtual-position rejection remains unchanged. Distance five makes +a candidate within radius two unique, but over-radius received words can still +lie within radius two of a different codeword. The direct and generic paths +preserve that BDD behavior; neither promises detection of all larger errors. + +### Integration and Coverage + +`WeakCandidateR4` in `src/reed_solomon/strong_weak_rs_product_code.cc` performs +three native table products per received symbol to accumulate four moments. +It performs no weak transposition or full candidate copy, and final weak +validation requests moments only, without locator solving. R=4 encoding uses +the existing `LCHEncoder(252,4)` on each information row before the unchanged +strong encode; no closed-form encoder speedup is claimed. + +The scheduler checks **every** proposed byte against the optional binary-image +gate (`popcount(delta)<=2`) and its protected column before committing **any** +byte in that row. Rejection is transactional for the entire candidate. Accepted +writes update exact bit/symbol counters, invalidate intersecting column validity, +and activate those columns. Generic reference scheduling now accepts one or two +weak repairs for R=4 while R=2 still accepts only one. + +Tests cover 65,280 single errors; all 32,640 position pairs with magnitudes +`(1,1)` and `(3,128)`; every normalized Artin-Schreier input, both soluble and +insoluble; explicit inconsistent rank-one and repeated-root cases; and 4,096 +random one-through-nine-injection comparisons against generic BDD. Exhaustive +position tests use independently scalar-encoded nonzero rows as their oracle. +All Artin-Schreier and randomized cases compare generic status, count, and full +output. Product scheduler differential coverage includes strong 4/2 and 256/224 +with weak 256/252, all gate combinations, caps 2/3/4/5/6/16, generic single and +batch paths, public direct correction, and every retained optimization bitset. +It compares final bytes, independent component validity, and every result field. +Dedicated tests reject both repairs if only one delta exceeds two bits or only +one target is protected, including data/parity targets. + +Native CLI dimension validation/help, atomic summary metadata, report/replay, +and plotting accept the new dimensions without changing schema or old defaults. +R=2 and R=4 reports never pool into one code group. Random-message native trials +already route nondefault dimensions through `code->Encode`, so the legacy +default-only R=2 encoder helper requires no change. Extended tests check this +route against all-zero trials, all gate combinations and saved-position replay. +Floyd remains thread-count reproducible; Fisher-Yates intentionally retains +worker-local permutations and is tested by saved-position replay instead. + +```cpp +gf2p8::rs::StrongWeakRSProductCode code(256, 224, 256, 252); +// code.Encode(block), then code.Correct(block, options), as before. +``` + +```bash +rs-product-monte-carlo --output /tmp/kilo/mc-r4 --n2 256 --k2 252 --seed 42 --batches 1 --batch-size 64 --threads 1 +``` + +### Matched R=4 Measurements + +AMD Ryzen 7 8845HS, WSL2, GCC 13.3.0, Google Benchmark 1.9.4, Release +`-O3 -DNDEBUG`. Native uses `-march=native`; AVX2-only disables GFNI and AVX-512. +Runs were CPU-0 pinned, sequential, three randomly interleaved repetitions of +at least 0.3s, with `CPP_JOBS=2`. Process inspection showed active editor indexing +and agent processes, but no Monte Carlo worker. Nothing was killed or re-pinned. +Host load and clocks are uncontrolled; these are bounded local measurements, +not Intel-host or multithread scaling evidence. + +Both rows use the same 64-block BSC(0.005) corpus, seed `0x5b5c0224`, strong +256/224, weak 256/252, cap 16, both gates enabled. Corpus setup compares every +output byte and result field with generic scalar/full-validation scheduling. +Timing includes correction and final validity, excludes encode/reset/accounting, +and normalizes to **224*252 = 56,448 information bytes per block**. + +| Profile | Generic batch CPU ms / 64 blocks | Direct CPU ms / 64 blocks | Generic MiB/s | Direct MiB/s | Throughput gain | Generic / direct CPU CV | +| --- | ---: | ---: | ---: | ---: | ---: | ---: | +| Native | 20.290 | 17.756 | 169.801 | 194.033 | 14.3% | 0.48% / 2.30% | +| AVX2-only | 23.501 | 20.852 | 146.600 | 165.230 | 12.7% | 3.38% / 0.34% | + +All measured corpus residual bits and message/block/validity/valid-wrong/pass-limit +failures were zero. Both paths had mean 3.734375 passes, 262.218750 strong lines, +258.765625 weak lines, 2570.328125 byte writes and 2615.734375 bit toggles; +weak-only writes/toggles were 103.453125/105.453125. These are corpus results, +not error-floor estimates. R=2 numbers elsewhere describe a different code and +are not the denominator of these speedup claims. No new instruction-count, +assembly-model, or isolated arithmetic throughput claim is made. + +Reproduction, with exported `CPP_JOBS=2 CPP_BENCH_CPU=0`, +`BENCHMARK_OUT_FORMAT=json` and distinct `BENCHMARK_OUT` paths: + +```bash +/home/user/.config/kilo/scripts/cpp-bench native matrix '^LCH/Owned/StrongWeakRSProductCode/R4/(Direct|Generic)256252$' 3 0.3s +/home/user/.config/kilo/scripts/cpp-bench avx2 matrix '^LCH/Owned/StrongWeakRSProductCode/R4/(Direct|Generic)256252$' 3 0.3s +``` + +Artifacts: `/tmp/kilo/product-r4-native.json` and +`/tmp/kilo/product-r4-avx2.json`. `/check quick` passed portable and native +unfiltered CTest suites; `/check avx2` passed its unfiltered suite. Focused +ASan/UBSan `/check sanitize` with +`GTEST_FILTER='ProductCode.*:WholeCodewordBatch.*'` passed (204.40s selected RS +tests, 268.02s total including CLI/legacy/plot integration). This sanitizer run +excludes unrelated GTest suites and is not described as a full sanitizer pass. +`/check cli` passed all seven tests, including the unrelated fragmenter CLI. +`/check fallback` passed; because that wrapper omits the product source, these +additional strict compilations also passed: + +```bash +g++ -std=c++20 -O2 -Wall -Wextra -Wpedantic -Werror -I include -I src -mno-avx2 -mno-ssse3 -mno-gfni -mno-avx512f -mno-avx512bw -c src/reed_solomon/strong_weak_rs_product_code.cc -o /tmp/kilo/product-r4-scalar.o +g++ -std=c++20 -O2 -Wall -Wextra -Wpedantic -Werror -I include -I src -mavx2 -mssse3 -mno-gfni -mno-avx512f -mno-avx512bw -c src/reed_solomon/strong_weak_rs_product_code.cc -o /tmp/kilo/product-r4-avx2.o +``` + +Formatting uses clang-format 18.1.3; `git diff --check` passed. No commits made. + +## Post-Profile Experiments (2026-09-10) + +All three requested experiments were implemented, differentially tested and +measured. Retained: sparse strong-mask traversal, coordinate-bit weak syndrome +reduction, and padding independent strong batch lanes. Public APIs, gates, +transactionality, stopping rules, counters, mother codes and R=2 policy are +unchanged. The earlier uncommitted direct-R2 implementation and the report below +are preserved. No commits were made. + +### Retained Techniques and Limits + +1. **Strong masks:** AVX2 compares 32 mask bytes against zero and extracts a + bitset. Set-bit traversal visits only verified repairs; old/new XOR popcount, + byte write counts, next-direction activation and clean-state invalidation + still occur for every committed byte. Real partial tails use the exact scalar + scan. Runtime availability is checked once per product correction, not once + per position row. Staging remains transactional per lane; failed lanes are + not committed. This is a modest, workload-dependent improvement, not removal + of the entire old 19-23% profiled region. +2. **Weak reduction:** for each complete 32-byte group, XOR symbols into eight + accumulators selected by the coordinate bits of weight `(j+2) XOR 1`. + Linearity gives `S0 = XOR_b MultiplyCantor(1<