mirror of
https://github.com/aicodix/code.git
synced 2026-04-27 14:30:36 +00:00
added Cauchy Reed Solomon erasure coding
Read more about the Cauchy matrix: https://en.wikipedia.org/wiki/Cauchy_matrix The Art of Computer Programming Volume 1: Fundamental Algorithms Donald E. Knuth - 1997 Paragraph 1.2.3: Sums and Products Page 37, Cauchy's matrix: $a_{ij} = \frac{1}{x_i + y_j}$ Page 38, Exercise 41: $b_{ij} = \frac{\prod_{k=1}^{n}{(x_j + y_k)(x_k + y_i)}}{(x_j + y_i)\prod_{k \ne j}^{n}{(x_j - x_k)}\prod_{k \ne i}^{n}{(y_i - y_k)}}$
This commit is contained in:
Ahmet Inan
parent
3fa6dde2a8
commit
eb6e002d52
2 changed files with 174 additions and 0 deletions
104
cauchy_reed_solomon_erasure_coding.hh
Normal file
104
cauchy_reed_solomon_erasure_coding.hh
Normal file
|
|
@ -0,0 +1,104 @@
|
|||
/*
|
||||
Cauchy Reed Solomon Erasure Coding
|
||||
|
||||
Copyright 2023 Ahmet Inan <inan@aicodix.de>
|
||||
*/
|
||||
|
||||
#pragma once
|
||||
|
||||
namespace CODE {
|
||||
|
||||
template <typename GF>
|
||||
struct CauchyReedSolomonEncoder
|
||||
{
|
||||
typedef typename GF::value_type value_type;
|
||||
typedef typename GF::ValueType ValueType;
|
||||
typedef typename GF::IndexType IndexType;
|
||||
// $a_{ij} = \frac{1}{x_i + y_j}$
|
||||
IndexType cauchy_matrix(int i, int j)
|
||||
{
|
||||
ValueType row(i), col(ValueType::N - j);
|
||||
return rcp(index(row + col));
|
||||
}
|
||||
void operator()(const ValueType *data, ValueType *block, int block_num, int block_len, int block_cnt)
|
||||
{
|
||||
assert(block_num <= ValueType::N - block_cnt);
|
||||
for (int k = 0; k < block_cnt; k++) {
|
||||
IndexType a_ik = cauchy_matrix(block_num, k);
|
||||
for (int j = 0; j < block_len; j++) {
|
||||
if (k)
|
||||
block[j] = fma(a_ik, data[block_len*k+j], block[j]);
|
||||
else
|
||||
block[j] = a_ik * data[block_len*k+j];
|
||||
}
|
||||
}
|
||||
}
|
||||
void operator()(const value_type *data, value_type *block, int block_num, int block_len, int block_cnt)
|
||||
{
|
||||
(*this)(reinterpret_cast<const ValueType *>(data), reinterpret_cast<ValueType *>(block), block_num, block_len, block_cnt);
|
||||
}
|
||||
void operator()(const void *data, void *block, int block_number, int block_bytes, int block_count)
|
||||
{
|
||||
assert(block_bytes % sizeof(value_type) == 0);
|
||||
(*this)(reinterpret_cast<const value_type *>(data), reinterpret_cast<value_type *>(block), block_number, block_bytes / sizeof(value_type), block_count);
|
||||
}
|
||||
};
|
||||
|
||||
template <typename GF>
|
||||
struct CauchyReedSolomonDecoder
|
||||
{
|
||||
typedef typename GF::value_type value_type;
|
||||
typedef typename GF::ValueType ValueType;
|
||||
typedef typename GF::IndexType IndexType;
|
||||
// $b_{ij} = \frac{\prod_{k=1}^{n}{(x_j + y_k)(x_k + y_i)}}{(x_j + y_i)\prod_{k \ne j}^{n}{(x_j - x_k)}\prod_{k \ne i}^{n}{(y_i - y_k)}}$
|
||||
IndexType inverse_cauchy_matrix(const ValueType *rows, int i, int j, int n)
|
||||
{
|
||||
ValueType col_i(ValueType::N - i);
|
||||
IndexType prod_xy(0);
|
||||
for (int k = 0; k < n; k++) {
|
||||
ValueType col_k(ValueType::N - k);
|
||||
prod_xy *= index(rows[j] + col_k) * index(rows[k] + col_i);
|
||||
}
|
||||
IndexType prod_x(0);
|
||||
for (int k = 0; k < n; k++) {
|
||||
if (k != j) {
|
||||
prod_x *= index(rows[j] + rows[k]);
|
||||
}
|
||||
}
|
||||
IndexType prod_y(0);
|
||||
for (int k = 0; k < n; k++) {
|
||||
if (k != i) {
|
||||
ValueType col_k(ValueType::N - k);
|
||||
prod_y *= index(col_i + col_k);
|
||||
}
|
||||
}
|
||||
return prod_xy / (index(rows[j] + col_i) * prod_x * prod_y);
|
||||
}
|
||||
void operator()(ValueType *data, const ValueType *blocks, const ValueType *block_nums, int block_len, int block_cnt)
|
||||
{
|
||||
for (int i = 0; i < block_cnt; i++) {
|
||||
for (int k = 0; k < block_cnt; k++) {
|
||||
IndexType b_ik = inverse_cauchy_matrix(block_nums, i, k, block_cnt);
|
||||
for (int j = 0; j < block_len; j++) {
|
||||
if (k)
|
||||
data[j] = fma(b_ik, blocks[block_len*k+j], data[j]);
|
||||
else
|
||||
data[j] = b_ik * blocks[block_len*k+j];
|
||||
}
|
||||
}
|
||||
data += block_len;
|
||||
}
|
||||
}
|
||||
void operator()(value_type *data, const value_type *blocks, const value_type *block_nums, int block_len, int block_cnt)
|
||||
{
|
||||
(*this)(reinterpret_cast<ValueType *>(data), reinterpret_cast<const ValueType *>(blocks), reinterpret_cast<const ValueType *>(block_nums), block_len, block_cnt);
|
||||
}
|
||||
void operator()(void *data, const void *blocks, const value_type *block_numbers, int block_bytes, int block_count)
|
||||
{
|
||||
assert(block_bytes % sizeof(value_type) == 0);
|
||||
(*this)(reinterpret_cast<value_type *>(data), reinterpret_cast<const value_type *>(blocks), block_numbers, block_bytes / sizeof(value_type), block_count);
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
70
tests/crs_regression_test.cc
Normal file
70
tests/crs_regression_test.cc
Normal file
|
|
@ -0,0 +1,70 @@
|
|||
/*
|
||||
Regression Test for the Cauchy Reed Solomon Encoder and Decoder
|
||||
|
||||
Copyright 2023 Ahmet Inan <inan@aicodix.de>
|
||||
*/
|
||||
|
||||
#include <cassert>
|
||||
#include <random>
|
||||
#include <iostream>
|
||||
#include <functional>
|
||||
#include "galois_field.hh"
|
||||
#include "cauchy_reed_solomon_erasure_coding.hh"
|
||||
|
||||
template <typename GF>
|
||||
void crs_test(int trials)
|
||||
{
|
||||
CODE::CauchyReedSolomonEncoder<GF> encode;
|
||||
CODE::CauchyReedSolomonDecoder<GF> decode;
|
||||
std::random_device rd;
|
||||
std::default_random_engine generator(rd());
|
||||
typedef std::uniform_int_distribution<int> distribution;
|
||||
auto rnd_cnt = std::bind(distribution(1, std::min(GF::Q / 2, 256)), generator);
|
||||
auto rnd_len = std::bind(distribution(1, 1 << 10), generator);
|
||||
auto rnd_dat = std::bind(distribution(0, 255), generator);
|
||||
while (--trials) {
|
||||
int block_count = rnd_cnt();
|
||||
int numbers_total = GF::Q - block_count;
|
||||
int block_bytes = rnd_len() * sizeof(typename GF::value_type);
|
||||
int data_bytes = block_count * block_bytes;
|
||||
uint8_t *orig = new uint8_t[data_bytes];
|
||||
for (int i = 0; i < data_bytes; ++i)
|
||||
orig[i] = rnd_dat();
|
||||
uint8_t *blocks = new uint8_t[data_bytes];
|
||||
auto numbers = new typename GF::value_type[numbers_total];
|
||||
for (int i = 0; i < numbers_total; ++i)
|
||||
numbers[i] = i;
|
||||
for (int i = 0; i < block_count; i++) {
|
||||
auto hat = std::bind(distribution(i, numbers_total - 1), generator);
|
||||
std::swap(numbers[i], numbers[hat()]);
|
||||
}
|
||||
for (int i = 0; i < block_count; ++i)
|
||||
encode(orig, blocks + block_bytes * i, numbers[i], block_bytes, block_count);
|
||||
uint8_t *data = new uint8_t[data_bytes];
|
||||
decode(data, blocks, numbers, block_bytes, block_count);
|
||||
for (int i = 0; i < data_bytes; ++i)
|
||||
assert(data[i] == orig[i]);
|
||||
delete[] numbers;
|
||||
delete[] blocks;
|
||||
delete[] orig;
|
||||
delete[] data;
|
||||
}
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
if (1) {
|
||||
typedef CODE::GaloisField<8, 0b100011101, uint8_t> GF;
|
||||
GF instance;
|
||||
crs_test<GF>(200);
|
||||
}
|
||||
if (1) {
|
||||
typedef CODE::GaloisField<16, 0b10001000000001011, uint16_t> GF;
|
||||
GF *instance = new GF();
|
||||
crs_test<GF>(100);
|
||||
delete instance;
|
||||
}
|
||||
std::cerr << "Cauchy Reed Solomon regression test passed!" << std::endl;
|
||||
return 0;
|
||||
}
|
||||
|
||||
Loading…
Add table
Add a link
Reference in a new issue