status-im
/
leopard
mirror of https://github.com/status-im/leopard.git
2
0
Fork 0
Code Issues Projects Releases Wiki Activity
leopard/tests/experiments.cpp

621 lines
17 KiB
C++
Raw Blame History

/*
Copyright (c) 2017 Christopher A. Taylor. All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
* Neither the name of Leopard-RS nor the names of its contributors may be
used to endorse or promote products derived from this software without
specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
*/
#include <string.h>
#include <time.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#define LEO_SHORT_FIELD
//#define LEO_EXPERIMENT_EXTRA_XOR
//#define LEO_EXPERIMENT_EXTRA_MULS
#define LEO_EXPERIMENT_CANTOR_BASIS
//------------------------------------------------------------------------------
// Debug
// Some bugs only repro in release mode, so this can be helpful
//#define LEO_DEBUG_IN_RELEASE
#if defined(_DEBUG) || defined(DEBUG) || defined(LEO_DEBUG_IN_RELEASE)
#define LEO_DEBUG
#ifdef _WIN32
#define LEO_DEBUG_BREAK __debugbreak()
#else
#define LEO_DEBUG_BREAK __builtin_trap()
#endif
#define LEO_DEBUG_ASSERT(cond) { if (!(cond)) { LEO_DEBUG_BREAK; } }
#else
#define LEO_DEBUG_BREAK ;
#define LEO_DEBUG_ASSERT(cond) ;
#endif
//------------------------------------------------------------------------------
// Platform/Architecture
// Compiler-specific C++11 restrict keyword
#define LEO_RESTRICT __restrict
// Compiler-specific force inline keyword
#ifdef _MSC_VER
#define LEO_FORCE_INLINE inline __forceinline
#else
#define LEO_FORCE_INLINE inline __attribute__((always_inline))
#endif
//------------------------------------------------------------------------------
// Field
#ifdef LEO_SHORT_FIELD
typedef uint8_t ffe_t;
static const unsigned kGFBits = 8;
static const unsigned kGFPolynomial = 0x11D;
ffe_t kGFBasis[kGFBits] = {
#ifdef LEO_EXPERIMENT_CANTOR_BASIS
1, 214, 152, 146, 86, 200, 88, 230 // Cantor basis
#else
1, 2, 4, 8, 16, 32, 64, 128 // Monomial basis
#endif
};
#else
typedef uint16_t ffe_t;
static const unsigned kGFBits = 16;
static const unsigned kGFPolynomial = 0x1002D;
ffe_t kGFBasis[kGFBits] = {
#ifdef LEO_EXPERIMENT_CANTOR_BASIS
0x0001, 0xACCA, 0x3C0E, 0x163E, // Cantor basis
0xC582, 0xED2E, 0x914C, 0x4012,
0x6C98, 0x10D8, 0x6A72, 0xB900,
0xFDB8, 0xFB34, 0xFF38, 0x991E
#else
1, 2, 4, 8, // Monomial basis
16, 32, 64, 128,
256, 512, 1024, 2048,
4096, 8192, 16384, 32768
#endif
};
#endif
/*
Cantor Basis introduced by:
D. G. Cantor, "On arithmetical algorithms over finite fields",
Journal of Combinatorial Theory, Series A, vol. 50, no. 2, pp. 285-300, 1989.
*/
static const unsigned kFieldSize = (unsigned)1 << kGFBits; //Field size
static const unsigned kModulus = kFieldSize - 1;
static ffe_t GFLog[kFieldSize];
static ffe_t GFExp[kFieldSize];
// Initialize GFLog[], GFExp[]
static void InitField()
{
unsigned state = 1;
for (unsigned i = 0; i < kModulus; ++i)
{
GFExp[state] = static_cast<ffe_t>(i);
state <<= 1;
if (state >= kFieldSize)
state ^= kGFPolynomial;
}
GFExp[0] = kModulus;
// Conversion to chosen basis:
GFLog[0] = 0;
for (unsigned i = 0; i < kGFBits; ++i)
{
const ffe_t basis = kGFBasis[i];
const unsigned width = (unsigned)(1UL << i);
for (unsigned j = 0; j < width; ++j)
GFLog[j + width] = GFLog[j] ^ basis;
}
for (unsigned i = 0; i < kFieldSize; ++i)
GFLog[i] = GFExp[GFLog[i]];
for (unsigned i = 0; i < kFieldSize; ++i)
GFExp[GFLog[i]] = i;
GFExp[kModulus] = GFExp[0];
}
//------------------------------------------------------------------------------
// Mod Q Field Operations
//
// Q is the maximum symbol value, e.g. 255 or 65535.
// z = x + y (mod Q)
static inline ffe_t AddModQ(ffe_t a, ffe_t b)
{
const unsigned sum = (unsigned)a + b;
// Partial reduction step, allowing for Q to be returned
return static_cast<ffe_t>(sum + (sum >> kGFBits));
}
// z = x - y (mod Q)
static inline ffe_t SubModQ(ffe_t a, ffe_t b)
{
const unsigned dif = (unsigned)a - b;
// Partial reduction step, allowing for Q to be returned
return static_cast<ffe_t>(dif + (dif >> kGFBits));
}
// return a*GFExp[b] over GF(2^r)
static ffe_t mulE(ffe_t a, ffe_t b)
{
if (a == 0)
return 0;
// if (b == 0)
// return a;
const ffe_t sum = static_cast<ffe_t>(AddModQ(GFLog[a], b));
return GFExp[sum];
}
//------------------------------------------------------------------------------
// Fast Walsh-Hadamard Transform (FWHT) Mod Q
//
// Q is the maximum symbol value, e.g. 255 or 65535.
// Define this to enable the optimized version of FWHT()
#define LEO_FWHT_OPTIMIZED
typedef ffe_t fwht_t;
// {a, b} = {a + b, a - b} (Mod Q)
static LEO_FORCE_INLINE void FWHT_2(fwht_t& LEO_RESTRICT a, fwht_t& LEO_RESTRICT b)
{
const fwht_t sum = AddModQ(a, b);
const fwht_t dif = SubModQ(a, b);
a = sum;
b = dif;
}
// Reference implementation
static void FWHT(fwht_t* data, const unsigned bits)
{
const unsigned size = (unsigned)(1UL << bits);
for (unsigned width = 1; width < size; width <<= 1)
for (unsigned i = 0; i < size; i += (width << 1))
for (unsigned j = i; j < (width + i); ++j)
FWHT_2(data[j], data[j + width]);
}
//------------------------------------------------------------------------------
// Formal Derivative
// Formal derivative of polynomial in the new basis
static void formal_derivative(ffe_t* cos, const unsigned size)
{
/*
Left to right xoring data ahead into data behind.
If the data ends in all zeroes, this can simply stop.
*/
for (unsigned i = 1; i < size; ++i)
{
const unsigned leng = ((i ^ (i - 1)) + 1) >> 1;
// If a large number of values are being XORed:
for (unsigned j = i - leng; j < i; ++j)
cos[j] ^= cos[j + leng];
}
// Doesn't seem to be needed
#ifdef LEO_EXPERIMENT_EXTRA_XOR
/*
Same here - Zeroes on the right are preserved
*/
for (unsigned i = size; i < kFieldSize; i <<= 1)
{
for (unsigned j = 0; j < size; ++j)
cos[j] ^= cos[j + i];
}
#endif
}
//------------------------------------------------------------------------------
// Fast Fourier Transform
static ffe_t skewVec[kModulus]; // twisted factors used in FFT
static LEO_FORCE_INLINE void ifft_butterfly(ffe_t& a, ffe_t& b, ffe_t skew)
{
b ^= a;
if (skew != kModulus)
a ^= mulE(b, skew);
}
// IFFT in the proposed basis
static void IFLT(ffe_t* data, const unsigned size, const unsigned index)
{
for (unsigned width = 1; width < size; width <<= 1)
{
for (unsigned j = width; j < size; j += (width << 1))
{
const ffe_t skew = skewVec[j + index - 1];
for (unsigned i = j - width; i < j; ++i)
ifft_butterfly(data[i], data[i + width], skew);
}
}
}
static LEO_FORCE_INLINE void fft_butterfly(ffe_t& a, ffe_t& b, ffe_t skew)
{
if (skew != kModulus)
a ^= mulE(b, skew);
b ^= a;
}
// FFT in the proposed basis
static void FLT(ffe_t* data, const unsigned size, const unsigned skewIndex, const unsigned output_elements)
{
for (unsigned width = (size >> 1); width > 0; width >>= 1)
{
const ffe_t* skewLUT = skewVec + width + skewIndex - 1;
for (unsigned j = 0; j < output_elements; j += (width << 1))
{
const ffe_t skew = skewLUT[j];
for (unsigned i = j; i < j + width; ++i)
fft_butterfly(data[i], data[i + width], skew);
}
}
}
//------------------------------------------------------------------------------
// FFT Initialization
#ifdef LEO_EXPERIMENT_EXTRA_MULS
static ffe_t B[kFieldSize >> 1]; // factors used in formal derivative
#endif
static fwht_t log_walsh[kFieldSize]; // factors used in the evaluation of the error locator polynomial
// Initialize skewVec[], B[], log_walsh[]
static void InitFieldOperations()
{
ffe_t temp[kGFBits - 1];
for (unsigned i = 1; i < kGFBits; ++i)
temp[i - 1] = (ffe_t)((unsigned)1 << i);
for (unsigned m = 0; m < (kGFBits - 1); ++m)
{
const unsigned step = (unsigned)1 << (m + 1);
skewVec[((unsigned)1 << m) - 1] = 0;
for (unsigned i = m; i < (kGFBits - 1); ++i)
{
const unsigned s = ((unsigned)1 << (i + 1));
for (unsigned j = ((unsigned)1 << m) - 1; j < s; j += step)
skewVec[j + s] = skewVec[j] ^ temp[i];
}
temp[m] = kModulus - GFLog[mulE(temp[m], GFLog[temp[m] ^ 1])];
for (unsigned i = m + 1; i < (kGFBits - 1); ++i)
{
const ffe_t sum = AddModQ(GFLog[temp[i] ^ 1], temp[m]);
temp[i] = mulE(temp[i], sum);
}
}
for (unsigned i = 0; i < kFieldSize; ++i)
skewVec[i] = GFLog[skewVec[i]];
#ifdef LEO_EXPERIMENT_EXTRA_MULS
temp[0] = kModulus - temp[0];
for (unsigned i = 1; i < (kGFBits - 1); ++i)
temp[i] = (kModulus - temp[i] + temp[i - 1]) % kModulus;
B[0] = 0;
for (unsigned i = 0; i < (kGFBits - 1); ++i)
{
const unsigned depart = ((unsigned)1 << i);
for (unsigned j = 0; j < depart; ++j)
B[j + depart] = (B[j] + temp[i]) % kModulus;
}
#endif
for (unsigned i = 0; i < kFieldSize; ++i)
log_walsh[i] = GFLog[i];
log_walsh[0] = 0;
FWHT(log_walsh, kGFBits);
}
//------------------------------------------------------------------------------
// Encoder
// Encoding alg for k/n<0.5: message is a power of two
static void encodeL(ffe_t* data, const unsigned k, ffe_t* codeword)
{
memcpy(codeword, data, sizeof(ffe_t) * k);
IFLT(codeword, k, 0);
for (unsigned i = k; i < kFieldSize; i += k)
{
memcpy(&codeword[i], codeword, sizeof(ffe_t) * k);
FLT(&codeword[i], k, i, k);
}
memcpy(codeword, data, sizeof(ffe_t) * k);
}
// Encoding alg for k/n>0.5: parity is a power of two.
// data: message array. parity: parity array. mem: buffer(size>= n-k)
static void encodeH(const ffe_t* data, const unsigned m, const unsigned original_count, ffe_t* parity, ffe_t* mem)
{
// Note: Assumes data is padded with zeroes out to the next multiple of m
memcpy(parity, data, m * sizeof(ffe_t));
IFLT(parity, m, m);
for (unsigned i = m; i < original_count; i += m)
{
memcpy(mem, data + i, m * sizeof(ffe_t));
IFLT(mem, m, m + i);
for (unsigned j = 0; j < m; ++j)
parity[j] ^= mem[j];
}
FLT(parity, m, 0, m);
}
//------------------------------------------------------------------------------
// Decoder
static void decode(ffe_t* codeword, const unsigned m, const unsigned original_count, const unsigned n, const bool* erasure)
{
fwht_t log_walsh2[kFieldSize];
// Compute the evaluations of the error locator polynomial
for (unsigned i = 0; i < kFieldSize; ++i)
log_walsh2[i] = erasure[i] ? 1 : 0;
FWHT(log_walsh2, kGFBits);
for (unsigned i = 0; i < kFieldSize; ++i)
log_walsh2[i] = ((unsigned)log_walsh2[i] * (unsigned)log_walsh[i]) % kModulus;
FWHT(log_walsh2, kGFBits);
// k2 can be replaced with k
//const unsigned k2 = kFieldSize;
//const unsigned k2 = k; // cannot actually be replaced with k. maybe for encodeL() only?
for (unsigned i = 0; i < m + original_count; ++i)
{
if (erasure[i])
{
codeword[i] = 0;
}
else
{
codeword[i] = mulE(codeword[i], log_walsh2[i]);
}
}
for (unsigned i = m + original_count; i < n; ++i)
codeword[i] = 0;
IFLT(codeword, n, 0);
// Note: This is not needed to recover successfully...
#ifdef LEO_EXPERIMENT_EXTRA_MULS
// formal derivative
// Note: Preserves zeroes on the right
for (unsigned i = 0; i < m + original_count; i += 2)
{
codeword[i] = mulE(codeword[i], kModulus - B[i >> 1]);
codeword[i + 1] = mulE(codeword[i + 1], kModulus - B[i >> 1]);
}
#endif
formal_derivative(codeword, n);
#ifdef LEO_EXPERIMENT_EXTRA_MULS
// Note: Preserves zeroes on the right
for (unsigned i = 0; i < m + original_count; i += 2)
{
codeword[i] = mulE(codeword[i], B[i >> 1]);
codeword[i + 1] = mulE(codeword[i + 1], B[i >> 1]);
}
#endif
FLT(codeword, n, 0, m + original_count);
for (unsigned i = 0; i < kFieldSize; ++i)
{
if (erasure[i])
{
codeword[i] = mulE(codeword[i], kModulus - log_walsh2[i]);
}
}
}
#ifdef _MSC_VER
#include <intrin.h>
#endif
// Returns highest bit index 0..63 where the first non-zero bit is found
// Precondition: x != 0
LEO_FORCE_INLINE unsigned LastNonzeroBit64(uint64_t x)
{
#ifdef _MSC_VER
#ifdef _WIN64
unsigned long index;
// Note: Ignoring result because x != 0
_BitScanReverse64(&index, x);
return (unsigned)index;
#else
unsigned long index;
if (0 != _BitScanReverse(&index, (uint32_t)x))
return (unsigned)index;
// Note: Ignoring result because x != 0
_BitScanReverse(&index, (uint32_t)(x >> 32));
return (unsigned)index + 32;
#endif
#else
// Note: Ignoring return value of 0 because x != 0
return 63 - (unsigned)__builtin_clzll(x);
#endif
}
//------------------------------------------------------------------------------
// Test Application
void test(unsigned original_count, unsigned recovery_count, unsigned seed)
{
unsigned m = 2UL << LastNonzeroBit64(recovery_count - 1);
unsigned n = 2UL << LastNonzeroBit64(m + original_count - 1);
srand(seed);
//-----------Generating message----------
// Message array
ffe_t data[kFieldSize] = {0};
// Filled with random numbers
for (unsigned i = m; i < m + original_count; ++i)
data[i] = (ffe_t)rand();
//---------encoding----------
ffe_t codeword[kFieldSize] = {};
// First m codewords are for the parity data
encodeH(data + m, m, original_count, data, codeword);
//encodeL(data, k, codeword); // does not seem to work with any input? what else needs to change?
memcpy(codeword, data, sizeof(ffe_t) * kFieldSize);
//--------erasure simulation---------
// Array indicating erasures
bool erasure[kFieldSize] = {
false
};
// Tag the first "recovery_count" elements as erasures
for (unsigned i = m; i < m + recovery_count; ++i)
erasure[i] = true;
// permuting the erasure array
for (unsigned i = m + original_count - 1; i > 0; --i)
{
unsigned pos = rand() % (i + 1);
if (i != pos)
{
bool tmp = erasure[i];
erasure[i] = erasure[pos];
erasure[pos] = tmp;
}
}
//---------main processing----------
decode(codeword, m, original_count, n, erasure);
// Check the correctness of the result
for (unsigned i = 0; i < kFieldSize; ++i)
{
if (erasure[i])
{
if (data[i] != codeword[i])
{
printf("Decoding Error with seed = %d!\n", seed);
LEO_DEBUG_BREAK;
return;
}
}
}
printf(":D ");
}
//------------------------------------------------------------------------------
// Entrypoint
int main(int argc, char **argv)
{
// Fill GFLog table and GFExp table
InitField();
// Compute factors used in erasure decoder
InitFieldOperations();
unsigned seed = (unsigned)time(NULL);
for (;;)
{
#ifdef LEO_SHORT_FIELD
const unsigned input_count = 100;
const unsigned recovery_count = 10;
#else // LEO_SHORT_FIELD
const unsigned input_count = 10000;
const unsigned recovery_count = 2000;
#endif // LEO_SHORT_FIELD
test(input_count, recovery_count, seed);
++seed;
}
return 0;
}
Reference in New Issue View Git Blame Copy Permalink