/* * 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 #include #include #include #include #include #include #include // // 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 #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_ -- compute the log base B of N at compile-time template struct log_ { enum { value = 1 + log_::value }; }; template struct log_<1, B> { enum { value = 0 }; }; template 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 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 cuTpair; typedef std::pair 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 int encode(std::vector& data) const { typedef typename std::make_unsigned::type uT; typedef std::pair uTpair; data.resize(data.size() + nroots()); return encode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size())); } template int encode(const std::vector& data, std::vector& parity) const { typedef typename std::make_unsigned::type uT; typedef std::pair cuTpair; typedef std::pair 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 int encode(std::array& data, int pad = 0) const { typedef typename std::make_unsigned::type uT; typedef std::pair uTpair; return encode(uTpair((uT*)&data.front() + pad, (uT*)&data.front() + data.size())); } virtual int encode(const std::pair& data) const = 0; virtual int encode(const std::pair& data, const std::pair& parity) const = 0; virtual int encode(const std::pair& data) const = 0; virtual int encode(const std::pair& data, const std::pair& parity) const = 0; virtual int encode(const std::pair& data) const = 0; virtual int encode(const std::pair& data, const std::pair& parity) const = 0; int decode(std::string& data, const std::vector& erasure = std::vector(), std::vector* position = 0) const { typedef uint8_t uT; typedef std::pair 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& erasure = std::vector(), std::vector* position = 0) const { typedef uint8_t uT; typedef std::pair uTpair; return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()), uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()), erasure, position); } template int decode(std::vector& data, const std::vector& erasure = std::vector(), std::vector* position = 0) const { typedef typename std::make_unsigned::type uT; typedef std::pair uTpair; return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()), erasure, position); } template int decode(std::vector& data, std::vector& parity, const std::vector& erasure = std::vector(), std::vector* position = 0) const { typedef typename std::make_unsigned::type uT; typedef std::pair uTpair; return decode(uTpair((uT*)&data.front(), (uT*)&data.front() + data.size()), uTpair((uT*)&parity.front(), (uT*)&parity.front() + parity.size()), erasure, position); } template int decode(std::array& data, int pad = 0, const std::vector& erasure = std::vector(), std::vector* position = 0) const { typedef typename std::make_unsigned::type uT; typedef std::pair uTpair; return decode(uTpair((uT*)&data.front() + pad, (uT*)&data.front() + data.size()), erasure, position); } virtual int decode(const std::pair& data, const std::vector& erasure = std::vector(), std::vector* position = 0) const = 0; virtual int decode(const std::pair& data, const std::pair& parity, const std::vector& erasure = std::vector(), std::vector* position = 0) const = 0; virtual int decode(const std::pair& data, const std::vector& erasure = std::vector(), std::vector* position = 0) const = 0; virtual int decode(const std::pair& data, const std::pair& parity, const std::vector& erasure = std::vector(), std::vector* position = 0) const = 0; virtual int decode(const std::pair& data, const std::vector& erasure = std::vector(), std::vector* position = 0) const = 0; virtual int decode(const std::pair& data, const std::pair& parity, const std::vector& erasure = std::vector(), std::vector* 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 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 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 alpha_to; static std::array index_of; static std::array 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; } ///

Finalizes a instance of the reed_solomon_tabs class.

~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< 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 class reed_solomon : public reed_solomon_tabs { public: typedef reed_solomon_tabs 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 genpoly; public: /** @brief Initializes a new instance of the reed_solomon class. */ reed_solomon() : reed_solomon_tabs() { // 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 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& data) const { return encode_mask(data.first, int(data.second - data.first) - NROOTS, data.second - NROOTS); } virtual int encode(const std::pair& data, const std::pair& 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& data) const { return encode_mask(data.first, int(data.second - data.first) - NROOTS, data.second - NROOTS); } virtual int encode(const std::pair& data, const std::pair& 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& data) const { return encode_mask(data.first, int(data.second - data.first) - NROOTS, data.second - NROOTS); } virtual int encode(const std::pair& data, const std::pair& 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 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 tmp; TYP msk = static_cast(~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(data); pariptr = reinterpret_cast(parity); encode(dataptr, len, pariptr); } return NROOTS; } using reed_solomon_base::decode; virtual int decode(const std::pair& data, const std::vector& erasure = std::vector(), std::vector* position = 0) const { return decode_mask(data.first, int(data.second - data.first), (uint8_t*)0, erasure, position); } virtual int decode(const std::pair& data, const std::pair& parity, const std::vector& erasure = std::vector(), std::vector* 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& data, const std::vector& erasure = std::vector(), std::vector* position = 0) const { return decode_mask(data.first, int(data.second - data.first), (uint16_t*)0, erasure, position); } virtual int decode(const std::pair& data, const std::pair& parity, const std::vector& erasure = std::vector(), std::vector* 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& data, const std::vector& erasure = std::vector(), std::vector* position = 0) const { return decode_mask(data.first, int(data.second - data.first), (uint32_t*)0, erasure, position); } virtual int decode(const std::pair& data, const std::pair& parity, const std::vector& erasure = std::vector(), std::vector* 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 int decode_mask(INP* data, int len, INP* parity = 0, const std::vector& erasure = std::vector(), std::vector* 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 tmp; TYP msk = static_cast(~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(data); pariptr = reinterpret_cast(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 _pos; std::vector& 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; typedef std::array typ_nroots_1; typedef std::array 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 int reed_solomon_tabs::iprim = -1; template std::array::NN + 1> reed_solomon_tabs::alpha_to; template std::array::NN + 1> reed_solomon_tabs::index_of; template std::array::MODS> reed_solomon_tabs::mod_of; template std::array::NROOTS + 1> reed_solomon::genpoly; } // namespace rs #endif