/* 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 #include #include #include #include //#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(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(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(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; const ffe_t sum = static_cast(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; 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]; if (skew != kModulus) { for (unsigned i = j - width; i < j; ++i) ifft_butterfly(data[i], data[i + width], skew); } else { for (unsigned i = j - width; i < j; ++i) data[i + width] ^= data[i]; } } } } static LEO_FORCE_INLINE void fft_butterfly(ffe_t& a, ffe_t& b, ffe_t skew) { 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]; if (skew != kModulus) { for (unsigned i = j; i < j + width; ++i) fft_butterfly(data[i], data[i + width], skew); } else { for (unsigned i = j; i < j + width; ++i) data[i + width] ^= data[i]; } } } } //------------------------------------------------------------------------------ // 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) temp[i] = mulE(temp[i], (GFLog[temp[i] ^ 1] + temp[m]) % kModulus); } 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 #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; }