Strong-weak RS product codes

CMakeLists.txt

@@ -62,6 +62,7 @@ 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)
 target_include_directories(gf256_core
@@ -142,7 +143,8 @@ if(GF256_BUILD_TESTS)
   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")
@@ -170,8 +172,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

benchmarks/strong_weak_rs_product_code_benchmarks.cc

@@ -0,0 +1,199 @@
+#include <algorithm>
+#include <bit>
+#include <cstdint>
+#include <random>
+#include <vector>
+
+#include "benchmark/benchmark.h"
+#include "reed_solomon/strong_weak_rs_product_code.h"
+#include "reed_solomon/product_code_internal.h"
+
+namespace {
+
+void BenchmarkProductCorrectionBSC(benchmark::State& state, int batch_passes = -1) {
+  using gf2p8::Element;
+  using gf2p8::rs::ProductCorrectionResult;
+  using gf2p8::rs::ProductTermination;
+  constexpr size_t kN = 256;
+  constexpr size_t kStrongK = 224;
+  constexpr size_t kWeakK = 254;
+  constexpr size_t kInformationBytes = kStrongK * kWeakK;
+  constexpr size_t kBlockBytes = kN * kN;
+  constexpr size_t kCorpusCount = 64;
+  constexpr size_t kPassLimit = 16;
+  constexpr uint32_t kSeed = 0x5b5c0224;
+  gf2p8::rs::StrongWeakRSProductCode code(kN, kStrongK, kN, kWeakK);
+  if (!code.Valid() || code.BlockSize() != kBlockBytes) {
+    state.SkipWithError("invalid product-code dimensions");
+    return;
+  }
+
+  std::mt19937 messages(kSeed);
+  std::mt19937 channel(kSeed ^ 0x9e3779b9U);
+  std::vector<std::vector<Element>> original(
+      kCorpusCount, std::vector<Element>(kBlockBytes));
+  auto corrupted = original;
+  auto work = original;
+  std::vector<ProductCorrectionResult> 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 * kN + col] = static_cast<Element>(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<uint32_t>(channel());
+        } while (draw >= 4294967200U);
+        if (draw % 200 == 0) {
+          value ^= static_cast<Element>(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;
+  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 = batch_passes < 0 ? code.Correct(work[sample], kPassLimit)
+        : gf2p8::rs::detail::ProductCorrectionAccess::Correct(
+              code, work[sample], kPassLimit, batch_passes);
+    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) {
+      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] = batch_passes < 0 ? code.Correct(work[sample], kPassLimit)
+          : gf2p8::rs::detail::ProductCorrectionAccess::Correct(
+                code, work[sample], kPassLimit, batch_passes);
+      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<unsigned>(
+            work[sample][pos] ^ original[sample][pos]));
+        block_bits += bits;
+        if (pos / kN < kStrongK && pos % kN < 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;
+    }
+    state.ResumeTiming();
+    if (state.skipped()) break;
+  }
+
+  const double sweeps = static_cast<double>(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<double>(channel_bits) / (kCorpusCount * kBlockBytes * 8);
+  state.counters["channel_flipped_bits"] = static_cast<double>(channel_bits);
+  state.counters["channel_corrupted_bytes"] = static_cast<double>(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<double>(max_passes);
+  state.counters["mean_strong_lines"] = strong_lines / blocks;
+  state.counters["mean_weak_lines"] = weak_lines / blocks;
+  state.counters["pass_limit"] = kPassLimit;
+  state.SetItemsProcessed(state.iterations() * kCorpusCount);
+  state.SetBytesProcessed(state.iterations() * kCorpusCount * kInformationBytes);
+  state.SetLabel("64 blocks/iteration; information bytes=224*254; "
+                  "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* kProductCorrectionSingle = benchmark::RegisterBenchmark(
+    "LCH/Owned/StrongWeakRSProductCode/Correct/BSC005/Single",
+    [](benchmark::State& state) { BenchmarkProductCorrectionBSC(state, 0); });
+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); });
+
+}  // namespace

include/reed_solomon/error_correction.h

@@ -0,0 +1,35 @@
+#pragma once
+
+#include "reed_solomon/lch_decoder.h"
+
+namespace gf2p8::rs {
+
+/** @brief Bounded-distance correction outcome. */
+enum class CorrectionStatus {
+  ok,
+  invalid_argument,
+  unsupported_dimensions,
+  uncorrectable,
+  reconstruction_failed,
+};
+
+/** @brief Status and number of repaired symbols (zero on failure). */
+struct CorrectionResult {
+  CorrectionStatus status = CorrectionStatus::invalid_argument;
+  size_t error_count = 0;
+};
+
+/**
+ * @brief Transactionally repairs data AND parity of one Cantor RS codeword.
+ * @param decoder Decoder defining K data and R recovery symbols.
+ * @param codeword Exactly K+R mutable symbols in [data][recovery] order.
+ * @return Status and corrected symbol count; zero count on an already valid word.
+ * @details Requires unshortened power-of-two N=K+R <=256 and power-of-two
+ * R<=K. Corrects up to floor(R/2) unknown symbol errors, with no erasures.
+ * On any non-ok status every input byte is unchanged. Success means the result
+ * is a codeword within that radius, NOT necessarily the transmitted codeword.
+ */
+CorrectionResult CorrectCodeword(const LCHDecoder& decoder,
+                                std::span<Element> codeword);
+
+}  // namespace gf2p8::rs

include/reed_solomon/strong_weak_rs_product_code.h

@@ -0,0 +1,91 @@
+#pragma once
+
+#include "reed_solomon/error_correction.h"
+#include "reed_solomon/lch_encoder.h"
+
+namespace gf2p8::rs {
+
+namespace detail { struct ProductCorrectionAccess; }
+
+/** @brief Reason iterative product correction stopped, independent of validity. */
+enum class ProductTermination { invalid_argument, no_change, pass_limit };
+
+/** @brief Product correction outcome; validity does not prove original content. */
+struct ProductCorrectionResult {
+  ProductTermination termination = ProductTermination::invalid_argument;
+  bool all_zero_syndromes = false;
+  size_t directional_passes = 0;
+  size_t strong_lines_visited = 0;
+  size_t weak_lines_visited = 0;
+  /** @brief Accepted symbol writes, counting repeated repairs separately. */
+  size_t changed_symbols = 0;
+};
+
+/**
+ * @brief Systematic Cantor RS product with strong columns and weak rows.
+ * @details Row-major block has Nstrong rows and Nweak columns. The top-left
+ * Kstrong by Kweak rectangle holds data; every column and every row, including
+ * parity regions, is a component codeword. Errors only: no erasures,
+ * shortening, or backtracking. Both N must be powers of two <=256; strong R
+ * must be a power of two with 2<=R<=K, and weak R must equal 2.
+ */
+class StrongWeakRSProductCode {
+ public:
+  /**
+   * @brief Constructs a product code; unsupported dimensions make Valid false.
+   * @param strong_n Number of rows.
+   * @param strong_k Number of systematic rows.
+   * @param weak_n Number of columns.
+   * @param weak_k Number of systematic columns.
+   */
+  StrongWeakRSProductCode(size_t strong_n = 256, size_t strong_k = 224,
+                          size_t weak_n = 256, size_t weak_k = 254);
+
+  /** @brief Reports supported dimensions. @return Whether operations are valid. */
+  bool Valid() const;
+  /** @brief Returns required row-major block size. @return Bytes, or zero if invalid. */
+  size_t BlockSize() const;
+
+  /**
+   * @brief Fills all product parity, preserving the systematic rectangle.
+   * @param block Exactly BlockSize() symbols, with data already in place.
+   * @param backend Owned encoder backend, scalar available for reference checks.
+   * @return Encoding status; on failure block is unchanged.
+   */
+  lch::Status Encode(std::span<Element> block,
+                     lch::Backend backend = lch::Backend::tuned) const;
+
+  /**
+   * @brief Applies alternating strong BDD and gated weak single-error passes.
+   * @param block Exactly BlockSize() mutable symbols.
+   * @param max_directional_passes Cap counting each direction separately; >=2.
+   * @return Termination, pass/work counts, and final all-component validity.
+   * @details Always starts with a full strong pass then a full weak pass, even
+   * if strong changes nothing. From then on visits only lines intersecting
+   * accepted byte changes in the preceding pass. Stops on a no-change pass
+   * (including the initial weak pass), or the cap; no-change takes precedence.
+   * Strong success, including zero syndromes, protects that column; failure
+   * leaves its bytes unchanged and unprotects it. Unvisited protection persists.
+   * Weak candidates commit only for one changed symbol with popcount(old XOR
+   * new)<=2 in an unprotected column; rejected rows are unchanged. Protection
+   * changes alone never activate lines. Accepted repairs are retained at exit;
+   * the entire iterative operation is not transactional. Invalid arguments
+   * leave the block untouched. Validity is checked separately at exit and does
+   * not count as a directional pass or guarantee the original message.
+   */
+  ProductCorrectionResult Correct(std::span<Element> block,
+                                   size_t max_directional_passes = 16) const;
+
+ private:
+  friend struct detail::ProductCorrectionAccess;
+  ProductCorrectionResult CorrectImpl(std::span<Element> block,
+                                      size_t max_directional_passes,
+                                       unsigned batch_passes,
+                                       bool tracked_validation = true) const;
+  size_t strong_n_, strong_k_, weak_n_, weak_k_;
+  bool valid_;
+  LCHEncoder strong_encoder_, weak_encoder_;
+  LCHDecoder strong_decoder_, weak_decoder_;
+};
+
+}  // namespace gf2p8::rs