1## Copyright (C) 2003 David Bateman
2##
3## This program is free software; you can redistribute it and/or modify it under
4## the terms of the GNU General Public License as published by the Free Software
5## Foundation; either version 3 of the License, or (at your option) any later
6## version.
7##
8## This program is distributed in the hope that it will be useful, but WITHOUT
9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
11## details.
12##
13## You should have received a copy of the GNU General Public License along with
14## this program; if not, see <http://www.gnu.org/licenses/>.
15
16## -*- texinfo -*-
17## @deftypefn  {Function File} {@var{p} =} bchpoly ()
18## @deftypefnx {Function File} {@var{p} =} bchpoly (@var{n})
19## @deftypefnx {Function File} {@var{p} =} bchpoly (@var{n}, @var{k})
20## @deftypefnx {Function File} {@var{p} =} bchpoly (@var{prim}, @var{k})
21## @deftypefnx {Function File} {@var{p} =} bchpoly (@var{n}, @var{k}, @var{prim})
22## @deftypefnx {Function File} {@var{p} =} bchpoly (@dots{}, @var{probe})
23## @deftypefnx {Function File} {[@var{p}, @var{f}] =} bchpoly (@dots{})
24## @deftypefnx {Function File} {[@var{p}, @var{f}, @var{c}] =} bchpoly (@dots{})
25## @deftypefnx {Function File} {[@var{p}, @var{f}, @var{c}, @var{par}] =} bchpoly (@dots{})
26## @deftypefnx {Function File} {[@var{p}, @var{f}, @var{c}, @var{par}, @var{t}] =} bchpoly (@dots{})
27##
28## Calculates the generator polynomials for a BCH coder. Called with no input
29## arguments @code{bchpoly} returns a list of all of the valid BCH codes for
30## the codeword length 7, 15, 31, 63, 127, 255 and 511. A three column matrix
31## is returned with each row representing a separate valid BCH code. The first
32## column is the codeword length, the second the message length and the third
33## the error correction capability of the code.
34##
35## Called with a single input argument, @code{bchpoly} returns the valid BCH
36## codes for the specified codeword length @var{n}. The output format is the
37## same as above.
38##
39## When called with two or more arguments, @code{bchpoly} calculates the
40## generator polynomial of a particular BCH code. The generator polynomial
41## is returned in @var{p} as a vector representation of a polynomial in
42## GF(2). The terms of the polynomial are listed least-significant term
43## first.
44##
45## The desired BCH code can be specified by its codeword length @var{n}
46## and its message length @var{k}. Alternatively, the primitive polynomial
47## over which to calculate the polynomial can be specified as @var{prim}.
48## If a vector representation of the primitive polynomial is given, then
49## @var{prim} can be specified as the first argument of two arguments,
50## or as the third argument. However, if an integer representation of the
51## primitive polynomial is used, then the primitive polynomial must be
52## specified as the third argument.
53##
54## When called with two or more arguments, @code{bchpoly} can also return the
55## factors @var{f} of the generator polynomial @var{p}, the cyclotomic coset
56## for the Galois field over which the BCH code is calculated, the parity
57## check matrix @var{par} and the error correction capability @var{t}. It
58## should be noted that the parity check matrix is calculated with
59## @code{cyclgen} and limitations in this function means that the parity
60## check matrix is only available for codeword length up to 63. For
61## codeword length longer than this @var{par} returns an empty matrix.
62##
63## With a string argument @var{probe} defined, the action of @code{bchpoly}
64## is to calculate the error correcting capability of the BCH code defined
65## by @var{n}, @var{k} and @var{prim} and return it in @var{p}. This is
66## similar to a call to @code{bchpoly} with zero or one argument, except that
67## only a single code is checked. Any string value for @var{probe} will
68## force this action.
69##
70## In general the codeword length @var{n} can be expressed as
71## @code{2^@var{m}-1}, where @var{m} is an integer. However, if
72## [@var{n},@var{k}] is a valid BCH code, then a shortened BCH code of
73## the form [@var{n}-@var{x},@var{k}-@var{x}] can be created with the
74## same generator polynomial
75##
76## @seealso{cyclpoly, encode, decode, cosets}
77## @end deftypefn
78
79function [p, f, c, par, t] = bchpoly (nn, k, varargin)
80
81  if (nargin < 0 || nargin > 4)
82    print_usage ();
83  endif
84
85  probe = 0;
86  prim = 0;    ## Set to zero to use default primitive polynomial
87  if (nargin == 0)
88    m = [3:9];
89    n = 2.^m - 1;
90    nn = n;
91  elseif (isscalar (nn))
92    m = ceil (log2 (nn+1));
93    n = 2.^m - 1;
94    if (! (n == fix (n) && n >= 7 && m == fix (m)))
95      error ("bchpoly: N must be a integer greater than 3");
96    endif
97  else
98    prim = bi2de (n);
99    if (!isprimitive (prim))
100      error ("bchpoly: PRIM must be a primitive polynomial of GF(2^M)");
101    endif
102    m = length (n) - 1;
103    n = 2^m - 1;
104  endif
105
106  if (nargin > 1 && ! (isscalar (k) && k == fix (k) && k <= n))
107    error ("bchpoly: K must be an integer less than N");
108  endif
109
110  for i = 1:length (varargin)
111    arg = varargin{i};
112    if (ischar (arg))
113      probe = 1;
114      if (nargout > 1)
115        error ("bchpoly: only one output argument allowed when probing valid codes");
116      endif
117    else
118      if (prim != 0)
119        error ("bchpoly: primitive polynomial already defined");
120      endif
121      prim = arg;
122      if (!isscalar (prim))
123        prim = bi2de (prim);
124      endif
125      if (! (prim == fix (prim) && prim >= 2^m && prim <= 2^(m+1)
126             && isprimitive (prim)))
127        error ("bchpoly: PRIM must be a primitive polynomial of GF(2^M)");
128      endif
129    endif
130  endfor
131
132  ## Am I using the right algo to calculate the correction capability?
133  if (nargin < 2)
134    if (nargout > 1)
135      error ("bchpoly: only one output argument allowed when probing valid codes");
136    endif
137
138    p = [];
139    for ni = 1:length (n)
140      c = cosets (m(ni), prim);
141      nc = length (c);
142      fc = zeros (1, nc);
143      f = [];
144
145      for t = 1:floor (n(ni)/2)
146        for i = 1:nc
147          if (fc(i) != 1)
148            cl = log (c{i});
149            for j = 2*(t-1)+1:2*t
150              if (find (cl == j))
151                f = [f, c{i}.x];
152                fc(i) = 1;
153                break;
154              endif
155            endfor
156          endif
157        endfor
158
159        k = nn(ni) - length (f);
160        if (k < 2)
161          break;
162        endif
163
164        if (!isempty (p) && (k == p(size (p, 1),2)))
165          p(size (p, 1),:) = [nn(ni), k, t];
166        else
167          p = [p; [nn(ni), k, t]];
168        endif
169      endfor
170    endfor
171  else
172    c = cosets (m, prim);
173    nc = length (c);
174    fc = zeros (1, nc);
175    f = [];
176    fl = 0;
177    f0 = [];
178    f1 = [];
179    t = 0;
180    do
181      t++;
182      f0 = f1;
183      for i = 1:nc
184        if (fc(i) != 1)
185          cl = log (c{i});
186          for j = 2*(t-1)+1:2*t
187            if (find (cl == j))
188              f1 = [f1, c{i}.x];
189              fc(i) = 1;
190              ptmp = gf ([c{i}(1), 1], m, prim);
191              for l = 2:length (c{i})
192                ptmp = conv (ptmp, [c{i}(l), 1]);
193              endfor
194              f = [f; [ptmp.x, zeros(1, m - length (ptmp) + 1)]];
195              fl = fl + length (ptmp);
196              break;
197            endif
198          endfor
199        endif
200      endfor
201    until (length (f1) > nn - k)
202    t--;
203
204    if (nn - length (f0) != k)
205      error ("bchpoly: could not find valid generator polynomial for parameters");
206    endif
207
208    if (probe)
209      p = [nn, k, t];
210    else
211
212      ## Have to delete a line from the list of minimum polynomials
213      ## since we've gone one past in calculating f1 above to be
214      ## sure or the error correcting capability
215      f = f(1:size (f, 1) - 1,:);
216
217      p = gf ([f0(1), 1], m, prim);
218      for i = 2:length (f0)
219        p = conv (p, [f0(i), 1]);
220      endfor
221      p = p.x;
222
223      if (nargout > 3)
224        if (n > 64)
225          warning ("bchpoly: could not create parity matrix");
226          par = [];
227        else
228          par = cyclgen (n, p);
229        endif
230      endif
231    endif
232  endif
233
234endfunction
235
236%% Test input validation
237%!error bchpoly (1)
238%!error bchpoly (1, 2, 3, 4, 5)
239%!error bchpoly (5, 10)
240