From bc43db1e02f26cc7bce9e9058c9ab574329aa61c Mon Sep 17 00:00:00 2001 From: Jonathan Naylor Date: Mon, 30 Sep 2024 18:27:00 +0100 Subject: [PATCH] Update the P25 Reed-Solomon FEC to be the same as the DVMM project. --- MMDVMHost.vcxproj | 5 +- MMDVMHost.vcxproj.filters | 15 +- Makefile | 2 +- Makefile.Pi.Adafruit | 2 +- Makefile.Pi.HD44780 | 2 +- Makefile.Pi.I2C | 2 +- Makefile.Pi.OLED | 2 +- Makefile.Pi.PCF8574 | 2 +- P25Data.cpp | 16 +- P25Data.h | 6 +- RS.h | 1107 ++++++++++++++++++++++++++++++++++ RS241213.cpp => RS634717.cpp | 379 ++++-------- RS241213.h => RS634717.h | 19 +- Version.h | 2 +- 14 files changed, 1249 insertions(+), 312 deletions(-) create mode 100644 RS.h rename RS241213.cpp => RS634717.cpp (62%) rename RS241213.h => RS634717.h (72%) diff --git a/MMDVMHost.vcxproj b/MMDVMHost.vcxproj index 0b29f6354..230e6e5cb 100644 --- a/MMDVMHost.vcxproj +++ b/MMDVMHost.vcxproj @@ -251,8 +251,9 @@ + - + @@ -359,7 +360,7 @@ - + diff --git a/MMDVMHost.vcxproj.filters b/MMDVMHost.vcxproj.filters index 3b74c3121..33b1adb6b 100644 --- a/MMDVMHost.vcxproj.filters +++ b/MMDVMHost.vcxproj.filters @@ -191,9 +191,6 @@ Header Files - - Header Files - Header Files @@ -350,6 +347,12 @@ Header Files + + Header Files + + + Header Files + @@ -511,9 +514,6 @@ Source Files - - Source Files - Source Files @@ -658,5 +658,8 @@ Source Files + + Source Files + \ No newline at end of file diff --git a/Makefile b/Makefile index 495a862e9..9132bc240 100644 --- a/Makefile +++ b/Makefile @@ -16,7 +16,7 @@ OBJECTS = \ Modem.o ModemPort.o ModemSerialPort.o Mutex.o NetworkInfo.o Nextion.o NullController.o NullDisplay.o NXDNAudio.o NXDNControl.o \ NXDNConvolution.o NXDNCRC.o NXDNFACCH1.o NXDNIcomNetwork.o NXDNKenwoodNetwork.o NXDNLayer3.o NXDNLICH.o NXDNLookup.o NXDNNetwork.o NXDNSACCH.o \ NXDNUDCH.o P25Audio.o P25Control.o P25Data.o P25LowSpeedData.o P25Network.o P25NID.o P25Trellis.o P25Utils.o PseudoTTYController.o POCSAGControl.o \ - POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS241213.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o TFTSurenoo.o Thread.o \ + POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS634717.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o TFTSurenoo.o Thread.o \ Timer.o UARTController.o UDPController.o UDPSocket.o UserDB.o UserDBentry.o Utils.o YSFControl.o YSFConvolution.o YSFFICH.o YSFNetwork.o YSFPayload.o all: MMDVMHost RemoteCommand diff --git a/Makefile.Pi.Adafruit b/Makefile.Pi.Adafruit index 42147587b..5442c425f 100644 --- a/Makefile.Pi.Adafruit +++ b/Makefile.Pi.Adafruit @@ -17,7 +17,7 @@ OBJECTS = \ MMDVMHost.o Modem.o ModemPort.o ModemSerialPort.o Mutex.o NetworkInfo.o Nextion.o NullController.o NullDisplay.o NXDNAudio.o \ NXDNControl.o NXDNConvolution.o NXDNCRC.o NXDNFACCH1.o NXDNIcomNetwork.o NXDNKenwoodNetwork.o NXDNLayer3.o NXDNLICH.o NXDNLookup.o NXDNNetwork.o \ NXDNSACCH.o NXDNUDCH.o P25Audio.o P25Control.o P25Data.o P25LowSpeedData.o P25Network.o P25NID.o P25Trellis.o P25Utils.o PseudoTTYController.o \ - POCSAGControl.o POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS241213.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o \ + POCSAGControl.o POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS634717.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o \ TFTSurenoo.o Thread.o Timer.o UARTController.o UDPController.o UDPSocket.o UserDB.o UserDBentry.o Utils.o YSFControl.o YSFConvolution.o YSFFICH.o \ YSFNetwork.o YSFPayload.o diff --git a/Makefile.Pi.HD44780 b/Makefile.Pi.HD44780 index ae98870ee..f857e52bc 100644 --- a/Makefile.Pi.HD44780 +++ b/Makefile.Pi.HD44780 @@ -16,7 +16,7 @@ OBJECTS = \ MMDVMHost.o Modem.o ModemPort.o ModemSerialPort.o Mutex.o NetworkInfo.o Nextion.o NullController.o NullDisplay.o NXDNAudio.o \ NXDNControl.o NXDNConvolution.o NXDNCRC.o NXDNFACCH1.o NXDNIcomNetwork.o NXDNKenwoodNetwork.o NXDNLayer3.o NXDNLICH.o NXDNLookup.o NXDNNetwork.o \ NXDNSACCH.o NXDNUDCH.o P25Audio.o P25Control.o P25Data.o P25LowSpeedData.o P25Network.o P25NID.o P25Trellis.o P25Utils.o PseudoTTYController.o \ - POCSAGControl.o POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS241213.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o \ + POCSAGControl.o POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS634717.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o \ TFTSurenoo.o Thread.o Timer.o UARTController.o UDPController.o UDPSocket.o UserDB.o UserDBentry.o Utils.o YSFControl.o YSFConvolution.o YSFFICH.o \ YSFNetwork.o YSFPayload.o diff --git a/Makefile.Pi.I2C b/Makefile.Pi.I2C index 4c3e2b4c0..43e1d4361 100644 --- a/Makefile.Pi.I2C +++ b/Makefile.Pi.I2C @@ -16,7 +16,7 @@ OBJECTS = \ Modem.o ModemPort.o ModemSerialPort.o Mutex.o NetworkInfo.o Nextion.o NullController.o NullDisplay.o NXDNAudio.o NXDNControl.o \ NXDNConvolution.o NXDNCRC.o NXDNFACCH1.o NXDNIcomNetwork.o NXDNKenwoodNetwork.o NXDNLayer3.o NXDNLICH.o NXDNLookup.o NXDNNetwork.o NXDNSACCH.o \ NXDNUDCH.o P25Audio.o P25Control.o P25Data.o P25LowSpeedData.o P25Network.o P25NID.o P25Trellis.o P25Utils.o PseudoTTYController.o POCSAGControl.o \ - POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS241213.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o TFTSurenoo.o Thread.o \ + POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS634717.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o TFTSurenoo.o Thread.o \ Timer.o UARTController.o UDPController.o UDPSocket.o UserDB.o UserDBentry.o Utils.o YSFControl.o YSFConvolution.o YSFFICH.o YSFNetwork.o YSFPayload.o all: MMDVMHost RemoteCommand diff --git a/Makefile.Pi.OLED b/Makefile.Pi.OLED index 8fdb0503f..ab93d3f86 100644 --- a/Makefile.Pi.OLED +++ b/Makefile.Pi.OLED @@ -20,7 +20,7 @@ OBJECTS = \ Modem.o ModemPort.o ModemSerialPort.o Mutex.o NetworkInfo.o Nextion.o NullController.o NullDisplay.o NXDNAudio.o NXDNControl.o \ NXDNConvolution.o NXDNCRC.o NXDNFACCH1.o NXDNIcomNetwork.o NXDNKenwoodNetwork.o NXDNLayer3.o NXDNLICH.o NXDNLookup.o NXDNNetwork.o NXDNSACCH.o \ NXDNUDCH.o OLED.o P25Audio.o P25Control.o P25Data.o P25LowSpeedData.o P25Network.o P25NID.o P25Trellis.o P25Utils.o PseudoTTYController.o \ - POCSAGControl.o POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS241213.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o \ + POCSAGControl.o POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS634717.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o \ TFTSurenoo.o Thread.o Timer.o UARTController.o UDPController.o UDPSocket.o UserDB.o UserDBentry.o Utils.o YSFControl.o YSFConvolution.o YSFFICH.o \ YSFNetwork.o YSFPayload.o diff --git a/Makefile.Pi.PCF8574 b/Makefile.Pi.PCF8574 index 86a3351a7..6e478ffb2 100644 --- a/Makefile.Pi.PCF8574 +++ b/Makefile.Pi.PCF8574 @@ -17,7 +17,7 @@ OBJECTS = \ MMDVMHost.o Modem.o ModemPort.o ModemSerialPort.o Mutex.o NetworkInfo.o Nextion.o NullController.o NullDisplay.o NXDNAudio.o \ NXDNControl.o NXDNConvolution.o NXDNCRC.o NXDNFACCH1.o NXDNIcomNetwork.o NXDNKenwoodNetwork.o NXDNLayer3.o NXDNLICH.o NXDNLookup.o NXDNNetwork.o \ NXDNSACCH.o NXDNUDCH.o P25Audio.o P25Control.o P25Data.o P25LowSpeedData.o P25Network.o P25NID.o P25Trellis.o P25Utils.o PseudoTTYController.o \ - POCSAGControl.o POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS241213.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o \ + POCSAGControl.o POCSAGNetwork.o QR1676.o RemoteControl.o RS129.o RS634717.o RSSIInterpolator.o SerialPort.o SMeter.o StopWatch.o Sync.o SHA256.o \ TFTSurenoo.o Thread.o Timer.o UARTController.o UDPController.o UDPSocket.o UserDB.o UserDBentry.o Utils.o YSFControl.o YSFConvolution.o YSFFICH.o \ YSFNetwork.o YSFPayload.o diff --git a/P25Data.cpp b/P25Data.cpp index 4a8aa58a0..21e71cc66 100644 --- a/P25Data.cpp +++ b/P25Data.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2016,2017,2023 by Jonathan Naylor G4KLX +* Copyright (C) 2016,2017,2023,2024 by Jonathan Naylor G4KLX * Copyright (C) 2018 by Bryan Biedenkapp N2PLL * * This program is free software; you can redistribute it and/or modify @@ -44,7 +44,7 @@ m_lcf(0x00U), m_emergency(false), m_srcId(0U), m_dstId(0U), -m_rs241213(), +m_rs(), m_trellis() { m_mi = new unsigned char[P25_MI_LENGTH_BYTES]; @@ -89,7 +89,7 @@ bool CP25Data::decodeHeader(const unsigned char* data) // decode RS (36,20,17) FEC try { - bool ret = m_rs241213.decode362017(rs); + bool ret = m_rs.decode362017(rs); if (!ret) return false; } catch (...) { @@ -133,7 +133,7 @@ void CP25Data::encodeHeader(unsigned char* data) rs[14U] = (m_dstId >> 0) & 0xFFU; // Talkgroup Address LSB // encode RS (36,20,17) FEC - m_rs241213.encode362017(rs); + m_rs.encode362017(rs); unsigned char raw[81U]; ::memset(raw, 0x00U, 81U); @@ -171,7 +171,7 @@ bool CP25Data::decodeLDU1(const unsigned char* data) decodeLDUHamming(raw, rs + 15U); try { - bool ret = m_rs241213.decode(rs); + bool ret = m_rs.decode241213(rs); if (!ret) return false; } catch (...) { @@ -234,7 +234,7 @@ void CP25Data::encodeLDU1(unsigned char* data) break; } - m_rs241213.encode(rs); + m_rs.encode241213(rs); unsigned char raw[5U]; encodeLDUHamming(raw, rs + 0U); @@ -284,7 +284,7 @@ bool CP25Data::decodeLDU2(const unsigned char* data) // decode RS (24,16,9) FEC try { - bool ret = m_rs241213.decode24169(rs); + bool ret = m_rs.decode24169(rs); if (!ret) return false; } catch (...) { @@ -324,7 +324,7 @@ void CP25Data::encodeLDU2(unsigned char* data) rs[11U] = (m_kId >> 0) & 0xFFU; // Key ID LSB // encode RS (24,16,9) FEC - m_rs241213.encode24169(rs); + m_rs.encode24169(rs); // encode Hamming (10,6,3) FEC and interleave for LC data unsigned char raw[5U]; diff --git a/P25Data.h b/P25Data.h index 026c5c4a9..ff5e8073f 100644 --- a/P25Data.h +++ b/P25Data.h @@ -1,5 +1,5 @@ /* -* Copyright (C) 2016,2017 by Jonathan Naylor G4KLX +* Copyright (C) 2016,2017,2024 by Jonathan Naylor G4KLX * Copyright (C) 2018 by Bryan Biedenkapp N2PLL * * This program is free software; you can redistribute it and/or modify @@ -20,7 +20,7 @@ #if !defined(P25Data_H) #define P25Data_H -#include "RS241213.h" +#include "RS634717.h" #include "P25Trellis.h" class CP25Data { @@ -81,7 +81,7 @@ class CP25Data { unsigned int m_srcId; unsigned int m_dstId; unsigned char m_serviceType; - CRS241213 m_rs241213; + CRS634717 m_rs; CP25Trellis m_trellis; void decodeLDUHamming(const unsigned char* raw, unsigned char* data); diff --git a/RS.h b/RS.h new file mode 100644 index 000000000..16cc53f28 --- /dev/null +++ b/RS.h @@ -0,0 +1,1107 @@ +/* + * Copyright (C) 2014 Hard Consulting Corporation + * Copyright (C) 2023 Bryan Biedenkapp, N2PLL + * Copyright (C) 2024 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 != NULL) + 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 != NULL) + data[loc[j] - pad] ^= cor; + } else if (loc[j] < NN) { + if (parity != NULL) + 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 != NULL) { + 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 diff --git a/RS241213.cpp b/RS634717.cpp similarity index 62% rename from RS241213.cpp rename to RS634717.cpp index 246da0047..372598d69 100644 --- a/RS241213.cpp +++ b/RS634717.cpp @@ -1,6 +1,6 @@ /* -* Copyright (C) 2016 by Jonathan Naylor G4KLX -* Copyright (C) 2018 by Bryan Biedenkapp N2PLL +* Copyright (C) 2016,2024 by Jonathan Naylor G4KLX +* Copyright (C) 2018,2023 by Bryan Biedenkapp N2PLL * * 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 @@ -17,13 +17,16 @@ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ -#include "RS241213.h" +#include "RS634717.h" +#include "RS.h" + +#include #include #include #include -const unsigned char ENCODE_MATRIX[12U][24U] = { +const unsigned char ENCODE_MATRIX_241213[12U][24U] = { {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 062, 044, 003, 025, 014, 016, 027, 003, 053, 004, 036, 047}, {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 011, 012, 011, 011, 016, 064, 067, 055, 001, 076, 026, 073}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 003, 001, 005, 075, 014, 006, 020, 044, 066, 006, 070, 066}, @@ -118,20 +121,69 @@ static void hex2Bin(unsigned char input, unsigned char* output, unsigned int off WRITE_BIT(output, offset + 5U, input & 0x01U); } -CRS241213::CRS241213() +#define __RS(TYPE, SYMBOLS, PAYLOAD, POLY, FCR, PRIM) \ + rs::reed_solomon::value, \ + (SYMBOLS) - (PAYLOAD), FCR, PRIM, \ + rs::gfpoly::value, POLY>> + +#define __RS_63(PAYLOAD) __RS(unsigned char, 63, PAYLOAD, 0x43, 1, 1) + +class RS6347 : public __RS_63(47) { +public: + RS6347() : __RS_63(47)() { /* stub */ } +}; +RS6347 rs634717; // 16 bit / 8 bit corrections max / 5 bytes total + +/** + * @brief Implements Reed-Solomon (24,12,13) + */ +class RS6351 : public __RS_63(51) { +public: + RS6351() : __RS_63(51)() { /* stub */ } +}; +RS6351 rs241213; // 12 bit / 6 bit corrections max / 3 bytes total + +/** + * @brief Implements Reed-Solomon (24,16,9) + */ +class RS6355 : public __RS_63(55) { +public: + RS6355() : __RS_63(55)() { /* stub */ } +}; +RS6355 rs24169; // 8 bit / 4 bit corrections max / 2 bytes total + +CRS634717::CRS634717() { } -CRS241213::~CRS241213() +CRS634717::~CRS634717() { } -bool CRS241213::decode(unsigned char* data) +bool CRS634717::decode241213(unsigned char* data) { - return decode(data, 24U, 39, 12); + assert(data != NULL); + + std::vector codeword(63, 0); + + unsigned int offset = 0U; + for (unsigned int i = 0U; i < 24U; i++, offset += 6U) + codeword[39U + i] = bin2Hex(data, offset); + + int ec = rs241213.decode(codeword); + + offset = 0U; + for (unsigned int i = 0U; i < 12U; i++, offset += 6U) + hex2Bin(codeword[39U + i], data, offset); + + if (ec == -1 || ec >= 6) + return false; + + return true; } -void CRS241213::encode(unsigned char* data) +void CRS634717::encode241213(unsigned char* data) { assert(data != NULL); @@ -143,7 +195,7 @@ void CRS241213::encode(unsigned char* data) unsigned int offset = 0U; for (unsigned int j = 0U; j < 12U; j++, offset += 6U) { unsigned char hexbit = bin2Hex(data, offset); - codeword[i] ^= gf6Mult(hexbit, ENCODE_MATRIX[j][i]); + codeword[i] ^= gf6Mult(hexbit, ENCODE_MATRIX_241213[j][i]); } } @@ -152,12 +204,29 @@ void CRS241213::encode(unsigned char* data) hex2Bin(codeword[i], data, offset); } -bool CRS241213::decode24169(unsigned char* data) +bool CRS634717::decode24169(unsigned char* data) { - return decode(data, 24U, 39, 8); + assert(data != NULL); + + std::vector codeword(63, 0); + + unsigned int offset = 0U; + for (unsigned int i = 0U; i < 24U; i++, offset += 6U) + codeword[39U + i] = bin2Hex(data, offset); + + int ec = rs24169.decode(codeword); + + offset = 0U; + for (unsigned int i = 0U; i < 16U; i++, offset += 6U) + hex2Bin(codeword[39U + i], data, offset); + + if (ec == -1 || ec >= 4) + return false; + + return true; } -void CRS241213::encode24169(unsigned char* data) +void CRS634717::encode24169(unsigned char* data) { assert(data != NULL); @@ -178,12 +247,29 @@ void CRS241213::encode24169(unsigned char* data) hex2Bin(codeword[i], data, offset); } -bool CRS241213::decode362017(unsigned char* data) +bool CRS634717::decode362017(unsigned char* data) { - return decode(data, 36U, 27, 16); + assert(data != NULL); + + std::vector codeword(63, 0); + + unsigned int offset = 0U; + for (unsigned int i = 0U; i < 36U; i++, offset += 6U) + codeword[27U + i] = bin2Hex(data, offset); + + int ec = rs634717.decode(codeword); + + offset = 0U; + for (unsigned int i = 0U; i < 20U; i++, offset += 6U) + hex2Bin(codeword[27U + i], data, offset); + + if (ec == -1 || ec >= 8) + return false; + + return true; } -void CRS241213::encode362017(unsigned char* data) +void CRS634717::encode362017(unsigned char* data) { assert(data != NULL); @@ -205,7 +291,7 @@ void CRS241213::encode362017(unsigned char* data) } // GF(2 ^ 6) multiply(for Reed - Solomon encoder) -unsigned char CRS241213::gf6Mult(unsigned char a, unsigned char b) const +unsigned char CRS634717::gf6Mult(unsigned char a, unsigned char b) const { unsigned char p = 0x00U; @@ -223,262 +309,3 @@ unsigned char CRS241213::gf6Mult(unsigned char a, unsigned char b) const return p; } - -bool CRS241213::decode(unsigned char* data, const unsigned int bitLength, const int firstData, const int roots) -{ - assert(data != NULL); - - //unsigned char HB[24U]; - unsigned char HB[63U]; - ::memset(HB, 0x00U, 63U); - - unsigned int offset = 0U; - for (unsigned int i = 0U; i < bitLength; i++, offset += 6) - HB[i] = bin2Hex(data, offset); - - //RS (63,63-nroots,nroots+1) decoder where nroots = number of parity bits - // rsDec(8, 39) rsDec(16, 27) rsDec(12, 39) - - const int nroots = roots; - int lambda[18]; // Err+Eras Locator poly - int S[17]; // syndrome poly - int b[18]; - int t[18]; - int omega[18]; - int root[17]; - int reg[18]; - int locn[17]; - - int i, j, count, r, el, SynError, DiscrR, q, DegOmega, tmp, num1, num2, den, DegLambda; - - // form the syndromes; i.e., evaluate HB(x) at roots of g(x) - for (i = 0; i <= nroots - 1; i++) { - S[i] = HB[0]; - } - - //for (j = 1; j <= 24; j++) { // XXX was 62 - //for (j = 1; j <= (int)(bitLength - 1); j++) { - for (j = 1; j <= 62; j++) { - for (i = 0; i <= nroots - 1; i++) { - if (S[i] == 0) { - S[i] = HB[j]; - } - else { - S[i] = HB[j] ^ rsGFexp[(rsGFlog[S[i]] + i + 1) % 63]; - } - } - } - - // convert syndromes to index form, checking for nonzero condition - SynError = 0; - - for (i = 0; i <= nroots - 1; i++) { - SynError = SynError | S[i]; - S[i] = rsGFlog[S[i]]; - } - - if (SynError == 0) { - // if syndrome is zero, rsData[] is a codeword and there are - // no errors to correct. So return rsData[] unmodified - count = 0; - return true; - } - - for (i = 1; i <= nroots; i++) { - lambda[i] = 0; - } - - lambda[0] = 1; - - for (i = 0; i <= nroots; i++) { - b[i] = rsGFlog[lambda[i]]; - } - - // begin Berlekamp-Massey algorithm to determine error+erasure - // locator polynomial - r = 0; - el = 0; - while (++r <= nroots) { - // r is the step number - //r = r + 1; - // compute discrepancy at the r-th step in poly-form - DiscrR = 0; - - for (i = 0; i <= r - 1; i++) { - if ((lambda[i] != 0) && (S[r - i - 1] != 63)) { - DiscrR = DiscrR ^ rsGFexp[(rsGFlog[lambda[i]] + S[r - i - 1]) % 63]; - } - } - - DiscrR = rsGFlog[DiscrR]; // index form - - if (DiscrR == 63) { - // shift elements upward one step - for (i = nroots; i >= 1; i += -1) { - b[i] = b[i - 1]; - } - - b[0] = 63; - } - else { - // t(x) <-- lambda(x) - DiscrR*x*b(x) - t[0] = lambda[0]; - - for (i = 0; i <= nroots - 1; i++) { - if (b[i] != 63) { - t[i + 1] = lambda[i + 1] ^ rsGFexp[(DiscrR + b[i]) % 63]; - } - else { - t[i + 1] = lambda[i + 1]; - } - } - - if (2 * el <= r - 1) { - el = r - el; - // b(x) <-- inv(DiscrR) * lambda(x) - - for (i = 0; i <= nroots; i++) { - if (lambda[i]) { - b[i] = (rsGFlog[lambda[i]] - DiscrR + 63) % 63; - } - else { - b[i] = 63; - } - } - } - else { - // shift elements upward one step - for (i = nroots; i >= 1; i += -1) { - b[i] = b[i - 1]; - } - - b[0] = 63; - } - - for (i = 0; i <= nroots; i++) { - lambda[i] = t[i]; - } - } - } /* end while() */ - - // convert lambda to index form and compute deg(lambda(x)) - DegLambda = 0; - for (i = 0; i <= nroots; i++) { - lambda[i] = rsGFlog[lambda[i]]; - - if (lambda[i] != 63) { - DegLambda = i; - } - } - - // Find roots of the error+erasure locator polynomial by Chien search - for (i = 1; i <= nroots; i++) { - reg[i] = lambda[i]; - } - - count = 0;// number of roots of lambda(x) - - for (i = 1; i <= 63; i++) { - q = 1;// lambda[0] is always 0 - - for (j = DegLambda; j >= 1; j += -1) { - if (reg[j] != 63) { - reg[j] = (reg[j] + j) % 63; - q = q ^ rsGFexp[reg[j]]; - } - } - - // it is a root - if (q == 0) { - // store root (index-form) and error location number - root[count] = i; - locn[count] = i - 40; - // if we have max possible roots, abort search to save time - count = count + 1; - - if (count == DegLambda) { - break; - } - } - } - - if (DegLambda != count) { - // deg(lambda) unequal to number of roots => uncorrectable error detected - return false; - } - - // compute err+eras evaluator poly omega(x) - // = s(x)*lambda(x) (modulo x**nroots). in index form. Also find deg(omega). - DegOmega = 0; - for (i = 0; i <= nroots - 1; i++) { - tmp = 0; - if (DegLambda < i) { - j = DegLambda; - } - else { - j = i; - } - - for ( /* j = j */; j >= 0; j += -1) { - if ((S[i - j] != 63) && (lambda[j] != 63)) { - tmp = tmp ^ rsGFexp[(S[i - j] + lambda[j]) % 63]; - } - } - - if (tmp) { - DegOmega = i; - } - - omega[i] = rsGFlog[tmp]; - } - - omega[nroots] = 63; - - //compute error values in poly-form: - // num1 = omega(inv(X(l))) - // num2 = inv(X(l))**(FCR - 1) - // den = lambda_pr(inv(X(l))) - for (j = count - 1; j >= 0; j += -1) { - num1 = 0; - - for (i = DegOmega; i >= 0; i += -1) { - if (omega[i] != 63) { - num1 = num1 ^ rsGFexp[(omega[i] + i * root[j]) % 63]; - } - } - - num2 = rsGFexp[0]; - den = 0; - - // lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] - if (DegLambda < nroots) { - i = DegLambda; - } - else { - i = nroots; - } - - for (i = i & ~1; i >= 0; i += -2) { - if (lambda[i + 1] != 63) { - den = den ^ rsGFexp[(lambda[i + 1] + i * root[j]) % 63]; - } - } - - if (den == 0) { - return false; - } - - // apply error to data - if (num1 != 0) { - if (locn[j] < firstData) - return false; - HB[locn[j]] = HB[locn[j]] ^ (rsGFexp[(rsGFlog[num1] + rsGFlog[num2] + 63 - rsGFlog[den]) % 63]); - } - } - - offset = 0U; - for (unsigned int i = 0U; i < (unsigned int)nroots; i++, offset += 6) - hex2Bin(HB[i], data, offset); - - return true; -} diff --git a/RS241213.h b/RS634717.h similarity index 72% rename from RS241213.h rename to RS634717.h index dcbeeed99..de1ebd7bd 100644 --- a/RS241213.h +++ b/RS634717.h @@ -1,6 +1,6 @@ /* -* Copyright (C) 2016 by Jonathan Naylor G4KLX -* Copyright (C) 2018 by Bryan Biedenkapp N2PLL +* Copyright (C) 2016,2024 by Jonathan Naylor G4KLX +* Copyright (C) 2018,2023 by Bryan Biedenkapp N2PLL * * 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 @@ -17,26 +17,25 @@ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ -#if !defined(RS241213_H) -#define RS241213_H +#if !defined(RS634717_H) +#define RS634717_H -class CRS241213 +class CRS634717 { public: - CRS241213(); - ~CRS241213(); + CRS634717(); + ~CRS634717(); - bool decode(unsigned char* data); + bool decode241213(unsigned char* data); bool decode24169(unsigned char* data); bool decode362017(unsigned char* data); - void encode(unsigned char* data); + void encode241213(unsigned char* data); void encode24169(unsigned char* data); void encode362017(unsigned char* data); private: unsigned char gf6Mult(unsigned char a, unsigned char b) const; - bool decode(unsigned char* data, const unsigned int bitLength, const int firstData, const int roots); }; #endif diff --git a/Version.h b/Version.h index ed23b164b..fcb4eb522 100644 --- a/Version.h +++ b/Version.h @@ -19,6 +19,6 @@ #if !defined(VERSION_H) #define VERSION_H -const char* VERSION = "20240923"; +const char* VERSION = "20240930"; #endif