berlekamp.c

rs

//#include <stdio.h>
#include "rs.h"




static int Lambda[MAXDEG];

static int Omega[MAXDEG];

static unsigned char ErrorLocs[128];
static int NErrors;

static unsigned char ErasureLocs[128];
static int NErasures;



static int compute_discrepancy(int lambda[], int S[], int L, int n);
static void init_gamma(int gamma[]);
static void compute_modified_omega (int NPAR);
static void mul_z_poly (int src[]);


void Modified_Berlekamp_Massey (int NPAR)
{	
  int n, L, L2, k, d, i;
  int psi[MAXDEG], psi2[MAXDEG], D[MAXDEG];
  int gamma[MAXDEG];
	
  init_gamma(gamma);

  copy_poly(D, gamma);
  mul_z_poly(D);
	
  copy_poly(psi, gamma);	
  k = -1; L = NErasures;
	
  for (n = NErasures; n < NPAR; n++) 
  {
	
    d = compute_discrepancy(psi, synBytes, L, n);
		
    if (d != 0) 
    {
		
      for (i = 0; i < MAXDEG; i++) psi2[i] = psi[i] ^ gmult(d, D[i]);
		
		
      if (L < (n-k))
	  {
	L2 = n-k;
	k = n-L;

	for (i = 0; i < MAXDEG; i++) D[i] = gmult(psi[i], ginv(d));
	L = L2;
      }
			
      for (i = 0; i < MAXDEG; i++) psi[i] = psi2[i];
    }
		
    mul_z_poly(D);
  }
	
  for(i = 0; i < MAXDEG; i++) Lambda[i] = psi[i];
  compute_modified_omega(NPAR);

	
}


void compute_modified_omega (int NPAR)
{
  int i;
  int product[MAXDEG*2];
	
  mult_polys(product, Lambda, synBytes);	
  zero_poly(Omega);
  for(i = 0; i < NPAR; i++) Omega[i] = product[i];

}

void mult_polys(int dst[],int pp1[],int pp2[])
{
  int i, j;
  int tmp1[MAXDEG*2];
	
  for (i=0; i < (MAXDEG*2); i++) dst[i] = 0;
	
  for (i = 0; i < MAXDEG; i++) {
    for(j=MAXDEG; j<(MAXDEG*2); j++) tmp1[j]=0;
		
    for(j=0; j<MAXDEG; j++) tmp1[j]=gmult(pp2[j], pp1[i]);

    for (j = (MAXDEG*2)-1; j >= i; j--) tmp1[j] = tmp1[j-i];
    for (j = 0; j < i; j++) tmp1[j] = 0;
		

    for(j=0; j < (MAXDEG*2); j++) dst[j] ^= tmp1[j];
  }
}


	
void init_gamma (int gamma[])
{
  int e, tmp[MAXDEG];
	
  zero_poly(gamma);
  zero_poly(tmp);
  gamma[0] = 1;
	
  for (e = 0; e < NErasures; e++) 
  {
    copy_poly(tmp, gamma);
    scale_poly(gexp[ErasureLocs[e]], tmp);
    mul_z_poly(tmp);
    add_polys(gamma, tmp);
  }
}
	
	
	

int compute_discrepancy (int lambda[], int S[], int L, int n)
{
  int i, sum=0;
	
  for (i = 0; i <= L; i++) 
    sum ^= gmult(lambda[i], S[n-i]);
  return (sum);
}


void add_polys (int dst[], int src[]) 
{
  int i;
  for (i = 0; i < MAXDEG; i++) dst[i] ^= src[i];
}

void copy_poly (int dst[], int src[]) 
{
  int i;
  for (i = 0; i < MAXDEG; i++) dst[i] = src[i];
}

void scale_poly (int k, int poly[]) 
{	
  int i;
  for (i = 0; i < MAXDEG; i++) poly[i] = gmult(k, poly[i]);
}


void zero_poly (int poly[]) 
{
  int i;
  for (i = 0; i < MAXDEG; i++) poly[i] = 0;
}


static void mul_z_poly (int src[])
{
  int i;
  for (i = MAXDEG-1; i > 0; i--) src[i] = src[i-1];
  src[0] = 0;
}




void  Find_Roots (int NPAR)
{
  int sum, r, k;	
  NErrors = 0;
  
  for (r = 1; r < 256; r++)
   {
    sum = 0;

    for (k = 0; k < NPAR+1; k++) 
    {
      sum ^= gmult(gexp[(k*r)%255], Lambda[k]);
    }
    if (sum == 0) 
      { 
	     ErrorLocs[NErrors] = (255-r); NErrors++; 
	     //fprintf(stderr, "Root found at r = %d, (255-r) = %d\n", r, (255-r));
      }
  }
}


int correct_errors_erasures (unsigned char *buffer, int nbytes,int csize,int nerasures, int erasures[],int NPAR,unsigned char *parity_point)
{
  int r, i, j, err;

  unsigned char codeword[150];

  for(i=0;i<nbytes;i++)
  {
    codeword[i]=*(buffer+i);
  }
  for(i=0;i<NPAR;i++)
  {
    codeword[i+nbytes]=*(parity_point+i);
  }

  NErasures = nerasures;
  for (i = 0; i < NErasures; i++) ErasureLocs[i] = erasures[i];

  Modified_Berlekamp_Massey(NPAR);
  Find_Roots(NPAR);
  

  if ((NErrors <= NPAR) && NErrors > 0) 
  { 

    for (r = 0; r < NErrors; r++) 
	  {
        if (ErrorLocs[r] >= csize) 
	    {     
    	 return -1;
      }
     }

    for (r = 0; r < NErrors; r++) 
	 {
      int num, denom;
      i = ErrorLocs[r];

      num = 0;
      for (j = 0; j < 2*NPAR; j++) 
	       num ^= gmult(Omega[j], gexp[((255-i)*j)%255]);
      
      denom = 0;
      for (j = 1; j < 2*NPAR; j += 2) 
	   {
        	denom ^= gmult(Lambda[j], gexp[((255-i)*(j-1)) % 255]);
      }
      
      err = gmult(num, ginv(denom));
             
	  if((csize-i-1)<nbytes)   *(buffer+(csize-i-1))^= err;
	  else *(parity_point+(csize-i-1-nbytes))^= err;
    }

    return 1;
  }
  else 
  {
	if ( NErrors)
		return -1;
	else
		return 1;
    
  }
}