Reed–Solomon codes for coders

From Wikiversity
Jump to: navigation, search

Reed–Solomon error correction is widely used in applications such as digital data storage (like CDs) and transmission. However, in these situations, the codes themselves are hidden inside an electronic device, so there's no opportunity to take a peek at them and see how they work. Some of the more complicated barcode designs also incorporate Reed–Solomon codes, which puts everything in plain view. This makes them an interesting subject for a hobbyist who wants to see firsthand how this kind of technology works.

In this essay, I will attempt to introduce the principles of Reed–Solomon codes from the point of view of a programmer rather than a mathematician. I will provide real-world examples taken from the popular QR code barcode system as well as working code samples. I chose to use Python for the samples (mainly because it looks pretty), but I will try to explain any non-obvious features for those who are not familiar with it. The math involved is advanced in the sense that it is not usually taught below the university level, but it should be understandable to someone with a good grasp of high-school algebra.

QR code structure[edit]

This section introduces the structure of QR codes in detail. The information in this section is deliberately incomplete. Only the more common features of the small 21×21 size symbols (also known as version 1) are presented here, but see the appendix for additional information.

Here is a QR symbol that will be used as an example. It consists of dark and light squares, known as modules in the barcoding world. The three square locator patterns in the corners are a visually distinctive feature of QR symbols.

QR Code Example.svg


A masking process is used to avoid features in the symbol that might confuse a scanner, such as misleading shapes that look like the locator patterns and large blank areas. Masking inverts certain modules (white becomes black and black becomes white) while leaving others alone.

In the diagram below, the red areas encode format information and use a fixed masking pattern. The data area (in black and white) is masked with a variable pattern. When the code is created, the encoder tries a number of different masks and chooses the one that minimizes undesirable features in the result. The chosen mask pattern is then indicated in the format information so that the decoder knows which one to use. The light gray areas are fixed patterns which do not encode any information. In addition to the obvious locator patterns, there are also timing patterns which contain alternating light and dark modules.

QR Code Masking Example.svg

The masking transformation is easily applied (or removed) using the exclusive-or operation (denoted by a caret ^ in many programming languages). The unmasking of the format information is shown below. Reading counter-clockwise around the upper-left locator pattern, we have the following sequence of bits. White modules represent 0 and black modules represent 1.

Input       101101101001011
Mask      ^ 101010000010010
Output      000111101011001

Format information[edit]

There are two identical copies of the format information, so that the symbol can still be decoded even if it is damaged. The second copy is broken in two pieces and placed around the other two locators, and is also read in a counter-clockwise direction (upwards in the lower-left corner, then left-to-right in the upper-right corner).

The first two bits of format information give the error correction level used for the message data. A QR symbol this size contains 26 bytes of information. Some of these are used to store the message and some are used for error correction, as shown in the table below. The left-hand column is simply a name given to that level.

Error Correction Level Level Indicator Error Correction Bytes Message Data Bytes
L 01 7 19
M 00 10 16
Q 11 13 13
H 10 17 9

The next three bits of format information select the masking pattern to be used in the data area. The patterns are illustrated below, including the mathematical formula that tells whether a module is black (i and j are the row and column numbers, respectively, and start with 0 in the upper-left hand corner).

QR Code Mask Patterns.svg

The remaining ten bits of format information are for correcting errors in the format itself. This will be explained in a later section.

Message data[edit]

Here is a larger diagram showing the "unmasked" QR code. Different regions of the symbol are indicated, including the boundaries of the message data bytes.

QR Code Unmasked.svg

Data bits are read starting from the lower-right corner and moving up the two right-hand columns in a zig-zag pattern. The first three bytes are 01000000 11010010 01110101. The next two columns are read in a downward direction, so the next byte is 01000111. Upon reaching the bottom, the two columns after that are read upward. Proceed in this up-and-down fashion all the way to the left side of the symbol (skipping over the timing pattern where necessary). Here is the complete message in hexadecimal notation.

Message data bytes: 40 d2 75 47 76 17 32 06 27 26 96 c6 c6 96 70 ec
Error correction bytes: bc 2a 90 13 6b af ef fd 4b e0


The final step is to decode the message bytes into something readable. The first four bits indicate how the message is encoded. QR codes use several different encoding schemes, so that different kinds of messages can be stored efficiently. These are summarized in the table below. After the mode indicator is a length field, which tells how many characters are stored. The size of the length field depends on the specific encoding.

Mode Name Mode Indicator Length Bits Data Bits
Numeric 0001 10 10 bits per 3 digits
Alphanumeric 0010 9 11 bits per 2 characters
Byte 0100 8 8 bits per character
Kanji 1000 8 13 bits per character

(The length field sizes above are valid only for smaller QR codes.)

Our sample message starts with 0100, indicating that there are 8 bits per character. The next 8 bits are the length field, 00001101, or 13 in decimal notation. After that are the actual characters of the message. The first two are 00100111 and 01010100 (the ASCII codes for apostrophe and T). Interested readers may want to decode the rest of the message for themselves.

After the last of the data bits is another 4-bit mode indicator. It can be different from the first one, allowing different encodings to be mixed within the same QR symbol. When there is no more data to store, the special end-of-message code 0000 is given. (Note that the standard allows the end-of-message code to be omitted if it wouldn't fit in the available number of data bytes.)

BCH codes[edit]

The format information is encoded with a BCH code which allows a certain number of bit-errors to be detected and corrected. BCH codes are a generalization of Reed–Solomon codes (all Reed–Solomon codes are also BCH codes). In the case of QR codes, the BCH code used for the format information is much simpler than the Reed–Solomon code used for the message data, so it makes sense to start here.

BCH error detection[edit]

The process for checking the encoded information is similar to long division, but uses exclusive-or instead of subtraction. The format code should produce a remainder of zero when it is is "divided" by the so-called generator of the code. QR format codes use the generator 10100110111. This process is demonstrated for the format information in the example code (000111101011001) below.

10100110111 ) 000111101011001
               ^ 10100110111 
                ^ 10100110111

Here is a Python function which implements this calculation.

def qr_check_format(fmt):
   g = 0x537 # = 0b10100110111 in python 2.6+
   for i in range(4,-1,-1):
      if fmt & (1 << (i+10)):
         fmt ^= g << i
   return fmt

Python note: The range function may not be clear to non-Python programmers. It produces a list of numbers counting down from 4 to 0. In C-derived languages, the for loop might be written as for (i = 4; i >= 0; i--); in Pascal-derived languages, for i := 4 downto 0.

Python note 2: The & operator performs bitwise and, while << is a left bit-shift. This is consistent with C-like languages.

This function can also be used to encode the 5-bit format information.

encoded_format = (format<<10) ^ qr_check_format(format<<10)

Readers may find it an interesting exercise to generalize this function to divide by different numbers. For example, larger QR codes contain six bits of version information with 12 error correction bits using the generator 1111100100101.

In mathematical formalism, these binary numbers are described as polynomials whose coefficients are integers mod 2. Each bit of the number is a coefficient of one term. For example:

10100110111 = 1 x10 + 0 x9 + 1 x8 + 0 x7 + 0 x6 + 1 x5 + 1 x4 + 0 x3 + 1 x2 + 1 x + 1 = x10 + x8 + x5 + x4 + x2 + x + 1

BCH error correction[edit]

If the remainder produced by qr_check_format is not zero, then the code has been damaged or misread. The next step is to determine which format code is most likely the one that was intended. Although sophisticated algorithms for decoding BCH codes exist, they are probably overkill in this case. Since there are only 32 possible format codes, it's much easier to simply try each one and pick the one that has the smallest number of bits different from the code in question (this is known as the Hamming distance).

def hamming_weight(x):
   weight = 0
   while x > 0:
      weight += x & 1
      x >>= 1
   return weight
def qr_decode_format(fmt):
   best_fmt = -1
   best_dist = 15
   for test_fmt in range(0,32):
      test_code = (test_fmt<<10) ^ qr_check_format(test_fmt<<10)
      test_dist = hamming_weight(fmt ^ test_code)
      if test_dist < best_dist:
         best_dist = test_dist
         best_fmt = test_fmt
      elif test_dist == best_dist:
         best_fmt = -1
   return best_fmt

The function qr_decode_format returns -1 if the format code could not be unambiguously decoded. This happens when two or more format codes have the same distance from the input.

To run this code in Python, first start IDLE, Python's integrated development environment. You should see a version message and the interactive input prompt >>>. Open a new window, copy the functions qr_check_format, hamming_weight, and qr_decode_format into it, and save as Return to the prompt and type the lines following >>> below.

>>> from qr import *
>>> qr_decode_format(int("000111101011001",2))  # no errors
>>> qr_decode_format(int("111111101011001",2))  # 3 bit-errors
>>> qr_decode_format(int("111011101011001",2))  # 4 bit-errors

You can also start Python by typing python at a command prompt.

Finite field arithmetic[edit]

Before discussing the Reed–Solomon codes used for the message, it will be useful to introduce a bit more math.

We'd like to define addition, subtraction, multiplication, and division for 8-bit bytes and always produce 8-bit bytes as a result. Naively, we might attempt to use the normal definitions for these operations, and then mod by 256 to keep results from overflowing. This works for everything but division: what's 7/5?

Instead, we throw out those definitions and introduce new ones, inspired by our polynomial analogy above. These new operations will let us define what's known as a finite field or Galois field.

Addition and Subtraction[edit]

Both addition and subtraction are replaced with exclusive-or. Thinking of our 8-bit values as polynomials with coefficients mod 2:

   0101 + 0110 = 0011

The same way

(x2 + 1) + (x2 + x) = 2 x2 + x + 1 = 0 x2 + x + 1 = x + 1

Since (a ^ a) = 0, every number is its own opposite, so (x - y) is the same as (x + y).


Multiplication is likewise based on polynomial multiplication. Simply write the inputs as polynomials and multiply them out using the distributive law as normal. As an example, 10001001 times 00101010 is calculated as follows.

(x7 + x3 + 1) (x5 + x3 + x) = x7 (x5 + x3 + x) + x3 (x5 + x3 + x) + 1 (x5 + x3 + x)
= x12 + x10 + 2 x8 + x6 + x5 + x4 + x3 + x
= x12 + x10 + x6 + x5 + x4 + x3 + x

The same result can be obtained by a modified version of the standard grade-school multiplication procedure, in which we replace addition with exclusive-or.

*      00101010
^   10001001
^ 10001001

Here is a Python function which implements this polynomial multiplication.

def poly_mult(x,y):
   z = 0
   i = 0
   while (y>>i) > 0:
      if y & (1<<i):
         z ^= x<<i
      i += 1
   return z

Of course, the result no longer fits in an 8-bit byte (in this example, it is 13 bits long), so we need to perform one more step before we are finished. The result is reduced modulo 100011101, using the long division process described previously. This produces the final answer 11000011 in our example.

^ 100011101
  ^ 100011101
    ^ 100011101

Why mod 100011101? The math is a little complicated here, but in short, 100011101 represents an 8th degree polynomial which is "irreducible" (meaning it can't represented as the product of two smaller polynomials). This is necessary for division to be well-behaved. There are other numbers we could have chosen, but they're all essentially the same, and 100011101 is what Reed-Solomon codes are defined to use.

Multiplication with logarithms[edit]

The procedure described above is not the most convenient way to implement Galois field multiplication. Multiplying two numbers takes up to eight iterations of the multiplication loop, followed by up to eight iterations of the division loop. However, we can multiply with no looping by using lookup tables. One solution would be to construct the entire multiplication table in memory, but that would require a bulky 64k table. The solution described below is much more compact.

First, notice that it is particularly easy to multiply by 00000010 (by convention, this number is referred to as α): simply left-shift by one place, then exclusive-or with the modulus 100011101 if necessary (why xor is sufficient for taking the mod in this case is an exercise left to the reader). Here are the first few powers of α.

α0 = 00000001 α4 = 00010000 α8  = 00011101 α12 = 11001101
α1 = 00000010 α5 = 00100000 α9  = 00111010 α13 = 10000111
α2 = 00000100 α6 = 01000000 α10 = 01110100 α14 = 00010011
α3 = 00001000 α7 = 10000000 α11 = 11101000 α15 = 00100110

If this table is continued in the same fashion, the powers of α do not repeat themselves until α255 = 00000001. Thus, every element of the field except zero is equal to some power of α. The element α is known as a primitive element or generator of the Galois field.

This observation suggests another way to implement multiplication: by adding the exponents of α.

10001001 * 00101010 = α74 * α142 = α74 + 142 = α216 = 11000011

The problem is, how do we find the power of α that corresponds to 10001001? This is known as the discrete logarithm problem, and no efficient general solution is known. However, since there are only 256 elements in this field, we can easily construct a table of logarithms. While we're at it, a corresponding table of antilogs (exponentials) will also be useful.

gf_exp = [0] * 512 # Create list of 512 elements. In Python 2.6+, consider using bytearray
gf_log = [0] * 256
gf_exp[0] = 1
x = 1
for i in range(1,255):
   x <<= 1
   if x & 0x100:
      x ^= 0x11d
   gf_exp[i] = x
   gf_log[x] = i
for i in range(255,512):
   gf_exp[i] = gf_exp[i-255]
gf_log[gf_exp[255]] = 255 # Set last missing elements in gf_log[]

Python note: The range operator's upper bound is exclusive, so gf_exp[255] is not set twice by the above.

The gf_exp table is oversized in order to simplify the multiplication function. This way, we don't have to check to make sure that gf_log[x] + gf_log[y] is within the table size.

def gf_mul(x,y):
   if x==0 or y==0:
      return 0
   return gf_exp[gf_log[x] + gf_log[y]]


Another advantage of the logarithm table approach is that it allows us to define division using the difference of logarithms. In the code below, 255 is added to make sure the difference isn't negative.

def gf_div(x,y):
   if y==0:
      raise ZeroDivisionError()
   if x==0:
      return 0
   return gf_exp[gf_log[x] + 255 - gf_log[y]]

Python note: The raise statement throws an exception and aborts execution of the gf_div function.

With this definition of division, gf_div(gf_mul(x,y),y)==x for any x and any nonzero y.

Readers who are more advanced programmers may find it interesting to write a class encapsulating Galois field arithmetic. Operator overloading can be used to replace calls to gf_mul and gf_div with the familiar operators * and /, but this can lead to confusion as to exactly what type of operation is being performed. Certain details can be generalized in ways that would make the class more widely useful. For example, Aztec codes use five different Galois fields with element sizes ranging from 4 to 12 bits.


Before moving on to Reed–Solomon codes, we need to define several operations on polynomials whose coefficients are Galois field elements. This is a potential source of confusion, since the elements themselves are described as polynomials; my advice is not to think about it too much. Adding to the confusion is the fact that x is still used as the placeholder. This x has nothing to do with the x mentioned previously, so don't mix them up.

The binary notation used previously for Galois field elements starts to become inconveniently bulky at this point, so I will switch to hexadecimal instead.

00000001 x4 + 00001111 x3 + 00110110 x2 + 01111000 x + 01000000 = 01 x4 + 0f x3 + 36 x2 + 78 x + 40

In Python, polynomials will be represented by a list of numbers in descending order of powers of x, so the polynomial above becomes [ 0x01, 0x0f, 0x36, 0x78, 0x40 ]. (The reverse order could have been used instead; both choices have their advantages and disadvantages.)

The first function multiplies a polynomial by a scalar.

def gf_poly_scale(p,x):
   r = [0] * len(p)
   for i in range(0, len(p)):
      r[i] = gf_mul(p[i], x)
   return r

Note to Python programmers: This function is not written in a "pythonic" style. It could be expressed quite elegantly as a list comprehension, but I have limited myself to language features that are easier to translate to other programming languages.

This function "adds" two polynomials (using exclusive-or, as usual).

def gf_poly_add(p,q):
   r = [0] * max(len(p),len(q))
   for i in range(0,len(p)):
      r[i+len(r)-len(p)] = p[i]
   for i in range(0,len(q)):
      r[i+len(r)-len(q)] ^= q[i]
   return r

The next function multiplies two polynomials.

def gf_poly_mul(p,q):
   r = [0] * (len(p)+len(q)-1)
   for j in range(0, len(q)):
      for i in range(0, len(p)):
         r[i+j] ^= gf_mul(p[i], q[j])
   return r

Finally, we need a function to evaluate a polynomial at a particular value of x, producing a scalar result. Horner's method is used to avoid explicitly calculating powers of x.

01 x4 + 0f x3 + 36 x2 + 78 x + 40 = (((01 x + 0f) x + 36) x + 78) x + 40
def gf_poly_eval(p,x):
   y = p[0]
   for i in range(1, len(p)):
      y = gf_mul(y,x) ^ p[i]
   return y

Reed–Solomon codes[edit]

Now that the preliminaries are out of the way, we are ready to begin looking at Reed–Solomon codes.

RS generator polynomial[edit]

Reed–Solomon codes use a generator polynomial similar to BCH codes. The generator is the product of factors (x - αn), starting with n=0 for QR codes. For example:

g4(x) = (x - α0) (x - α1) (x - α2) (x - α3) = 01 x4 + 0f x3 + 36 x2 + 78 x + 40

Here is a function that computes the generator polynomial for a given number of error correction symbols.

def rs_generator_poly(nsym):
   g = [1]
   for i in range(0,nsym):
      g = gf_poly_mul(g, [1, gf_exp[i]])
   return g

This function is somewhat inefficient in that it allocates successively larger arrays for g. While this is unlikely to be a performance problem in practice, readers who are inveterate optimizers may find it interesting to rewrite it so that g is only allocated once.

RS encoding[edit]

Like BCH codes, Reed–Solomon codes are encoded with a process analogous to long division. This example shows the calculation for the message 12 34 56.

                             12 da df
01 0f 36 78 40 ) 12 34 56 00 00 00 00
               ^ 12 ee 2b 23 f4
                    da 7d 23 f4 00
                  ^ da a2 85 79 84
                       df a6 8d 84 00
                     ^ df 91 6b fc d9
                          37 e6 78 d9

The remainder is concatenated with the message, so the encoded message is 12 34 56 37 e6 78 d9. Here is a function that implements this encoding procedure.

def rs_encode_msg(msg_in, nsym):
   gen = rs_generator_poly(nsym)
   msg_out = [0] * (len(msg_in) + nsym)
   for i in range(0, len(msg_in)):
      msg_out[i] = msg_in[i]
   for i in range(0, len(msg_in)):
      coef = msg_out[i]
      if coef != 0:
         for j in range(0, len(gen)):
            msg_out[i+j] ^= gf_mul(gen[j], coef)
   for i in range(0, len(msg_in)):
      msg_out[i] = msg_in[i]
   return msg_out

This example shows the encode function applied to the message in the sample QR code introduced earlier. The calculated error correction symbols (on the second line) match the values decoded from the QR code.

>>> msg_in = [ 0x40, 0xd2, 0x75, 0x47, 0x76, 0x17, 0x32, 0x06,
...            0x27, 0x26, 0x96, 0xc6, 0xc6, 0x96, 0x70, 0xec ]
>>> msg = rs_encode_msg(msg_in, 10)
>>> for i in range(0,len(msg)):
...    print(hex(msg[i]), end=' ')
0x40 0xd2 0x75 0x47 0x76 0x17 0x32 0x6 0x27 0x26 0x96 0xc6 0xc6 0x96 0x70 0xec
0xbc 0x2a 0x90 0x13 0x6b 0xaf 0xef 0xfd 0x4b 0xe0

Python version note: The syntax for the print function has changed, and this example uses the Python 3.0+ version. In previous versions of Python, replace the print line with print hex(msg[i]), (including the final comma).

Syndrome calculation[edit]

Decoding a Reed–Solomon message involves several steps. The first step is to calculate the "syndrome" of the message. Treat the message as a polynomial and evaluate it at α0, α1, α2, ..., αn. Since these are the zeros of the generator polynomial, the result should be zero if the scanned message is undamaged. If not, the syndromes contain all the information necessary to determine the correction that should be made. It is simple to write a function to calculate the syndromes.

def rs_calc_syndromes(msg, nsym):
   synd = [0] * nsym
   for i in range(0, nsym):
      synd[i] = gf_poly_eval(msg, gf_exp[i])
   return synd

Continuing the example, we see that the syndromes of the message are indeed zero. Introducing a change into the message gives nonzero syndromes.

>>> synd = rs_calc_syndromes(msg, 10)
>>> print(synd)
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
>>> msg[0] = 0  # deliberately damage the message
>>> synd = rs_calc_syndromes(msg, 10)
>>> print(synd)
[64, 192, 93, 231, 52, 92, 228, 49, 83, 245]

Erasure correction[edit]

It is simplest to correct mistakes in the code if the locations of the mistakes are already known. This is known as erasure correction. It is possible to correct one erased symbol for each error-correction symbol added to the code. If the error locations are not known, two EC symbols are needed for each symbol error. This makes erasure correction useful in practice if part of the QR code being scanned is covered or physically torn away. It may be difficult for a scanner to determine that this has happened, though, so not all QR code scanners can perform erasure correction.

The Forney algorithm is used to calculate the correction values. It is implemented in the function below.

def rs_correct_errata(msg, synd, pos):
   # calculate error locator polynomial
   q = [1]
   for i in range(0, len(pos)):
      x = gf_exp[len(msg)-1-pos[i]]
      q = gf_poly_mul(q, [x,1])
   # calculate error evaluator polynomial
   p = synd[0:len(pos)]
   p = gf_poly_mul(p, q)
   p = p[len(p)-len(pos):len(p)]
   # formal derivative of error locator eliminates even terms
   qprime = q[len(q)&1:len(q):2]
   # compute corrections
   for i in range(0, len(pos)):
      x = gf_exp[pos[i]+256-len(msg)]
      y = gf_poly_eval(p, x)
      z = gf_poly_eval(qprime, gf_mul(x,x))
      msg[pos[i]] ^= gf_div(y, gf_mul(x,z))

Math note: The denominator of the expression for the error value is the formal derivative of the error locator polynomial q. This is calculated by the usual procedure of replacing each term cn xn with n cn xn-1. Since we're working in a field of characteristic two, n cn is equal to cn when n is odd, and 0 when n is even. Thus, we can simply remove the even coefficients (resulting in the polynomial qprime) and evaluate qprime(x2).

Python note: This function uses array slicing to extract parts of an array. The expression synd[0:len(pos)] returns the first few elements of synd, while p[len(p)-len(pos):len(p)] returns the last few elements of p. The more complicated expression q[len(q)&1:len(q):2] returns every second element of q, skipping the first element if the length of q is odd.

Continuing the example, here we use rs_correct_errata to restore the first byte of the message.

>>> msg[0] = 0
>>> synd = rs_calc_syndromes(msg, 10)
>>> rs_correct_errata(msg, synd, [0])
>>> print(hex(msg[0]))

Error correction[edit]

In the more likely situation where the error locations are unknown, the first step is to find them. The Berlekamp–Massey algorithm is used to calculate the error locator polynomial. Then we simply use trial substitution to find the zeros of this polynomial, which identifies the error locations.

def rs_find_errors(synd, nmess):
   # find error locator polynomial with Berlekamp-Massey algorithm
   err_poly = [1]
   old_poly = [1]
   for i in range(0,len(synd)):
      delta = synd[i]
      for j in range(1,len(err_poly)):
         delta ^= gf_mul(err_poly[len(err_poly)-1-j], synd[i-j])
      if delta != 0:
         if len(old_poly) > len(err_poly):
            new_poly = gf_poly_scale(old_poly, delta)
            old_poly = gf_poly_scale(err_poly, gf_div(1,delta))
            err_poly = new_poly
         err_poly = gf_poly_add(err_poly, gf_poly_scale(old_poly, delta))
   errs = len(err_poly)-1
   if errs*2 > len(synd):
      return None    # too many errors to correct
   # find zeros of error polynomial
   err_pos = []
   for i in range(0, nmess):
      if gf_poly_eval(err_poly, gf_exp[255-i]) == 0:
   if len(err_pos) != errs:
      return None    # couldn't find error locations
   return err_pos

Math note: When the error locator polynomial is linear (err_poly has length 2), it can be solved easily without resorting to a brute-force approach. The function presented above does not take advantage of this fact, but the interested reader may wish to implement the more efficient solution. Similarly, when the error locator is quadratic, it can be solved by using a generalization of the quadratic formula. A more ambitious reader may wish to implement this procedure as well.

Here is an example where three errors in the message are corrected.

>>> print(hex(msg[10]))
>>> msg[0] = 6
>>> msg[10] = 7
>>> msg[20] = 8
>>> synd = rs_calc_syndromes(msg, 10)
>>> pos = rs_find_errors(synd, len(msg))
>>> print(pos)
[20, 10, 0]
>>> rs_correct_errata(msg, synd, pos)
>>> print(hex(msg[10]))

Error and erasure correction[edit]

In order to correct both errors and erasures, we must prevent the erasures from interfering with the error location process. This can be done by calculating the Forney syndromes, as follows.

def rs_forney_syndromes(synd, pos, nmess):
   fsynd = list(synd)      # make a copy
   for i in range(0, len(pos)):
      x = gf_exp[nmess-1-pos[i]]
      for j in range(0,len(fsynd)-1):
         fsynd[j] = gf_mul(fsynd[j], x) ^ fsynd[j+1]
   return fsynd

The Forney syndromes can then be used in place of the regular syndromes in the error location process.

The function rs_correct_msg below brings the complete procedure together. Erasures are indicated by negative values in msg_in.

def rs_correct_msg(msg_in, nsym):
   msg_out = list(msg_in)     # copy of message
   # find erasures
   erase_pos = []
   for i in range(0, len(msg_out)):
      if msg_out[i] < 0:
         msg_out[i] = 0
   if len(erase_pos) > nsym:
      return None     # too many erasures to correct
   synd = rs_calc_syndromes(msg_out, nsym)
   if max(synd) == 0:
      return msg_out  # no errors
   fsynd = rs_forney_syndromes(synd, erase_pos, len(msg_out))
   err_pos = rs_find_errors(fsynd, len(msg_out))
   if err_pos == None:
      return None    # error location failed
   rs_correct_errata(msg_out, synd, erase_pos + err_pos)
   synd = rs_calc_syndromes(msg_out, nsym)
   if max(synd) > 0:
      return None     # message is still not right
   return msg_out

Python note: The lists erase_pos and err_pos are concatenated with the + operator.

This is the last piece needed for a fully operational error-and-erasure correcting Reed–Solomon decoder.


The basic principles of Reed–Solomon codes has been presented in this essay. Working Python code for a particular implementation (QR codes) has been included. Many variations and refinements of these ideas are possible, since coding theory is a very rich field of study.