Subject: Re: Reed Solomon code needed

From: [email-masked] (Risto Rajala NTC/SWS/CPC Kalkki2 puh 90-511 5843)
Date: 28 Nov 94 15:44:02 EET

In article <[email-masked]>,
Andrew Powell <[email-masked]> writes:

> I am looking for source code for the Reed Solomon ECC algorithm. 

Here is R-S-Code articles collected and posted by Tom Moore two years ago.

        Risto Rajala, [email-masked]

--------------------------------------------------------------------------------

Subject: REED-SOLOMON PACKAGE
From: [email-masked] (Tom Moore)

The collection of Reed-Solomon goodies.

The following collection should help you with the theory and
actual code behind REED-SOLOMON Codes.

The included files are responses sent back to me by request.
If you use anything, please give credit to the proper
authors.  As of yet, I have NOT written anything to do
with RS (But I will shortly........)

Tom Moore  			{...uunet!isc-br!ewu!tmoore}
Eastern Washington University   [email-masked]}
Head Consultant, Computer Graphics Lab
(509) 623-4235

---------------------------------------------------------------------------
# This is a shell archive.  Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by cg5!tmoore on Fri Nov 15 14:24:06 PST 1991
# Contents:  rs.galois rs.rply.1 rs.rply.2 rs.rply.3
 
echo x - rs.galois
sed 's/^@//' > "rs.galois" <<'@//E*O*F rs.galois//'
@From: [email-masked] (Colin Plumb)
Subject: REED/SOLOMON -Galois fields
Date: 13 Nov 91 02:31:18 GMT
Organization: Array Systems Computing, Inc., Toronto, Ontario, CANADA

[email-masked] writes:

>There has got to be someone out there who has messed with REED/SOLOMON
>error correcting codes.  I know, I know, you all are interested in
>getting avery little bit of redundancy out of text, and here I come
>wanting to put some back in.  It may not sit will with you, but hey,
>I need some help here, so how about a little mail to me.

Okay, time for the old tutorial...

First, Galois fields.  In number theory, we have a thing called a
"group".  This is a bunch of things I'll call numbers, and
an operation, usually written as addition or multiplication.
The axioms of a group are, roughly:

For each a, and b in the group, a+b is in the group
There is an identity element, often written e, such that a+e = e + a = a.
For each a in the group, there is an inverse -a so a + -a = 0.
For all a, b and c, (a+b)+c = a+(b+c)

Some groups are commutative, or Abelian (after the mathematician Abel),
in which case you also have a+b = b+a.  Everything I'll work with here
is, in fact, Abelian, although there are structures that aren't.

An example of a group is the integers, or the integers modulo 10,
with normal addition.  5+5 = 0, so -5 = 5, and so on.

Fancier than that is a ring, which has addition and multiplication.
Multiplication also stays within the ring, and it distributes:
a*(b+c) = (a*b)+(a*c).  There's also a multiplicative identity,
usually written as 1.  a*1 = 1*a = a.
The integers are a ring, as are the integers modulo 10.

Fancier still is a field, which has multiplicative inverses.  For every
a other than 0 (the additive inverse), there is a multiplicative
inverse a', such that a*a' = a'*a = 1.  The integers are not a field,
although the rationals are.  So are the integers modulo any prime
number.  Modulo 10, 5 does not have a multiplicative inverse: times
anything, it's either 5 or 0.  Modulo 9 (also not a prime), 3 does not
have a multiplicative inverse.  But modulo 11, everything has a
multiplicative inverse.  (This property comes up in the prime-testing
algorithms that public-key cryptography has produced so much excitement
in, but I'm drifting...)

Anyway, it turns out that given a ring, you can define primes.
And a ring, modulo a prime, produces a field.  (There are actually
come complications involved in prime vs. irreducible, but I'll ignore
them - this theory is *not* designed to let you pass your exam
tomorrow; it's the bare minimum motivation for what's to follow.)

So, anyway, we now have the material for producing rings of "nice"
sizes.  2 is prime, so we can make a field of size 2, but it's not
big enough to do much with.  And no other prime is a power of 2,
as is suitable to work with on binary computers.  Galois to the rescue!

Given a ring, we can define polynomials over that ring.  You know, like
3x^2 + 2x + 1.  Basically, it's just a list of numbers in the ring
(e.g. 1, 2, 3 for the above - equivalent to 1, 2, 3, 0, 0, 0,...).
There are a few more things you can do if the numbers are from a field,
and it turns out we'll be using a field, so I'll stop making the
distinction.

Now, we can multiply polynomials, and add them, and so on.  Furthermore,
we have an additive identity, additive inverses, a multiplicative
identity, and so on.  It turns out to be a ring.

So, let's take as our basic field the integers modulo 2.  I.e.
0 and 1.  0+0 = 1+1 = 0, 0+1 = 1+0 = 1, 0*0 = 0*1 = 1*0 = 0, 1*1 = 1.

We can now form polynomials over the integers modulo 2, in which
each power of x (like, say, x^12) is either present (coefficient of 1)
or absent (coefficient of 0).  These can be expressed quite nicely
as bit strings.  E.g. 10001000000100001 is the polynomial x^16+x^12+x^5+1.
This polynomial can be factored as 1111000000011111 * 11,
x^15+x^14+x^13+x^12+x^4+x^3+x^2+x+1 * x+1.  Both of these are prime -
they cannot be expressed as the product of smaller polynomials.

Now, observe that adding polynomials involves adding corresponding
coefficients, modulo 2.  This is also known as the xor operation.
Also observe that each polynomial is its own additive inverse, so
subtracting is the same as adding.  And multiplication is just
shift-and-add, like usual, except it's shift-and-xor.

Now, given primes, we can form the polynomials modulo a prime.
Let's take 1011, x^3+x+1, a polynomial which happens to be prime.
We can take any given polynomial (say, 11111, x^4+x^3+x^2+x+1), and
subtract multiples of 1011 until it has degree less than 3.

First, subtract x times 1011, i.e. 10 times 1011, or 10110:

11111
10110
=====
 1001

Then, subtract 1 times 1011:
 1001
 1011
=====
 0010

Thus, x^4+x^3+x^2+x+1 modulo x^3+x+1 is x.

If you think for a bit, you'll realize that exactly the polynomials with
degree less than the modulus are the possible remainders.  Other than
that, any coefficients are possible.  So if the modulus is of degree n
(highest coefficient is x^n), then there are 2^n possible remainders,
corresponding to all possible combinations of x^0, x^1, ..., x^(n-1).
These are written as all possible n-bit strings.

Now, we're getting somewhere!  We have numbers we can work with (it
turns out that multiplicative inverses exist, if the modulus is
irreducible, which is prime plus a little bit that I won't get into),
and we have a nice equivalnece with bit-strings.

So, let's take some degree-8 polynomial and use it as the modulus,
so all possible bytes are our numbers.

Now, consider what would happen if we were to store, after a bunch of
bytes a0, a1, ... a(n-1), their sum s (using this field, which means we
xor them).  The total would be a0+...+a(n-1)+s = a0+...+a(n-1)+a0+...+a(n-1)
= 0.  This is basically parity.

If an error hit one of the bytes, say an error e in position k, so the
received data is b0 = a0, b1 = a1, ..., bk = ak+e, ... b(n-1) = a(n-1),
bn = s.  Then the sum wouldn't add up, and we'd get, instead of 0, e.

Now, this is an error-detection code, but isn't much better than parity.
It is, in fact, 8-way interleaved parity.  But we can get better.

Now, let's come up with a second sum character, t.
t = a0*x^0 + a1*x^1 + a2*x^2 + ... + a(n-1)*x^(n-1).
Now, if an error occurs, then summing t on the received side will
give us, not the received t, but t+e*x^k.  We cah subtract off the
received t to get e*x^k.

Now, we have e and e*x^k.  Since e is not 0 (if it was, there'd be
no error), it has a multiplicative inverse, e'.  We can multiply
this by e*x^k to get x^k.  As long as n is low enough that each
possible k has a different x^k (there are 255 such possibilities for
an 8-bit code), then from this we can recover the position, k, of
the error.  Since we know the error pattern and the location,
we can fix it.  Bingo, error correction!

(If the error is in the user data, then x^k is non-zero, so both
e and e*x^k are non-zero.  If one or the other is zero, then the
error is in the coresponding check byte and can be ignored.  If
both are zero, there was no error.)

This is the simplest single-error-correcting Reed-Solomon code.
It turns out that by adding more check bytes, computed by
u = a0*x^0+a1*x^2+a2*x^4+..., v = a0*x^0+a1*x^3+a2*x^6+...,
etc. and solving more complex equations to find the solutions,
we can determine the patterns and locations of multiple errors.

For example, suppose we have 4 check bytes and 2 errors.  e1 at
location k1 and e2 at location k2.

We can derive, from the check bytes, e1+e2, e1*x^k1+e2*x^k2,
e1*x^(2*e1)+e2*x^(2*k2) and e1*x^(3*e1)+e2*x^(3*k2).  From
these four equations, we can derive the four unknowns e1, k1, e2, and
k2 we need to fix the error.  See various journal articles and
textbooks for efficient ways to solve these systems of equations.  With
r check bytes, you can detect r errors, correct r/2 errors, or any
combination.  (Say, correct up to r/3 errors and detect up to 2r/3.) In
fact, even with more errors, you'll *probably* detect them, but it's
just not absolutely guaranteed.

What CD's do is a bit more complicated... they have two levels of
Reed-Solomon ECC.  First, they scatter the data around on the disc so a
burst error will take a bit out of a bunch of blocks rather than
totally destroying a few blocks.  Then, they arrange the block in a
grid and apply error correction codes to the rows.  If there are just
one or two errors, they repair them, but if there are a bunch, they
mark the whole row as bad, and don't try to correct it.

Next, they apply error correction down the columns.  Because any rows
with errors are already marked bad, we don't have to derive the error
locations k1, k2, k3, etc. and so with the number of check bytes we
have (which determines the number of equations), we can have more
errors e1, e2,... without having too many unknowns.  (It's also easier
to solve the equations for the errors than the locations, so this makes
the implementation easier.)  This results in really spectacular
error-correction ability, over a large block, without having to deal
with hundreds of equations in hundreds of unknowns, which would make CD
players ridiculously expensive.  (Remember, you can't slow down the
music when an error occurs - that's just as bad as no correction at
all.)

That was a *really* quick summary, but the essentials are there.
-- 
	-Colin

@//E*O*F rs.galois//
chmod u=rw,g=,o= rs.galois
 
echo x - rs.rply.1
sed 's/^@//' > "rs.rply.1" <<'@//E*O*F rs.rply.1//'
@From: [email-masked]
To: [email-masked]
Subject: REED/SOLOMON
Cc: [email-masked]

If you have no experience with error correcting codes, then
I think it may be a little difficult to start with RS codes.
But anyway,

RS codes are a general case of BCH codes (i.e., RS codes were
discovered first and then they were generalized to BCH codes).
So, any literature you find about nonbinary BCH codes can be 
applied to RS codes.  A good starter book on coding theory is:
ERROR CONTROL CODING: FUNDAMENTALS AND APPLICATIONS written by
Shu Lin and Daniel J. Costello, Jr.  In particular, Chapter 6,
discusses the encoding and decoding of binary and nonbinary BCH
codes along with RS codes.  It should be very simple to write the
encoder.  The hard part is determining the type of decoder you want
(i.e., do I want a very fast algorithm, am I worried about memory,
how large are the codewords, how much parity do I need, and do I 
need to correct erasures).

If you have any questions, let me know.   

Jim Kosmach
Research Engineer
Motorola Corporate Research Center
[email-masked]

@//E*O*F rs.rply.1//
chmod u=rw,g=,o= rs.rply.1
 
echo x - rs.rply.2
sed 's/^@//' > "rs.rply.2" <<'@//E*O*F rs.rply.2//'
@From: [email-masked] (Rob Warnock)
Subject: Re: REED/SOLOMON

Rob Warnock, MS-1L/515		[email-masked]
Silicon Graphics, Inc.		(415)335-1673
2011 N. Shoreline Blvd.
Mountain View, CA  94039-7311

=========================================================================
============= attachment 1 of 2 =========================================
=========================================================================
@From: [email-masked]
Subject: ECC bibliocraphy 
Organization: Silicon Graphics Inc., Mountain View, CA

Here's my bibliography for ECC:

        - Lin & Costello, "Error Control Coding: Fundamentals and
          Applications" (Prentice-Hall)
        - Peterson & Weldon, "Error Correcting Codes" (MIT Press)
        - Blahut, "Theory and Practice of Error Control Codes" (Addison-Wesley)
        - Berlekamp, "Algebraic Coding Theory" (Aegean Park Press)

But for beginners, good books are:

        - Hamming, "Coding and Information Theory" (Prentice-Hall)
        - Arazi, "A Commonsense Approach to the Theory of Error-Correcting
          Codes" (MIT Press)

Hamming's little book is an easy start, but only goes through single-error
correction. Lin & Costello is a pretty good cross-section, and has some
good stuff about Reed-Solomon codes used in disks & tapes. For example,
on pp.525-531 they discuss the particularly simple single-error correcting
shortened (174,171) R-S code over GF(256) used in the IBM 3370 disk. [You
can also use it in the unshortened (255,252) form.] This R-S code can easily
be done in software, with small table lookups. [I posted some rough, untested
pseudo-C for it in this newsgroup back in July...]

=========================================================================
============= attachment 2 of 2 =========================================
=========================================================================
@From: [email-masked]
Subject: Re: Compression is great, but, robustness?
Organization: Silicon Graphics Inc., Mountain View, CA

[email-masked] (Mark Adler) writes:
+---------------
| Bill Davidsen notes:
| >>   You can, of course, put error detection and correction in the data,
| >> using something simple like Hamming code, or larger and more complex
| >> like Fire codes.
+---------------

Fire codes are totally obsolete these days. Their miscorrection rates are
quite high compared to the available alternatives.

+---------------
| Better still are Reed-Solomon codes.  The algebra is a little complicated,
| but the run-time of a good implementation is pretty fast.  They can correct
| rather long error bursts and are naturals for computer data streams--you
| can use eight bit symbols in the algebra which correspond to the bytes in
| the data stream.
+---------------

Particularly simple is the shortened (174,171) R-S code over GF(256) used in
the IBM 3370 disk. [You can also use it in the unshortened (255,252) form.]
Each of the three checksum bytes is computed independently, rather than as
the remainder of one large polynomial, computed (roughly) as:

        for (s0 = s1 = sm1 = i = 0; i < BLOCK_SIZE; ++i) {
                s0  =                           s0 ^ input[i];
                s1  =         GF_mult_by_alpha[ s1 ^ input[i] ];
                sm1 = GF_mult_by_alpha_inverse[sm1 ^ input[i] ];
        }

where "GF_mult_by_alpha" is a lookup table to multiply an element of GF(256)
by the constant "alpha", the chosen primitive element. The two tables are
256 bytes each. [These checksums can also be computed in hardware at the
speed of one clock per byte, with a simple PAL per checksum byte (3).]

Correction is just as simple, but needs more memory to store the division
table (since both operands are variables). You start by computing the
"syndrome" bytes the same was as the checksum bytes, but this time the
transmitted checksum bytes are folded in too. then, assuming just one
byte in error, the location of the byte in error is given by:

        err_i = GF_log_base_alpha[ GF_divide[s1][s0] ];

where "GF_log_base_alpha" is another 256-byte table and "GF_divide" is a
65,536-byte table. Given the location "err_i" of the error, correction is
then simply a matter of:

        input[err_i] ^= s0;

(The "sm1" term is actually redundant, but can be used to help protect
against miscorrections that can occur when there are multiple errors.
It's cheap insurance.)

Although this code corrects only one symbol (byte) in error, interleaving
multiple blocks can improve the burst performance almost as well as one
likes. [The IBM 3370 uses 3:1 interleaving to correct any burst confined
to three consecutive bytes.] Still, it can only correct one burst. But the
cost is really low...

+---------------
| The guys from JPL would tell you this for no beers at all!
+---------------

Reference:

        Lin & Costello, "Error Control Coding: Fundamentals & Applications",
        (Prentice-Hall 1983), pp.525-531.

@//E*O*F rs.rply.2//
chmod u=rw,g=,o= rs.rply.2
 
echo x - rs.rply.3
sed 's/^@//' > "rs.rply.3" <<'@//E*O*F rs.rply.3//'
@From: "[email-masked] -- Advanced Media Laboratory of Samsung Electronics 1009 Lenox Dr. Lawrenceville" <[email-masked]>
Subject: Reed Solomon Code

Here is some code that came down the net some time ago.
I works really well. I have put this code into a Khoros routine.

Steve Jaffe
[email-masked]

Received: from munnari.OZ.AU by tigger.jvnc.net (5.65/1.34)
	id AA17750; Wed, 26 Jun 91 02:11:52 -0400
Received: from augean.ua.oz.au by munnari.oz.au with SMTP (5.64+1.3.1+0.50)
	id AA22760; Wed, 26 Jun 1991 16:11:30 +1000 (from [email-masked])
Date: Wed, 26 Jun 1991 16:11:30 +1000
@From: [email-masked]
Message-Id: <[email-masked]>
Received: by augean (5.61+IDA+MU/4.8N)
	id AA27938; Wed, 26 Jun 1991 15:32:11 +0930
To: [email-masked], [email-masked], [email-masked],
        [email-masked], [email-masked],
        [email-masked], [email-masked], [email-masked],
        [email-masked], [email-masked],
        [email-masked], [email-masked],
        [email-masked], [email-masked],
        [email-masked], [email-masked],
        [email-masked], [email-masked],
        [email-masked], [email-masked], [email-masked],
        [email-masked], [email-masked],
        [email-masked], [email-masked],
        [email-masked], [email-masked], [email-masked],
        [email-masked], [email-masked], [email-masked]
Subject: Re:  Reed-Solomon code

To all those people who requested the Reed-Solomon encoder/decoder program,
here it is.  I had to dig it out of my archive, and having posted that I had
it without checking it was all okay, discovered to my chagrin that it dumped
core when run after being compiled using cc -O on our new Sparc 2. (Okay on cc
with no optimize, and also on gcc with or without optimize, and ran heaps of
times on our old Dec and Sun boxes - funny how these bugs can hide for years!).
So after a night of reviewing my coding theory so I could follow the program
(it's been a while...), I found and fixed the bug (you guessed it, reading one
past the end of an array to give a compiler dependent segmentation violation!).
At the same time I put in a couple more comments so it hopefully won't be too
hard to follow with the aid of a reference text such as those mentioned in the
comments (can't recall the titles off-hand, just the authors, but they are
standard texts).
Anyhow, do what you will with it. It may or may not suit your needs, but its
one way of doing the job (I have no idea of it's speed compared to other
techniques or implementations). Hopefully it may save someone some work!

                                                Simon

----------------------- Cut here -  code follows ------------------------------

/*             rs.c        */
/* This program is an encoder/decoder for Reed-Solomon codes. Encoding is in
   systematic form, decoding via the Berlekamp iterative algorithm.
   In the present form , the constants mm, nn, tt, and kk=nn-2tt must be
   specified  (the double letters are used simply to avoid clashes with
   other n,k,t used in other programs into which this was incorporated!)
   Also, the irreducible polynomial used to generate GF(2**mm) must also be
   entered -- these can be found in Lin and Costello, and also Clark and Cain.

   The representation of the elements of GF(2**m) is either in index form,
   where the number is the power of the primitive element alpha, which is
   convenient for multiplication (add the powers modulo 2**m-1) or in
   polynomial form, where the bits represent the coefficients of the
   polynomial representation of the number, which is the most convenient form
   for addition.  The two forms are swapped between via lookup tables.
   This leads to fairly messy looking expressions, but unfortunately, there
   is no easy alternative when working with Galois arithmetic.

   The code is not written in the most elegant way, but to the best
   of my knowledge, (no absolute guarantees!), it works.
   However, when including it into a simulation program, you may want to do
   some conversion of global variables (used here because I am lazy!) to
   local variables where appropriate, and passing parameters (eg array
   addresses) to the functions  may be a sensible move to reduce the number
   of global variables and thus decrease the chance of a bug being introduced.

   This program does not handle erasures at present, but should not be hard
   to adapt to do this, as it is just an adjustment to the Berlekamp-Massey
   algorithm. It also does not attempt to decode past the BCH bound -- see
   Blahut "Theory and practice of error control codes" for how to do this.

              Simon Rockliff, University of Adelaide   21/9/89

   26/6/91 Slight modifications to remove a compiler dependent bug which hadn't
           previously surfaced. A few extra comments added for clarity.
           Appears to all work fine, ready for posting to net!

                  Notice
                 --------
   This program may be freely modified and/or given to whoever wants it.
   A condition of such distribution is that the author's contribution be
   acknowledged by his name being left in the comments heading the program,
   however no responsibility is accepted for any financial or other loss which
   may result from some unforseen errors or malfunctioning of the program
   during use.
                                 Simon Rockliff, 26th June 1991
*/

#include <math.h>
#include <stdio.h>
#define mm  4           /* RS code over GF(2**4) - change to suit */
#define nn  15          /* nn=2**mm -1   length of codeword */
#define tt  3           /* number of errors that can be corrected */
#define kk  9           /* kk = nn-2*tt  */

int pp [mm+1] = { 1, 1, 0, 0, 1} ; /* specify irreducible polynomial coeffts */
int alpha_to [nn+1], index_of [nn+1], gg [nn-kk+1] ;
int recd [nn], data [kk], bb [nn-kk] ;

void generate_gf()
/* generate GF(2**mm) from the irreducible polynomial p(X) in pp[0]..pp[mm]
   lookup tables:  index->polynomial form   alpha_to[] contains j=alpha**i;
                   polynomial form -> index form  index_of[j=alpha**i] = i
   alpha=2 is the primitive element of GF(2**mm)
*/
 {
   register int i, mask ;

  mask = 1 ;
  alpha_to[mm] = 0 ;
  for (i=0; i<mm; i++)
   { alpha_to[i] = mask ;
     index_of[alpha_to[i]] = i ;
     if (pp[i]!=0)
       alpha_to[mm] ^= mask ;
     mask <<= 1 ;
   }
  index_of[alpha_to[mm]] = mm ;
  mask >>= 1 ;
  for (i=mm+1; i<nn; i++)
   { if (alpha_to[i-1] >= mask)
        alpha_to[i] = alpha_to[mm] ^ ((alpha_to[i-1]^mask)<<1) ;
     else alpha_to[i] = alpha_to[i-1]<<1 ;
     index_of[alpha_to[i]] = i ;
   }
  index_of[0] = -1 ;
 }

void gen_poly()
/* Obtain the generator polynomial of the tt-error correcting, length
  nn=(2**mm -1) Reed Solomon code  from the product of (X+alpha**i), i=1..2*tt
*/
 {
   register int i,j ;

   gg[0] = 2 ;    /* primitive element alpha = 2  for GF(2**mm)  */
   gg[1] = 1 ;    /* g(x) = (X+alpha) initially */
   for (i=2; i<=nn-kk; i++)
    { gg[i] = 1 ;
      for (j=i-1; j>0; j--)
        if (gg[j] != 0)  gg[j] = gg[j-1]^ alpha_to[(index_of[gg[j]]+i)%nn] ;
        else gg[j] = gg[j-1] ;
      gg[0] = alpha_to[(index_of[gg[0]]+i)%nn] ;     /* gg[0] can never be zero */
    }
   /* convert gg[] to index form for quicker encoding */
   for (i=0; i<=nn-kk; i++)  gg[i] = index_of[gg[i]] ;
 }

void encode_rs()
/* take the string of symbols in data[i], i=0..(k-1) and encode systematically
   to produce 2*tt parity symbols in bb[0]..bb[2*tt-1]
   data[] is input and bb[] is output in polynomial form.
   Encoding is done by using a feedback shift register with appropriate
   connections specified by the elements of gg[], which was generated above.
   Codeword is   c(X) = data(X)*X**(nn-kk)+ b(X)          */
 {
   register int i,j ;
   int feedback ;

   for (i=0; i<nn-kk; i++)   bb[i] = 0 ;
   for (i=kk-1; i>=0; i--)
    {  feedback = index_of[data[i]^bb[nn-kk-1]] ;
       if (feedback != -1)
        { for (j=nn-kk-1; j>0; j--)
            if (gg[j] != -1)
              bb[j] = bb[j-1]^alpha_to[(gg[j]+feedback)%nn] ;
            else
              bb[j] = bb[j-1] ;
          bb[0] = alpha_to[(gg[0]+feedback)%nn] ;
        }
       else
        { for (j=nn-kk-1; j>0; j--)
            bb[j] = bb[j-1] ;
          bb[0] = 0 ;
        } ;
    } ;
 } ;

void decode_rs()
/* assume we have received bits grouped into mm-bit symbols in recd[i],
   i=0..(nn-1),  and recd[i] is index form (ie as powers of alpha).
   We first compute the 2*tt syndromes by substituting alpha**i into rec(X) and
   evaluating, storing the syndromes in s[i], i=1..2tt (leave s[0] zero) .
   Then we use the Berlekamp iteration to find the error location polynomial
   elp[i].   If the degree of the elp is >tt, we cannot correct all the errors
   and hence just put out the information symbols uncorrected. If the degree of
   elp is <=tt, we substitute alpha**i , i=1..n into the elp to get the roots,
   hence the inverse roots, the error location numbers. If the number of errors
   located does not equal the degree of the elp, we have more than tt errors
   and cannot correct them.  Otherwise, we then solve for the error value at
   the error location and correct the error.  The procedure is that found in
   Lin and Costello. For the cases where the number of errors is known to be too
   large to correct, the information symbols as received are output (the
   advantage of systematic encoding is that hopefully some of the information
   symbols will be okay and that if we are in luck, the errors are in the
   parity part of the transmitted codeword).  Of course, these insoluble cases
   can be returned as error flags to the calling routine if desired.   */
 {
   register int i,j,u,q ;
   int elp[nn-kk+2][nn-kk], d[nn-kk+2], l[nn-kk+2], u_lu[nn-kk+2], s[nn-kk+1] ;
   int count=0, syn_error=0, root[tt], loc[tt], z[tt+1], err[nn], reg[tt+1] ;

/* first form the syndromes */
   for (i=1; i<=nn-kk; i++)
    { s[i] = 0 ;
      for (j=0; j<nn; j++)
        if (recd[j]!=-1)
          s[i] ^= alpha_to[(recd[j]+i*j)%nn] ;      /* recd[j] in index form */
/* convert syndrome from polynomial form to index form  */
      if (s[i]!=0)  syn_error=1 ;        /* set flag if non-zero syndrome => error */
      s[i] = index_of[s[i]] ;
    } ;

   if (syn_error)       /* if errors, try and correct */
    {
/* compute the error location polynomial via the Berlekamp iterative algorithm,
   following the terminology of Lin and Costello :   d[u] is the 'mu'th
   discrepancy, where u='mu'+1 and 'mu' (the Greek letter!) is the step number
   ranging from -1 to 2*tt (see L&C),  l[u] is the
   degree of the elp at that step, and u_l[u] is the difference between the
   step number and the degree of the elp.
*/
/* initialise table entries */
      d[0] = 0 ;           /* index form */
      d[1] = s[1] ;        /* index form */
      elp[0][0] = 0 ;      /* index form */
      elp[1][0] = 1 ;      /* polynomial form */
      for (i=1; i<nn-kk; i++)
        { elp[0][i] = -1 ;   /* index form */
          elp[1][i] = 0 ;   /* polynomial form */
        }
      l[0] = 0 ;
      l[1] = 0 ;
      u_lu[0] = -1 ;
      u_lu[1] = 0 ;
      u = 0 ;

      do
      {
        u++ ;
        if (d[u]==-1)
          { l[u+1] = l[u] ;
            for (i=0; i<=l[u]; i++)
             {  elp[u+1][i] = elp[u][i] ;
                elp[u][i] = index_of[elp[u][i]] ;
             }
          }
        else
/* search for words with greatest u_lu[q] for which d[q]!=0 */
          { q = u-1 ;
            while ((d[q]==-1) && (q>0)) q-- ;
/* have found first non-zero d[q]  */
            if (q>0)
             { j=q ;
               do
               { j-- ;
                 if ((d[j]!=-1) && (u_lu[q]<u_lu[j]))
                   q = j ;
               }while (j>0) ;
             } ;

/* have now found q such that d[u]!=0 and u_lu[q] is maximum */
/* store degree of new elp polynomial */
            if (l[u]>l[q]+u-q)  l[u+1] = l[u] ;
            else  l[u+1] = l[q]+u-q ;

/* form new elp(x) */
            for (i=0; i<nn-kk; i++)    elp[u+1][i] = 0 ;
            for (i=0; i<=l[q]; i++)
              if (elp[q][i]!=-1)
                elp[u+1][i+u-q] = alpha_to[(d[u]+nn-d[q]+elp[q][i])%nn] ;
            for (i=0; i<=l[u]; i++)
              { elp[u+1][i] ^= elp[u][i] ;
                elp[u][i] = index_of[elp[u][i]] ;  /*convert old elp value to index*/
              }
          }
        u_lu[u+1] = u-l[u+1] ;

/* form (u+1)th discrepancy */
        if (u<nn-kk)    /* no discrepancy computed on last iteration */
          {
            if (s[u+1]!=-1)
                   d[u+1] = alpha_to[s[u+1]] ;
            else
              d[u+1] = 0 ;
            for (i=1; i<=l[u+1]; i++)
              if ((s[u+1-i]!=-1) && (elp[u+1][i]!=0))
                d[u+1] ^= alpha_to[(s[u+1-i]+index_of[elp[u+1][i]])%nn] ;
            d[u+1] = index_of[d[u+1]] ;    /* put d[u+1] into index form */
          }
      } while ((u<nn-kk) && (l[u+1]<=tt)) ;

      u++ ;
      if (l[u]<=tt)         /* can correct error */
       {
/* put elp into index form */
         for (i=0; i<=l[u]; i++)   elp[u][i] = index_of[elp[u][i]] ;

/* find roots of the error location polynomial */
         for (i=1; i<=l[u]; i++)
           reg[i] = elp[u][i] ;
         count = 0 ;
         for (i=1; i<=nn; i++)
          {  q = 1 ;
             for (j=1; j<=l[u]; j++)
              if (reg[j]!=-1)
                { reg[j] = (reg[j]+j)%nn ;
                  q ^= alpha_to[reg[j]] ;
                } ;
             if (!q)        /* store root and error location number indices */
              { root[count] = i;
                loc[count] = nn-i ;
                count++ ;
              };
          } ;
         if (count==l[u])    /* no. roots = degree of elp hence <= tt errors */
          {
/* form polynomial z(x) */
           for (i=1; i<=l[u]; i++)        /* Z[0] = 1 always - do not need */
            { if ((s[i]!=-1) && (elp[u][i]!=-1))
                 z[i] = alpha_to[s[i]] ^ alpha_to[elp[u][i]] ;
              else if ((s[i]!=-1) && (elp[u][i]==-1))
                      z[i] = alpha_to[s[i]] ;
                   else if ((s[i]==-1) && (elp[u][i]!=-1))
                          z[i] = alpha_to[elp[u][i]] ;
                        else
                          z[i] = 0 ;
              for (j=1; j<i; j++)
                if ((s[j]!=-1) && (elp[u][i-j]!=-1))
                   z[i] ^= alpha_to[(elp[u][i-j] + s[j])%nn] ;
              z[i] = index_of[z[i]] ;         /* put into index form */
            } ;

  /* evaluate errors at locations given by error location numbers loc[i] */
           for (i=0; i<nn; i++)
             { err[i] = 0 ;
               if (recd[i]!=-1)        /* convert recd[] to polynomial form */
                 recd[i] = alpha_to[recd[i]] ;
               else  recd[i] = 0 ;
             }
           for (i=0; i<l[u]; i++)    /* compute numerator of error term first */
            { err[loc[i]] = 1;       /* accounts for z[0] */
              for (j=1; j<=l[u]; j++)
                if (z[j]!=-1)
                  err[loc[i]] ^= alpha_to[(z[j]+j*root[i])%nn] ;
              if (err[loc[i]]!=0)
               { err[loc[i]] = index_of[err[loc[i]]] ;
                 q = 0 ;     /* form denominator of error term */
                 for (j=0; j<l[u]; j++)
                   if (j!=i)
                     q += index_of[1^alpha_to[(loc[j]+root[i])%nn]] ;
                 q = q % nn ;
                 err[loc[i]] = alpha_to[(err[loc[i]]-q+nn)%nn] ;
                 recd[loc[i]] ^= err[loc[i]] ;  /*recd[i] must be in polynomial form */
               }
            }
          }
         else    /* no. roots != degree of elp => >tt errors and cannot solve */
           for (i=0; i<nn; i++)        /* could return error flag if desired */
               if (recd[i]!=-1)        /* convert recd[] to polynomial form */
                 recd[i] = alpha_to[recd[i]] ;
               else  recd[i] = 0 ;     /* just output received codeword as is */
       }
     else         /* elp has degree has degree >tt hence cannot solve */
       for (i=0; i<nn; i++)       /* could return error flag if desired */
          if (recd[i]!=-1)        /* convert recd[] to polynomial form */
            recd[i] = alpha_to[recd[i]] ;
          else  recd[i] = 0 ;     /* just output received codeword as is */
    }
   else       /* no non-zero syndromes => no errors: output received codeword */
    for (i=0; i<nn; i++)
       if (recd[i]!=-1)        /* convert recd[] to polynomial form */
         recd[i] = alpha_to[recd[i]] ;
       else  recd[i] = 0 ;
 }

main()
{
  register int i;

/* generate the Galois Field GF(2**mm) */
  generate_gf() ;
  printf("Look-up tables for GF(2**%2d)\n",mm) ;
  printf("  i   alpha_to[i]  index_of[i]\n") ;
  for (i=0; i<=nn; i++)
   printf("%3d      %3d          %3d\n",i,alpha_to[i],index_of[i]) ;
  printf("\n\n") ;

/* compute the generator polynomial for this RS code */
  gen_poly() ;

/* for known data, stick a few numbers into a zero codeword. Data is in
   polynomial form.
*/
for  (i=0; i<kk; i++)   data[i] = 0 ;

/* for example, say we transmit the following message (nothing special!) */
data[0] = 8 ;
data[1] = 6 ;
data[2] = 8 ;
data[3] = 1 ;
data[4] = 2 ;
data[5] = 4 ;
data[6] = 15 ;
data[7] = 9 ;
data[8] = 9 ;

/* encode data[] to produce parity in bb[].  Data input and parity output
   is in polynomial form
*/
  encode_rs() ;

/* put the transmitted codeword, made up of data plus parity, in recd[] */
  for (i=0; i<nn-kk; i++)  recd[i] = bb[i] ;
  for (i=0; i<kk; i++) recd[i+nn-kk] = data[i] ;

/* if you want to test the program, corrupt some of the elements of recd[]
   here. This can also be done easily in a debugger. */
/* Again, lets say that a middle element is changed */
  data[nn-nn/2] = 3 ;

  for (i=0; i<nn; i++)
     recd[i] = index_of[recd[i]] ;          /* put recd[i] into index form */

/* decode recv[] */
  decode_rs() ;         /* recd[] is returned in polynomial form */

/* print out the relevant stuff - initial and decoded {parity and message} */
  printf("Results for Reed-Solomon code (n=%3d, k=%3d, t= %3d)\n\n",nn,kk,tt) ;
  printf("  i  data[i]   recd[i](decoded)   (data, recd in polynomial form)\n");
  for (i=0; i<nn-kk; i++)
    printf("%3d    %3d      %3d\n",i, bb[i], recd[i]) ;
  for (i=nn-kk; i<nn; i++)
    printf("%3d    %3d      %3d\n",i, data[i-nn+kk], recd[i]) ;
}

@//E*O*F rs.rply.3//
chmod u=rw,g=,o= rs.rply.3
 
exit 0
--------------------------------------------------------------------------------