/*
 * File:     general_RS.c
 * Title:    Encoder/decoder for RS codes in C
 * Authors:  Robert Morelos-Zaragoza (robert@spectra.eng.hawaii.edu)
 *           and Hari Thirumoorthy (harit@spectra.eng.hawaii.edu)
 * Date:     Aug 1995
 *
 *
 * ===============  Encoder/Decoder for RS codes in C =================
 *
 *
 * The encoding and decoding methods used in this program are based on the
 * book "Error Control Coding: Fundamentals and Applications", by Lin and
 * Costello, Prentice Hall, 1983. Most recently on "Theory and Practice of 
 * Error Control Codes", by R.E. Blahut.
 *
 *
 * NOTE:    
 *          The authors are not responsible for any malfunctioning of
 *          this program, nor for any damage caused by it. Please include the
 *          original program along with these comments in any redistribution.
 *
 * Portions of this program are from a Reed-Solomon encoder/decoder
 * in C, written by Simon Rockliff (simon@augean.ua.oz.au) on 21/9/89.
 *  For more information, suggestions, or other ideas on implementing error
 *  correcting codes, please contact the authors at:
 *
 *        Robert Morelos-Zaragoza         Hari Thirumoorthy
 *        770 S. Post Oak Ln. #200        2540 Dole Street, # 483
 *        Houston, TX 77056               Dept. of Electrical Engineering,
 *                                        University of Hawaii
 *                                        Honolulu, HI 96822
 *
 * COPYRIGHT NOTICE: This computer program is free for non-commercial purposes.
 * You may implement this program for any non-commercial application. You may 
 * also implement this program for commercial purposes, provided that you
 * obtain my written permission. Any modification of this program is covered
 * by this copyright.
 *
 * Copyright (c) 1995.  Robert Morelos-Zaragoza and Hari Thirumoorthy.
 *                      All rights reserved.
 *
 */

/*
  Program computes the generator polynomial of a RS code. Also performs 
encoding and decoding of the RS code or a shortened RS code. Compile
using one of the following options:

  In the commands below, using  the [-DNO_PRINT] option will ensure that
 the program runs without printing a lot of detail that is generally
 unnecessary. If the program is compiled without using the option a lot
 of detail will be printed. Some other options, for example changing the
 probability of symbol error are available by directly editing the program 
 itself. 

 1. cc general_RS.c -DNO_PRINT -lm -o rsgp
        "rsgp" computes the generator polynomial of a RS code whose
        length is "nn", dimension is "kk" and error correcting capability
        is "tt". The generator polynomial is displayed in two formats.
        In one the coefficients of the generator polynomial are given
        as exponents of the primitive element @ of the field. In the other
        the coefficients are given by their representations in the
        basis {@^7,...,@,1}.

 2. cc general_RS.c -DDEBUG -DNO_PRINT -lm -o rsfield
        "rsfield" does all that rsgp does and in addition it also
        displays the field GF(n) and the two representations of the
        field elements namely the "power of @" representation and
        the "Basis" representation.

 3. cc general_RS.c -DENCODE -DNO_PRINT -lm -o rsenc
        "rsenc" does all that rsgp does and in addition it prompts the
        user for the following data. (1) The length of the shortened
        RS code. Of course, if shortening is not desired then, the integer
        nn is entered as the length. (2) The number of codewords that are
        to be generated. (3) The name of the file into which they are to
        be stored. (4) Random number seeds that are needed to ensure true
        randomness in generation of the codewords. (The DowJones closing
        index is a good number to enter!).  These codewords can then be
        transmitted over a noisy channel and decoded using the rsdec program.

        The error correcting capability of the code cannot be set by running
        the program. It must be set in this file "general_RS.c" by changing
        the "tt" parameter. Each time the "tt" is
        changed, the dimension "kk" of the unshortened RS code changes and
        this is given by kk = nn - 2 * tt. This value of kk must be set in
        the program below. After both "tt" and "kk" are set, the program must
        be recompiled.

 
 4.  cc general_RS.c -DDECODE -DNO_PRINT -lm -o rsdec
        "rsdec" does all that rsgp does and in addition it prompts the user
        for the following data. (1) The length of the RS code. (2) The
        file from which the codewords are to be read. (3) The number of code
        words to read. (4) The probability of symbol error. The program then
        reads in codewords and transmits each codeword over a discrete
        memoryless channel that has a symbol error probability specified
        above. The received corrupted word is then decoded and the results
        of decoding displayed. This basically verifies that the decoding
        procedure is capable of producing the transmitted codeword provided
        the channel caused at most tt errors. Almost always when the channel
        causes more than tt errors, the decoder fails to produce a codeword.
        In such a case the decoder outputs the uncorrected information bytes
        as they were received from the channel.

5.  cc general_RS.c -DDECODE -DDECODER_DEBUG -lm -o rsddebug
        "rsddebug" is similar to rsdec and can be used to debug the
        the decoder and verify its operation. It reads in one codeword
        from a file specified by the user and allows the user to specify
        the number of errors to be caused and their locations. Then the
        decoder transmits the codeword and causes the specified number
        of errors at the specified locations. The actual value of the
        error byte is randomnly generated. The received word is the sum
        of the transmitted codeword and the error word. The decoder then
        decodes this received word and displays all the intermediate
        results and the final result as explained below.
 
        (1) The syndrome of the received word is shown. This consists of
        2*tt bytes {S(1),S(2), ... ,S(2*tt)}. Each S(i) is an element in
        the field GF(n). NOTE: Iff all the S(i)'s are zero, the received
        word is a codeword in the code. In such a case, the channel
        presumably caused no errors.
 
  INPUT:
	Change the Galois Field and error correcting capability by 
	setting the appropriate global variables below including
	1. primitive polynomial 2. error correcting capability 
	3. dimension 4. length.
  
  FUNCTIONS:
		generate_gf() generates the field.
		gen_poly() generates the generator polynomial.
		encode_rs() encodes in systematic form.
		decode_rs() errors-only decodes a vector assumed to be encoded in systematic form.
	        eras_dec_rs() errors+erasures decoding of a vector assumed to be encoded in 
		 	      systematic form.
*/

#include <math.h>
#include <stdio.h>
#define maxrand 0x7fffffff
#define mm  8           /* RS code over GF(2**mm) - change to suit */
#define n   256   	/* n = size of the field */
#define nn  255         /* nn=2**mm -1   length of codeword */
#define tt  16          /* number of errors that can be corrected */
#define kk  223         /* kk = nn-2*tt  */ /* Degree of g(x) = 2*tt */

/**** Primitive polynomial ****/
int pp [mm+1] = { 1, 0, 1, 1, 1, 0, 0, 0, 1}; /* 1+x^2+x^3+x^4+x^8 */

/* int pp [mm+1] = { 1, 1, 0, 0, 1}; */ /* 1 + x + x^4 */

/* int pp[mm+1] = {1,1,0,1}; */ /* 1 + x + x^3 */

/* generator polynomial, tables for Galois field */
int alpha_to [nn+1], index_of [nn+1], gg [nn-kk+1];
int b0;  /* g(x) has roots @**b0, @**(b0+1), ... ,@^(b0+2*tt-1) */

/* data[] is the info vector, bb[] is the parity vector, recd[] is the 
  noise corrupted received vector  */
int recd[nn], data[kk], bb[nn-kk] ;

/* generate GF(2**m) from the irreducible polynomial p(X) in p[0]..p[m]
   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**m)
   HARI's COMMENT: (4/13/94) alpha_to[] can be used as follows:
        Let @ represent the primitive element commonly called "alpha" that
   is the root of the primitive polynomial p(x). Then in GF(2^m), for any
   0 <= i <= 2^m-2,
        @^i = a(0) + a(1) @ + a(2) @^2 + ... + a(m-1) @^(m-1)
   where the binary vector (a(0),a(1),a(2),...,a(m-1)) is the representation
   of the integer "alpha_to[i]" with a(0) being the LSB and a(m-1) the MSB. Thus for
   example the polynomial representation of @^5 would be given by the binary
   representation of the integer "alpha_to[5]".
                   Similarily, index_of[] can be used as follows:
        As above, let @ represent the primitive element of GF(2^m) that is
   the root of the primitive polynomial p(x). In order to find the power
   of @ (alpha) that has the polynomial representation
        a(0) + a(1) @ + a(2) @^2 + ... + a(m-1) @^(m-1)
   we consider the integer "i" whose binary representation with a(0) being LSB
   and a(m-1) MSB is (a(0),a(1),...,a(m-1)) and locate the entry
   "index_of[i]". Now, @^index_of[i] is that element whose polynomial 
    representation is (a(0),a(1),a(2),...,a(m-1)).
   NOTE:
        The element alpha_to[2^m-1] = 0 always signifying that the
   representation of "@^infinity" = 0 is (0,0,0,...,0).
        Similarily, the element index_of[0] = -1 always signifying
   that the power of alpha which has the polynomial representation
   (0,0,...,0) is "infinity".
 
*/

void generate_gf()
 {
   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) /* If pp[i] == 1 then, term @^i occurs in poly-repr of @^mm */
       alpha_to[mm] ^= mask ;  /* Bit-wise EXOR operation */
     mask <<= 1 ; /* single left-shift */
   }
  index_of[alpha_to[mm]] = mm ;
  /* Have obtained poly-repr of @^mm. Poly-repr of @^(i+1) is given by 
     poly-repr of @^i shifted left one-bit and accounting for any @^mm 
     term that may occur when poly-repr of @^i is shifted. */
  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+@**(b0+i)), i = 0, ... ,(2*tt-1)
  Examples: 	If b0 = 1, tt = 1. deg(g(x)) = 2*tt = 2.
 	g(x) = (x+@) (x+@**2)
		If b0 = 0, tt = 2. deg(g(x)) = 2*tt = 4.
	g(x) = (x+1) (x+@) (x+@**2) (x+@**3)  	
*/
 {
   register int i,j ;

   gg[0] = alpha_to[b0];
   gg[1] = 1 ;    /* g(x) = (X+@**b0) initially */
   for (i=2; i <= nn-kk; i++)
    { gg[i] = 1 ;
      /* Below multiply (gg[0]+gg[1]*x + ... +gg[i]x^i) by (@**(b0+i-1) + x) */
      for (j=i-1; j>0; j--)
        if (gg[j] != 0)  gg[j] = gg[j-1]^ alpha_to[((index_of[gg[j]])+b0+i-1)%nn] ;
        else gg[j] = gg[j-1] ;
      gg[0] = alpha_to[((index_of[gg[0]])+b0+i-1)%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) /* feedback term is non-zero */
        { 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 /* feedback term is zero. encoder becomes a single-byte shifter */
        { for (j=nn-kk-1; j>0; j--)
            bb[j] = bb[j-1] ;
          bb[0] = 0 ;
        } ;
    } ;
 }


/* Performs ERRORS-ONLY decoding of RS codes. If decoding is successful,
 writes the codeword into recd[] itself. Otherwise recd[] is unaltered. 
 If channel caused no more than "tt" errors, the tranmitted codeword will
 be recovered.
*/ 

int decode_rs()
 {
   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] ;

/* recd[] is in polynomial form, convert to index form */
   for (i=0;i < nn;i++) recd[i] = index_of[recd[i]];

/* first form the syndromes; i.e., evaluate recd(x) at roots of g(x) namely
 @**(b0+i), i = 0, ... ,(2*tt-1) */
   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]+(b0+i-1)*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]] ;
    };

#ifdef DECODER_DEBUG
    printf("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
    printf("The Syndromes of the received vector in polynomial-form are:\n");
    for (i=1;i <= nn-kk;i++){
        if (s[i] == -1)
           printf("s(%d) = %#x ",i,0);
	else
           printf("s(%d) = %#x ",i,s[i]);
    }
    printf("\n NOTE: The basis is {@^%d,...,@,1}\n",mm-1);
#endif   


#ifdef DECODER_DEBUG
    printf("The Syndromes of the received vector in index-form are:\n");
    for (i=1;i <= nn-kk;i++){
        printf("s(%d) = %d ",i,s[i]);
    } 
    printf("\n NOTE: Since 0 is not expressible as a power of the primitive element\n");
    printf("alpha (denoted @),  -1 is printed if the corresponding syndrome is zero.\n");
    printf("---------------------------------------------------------------\n");
#endif

   if (syn_error)       /* if errors, try and correct */
    {
#ifdef DECODER_DEBUG
    printf("Received vector is not a codeword. Channel has caused at least one error\n");
#endif

/* 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.
  
   The notation is the same as that in Lin and Costello's book; pages 155-156 and 175. 

*/
/* 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 ;