close all; clear all;
%reed-solomon decoder simulation
%by:yumi


m=8;        %number of bits each symbol
n=2^m-1;    %code length
t=8;         %error correction ability
k=n-2*t;      %information length
pp = 301;    %primitive polynomial

msg = gf([1:k], m, pp);
code = rsenc(msg, n, k);
%code = gf(zeros(1, n), m, pp);


alpha = gf(2, m, pp);
zero = gf(0, m, pp);
one = gf(1, m, pp);

%reccode=zeros(1, n);
%reccode(7) = 5;
%reccode(5) = 2;
%reccode(2) = 11;

%generate integers between 
reccode=zeros(1, n);
reccode(1:8)=255;
%reccode(1)=255;
%reccode(2)=2;
%reccode(3)=3;
%reccode(15)=72;
%reccode(20)=57;
%reccode(32)=135;
%reccode(90)=172;
%reccode(156)=145;
%reccode(170)=74;
%reccode(255)=255;

reccode = gf(reccode, m, pp);

reccode = code + reccode;

%functional simulation of decoding
%creating alpha array
alpha_tb=gf(zeros(1, 2*t), m, pp);
for i=1:2*t;
    alpha_tb(i)=alpha^i;
end;

%syndrome
syndrome_o = gf(zeros(1, 2*t), m, pp);
syndrome = gf(zeros(1, 2*t), m, pp);
for i = 1:n,
    syndrome_o = syndrome_o.*alpha_tb +reccode(i);
    cat = 1;
end;
%reverse s(2t-1)...s(0)
for i=1:2*t,
    syndrome(i)=syndrome_o(2*t+1-i);
end;

%ibma
lamda = gf([1, zeros(1, t)], m, pp);
lamda0 = gf([1, zeros(1, t)], m, pp);
b = gf([0, 1, zeros(1, t)], m, pp);
k = 0;
gamma = one;
delta = zero;
syndrome_array1 = gf(zeros(1, t+1), m, pp);

for r = 1:2*t,
    %step 1:
    r1 = 2*t-r+1;
    r2 = min(r1+t, 2*t);
    num = r2-r1+1;
    syndrome_array1(1:num) = syndrome(r1:r2);
    delta = syndrome_array1*lamda';
    
    %step 2:
    lamda0 = lamda;
    lamda = gamma*lamda-delta*b(1:t+1);
    
    %step 3:
    if((delta ~= zero) && (k>=0))
        b(2:2+t)=lamda0;
        gamma = delta;
        k=-k-1;
    else
        b(2:2+t) = b(1:t+1);
        gamma = gamma;
        k=k+1;
    end
    
   mini=1;
end
 
omega = gf(zeros(1, t), m, pp);
syndrome_array2 = gf(zeros(1, t), m, pp);
for i= 1:t,
    i1 = 2*t-i+1;
    syndrome_array2(1:i) = syndrome(i1: 2*t);
    omega(i) = syndrome_array2*lamda(1:t)';    
end;


%chien search algorithm and forney algorithm in paralle
%seperate the even and the odd parts of lamda and omega
even=floor(t/2)*2;
if(even==t)
    odd=t-1;
else
    odd=t;
end
%lamda_even=lamda(1:2:even+1);
%lamda_odd=lamda(2:2:odd+1);

%inverstable
inverse_tb = gf(zeros(1, t+1), m, pp);
for i=1:t+1,
    inverse_tb(i) = alpha^(-i+1);
end;
%inverse_omega=inverse_tb(1:t);

%iteratively compute omega and 
lamda_v=gf(0, m, pp);
lamda_ov=gf(0, m, pp);
omega_v=gf(0, m, pp);
accu_tb=gf(ones(1, t+1), m,pp);
accu_tb1=gf(ones(1, t), m, pp);
for i=1:n,
    root=alpha^(1-i);
    lamda_v=lamda*accu_tb';
    lamda_ov=lamda(2:2:odd+1)*accu_tb(2:2:odd+1)';
    omega_v=omega*accu_tb1';
    
    accu_tb=accu_tb.*inverse_tb;
    accu_tb1=accu_tb1.*inverse_tb(1:t);
    
    if(lamda_v==zero)
        error(1,n-i+1)=1;
        ev=(omega_v/lamda_ov)*root;
        evv=ev.x;
        error(2, n-i+1)=double(evv);
    else
        error(1, n-i+1)=0;
        error(2, n-i+1)=0;
    end
    kk=1;
end


found = find(error(1,:)~=0)
error(2, found)