g4klx/MMDVMHost
1
0
Fork 0
mirror of https://github.com/g4klx/MMDVMHost.git synced 2025-12-06 05:32:00 +01:00
Code Issues Actions Packages Projects Releases Wiki Activity
MMDVMHost/RS.h
Jonathan Naylor cf15f42a0f Convert to C++ nullptr from NULL.
2025-03-14 12:13:11 +00:00

1108 lines
51 KiB
C++
Raw Blame History

/*
* Copyright (C) 2014 Hard Consulting Corporation
* Copyright (C) 2023 Bryan Biedenkapp, N2PLL
* Copyright (C) 2024,2025 Jonathan Naylor, G4KLX
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
#if !defined(RS_H)
#define RS_H
/*
* Ezpwd Reed-Solomon -- Reed-Solomon encoder / decoder library
*
* The core Reed-Solomon codec implementation in c++/ezpwd/rs_base is by Phil Karn, converted to C++
* by Perry Kundert (perry@hardconsulting.com), and may be used under the terms of the LGPL. Here
* is the terms from Phil's README file (see phil-karn/fec-3.0.1/README):
*
* COPYRIGHT
*
* This package is copyright 2006 by Phil Karn, KA9Q. It may be used
* under the terms of the GNU Lesser General Public License (LGPL). See
* the file "lesser.txt" in this package for license details.
*
* The c++/ezpwd/rs_base file is, therefore, redistributed under the terms of the LGPL, while the
* rest of Ezpwd Reed-Solomon is distributed under either the GPL or Commercial licenses.
* Therefore, even if you have obtained Ezpwd Reed-Solomon under a Commercial license, you must make
* available the source code of the c++/ezpwd/rs_base file with your product. One simple way to
* accomplish this is to include the following URL in your code or documentation:
*
* https://github.com/pjkundert/ezpwd-reed-solomon/blob/master/c++/ezpwd/rs_base
*
*
* The Linux 3.15.1 version of lib/reed_solomon was also consulted as a cross-reference, which (in
* turn) is basically verbatim copied from Phil Karn's LGPL implementation, to ensure that no new
* defects had been found and fixed; there were no meaningful changes made to Phil's implementation.
* I've personally been using Phil's implementation for years in a heavy industrial use, and it is
* rock-solid.
*
* However, both Phil's and the Linux kernel's (copy of Phil's) implementation will return a
* "corrected" decoding with impossible error positions, in some cases where the error load
* completely overwhelms the R-S encoding. These cases, when detected, are rejected in this
* implementation. This could be considered a defect in Phil's (and hence the Linux kernel's)
* implementations, which results in them accepting clearly incorrect R-S decoded values as valid
* (corrected) R-S codewords. We chose the report failure on these attempts.
*
*/
#include <algorithm>
#include <array>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <iomanip>
#include <type_traits>
#include <vector>
//
// Preprocessor defines available:
//
// EZPWD_NO_EXCEPTS -- define to use no exceptions; return -1, or abort on catastrophic failures
// EZPWD_NO_MOD_TAB -- define to force no "modnn" Galois modulo table acceleration
// EZPWD_ARRAY_SAFE -- define to force usage of bounds-checked arrays for most tabular data
// EZPWD_ARRAY_TEST -- define to force erroneous sizing of some arrays for non-production testing
//
#if defined(EZPWD_NO_EXCEPTS)
#include <cstdio>
#define EZPWD_RAISE_OR_ABORT(typ, str) do { \
std::fputs((str), stderr); std::fputc('\n', stderr);\
abort(); \
} while (false)
#define EZPWD_RAISE_OR_RETURN(typ, str, ret) return (ret)
#else
#define EZPWD_RAISE_OR_ABORT(typ, str) throw (typ)(str)
#define EZPWD_RAISE_OR_RETURN(typ, str, ret) throw (typ)(str)
#endif
namespace rs
{
// ezpwd::log_<N,B> -- compute the log base B of N at compile-time
template <size_t N, size_t B = 2> struct log_ { enum { value = 1 + log_<N / B, B>::value }; };
template <size_t B> struct log_<1, B> { enum { value = 0 }; };
template <size_t B> struct log_<0, B> { enum { value = 0 }; };
// ---------------------------------------------------------------------------
// Class Declaration
// ---------------------------------------------------------------------------
/**
* @brief Reed-Solomon codec generic base class.
* @ingroup edac_rs
*/
class reed_solomon_base {
public:
/** @brief A data element's bits. */
virtual size_t datum() const = 0;
/** @brief A symbol's bits. */
virtual size_t symbol() const = 0;
/** @brief R-S block size (maximum total symbols). */
virtual int size() const = 0;
/** @brief R-S roots (parity symbols). */
virtual int nroots() const = 0;
/** @brief R-S net payload (data symbols). */
virtual int load() const = 0;
/** @brief Initializes a new instance of the reed_solomon_base class. */
reed_solomon_base() = default;
/** @brief Finalizes a instance of the reed_solomon_base class. */
virtual ~reed_solomon_base() = default;
/** @brief */
virtual std::ostream& output(std::ostream& lhs) const { return lhs << "RS(" << this->size() << "," << this->load() << ")"; }
//
// {en,de}code -- Compute/Correct errors/erasures in a Reed-Solomon encoded container
//
// The encoded parity symbols may be included in 'data' (len includes nroots() parity
// symbols), or may (optionally) supplied separately in (at least nroots()-sized)
// 'parity'.
//
// For decode, optionally specify some known erasure positions (up to nroots()). If
// non-empty 'erasures' is provided, it contains the positions of each erasure. If a
// non-zero pointer to a 'position' vector is provided, its capacity will be increased to
// be capable of storing up to 'nroots()' ints; the actual deduced error locations will be
// returned.
//
// RETURN VALUE
//
// Return -1 on error. The encode returns the number of parity symbols produced;
// decode returns the number of symbols corrected. Both errors and erasures are included,
// so long as they are actually different than the deduced value. In other words, if a
// symbol is marked as an erasure but it actually turns out to be correct, it's index will
// NOT be included in the returned count, nor the modified erasure vector!
//
int encode(std::string& data) const
{
typedef uint8_t uT;
typedef std::pair<uT*, uT*> uTpair;
data.resize(data.size() + nroots());
return encode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()));
}
int encode(const std::string& data, std::string& parity) const
{
typedef uint8_t uT;
typedef std::pair<const uT*, const uT*> cuTpair;
typedef std::pair<uT*, uT*> uTpair;
parity.resize(nroots());
return encode(cuTpair((const uT*)&data.front(), (const uT*)&data.front() + data.size()),
uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()));
}
template<typename T> int encode(std::vector<T>& data) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
data.resize(data.size() + nroots());
return encode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()));
}
template<typename T> int encode(const std::vector<T>& data, std::vector<T>& parity) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<const uT*, const uT*> cuTpair;
typedef std::pair<uT*, uT*> uTpair;
parity.resize(nroots());
return encode(cuTpair((uT*)&data.front(), (uT*)&data.front() + data.size()),
uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()));
}
template<typename T, size_t N> int encode(std::array<T, N>& data, int pad = 0) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
return encode(uTpair((uT*)&data.front() + pad, (uT*)&data.front() + data.size()));
}
virtual int encode(const std::pair<uint8_t*, uint8_t*>& data) const = 0;
virtual int encode(const std::pair<const uint8_t*, const uint8_t*>& data,
const std::pair<uint8_t*, uint8_t*>& parity) const = 0;
virtual int encode(const std::pair<uint16_t*, uint16_t*>& data) const = 0;
virtual int encode(const std::pair<const uint16_t*, const uint16_t*>& data,
const std::pair<uint16_t*, uint16_t*>& parity) const = 0;
virtual int encode(const std::pair<uint32_t*, uint32_t*>& data) const = 0;
virtual int encode(const std::pair<const uint32_t*, const uint32_t*>& data,
const std::pair<uint32_t*, uint32_t*>& parity) const = 0;
int decode(std::string& data, const std::vector<int>& erasure = std::vector<int>(),
std::vector<int>* position = 0) const
{
typedef uint8_t uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()), erasure, position);
}
int decode(std::string& data, std::string& parity, const std::vector<int>& erasure = std::vector<int>(),
std::vector<int>* position = 0) const
{
typedef uint8_t uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()),
uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()), erasure, position);
}
template<typename T> int decode(std::vector<T>& data, const std::vector<int>& erasure = std::vector<int>(),
std::vector<int>* position = 0) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()), erasure, position);
}
template<typename T> int decode(std::vector<T>& data, std::vector<T>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()),
uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()), erasure, position);
}
template<typename T, size_t N> int decode(std::array<T, N>& data, int pad = 0,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
typedef typename std::make_unsigned<T>::type uT;
typedef std::pair<uT*, uT*> uTpair;
return decode(uTpair((uT*)&data.front() + pad, (uT*)&data.front() + data.size()), erasure, position);
}
virtual int decode(const std::pair<uint8_t*, uint8_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint8_t*, uint8_t*>& data, const std::pair<uint8_t*, uint8_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint16_t*, uint16_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint16_t*, uint16_t*>& data, const std::pair<uint16_t*, uint16_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint32_t*, uint32_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
virtual int decode(const std::pair<uint32_t*, uint32_t*>& data, const std::pair<uint32_t*, uint32_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const = 0;
};
//
// std::ostream << edac::rs::reed_solomon<...>
//
// Output a R-S codec description in standard form eg. RS(255,253)
//
inline std::ostream& operator<<(std::ostream& lhs, const rs::reed_solomon_base& rhs) { return rhs.output(lhs); }
// ---------------------------------------------------------------------------
// Structure Declaration
// ---------------------------------------------------------------------------
/**
* @brief Default field polynomial generator functor.
* @tparam SYM
* @tparam PLY
* @ingroup edac_rs
*/
template<int SYM, int PLY>
struct gfpoly {
int operator() (int sr) const
{
if (sr == 0) {
sr = 1;
} else {
sr <<= 1;
if (sr & (1 << SYM))
sr ^= PLY;
sr &= ((1 << SYM) - 1);
}
return sr;
}
};
// ---------------------------------------------------------------------------
// Class Declaration
// ---------------------------------------------------------------------------
/**
* @brief R-S tables common to all RS(NN,*) with same SYM, PRM and PLY.
* @tparam TYP
* @tparam SYM
* @tparam PRM
* @tparam PLY
* @ingroup edac_rs
*/
template <typename TYP, int SYM, int PRM, class PLY>
class reed_solomon_tabs : public reed_solomon_base {
public:
typedef TYP symbol_t;
/** @brief Bits / TYP */
static const size_t DATUM = 8 * sizeof TYP();
/** @brief Bits / Symbol */
static const size_t SYMBOL = SYM;
static const int MM = SYM;
static const int SIZE = (1 << SYM) - 1; // maximum symbols in field
static const int NN = SIZE;
static const int A0 = SIZE;
// modulo table: 1/2 the symbol size squared, up to 4k
static const int MODS = SYM > 8 ? (1 << 12) : (1 << SYM << SYM / 2);
static int iprim; // initialized to -1, below
protected:
static std::array<TYP, NN + 1> alpha_to;
static std::array<TYP, NN + 1> index_of;
static std::array<TYP, MODS> mod_of;
/** @brief Initializes a new instance of the reed_solomon_tabs class. */
reed_solomon_tabs() : reed_solomon_base()
{
// Do init if not already done. We check one value which is initialized to -1; this is
// safe, 'cause the value will not be set 'til the initializing thread has completely
// initialized the structure. Worst case scenario: multiple threads will initialize
// identically. No mutex necessary.
if (iprim >= 0)
return;
#if DEBUG_RS
LogDebug(LOG_HOST, "reed_solomon_tabs::reed_solomon_tabs() RS(%d,*): initialized for %d symbols size, %d modulo table", SIZE, NN, MODS);
#endif
// Generate Galois field lookup tables
index_of[0] = A0; // log(zero) = -inf
alpha_to[A0] = 0; // alpha**-inf = 0
PLY poly;
int sr = poly(0);
for (int i = 0; i < NN; i++) {
index_of[sr] = i;
alpha_to[i] = sr;
sr = poly(sr);
}
// If it's not primitive, raise exception or abort
if (sr != alpha_to[0])
EZPWD_RAISE_OR_ABORT(std::runtime_error, "reed-solomon: Galois field polynomial not primitive");
// Generate modulo table for some commonly used (non-trivial) values
for (int x = NN; x < NN + MODS; ++x)
mod_of[x - NN] = _modnn(x);
// Find prim-th root of 1, index form, used in decoding.
int iptmp = 1;
while (iptmp % PRM != 0)
iptmp += NN;
iprim = iptmp / PRM;
}
/// <summary>Finalizes a instance of the reed_solomon_tabs class.</summary>
~reed_solomon_tabs() override = default;
//
// modnn -- modulo replacement for galois field arithmetics, optionally w/ table acceleration
//
// @x: the value to reduce (will never be -'ve)
//
// where
// MM = number of bits per symbol
// NN = (2^MM) - 1
//
// Simple arithmetic modulo would return a wrong result for values >= 3 * NN
//
TYP _modnn(int x) const
{
while (x >= NN) {
x -= NN;
x = (x >> MM) + (x & NN);
}
return x;
}
TYP modnn(int x) const
{
while (x >= NN + MODS) {
x -= NN;
x = (x >> MM) + (x & NN);
}
if (MODS && x >= NN)
x = mod_of[x - NN];
return x;
}
};
// ---------------------------------------------------------------------------
// Class Declaration
// ---------------------------------------------------------------------------
/**
* @brief Reed-Solomon codec.
* @tparam TYP A symbol datum; {en,de}code operates on arrays of these
* @tparam SYM Bits per symbol
* @tparam RTS
* @tparam FCR First consecutive root, index form
* @tparam PRM Primitive element, index form
* @tparam PLY The primitive generator polynominal functor
* @ingroup edac_rs
*/
/*
* @TYP: A symbol datum; {en,de}code operates on arrays of these
* @DATUM: Bits per datum (a TYP())
* @SYM{BOL}, MM: Bits per symbol
* @NN: Symbols per block (== (1<<MM)-1)
* @alpha_to: log lookup table
* @index_of: Antilog lookup table
* @genpoly: Generator polynomial
* @NROOTS: Number of generator roots = number of parity symbols
* @FCR: First consecutive root, index form
* @PRM: Primitive element, index form
* @iprim: prim-th root of 1, index form
* @PLY: The primitive generator polynominal functor
*
* All reed_solomon<T, ...> instances with the same template type parameters share a common
* (static) set of alpha_to, index_of and genpoly tables. The first instance to be constructed
* initializes the tables.
*
* Each specialized type of reed_solomon implements a specific encode/decode method
* appropriate to its datum 'TYP'. When accessed via a generic reed_solomon_base pointer, only
* access via "safe" (size specifying) containers or iterators is available.
*/
template<typename TYP, int SYM, int RTS, int FCR, int PRM, class PLY>
class reed_solomon : public reed_solomon_tabs<TYP, SYM, PRM, PLY> {
public:
typedef reed_solomon_tabs<TYP, SYM, PRM, PLY> tabs_t;
using tabs_t::DATUM;
using tabs_t::SYMBOL;
using tabs_t::MM;
using tabs_t::SIZE;
using tabs_t::NN;
using tabs_t::A0;
using tabs_t::iprim;
using tabs_t::alpha_to;
using tabs_t::index_of;
using tabs_t::modnn;
static const int NROOTS = RTS;
static const int LOAD = SIZE - NROOTS; // maximum non-parity symbol payload
protected:
static std::array<TYP, NROOTS + 1> genpoly;
public:
/** @brief Initializes a new instance of the reed_solomon class. */
reed_solomon() : reed_solomon_tabs<TYP, SYM, PRM, PLY>()
{
// We check one element of the array; this is safe, 'cause the value will not be
// initialized 'til the initializing thread has completely initialized the array. Worst
// case scenario: multiple threads will initialize identically. No mutex necessary.
if (genpoly[0])
return;
#if DEBUG_RS
LogDebug(LOG_HOST, "reed_solomon::reed_solomon() RS(%d,%d): initialized for %d roots", SIZE, LOAD, NROOTS);
#endif
std::array<TYP, NROOTS + 1> tmppoly; // uninitialized
// Form RS code generator polynomial from its roots. Only lower-index entries are
// consulted, when computing subsequent entries; only index 0 needs initialization.
tmppoly[0] = 1;
for (int i = 0, root = FCR * PRM; i < NROOTS; i++, root += PRM) {
tmppoly[i + 1] = 1;
// Multiply tmppoly[] by @**(root + x)
for (int j = i; j > 0; j--) {
if (tmppoly[j] != 0)
tmppoly[j] = tmppoly[j - 1] ^ alpha_to[modnn(index_of[tmppoly[j]] + root)];
else
tmppoly[j] = tmppoly[j - 1];
}
// tmppoly[0] can never be zero
tmppoly[0] = alpha_to[modnn(index_of[tmppoly[0]] + root)];
}
// convert NROOTS entries of tmppoly[] to genpoly[] in index form for quicker encoding,
// in reverse order so genpoly[0] is last element initialized.
for (int i = NROOTS; i >= 0; --i)
genpoly[i] = index_of[tmppoly[i]];
}
/** @brief Finalizes a instance of the reed_solomon class. */
virtual ~reed_solomon() = default;
/** @brief A data element's bits. */
virtual size_t datum() const { return DATUM; }
/** @brief A symbol's bits. */
virtual size_t symbol() const { return SYMBOL; }
/** @brief R-S block size (maximum total symbols). */
virtual int size() const { return SIZE; }
/** @brief R-S roots (parity symbols). */
virtual int nroots() const { return NROOTS; }
/** @brief R-S net payload (data symbols). */
virtual int load() const { return LOAD; }
using reed_solomon_base::encode;
virtual int encode(const std::pair<uint8_t*, uint8_t*>& data) const { return encode_mask(data.first, int(data.second - data.first) - NROOTS, data.second - NROOTS); }
virtual int encode(const std::pair<const uint8_t*, const uint8_t*>& data,
const std::pair<uint8_t*, uint8_t*>& parity) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return encode_mask(data.first, int(data.second - data.first), parity.first);
}
virtual int encode(const std::pair<uint16_t*, uint16_t*>& data) const { return encode_mask(data.first, int(data.second - data.first) - NROOTS, data.second - NROOTS); }
virtual int encode(const std::pair<const uint16_t*, const uint16_t*>& data,
const std::pair<uint16_t*, uint16_t*>& parity) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return encode_mask(data.first, int(data.second - data.first), parity.first);
}
virtual int encode(const std::pair<uint32_t*, uint32_t*>& data) const { return encode_mask(data.first, int(data.second - data.first) - NROOTS, data.second - NROOTS); }
virtual int encode(const std::pair<const uint32_t*, const uint32_t*>& data,
const std::pair<uint32_t*, uint32_t*>& parity) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return encode_mask(data.first, int(data.second - data.first), parity.first);
}
template<typename INP>
int encode_mask(const INP* data, int len, INP* parity) const
{
if (len < 1)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: must provide space for all parity and at least one non-parity symbol", -1);
const TYP* dataptr;
TYP* pariptr;
const size_t INPUT = 8 * sizeof(INP);
if (DATUM != SYMBOL || DATUM != INPUT) {
// Our DATUM (TYP) size (eg. uint8_t ==> 8, uint16_t ==> 16, uint32_t ==> 32)
// doesn't exactly match our R-S SYMBOL size (eg. 6), or our INP size; Must mask and
// copy. The INP data must fit at least the SYMBOL size!
if (SYMBOL > INPUT)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: output data type too small to contain symbols", -1);
std::array<TYP, SIZE> tmp;
TYP msk = static_cast<TYP>(~0UL << SYMBOL);
for (int i = 0; i < len; ++i)
tmp[LOAD - len + i] = data[i] & ~msk;
dataptr = &tmp[LOAD - len];
pariptr = &tmp[LOAD];
encode(dataptr, len, pariptr);
// we copied/masked data; copy the parity symbols back (may be different sizes)
for (int i = 0; i < NROOTS; ++i)
parity[i] = pariptr[i];
} else {
// Our R-S SYMBOL size, DATUM size and INP type size exactly matches; use in-place.
dataptr = reinterpret_cast<const TYP*>(data);
pariptr = reinterpret_cast<TYP*>(parity);
encode(dataptr, len, pariptr);
}
return NROOTS;
}
using reed_solomon_base::decode;
virtual int decode(const std::pair<uint8_t*, uint8_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
return decode_mask(data.first, int(data.second - data.first), (uint8_t*)0, erasure, position);
}
virtual int decode(const std::pair<uint8_t*, uint8_t*>& data, const std::pair<uint8_t*, uint8_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return decode_mask(data.first, int(data.second - data.first), parity.first, erasure, position);
}
virtual int decode(const std::pair<uint16_t*, uint16_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
return decode_mask(data.first, int(data.second - data.first), (uint16_t*)0, erasure, position);
}
virtual int decode(const std::pair<uint16_t*, uint16_t*>& data, const std::pair<uint16_t*, uint16_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return decode_mask(data.first, int(data.second - data.first), parity.first, erasure, position);
}
virtual int decode(const std::pair<uint32_t*, uint32_t*>& data,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
return decode_mask(data.first, int(data.second - data.first), (uint32_t*)0, erasure, position);
}
virtual int decode(const std::pair<uint32_t*, uint32_t*>& data, const std::pair<uint32_t*, uint32_t*>& parity,
const std::vector<int>& erasure = std::vector<int>(), std::vector<int>* position = 0) const
{
if (parity.second - parity.first != NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity length incompatible with number of roots", -1);
return decode_mask(data.first, int(data.second - data.first), parity.first, erasure, position);
}
//
// decode_mask -- mask INP data into valid SYMBOL data
//
// Incoming data may be in a variety of sizes, and may contain information beyond the
// R-S symbol capacity. For example, we might use a 6-bit R-S symbol to correct the lower
// 6 bits of an 8-bit data character. This would allow us to correct common substitution
// errors (such as '2' for '3', 'R' for 'T', 'n' for 'm').
//
template<typename INP>
int decode_mask(INP* data, int len, INP* parity = 0, const std::vector<int>& erasure = std::vector<int>(),
std::vector<int>* position = 0) const {
if (len < (parity ? 0 : NROOTS) + 1)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: must provide all parity and at least one non-parity symbol", -1);
if (!parity) {
len -= NROOTS;
parity = data + len;
}
TYP* dataptr;
TYP* pariptr;
const size_t INPUT = 8 * sizeof(INP);
std::array<TYP, SIZE> tmp;
TYP msk = static_cast<TYP>(~0UL << SYMBOL);
const bool cpy = DATUM != SYMBOL || DATUM != INPUT;
if (cpy) {
// Our DATUM (TYP) size (eg. uint8_t ==> 8, uint16_t ==> 16, uint32_t ==> 32)
// doesn't exactly match our R-S SYMBOL size (eg. 6), or our INP size; Must copy.
// The INP data must fit at least the SYMBOL size!
if (SYMBOL > INPUT)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: input data type too small to contain symbols", -1);
for (int i = 0; i < len; ++i)
tmp[LOAD - len + i] = data[i] & ~msk;
dataptr = &tmp[LOAD - len];
for (int i = 0; i < NROOTS; ++i) {
if (TYP(parity[i]) & msk)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: parity data contains information beyond R-S symbol size", -1);
tmp[LOAD + i] = (TYP)parity[i];
}
pariptr = &tmp[LOAD];
} else {
// Our R-S SYMBOL size, DATUM size and INPUT type sizes exactly matches
dataptr = reinterpret_cast<TYP*>(data);
pariptr = reinterpret_cast<TYP*>(parity);
}
int corrects;
if (!erasure.size() && !position) {
// No erasures, and error position info not wanted.
corrects = decode(dataptr, len, pariptr);
} else {
// Either erasure location info specified, or resultant error position info wanted;
// Prepare pos (a temporary, if no position vector provided), and copy any provided
// erasure positions. After number of corrections is known, resize the position
// vector. Thus, we use any supplied erasure info, and optionally return any
// correction position info separately.
std::vector<int> _pos;
std::vector<int>& pos = position ? *position : _pos;
pos.resize(std::max(size_t(NROOTS), erasure.size()));
std::copy(erasure.begin(), erasure.end(), pos.begin());
corrects = decode(dataptr, len, pariptr, &pos.front(), int(erasure.size()));
if (corrects > int(pos.size()))
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: FATAL: produced too many corrections; possible corruption!", -1);
pos.resize(std::max(0, corrects));
}
if (cpy && corrects > 0) {
for (int i = 0; i < len; ++i) {
data[i] &= msk;
data[i] |= tmp[LOAD - len + i];
}
for (int i = 0; i < NROOTS; ++i)
parity[i] = tmp[LOAD + i];
}
return corrects;
}
int encode(const TYP* data, int len, TYP* parity) const
{
// Check length parameter for validity
int pad = NN - NROOTS - len;
if (pad < 0 || pad >= NN)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: data length incompatible with block size and error correction symbols", -1);
for (int i = 0; i < NROOTS; i++)
parity[i] = 0;
for (int i = 0; i < len; i++) {
TYP feedback = index_of[data[i] ^ parity[0]];
if (feedback != A0) {
for (int j = 1; j < NROOTS; j++)
parity[j] ^= alpha_to[modnn(feedback + genpoly[NROOTS - j])];
}
std::rotate(parity, parity + 1, parity + NROOTS);
if (feedback != A0)
parity[NROOTS - 1] = alpha_to[modnn(feedback + genpoly[0])];
else
parity[NROOTS - 1] = 0;
}
return NROOTS;
}
int decode(TYP* data, int len, TYP* parity, int* eras_pos = 0, int no_eras = 0, TYP* corr = 0) const
{
typedef std::array<TYP, NROOTS> typ_nroots;
typedef std::array<TYP, NROOTS + 1> typ_nroots_1;
typedef std::array<int, NROOTS> int_nroots;
typ_nroots_1 lambda{ {0} };
typ_nroots syn;
typ_nroots_1 b;
typ_nroots_1 t;
typ_nroots_1 omega;
int_nroots root;
typ_nroots_1 reg;
int_nroots loc;
int count = 0;
// Check length parameter and erasures for validity
int pad = NN - NROOTS - len;
if (pad < 0 || pad >= NN)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: data length incompatible with block size and error correction symbols", -1);
if (no_eras) {
if (no_eras > NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: number of erasures exceeds capacity (number of roots)", -1);
for (int i = 0; i < no_eras; ++i) {
if (eras_pos[i] < 0 || eras_pos[i] >= len + NROOTS)
EZPWD_RAISE_OR_RETURN(std::runtime_error, "reed-solomon: erasure positions outside data+parity", -1);
}
}
// form the syndromes; i.e., evaluate data(x) at roots of g(x)
for (int i = 0; i < NROOTS; i++)
syn[i] = data[0];
for (int j = 1; j < len; j++) {
for (int i = 0; i < NROOTS; i++) {
if (syn[i] == 0)
syn[i] = data[j];
else
syn[i] = data[j] ^ alpha_to[modnn(index_of[syn[i]] + (FCR + i) * PRM)];
}
}
for (int j = 0; j < NROOTS; j++) {
for (int i = 0; i < NROOTS; i++) {
if (syn[i] == 0)
syn[i] = parity[j];
else
syn[i] = parity[j] ^ alpha_to[modnn(index_of[syn[i]] + (FCR + i) * PRM)];
}
}
// Convert syndromes to index form, checking for nonzero condition
TYP syn_error = 0;
for (int i = 0; i < NROOTS; i++) {
syn_error |= syn[i];
syn[i] = index_of[syn[i]];
}
int deg_lambda = 0;
int deg_omega = 0;
int r = no_eras;
int el = no_eras;
if (!syn_error) {
// if syndrome is zero, data[] is a codeword and there are no errors to correct.
count = 0;
goto finish; // ewww; gotos!
}
lambda[0] = 1;
if (no_eras > 0) {
// Init lambda to be the erasure locator polynomial. Convert erasure positions
// from index into data, to index into Reed-Solomon block.
lambda[1] = alpha_to[modnn(PRM * (NN - 1 - (eras_pos[0] + pad)))];
for (int i = 1; i < no_eras; i++) {
TYP u = modnn(PRM * (NN - 1 - (eras_pos[i] + pad)));
for (int j = i + 1; j > 0; j--) {
TYP tmp = index_of[lambda[j - 1]];
if (tmp != A0)
lambda[j] ^= alpha_to[modnn(u + tmp)];
}
}
}
#if DEBUG_RS
// Test code that verifies the erasure locator polynomial just constructed
// Needed only for decoder debugging.
// find roots of the erasure location polynomial
for (int i = 1; i <= no_eras; i++) {
reg[i] = index_of[lambda[i]];
}
count = 0;
for (int i = 1, k = iprim - 1; i <= NN; i++, k = modnn(k + iprim)) {
TYP q = 1;
for (int j = 1; j <= no_eras; j++) {
if (reg[j] != A0) {
reg[j] = modnn(reg[j] + j);
q ^= alpha_to[reg[j]];
}
}
if (q != 0) {
continue;
}
// store root and error location number indices
root[count] = i;
loc[count] = k;
count++;
}
if (count != no_eras) {
LogDebug(LOG_HOST, "reed_solomon::decode(): count = %d, no_eras = %d, lambda(x) is WRONG", count, no_eras);
count = -1;
goto finish;
}
if (count) {
std::stringstream ss;
ss << "reed_solomon::decode(): Erasure positions as determined by roots of Eras Loc Poly: ";
for (int i = 0; i < count; i++) {
ss << loc[i] << ' ';
}
LogDebug(LOG_HOST, "%s", ss.str().c_str());
ss.clear();
ss << "reed_solomon::decode(): Erasure positions as determined by roots of eras_pos array: ";
for (int i = 0; i < no_eras; i++) {
ss << eras_pos[i] << ' ';
}
LogDebug(LOG_HOST, "%s", ss.str().c_str());
}
#endif
for (int i = 0; i < NROOTS + 1; i++)
b[i] = index_of[lambda[i]];
//
// Begin Berlekamp-Massey algorithm to determine error+erasure locator polynomial
//
while (++r <= NROOTS) {
// r is the step number
// Compute discrepancy at the r-th step in poly-form
TYP discr_r = 0;
for (int i = 0; i < r; i++) {
if ((lambda[i] != 0) && (syn[r - i - 1] != A0))
discr_r ^= alpha_to[modnn(index_of[lambda[i]] + syn[r - i - 1])];
}
discr_r = index_of[discr_r]; // Index form
if (discr_r == A0) {
// 2 lines below: B(x) <-- x*B(x)
// Rotate the last element of b[NROOTS+1] to b[0]
std::rotate(b.begin(), b.begin() + NROOTS, b.end());
b[0] = A0;
} else {
// 7 lines below: T(x) <-- lambda(x)-discr_r*x*b(x)
t[0] = lambda[0];
for (int i = 0; i < NROOTS; i++) {
if (b[i] != A0)
t[i + 1] = lambda[i + 1] ^ alpha_to[modnn(discr_r + b[i])];
else
t[i + 1] = lambda[i + 1];
}
if (2 * el <= r + no_eras - 1) {
el = r + no_eras - el;
// 2 lines below: B(x) <-- inv(discr_r) * lambda(x)
for (int i = 0; i <= NROOTS; i++)
b[i] = ((lambda[i] == 0) ? A0 : modnn(index_of[lambda[i]] - discr_r + NN));
} else {
// 2 lines below: B(x) <-- x*B(x)
std::rotate(b.begin(), b.begin() + NROOTS, b.end());
b[0] = A0;
}
lambda = t;
}
}
// Convert lambda to index form and compute deg(lambda(x))
for (int i = 0; i < NROOTS + 1; i++) {
lambda[i] = index_of[lambda[i]];
if (lambda[i] != NN)
deg_lambda = i;
}
// Find roots of error+erasure locator polynomial by Chien search
reg = lambda;
count = 0; // Number of roots of lambda(x)
for (int i = 1, k = iprim - 1; i <= NN; i++, k = modnn(k + iprim)) {
TYP q = 1; // lambda[0] is always 0
for (int j = deg_lambda; j > 0; j--) {
if (reg[j] != A0) {
reg[j] = modnn(reg[j] + j);
q ^= alpha_to[reg[j]];
}
}
if (q != 0)
continue; // Not a root
// store root (index-form) and error location number
#if DEBUG_RS
LogDebug(LOG_HOST, "reed_solomon::decode(): count = %d, root = %d, loc = %d", count, i, k);
#endif
root[count] = i;
loc[count] = k;
// If we've already found max possible roots, abort the search to save time
if (++count == deg_lambda)
break;
}
if (deg_lambda != count) {
// deg(lambda) unequal to number of roots => uncorrectable error detected
count = -1;
goto finish;
}
//
// Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo x**NROOTS). in
// index form. Also find deg(omega).
//
deg_omega = deg_lambda - 1;
for (int i = 0; i <= deg_omega; i++) {
TYP tmp = 0;
for (int j = i; j >= 0; j--) {
if ((syn[i - j] != A0) && (lambda[j] != A0)) {
tmp ^= alpha_to[modnn(syn[i - j] + lambda[j])];
}
}
omega[i] = index_of[tmp];
}
//
// Compute error values in poly-form. num1 = omega(inv(X(l))), num2 = inv(X(l))**(fcr-1)
// and den = lambda_pr(inv(X(l))) all in poly-form
//
for (int j = count - 1; j >= 0; j--) {
TYP num1 = 0;
for (int i = deg_omega; i >= 0; i--) {
if (omega[i] != A0)
num1 ^= alpha_to[modnn(omega[i] + i * root[j])];
}
TYP num2 = alpha_to[modnn(root[j] * (FCR - 1) + NN)];
TYP den = 0;
// lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i]
for (int i = std::min(deg_lambda, NROOTS - 1) & ~1; i >= 0; i -= 2) {
if (lambda[i + 1] != A0)
den ^= alpha_to[modnn(lambda[i + 1] + i * root[j])];
}
#if DEBUG_RS
if (den == 0) {
LogDebug(LOG_HOST, "reed_solomon::decode(): ERROR: denominator = 0");
count = -1;
goto finish;
}
#endif
// Apply error to data. Padding ('pad' unused symbols) begin at index 0.
if (num1 != 0) {
if (loc[j] < pad) {
// If the computed error position is in the 'pad' (the unused portion of the
// R-S data capacity), then our solution has failed -- we've computed a
// correction location outside of the data and parity we've been provided!
#if DEBUG_RS
std::stringstream ss;
ss << "reed_solomon::decode(): ERROR: RS(" << SIZE << "," << LOAD << ") computed error location: " << loc[j] <<
" within " << pad << " pad symbols, not within " << LOAD - pad << " data or " << NROOTS << " parity";
LogDebug(LOG_HOST, "%s", ss.str().c_str());
#endif
count = -1;
goto finish;
}
TYP cor = alpha_to[modnn(index_of[num1] + index_of[num2] + NN - index_of[den])];
// Store the error correction pattern, if a correction buffer is available
if (corr != nullptr)
corr[j] = cor;
// If a data/parity buffer is given and the error is inside the message or
// parity data, correct it
if (loc[j] < (NN - NROOTS)) {
if (data != nullptr)
data[loc[j] - pad] ^= cor;
} else if (loc[j] < NN) {
if (parity != nullptr)
parity[loc[j] - (NN - NROOTS)] ^= cor;
}
}
}
finish:
#if DEBUG_RS
if (count > NROOTS) {
LogDebug(LOG_HOST, "reed_solomon::decode(): ERROR: number of corrections %d exceeds NROOTS %d", count, NROOTS);
}
if (count > 0) {
std::string errors(2 * (len + NROOTS), '.');
for (int i = 0; i < count; ++i) {
errors[2 * (loc[i] - pad) + 0] = 'E';
errors[2 * (loc[i] - pad) + 1] = 'E';
}
for (int i = 0; i < no_eras; ++i) {
errors[2 * (eras_pos[i]) + 0] = 'e';
errors[2 * (eras_pos[i]) + 1] = 'e';
}
std::stringstream ss;
ss << "reed_solomon::decode(): e)rase, E)rror; count = " << count << ": " << std::endl << errors;
LogDebug(LOG_HOST, "%s", ss.str().c_str());
}
#endif
if (eras_pos != nullptr) {
for (int i = 0; i < count; i++)
eras_pos[i] = loc[i] - pad;
}
return count;
}
};
//
// Define the static reed_solomon...<...> members; allowed in header for template types.
//
// The reed_solomon_tags<...>::iprim < 0 is used to indicate to the first instance that the
// static tables require initialization.
//
template<typename TYP, int SYM, int PRM, class PLY> int reed_solomon_tabs<TYP, SYM, PRM, PLY>::iprim = -1;
template<typename TYP, int SYM, int PRM, class PLY> std::array<TYP, reed_solomon_tabs<TYP, SYM, PRM, PLY>::NN + 1> reed_solomon_tabs<TYP, SYM, PRM, PLY>::alpha_to;
template<typename TYP, int SYM, int PRM, class PLY> std::array<TYP, reed_solomon_tabs<TYP, SYM, PRM, PLY>::NN + 1> reed_solomon_tabs<TYP, SYM, PRM, PLY>::index_of;
template<typename TYP, int SYM, int PRM, class PLY> std::array<TYP, reed_solomon_tabs< TYP, SYM, PRM, PLY>::MODS> reed_solomon_tabs<TYP, SYM, PRM, PLY>::mod_of;
template<typename TYP, int SYM, int RTS, int FCR, int PRM, class PLY> std::array<TYP, reed_solomon< TYP, SYM, RTS, FCR, PRM, PLY>::NROOTS + 1> reed_solomon<TYP, SYM, RTS, FCR, PRM, PLY>::genpoly;
} // namespace rs
#endif
Reference in a new issue View git blame Copy permalink