rs_eedec_euc.c

The Reed-Somolon code is specified by the finite field, the length (length <= 2^m-1), the numbe

#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);
}