Berlekamp Massey Algorithm for BCH simplified binary version

Question

I am trying to follow Lin, Costello's explanation of the simplified BM algorithm for the binary case in page 210 of chapter 6 with no success on finding the error locator polynomial. I'm trying to implement it in MATLAB like this:

function [locator_polynom] = compute_error_locator(syndrome, t, m, field, alpha_powers)
% 
    % Initial conditions for the BM algorithm
    polynom_length = 2*t;
    syndrome = [syndrome; zeros(3, 1)];
    % Delta matrix storing the powers of alpha in the corresponding place
    delta_rho     = uint32(zeros(polynom_length, 1)); delta_rho(1)=1;
    delta_next    = uint32(zeros(polynom_length, 1));
   
    % Premilimnary values
    n_max = uint32(2^m - 1);

    % Initialize step mu = 1
    delta_next(1) = 1; delta_next(2) = syndrome(1); % 1 + S1*X
    % The discrepancy is stored in polynomial representation as uint32 numbers
    value = gf_mul_elements(delta_next(2), syndrome(2), field, alpha_powers, n_max);
    discrepancy_next = bitxor(syndrome(3), value);
    % The degree of the locator polynomial
    locator_degree_rho  = 0;
    locator_degree_next = 1;
    
    % Update all values
    locator_polynom     = delta_next;
    delta_current       = delta_next;
    discrepancy_rho     = syndrome(1);
    discrepancy_current = discrepancy_next;
    locator_degree_current     = locator_degree_next;
    rho = 0; % The row with the maximum value of 2mu - l starts at 1
    
    for i = 1:t % Only the even steps are needed (so make t out of 2*t)
        if discrepancy_current ~= 0
            % Compute the correction factor 
            correction_factor = uint32(zeros(polynom_length, 1));
            x_exponent = 2*(i - rho);
            if (discrepancy_current == 1 || discrepancy_rho == 1)
                d_mu_times_rho = discrepancy_current * discrepancy_rho;
            else
                alpha_discrepancy_mu        = alpha_powers(discrepancy_current);
                alpha_discrepancy_rho       = alpha_powers(discrepancy_rho);
                alpha_inver_discrepancy_rho = n_max - alpha_discrepancy_rho;
                % The alpha power for dmu * drho^-1 is
                alpha_d_mu_times_rho  = alpha_discrepancy_mu + alpha_inver_discrepancy_rho;
                % Equivalent to aux mod(2^m - 1)
                alpha_d_mu_times_rho  = alpha_d_mu_times_rho - ...
                                  n_max * uint32(alpha_d_mu_times_rho > n_max);
                d_mu_times_rho = field(alpha_d_mu_times_rho);          
            end
            correction_factor(x_exponent+1) = d_mu_times_rho;
            correction_factor = gf_mul_polynoms(correction_factor,...
                                                delta_rho,...
                                                field, alpha_powers, n_max);
            % Finally we add the correction factor to get the new delta
            delta_next = bitxor(delta_current, correction_factor(1:polynom_length));
            
            % Update used data
            l = polynom_length;
            while delta_next(l) == 0 && l>0
                l = l - 1;
            end
            locator_degree_next = l-1;
            % Update previous maximum if the degree is higher than recorded
            if (2*i - locator_degree_current) > (2*rho - locator_degree_rho)
                locator_degree_rho  = locator_degree_current;
                delta_rho           = delta_current;
                discrepancy_rho     = discrepancy_current;
                rho = i;
            end
        else
            % If the discrepancy is 0, the locator polynomial for this step
            % is passed to the next one. It satifies all newtons' equations
            % until now.
            delta_next = delta_current; 
        end

        % Compute the discrepancy for the next step
        syndrome_start_index = 2 * i + 3;
        discrepancy_next     = syndrome(syndrome_start_index); % First value
        for k = 1:locator_degree_next
            value = gf_mul_elements(delta_next(k + 1), ...
                                    syndrome(syndrome_start_index - k), ...
                                    field, alpha_powers, n_max);
            discrepancy_next = bitxor(discrepancy_next, value);
        end
        
        % Update all values
        locator_polynom = delta_next;
        delta_current = delta_next;
        discrepancy_current = discrepancy_next;
        locator_degree_current = locator_degree_next;
    end
end

I'm trying to see what's wrong but I can't. It works for the examples in the book, but not always. As an aside, to compute the discrepancy S_2mu+3 is needed, but when I have only 24 syndrome coefficients how is it computed on step 11 where 2*11 + 3 is 25? Thanks in advance!

Answer 1

It turns out the code is ok. I made a different implementation from Error Correction and Coding. Mathematical Methods and gives the same result. My problem is at the Chien Search.

Code for the interested:

function [c] = compute_error_locator_v2(syndrome, m, field, alpha_powers)
% 
    % Initial conditions for the BM algorithm
    % Premilimnary values
    N = length(syndrome);
    n_max = uint32(2^m - 1);
    polynom_length = N/2 + 1;
    L = 0; % The curent length of the LFSR
    % The current connection polynomial
    c = uint32(zeros(polynom_length, 1)); c(1) = 1;
    % The connection polynomial before last length change
    p = uint32(zeros(polynom_length, 1)); p(1) = 1;
    l = 1; % l is k - m, the amount of shift in update
    dm = 1; % The previous discrepancy
    for k = 1:2:N % For k = 1 to N in steps of 2
        % ========= Compute discrepancy ==========
        d = syndrome(k);
        for i = 1:L
            aux = gf_mul_elements(c(i+1), syndrome(k-i), field, alpha_powers, n_max);
            d   = bitxor(d, aux);
        end
        
        if d == 0 % No change in polynomial
            l = l + 1;
        else
            % ======== Update c ========
            t = c;
            % Compute the correction factor 
            correction_factor = uint32(zeros(polynom_length, 1));
            % This is d *  dm^-1
            dd_sum = modulo(alpha_powers(d) + n_max - alpha_powers(dm), n_max);
            for i = 0:polynom_length - 1
                if p(i+1) ~= 0
                    % Here we compute d*d^-1*p(x_i)
                    ddp_sum = modulo(dd_sum + alpha_powers(p(i+1)), n_max);
                    if ddp_sum == 0
                        correction_factor(i + l + 1) = 1;
                    else
                        correction_factor(i + l + 1) = field(ddp_sum);
                    end
                end
            end
            % Finally we add the correction factor to get the new locator
            c = bitxor(c, correction_factor);
            
            if (2*L >= k) % No length change in update
                l = l + 1;
            else
                p   = t;
                L   = k - L;
                dm  = d;
                l = 1;
            end
        end
        l = l + 1;
    end
end

The code comes from this implementation of the Massey algorithm