// ------------------------------------------------------------------
// File:    rs_eedec_euc.c
// Author:  Robert Morelos-Zaragoza
// Date:    March 20, 2002
//
// The Reed-Somolon code is specified by the finite field, the length 
// (length <= 2^m-1), the number of redundant symbols (length-k), and 
// the initial zero of the code, init_zero, such that the zeros are:
//    init_zero, init_zero+1, ..., init_zero+length-k-1
//
// Update 8/17/03: Added a line to check for zero in computation of
// erasure polynomial. Many thanks to Matteo Albanese, a graduate
// student at Politecnico di Milano, for pointing this out.
//
// ERRORS AND ERASURES CORRECTION WITH EXTENDED EUCLIDEAN ALGORITHM
//
// Copyright 2002 (c) Robert H. Morelos-Zaragoza. All rights reserved
// ------------------------------------------------------------------

#include <math.h>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#define MAX_RANDOM LONG_MAX    // Maximum value of random() 

int i;
int m, n, length, k, t, t2, t21, d, red;
int init_zero;
int p[10];
int alpha_to[1024], index_of[1024], g[1024];
int recd[1024], data[1024], b[1024];
int numerr, errpos[512], errval[512], decerror;
int biterror, error;
char filename[40], name2[40];
int numera;
int era[512], eraval[512];

void read_p(void);
void generate_gf(void);
void gen_poly(void);
void encode_rs(void);
void bpsk_awgn(void);
void decode_rs(void);
int weight(int word);

main(int argc, char *argv[])
{
  // Command line processing
  if (argc != 5)
  {
  printf("\nSimulation of RS codes \n");
  printf("Copyright 2002 (c) Robert Morelos-Zaragoza\n\n");
  printf("Usage: %s m length red init_zero \n", argv[0]);
  printf(" - m is the order of GF(2^m)\n");
  printf(" - length is the length of the RS code\n");
  printf(" - red = length - dimension\n");
  printf(" - init_zero = inital zero of the code\n");
  exit(0);
  }

  sscanf(argv[1],"%d", &m);
  sscanf(argv[2],"%d", &length);
  sscanf(argv[3],"%d", &red);
  sscanf(argv[4],"%d", &init_zero);

  k = length - red;
  t = red/2;
  t2 = 2*t;
  t21 = t2+1;

  read_p();        // Read m */
  generate_gf();   // Construct the Galois Field GF(2^m) */
  printf("--> This is an RS(%d,%d,%d) code over GF(2^%d), t = %d\n",
          length, k, red+1, m, t);
  printf("    with zeros: ");
  for (i=0;i<t2;i++) printf("%d ", init_zero+i);
  printf("\n\n");

  gen_poly();      // Compute the generator polynomial of RS code */

  // Randomly generate DATA */
  for (i = 0; i < k; i++)
    data[i] = (random()>>10) % n;

  encode_rs();    // encode data */

  // Codeword is c(X) = data(X)*X**(length-k)+ b(X)
  for (i = 0; i < length - k; i++) recd[i] = b[i];
  for (i = 0; i < k; i++) recd[i + length - k] = data[i];

  for (i = 0; i < length; i++) recd[i] = 0;

  // Introduce errors and erasures at will ...

  /* CHANNEL NOISE: Read errors and erasures */
  printf("Number of errors, e: "); scanf("%d", &numerr);
  if (numerr)
    for (i=0; i<numerr; i++)
      {
      printf("position: "); scanf("%d", &errpos[i]); // as power of alpha
      printf("value:    "); scanf("%d", &errval[i]); // as vector
      recd[errpos[i]] ^= errval[i];
      }
  printf("Number of erasures, b: "); scanf("%d", &numera);
  if (numera)
    for (i=0; i<numera; i++)
      {
      printf("position: "); scanf("%d", &era[i]); // as power of alpha
      printf("value:    "); scanf("%d", &eraval[i]); // as vector
      recd[era[i]] ^= eraval[i];
      }


  printf("\n\nRec =");
  for (i=0; i<length; i++) {
     printf("%4d ", index_of[recd[i]]);
     }
  printf("\n     ");
  for (i=0; i<length; i++) {
     printf("%4d ", recd[i]);
     }
  printf("\n");

  decode_rs();   // DECODE received codeword recv[] */

  // DECODING ERRORS? we compare only the data portion
  decerror = 0;
  biterror = 0;
  for (i = length-k; i < length; i++)
    if (data[i-length+k] != recd[i])
      {
      decerror++;
      biterror += weight(data[i-length+k]- recd[i]);
      }
}


void read_p()
//      Read m, the degree of a primitive polynomial p(x) used to compute the
//      Galois field GF(2**m). Get precomputed coefficients p[] of p(x). Read
//      the code length.
{
  int i, ninf;

  printf("\nSimulation of RS codes \n");
  printf("Copyright 2002 (c)Robert Morelos-Zaragoza. All rights reserved.\n\n");
  for (i=1; i<m; i++)
    p[i] = 0;
  p[0] = p[m] = 1;
  if (m == 2)             p[1] = 1;
  else if (m == 3)        p[1] = 1;
  else if (m == 4)        p[3] = 1;
  // else if (m == 4)        p[1] = 1;  // Commented out to match example p. 68
  else if (m == 5)        p[2] = 1;
  else if (m == 6)        p[1] = 1;
  else if (m == 7)        p[1] = 1;
  else if (m == 8)        p[4] = p[5] = p[6] = 1;
  else if (m == 9)        p[4] = 1;
  else if (m == 10)       p[3] = 1;
  else if (m == 11)       p[2] = 1;
  else if (m == 12)       p[3] = p[4] = p[7] = 1;
  else if (m == 13)       p[1] = p[3] = p[4] = 1;
  else if (m == 14)       p[1] = p[11] = p[12] = 1;
  else if (m == 15)       p[1] = 1;
  else if (m == 16)       p[2] = p[3] = p[5] = 1;
  else if (m == 17)       p[3] = 1;
  else if (m == 18)       p[7] = 1;
  else if (m == 19)       p[1] = p[5] = p[6] = 1;
  else if (m == 20)       p[3] = 1;
  printf("Primitive polynomial of GF(2^%d), (LSB first)   p(x) = ",m);
  n = 1;
  for (i = 0; i <= m; i++) 
    {
      n *= 2;
      printf("%1d", p[i]);
    }
  printf("\n");
  n = n / 2 - 1;
}



void generate_gf()
// generate GF(2^m) from the irreducible polynomial p(X) in p[0]..p[m]
//
// lookup tables:  log->vector form           alpha_to[] contains j=alpha**i;
//                 vector form -> log form    index_of[j=alpha**i] = i
// alpha=2 is the primitive element of GF(2^m)
{
register int i, mask; 
  mask = 1; 
  alpha_to[m] = 0; 
  for (i=0; i<m; i++)
   { 
     alpha_to[i] = mask; 
     index_of[alpha_to[i]] = i; 
     if (p[i]!=0)
       alpha_to[m] ^= mask; 
     mask <<= 1; 
   }
  index_of[alpha_to[m]] = m; 
  mask >>= 1; 
  for (i=m+1; i<n; i++)
   { 
     if (alpha_to[i-1] >= mask)
        alpha_to[i] = alpha_to[m] ^ ((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; 

#ifdef PRINT_GF
printf("Table of GF(%d):\n",n+1);
printf("----------------------\n");
printf("   i\tvector \tlog\n");
printf("----------------------\n");
for (i=0; i<=n; i++)
printf("%4d\t%4x\t%4d\n", i, alpha_to[i], index_of[i]);
#endif
 }


void gen_poly()
// Compute the generator polynomial of the t-error correcting, length
// n=(2^m -1) Reed-Solomon code from the product of (X+alpha^i), for
// i = init_zero, init_zero + 1, ..., init_zero+length-k-1
{
register int i,j; 

   g[0] = alpha_to[init_zero];  //  <--- vector form of alpha^init_zero
   g[1] = 1;     // g(x) = (X+alpha^init_zero)
   for (i=2; i<=length-k; i++)
    { 
      g[i] = 1;
      for (j=i-1; j>0; j--)
        if (g[j] != 0)  
          g[j] = g[j-1]^ alpha_to[(index_of[g[j]]+i+init_zero-1)%n]; 
        else 
          g[j] = g[j-1]; 
      g[0] = alpha_to[(index_of[g[0]]+i+init_zero-1)%n];
    }
   // convert g[] to log form for quicker encoding 
   for (i=0; i<=length-k; i++)  g[i] = index_of[g[i]]; 
#ifdef PRINT_POLY
printf("Generator polynomial (independent term first):\ng(x) = ");
for (i=0; i<=length-k; i++) printf("%5d", g[i]);
printf("\n");
#endif
}


void encode_rs()
// Compute the 2t parity symbols in b[0]..b[2*t-1]
// data[] is input and b[] is output in polynomial form.
// Encoding is done by using a feedback shift register with connections
// specified by the elements of g[].
{
   register int i,j; 
   int feedback; 

   for (i=0; i<length-k; i++)   
     b[i] = 0; 
   for (i=k-1; i>=0; i--)
    {
    feedback = index_of[data[i]^b[length-k-1]]; 
      if (feedback != -1)
        { 
        for (j=length-k-1; j>0; j--)
          if (g[j] != -1)
            b[j] = b[j-1]^alpha_to[(g[j]+feedback)%n]; 
          else
            b[j] = b[j-1]; 
          b[0] = alpha_to[(g[0]+feedback)%n]; 
        }
       else
        { 
        for (j=length-k-1; j>0; j--)
          b[j] = b[j-1]; 
        b[0] = 0; 
        }
    }
}



void decode_rs()
{
   register int i,j,u,q;
   int qt[513], r[129][513];
   int b[12][513];
   int degr[129], degb[129];
   int elp[129][1024], l[1], s[1025], forney[1025];
   int count=0, syn_error=0, tau[512], root[512], loc[512], z[513]; 
   int err[1024], reg[513], aux[513], omega[1025], phi[1025], phiprime[1025];
   int degphi, ell, temp;

   // Compute the syndromes

#ifdef PRINT_SYNDROME
   printf("\ns =         0 ");
#endif

   for (i=1; i<=t2; i++)
    { 
      s[i] = 0; 
      for (j=0; j<length; j++)
        if (recd[j]!=0)
          s[i] ^= alpha_to[(index_of[recd[j]]+(i+init_zero-1)*j)%n];
      // convert syndrome from vector form to log form  */
      if (s[i]!=0)  
        syn_error=1;         // set flag if non-zero syndrome => error
      //
      // Note:    If the code is used only for ERROR DETECTION, then
      //          exit program here indicating the presence of errors.
      //
      s[i] = index_of[s[i]]; 

#ifdef PRINT_SYNDROME
   printf("%4d ", s[i]);
#endif

    }

   if (syn_error)       // if syndromes are nonzero then try to correct
    {

     s[0] = 0;  // S(x) = 1 + s_1x + ...

      // TO HANDLE ERASURES

      if (numera)
      // if erasures are present, compute the erasure locator polynomial, tau(x)
        {
          for (i=0; i<=t2; i++) 
            { tau[i] = 0; aux[i] = 0;}
          
          aux[1] = alpha_to[era[0]];
          aux[0] = 1;       // (X + era[0]) 

          if (numera>1)

            for (i=1; i<numera; i++)
              {
              p[1] = era[i];
              p[0] = 0;
              for (j=0; j<2; j++)
                for (ell=0; ell<=i; ell++)
                  // Line below added 8/17/2003
                  if ((p[j] !=-1) && (aux[ell]!=0))
                    tau[j+ell] ^= alpha_to[(p[j]+index_of[aux[ell]])%n];
              if (i != (numera-1))
              for (ell=0; ell<=(i+1); ell++)
                {
                aux[ell] = tau[ell];
                tau[ell] = 0;
                }
              }

          else {
            tau[0] = aux[0]; tau[1] = aux[1];
          }

        // Put in index (log) form
        for (i=0; i<=numera; i++)
          tau[i] = index_of[tau[i]]; /* tau in log form */

#ifdef PRINT_SYNDROME
printf("\ntau =    ");
for (i=0; i<=numera; i++)
  printf("%4d ", tau[i]);
printf("\nforney = ");
#endif

        // Compute FORNEY modified syndrome:
        //            forney(x) = [ s(x) tau(x) + 1 ] mod x^{t2}

        for (i=0; i<=n-k; i++) forney[i] = 0;
    
        for (i=0; i<=n-k; i++)
          for (j=0; j<=numera; j++)
            if (i+j <= (n-k)) // mod x^{n-k+1}
              if ((s[i]!=-1)&&(tau[j]!=-1))
                forney[i+j] ^= alpha_to[(s[i]+tau[j])%n];

        forney[0] ^= 1; 

        for (i=0; i<=n-k; i++)
          forney[i] = index_of[forney[i]];

#ifdef PRINT_SYNDROME
for (i=0; i<=n-k; i++)
  printf("%4d ", forney[i]);
#endif

   }
   else // No erasures
   {
     tau[0]=0;     // i.e., tau(x) = 1
     for (i=1; i<=n-k; i++) forney[i] = s[i];
   }

#ifdef PRINT_SYNDROME
printf("\n");
#endif


  // --------------------------------------------------------------
  //    THE EUCLIDEAN ALGORITHM FOR ERRORS AND ERASURES
  // --------------------------------------------------------------

  for (i=0; i<= t21; i++) {
    r[0][i] = 0;
    r[1][i] = 0;
    b[0][i] = 0;
    b[1][i] = 0;
    qt[i] = 0;
   }

  b[1][0] = 1; degb[0] = 0; degb[1] = 0;

  r[0][t21] = 1; // x^{2t+1}
  degr[0] = t21;

  for (i=0; i<=t2; i++)
    {
    if (forney[i]!=-1) {
      r[1][i] = alpha_to[forney[i]];    // Forney's syndromes, T(X)
      degr[1] = i;
      }
    else
      r[1][i] = 0;
    }


  j = 1;

  if (numera<t2) {  // If errors can be corrected

printf("\nEuclidean algorithm:\n");
printf("r[0] = "); for (i=0; i<=degr[0];i++) printf("%4d ", index_of[r[0][i]]);
printf("\nr[1] = "); for (i=0; i<=degr[1];i++) printf("%4d ", index_of[r[1][i]]);
printf("\n");


  if( (degr[0]-degr[1]) < ((t2-numera)/2+1) ) {

  do {

    j++;

printf("\n************************  j=%3d\n", j);
    // ----------------------------------------
    // Apply long division to r[j-2] and r[j-1]
    // ----------------------------------------

    // Clean r[j]
    for (i=0; i<=t21; i++) r[j][i] = 0;

    for (i=0;i<=degr[j-2];i++) 
      r[j][i] = r[j-2][i]; 
    degr[j] = degr[j-2];

    temp = degr[j-2]-degr[j-1];
    for (i=temp; i>=0; i--) {
      u = degr[j-1]+i;
      if (degr[j] == u)
        {
        if ( r[j][degr[j]] && r[j-1][degr[j-1]] )
          qt[i] = alpha_to[(index_of[r[j][degr[j]]]
                                       +n-index_of[r[j-1][degr[j-1]]])%n];

//printf("r[j][degr[j]]] = %d,  r[j-1][degr[j-1]] = %d\n",
//index_of[r[j][degr[j]]], index_of[r[j-1][degr[j-1]]]);
printf("\nqt[%d]=%d\n", i, index_of[qt[i]]);

        for (u=0; u<=t21; u++) aux[u] = 0;

        temp = degr[j-1];
        for (u=0; u<=temp; u++)
          if ( qt[i] && r[j-1][u] )
            aux[u+i] = alpha_to[(index_of[qt[i]]+index_of[r[j-1][u]])%n];
          else
            aux[u+i] = 0;

printf("r    = ");
for (u=0; u<=degr[j]; u++) printf("%4d ", index_of[r[j][u]]);
printf("\n");
printf("aux  = ");
for (u=0; u<=degr[j-1]+i; u++) printf("%4d ", index_of[aux[u]]);
printf("\n");

        for (u=0; u<=degr[j]; u++)
          r[j][u] ^= aux[u];
        u = t21;
        while ( !r[j][u] && (u>0)) u--;
        degr[j] = u;
        }
      else
        qt[i] = 0;

printf("r    = ");
for (u=0; u<=degr[j]; u++) printf("%4d ", index_of[r[j][u]]);
printf("\n");

      }

printf("\nqt = ",j);
temp = degr[j-2]-degr[j-1];
for (i=0; i<=temp; i++) printf("%4d ", index_of[qt[i]]);
printf("\nr = ");
for (i=0; i<=degr[j]; i++) printf("%4d ", index_of[r[j][i]]);
printf("\nb = ");

    // Compute b(x)

    for (i=0; i<=t21; i++) 
      aux[i] = 0; 

    temp = degr[j-2]-degr[j-1];
    for (i=0; i<=temp; i++)
      for (u=0; u<=degb[j-1]; u++)
        if ( qt[i] && b[j-1][u] )
          aux[i+u] ^= alpha_to[(index_of[qt[i]]+index_of[b[j-1][u]])%n];

    for (i=0; i<=t21; i++) 
      b[j][i] = b[j-2][i] ^ aux[i];

    u = t21;
    while ( !b[j][u] && (u>0) ) u--;
    degb[j] = u;

for (i=0; i<=degb[j]; i++) printf("%4d ", index_of[b[j][i]]);
printf("\n");

  } while (degr[j]>t+numera/2); 

  }

  } // end if (numera<t2+1)



  u=1;
  temp = degb[j];
  // Normalize error locator polynomial
  for (i=0;i<=temp;i++) {
    elp[u][i] = alpha_to[(index_of[b[j][i]]-index_of[b[j][0]]+n)%n];
    }
  l[u] = temp;


  

printf("\nEuclidean algorithm, after %d iterations:\nsigma = ", (j-1));
for (i=0; i<=l[u]; i++)   printf("%4d ", index_of[elp[u][i]]);
printf("\n");

  // FIND ROOTS OF SIGMA

      if (l[u]<=t-numera/2)         // can correct errors
       {
         // 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<=n; i++)
          {  
             q = 1; 
             for (j=1; j<=l[u]; j++)
              if (reg[j]!=-1)
                { 
                  reg[j] = (reg[j]+j)%n; 
                  q ^= alpha_to[reg[j]]; 
                }
             if (!q)        // store root and error location number indices 
              { 
                root[count] = i;
                loc[count] = n-i; 
#ifdef PRINT_SYNDROME
printf("loc[%4d] = %4d\n", count, loc[count]);
#endif
                count++;
              }
           }

         if (count==l[u])    // no. roots = degree of elp hence <= t errors 
          {

            // Compute the errata evaluator polynomial, omega(x)

            forney[0] = 0;  // as a log, to construct 1+T(x)
            for (i=1; i<=t2; i++) 
              omega[i] = 0;
            for (i=0; i<=t2; i++)
              {
              for (j=0; j<=l[u]; j++)
                {
                if (i+j <= t2) // mod x^{t2}
                  if ((forney[i]!=-1) && (elp[u][j]!=-1))
                    omega[i+j] ^= alpha_to[(forney[i]+elp[u][j])%n];
                }
              }
            for (i=0; i<=t2; i++)
              omega[i] = index_of[omega[i]];

#ifdef PRINT_SYNDROME
printf("\nomega =    ");
for (i=0; i<=t2; i++)   printf("%4d ", omega[i]);
printf("\n");
#endif


           // Compute the errata locator polynomial, phi(x)

           degphi = numera+l[u];

           for (i=1; i<=degphi; i++) phi[i] = 0;
           for (i=0; i<=numera; i++)
             for (j=0; j<=l[u]; j++)
               if ((tau[i]!=-1)&&(elp[u][j]!=-1))
                 phi[i+j] ^= alpha_to[(tau[i]+elp[u][j])%n];
           for (i=0; i<=degphi; i++)
             phi[i] = index_of[phi[i]];

#ifdef PRINT_SYNDROME
printf("phi =      ");
for (i=0; i<=degphi; i++)   printf("%4d ", phi[i]);
printf("\n");
#endif


           // Compute the "derivative" of phi(x): phiprime

           for (i=0; i<=degphi; i++) phiprime[i] = -1; // as a log
           for (i=0; i<=degphi; i++)
             if (i%2)  // Odd powers of phi(x) give terms in phiprime(x)
               phiprime[i-1] = phi[i];

#ifdef PRINT_SYNDROME
printf("phiprime = ");
for (i=0; i<=degphi; i++)   printf("%4d ", phiprime[i]);
printf("\n\n");
#endif

for (i=0; i<count; i++)
printf("loc[%d] = %d\n", i, loc[i]);
printf("\n");

            if (numera)
            // Add erasure positions to error locations
            for (i=0; i<numera; i++) {
              loc[count+i] = era[i];
              root[count+i] = (n-era[i])%n;
              }


           // evaluate errors at locations given by errata locations, loc[i]
           //                             for (i=0; i<l[u]; i++)    
           for (i=0; i<degphi; i++)    
            { 
              // compute numerator of error term  
              err[loc[i]] = 0;
              for (j=0; j<=t2; j++)
                if ((omega[j]!=-1)&&(root[i]!=-1))
                  err[loc[i]] ^= alpha_to[(omega[j]+j*root[i])%n];

              // -------  The term loc[i]^{2-init_zero}
              if ((err[loc[i]]!=0)&&(loc[i]!=-1))
                 err[loc[i]] = alpha_to[(index_of[err[loc[i]]]
                                              +loc[i]*(2-init_zero+n))%n];

              if (err[loc[i]]!=0)
               { 
                 err[loc[i]] = index_of[err[loc[i]]];
                 // compute denominator of error term 
                 q = 0;             
                 for (j=0; j<=degphi; j++)
                   if ((phiprime[j]!=-1)&&(root[i]!=-1))
                     q ^= alpha_to[(phiprime[j]+j*root[i])%n];

                 // Division by q
                 err[loc[i]] = alpha_to[(err[loc[i]]-index_of[q]+n)%n];

#ifdef PRINT_SYNDROME
printf("errata[%4d] = %4d (%4d) \n",loc[i],index_of[err[loc[i]]],err[loc[i]]);
#endif

                 recd[loc[i]] ^= err[loc[i]];  


               }

            }


          printf("\nCor =");
          for (i=0; i<length; i++) {
             printf("%4d ", index_of[recd[i]]);
             }
          printf("\n     ");
          for (i=0; i<length; i++) {
             printf("%4d ", recd[i]);
             }
          printf("\n");

          }
         else    // no. roots != degree of elp => >t errors and cannot solve
          ;
       }
     else         // elp has degree has degree >t hence cannot solve 
       ;
    }
   else       // no non-zero syndromes => no errors: output received codeword 
     ;
}



int weight(int word)
{
  int i, res;
  res = 0;
  for (i=0; i<m; i++)
    if ( (word>>i) & 0x01 ) res++;
  return(res);
}