## Berlekamp Welch Algorithm Beispiel Essay

Error correcting codes are a signal processing technique to correct errors. They are nowadays ubiquitous, such as in communications (mobile phone, internet), data storage and archival (hard drives, optical discs CD/DVD/BluRay, archival tapes), warehouse management (barcodes) and advertisement (QR codes).

Reed–Solomon error correction is a specific type of error correction code. It is one of the oldest but it is still widely used, as it is very well defined and several efficient algorithms are now available under the public domain.

Usually, error correction codes are hidden and most users do not even know about them, nor when they are used. Yet, they are a critical component for some applications to be viable, such as communication or data storage. Indeed, a hard drive that would randomly lose data every few days would be useless, and a phone being able to call only on days with a cloud-less weather would be seldom used. Using error correction codes allows to recover a corrupted message into the full original message.

Barcodes and QR codes are interesting applications to study, as they have the specificity of displaying visually the error correction code, rendering these codes readily accessible to the curious user.

In this essay, we will attempt to introduce the principles of Reed–Solomon codes from the point of view of a programmer rather than a mathematician, which means that we will focus more on the practice than the theory, although we will also explain the theory, but only the necessary knowledge for intuition and implementation. Notable references in the domain will be provided, so that the interested reader can dig deeper into the mathematical theory at will. We will provide real-world examples taken from the popular QR code barcode system as well as working code samples. We chose to use Python for the samples (mainly because it looks pretty and similar to pseudocode), but we will try to explain any non-obvious features for those who are not familiar with it. The mathematics 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.

We will first gently introduce the intuitions behind error correction codes principles, then in a second section we will introduce the structural design of QR codes, in other words how information is stored in a QR code and how to read and produce it, and in a third section we will study error correction codes via the implementation of a Reed–Solomon decoder, with a quick introduction of the bigger BCH codes family, in order to reliably read damaged QR codes.

Note for the curious readers that extended information can be found in the appendix and on the discussion page.

## Principles of error correction

Before detailing the code, it might be useful to understand the intuition behind error correction. Indeed, although error correcting codes may seem daunting mathematically-wise, most of the mathematical operations are high school grade (with the exception of Galois Fields, but which are in fact easy and common for any programmer: it's simply doing operations on integers modulo a number). However, the complexity of the mathematical ingenuity behind error correction codes hide the quite intuitive goal and mechanisms at play.

Error correcting codes might seem like a difficult mathematical concept, but they are in fact based on an intuitive idea with an ingenious mathematical implementation: let's make the data structured, in a way that we can "guess" what the data was if it gets corrupted, just by "fixing" the structure. Mathematically-wise, we use polynomials from the Galois Field to implement this structure.

Let's take a more practical analogy: let's say you want to communicate messages to someone else, but these messages can get corrupted along the way. The main insight of error correcting codes is that, instead of using a whole dictionary of words, we can use a smaller set of carefully selected words, a "reduced dictionary", so that each word is as different as any other. This way, when we get a message, we just have to lookup inside our reduced dictionary to 1) detect which words are corrupted (as they are not in our reduced dictionary); 2) correct corrupted words by finding the most similar word in our dictionary.

Let's take a simple example: we have a reduced dictionary with only three words of 4 letters: , and . Let's say we receive a corrupted word: , where is an erasure. Since we have only 3 words in our dictionary, we can easily compare our received word with our dictionary to find the word that is the closest. In this case, it's . Thus the missing letters are .

Now let's say we receive the word . Here the problem is that we have two words in our dictionary that match the received word: and . In this case, we cannot be sure which one it is, and thus we cannot decode. This means that our dictionary is not very good, and we should replace with another more different word, such as to maximize the difference between each word. This difference, or more precisely the minimum number of different letters between any 2 words of our dictionary, is called the maximum Hamming distance of our dictionary. Making sure that any 2 words of the dictionary share only a minimum number of letters at the same position is called maximum separability.

The same principle is used for most error correcting codes: we generate only a reduced dictionary containing only words with maximum separability (we will detail more how to do that in the third section), and then we communicate only with the words of this reduced dictionary. What Galois Fields provide is the structure (ie, reduced dictionary basis), and Reed–Solomon is a way to automatically create a suitable structure (make a reduced dictionary with maximum separability tailored for a dataset), as well as provide the automated methods to detect and correct errors (ie, lookups in the reduced dictionary). To be more precise, Galois Fields are the structure (thanks to their cyclic nature, the modulo an integer) and Reed–Solomon is the codec (encoder/decoder) based on Galois Fields.

If a word gets corrupted in the communication, that's no big deal since we can easily fix it by looking inside our dictionary and find the closest word, which is probably the correct one (there is however a chance of choosing a wrong one if the input message is too heavily corrupted, but the probability is very small). Also, the longer our words are, the more separable they are, since more characters can be corrupted without any impact.

The simplest way to generate a dictionary of maximally separable words is to make words longer than they really are.

Let's take again our example:

t h i s t h a t c o r n

Append a unique set of characters so that there are no duplicated characters at any of the appended positions, and add one more word to help with the explanation:

t h i s a b c d t h a t b c d e c o r n c d e f t h i n d e f g

Note that each word in this dictionary differs from every other word by at least 5 characters, so the distance is 5. This allows up to 4 errors in known positions, which are called erasures or 2 errors in unknown positions to be corrected.

Assume that 4 erasures occur:

t * * * a b * d

Then a search of the dictionary for the 4 non-erased characters can be done to find the only entry that matches those 4 characters, since the distance is 5.

Assume that 2 errors occur as in one of these patterns:

t h i x a x c d x h i x a b c d t h x x a b c d

The issue here is the location of the errors is unknown. There are 28 possible sub-sets of 6 characters out of 8, so a search using each of the 28 subsets of 6 characters is done, and again there will only be a single sub-set that matches 6 characters since the distance is 5 (and assuming that 2 or fewer errors occurred).

With these examples, one can see the advantage of redundancy in recovering lost information: redundant characters help you recover your original data. The previous examples show how a crude error correcting scheme could work. Reed–Solomon's core idea is similar, append redundant data to a message based on Galois Field mathematics. The original error correcting decoder was similar to the error example above, search sub-sets of a received message that correspond to a valid message, and choose the one with the most matches as the corrected message. This isn't practical for larger messages, so mathematical algorithms were developed to perform error correction in a reasonable time.

## QR code structure

This section introduces the structure of QR codes, which is how data is stored in a QR code. The information in this section is deliberately incomplete. Only the most 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.

### Masking

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.

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

### Formatting information

There are two identical copies of the formatting 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 formatting 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 LevelLevel IndicatorError Correction BytesMessage Data Bytes
L01719
M001016
Q111313
H10179

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).

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

### Message data

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.

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

### Decoding

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 NameMode IndicatorLength BitsData Bits
Numeric00011010 bits per 3 digits
Alphanumeric0010911 bits per 2 characters
Byte010088 bits per character
Kanji1000813 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.)

At this point, we know how to decode, or read, a whole QR code. However, in real life conditions, a QR code is rarely whole: usually, it is scanned via a phone's camera, under potentially poor lighting conditions, or on a scratched surface where part of the QR code was ripped, or on a stained surface, etc.

To make our QR code decoder **reliable**, we need to be able to **correct** errors. The next part of this article will describe how to correct errors, by constructing a BCH decoder, and more specifically a Reed–Solomon decoder.

## BCH codes

In this section, we introduce a general class of error correction codes: the BCH codes, the parent family of modern Reed–Solomon codes, and the basic detection and correction mechanisms.

The formatting 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 (modern Reed–Solomon codes are 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 with the BCH code for format information.

### BCH error detection

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.

00011 10100110111 ) 000111101011001 ^ 10100110111 010100110111 ^ 10100110111 00000000000

Here is a Python function which implements this calculation.

defqr_check_format(fmt):g=0x537# = 0b10100110111 in python 2.6+foriinrange(4,-1,-1):iffmt&(1<<(i+10)):fmt^=g<<ireturnfmt

Python note: The 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 ; in Pascal-derived languages, .

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

If the remainder produced by 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 (ie, lookup in our reduced dictionary).

### BCH error correction

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 (the number of different bits is known as the Hamming distance). This method of finding the closest code is known as exhaustive search, and is possible only because we have very few codes (a code is a valid message, and here there are only 32, all other binary numbers aren't correct).

(Note that Reed–Solomon is also based on this principle, but since the number of possible codewords is simply too big, we can't afford to do an exhaustive search, and that's why clever but complicated algorithms have been devised, such as Berlekamp-Massey.)

defhamming_weight(x):weight=0whilex>0:weight+=x&1x>>=1returnweightdefqr_decode_format(fmt):best_fmt=-1best_dist=15fortest_fmtinrange(0,32):test_code=(test_fmt<<10)^qr_check_format(test_fmt<<10)test_dist=hamming_weight(fmt^test_code)iftest_dist<best_dist:best_dist=test_distbest_fmt=test_fmteliftest_dist==best_dist:best_fmt=-1returnbest_fmt

The function 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 , , and 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 3 >>> qr_decode_format(int("111111101011001",2)) # 3 bit-errors 3 >>> qr_decode_format(int("111011101011001",2)) # 4 bit-errors -1

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

In the next sections, we will study Finite Field Arithmetics and Reed–Solomon code, which is a subtype of BCH codes. The basic idea (ie, using a limited words dictionary with maximum separability) is the same, but since we will encode longer words (256 bytes instead of 2 bytes), with more symbols available (encoded on all 8bits, thus 256 different possible values), we cannot use this naive, exhaustive approach, because it would take way too much time: we need to use cleverer algorithms, and Finite Field mathematics will help us do just that, by giving us a structure.

## Finite field arithmetic

### Introduction to mathematical fields

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

We'd like to define addition, subtraction, multiplication, and division for 8-bit bytes and always produce 8-bit bytes as a result, so as to avoid any overflow. Naively, we might attempt to use the normal definitions for these operations, and then mod by 256 to keep results from overflowing. And this is exactly what we will be doing, and is what is called a Galois Field 2^8. You can easily imagine why it works for everything, except for division: what is 7/5 ?

Here's a brief introduction to Galois Fields: a finite field is a set of numbers, and a field needs to have six properties: Closure, Associative, Commutative, Distributive, Identity and Inverse. More simply put, using a field allow to study the relationship between numbers of this field, and apply the result to any other field that follows the same properties. For example, the set of reals ℝ is a field. In other words, mathematical fields studies the structure of a set of numbers.

However, integers ℤ aren't a field, because as we said above, not all divisions are defined (such as 7/5), which violates multiplicative inverse property (x such as 7*x=5 does not exist). One simple way to fix that is to use modulo using a prime number, such as 2: in this way, we are guaranteed that 7*x=5 exists since we will just wrap around. ℤ modulo 2 is called a Galois Field, and any number divisible by 2 is a Galois Field (since we need to modulo using a prime number), thus 256, the value of an 8-bit symbol, can be reduced to 2^8, and thus we say that we use a Galois Field of 2^8, or GF(2^8). In spoken language, 2 is the characteristic of the field, 8 is the exponent, and 256 is the field's cardinality. More information on finite fields can be found here.

Here we will define the usual mathematical operations that you are used to doing with integers, but adapted to GF(2^8), which is basically doing usual operations but modulo 2^8.

Another way to consider the link between GF(2) and GF(2^8) is to think that GF(2^8) represents a polynomial of 8 binary coefficients. For example, in GF(2^8), 170 is equivalent to . Both representations are equivalent, it's just that in the first case, 170, the representation is decimal, and in the other case it's binary, which can be thought as representing a polynomial by convention (only used in GF(2^p) as explained here). The latter is often the representation used in academic books and in hardware implementations (because of logical gates and registers, which work at the binary level). For a software implementation, the decimal representation can be preferred for clearer and more close-to-the-mathematics code (this is what we will use for the code in this tutorial, except for some examples that will use the binary representation).

In any case, try to not confuse the polynomial representing a single GF(2^p) symbol (each coefficient is a bit/boolean: either 0 or 1), and the polynomial representing a list of GF(2^p) symbols (in this case the polynomial is equivalent to the message+RScode, each coefficient is a value between 0 and 2^p and represent one character of the message+RScode). We will first describe operations on single symbol, then polynomial operations on a list of symbols.

### Addition and Subtraction

Both addition and subtraction are replaced with exclusive-or in Galois Field base 2. This is logical: addition modulo 2 is exactly like an XOR, and subtraction modulo 2 is exactly the same as addition modulo 2. This is possible because additions and subtractions in this Galois Field are carry-less.

Thinking of our 8-bit values as polynomials with coefficients mod 2:

0101 + 0110 = 0101 - 0110 = 0101 XOR 0110 = 0011

The same way (in binary representation of two single GF(2^8) integers):

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

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

Note that in books, you will find additions and subtractions to define some mathematical operations on GF integers, but in practice, you can just XOR (as long as you are in a Galois Field, this is not true in other fields).

Here is the equivalent Python code:

defgf_add(x,y):returnx^ydefgf_sub(x,y):returnx^y# in binary galois field, subtraction is just the same as addition (since we mod 2)

### Multiplication

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.

10001001 * 00101010 10001001 ^ 10001001 ^ 10001001 1010001111010

Note: the XOR multiplication here is carry-less! If you do it with-carry, you will get the wrong result 1011001111010 with the extra term x9 instead of the correct result 1010001111010.

Here is a Python function which implements this polynomial multiplication on single GF(2^8) integers.

Note: this function (and some other functions below) use a lot of bitwise operators such as >> and <<, because they are both faster and more concise to do what we want to do. These operators are available in most languages, they are not specific to Python, and you can get more information about them here.

defcl_mul(x,y):'''Bitwise carry-less multiplication on integers'''z=0i=0while(y>>i)>0:ify&(1<<i):z^=x<<ii+=1returnz

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 (the choice of this number is explained below the code), using the long division process described previously. In this instance, this is called "modular reduction", because basically what we do is that we divide and keep only the remainder, using a modulo. This produces the final answer 11000011 in our example.

1010001111010 ^ 100011101 0010110101010 ^ 100011101 00111011110 ^ 100011101 011000011

Here is the Python code to do the whole Galois Field multiplication with modular reduction:

defgf_mult_noLUT(x,y,prim=0):'''Multiplication in Galois Fields without using a precomputed look-up table (and thus it's slower) by using the standard carry-less multiplication + modular reduction using an irreducible prime polynomial'''### Define bitwise carry-less operations as inner functions ###defcl_mult(x,y):'''Bitwise carry-less multiplication on integers'''z=0i=0while(y>>i)>0:ify&(1<<i):z^=x<<ii+=1returnzdefbit_length(n):'''Compute the position of the most significant bit (1) of an integer. Equivalent to int.bit_length()'''bits=0whilen>>bits:bits+=1returnbitsdefcl_div(dividend,divisor=None):'''Bitwise carry-less long division on integers and returns the remainder'''# Compute the position of the most significant bit for each integersdl1=bit_length(dividend)dl2=bit_length(divisor)# If the dividend is smaller than the divisor, just exitifdl1<dl2:returndividend# Else, align the most significant 1 of the divisor to the most significant 1 of the dividend (by shifting the divisor)foriinrange(dl1-dl2,-1,-1):# Check that the dividend is divisible (useless for the first iteration but important for the next ones)ifdividend&(1<<i+dl2-1):# If divisible, then shift the divisor to align the most significant bits and XOR (carry-less subtraction)dividend^=divisor<<ireturndividend### Main GF multiplication routine #### Multiply the gf numbersresult=cl_mult(x,y)# Then do a modular reduction (ie, remainder from the division) with an irreducible primitive polynomial so that it stays inside GF boundsifprim>0:result=cl_div(result,prim)returnresult

Result:

>>> a = 0b10001001 >>> b = 0b00101010 >>> print bin(gf_mult_noLUT(a, b, 0)) # multiplication only 0b1010001111010 >>> print bin(gf_mult_noLUT(a, b, 0x11d)) # multiplication + modular reduction 0b11000011

Why mod 100011101 (in hexadecimal: 0x11d)? The mathematics is a little complicated here, but in short, 100011101 represents an 8th degree polynomial which is "irreducible" (meaning it can't be represented as the product of two smaller polynomials). This number is called a primitive polynomial or irreducible polynomial, or prime polynomial (we will mainly use this latter name for the rest of this tutorial). This is necessary for division to be well-behaved, which is to stay in the limits of the Galois Field, but without duplicating values. There are other numbers we could have chosen, but they're all essentially the same, and 100011101 (0x11d) is a common primitive polynomial for Reed–Solomon codes. If you are curious to know how to generate those prime numbers, please see the appendix.

Additional infos on the prime polynomial (you can skip): What is a prime polynomial? It is the equivalent of a prime number, but in the Galois Field. Remember that a Galois Field uses values that are multiples of 2 as the generator. Of course, a prime number cannot be a multiple of two in standard arithmetics, but in a Galois Field it is possible. Why do we need a prime polynomial? Because to stay in the bound of the field, we need to compute the modulo of any value above the Galois Field. Why don't we just modulo with the Galois Field size? Because we will end up with lots of duplicate values, and we want to have as many unique values as possible in the field, so that a number always has one and only projection when doing a modulo or a XOR with the prime polynomial.

Note for the interested reader: as an example of what you can achieve with clever algorithms, here is another way to achieve multiplication of GF numbers in a more concise and faster way, using the Russian Peasant Multiplication algorithm:

defgf_mult_noLUT(x,y,prim=0,field_charac_full=256,carryless=True):'''Galois Field integer multiplication using Russian Peasant Multiplication algorithm (faster than the standard multiplication + modular reduction). If prim is 0 and carryless=False, then the function produces the result for a standard integers multiplication (no carry-less arithmetics nor modular reduction).'''r=0whiley:# while y is above 0ify&1:r=r^xifcarrylesselser+x# y is odd, then add the corresponding x to r (the sum of all x's corresponding to odd y's will give the final product). Note that since we're in GF(2), the addition is in fact an XOR (very important because in GF(2) the multiplication and additions are carry-less, thus it changes the result!).y=y>>1# equivalent to y // 2x=x<<1# equivalent to x*2ifprim>0andx&field_charac_full:x=x^prim# GF modulo: if x >= 256 then apply modular reduction using the primitive polynomial (we just subtract, but since the primitive number can be above 256 then we directly XOR).returnr

Note that using this last function with parameters prim=0 and carryless=False will return the result for a standard integers multiplication (and thus you can see the difference between carryless and with-carry addition and its impact on multiplication).

### Multiplication with logarithms

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 2=00000010 (by convention, this number is referred to as α or the generator number): 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 α, that we define, 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. More mathematical information about this trick can be found here.

gf_exp=[0]*512# Create list of 512 elements. In Python 2.6+, consider using bytearraygf_log=[0]*256definit_tables(prim=0x11d):'''Precompute the logarithm and anti-log tables for faster computation later, using the provided primitive polynomial.'''# prim is the primitive (binary) polynomial. Since it's a polynomial in the binary sense,# it's only in fact a single galois field value between 0 and 255, and not a list of gf values.globalgf_exp,gf_loggf_exp=[0]*512# anti-log (exponential) tablegf_log=[0]*256# log table# For each possible value in the galois field 2^8, we will pre-compute the logarithm and anti-logarithm (exponential) of this valuex=1foriinrange(0,255):gf_exp[i]=x# compute anti-log for this value and store it in a tablegf_log[x]=i# compute log at the same timex=gf_mult_noLUT(x,2,prim)# If you use only generator==2 or a power of 2, you can use the following which is faster than gf_mult_noLUT():#x <<= 1 # multiply by 2 (change 1 by another number y to multiply by a power of 2^y)#if x & 0x100: # similar to x >= 256, but a lot faster (because 0x100 == 256)#x ^= prim # substract the primary polynomial to the current value (instead of 255, so that we get a unique set made of coprime numbers), this is the core of the tables generation# Optimization: double the size of the anti-log table so that we don't need to mod 255 to# stay inside the bounds (because we will mainly use this table for the multiplication of two GF numbers, no more).foriinrange(255,512):gf_exp[i]=gf_exp[i-255]return[gf_log,gf_exp]

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

The table is oversized in order to simplify the multiplication function. This way, we don't have to check to make sure that is within the table size.

defgf_mul(x,y):ifx==0ory==0:return0returngf_exp[gf_log[x]+gf_log[y]]# should be gf_exp[(gf_log[x]+gf_log[y])%255] if gf_exp wasn't oversized

### Division

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.

defgf_div(x,y):ify==0:raiseZeroDivisionError()ifx==0:return0returngf_exp[(gf_log[x]+255-gf_log[y])%255]

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

With this definition of division, for any and any nonzero .

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 and 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.

### Power and Inverse

The logarithm table approach will once again simplify and speed up our calculations when computing the power and the inverse:

defgf_pow(x,power):returngf_exp[(gf_log[x]*power)%255]defgf_inverse(x):returngf_exp[255-gf_log[x]]# gf_inverse(x) == gf_div(1, x)

### Polynomials

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 = x4 + x3 + x2 + x +

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

The first function multiplies a polynomial by a scalar.

defgf_poly_scale(p,x):r=[0]*len(p)foriinrange(0,len(p)):r[i]=gf_mul(p[i],x)returnr

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).

defgf_poly_add(p,q):r=[0]*max(len(p),len(q))foriinrange(0,len(p)):r[i+len(r)-len(p)]=p[i]foriinrange(0,len(q)):r[i+len(r)-len(q)]^=q[i]returnr

The next function multiplies two polynomials.

defgf_poly_mul(p,q):'''Multiply two polynomials, inside Galois Field'''# Pre-allocate the result arrayr=[0]*(len(p)+len(q)-1)# Compute the polynomial multiplication (just like the outer product of two vectors,# we multiply each coefficients of p with all coefficients of q)forjinrange(0,len(q)):foriinrange(0,len(p)):r[i+j]^=gf_mul(p[i],q[j])# equivalent to: r[i + j] = gf_add(r[i+j], gf_mul(p[i], q[j]))# -- you can see it's your usual polynomial multiplicationreturnr

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. Horner's method works by factorizing consecutively the terms, so that we always deal with x^1, iteratively, avoiding the computation of higher degree terms:

x4 + x3 + x2 + x + = (((x + ) x + ) x + ) x +
defgf_poly_eval(poly,x):'''Evaluates a polynomial in GF(2^p) given the value for x. This is based on Horner's scheme for maximum efficiency.'''y=poly[0]foriinrange(1,len(poly)):y=gf_mul(y,x)^poly[i]returny

There's still one missing polynomial operation that we will need: polynomial division. This is more complicated than the other operations on polynomial, so we will study it in the next chapter, along with Reed–Solomon encoding.

## Reed–Solomon codes

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

### Insight of the coding theory

But first, why did we have to learn about finite fields and polynomials? Because this is the main insight of error-correcting codes like Reed–Solomon: instead of just seeing a message as a series of (ASCII) numbers, we see it as a polynomial following the very well-defined rules of finite field arithmetic. In other words, by representing the data using polynomials and finite fields arithmetic, we added a structure to the data. The values of the message are still the same, but this conceptual structure now allow us to operate on the message, on its characters values, using well defined mathematical rules. This structure, that we always know because it's outside and independent of the data, is what allows us to repair a corrupted message.

Thus, even if in your code implementation you may choose to not explicitly represent the polynomials and the finite field arithmetic, these notions are essential for the error-correcting codes to work, and you will find these notions to underlie (even if implicitly) any implementation.

And now we will put these notions into practice!

### RS encoding

#### Encoding outline

Like BCH codes, Reed–Solomon codes are encoded by dividing the polynomial representing the message by an irreducible generator polynomial, and then the remainder is the RS code, which we will just append to the original message.

Why? We previously said that the principle behind BCH codes, and most other error correcting codes, is to use a reduced dictionary with very different words as to maximize the distance between words, and that longer words have greater distance: here it's the same principle, first because we lengthen the original message with additional symbols (the remainder) which raise the distance, and secondly because the remainder is almost unique (thanks to the carefully designed irreducible generator polynomial), so that it can be exploited by clever algorithms to deduce parts of the original message.

To summary, with an approximated analogy to encryption: our generator polynomial is our encoding dictionary, and polynomial division is the operator to convert our message using the dictionary (the generator polynomial) into a RS code.

#### Exception management

To manage errors and cases where we can't correct a message, we will display a meaningful error message, by raising an exception. We will make our own custom exception so that users can easily catch and manage them:

classReedSolomonError(Exception):pass

To display an error by raising our custom exception, we can then simply do the following:

raiseReedSolomonError("Some error message")

And you can easily catch this exception to manage it by using a try/except block:

try:raiseReedSolomonError("Some error message")exceptReedSolomonError,e:pass# do something here

#### RS generator polynomial

Reed–Solomon codes use a generator polynomial similar to BCH codes (not to be confused with the generator number alpha). 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) = x4 + x3 + x2 + x +

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

defrs_generator_poly(nsym):'''Generate an irreducible generator polynomial (necessary to encode a message into Reed-Solomon)'''g=[1]foriinrange(0,nsym):g=gf_poly_mul(g,[1,gf_pow(2,i)])returng

This function is somewhat inefficient in that it allocates successively larger arrays for . 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 is only allocated once, or you can compute once and memorize g since it is fixed for a given nsym, so you can reuse g.

#### Polynomial division

Several algorithms for polynomial division exist, the simplest one that is often taught in high school is long division. This example shows the calculation for the message .

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

Note: This may seem complicated at first, but it works exactly the usual way when you do polynomial long division with integers in base 10, it's just that here it's in base 16 (hexadecimal). If you have some trouble understanding this example, you can try to convert the numbers into base 10, and you will probably quickly see that the result is natural.

The remainder is concatenated with the message, so the encoded message is .

However, long division is quite slow as it requires a lot of recursive iterations to terminate. More efficient strategies can be devised, such as using synthetic division (also called Horner's method, a good tutorial video can be found on Khan Academy). Here is a function that implements extended synthetic division of GF(2^p) polynomials (extended because the divisor is a polynomial instead of a monomial):

defgf_poly_div(dividend,divisor):'''Fast polynomial division by using Extended Synthetic Division and optimized for GF(2^p) computations (doesn't work with standard polynomials outside of this galois field, see the Wikipedia article for generic algorithm).'''# CAUTION: this function expects polynomials to follow the opposite convention at decoding:# the terms must go from the biggest to lowest degree (while most other functions here expect# a list from lowest to biggest degree). eg: 1 + 2x + 5x^2 = [5, 2, 1], NOT [1, 2, 5]msg_out=list(dividend)# Copy the dividend#normalizer = divisor[0] # precomputing for performanceforiinrange(0,len(dividend)-(len(divisor)-1)):#msg_out[i] /= normalizer # for general polynomial division (when polynomials are non-monic), the usual way of using# synthetic division is to divide the divisor g(x) with its leading coefficient, but not needed here.coef=msg_out[i]# precachingifcoef!=0:# log(0) is undefined, so we need to avoid that case explicitly (and it's also a good optimization).forjinrange(1,len(divisor)):# in synthetic division, we always skip the first coefficient of the divisior,# because it's only used to normalize the dividend coefficientifdivisor[j]!=0:# log(0) is undefinedmsg_out[i+j]^=gf_mul(divisor[j],coef)# equivalent to the more mathematically correct# (but xoring directly is faster): msg_out[i + j] += -divisor[j] * coef# The resulting msg_out contains both the quotient and the remainder, the remainder being the size of the divisor# (the remainder has necessarily the same degree as the divisor -- not length but degree == length-1 -- since it's# what we couldn't divide from the dividend), so we compute the index where this separation is, and return the quotient and remainder.separator=-(len(divisor)-1)returnmsg_out[:separator],msg_out[separator:]# return quotient, remainder.

#### Encoding main function

And now, here's how to encode a message to get its RS code:

defrs_encode_msg(msg_in,nsym):'''Reed-Solomon main encoding function'''gen=rs_generator_poly(nsym)# Pad the message, then divide it by the irreducible generator polynomial_,remainder=gf_poly_div(msg_in+[0]

The Berlekamp–Welch algorithm, also known as the Welch–Berlekamp algorithm, is named for Elwyn R. Berlekamp and Lloyd R. Welch. This is a decoder algorithm that efficiently corrects errors in Reed–Solomon codes for an RS(n, k), code based on the Reed Solomon original view where a message is used as coefficients of a polynomial or used with Lagrange interpolation to generate the polynomial of degree < k for inputs and then is applied to to create an encoded codeword .

The goal of the decoder is to recover the original encoding polynomial , using the known inputs and received codeword with possible errors. It also computes an error polynomial where corresponding to errors in the received codeword.

## The key equations

Defining e = number of errors, the key set of n equations is

Where E(ai) = 0 for the e cases when bi ≠ F(ai), and E(ai) ≠ 0 for the n - e non error cases where bi = F(ai) . These equations can't be solved directly, but by defining Q() as the product of E() and F():

and adding the constraint that the most significant coefficient of E(ai) = ee = 1, the result will lead to a set of equations that can be solved with linear algebra.

where q = n - e - 1. Since ee is constrained to be 1, the equations become:

resulting in a set of equations which can be solved using linear algebra, with time complexity O(n^3).

The algorithm begins assuming the maximum number of errors e = ⌊ (n-k)/2 ⌋. If the equations can not be solved (due to redundancy), e is reduced by 1 and the process repeated, until the equations can be solved or e is reduced to 0, indicating no errors. If Q()/E() has remainder = 0, then F() = Q()/E() and the code word values F(ai) are calculated for the locations where E(ai) = 0 to recover the original code word. If the remainder ≠ 0, then an uncorrectable error has been detected.

## Example

Consider RS(7,3) (n = 7, k = 3) defined in GF(7) with α = 3 and input values: ai = i-1 : {0,1,2,3,4,5,6}. The message to be systematically encoded is {1,6,3}. Using Lagrange interpolation, F(ai) = 3 x2 + 2 x + 1, and applying F(ai) for a4 = 3 to a7 = 6, results in the code word {1,6,3,6,1,2,2}. Assume errors occur at c2 and c5 resulting in the received code word {1,5,3,6,3,2,2}. Start off with e = 2 and solve the linear equations:

Starting from the bottom of the right matrix, and the constraint e2 = 1:

with remainder = 0.

E(ai) = 0 at a2 = 1 and a5 = 4 Calculate F(a2 = 1) = 6 and F(a5 = 4) = 1 to produce corrected code word {1,6,3,6,1,2,2}.

## External links

• MIT Lecture Notes on Essential Coding Theory – Dr. Madhu Sudan
• University at Buffalo Lecture Notes on Coding Theory – Dr. Atri Rudra
• Algebraic Codes on Lines, Planes and Curves, An Engineering Approach – Richard E. Blahut
• Welch Berlekamp Decoding of Reed–Solomon Codes – L. R. Welch
• US 4,633,470, Welch, Lloyd R. & Elwyn R. Berlekamp, "Error Correction for Algebraic Block Codes", published September 27, 1983, issued December 30, 1986  – The patent by Lloyd R. Welch and Elewyn R. Berlekamp