diff --git a/EDAC/RS/GaloisField.cs b/EDAC/RS/GaloisField.cs new file mode 100644 index 0000000..4f4a255 --- /dev/null +++ b/EDAC/RS/GaloisField.cs @@ -0,0 +1,624 @@ +/** +* Digital Voice Modem - Fixed Network Equipment +* AGPLv3 Open Source. Use is subject to license terms. +* DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. +* +* @package DVM / Fixed Network Equipment +* +*/ +// +// Based on code from the ErrorCorrection project. (https://github.com/antiduh/ErrorCorrection) +// Licensed under the BSD 2-clause License (https://opensource.org/license/bsd-2-clause/) +// +/* +* Copyright (C) 2016 by Kevin Thompson +* +* Redistribution and use in source and binary forms, with or without +* modification, are permitted provided that the following conditions are met: +* +* 1. Redistributions of source code must retain the above copyright notice, this +* list of conditions and the following disclaimer. +* 2. Redistributions in binary form must reproduce the above copyright notice, +* this list of conditions and the following disclaimer in the documentation +* and/or other materials provided with the distribution. +* +* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR +* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ +using System; +using System.Text; + +namespace fnecore.EDAC.RS +{ + /// + /// Implements polynomial math for polynomials over a galois field with characteristic 2 (a binary + /// galois field), with restrictions on the practical size of the field. + /// + /// + /// This paper has a Reed Solomon-relevant intro to GFs: + /// http://downloads.bbc.co.uk/rd/pubs/whp/whp-pdf-files/WHP031.pdf + /// + /// A Galois field is an arithmetical field that is: + /// - Finite: Arithmetic operations are defined for only the elements defined to be in the field, + /// and there are a finite number of elements. + /// - Closed: Arithmetic operations on elements of the field return a value also in the field. + /// + /// A Galois field's size is a parameter of the field, but only certain sizes can exist. The + /// field must be a size = p^k where p is prime and k is a positive integer. p is also called the + /// characteristic of the field. + /// + /// A Galois field is composed of an additive group and a multiplicative group - what it means to + /// perform addition or multiplication in the field depends on the definitions of these groups for + /// the field in question - the choice of what meaning to use is up to the designer of the field. + /// However, both groups must meet the requirements for operations in a Galois field - they must be + /// finite and closed. + /// + /// One such possible definition for addition over GF(p) (where p is prime) could be `z = (x + y) mod p`. + /// Another might be `z = a XOR b` where XOR represents the bitwise exclusive-or operation. Both of these + /// definitions are sufficient for defining the meaning of addition in a Galois field. + /// + /// The same is true for defining multiplication - in GF(p), multiplication could be defined as + /// `z = (x * y ) mod p` - such a definition is finite and closed. + /// + /// A Galois field need not be defined simply over a prime number of elements - it is possible + /// to define a Galois field GF(p^n) (where p is prime and n is > 1), however, providing definitions for + /// addition and multiplication prove a bit tricky - particularly, multiplication is the trouble maker. + /// Our go-to definition for multiplication over such a field might be `z = (x * y) mod p^k`, however + /// such a definition is not sufficient - the ring of integers modulo p^k doesn't meet the definition of + /// a field. + /// + /// Instead, however, we can define the multiplicative group for our finite field in a different manner. + /// What if we treated the elements of GF(p^n) as polynomials over the field GF(p)? Each element in GF(p^n) + /// would be decomposed into a set of elements in GF(p), and those elements would be used as coefficients + /// for a polynomial over GF(p). + /// + /// Consider a galois field that uses characteristic '2', such as GF(2^4). Such a field has 16 elements, and + /// we could consider each element as a set of elements from GF(2) - we could consider them by their bits. + /// + /// For instance, the element '14' from GF(2^4) could be decomposed into binary '1110'. Now we could + /// represent the element '14' from GF(2^4) as the polynomial '1x^3 + 1x^2 + 1x^1 + 0x^0' from GF(2). + /// + /// Let's keep that in our pocket for now. + /// + /// In Galois Fields, non-zero elements of the field form a multiplicative group that is also a cyclic + /// group. For certain elements of the field, successive powers of these elements form a cyclic group + /// that contains the entire set of non-zero elements in the field. These elements are called primitive + /// elements, usually denoted as some value 'a'. Since successive powers of these elements forms + /// a complete cyclic group containing all non-zero elements of the field - all elements are + /// represented, and represented only once - we can represent the elements of the field indexed + /// by the power of some primitive element a, since we can guarantee that there is a one-to-one + /// mapping from some power a^i to some element in the field. The field can be represented in + /// an alternate representation { 0, a^0, a^1, ..., a^(n-2) }, which is particularly useful when + /// defining multiplication in some GF(p^n) as polynomials over GF(p). + /// + /// For example, consider GF(5) - this is a prime group, so we can choose to define arithmetic as + /// simply being addition modulo p and multiplication modulo p. + /// + /// The elements of GF(5) are {0, 1, 2, 3, 4}. 2 + 4 = 6 mod 5 = 1. 2 * 4 = 8 mod 5 = 3. + /// + /// Lets test to see if the element 2 in GF(5) is also a primitive element: + /// - 2^1 = 2 + /// - 2^2 = 4 + /// - 2^3 = 8 mod 5 = 3 + /// - 2^4 = 16 mod 5 = 1 + /// + /// Each element was generated, and generated only once, thus the element 2 in GF(5) is primitive in GF(5). + /// We can now write the elements of GF(5) indexed by a=2: + /// {0, a^0, a^1, a^2, a^3} = + /// {0 1 2, 4, 3} + /// + /// Consider the field GF(7) equipped with addition/multiplication modulo 7, and see if 2 is a primitive element: + /// - 2^0 = 1 + /// - 2^1 = 2 + /// - 2^2 = 4 + /// - 2^3 = 8 mod 7 = 1 xxx + /// + /// 2 in GF(7) forms a cyclic group {2, 4, 1}, but isn't a complete cyclic group, + /// and thus is not a primitive element in GF(7). + /// + /// -- + /// + /// Recall that we don't currently have a definition for how to perform arithmetic on elements + /// from some arbitrary field GF(p^k), since our go-to operator modulus doesn't fit the bill - + /// the ring of integers modulo p^k doesn't meet the definitions for a field. + /// + /// However, we *can* use modulus to define arthimetic on Galois Fields of the form GF(p). Then, we + /// can use polynomials over GF(p), which we are now equiped to compute, as the basis for defining + /// arithmetic over some arbitrary GF(p^k). + /// + /// A polynomial over some arbitrary GF(x) has polynomial coefficients that are elements of GF(x). + /// Consider the following equation of polynomials over the GF(5) that is using modulus to define + /// arithmetic: + /// ( 3x^2 + x + 2 ) + ( 3x^2 + x + 3) = ( 6x^2 + 2x + 5 ) = ( x^2 + 2x + 0 ) + /// + /// We can represent *elements* of GF(p^k) as *polynomials* over GF(p). + /// For instance, element '7' in GF(2^4) would be the following polynomial over GF(2): + /// 0x^3 + 1x^2 + 1x + 1. + /// + /// And element '10' would be the polynomial: + /// 1x^3 + 0x^2 + 1x + 0 + /// + /// We can then add those two polynomials together according to our definition of arithmetic in GF(2): + /// 1x^3 + 1x^2 + 2x + 1 which reduces modulo 2 to + /// 1x^3 + 1x^2 + 0x + 1 + /// + /// Which, when converted back to an element in GF(2^4) would be 13. + /// Thus, we've defined elemental addition in GF(2^4) to be polynomial addition in + /// the system of GF(2) equipped with elemental addition modulu 2. + /// + /// We can similarly use GF(2) with modulu to define multiplication in GF(2^4), however, we + /// run into one little snag - we end up with polynomials with higher-power terms than defined + /// in GF(2^4): + /// + /// Multiply 7 * 6 in GF(2^4): + /// (x^2 + x + 1) * (x^2 + x) = + /// (x^4 + x^3) + (x^3 + x^2) + (x^2 + x) = + /// x^4 + 2x^3 + 2x^2 + x + 0 = + /// x^4 + 0x^3 + 0x^2 + x + 0 = + /// x^4 + x + /// + /// We're left with this x^4 term that has no representation in GF(2^4). + /// + /// However, if we can find a irreducible polynomial in GF(p) that has degree k, then we can + /// perform polynomial division and use the remainder as the result, effectively, implmenting + /// polynomial modulus. The irreducible polynomial must have a non-zero term of degree k (it must be monic) + /// for it to have the reducing strength required to remove polynomial terms that are too high order. + /// + /// For GF(2^4), such a polynomial might be x^4 + x + 1 + /// + /// Thus, we could take our earlier result x^4 + x and reduce it x^4 + x + 1 via polynomial division: + /// 1 + /// __________ + /// x^4 + x + 1 | x^4 + x + 0 + /// - x^4 + x + 1 + /// 0 + 0 + 1 + /// + /// And thus, we're left with just 1. This was a polynomial in GF(2), which + /// we can convert back to an element in GF(2^4), which is simply element 1. + /// + /// Thus, for the field GF(2^4), with reducing polynomial x^4 + x + 1, + /// 7 * 6 = 1. + /// + /// We call this reducing polynomial the field generator polynomial p(x). + /// + /// -- + /// + /// The primitive elements of the field GF(p^k) also form roots of all generator polynomials. + /// For instance, the generator polynomial we chose above: + /// p(x) = x^4 + x + 1 --> + /// p(a) = a^4 + a + 1 = 0 ---> + /// a^4 = a + 1 + /// + /// We have defined addition and multiplication for non-prime field of the form GF(p^k) + /// by using polynomal representation in GF(p) equipped with modular arithmetic. + /// + /// We can use this fact to generate the field of GF(p^k) in terms of the elements { 0, a^0, .., a^(n-2) } + /// + /// Suppose we pick 2 as our primitive element for GF(2^4), and use reducing polynomial x^4 + x + 1, and thus + /// know that a^4 = a + 1 + /// + /// + /// a^0 = = 1 + /// a^1 = = 2 + /// a^2 = = 4 + /// a^3 = = 8 + /// a^4 = a + 1 = 3 + /// a^5 = a(a^4 + 1) = a(a+1) = a^2 + a = 6 + /// ... + /// + /// And thus, we can represent the entire field indexed by the power of our primitive element, and implement + /// multiplication and division using logarithms. + /// + /// -- + /// + /// This design uses the following optimizating restrictions: + /// - The field characteristic must be 2; that is the field must be of the form GF(2^m). + /// - The generator polynomial is specified in the form of a 32-bit 'int' variable, + /// limiting the length of the polynomial to 32 coefficients. This limit will never be + /// practically hit, due to the next restriction. + /// - Multiplication and division are implemented using a look-up table, for performance. + /// Large galois fields require n^2 memory to store the table. For GF(2^8), this is + /// 256 * 256 == 65536 entries, times 4 bytes per entry = 256 kiB. + /// GF(2^16) would require 16 GiB of memory. + /// + public sealed class GaloisField + { + // Note: Most examples in the comments for this class are made in GF(2^4) with p(x) = x^4 + x + 1 + // The field looks like: + // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 + // { 0, 1, 2, 4, 8, 3, 6, 12, 11, 5, 10, 7, 14, 15 13, 9} + // + // Eg: + // 0 = field[0] = 0 + // a^0 = field[1] = 1; + // a^1 = field[2] = 2; + // a^2 = feild[3] = 4; + // a^3 = field[4] = 8; + // ... + // a^7 = field[8] = 11. + + /// + /// The elements of the field, ordered as generated by the field generator polynomial. + /// + public readonly int[] Field; + + /// + /// Stores the multiplicative inverses of the elements of the field, eg + /// the value of 1/a^9 is stored + /// + public readonly int[] Inverses; + + /// + /// Stores the values of the field indexed by their multiplicative logarithm. + /// + public readonly int[] Logarithms; + + /// + /// The primitive polynomial used to generate the elements of the field by taking successive powers of + /// the polynomial. + /// + private int fieldGenPoly; + + /// + /// Caches multiplication values for field elements. + /// + private int[,] multTable; + + /// + /// Stores the total number of elements in the field. As a consequence of the structure of + /// Galois Feilds, this value must be of the form p^k where p is a prime number and k is a + /// positive integer. + /// + private int size; + + /* + ** Methods + */ + + /// + /// Initializes a new instance of the class. + /// + /// The total number of elements in the field. Must be a power of two. + /// A primitive element of the field, represented in polynomial form, + /// that is used to generate the elements of the field in an ordered manner. The bits of the value + /// indicate the coefficients of the polynomial, eg, 0b1011 indicate 1*x^3 + 0*x^2 + 1*x^1 + 1*x^0. + public GaloisField(int size, int fieldGenPoly) + { + this.size = size; + this.fieldGenPoly = fieldGenPoly; + + this.Field = new int[size]; + this.Logarithms = new int[this.Field.Length]; + this.Inverses = new int[this.Field.Length]; + + BuildField(); + BuildLogarithms(); + BuildMultTable(); + BuildInverses(); + } + + /// + /// Pretty-prints a polynomial with the given coefficients. + /// + /// + /// + public static string PolyPrint(int[] poly) + { + StringBuilder builder = new StringBuilder(poly.Length * 3); + for (int i = poly.Length - 1; i >= 0; i--) + { + if (i > 1) + builder.Append(poly[i]).Append("x^").Append(i).Append(" + "); + else if (i == 1) + builder.Append(poly[i]).Append("x").Append(" + "); + else + builder.Append(poly[i]); + } + + return builder.ToString(); + } + + /// + /// Performs division, as defined by by this Galois field implementation, on the + /// given operands. + /// + /// The value to act as the dividend. + /// The value to act as the divisor. + /// + public int Divide(int dividend, int divisor) + { + // Using the original computation is the same speed as the multiplication table. + // I don't know why. + return this.multTable[dividend, Inverses[divisor]]; + } + + /// + /// Performs the multiplication operation, as defined by this Galois field implemention, + /// on the given operands. + /// + /// The first value to multiply. + /// The second value to multiply. + /// + public int Multiply(int left, int right) + { + // Using the multiplication table is a lot faster than the original computation. + return this.multTable[left, right]; + } + + /// + /// Evaluates the polynomial for the given value of x + /// + /// The polynomial to evaluate + /// The value with which to evaluate the polynomail. + /// + public int PolyEval(int[] poly, int x) + { + int sum; + int xLog; + int coeffLog; + int power; + + // The constant term in the poly requires no multiplication. + sum = poly[0]; + + xLog = Logarithms[x]; + + for (int i = 1; i < poly.Length; i++) + { + // The polynomial at this spot has some coefficent, call it a^j. + // x itself has some value, call it a^k. + // x is raised to some power, which is 'i' in the loop. + // a^j * (a^k)^i + // == a^j * a^(k*i) + // == a^(j+k+i) + // Remember that exponents are added modulo 'size - 1', because the field elements + // are {0, a^0, a^1, ... } - we use size - 1 because we don't use 0. + + // If the coeff is 0, then this monomial contributes no value to the sum. + if (poly[i] == 0) + continue; + + coeffLog = Logarithms[poly[i]]; + + // Add the powers together, then lookup which a^L you ended up with. + power = (coeffLog + xLog * i) % (size - 1); + sum ^= Field[power + 1]; + } + + return sum; + } + + /// + /// Multiplies two polynomials of degrees X and Y to produce a single polynomial + /// of degree X + Y, where 'degree' means the exponent of the highest order monomial. + /// Eg, x^2 + x + 1 is of degree 2, and has 3 monomials. + /// + /// + /// + /// + public int[] PolyMult(int[] left, int[] right) + { + int[] result; + + // Lets say we have the following two polynomials: + // 3x^2 + 0x + 1 == int[]{1, 0, 3} + // 1x^2 + 1x + 0 == int[]{0, 1, 1} + // + // The naive result (ignoring galois fields) will be: + // (3x^4 + 3x^3 + 0) + (0x^3 + 0x^2 + 0) + (1x^2 + 1x + 0) == + // 3x^4 + (3 + 0)x^3 + (0+1)x^2 + 1x + 0 == + // 3x^4 + 3x^3 + 1x^2 + 1x + 0 == {3, 3, 1, 1, 0} + + // The result is of degree X + Y == 2 + 2 = 4. The number of monomials, and thus the number of + // coefficients is degree + 1. Thus, for two polynomials with degree X and Y, the resultant + // polynomial has degree X + Y but has X + Y - 1 coefficents. + // This establishes our basis for the degree of the resultant polynomial. + + // Now, in our galois fields, coefficient multiplication and addition are different. + // In our GF(2^x), addition is simple XOR and multiplication is best represented as being logarithm based. + // If + // a^9 = 10 and a^13 = 13 in GF(16) p(x) = x^4 + x + 1, + // Then + // 10 * 13 == + // a^9 * a^13 == + // a^((9+13) mod 15) == + // a^(22 mod 15) == a^7 == + // 11 + // + // Taking another example: + // 5x + 1 == {1, 5} + // 7x + 10 == {10, 7} + // + // (5x^1 + 2x^0)(7x^1 + 10x^0) == + // [(5*7)x^2 + (5*10)x^1 ] + [(2*7)x^1 + (2*10)x^0] + // + // --> Perform multiplications: + // 5*7 = a^8 * a^10 = a^18 = a^3 = 8 + // 5*10 = a^8 * a^9 = a^17 = a^2 = 4 + // 2*7 = a^1 * a^10 = a^11 = a^11 = 14 + // 2*10 = a^1 * a^9 = a^10 = a^10 = 7 + // [8x^2 + 4x^1 ] + [14x^1 + 7x^0] + // + // --> Combine like terms + // 8x^2 + (4+14)x^1 + 7x^0 + // + // --> Perform sums + // 4+14 = 4 XOR 14 = 10 + // 8x^2 + 10x^1 + 7x^0 + // + // Done. + + int coeff; + result = new int[left.Length + right.Length - 1]; + + for (int i = 0; i < left.Length; i++) + { + for (int j = 0; j < right.Length; j++) + { + coeff = InternalMult(left[i], right[j]); + result[i + j] = result[i + j] ^ coeff; + } + } + + return result; + } + + /// + /// + /// + private void BuildField() + { + int next; + int last; + + this.Field[0] = 0; + this.Field[1] = 1; + last = 1; + + for (int i = 2; i < this.Field.Length; i++) + { + next = last << 1; + if (next >= size) + next = next ^ fieldGenPoly; + + this.Field[i] = next; + last = next; + } + } + + /// + /// + /// + private void BuildInverses() + { + // Build the mulitiplicative inverses. + // if a^9 = 10, then what is 1/a^9 aka a^(-9) ? + this.Inverses[0] = 0; + for (int i = 1; i < this.Inverses.Length; i++) + this.Inverses[this.Field[i]] = InternalDivide(1, this.Field[i]); + } + + /// + /// + /// + private void BuildLogarithms() + { + // In GF(2^8) with p(x) = x^4 + x + 1, the field has elements 0, a^0, .., a^15: + // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 + // { 0, 1, 2, 4, 8, 3, 6, 12, 11, 5, 10, 7, 14, 15 13, 9} + + // In the above, we have the elements of the field by their index. + // What about going from the element to its power of a, for multiplying? + // In above, we have element 15 at index 13, but 15 is a^12 - it's inverse is one higher + // than its power, because 0 is in there. + + // field[13] = 15; + // logarithms[15] = 13 - 1; + + // logarithms[ field[13] ] = 13 - 1; + // logarithms[ 15 ] = 13 - 1; + // logarithms[ 15 ] = 12; + + // This means that zero will be stored with a logarithm of -1. + // This is intentional, but we have to be careful to handle it specially when we actually + // do multiplication. + for (int i = 0; i < this.Field.Length; i++) + this.Logarithms[this.Field[i]] = i - 1; + } + + /// + /// + /// + private void BuildMultTable() + { + this.multTable = new int[this.size, this.size]; + for (int left = 0; left < size; left++) + for (int right = 0; right < size; right++) + this.multTable[left, right] = InternalMult(left, right); + } + + /// + /// + /// + /// + /// + /// + private int InternalDivide(int dividend, int divisor) + { + // Same general concept as Mult. Convert to logarithms and subtract. + if (dividend == 0) + return 0; + + dividend = Logarithms[dividend]; + divisor = Logarithms[divisor]; + + // Note the extra '... + size - 1' term. This is to push the subtraction above + // zero, so that the modulus operator will do the right thing. + // (10 - 11) % 5 == -1 wrong + // ((10 - 11) + 5) % 5 == 4 right + dividend = (dividend - divisor + (size - 1)) % (size - 1); + + return this.Field[dividend + 1]; + } + + /// + /// + /// + /// + /// + /// + private int InternalMult(int left, int right) + { + // Conceptual notes: + // If + // a^9 = 10 and a^13 = 13 in GF(16) p(x) = x^4 + x + 1, + // Then + // 10 * 13 == + // a^9 * a^13 == + // a^((9+13) mod 15) == + // a^(22 mod 15) == a^7 == + // 11 + // + // Implementation notes: + // Our logarithms table stores a^9 at logarithms[9] = 10; + // Our field table stores a^9 at field[10] (0 is the first element, a^0 is the second). + // a^i is stored at field[i+1]; + // + // Plan: + // Convert each field element to its logarithm: + // left = a^i --> i + // right = a^j --> j + // + // Sum the logarithms to perform the multiplication. + // Modulus the sum to convert from, eg, a^15 back to a^0 + // k = (i + j) mod (size-1). + // + // Convert k back to a^k. Remember that a^k is stored at field[k+1]; + // a^k = field[k+1]; + // + // Return a^k. + + // Handle the special case 0 + if (left == 0 || right == 0) + return 0; + + // Convert each to their logarithm; + left = Logarithms[left]; + right = Logarithms[right]; + + // Sum the logarithms, using left to store the result. + left = (left + right) % (size - 1); + + // Convert the logarithm back to the field value. + return this.Field[left + 1]; + } + } // public sealed class GaloisField +} // namespace fnecore.EDAC.RS diff --git a/EDAC/RS/ReedSolomonAlgorithm.cs b/EDAC/RS/ReedSolomonAlgorithm.cs index df3920e..2051940 100644 --- a/EDAC/RS/ReedSolomonAlgorithm.cs +++ b/EDAC/RS/ReedSolomonAlgorithm.cs @@ -30,7 +30,7 @@ namespace fnecore.EDAC.RS /// /// (63,47,17) /// - ReedSolomon_634717, + ReedSolomon_362017, /// /// (24,12,13) /// @@ -62,7 +62,7 @@ namespace fnecore.EDAC.RS byte[] output = new byte[message.Length]; Buffer.BlockCopy(message, 0, output, 0, message.Length); - if (eccType == ErrorCorrectionCodeType.ReedSolomon_634717) + if (eccType == ErrorCorrectionCodeType.ReedSolomon_362017) { reedSolomonEncoder.encode362017(ref output); return output; @@ -89,7 +89,92 @@ namespace fnecore.EDAC.RS /// Returns the repaired message, or null if it cannot be repaired. public static byte[] Decode(byte[] message, ErrorCorrectionCodeType eccType) { - throw new NotImplementedException(); + byte[] codeword = new byte[63]; + + if (eccType == ErrorCorrectionCodeType.ReedSolomon_362017) + { + uint offset = 0U; + for (uint i = 0U; i < 36U; i++, offset += 6) + codeword[27 + i] = bin2Hex(message, offset); + + ReedSolomonDecoder rs362017 = new ReedSolomonDecoder(64, 47, 16, 0x43); + byte[] codewordEC = rs362017.DecodeEx(codeword); + + byte[] messageOut = new byte[message.Length]; + offset = 0U; + for (uint i = 0U; i < 20U; i++, offset += 6) + hex2Bin(codewordEC[27 + i], ref messageOut, offset); + return messageOut; + } + else if (eccType == ErrorCorrectionCodeType.ReedSolomon_241213) + { + uint offset = 0U; + for (uint i = 0U; i < 24U; i++, offset += 6) + codeword[39 + i] = bin2Hex(message, offset); + + ReedSolomonDecoder rs241213 = new ReedSolomonDecoder(64, 51, 12, 0x43); + byte[] codewordEC = rs241213.DecodeEx(codeword); + + byte[] messageOut = new byte[message.Length]; + offset = 0U; + for (uint i = 0U; i < 12U; i++, offset += 6) + hex2Bin(codewordEC[39 + i], ref messageOut, offset); + return messageOut; + } + else if (eccType == ErrorCorrectionCodeType.ReedSolomon_24169) + { + uint offset = 0U; + for (uint i = 0U; i < 24U; i++, offset += 6) + codeword[39 + i] = bin2Hex(message, offset); + + ReedSolomonDecoder rs24169 = new ReedSolomonDecoder(64, 55, 8, 0x43); + byte[] codewordEC = rs24169.DecodeEx(codeword); + + byte[] messageOut = new byte[message.Length]; + offset = 0U; + for (uint i = 0U; i < 16U; i++, offset += 6) + hex2Bin(codewordEC[39 + i], ref messageOut, offset); + return messageOut; + } + else + throw new ArgumentException($"Invalid '{nameof(eccType)}' argument.", nameof(eccType)); + } + + /// + /// + /// + /// + /// + /// + internal static byte bin2Hex(byte[] input, uint offset) + { + byte output = 0x00; + + output |= (byte)(FneUtils.ReadBit(input, offset + 0U) ? 0x20U : 0x00U); + output |= (byte)(FneUtils.ReadBit(input, offset + 1U) ? 0x10U : 0x00U); + output |= (byte)(FneUtils.ReadBit(input, offset + 2U) ? 0x08U : 0x00U); + output |= (byte)(FneUtils.ReadBit(input, offset + 3U) ? 0x04U : 0x00U); + output |= (byte)(FneUtils.ReadBit(input, offset + 4U) ? 0x02U : 0x00U); + output |= (byte)(FneUtils.ReadBit(input, offset + 5U) ? 0x01U : 0x00U); + + return output; + } + + /// + /// + /// + /// + /// + /// + /// + internal static void hex2Bin(byte input, ref byte[] output, uint offset) + { + FneUtils.WriteBit(ref output, offset + 0U, (input & 0x20U) == 0x20U); + FneUtils.WriteBit(ref output, offset + 1U, (input & 0x10U) == 0x10U); + FneUtils.WriteBit(ref output, offset + 2U, (input & 0x08U) == 0x08U); + FneUtils.WriteBit(ref output, offset + 3U, (input & 0x04U) == 0x04U); + FneUtils.WriteBit(ref output, offset + 4U, (input & 0x02U) == 0x02U); + FneUtils.WriteBit(ref output, offset + 5U, (input & 0x01U) == 0x01U); } } // public static class ReedSolomonAlgorithm } // namespace fnecore.EDAC.RS diff --git a/EDAC/RS/ReedSolomonDecoder.cs b/EDAC/RS/ReedSolomonDecoder.cs new file mode 100644 index 0000000..9611ce9 --- /dev/null +++ b/EDAC/RS/ReedSolomonDecoder.cs @@ -0,0 +1,463 @@ +/** +* Digital Voice Modem - Fixed Network Equipment +* AGPLv3 Open Source. Use is subject to license terms. +* DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. +* +* @package DVM / Fixed Network Equipment +* +*/ +// +// Based on code from the MMDVMHost project. (https://github.com/g4klx/MMDVMHost) +// Licensed under the GPLv2 License (https://opensource.org/licenses/GPL-2.0) +// +/* +* Copyright (C) 2016 by Jonathan Naylor G4KLX +* Copyright (C) 2017-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 +* 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. +*/ + +using System; +using System.Linq; + +namespace fnecore.EDAC.RS +{ + /// + /// Implements Reed-Solomon decoding, as the name implies. + /// + /// This does not do error correction, only detection. + internal sealed class ReedSolomonDecoder + { + private readonly int fieldSize; + private readonly int messageSymbols; + private readonly int paritySymbols; + + private readonly int fieldGenPoly; + + private readonly GaloisField gf; + + private readonly int[] syndroms; + + private readonly int[] lambda; + private readonly int[] corrPoly; + private readonly int[] lambdaStar; + + private readonly int[] lambdaPrime; + + private readonly int[] omega; + + private readonly int[] errorIndexes; + + private readonly int[] chienCache; + + /* + ** Properties + */ + + /// + /// The number of symbols that make up an entire encoded message. An encoded message is composed of the + /// original data symbols plus parity symbols. + /// + public int BlockSize { get; private set; } + + /// + /// The number of symbols per block that store original message symbols. + /// + public int MessageSize + { + get { return this.messageSymbols; } + } + + /* + ** Methods + */ + + /// + /// Initializes a new instance of the class. + /// + /// The size of the Galois field to create. Must be a value that is + /// a power of two. The length of the output block is set to `fieldSize - 1`. + /// The number of original message symbols per block. + /// The number of parity symbols per block. + /// A value representing the field generator polynomial, + /// which must be order N for a field GF(2^N). + /// BlockSize is equal to `fieldSize - 1`. messageSymbols plus paritySymbols must equal BlockSize. + public ReedSolomonDecoder(int fieldSize, int messageSymbols, int paritySymbols, int fieldGenPoly) + { + if (fieldSize - 1 != messageSymbols + paritySymbols) + throw new ArgumentOutOfRangeException("Invalid reed-solomon block parameters were provided - the number of message symbols plus the number of parity symbols " + + "does not add up to the size of a block"); + + this.fieldSize = fieldSize; + this.messageSymbols = messageSymbols; + this.paritySymbols = paritySymbols; + this.BlockSize = fieldSize - 1; + + this.fieldGenPoly = fieldGenPoly; + + this.gf = new GaloisField(fieldSize, fieldGenPoly); + + // Syndrom calculation buffers + this.syndroms = new int[paritySymbols]; + + // Lamda calculation buffers + this.lambda = new int[paritySymbols - 1]; + this.corrPoly = new int[paritySymbols - 1]; + this.lambdaStar = new int[paritySymbols - 1]; + + // LambdaPrime calculation buffers + this.lambdaPrime = new int[paritySymbols - 2]; + + // Omega calculation buffers + this.omega = new int[paritySymbols - 2]; + + // Error position calculation + this.errorIndexes = new int[fieldSize - 1]; + + // Cache of the lookup used in the ChienSearch process. + this.chienCache = new int[fieldSize - 1]; + + for (int i = 0; i < this.chienCache.Length; i++) + this.chienCache[i] = gf.Inverses[gf.Field[i + 1]]; + } + + /// + /// Discovers and corrects any errors in the block provided. + /// + /// A block containing a reed-solomon encoded message. + public byte[] DecodeEx(byte[] message) + { + int[] received = message.Select(x => (int)x).ToArray(); + Decode(ref received); + + return received.Take(message.Length).Select(x => (byte)x).ToArray(); + } + + /// + /// Discovers and corrects any errors in the block provided. + /// + /// A block containing a reed-solomon encoded message. + public void Decode(ref int[] message) + { + if (message.Length != this.BlockSize) + throw new ArgumentException("The provided message's size was not the size of a block"); + + CalcSyndromPoly(message); + CalcLambda(); + CalcLambdaPrime(); + CalcOmega(); + + ChienSearch(); + + RepairErrors(ref message, errorIndexes, omega, lambdaPrime); + } + + /// + /// + /// + /// + /// + /// + /// + private void RepairErrors(ref int[] message, int[] errorIndexes, int[] omega, int[] lp) + { + int top, bottom, x, xInverse; + int messageLen = message.Length; + + for (int i = 0; i < messageLen; i++) + { + // If i = 2, then use a^2 in the evaluation. + // remember that field[i + 1] = a^i, since field[0] = 0 and field[1] = a^0. + + // i = 2 + // a^2 is the element. a^2 = 4 + // + // 13 = 4 * [6*(4^-1)^2 + 15] / 14 + + // 8 ^ 15 = 7 + // 6 * a^-2 = 8 + // 8 + 15 = 7 + + if (errorIndexes[i] == 0) + { + // This spot has an error. + // Equation: + // value = X_J * [ Omega(1/X_J) / L'(1/X_J) ]; + // Break it down to: + // top = eval(omega, 1/X_J); + // top = X_J * top; + // bottom = eval( lambaPrime, 1/X_J ); + // errorMagnitude = top / bottom; + // + // To repair the message, we add the error magnitude to the received value. + // message[i] = message[i] ^ errorMagnitude + + x = gf.Field[i + 1]; + + xInverse = gf.Inverses[x]; + + top = gf.PolyEval(omega, xInverse); + top = gf.Multiply(top, x); + bottom = gf.PolyEval(lp, xInverse); + + message[i] ^= gf.Divide(top, bottom); + } + } + } + + /// + /// + /// + private void CalcLambda() + { + // Explanation of terms: + // S = S(x) - syndrom polynomial + // C = C(x) - correction polynomial + // D = 'lamba'(x) - error locator estimate polynomial. + // D* = 'lambda-star'(x) - a new error locator estimate polynomial. + // S_x = the element at index 'x' in S, eg, if S = {5,6,7,8}, then S_0 = 5, S_1 = 6, etc. + // 2T = the number of error correction symbols, eg, numCheckBytes. + // T must be >= 1, so 2T is guarenteed to be at least 2. + // + // Start with + // K = 1; + // L = 0; + // C = 0x^n + ... + 0x^2 + x + 0 aka {0,1,0, ...}; + // D = 0x^n + ... + 0x^2 + 0x + 1 aka {1,0,0,...}; + // Both C and D are guarenteed to be at least 2 elements, which is why they can have + // hardcoded initial values. + + // Step 1: Calculate e. + // -------------- + // e = S_(K-1) + sum(from i=1 to L: D_i * S_(K-1-i) + // + // Example + // 0 1 2 0 1 2 3 + // K = 4, L = 2; D = {1, 11, 15}; S = {15, 3, 4, 12} + // + // e = S_3 + sum(from i = 1 to 2: D_i * S_(3- i) + // = 12 + D_1 * S_2 + D_2 * S_1 + // = 12 + 11 * 4 + 15 * 3 + // = 12 + 10 + 2 + // = 12 XOR 10 XOR 2 + // = 4 + + // Step 2: Update estimate if e != 0 + // + // If e != 0 { + // D*(x) = D(x) + e * C(X) -- Note that this just assigns D(x) to D*(x) if e is zero. + // If 2L < k { + // L = K - L + // C(x) = D(x) * e^(-1) -- Multiply D(x) times the multiplicative inverse of e. + // } + // } + + // Step 3: Advance C(x): + // C(x) = C(x) * x + // This just shifts the coeffs down; eg, x + 1 {1, 1, 0} turns into x^2 + x {0, 1, 1} + // + // D(x) = D*(x) (only if a new D*(x) was calulated) + // + // K = K + 1 + + // Step 4: Compute end conditions + // If K <= 2T goto 1 + // Else, D(x) is the error locator polynomial. + + int k, l, e; + int eInv; // temp to store calculation of 1 / e aka e^(-1) + + // --- Initial conditions ---- + // Need to clear lambda and corrPoly, but not lambdaStar. lambda and corrPoly + // are used and initialized iteratively in the algorithm, whereas lambdaStar isn't. + Array.Clear(corrPoly, 0, corrPoly.Length); + Array.Clear(lambda, 0, lambda.Length); + k = 1; + l = 0; + corrPoly[1] = 1; + lambda[0] = 1; + + while (k <= paritySymbols) + { + // --- Calculate e --- + e = syndroms[k - 1]; + + for (int i = 1; i <= l; i++) + e ^= gf.Multiply(lambda[i], syndroms[k - 1 - i]); + + // --- Update estimate if e != 0 --- + if (e != 0) + { + // D*(x) = D(x) + e * C(x); + for (int i = 0; i < lambdaStar.Length; i++) + lambdaStar[i] = lambda[i] ^ gf.Multiply(e, corrPoly[i]); + + if (2 * l < k) + { + // L = K - L; + l = k - l; + + // C(x) = D(x) * e^(-1); + eInv = gf.Inverses[e]; + for (int i = 0; i < corrPoly.Length; i++) + corrPoly[i] = gf.Multiply(lambda[i], eInv); + } + } + + // --- Advance C(x) --- + + // C(x) = C(x) * x + for (int i = corrPoly.Length - 1; i >= 1; i--) + corrPoly[i] = corrPoly[i - 1]; + corrPoly[0] = 0; + + if (e != 0) + { + // D(x) = D*(x); + Array.Copy(lambdaStar, lambda, lambda.Length); + } + + k += 1; + } + } + + /// + /// + /// + private void CalcLambdaPrime() + { + // Forney's says that we can just set even powers to 0 and then take the rest and + // divide it by x (shift it down one). + + // No need to clear this.lambdaPrime between calls; full assignment is done every call. + for (int i = 0; i < lambdaPrime.Length; i++) + if ((i & 0x1) == 0) + lambdaPrime[i] = lambda[i + 1]; + else + lambdaPrime[i] = 0; + } + + /// + /// + /// + private void CalcOmega() + { + // O(x) is shorthand for Omega(x). + // L(x) is shorthand for Lambda(x). + // + // O_i is the coefficient of the term in omega with degree i. Ditto for L_i. + // Eg, O(x) = 6x + 15; O_0 = 15, O_1 = 6 + // + // From the paper: + // O_0 = S_b + // ---> b in our implementation is 0. + // O_1 = S_{b+1} + S_b * L_1 + // + // O_{v-1} = S_{b+v-1} + S_{b+v-2} * L_1 + ... + S_b * L_{v-1}. + // O_i = S_{b+i} + S_{b+ i-1} * L_1 + ... + S_{b+0} * L_i + + // Lets say : + // L(x) = 14x^2 + 14x + 1 aka {1, 14, 14}. + // S(x) = 12x^3 + 4x^2 + 3x + 15 aka {15, 3, 4, 12} + // b = 0; + // v = 2 because the power of the highest monomial in L(x), 14x^2, is 2. + // + // O_0 = S_{b} = S_0 = 15 + // O_1 = S_{b+1} + S_b * L_1 = S_1 + S_0 * L_1 = 3 + 15 * 14 = 6. + // + // O(x) = 6x + 15. + + // Lets make up another example (these are completely made up so they may not work): + // L(x) = 10x^3 + 9x^2 + 8x + 7 aka { 7, 8, 9, 10} + // S(x) = 2^4 + 3x^3 + 4x^2 + 5x + 6 aka { 6, 5, 4, 3, 2} + // b = 0 (design parameter) + // v = 3 + + // O_i for i = 0 .. v - 1 = 2. Thus, O has form ax^2 + bx^1 + cx^0 + // Compute O_0, O_1, O_2 + + // O_0 = S_{b+0} + // = S_0 + // + // O_1 = S_{b+1} + S_{b+0} * L_1 + // = S_1 + S_0 * L_1 + // + // O_2 = S_{b+2} + S_{b+1} * L_1 + S_{b+0} * L_2 + // + S_2 + S_1 * L_1 + S_0 * L_2 + + // Don't need to zero this.omega first - it's assigned to before we use it. + + for (int i = 0; i < omega.Length; i++) + { + omega[i] = syndroms[i]; + for (int lIter = 1; lIter <= i; lIter++) + omega[i] ^= gf.Multiply(syndroms[i - lIter], lambda[lIter]); + } + } + + /// + /// + /// + private void ChienSearch() + { + // The cheap chien search evaluates the lamba polynomial for the multiplicate inverse + // each element in the field other than 0. + + // Don't need to zero this.errorIndexes first - it's not used before its assigned to. + + /* + * This loop uses an optimization where I precalculate one lookup. + * The code before the optimization was: + * for( int i = 0; i < errorIndexes.Length; i++ ) + * { + * errorIndexes[i] = gf.PolyEval( + * lambda, + * gf.Inverses[ gf.Field[ i + 1] ] + * ); + * } + * We precompute the lookup gf.Inverses[ gf.Field[ i + 1] ] when creating the decoder. + */ + + for (int i = 0; i < errorIndexes.Length; i++) + errorIndexes[i] = gf.PolyEval(lambda, chienCache[i]); + } + + /// + /// + /// + /// + private void CalcSyndromPoly(int[] message) + { + int syndrome, root; + + // Don't need to zero this.syndromes first - it's not used before its assigned to. + for (int synIndex = 0; synIndex < syndroms.Length; synIndex++) + { + // EG, if g(x) = (x+a^0)(x+a^1)(x+a^2)(x+a^3) + // = (x+1)(x+2)(x+4)(x+8), + // Then for the first syndrom S_0, we would provide root = a^0 = 1. + // S_1 --> root = a^1 = 2 etc. + + root = gf.Field[synIndex + 1]; + syndrome = 0; + + for (int i = message.Length - 1; i > 0; i--) + syndrome = gf.Multiply((syndrome ^ message[i]), root); + + syndroms[synIndex] = syndrome ^ message[0]; + } + } + } // internal sealed class ReedSolomonDecoder +} // namespace fnecore.EDAC.RS diff --git a/EDAC/RS/ReedSolomonEncoder.cs b/EDAC/RS/ReedSolomonEncoder.cs index f3c3feb..5222026 100644 --- a/EDAC/RS/ReedSolomonEncoder.cs +++ b/EDAC/RS/ReedSolomonEncoder.cs @@ -117,14 +117,14 @@ namespace fnecore.EDAC.RS uint offs = 0U; for (uint j = 0U; j < 12U; j++, offs += 6U) { - byte hexbit = bin2Hex(data, offs); + byte hexbit = ReedSolomonAlgorithm.bin2Hex(data, offs); codeword[i] ^= gf6Mult(hexbit, ENCODE_MATRIX[j][i]); } } uint offset = 0U; for (uint i = 0U; i < 24U; i++, offset += 6U) - hex2Bin(codeword[i], ref data, offset); + ReedSolomonAlgorithm.hex2Bin(codeword[i], ref data, offset); } /// @@ -145,14 +145,14 @@ namespace fnecore.EDAC.RS uint offs = 0U; for (uint j = 0U; j < 16U; j++, offs += 6U) { - byte hexbit = bin2Hex(data, offs); + byte hexbit = ReedSolomonAlgorithm.bin2Hex(data, offs); codeword[i] ^= gf6Mult(hexbit, ENCODE_MATRIX_24169[j][i]); } } uint offset = 0U; for (uint i = 0U; i < 24U; i++, offset += 6U) - hex2Bin(codeword[i], ref data, offset); + ReedSolomonAlgorithm.hex2Bin(codeword[i], ref data, offset); } /// @@ -173,51 +173,14 @@ namespace fnecore.EDAC.RS uint offs = 0U; for (uint j = 0U; j < 20U; j++, offs += 6U) { - byte hexbit = bin2Hex(data, offs); + byte hexbit = ReedSolomonAlgorithm.bin2Hex(data, offs); codeword[i] ^= gf6Mult(hexbit, ENCODE_MATRIX_362017[j][i]); } } uint offset = 0U; for (uint i = 0U; i < 36U; i++, offset += 6U) - hex2Bin(codeword[i], ref data, offset); - } - - /// - /// - /// - /// - /// - /// - private byte bin2Hex(byte[] input, uint offset) - { - byte output = 0x00; - - output |= (byte)(FneUtils.ReadBit(input, offset + 0U) ? 0x20U : 0x00U); - output |= (byte)(FneUtils.ReadBit(input, offset + 1U) ? 0x10U : 0x00U); - output |= (byte)(FneUtils.ReadBit(input, offset + 2U) ? 0x08U : 0x00U); - output |= (byte)(FneUtils.ReadBit(input, offset + 3U) ? 0x04U : 0x00U); - output |= (byte)(FneUtils.ReadBit(input, offset + 4U) ? 0x02U : 0x00U); - output |= (byte)(FneUtils.ReadBit(input, offset + 5U) ? 0x01U : 0x00U); - - return output; - } - - /// - /// - /// - /// - /// - /// - /// - private void hex2Bin(byte input, ref byte[] output, uint offset) - { - FneUtils.WriteBit(ref output, offset + 0U, (input & 0x20U) == 0x20U); - FneUtils.WriteBit(ref output, offset + 1U, (input & 0x10U) == 0x10U); - FneUtils.WriteBit(ref output, offset + 2U, (input & 0x08U) == 0x08U); - FneUtils.WriteBit(ref output, offset + 3U, (input & 0x04U) == 0x04U); - FneUtils.WriteBit(ref output, offset + 4U, (input & 0x02U) == 0x02U); - FneUtils.WriteBit(ref output, offset + 5U, (input & 0x01U) == 0x01U); + ReedSolomonAlgorithm.hex2Bin(codeword[i], ref data, offset); } ///