% File: egcodepol.m
% Description:
% Determine the generator polynomial of a binary EG code
%
% Copyright (c) 2006. Robert Morelos-Zaragoza. All rights reserved.

clear all
warning off all

OK = 0;
while ~OK
    mu = input('Enter the value of mu: ');
    s = input('Enter the value of s: ');
    m = input('Enter the value of m: ');

    ms = m*s;
    ms2 = 2^ms;
    s2 = 2^s;

    % Generate GF(2^ms)
    if (ms == 2)
        pol = [1 1 1];              % Primitve polynomial x^2+x+1
        pold = 7;                   % In decimal
        OK = 1;
    elseif (ms == 3)
        pol = [1 0 1 1];            % Primitve polynomial x^3+x+1
        pold = 11;                  % In decimal
        OK = 1;
    elseif (ms == 4)
        pol = [1 0 0 1 1];          % Primitve polynomial x^4+x+1
        pold = 19;                  % In decimal
        OK = 1;
    elseif (ms == 5)
        pol = [1 0 0 1 0 1];        % Primitve polynomial x^5+x^2+1
        pold = 37;                  % In decimal
        OK = 1;
    elseif (ms == 6)
        pol = [1 0 0 0 0 1 1];      % Primitve polynomial x^6+x+1
        pold = 67;                  % In decimal
        OK = 1;
    elseif (ms == 7)
        pol = [1 0 0 0 0 0 1 1];    % Primitve polynomial x^7+x+1
        pold = 131;                 % In decimal
        OK = 1;
    elseif (ms == 8)
        pol = [1 0 0 0 1 1 1 0 1];  % Primitve polynomial x^8+x^4+x^3+x^2+1
        pold = 285;                 % In decimal
        OK = 1;
    elseif (ms == 9)
        pol = [1 0 0 0 0 1 0 0 0 1];% Primitve polynomial x^9+x^4+1
        pold = 529;                  % In decimal
        OK = 1;
    elseif (ms == 10)
        pol = [1 0 0 0 0 0 0 1 0 0 1];% Primitve polynomial x^10+x^3+1
        pold = 1033;                     % In decimal
        OK = 1;
    else
        fprintf('Wrong paraneters!\n\n');
    end
end

fprintf('Determining the generator polynomial of a (%d,%d)-order', mu,s);
fprintf(' EG code of length %d\n', ms2-1);

alpha = gf(2,ms,pold)                   % Primitive element in GF(2^m)

% Construct the set of zeros powers by checking condition
for h=1:ms2-2
    max(h) = -inf;
    for l=0:s-1
        vec = [];
        hl = mod(h*(2^l),ms2-1);
        % Compute the 2^s-weight of hl
        hnew = hl;
        for i=1:m
            vec(i) = mod(hnew,s2);
            hnew = floor(hnew/s2);
        end
        if sum(vec(:)) > max(h)
            max(h) = sum(vec(:));
        end
    end
end
zero = [];
limit = (m-mu-1)*(s2-1);
for h=1:ms2-2
    if ((max(h) > 0) && (max(h) <= limit))
        zero = [ zero h];
    end
end

fprintf('Set of zeros powers: S = { %d', zero(1));
for i=2:length(zero)
    fprintf(', %d', zero(i));
end
fprintf('} \n');

% Compute the generator polynomial of the code
g = [1 alpha^zero(1)];
for i=2:length(zero)
    g = conv(g,[1 alpha^zero(i)]);
end

fprintf('Generator polynomial: g(x) = 1');      % g(0) = 1 always
for i=length(zero):-1:1
    if g(i) ~= 0
        fprintf(' + x^%d',length(zero)-i+1);
    end
end
fprintf('\n');
fprintf('This is a binary cyclic (%d,%d) EG code.\n', ms2-1, ms2-1-length(zero));