function [GENPOLY, FACTORS, cs, h, t] = bchpoly(n, k)
%BCHPOLY
%
%WARNING: This is an obsolete function and may be removed in the future.
%         Please use BCHGENPOLY instead.


%BCHPOLY Produce parameters or generator polynomial for binary BCH code.
%   BCHPOLY lists the codeword length, message length, and error-correction
%   capability for BCH code in a figure. The codeword lengths listed are 7, 15,
%   31, 63, 127, 255, and 511.
%
%   GENPOLY = BCHPOLY outputs the list to a vector GENPOLY. The first column is the code
%   word length, the second column is the message length, and the third column
%   is the error-correction capability.
%
%   GENPOLY = BCHPOLY(N) outputs the code word length equal or larger than N, its
%   message lengths, and error-correction capability on the first, second and
%   third column of GENPOLY, respectively. When N > 1023, the error-correction
%   capability is not listed.
%
%   GENPOLY = BCHPOLY(N, K) produces BCH code generator polynomial GENPOLY based on two
%   positive numbers: the codeword length N and the message length K. The
%   codeword length must be 2^M - 1 for some integer M greater than 2. K
%   should be a valid message length, which can be found by using
%   GENPOLY = BCHPOLY(N) to get a list.  When N is a binary vector instead of a
%   positive number, N is a primitive polynomial generating the field GF(2^DEG),
%   where DEG is the degree of polynomial N(X). The codeword length in this case
%   should be 2^DEG - 1.
%
%   [GENPOLY, FACTORS] = BCHPOLY(N, K) outputs minimal polynomials FACTORS in the field
%   GF(2^M). Each row of FACTORS is a minimal polynomial.
%
%   [GENPOLY, FACTORS, CST] = BCHPOLY(N, K) also lists the cosets CST that are the roots of the
%   polynomial GENPOLY. The cosets are in the same format as the output from the
%   GFCOSETS function.
%
%   [GENPOLY, FACTORS, CST, H] = BCHPOLY(N, K) outputs the parity-check matrix H.
%
%   [GENPOLY, FACTORS, CST, H, T] = BCHPOLY(N, K) outputs the error-correction capability
%   of the BCH code.
%
%   See also CYCLGEN, GFMINPOL.

%   Copyright 1996-2002 The MathWorks, Inc.
%   $Revision: 1.1.6.1 $

%   Reference: Peterson, W. Wesley and E. J. Weldon, Jr. Error-correcting
%   Codes, 2nd ed. Cambridge, Mass.: MIT Press, 1972.

% parameter assignment
n_lst = [3    7   1;
         4   15   3;
         5   31   5;
         6   63  11;
         7  127  17;
         8  255  33;
         9  511  57;
        10 1023 105];

% known error-correction capability
t_lst = [1,...
         [1:3],...
         [1:3],5 7,...
         [1:7],10 11 13 15,...
         [1:7], 9,10 11 13 14 15 21 23 27 31,...
         [1:15], 18 19 [21:23],[25:27],[29:31],42 43 45 47 55 59 63,...
         [1:16],[18:23],[25:31],[36:39],[41:43],[45:47],51,[53:55],58 59,[61:63],...
            85 87 91 93 95 109 111 119 121,...
         [1:31],[34:39],[41:47],[49:55],[57:63],[73:75],[77:79],82,83,85+[0:2],...
         89+[0:2],93+[0:2],102,103,106,107,109,110,111,115,117+[0:2],122,123,...
        125+[0:2],170,171,173,175,181,183,187,189,191,219,223,239,247,255];

% give a list of all BCH code
% code word length / message length / error-correction capability
% When there is no output parameter, draw a figure.
% When there is output parameter, output the line

if nargin < 1
    if nargout < 1
       % title of the figure.
        yy = 'BCH Code Generated by Default Primitive Polynomial';

        % In case the figure exist
        hand = get(0, 'Child');
        if ~isempty(hand)
            for i = 1:length(hand)
                if strcmp(yy, get(hand(i), 'Name'))
                    set(0,'Current',hand(i));
                    if length(get(get(hand(i),'Child'),'Child')) == 278
                        set(hand(i), 'Position', [10 10 700 700]);
                        return;
                    else
                        close(hand(i));
                    end
                end;
            end;
        end;

        % Create figure
        h = figure('Position',    [10 10 700 700],...
                   'NumberTitle', 'off',....
                   'Clipping',    'off',...
                   'Name',        yy);
        hght = 43;
        x = axes('Xlim', [0 9], 'Ylim', [0 hght], 'Visible','off',...
           'Clipping','off');
        set(x, 'visible','off');
        text(0, -2, 'N: code word length; K: message length; T: error-correction capability');
        for i=1:3
            tmp = text(3*(i-1),  hght+1,'   N');
            tmp = text(3*(i-1)+1,hght+1,'   K');
	    set(tmp, 'color', [0 .5 .5]);
            tmp = text(3*(i-1)+2,hght+1,'   T');
	    set(tmp, 'color', [1 0 0]);
        end;
        XData=[0 9;0 9;0 9; 2.5 2.5;5.5 5.5];
        YData=[hght hght;[hght hght]+2;-1 -1;-1 hght+2;-1 hght+2];
        for inlp=1:size(XData,1),
       	  line('Parent',x            , ...
               'XData' ,XData(inlp,:), ...
               'YData' ,YData(inlp,:), ...
               'Color' ,'b'         , ...
               'LineStyle','-');
        end % for
    end;

    % list all numbers in the figure or output to variable GENPOLY.
    next_coset = 1;
    n_m = 0;
    n_base = 0;

    if ~(nargout == 0),
       GENPOLY = [];
    end

    for i = 0 : 126
       if (nargout < 1),
          y_p = rem(i,hght);
          x_p = floor(i / hght) * 3;

          if ((next_coset) | (y_p == 0)),
            drawnow;
            if next_coset
                xx = num2str(n_lst(n_m+1, 2));
            else
                xx = num2str(n_lst(n_m, 2));
            end;
            xx = [char(ones(1, 4-length(xx))*32) xx];
            xx=text(x_p, hght-y_p-1, xx);
            set(xx,'Color',[0 0 0]);
          end
        end
        if next_coset
            n_m = n_m + 1;
            cs = gfcosets(n_lst(n_m, 1));
            next_coset = 0;
        end;
        if (nargout < 1)
           xx = num2str(n_lst(n_m, 2) - sum(sum(~isnan(cs(2:2+i-n_base,:)))));
           xx = [char(ones(1, 4-length(xx))*32) xx];
           xx = text(x_p+1, hght-y_p-1, xx);
           set(xx,'Color',[0 .5 .5]);
           xx = num2str(t_lst(i+1));
           xx = [char(ones(1, 4-length(xx))*32) xx];
           xx = text(x_p+2, hght-y_p-1, xx);
           set(xx,'Color',[1 0 0]);
        else
           GENPOLY = [GENPOLY;
                [n_lst(n_m, 2),...
                    n_lst(n_m, 2)-sum(sum(~isnan(cs(2:2+i-n_base,:)))),...
                    t_lst(i+1)]];
        end;
        if n_base + n_lst(n_m, 3) == i + 1
            next_coset = 1;
            n_base = n_base + n_lst(n_m, 3);
        end;
    end;
    if (nargout < 1),
       set(h, 'HandleVisibility','callback');
    end

    return;
end;

% When specific code word length is specified, output a list of the
% code word length / message length / error-correction capability

if length(n) > 1
    m = gftrunc(n);
    n = 2^length(m) - 1;
    df_flag = 0;
else
    df_flag = 1;
end;
indx = find(n <= n_lst(:, 2));
if isempty(indx)
    % beyond the range
    dim = 0;
    k = 11;
    while (dim == 0)
        if (2^k-1 >= n)
            dim = k;
        else
            k = k + 1;
        end;
    end;
else
    dim = n_lst(indx(1), 1);
end;
n = 2^dim - 1;


if nargin < 2
   GENPOLY = [];
   cs = gfcosets(dim);
   [n_cs, m_cs] = size(cs);
   if nargout < 1
      yy = ['BCH Code Generated by ', num2str(dim), 'th Order Primitive Polynomial'];

      % in case the figure exists
      hand = get(0, 'Child');
      if ~isempty(hand)
         for i = hand(:)'
            if strcmp(yy, get(i, 'Name'))
               set(0, 'Current', i);
               sz = get(i, 'UserData');
               if length(get(get(i, 'Child'), 'Child')) == sz(1)
                  set(i, 'Position', sz(2:5));
                        return;
                     else
                        close(i);
                     end;
                  end;
               end;
            end;
        % figure frame construction

        x_a = 1 + (dim>=8) + (dim>10);
        if x_a < 3
            hght = ceil((n_cs - 2) / x_a);
        else
            hght = ceil((n_cs - 2) / x_a / 2);
        end;
        h = figure('Position', [10, 10, 240*x_a, hght*12+100],...
                   'NumberTitle', 'off',...
                   'Clipping',    'off',...
                   'Name',        yy);
        x = axes('Xlim', [0 2*x_a], 'Ylim', [-1 hght+2], 'Visible','off',...
            'Clipping','off');
        set(x, 'visible','off');
        if n <= 1023
            text(0, -1, 'K: message length; T: error-correction capability.');
        else
            text(0, -1, 'K: message length.');
        end;

        for i=1 : x_a
            tmp = text(2*(i-1),hght+1,'   K');
	    set(tmp, 'color', [0 .5 .5]);
            if (dim <= 10)
                tmp = text(2*(i-1)+1,hght+1,'   T');
                set(tmp, 'color', [1 0 0]);
            else
                tmp = text(2*(i-1)+1,hght+1,'   K');
                set(tmp, 'color', [0 .5 .5]);
            end;
        end;
        text(0, hght+2, ['Code word length', num2str(n)]);

        plot([0, 2*x_a],   [hght, hght]+.5,  'b-');
        plot([0, 2*x_a],   [hght, hght]+2.5,'b-');
        plot([0, 2*x_a],   [-.5 -.5],       'b-');

        if x_a > 1
            plot([1.5,1.5],[-1,hght+2],   'b-');
        end;
        if x_a > 2
            plot([3.5,3.5],[-1,hght+2],   'b-');
        end;

        lengt = length(char(n));
    end;

    ind = find(n_lst(:,1) < dim);
    if ~isempty(ind)
        base_i = sum(n_lst(ind, 3));
    else
        base_i = 0;
    end;

    for i = 1:n_cs-2
        if (nargout < 1)
             xx = num2str(n - sum( sum( ~isnan( cs(1 : i+1, :))))+1);
            xx = [char(ones(1, lengt-length(xx))*32) xx];
            if (x_a <= 1) | ((x_a <= 2) & (i <= hght)) | ((x_a <= 3) & (i<=hght))
                % the first column
                xx=text(0, hght-i, xx);
            elseif (x_a<=2) | ((x_a<=3) & (i<=2*hght))
                % the second column
                xx=text(4-x_a, 2*hght-i, xx);
            else
                if (i <= 3*hght)
                    % the third column
                    xx=text(2, 3*hght-i, xx);
                elseif (i <= 4*hght)
                    % the fourth column
                    xx=text(3, 4*hght-i, xx);
                elseif (i <= 5*hght)
                    % the fifth column
                    xx=text(4, 5*hght-i, xx);
                else
                    % the sixth column
                    xx=text(5, 6*hght-i, xx);
                end;
            end;
            set(xx,'Color',[0 0.5 .5]);
            if dim <= 10
                xx = num2str(t_lst(i+base_i));
                xx = [char(ones(1, 4-length(xx))*32) xx];
                if (x_a<=1) | ((x_a<=2) & (i*2<=n_cs)) | ((x_a<=3) & (i*3<=n_cs))
                    % the first column
                    xx=text(1, hght-i, xx);
                elseif (x_a<=2) | ((x_a<=3) & (i*3<=2*n_cs))
                    % the second column
                    xx=text(3, 2*hght-i, xx);
                else
                    % the third column
                    xx=text(5, 3*hght-i, xx);
                end;
                set(xx,'Color',[1 0 0]);
            end;
            usda = [length(get(x, 'child')), [10, 10, 240*x_a, hght*12+100]];
            set(h, 'UserData', usda);
        else
            if dim <= 10
                GENPOLY = [GENPOLY;
                  [n, ...
                   n - sum( sum( ~isnan( cs(1 : i+1, :))))+1,...
                   t_lst(i+base_i)]];
            else
                GENPOLY = [GENPOLY;
                  [n, ...
                   n - sum( sum( ~isnan( cs(1 : i+1, :))))+1]];
            end;
        end;

    end;
    return;
end;

if (dim<1) | (k<1) | (floor(dim) ~= dim) | (floor(k) ~= k)
    error('Input variable for BCHPOLY is not valid');
end

if df_flag
    m = gfprimdf(dim);
end;

% compute the minimal polynomial
if k == n - dim
    % It is a trivial case, so simply assign the primitive polynomial
    % as BCH polynomial.
    FACTORS = m;
    GENPOLY = m;
    t = 1;
    if (nargout > 2)
        cs = gfcosets(dim);
    end;
    n_i = 1;
else
    % list cosets and take off the first row.
    cs = gfcosets(dim);
    [n_cs, m_cs] = size(cs);

    pl = gfminpol(cs(2:n_cs,1), m);

    n_terms = 0;
    n_i = 0;
    FACTORS = [];
    while ((n_terms < (n - k)) & (n_i+1 < n_cs))
        n_i = n_i + 1;
        n_terms = n_terms + sum(~isnan( cs(n_i+1, :)));
        FACTORS = [FACTORS; pl(n_i, :)];
    end;

    if ((n_i+1 >= n_cs) & (n_terms < (n - k)))
        warning(['The maximum message length is ', num2str(n-n_terms)]);
    end;

    GENPOLY = gftrunc(FACTORS(1,:));
    [n_FACTORS, m_FACTORS] = size(FACTORS);
    if n_FACTORS > 1
        for i = 2 : n_FACTORS
            GENPOLY = gfconv(GENPOLY, gftrunc(FACTORS(i,:)));
        end;
    end;
end;

% construct the parity-check matrix
if nargout > 3
    tp = gftuple([-1:n-1]', m);
    h = tp(2:n+1,:)';
    if n_i > 1
        for i = 2 : n_i

            % build the rows led by new coset.
            ad = [];
            seed = cs(i+1,1);
            for j = 0:n-1
                ad = [ad tp(mod(j*seed,n+1)+1,:)'];
            end;

            % remove unnecessary rows.
            for j = dim : -1 : 1
                % remove the all-zero rows.
                if max(ad(j, :)) == 0
                    ad(j,:) = [];
                else
                    % remove the repeated rows.
                    run_flag = 1;
                    jj = j - 1;
                    while run_flag & (jj > 0)
                        if max(rem(ad(j,:)+ad(jj, :), 2)) == 0
                            ad(jj,:)=[];
                            run_flag = 0;
                        end;
                        jj = jj - 1;
                    end;
                end;
            end;

            % add the new rows to the parity-check matrix.
            h = [h; ad];
        end;
    end;
end;

if (~df_flag) & (dim <= 10)
    if (max(rem(m - gfprimdf(dim), 2)) == 0)
        df_flag == 0;
    end;
end;

% detect error-correction capability.
if df_flag & (dim <= 10)
    % look up the data base
    ind = find(n_lst(:,1) < dim);
    if ~isempty(ind)
        base_i = sum(n_lst(ind, 3));
    else
        base_i = 0;
    end;
    t = t_lst(base_i + n_i);
else
    % designed distance. By Corollary 9.2 in Peterson & Weldon
    %    t = floor((n - k) / dim);
    t = n_i;
    %    dis_min = 2 * n_i + 1;  % the minimum by theory.
    tt = 1;
    for i = 1:t+1
        tt = tt * i;
    end;

    % by Farr Theorem (Theorem 9.6) in Peterson & Weldon, that's the solution.
    if (dim > 1 + log(tt) / log(2)) & (t < dim/2)
        return;
    end;

    if dim > 3
        disp('It takes a long time to process the minimum distance calculation');
    end;

    dis_min = 2 * t + 1;

    % The maximum distance will not be more than the following value
    % by Theory 9.8 of Peterson and Weldon.
    dis = 2 * dis_min-1;
    h = h';
    [x1, x2] = size(h);
    dis = min(dis, x2);
    ind = [0 1 zeros(1, x1 - 2)];
    i = 3;
    px1 = 2^x1 - 1;
    tt = floor((dis - 1) / 2);
    fprintf(['The error-correction capability ',num2str(t),' and ',num2str(tt),'.\n']);
    fprintf([num2str(px1), ' number to be processed. \nPlease wait....']);
    count = 0;
    while (i < px1) & (tt > t)
        n_j = 1;
        i = i + 1;
        n_k = 1;
        while n_j
            if ind(n_k)
                ind(n_k) = 0;
                n_k = n_k + 1;
            else
                ind(n_k) = 1;
                n_j = 0;
            end;
        end;
        sum_ind = sum(ind);
        if (sum_ind >= dis_min)
            if (sum_ind < dis)
                count = count + 1;
                if count > 1000
                    fprintf('.');
                    count = 0;
                end;
                idx = find(ind == 1);
                if max(rem(sum(h(idx, :)), 2)) == 0
                    dis = sum_ind;
                    tt = floor((dis-1)/2);
                    if tt ~= t
                        fprintf(['\nThe error-correction capability is, now, between ', num2str(t), ' and ',num2str(tt), '\n']);
                    end;
                end;
            end;
        end;
    end;
    fprintf('..done.\n');
    fprintf(['The error-correction capability is found to be ',num2str(tt), '\n']);
    t = floor((dis - 1) / 2);
    h = h';
end;

%--end of bchpoly.m--