rs_eedec_bm_c.c

很不错的RS编解码程序

        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;
    for (i=1; i<=n-k; i++) forney[i] = s[i];
    }

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


  // --------------------------------------------------------------
  //    THE BERLEKAMP-MASSEY ALGORITHM FOR ERRORS AND ERASURES
  // --------------------------------------------------------------

      // initialize table entries
      d[0] = 0;                // log form
      d[1] = forney[numera+1]; // log form
      elp[0][0] = 0;           // log form 
      elp[1][0] = 1;           // vector form 
      for (i=1; i<t2; i++)
        { 
          elp[0][i] = -1;   // log form
          elp[1][i] = 0;    // vector form 
        }
      l[0] = 0; 
      l[1] = 0; 
      u_lu[0] = -1; 
      u_lu[1] = 0; 
      u = 0; 

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

      do
      {
        u++; 
        if (d[u]==-1)
          { 
#ifdef PRINT_SYNDROME
printf("d[%d] is zero\n",u);
#endif
            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); 
             }

#ifdef PRINT_SYNDROME
printf("u = %4d, q = %4d, d[q] = %4d d[u] = %4d\n", u, q, d[q],d[u]);
#endif

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

#ifdef PRINT_SYNDROME
printf("l[q] = %4d, l[u] = %4d\n", l[q], l[u]);
#endif

            // compute new elp(x) 
            // for (i=0; i<t2-numera; i++)    elp[u+1][i] = 0; 
            for (i=0; i<t2; 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]+n-d[q]+elp[q][i])%n]; 
            for (i=0; i<=l[u]; i++)
              { 
                elp[u+1][i] ^= elp[u][i]; 
                elp[u][i] = index_of[elp[u][i]];   
              }

#ifdef PRINT_SYNDROME
printf("l[u+1] = %4d, elp[u+1] = ", l[u+1]);
for (i=0;  i<=l[u+1]; i++) printf("%4d ",index_of[elp[u+1][i]]);
printf("\n");
#endif

          }
        u_lu[u+1] = u-l[u+1]; 
        // compute (u+1)th discrepancy 
        // if (u<t2)    // no discrepancy computed on last iteration 
        if (u<(t2-numera)) // no discrepancy computed on last iteration 
        // --- if ( u < (l[u+1]+t-1-(numera/2)) ) 
          {
            // if (s[u+1]!=-1)
            if (forney[numera+u+1]!=-1)
              d[u+1] = alpha_to[forney[numera+u+1]]; 
            else
              d[u+1] = 0; 
#ifdef PRINT_SYNDROME
printf("discrepancy for u = %d: d[u+1] = %4d\n", u, index_of[d[u+1]]);
#endif

            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[numera+u+1-i]
              if ((forney[numera+u+1-i]!=-1) && (elp[u+1][i]!=0)) {
                 d[u+1] ^= alpha_to[(forney[numera+u+1-i]
                                      + index_of[elp[u+1][i]])%n]; 
#ifdef PRINT_SYNDROME
printf("i=%d, forney[%d] = %4d, d[u+1] = %4d\n",i,numera+u+1-i,
                   forney[numera+u+1-i],index_of[d[u+1]]);
#endif
                }
            d[u+1] = index_of[d[u+1]];     // put d[u+1] into index form 

#ifdef PRINT_SYNDROME
printf("d[u+1] = %4d\n", d[u+1]);
#endif

          }
      // } while ((u<t2) && (l[u+1]<=t)); 
      } while ((u<(t2-numera)) && (l[u+1]<=((t2-numera)/2))); 
      }

      // else // case of 2t erasures
      // {
      //   elp[1][0] = 0;
      //   count = 0;
      // }

      u++; 

      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]]; 

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

         // 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=0; 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=0; 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

printf("**** count = %d\n", count);

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