Mercurial > mxe-octave
view src/of-communications-5-fixes.patch @ 4733:1ecb1e67eaa1
* src/src-msys-dos2unix.mk: update checksum
| author | John Donoghue <john.donoghue@ieee.org> |
|---|---|
| date | Sat, 09 Jun 2018 10:13:46 -0400 |
| parents | 00e61c4a5657 |
| children |
line wrap: on
line source
diff -uNr a/src/base-lu.cc b/src/base-lu.cc --- a/src/base-lu.cc 2018-04-09 13:25:42.884981069 -0400 +++ b/src/base-lu.cc 2018-04-09 13:53:52.501958203 -0400 @@ -125,7 +125,7 @@ for (octave_idx_type i = 0; i < a_nr; i++) pvt.xelem (i) = i; - for (octave_idx_type i = 0; i < ipvt.length (); i++) + for (octave_idx_type i = 0; i < ipvt.numel (); i++) { octave_idx_type k = ipvt.xelem (i); diff -uNr a/src/base-lu.cc~ b/src/base-lu.cc~ --- a/src/base-lu.cc~ 1969-12-31 19:00:00.000000000 -0500 +++ b/src/base-lu.cc~ 2018-04-09 13:26:27.478896947 -0400 @@ -0,0 +1,187 @@ +/* + +Copyright (C) 1996-2015 John W. Eaton +Copyright (C) 2009 VZLU Prague + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#include "base-lu.h" + +template <class lu_type> +base_lu<lu_type>::base_lu (const lu_type& l, const lu_type& u, + const PermMatrix& p) + : a_fact (u), l_fact (l), ipvt (p.transpose ().col_perm_vec ()) +{ + if (l.columns () != u.rows ()) + (*current_liboctave_error_handler) ("lu: dimension mismatch"); +} + +template <class lu_type> +bool +base_lu <lu_type> :: packed (void) const +{ + return l_fact.dims () == dim_vector (); +} + +template <class lu_type> +void +base_lu <lu_type> :: unpack (void) +{ + if (packed ()) + { + l_fact = L (); + a_fact = U (); // FIXME: sub-optimal + ipvt = getp (); + } +} + +template <class lu_type> +lu_type +base_lu <lu_type> :: L (void) const +{ + if (packed ()) + { + octave_idx_type a_nr = a_fact.rows (); + octave_idx_type a_nc = a_fact.cols (); + octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); + + lu_type l (a_nr, mn, lu_elt_type (0.0)); + + for (octave_idx_type i = 0; i < a_nr; i++) + { + if (i < a_nc) + l.xelem (i, i) = 1.0; + + for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++) + l.xelem (i, j) = a_fact.xelem (i, j); + } + + return l; + } + else + return l_fact; +} + +template <class lu_type> +lu_type +base_lu <lu_type> :: U (void) const +{ + if (packed ()) + { + octave_idx_type a_nr = a_fact.rows (); + octave_idx_type a_nc = a_fact.cols (); + octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); + + lu_type u (mn, a_nc, lu_elt_type (0.0)); + + for (octave_idx_type i = 0; i < mn; i++) + { + for (octave_idx_type j = i; j < a_nc; j++) + u.xelem (i, j) = a_fact.xelem (i, j); + } + + return u; + } + else + return a_fact; +} + +template <class lu_type> +lu_type +base_lu <lu_type> :: Y (void) const +{ + if (! packed ()) + (*current_liboctave_error_handler) + ("lu: Y () not implemented for unpacked form"); + return a_fact; +} + +template <class lu_type> +Array<octave_idx_type> +base_lu <lu_type> :: getp (void) const +{ + if (packed ()) + { + octave_idx_type a_nr = a_fact.rows (); + + Array<octave_idx_type> pvt (dim_vector (a_nr, 1)); + + for (octave_idx_type i = 0; i < a_nr; i++) + pvt.xelem (i) = i; + + for (octave_idx_type i = 0; i < ipvt.length (); i++) + { + octave_idx_type k = ipvt.xelem (i); + + if (k != i) + { + octave_idx_type tmp = pvt.xelem (k); + pvt.xelem (k) = pvt.xelem (i); + pvt.xelem (i) = tmp; + } + } + + return pvt; + } + else + return ipvt; +} + +template <class lu_type> +PermMatrix +base_lu <lu_type> :: P (void) const +{ + return PermMatrix (getp (), false); +} + +template <class lu_type> +ColumnVector +base_lu <lu_type> :: P_vec (void) const +{ + octave_idx_type a_nr = a_fact.rows (); + + ColumnVector p (a_nr); + + Array<octave_idx_type> pvt = getp (); + + for (octave_idx_type i = 0; i < a_nr; i++) + p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1); + + return p; +} + +template <class lu_type> +bool +base_lu<lu_type>::regular (void) const +{ + bool retval = true; + + octave_idx_type k = std::min (a_fact.rows (), a_fact.columns ()); + + for (octave_idx_type i = 0; i < k; i++) + { + if (a_fact(i, i) == lu_elt_type ()) + { + retval = false; + break; + } + } + + return retval; +} diff -uNr a/src/cyclgen.cc b/src/cyclgen.cc --- a/src/cyclgen.cc 2015-04-04 12:28:43.942510204 -0400 +++ b/src/cyclgen.cc 2018-04-09 13:51:41.852070262 -0400 @@ -29,7 +29,7 @@ const Array<int>& x, const int& n) { - int x_len = x.length (); + int x_len = x.numel (); Array<int> si (dim_vector (n, 1), 0); Array<int> y (dim_vector (x_len, 1), 0); diff -uNr a/src/cyclgen.cc~ b/src/cyclgen.cc~ --- a/src/cyclgen.cc~ 1969-12-31 19:00:00.000000000 -0500 +++ b/src/cyclgen.cc~ 2015-04-04 12:28:43.942510204 -0400 @@ -0,0 +1,278 @@ +//Copyright (C) 2003 David Bateman +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, see +// <http://www.gnu.org/licenses/>. +// +// In addition to the terms of the GPL, you are permitted to link this +// program with any Open Source program, as defined by the Open Source +// Initiative (www.opensource.org) + +#include <string> + +#include <octave/oct.h> + +// A simplified version of the filter function for specific lengths of a and b +// in the Galois field GF(2) +Array<int> +filter_gf2 (const Array<int>& b, const Array<int>& a, + const Array<int>& x, const int& n) +{ + + int x_len = x.length (); + Array<int> si (dim_vector (n, 1), 0); + Array<int> y (dim_vector (x_len, 1), 0); + + for (int i = 0; i < x_len; i++) + { + y(i) = si(0); + if (b(0) && x(i)) + y(i) ^= 1; + + for (int j = 0; j < n - 1; j++) + { + si(j) = si(j+1); + if (a(j+1) && y(i)) + si(j) ^= 1; + if (b(j+1) && x(i)) + si(j) ^= 1; + } + si(n-1) = 0; + if (a(n) && y(i)) + si(n-1) ^= 1; + if (b(n) && x(i)) + si(n-1) ^= 1; + } + + return y; +} + +// Cyclic polynomial is irreducible. I.E. it divides into x^n-1 +// without remainder There must surely be an easier way of doing this +// as the polynomials are over GF(2). +static bool +do_is_cyclic_polynomial (const Array<int>& a, const int& n, const int& m) +{ + Array<int> y (dim_vector (n+1, 1), 0); + Array<int> x (dim_vector (n-m+2, 1), 0); + y(0) = 1; + y(n) = 1; + x(0) = 1; + + Array<int> b = filter_gf2 (y, a, x, n); + b.resize (dim_vector (n+1, 1), 0); + Array<int> p (dim_vector (m+1, 1), 0); + p(0) = 1; + Array<int> q = filter_gf2 (a, p, b, m); + + for (int i = 0; i < n+1; i++) + if (y(i) ^ q(i)) + return false; + + return true; +} + +DEFUN_DLD (cyclgen, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{h} =} cyclgen (@var{n}, @var{p})\n\ +@deftypefnx {Loadable Function} {@var{h} =} cyclgen (@var{n}, @var{p}, @var{typ})\n\ +@deftypefnx {Loadable Function} {[@var{h}, @var{g}] =} cyclgen (@dots{})\n\ +@deftypefnx {Loadable Function} {[@var{h}, @var{g}, @var{k}] =} cyclgen (@dots{})\n\ +Produce the parity check and generator matrix of a cyclic code. The parity\n\ +check matrix is returned as a @var{m} by @var{n} matrix, representing the\n\ +[@var{n},@var{k}] cyclic code. @var{m} is the order of the generator\n\ +polynomial @var{p} and the message length @var{k} is given by\n\ +@code{@var{n} - @var{m}}.\n\ +\n\ +The generator polynomial can either be a vector of ones and zeros,\n\ +and length @var{m} representing,\n\ +@tex\n\ +$$ p_0 + p_1 x + p_2 x^2 + \\cdots + p_m x^{m-1} $$\n\ +@end tex\n\ +@ifnottex\n\ +\n\ +@example\n\ +@var{p}(1) + @var{p}(2) * x + @var{p}(3) * x^2 + ... + @var{p}(@var{m}) * x^(m-1)\n\ +@end example\n\ +@end ifnottex\n\ +\n\ +The terms of the polynomial are stored least-significant term first.\n\ +Alternatively, @var{p} can be an integer representation of the same\n\ +polynomial.\n\ +\n\ +The form of the parity check matrix is determined by @var{typ}. If\n\ +@var{typ} is 'system', a systematic parity check matrix is produced. If\n\ +@var{typ} is 'nosys' and non-systematic parity check matrix is produced.\n\ +\n\ +If requested @code{cyclgen} also returns the @var{k} by @var{n} generator\n\ +matrix @var{g}.\ +\n\ +@seealso{hammgen, gen2par, cyclpoly}\n\ +@end deftypefn") +{ + octave_value_list retval; + int nargin = args.length (); + unsigned long long p = 0; + int n, m, k, mm; + bool system = true; + Array<int> pp; + + if (nargin < 2 || nargin > 3) + { + print_usage (); + return retval; + } + + n = args(0).int_value (); + m = 1; + while (n > (1<<(m+1))) + m++; + pp.resize (dim_vector (n+1, 1), 0); + + if (args(1).is_scalar_type ()) + { + p = (unsigned long long)(args(1).int_value ()); + mm = 1; + while (p > ((unsigned long long)1<<(mm+1))) + mm++; + for (int i = 0; i < mm+1; i++) + pp(i) = (p & (1<<i) ? 1 : 0); + } + else + { + Matrix tmp = args(1).matrix_value (); + if ((tmp.rows () != 1) && (tmp.columns () != 1)) + { + error ("cyclgen: generator polynomial must be a vector"); + return retval; + } + + if (tmp.rows () == 1) + { + mm = tmp.columns (); + for (int j = 0; j < mm; j++) { + if (tmp(0, j) == 1) { + p |= ((unsigned long long)1 << j); + pp(j) = 1; + } + else if (tmp(0, j) != 0) { + error ("cyclgen: illegal generator polynomial"); + return retval; + } + } + } + else + { + mm = tmp.rows (); + for (int i = 0; i < mm; i++) + { + if (tmp(i, 0) == 1) + { + p |= ((unsigned long long)1 << i); + pp(i) = 1; + } + else if (tmp(i, 0) != 0) + { + error ("cyclgen: illegal generator polynomial"); + return retval; + } + } + } + mm = mm - 1; + } + k = n - mm; + + if (nargin > 2) + { + if (args(2).is_string ()) + { + std::string s_arg = args(2).string_value (); + + if (s_arg == "system") + system = true; + else if (s_arg == "nosys") + system = false; + else + { + error ("cyclgen: illegal argument"); + return retval; + } + } + else + { + error ("cyclgen: illegal argument"); + return retval; + } + } + + // Haven't implemented this since I'm not sure what matlab wants here + if (!system) + { + error ("cyclgen: non-systematic generator matrices not implemented"); + return retval; + } + + if (!do_is_cyclic_polynomial (pp, n, mm)) + { + error ("cyclgen: generator polynomial does not produce cyclic code"); + return retval; + } + + unsigned long long mask = 1; + unsigned long long *alpha_to = + (unsigned long long *)malloc (sizeof (unsigned long long) * n); + for (int i = 0; i < n; i++) + { + alpha_to[i] = mask; + mask <<= 1; + if (mask & ((unsigned long long)1<<mm)) + mask ^= p; + } + + Matrix parity (mm, n, 0); + for (int i = 0; i < n; i++) + for (int j = 0; j < mm; j++) + if (alpha_to[i] & ((unsigned long long)1<<j)) + parity(j, i) = 1; + + free (alpha_to); + retval(0) = octave_value (parity); + + if (nargout > 1) + { + Matrix generator (k, n, 0); + + for (int i = 0; i < (int)k; i++) + for (int j = 0; j < (int)mm; j++) + generator(i, j) = parity(j, i+mm); + for (int i = 0; i < (int)k; i++) + generator(i, i+mm) = 1; + + retval(1) = octave_value (generator); + retval(2) = octave_value ((double)k); + } + return retval; +} + +/* +%% Test input validation +%!error cyclgen () +%!error cyclgen (1) +%!error cyclgen (1, 2, 3, 4) +*/ + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ diff -uNr a/src/cyclpoly.cc b/src/cyclpoly.cc --- a/src/cyclpoly.cc 2015-04-04 12:28:43.942510204 -0400 +++ b/src/cyclpoly.cc 2018-04-09 13:52:48.344959602 -0400 @@ -38,7 +38,7 @@ const Array<int>& x, const int& n) { - int x_len = x.length (); + int x_len = x.numel (); Array<int> si (dim_vector (n, 1), 0); Array<int> y (dim_vector (x_len, 1), 0); @@ -217,8 +217,8 @@ for (unsigned long long i = (1UL<<m)+1; i < (1UL<<(1+m)); i+=2) if (do_is_cyclic_polynomial (i, n, m)) { - cyclic_polys.resize (cyclic_polys.length ()+1); - cyclic_polys(cyclic_polys.length ()-1) = (double)i; + cyclic_polys.resize (cyclic_polys.numel ()+1); + cyclic_polys(cyclic_polys.numel ()-1) = (double)i; } break; case CYCLIC_POLY_L: @@ -233,8 +233,8 @@ { if (do_is_cyclic_polynomial (i, n, m)) { - cyclic_polys.resize (cyclic_polys.length ()+1); - cyclic_polys(cyclic_polys.length ()-1) = (double)i; + cyclic_polys.resize (cyclic_polys.numel ()+1); + cyclic_polys(cyclic_polys.numel ()-1) = (double)i; } } } @@ -244,7 +244,7 @@ break; } - if (cyclic_polys.length () == 0) + if (cyclic_polys.numel () == 0) { octave_stdout << "cyclpoly: no generator polynomial statifies constraints" << std::endl; @@ -254,8 +254,8 @@ { if (polyrep) { - Matrix polys (cyclic_polys.length (), m+1, 0); - for (int i = 0 ; i < cyclic_polys.length (); i++) + Matrix polys (cyclic_polys.numel (), m+1, 0); + for (int i = 0 ; i < cyclic_polys.numel (); i++) for (int j = 0; j < m+1; j++) if ((unsigned long long)cyclic_polys(i) & (1<<j)) polys(i, j) = 1; diff -uNr a/src/cyclpoly.cc~ b/src/cyclpoly.cc~ --- a/src/cyclpoly.cc~ 1969-12-31 19:00:00.000000000 -0500 +++ b/src/cyclpoly.cc~ 2018-04-09 13:52:05.442966638 -0400 @@ -0,0 +1,282 @@ +//Copyright (C) 2003 David Bateman +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, see +// <http://www.gnu.org/licenses/>. +// +// In addition to the terms of the GPL, you are permitted to link this +// program with any Open Source program, as defined by the Open Source +// Initiative (www.opensource.org) + +#include <iostream> +#include <string> + +#include <octave/oct.h> + +enum cyclic_poly_type +{ + CYCLIC_POLY_MIN=0, + CYCLIC_POLY_MAX, + CYCLIC_POLY_ALL, + CYCLIC_POLY_L +}; + +// A simplified version of the filter function for specific lengths of +// a and b in the Galois field GF(2) +Array<int> +filter_gf2 (const Array<int>& b, const Array<int>& a, + const Array<int>& x, const int& n) +{ + + int x_len = x.numel (); + Array<int> si (dim_vector (n, 1), 0); + Array<int> y (dim_vector (x_len, 1), 0); + + for (int i = 0; i < x_len; i++) + { + y(i) = si(0); + if (b(0) && x(i)) + y(i) ^= 1; + + for (int j = 0; j < n - 1; j++) + { + si(j) = si(j+1); + if (a(j+1) && y(i)) + si(j) ^= 1; + if (b(j+1) && x(i)) + si(j) ^= 1; + } + si(n-1) = 0; + if (a(n) && y(i)) + si(n-1) ^= 1; + if (b(n) && x(i)) + si(n-1) ^= 1; + } + + return y; +} + +// Cyclic polynomial is irreducible. I.E. it divides into x^n-1 +// without remainder There must surely be an easier way of doing this +// as the polynomials are over GF(2). +static bool +do_is_cyclic_polynomial (const unsigned long long& a1, const int& n, + const int& m) +{ + Array<int> a (dim_vector (n+1, 1), 0); + Array<int> y (dim_vector (n+1, 1), 0); + Array<int> x (dim_vector (n-m+2, 1), 0); + y(0) = 1; + y(n) = 1; + x(0) = 1; + for (int i=0; i < m+1; i++) + a(i) = (a1 & (1UL << i) ? 1 : 0); + + Array<int> b = filter_gf2 (y, a, x, n); + b.resize(dim_vector (n+1, 1), 0); + Array<int> p (dim_vector (m+1, 1), 0); + p(0) = 1; + Array<int> q = filter_gf2 (a, p, b, m); + + for (int i=0; i < n+1; i++) + if (y(i) ^ q(i)) + return false; + + return true; +} + +DEFUN_DLD (cyclpoly, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} cyclpoly (@var{n}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{y} =} cyclpoly (@var{n}, @var{k}, @var{opt})\n\ +@deftypefnx {Loadable Function} {@var{y} =} cyclpoly (@var{n}, @var{k}, @var{opt}, @var{rep})\n\ +This function returns the cyclic generator polynomials of the code\n\ +[@var{n},@var{k}]. By default the polynomial with the smallest weight\n\ +is returned. However this behavior can be overridden with the @var{opt}\n\ +flag. Valid values of @var{opt} are:\n\ +\n\ +@table @asis\n\ +@item @code{\"all\"}\n\ +Returns all of the polynomials of the code [@var{n},@var{k}]\n\ +@item @code{\"min\"}\n\ +Returns the polynomial of minimum weight of the code [@var{n},@var{k}]\n\ +@item @code{\"max\"}\n\ +Returns the polynomial of the maximum weight of the code [@var{n},@var{k}]\n\ +@item @var{l}\n\ +Returns the polynomials having exactly the weight @var{l}\n\ +@end table\n\ +\n\ +The polynomials are returns as row-vectors in the variable @var{y}. Each\n\ +row of @var{y} represents a polynomial with the least-significant term\n\ +first. The polynomials can be returned with an integer representation\n\ +if @var{rep} is @code{\"integer\"}. The default behavior is given if @var{rep}\n\ +is @code{\"polynomial\"}.\n\ +@seealso{gf, isprimitive}\n\ +@end deftypefn") +{ + octave_value retval; + int nargin = args.length (); + bool polyrep = true; + enum cyclic_poly_type type = CYCLIC_POLY_MIN; + RowVector cyclic_polys; + int l=0; + + if (nargin < 2 || nargin > 4) + { + print_usage (); + return retval; + } + + int n = args(0).int_value (); + int k = args(1).int_value ();; + + if (n < 1) + { + error ("cyclpoly: n must be 1 or greater"); + return retval; + } + + if (n <= k) + { + error ("cyclpoly: k must be less than n"); + return retval; + } + + for (int i = 2; i < nargin; i++) + { + if (args(i).is_scalar_type ()) + { + l = args(i).int_value (); + type = CYCLIC_POLY_L; + } + else if (args(i).is_string ()) + { + std::string s_arg = args(i).string_value (); + + if (s_arg == "integer") + polyrep = false; + else if (s_arg == "polynomial") + polyrep = true; + else if (s_arg == "min") + type = CYCLIC_POLY_MIN; + else if (s_arg == "max") + type = CYCLIC_POLY_MAX; + else if (s_arg == "all") + type = CYCLIC_POLY_ALL; + else + { + error ("cyclpoly: invalid argument"); + return retval; + } + } + else + { + error ("cyclpoly: incorrect argument type"); + return retval; + } + } + + int m = n - k; + + // Matlab code seems to think that 1+x+x^3 is of larger weight than + // 1+x^2+x^3. So for matlab compatiability the list of polynomials + // should be reversed by replacing "i+=2" with "i-=2" and visa-versa. + // Thats not going to happen!!! + + switch (type) + { + case CYCLIC_POLY_MIN: + cyclic_polys.resize (1); + for (unsigned long long i = (1UL<<m)+1; i < (1UL<<(1+m)); i+=2) + if (do_is_cyclic_polynomial (i, n, m)) + { + cyclic_polys(0) = (double)i; + break; + } + break; + case CYCLIC_POLY_MAX: + cyclic_polys.resize (1); + for (unsigned long long i = (1UL<<(m+1))-1; i > (1UL<<m); i-=2) + if (do_is_cyclic_polynomial (i, n, m)) + { + cyclic_polys(0) = (double)i; + break; + } + break; + case CYCLIC_POLY_ALL: + for (unsigned long long i = (1UL<<m)+1; i < (1UL<<(1+m)); i+=2) + if (do_is_cyclic_polynomial (i, n, m)) + { + cyclic_polys.resize (cyclic_polys.length ()+1); + cyclic_polys(cyclic_polys.length ()-1) = (double)i; + } + break; + case CYCLIC_POLY_L: + for (unsigned long long i = ((unsigned long long)1<<m)+1; + i < ((unsigned long long)1<<(1+m)); i+=2) + { + int li = 0; + for (int j=0; j < m+1; j++) + if (i & ((unsigned long long)1 << j)) + li++; + if (li == l) + { + if (do_is_cyclic_polynomial (i, n, m)) + { + cyclic_polys.resize (cyclic_polys.length ()+1); + cyclic_polys(cyclic_polys.length ()-1) = (double)i; + } + } + } + break; + default: + error ("cyclpoly: impossible"); + break; + } + + if (cyclic_polys.length () == 0) + { + octave_stdout << + "cyclpoly: no generator polynomial statifies constraints" << std::endl; + retval = octave_value (Matrix (0, 0)); + } + else + { + if (polyrep) + { + Matrix polys (cyclic_polys.length (), m+1, 0); + for (int i = 0 ; i < cyclic_polys.length (); i++) + for (int j = 0; j < m+1; j++) + if ((unsigned long long)cyclic_polys(i) & (1<<j)) + polys(i, j) = 1; + retval = octave_value (polys); + } + else + retval = octave_value (cyclic_polys); + } + + return retval; +} + +/* +%% Test input validation +%!error cyclpoly () +%!error cyclpoly (1) +%!error cyclpoly (1, 2, 3, 4, 5) +*/ + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ diff -uNr a/src/galois.cc b/src/galois.cc --- a/src/galois.cc 2018-04-09 13:25:42.880981256 -0400 +++ b/src/galois.cc 2018-04-09 13:53:27.547125644 -0400 @@ -19,7 +19,7 @@ // Initiative (www.opensource.org) #include <octave/error.h> -#include <octave/gripes.h> +#include <octave/errwarn.h> #include <octave/mx-op-defs.h> #include "galois.h" @@ -215,7 +215,7 @@ if (nr != a_nr || nc != a_nc) { - gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc); + octave::err_nonconformant ("operator +=", nr, nc, a_nr, a_nc); return *this; } @@ -251,7 +251,7 @@ if (nr != a_nr || nc != a_nc) { - gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc); + octave::err_nonconformant ("operator -=", nr, nc, a_nr, a_nc); return *this; } @@ -517,7 +517,7 @@ { if (a_nr != b_nr || a_nc != b_nc) { - gripe_nonconformant ("operator .^", a_nr, a_nc, a_nr, a_nc); + octave::err_nonconformant ("operator .^", a_nr, a_nc, a_nr, a_nc); return galois (); } @@ -548,7 +548,7 @@ if (a_nr != b_nr || a_nc != b_nc) { - gripe_nonconformant ("operator .^", a_nr, a_nc, b_nr, b_nc); + octave::err_nonconformant ("operator .^", a_nr, a_nc, b_nr, b_nc); return galois (); } @@ -757,7 +757,7 @@ return product (a, b); else if (a_nc != b_nr) { - gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); + octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); return galois (); } else @@ -1302,7 +1302,7 @@ // Apply row interchanges to the right hand sides. //for (int j = 0; j < IP.length (); j++) - for (int j = IP.length ()-1; j >= 0; j--) + for (int j = IP.numel ()-1; j >= 0; j--) { int piv = IP(j); for (int i = 0; i < b_nc; i++) @@ -1334,7 +1334,7 @@ Array<int> IP (fact.ipvt); // Apply row interchanges to the right hand sides. - for (int j = 0; j < IP.length (); j++) + for (int j = 0; j < IP.numel (); j++) { int piv = IP(j); for (int i = 0; i < b_nc; i++) @@ -1419,7 +1419,7 @@ int a_nr = a.rows (); int b_nr = b.rows (); - gripe_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc); + octave::err_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc); return galois (); } @@ -1463,7 +1463,7 @@ int a_nc = a.cols (); int b_nc = b.cols (); - gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc); + octave::err_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc); return galois (); } diff -uNr a/src/galois.cc~ b/src/galois.cc~ --- a/src/galois.cc~ 1969-12-31 19:00:00.000000000 -0500 +++ b/src/galois.cc~ 2018-04-09 13:37:49.607001400 -0400 @@ -0,0 +1,1494 @@ +//Copyright (C) 2003 David Bateman +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, see +// <http://www.gnu.org/licenses/>. +// +// In addition to the terms of the GPL, you are permitted to link this +// program with any Open Source program, as defined by the Open Source +// Initiative (www.opensource.org) + +#include <octave/error.h> +#include <octave/errwarn.h> +#include <octave/mx-op-defs.h> + +#include "galois.h" +#include "galoisfield.h" +#include "galois-def.h" + +#include "base-lu.cc" + +galois_field_list stored_galois_fields; + +// galois class + +galois::galois (const Array<int>& a, const int& _m, + const int& _primpoly) : MArray<int> (a.dims ()), field (NULL) +{ + int _n = (1<<_m) - 1; + + // Check the validity of the data in the matrix + for (int i = 0; i < rows (); i++) + { + for (int j = 0; j < columns (); j++) + { + if ((a(i, j) < 0) || (a(i, j) > _n)) + { + gripe_range_galois (_m); + return; + } + xelem(i, j) = (int)a(i, j); + } + } + + field = stored_galois_fields.create_galois_field (_m, _primpoly); +} + +galois::galois (const MArray<int>& a, const int& _m, + const int& _primpoly) : MArray<int> (a.dims ()), field (NULL) +{ + int _n = (1<<_m) - 1; + + // Check the validity of the data in the matrix + for (int i = 0; i < rows (); i++) + { + for (int j = 0; j < columns (); j++) + { + if ((a(i, j) < 0) || (a(i, j) > _n)) + { + gripe_range_galois (_m); + return; + } + xelem(i, j) = (int)a(i, j); + } + } + + field = stored_galois_fields.create_galois_field (_m, _primpoly); +} + +galois::galois (const Matrix& a, const int& _m, + const int& _primpoly) : MArray<int> (a.dims ()), field (NULL) +{ + int _n = (1<<_m) - 1; + + // Check the validity of the data in the matrix + for (int i = 0; i < rows (); i++) + { + for (int j = 0; j < columns (); j++) + { + if ((a(i, j) < 0) || (a(i, j) > _n)) + { + gripe_range_galois (_m); + return; + } + if ((a(i, j) - (double)((int)a(i, j))) != 0.) + { + gripe_integer_galois (); + return; + } + xelem(i, j) = (int)a(i, j); + } + } + + field = stored_galois_fields.create_galois_field (_m, _primpoly); +} + +galois::galois (int nr, int nc, const int& val, const int& _m, + const int& _primpoly) + : MArray<int> (dim_vector (nr, nc), val), field (NULL) +{ + int _n = (1<<_m) - 1; + + // Check the validity of the data in the matrix + if ((val < 0) || (val > _n)) + { + gripe_range_galois (_m); + return; + } + + field = stored_galois_fields.create_galois_field (_m, _primpoly); +} + +galois::galois (int nr, int nc, double val, const int& _m, + const int& _primpoly) + : MArray<int> (dim_vector (nr, nc), (int)val), field (NULL) +{ + int _n = (1<<_m) - 1; + + // Check the validity of the data in the matrix + if ((val < 0) || (val > _n)) + { + gripe_range_galois (_m); + return; + } + + if ((val - (double)((int)val)) != 0.) + { + gripe_integer_galois (); + return; + } + + field = stored_galois_fields.create_galois_field (_m, _primpoly); +} + +galois::galois (const galois& a) : MArray<int> (a) +{ + + if (!a.have_field ()) + { + gripe_copy_invalid_galois (); + field = NULL; + return; + } + + // This call to create_galois_field will just increment the usage counter + field = stored_galois_fields.create_galois_field (a.m (), a.primpoly ()); +} + +galois::~galois (void) +{ + stored_galois_fields.delete_galois_field (field); + field = NULL; +} + +galois& +galois::operator = (const galois& t) +{ + if (!t.have_field ()) + { + gripe_copy_invalid_galois (); + if (have_field ()) + stored_galois_fields.delete_galois_field (field); + field = NULL; + return *this; + } + + if (have_field ()) + { + if ((m () != t.m ()) || (primpoly () != t.primpoly ())) + { + stored_galois_fields.delete_galois_field (field); + field = stored_galois_fields.create_galois_field (t.m (), t.primpoly ()); + } + } + else + field = stored_galois_fields.create_galois_field (t.m (), t.primpoly ()); + + // Copy the data + MArray<int>::operator = (t); + + return *this; +} + +galois& +galois::operator += (const galois& a) +{ + int nr = rows (); + int nc = cols (); + + int a_nr = a.rows (); + int a_nc = a.cols (); + + if (have_field () && a.have_field ()) + { + if ((m () != a.m ()) || (primpoly () != a.primpoly ())) + { + gripe_differ_galois (); + return *this; + } + } + else + { + gripe_invalid_galois (); + return *this; + } + + if (nr != a_nr || nc != a_nc) + { + octave::err_nonconformant ("operator +=", nr, nc, a_nr, a_nc); + return *this; + } + + for (int i = 0; i < rows (); i++) + for (int j = 0; j < columns (); j++) + xelem(i, j) ^= a (i, j); + + return *this; +} + +galois& +galois::operator -= (const galois& a) +{ + int nr = rows (); + int nc = cols (); + + int a_nr = a.rows (); + int a_nc = a.cols (); + + if (have_field () && a.have_field ()) + { + if ((m () != a.m ()) || (primpoly () != a.primpoly ())) + { + gripe_differ_galois (); + return *this; + } + } + else + { + gripe_invalid_galois (); + return *this; + } + + if (nr != a_nr || nc != a_nc) + { + octave::err_nonconformant ("operator -=", nr, nc, a_nr, a_nc); + return *this; + } + + for (int i = 0; i < rows (); i++) + for (int j = 0; j < columns (); j++) + xelem(i, j) ^= a (i, j); + + return *this; +} + +galois +galois::index (idx_vector& i, int resize_ok, const int& rfv) const +{ + galois retval (MArray<int>::index(i, resize_ok, rfv), m (), primpoly ()); + + return retval; +} + +galois +galois::index (idx_vector& i, idx_vector& j, int resize_ok, + const int& rfv) const +{ + galois retval (MArray<int>::index(i, j, resize_ok, rfv), m (), primpoly ()); + + return retval; +} + +galois +galois::concat (const galois& rb, const Array<int>& ra_idx) +{ + if (rb.numel () > 0) + insert (rb, ra_idx(0), ra_idx(1)); + return *this; +} + +galois +galois::concat (const Matrix& rb, const Array<int>& ra_idx) +{ + if (numel () == 1) + return *this; + + galois tmp (0, 0, 0, m (), primpoly ()); + int _n = (1<<m ()) - 1; + int r = rb.rows (); + int c = rb.columns (); + tmp.resize (dim_vector (r, c)); + + // Check the validity of the data in the matrix + for (int i = 0; i < r; i++) + { + for (int j = 0; j < c; j++) + { + if ((rb(i, j) < 0) || (rb(i, j) > _n)) + { + gripe_range_galois (m ()); + return *this; + } + if ((rb(i, j) - (double)((int)rb(i, j))) != 0.) + { + gripe_integer_galois (); + return *this; + } + tmp(i, j) = (int)rb(i, j); + } + } + + insert (tmp, ra_idx(0), ra_idx(1)); + return *this; +} + +galois +concat (const Matrix& ra, const galois& rb, const Array<int>& ra_idx) +{ + galois retval (0, 0, 0, rb.m (), rb.primpoly ()); + int _n = (1<<rb.m ()) - 1; + int r = ra.rows (); + int c = ra.columns (); + retval.resize (dim_vector (r, c)); + if (ra.numel () < 1) + return retval; + + // FIXME: + // Check the validity of the data in the matrix. This is problematic + // as "ra" is not initialized on the initial resize and so contains + // random data that will be replaced. Humm, disable for now + for (int i = 0; i < r; i++) + { + for (int j = 0; j < c; j++) + { +#if 0 + if ((ra(i, j) < 0) || (ra(i, j) > _n)) + { + gripe_range_galois (rb.m ()); + return retval; + } + if ((ra(i, j) - (double)((int)ra(i, j))) != 0.) + { + gripe_integer_galois (); + return retval; + } + retval(i, j) = (int)ra(i, j); +#else + int tmp = (int)ra(i, j); + if (tmp < 0) + retval(i, j) = 0; + else if (tmp > _n) + retval(i, j) = _n; + else + retval(i, j) = tmp; +#endif + } + } + + retval.insert (rb, ra_idx(0), ra_idx(1)); + return retval; +} + +galois& +galois::insert (const galois& t, int r, int c) +{ + if ((m () != t.m ()) || (primpoly () != t.primpoly ())) + (*current_liboctave_error_handler) ("inserted galois variable must " + "be in the same field"); + else + Array<int>::insert (t, r, c); + return *this; +} + +galois +galois::diag (void) const +{ + return diag (0); +} + +galois +galois::diag (int k) const +{ + int nnr = rows (); + int nnc = cols (); + galois retval (0, 0, 0, m (), primpoly ()); + + if (k > 0) + nnc -= k; + else if (k < 0) + nnr += k; + + if (nnr > 0 && nnc > 0) + { + int ndiag = (nnr < nnc) ? nnr : nnc; + retval.resize (dim_vector (ndiag, 1)); + + if (k > 0) + { + for (int i = 0; i < ndiag; i++) + retval(i, 0) = xelem (i, i+k); + } + else if ( k < 0) + { + for (int i = 0; i < ndiag; i++) + retval(i, 0) = xelem (i-k, i); + } + else + { + for (int i = 0; i < ndiag; i++) + retval(i, 0) = xelem (i, i); + } + } + else + error ("diag: requested diagonal out of range"); + + return retval; +} + +// unary operations + +boolMatrix +galois::operator ! (void) const +{ + int nr = rows (); + int nc = cols (); + + boolMatrix b (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + b (i, j) = ! xelem (i, j); + + return b; +} + +galois +galois::transpose (void) const +{ + galois a (Matrix (0, 0), m (), primpoly ()); + int d1 = rows (); + int d2 = cols (); + + a.resize (dim_vector (d2, d1)); + for (int j = 0; j < d2; j++) + for (int i = 0; i < d1; i++) + a (j, i) = xelem (i, j); + + return a; +} + +static inline int +modn (int x, int m, int n) +{ + while (x >= n) + { + x -= n; + x = (x >> m) + (x & n); + } + return x; +} + +galois +elem_pow (const galois& a, const galois& b) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + galois result (a_nr, a_nc, 0, a.m (), a.primpoly ()); + + int b_nr = b.rows (); + int b_nc = b.cols (); + + if (a.have_field () && b.have_field ()) + { + if ((a.m () != b.m ()) || (a.primpoly () != b.primpoly ())) + { + gripe_differ_galois (); + return galois (); + } + } + else + { + gripe_invalid_galois (); + return galois (); + } + + if (a_nr == 1 && a_nc == 1) + { + result.resize (dim_vector (b_nr, b_nc), 0); + int tmp = a.index_of (a(0, 0)); + for (int j = 0; j < b_nc; j++) + for (int i = 0; i < b_nr; i++) + if (b(i, j) == 0) + result(i, j) = 1; + else if (a(0, 0) != 0) + result(i, j) = a.alpha_to (modn (tmp * b(i, j), a.m (), a.n ())); + } + else if (b_nr == 1 && b_nc == 1) + { + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < a_nr; i++) + if (b(0, 0) == 0) + result(i, j) = 1; + else if (a(i, j) != 0) + result(i, j) = a.alpha_to (modn (a.index_of (a(i, j)) * + b(0, 0), a.m (), a.n ())); + } + else + { + if (a_nr != b_nr || a_nc != b_nc) + { + octave::err_nonconformant ("operator .^", a_nr, a_nc, a_nr, a_nc); + return galois (); + } + + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < a_nr; i++) + if (b(i, j) == 0) + result(i, j) = 1; + else if (a(i, j) != 0) + result(i, j) = a.alpha_to (modn (a.index_of (a(i, j)) * + b(i, j), a.m (), a.n ())); + } + + return result; +} + +galois +elem_pow (const galois& a, const Matrix& b) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + galois result (a_nr, a_nc, 0, a.m (), a.primpoly ()); + + int b_nr = b.rows (); + int b_nc = b.cols (); + + if (b_nr == 1 && b_nc == 1) + return elem_pow (a, b(0, 0)); + + if (a_nr != b_nr || a_nc != b_nc) + { + octave::err_nonconformant ("operator .^", a_nr, a_nc, b_nr, b_nc); + return galois (); + } + + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < a_nr; i++) + { + int tmp = (int)b(i, j); + while (tmp < 0) + tmp += a.n (); + if (tmp == 0) + result(i, j) = 1; + else if (a(i, j) != 0) + result(i, j) = a.alpha_to (modn (a.index_of (a(i, j)) * tmp, + a.m (), a.n ())); + } + return result; +} + +galois +elem_pow (const galois& a, double b) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + galois result (a_nr, a_nc, 0, a.m (), a.primpoly ()); + int bi = (int) b; + + if ((double)bi != b) + { + gripe_integer_galois (); + return galois (); + } + + while (bi < 0) + bi += a.n (); + + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < a_nr; i++) + { + if (bi == 0) + result(i, j) = 1; + else if (a(i, j) != 0) + result(i, j) = a.alpha_to (modn (a.index_of (a(i, j)) * + bi, a.m (), a.n ())); + } + return result; +} + +galois +elem_pow (const galois &a, int b) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + galois result (a_nr, a_nc, 0, a.m (), a.primpoly ()); + + while (b < 0) + b += a.n (); + + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < a_nr; i++) + { + if (b == 0) + result(i, j) = 1; + else if (a(i, j) != 0) + result(i, j) = a.alpha_to (modn (a.index_of (a(i, j)) * b, + a.m (), a.n ())); + } + return result; +} + +galois +pow (const galois& a, double b) +{ + int bi = (int)b; + if ((double)bi != b) + { + gripe_integer_power_galois (); + return galois (); + } + + return pow (a, bi); +} + +galois +pow (const galois& a, const galois& b) +{ + int nr = b.rows (); + int nc = b.cols (); + + if (a.have_field () && b.have_field ()) + { + if ((a.m () != b.m ()) || (a.primpoly () != b.primpoly ())) + { + gripe_differ_galois (); + return galois (); + } + } + else + { + gripe_invalid_galois (); + return galois (); + } + + if (nr != 1 || nc != 1) + { + gripe_square_galois (); + return galois (); + } + else + return pow (a, b(0, 0)); +} + +galois +pow (const galois& a, int b) +{ + galois retval; + int nr = a.rows (); + int nc = a.cols (); + + if (!a.have_field ()) + { + gripe_invalid_galois (); + return retval; + } + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_galois (); + else if (b == 0) + { + retval = galois (nr, nc, 0, a.m (), a.primpoly ()); + for (int i = 0; i < nr; i++) + retval(i, i) = 1; + } + else + { + galois atmp; + + if (b < 0 ) + { + atmp = a.inverse (); + b = abs (b); + } + else + atmp = a; + + retval = atmp; + b--; + while (b > 0) + { + if (b & 1) + retval = retval * atmp; + + b >>= 1; + + if (b > 0) + atmp = atmp * atmp; + } + } + + return retval; +} + +galois +operator * (const Matrix& a, const galois& b) +{ + galois tmp (a, b.m (), b.primpoly ()); + + OCTAVE_QUIT; + + return tmp * b; +} + +galois +operator * (const galois& a, const Matrix& b) +{ + galois tmp (b, a.m (), a.primpoly ()); + + OCTAVE_QUIT; + + return a * tmp; +} + +galois +operator * (const galois& a, const galois& b) +{ + if (a.have_field () && b.have_field ()) + { + if ((a.m () != b.m ()) || (a.primpoly () != b.primpoly ())) + { + gripe_differ_galois (); + return galois (); + } + } + else + { + gripe_invalid_galois (); + return galois (); + } + + int a_nr = a.rows (); + int a_nc = a.cols (); + + int b_nr = b.rows (); + int b_nc = b.cols (); + + if ((a_nr == 1 && a_nc == 1) || (b_nr == 1 && b_nc == 1)) + return product (a, b); + else if (a_nc != b_nr) + { + octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); + return galois (); + } + else + { + galois retval (a_nr, b_nc, 0, a.m (), a.primpoly ()); + if (a_nr != 0 && a_nc != 0 && b_nc != 0) + { + // This is not optimum for referencing b, but can use vector + // to represent index(a(k,j)). Seems to be the fastest. + galois c (a_nr, 1, 0, a.m (), a.primpoly ()); + for (int j = 0; j < b_nr; j++) + { + for (int k = 0; k < a_nr; k++) + c(k, 0) = a.index_of (a(k, j)); + + for (int i = 0; i < b_nc; i++) + if (b(j, i) != 0) + { + int tmp = a.index_of (b(j, i)); + for (int k = 0; k < a_nr; k++) + { + if (a(k, j) != 0) + retval(k, i) = retval(k, i) + ^ a.alpha_to (modn (tmp + c(k, 0), + a.m (), a.n ())); + } + } + } + } + return retval; + } +} + +// Other operators +boolMatrix +galois::all (int dim) const +{ + return do_mx_red_op<bool, int> (*this, dim, mx_inline_all); +} + +boolMatrix +galois::any (int dim) const +{ + return do_mx_red_op<bool, int> (*this, dim, mx_inline_any); +} + +galois +galois::prod (int dim) const +{ + if (!have_field ()) + { + gripe_invalid_galois (); + return galois (); + } + + galois retval (0, 0, 0, m (), primpoly ()); + +#define ROW_EXPR \ + if ((retval(i, 0) == 0) || (elem (i, j) == 0)) \ + retval(i, 0) = 0; \ + else \ + retval(i, 0) = alpha_to (modn (index_of (retval(i, 0)) + \ + index_of (elem (i, j)), m (), n ())); + +#define COL_EXPR \ + if ((retval(0, j) == 0) || (elem (i, j) == 0)) \ + retval(0, j) = 0; \ + else \ + retval(0, j) = alpha_to (modn (index_of (retval(0, j)) + \ + index_of (elem (i, j)), m (), n ())); + + GALOIS_REDUCTION_OP (retval, ROW_EXPR, COL_EXPR, 1, 1); + return retval; + +#undef ROW_EXPR +#undef COL_EXPR +} + +galois +galois::sum (int dim) const +{ + if (!have_field ()) + { + gripe_invalid_galois (); + return galois (); + } + + galois retval (0, 0, 0, m (), primpoly ()); + + +#define ROW_EXPR \ + retval(i, 0) ^= elem (i, j); + +#define COL_EXPR \ + retval(0, j) ^= elem (i, j); + + GALOIS_REDUCTION_OP (retval, ROW_EXPR, COL_EXPR, 0, 0); + return retval; + +#undef ROW_EXPR +#undef COL_EXPR +} + +galois +galois::sumsq (int dim) const +{ + if (!have_field ()) + { + gripe_invalid_galois (); + return galois (); + } + + galois retval (0, 0, 0, m (), primpoly ()); + +#define ROW_EXPR \ + if (elem (i, j) != 0) \ + retval(i, 0) ^= alpha_to (modn (2*index_of (elem (i, j)), m (), n ())); + +#define COL_EXPR \ + if (elem (i, j) != 0) \ + retval(0, j) ^= alpha_to (modn (2*index_of (elem (i, j)), m (), n ())); + + GALOIS_REDUCTION_OP (retval, ROW_EXPR, COL_EXPR, 0, 0); + return retval; + +#undef ROW_EXPR +#undef COL_EXPR +} + +galois +galois::sqrt (void) const +{ + galois retval (*this); + int nr = rows (); + int nc = cols (); + + for (int j = 0; j < nc; j++) + { + for (int i = 0; i < nr; i++) + if (retval.index_of (retval(i, j)) & 1) + retval(i, j) = retval.alpha_to ((retval.index_of (retval(i, j)) + + retval.n ()) / 2); + else + retval(i, j) = retval.alpha_to (retval.index_of (retval(i, j)) + / 2); + } + return retval; +} + +galois +galois::log (void) const +{ + bool warned = false; + if (!have_field ()) + { + gripe_invalid_galois (); + return galois (); + } + + galois retval (*this); + int nr = rows (); + int nc = cols (); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + if (retval(i, j) == 0) + { + if (!warned) + { + warning ("log of zero undefined in Galois field"); + warned = true; + } + // How do I flag a NaN without either + // 1) Having to check everytime that the data is valid + // 2) Causing overflow in alpha_to or index_of!! + retval(i, j) = retval.index_of (retval(i, j)); + } + else + retval(i, j) = retval.index_of (retval(i, j)); + } + return retval; +} + +galois +galois::exp (void) const +{ + bool warned = false; + if (!have_field ()) + { + gripe_invalid_galois (); + return galois (); + } + + galois retval (*this); + int nr = rows (); + int nc = cols (); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + if (retval(i, j) == n ()) + { + if (!warned) + { + warning ("warning: exp of 2^m-1 undefined in Galois field"); + warned = true; + } + // How do I flag a NaN without either + // 1) Having to check everytime that the data is valid + // 2) Causing overflow in alpha_to or index_of!! + retval(i, j) = retval.alpha_to (retval(i, j)); + } + else + retval(i, j) = retval.alpha_to (retval(i, j)); + } + return retval; +} + +template class base_lu <galois>; + +void +galoisLU::factor (const galois& a, const pivot_type& typ) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + int mn = (a_nr > a_nc ? a_nc : a_nr); + + ptype = typ; + info = 0; + ipvt.resize (dim_vector (mn, 1)); + + a_fact = a; + + for (int j = 0; j < mn; j++) + { + int jp = j; + + // Find the pivot and test for singularity + if (ptype == galoisLU::ROW) + { + for (int i = j+1; i < a_nr; i++) + if (a_fact(i, j) > a_fact(jp, j)) + jp = i; + } + else + { + for (int i = j+1; i < a_nc; i++) + if (a_fact(j, i) > a_fact(j, jp)) + jp = i; + } + + ipvt(j) = jp; + + if (a_fact(jp, j) != 0) + { + if (ptype == galoisLU::ROW) + { + // Apply the interchange to columns 1:NC. + if (jp != j) + for (int i = 0; i < a_nc; i++) + { + int tmp = a_fact(j, i); + a_fact(j, i) = a_fact(jp, i); + a_fact(jp, i) = tmp; + } + } + else + { + // Apply the interchange to rows 1:NR. + if (jp != j) + for (int i = 0; i < a_nr; i++) + { + int tmp = a_fact(i, j); + a_fact(i, j) = a_fact(i, jp); + a_fact(i, jp) = tmp; + } + } + + // Compute elements J+1:M of J-th column. + if ( j < a_nr-1) + { + int idxj = a_fact.index_of (a_fact(j, j)); + for (int i = j+1; i < a_nr; i++) + { + if (a_fact(i, j) == 0) + a_fact(i, j) = 0; + else + a_fact(i, j) = a_fact.alpha_to (modn (a_fact.index_of (a_fact(i, j)) + - idxj + a_fact.n (), a_fact.m (), + a_fact.n ())); + } + } + } + else + { + info = 1; + } + + if (j < mn-1) + { + // Update trailing submatrix. + for (int i = j+1; i < a_nr; i++) + { + if (a_fact(i, j) != 0) + { + int idxi = a_fact.index_of (a_fact(i, j)); + for (int k = j+1; k < a_nc; k++) + { + if (a_fact(j, k) != 0) + a_fact(i, k) ^= a_fact.alpha_to (modn (a_fact.index_of (a_fact(j, k)) + + idxi, a_fact.m (), + a_fact.n ())); + } + } + } + } + } +} + +galois +galoisLU::L (void) const +{ + int a_nr = a_fact.rows (); + int a_nc = a_fact.cols (); + int mn = (a_nr < a_nc ? a_nr : a_nc); + + galois l (a_nr, mn, 0, a_fact.m (), a_fact.primpoly ()); + + for (int i = 0; i < mn; i++) + l(i, i) = 1; + + for (int j = 0; j < mn; j++) + for (int i = j+1; i < a_nr; i++) + l(i, j) = a_fact (i, j); + + return l; +} + +galois +galoisLU::U (void) const +{ + int a_nr = a_fact.rows (); + int a_nc = a_fact.cols (); + int mn = (a_nr < a_nc ? a_nr : a_nc); + + galois u (mn, a_nc, 0, a_fact.m (), a_fact.primpoly ()); + + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < (j+1 > mn ? mn : j+1); i++) + u (i, j) = a_fact (i, j); + return u; +} + +galois +galois::inverse (void) const +{ + int info; + return inverse (info); +} + +galois +galois::inverse (int& info, int force) const +{ + int nr = rows (); + int nc = cols (); + info = 0; + + if (nr != nc || nr == 0 || nc == 0) + { + (*current_liboctave_error_handler) ("inverse requires square matrix"); + return galois (); + } + else + { + int info = 0; + + // Solve with identity matrix to find the inverse. + galois btmp (nr, nr, 0, m (), primpoly ()); + for (int i = 0; i < nr; i++) + btmp(i, i) = 1; + + galois retval = solve (btmp, info, 0); + + if (info == 0) + return retval; + else + return galois (); + } +} + +galois +galois::determinant (void) const +{ + int info; + return determinant (info); +} + +galois +galois::determinant (int& info) const +{ + galois retval (1, 1, 0, m (), primpoly ()); + + int nr = rows (); + int nc = cols (); + info = -1; + + if (nr == 0 || nc == 0) + { + info = 0; + retval(0, 0) = 1; + } + else + { + galoisLU fact (*this); + + if ( ! fact.singular ()) + { + galois A (fact.a_fact); + info = 0; + + retval(0, 0) = A(0, 0); + for (int i = 1; i < nr; i++) + { + if ((retval(0, 0) == 0) || (A(i, i) == 0)) + { + retval(0, 0) = 0; + error ("What the hell are we doing here!!!"); + } + else + retval(0, 0) = alpha_to (modn (index_of (retval(0, 0)) + + index_of (A(i, i)), m (), n ())); + } + } + } + + return retval; +} + +galois +galois::solve (const galois& b) const +{ + int info; + return solve (b, info); +} + +galois +galois::solve (const galois& b, int& info) const +{ + return solve (b, info, 0); +} + +galois +galois::solve (const galois& b, int& info, + solve_singularity_handler sing_handler) const +{ + galois retval (b); + + if (!have_field () || !b.have_field ()) + { + gripe_invalid_galois (); + return galois (); + } + else if ((m () != b.m ()) || (primpoly () != b.primpoly ())) + { + gripe_differ_galois (); + return galois (); + } + + int nr = rows (); + int nc = cols (); + int b_nr = b.rows (); + int b_nc = b.cols (); + galois c (nr, 1, 0, m (), primpoly ()); + + // if (nr == 0 || nc == 0 || nr != nc || nr != b_nr) + if (nr == 0 || nc == 0 || nr != b_nr) + { + (*current_liboctave_error_handler) + ("matrix dimension mismatch solution of linear equations"); + return galois (); + } + else if (nc > nr) + { + // Under-determined system, use column interchanges. + galoisLU fact ((*this), galoisLU::COL); + + if (fact.singular ()) + { + info = -1; + if (sing_handler) + sing_handler (0.0); + else + (*current_liboctave_error_handler)("galois matrix singular"); + + return galois (); + } + else + { + galois A (fact.a_fact); + Array<int> IP (fact.ipvt); + + // Resize the number of solution rows if needed + if (nc > nr) + retval.resize (dim_vector (b_nr+nc-nr, b_nc), 0); + + //Solve L*X = B, overwriting B with X. + int mn = (nc < nr ? nc : nr); + for (int k = 0; k < mn; k++) + { + for (int i = k+1; i < nr; i++) + c(i, 0) = index_of (A(i, k)); + + for (int j = 0; j < b_nc; j++) + if (retval(k, j) != 0) + { + int idx = index_of (retval(k, j)); + for (int i = k+1; i < nr; i++) + if (A(i, k) != 0) + retval(i, j) ^= alpha_to (modn (c(i, 0) + idx, m (), n ())); + } + } + + // Solve U*X = B, overwriting B with X. + for (int k = (nc < nr ? nc-1 : nr-1); k >= 0; k--) + { + int mn = k+1 < nr ? k+1 : nr; + for (int i = 0; i < mn; i++) + c(i, 0) = index_of (A(i, k)); + mn = k < nr ? k : nr; + for (int j = 0; j < b_nc; j++) + if (retval(k, j) != 0) + { + retval(k, j) = alpha_to (modn (index_of (retval(k, j)) - + c(k, 0) + n (), m (), n ())); + int idx = index_of (retval(k, j)); + for (int i = 0; i < mn; i++) + if (A(i, k) != 0) + retval(i, j) ^= alpha_to (modn (c(i, 0) + idx, m (), n ())); + } + } + + // Apply row interchanges to the right hand sides. + //for (int j = 0; j < IP.length (); j++) + for (int j = IP.length ()-1; j >= 0; j--) + { + int piv = IP(j); + for (int i = 0; i < b_nc; i++) + { + int tmp = retval(j, i); + retval(j, i) = retval(piv, i); + retval(piv, i) = tmp; + } + } + } + } + else + { + galoisLU fact (*this); + + if (fact.singular ()) + { + info = -1; + if (sing_handler) + sing_handler (0.0); + else + (*current_liboctave_error_handler)("galois matrix singular"); + + return galois (); + } + else + { + galois A (fact.a_fact); + Array<int> IP (fact.ipvt); + + // Apply row interchanges to the right hand sides. + for (int j = 0; j < IP.length (); j++) + { + int piv = IP(j); + for (int i = 0; i < b_nc; i++) + { + int tmp = retval(j, i); + retval(j, i) = retval(piv, i); + retval(piv, i) = tmp; + } + } + + //Solve L*X = B, overwriting B with X. + int mn = (nc < nr ? nc : nr); + for (int k = 0; k < mn; k++) + { + for (int i = k+1; i < nr; i++) + c(i, 0) = index_of (A(i, k)); + for (int j = 0; j < b_nc; j++) + if (retval(k, j) != 0) + { + int idx = index_of (retval(k, j)); + for (int i = k+1; i < nr; i++) + if (A(i, k) != 0) + retval(i, j) ^= alpha_to (modn (c(i, 0) + idx, m (), n ())); + } + } + + // Solve U*X = B, overwriting B with X. + for (int k = (nc < nr ? nc-1 : nr-1); k >= 0; k--) + { + int mn = k+1 < nr ? k+1 : nr; + for (int i = 0; i < mn; i++) + c(i, 0) = index_of (A(i, k)); + mn = k < nr ? k : nr; + for (int j = 0; j < b_nc; j++) + if (retval(k, j) != 0) + { + retval(k, j) = alpha_to (modn (index_of (retval(k, j)) - + c(k, 0) + n (), m (), n ())); + int idx = index_of (retval(k, j)); + for (int i = 0; i < mn; i++) + if (A(i, k) != 0) + retval(i, j) ^= alpha_to (modn (c(i, 0) + idx, m (), n ())); + } + } + + // Resize the number of solution rows if needed + if (nc < nr) + retval.resize (dim_vector (b_nr+nc-nr, b_nc)); + + } + } + + return retval; +} + +galois +xdiv (const galois& a, const Matrix& b) +{ + galois btmp (b, a.m (), a.primpoly ()); + + return xdiv (a, btmp); +} + +galois +xdiv (const Matrix& a, const galois& b) +{ + galois atmp (a, b.m (), b.primpoly ()); + + return xdiv (atmp, b); +} + +galois +xdiv (const galois& a, const galois& b) +{ + int info = 0; + int a_nc = a.cols (); + int b_nc = b.cols (); + + // if ((a_nc != b_nc) || (b.rows () != b.cols ())) + if (a_nc != b_nc) + { + int a_nr = a.rows (); + int b_nr = b.rows (); + + octave::err_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc); + return galois (); + } + + galois atmp = a.transpose (); + galois btmp = b.transpose (); + galois result = btmp.solve (atmp, info, 0); + + if (info == 0) + return galois (result.transpose ()); + else + return galois (); +} + + +galois +xleftdiv (const galois& a, const Matrix& b) +{ + galois btmp (b, a.m (), a.primpoly ()); + + return xleftdiv (a, btmp); +} + +galois +xleftdiv (const Matrix& a, const galois& b) +{ + galois atmp (a, b.m (), b.primpoly ()); + + return xleftdiv (atmp, b); +} + +galois +xleftdiv (const galois& a, const galois& b) +{ + int info = 0; + int a_nr = a.rows (); + int b_nr = b.rows (); + + // if ((a_nr != b_nr) || (a.rows () != a.columns ())) + if (a_nr != b_nr) + { + int a_nc = a.cols (); + int b_nc = b.cols (); + + octave::err_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc); + return galois (); + } + + galois result = a.solve (b, info, 0); + + if (info == 0) + return result; + else + return galois (); +} + +MM_BIN_OPS1 (galois, galois, galois, 1, 2, GALOIS) +MM_BIN_OPS1 (galois, galois, Matrix, 1, 2, MATRIX) +MM_BIN_OPS1 (galois, Matrix, galois, 2, 1, MATRIX) + +MM_CMP_OPS1 (galois, , galois, , 1, 2, GALOIS) +MM_CMP_OPS1 (galois, , Matrix, , 1, 2, MATRIX) +MM_CMP_OPS1 (Matrix, , galois, , 2, 1, MATRIX) + +MM_BOOL_OPS1 (galois, galois, 0.0, 1, 2, GALOIS) +MM_BOOL_OPS1 (galois, Matrix, 0.0, 1, 2, MATRIX) +MM_BOOL_OPS1 (Matrix, galois, 0.0, 2, 1, MATRIX) + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ diff -uNr a/src/galois-def.h b/src/galois-def.h --- a/src/galois-def.h 2015-04-04 12:28:43.942510204 -0400 +++ b/src/galois-def.h 2018-04-09 13:37:40.547425139 -0400 @@ -137,7 +137,7 @@ r(i, j) = (int)m1(i, j) ^ (int)m2(0, 0); \ } \ else \ - gripe_nonconformant (#OP, m1_nr, m1_nc, m2_nr, m2_nc); \ + octave::err_nonconformant (#OP, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ else \ { \ @@ -221,7 +221,7 @@ } \ } \ else \ - gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ + octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ else \ if (m1_nr > 0 && m1_nc > 0) \ @@ -289,7 +289,7 @@ r(i, j) = C1 (m1(i, j)) OP C2 (m2(0, 0)); \ } \ else \ - gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ + octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ \ return r; \ @@ -350,7 +350,7 @@ OP (m2(0, 0) != ZERO); \ } \ else if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ - gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ + octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ \ return r; \ diff -uNr a/src/genqamdemod.cc b/src/genqamdemod.cc --- a/src/genqamdemod.cc 2015-04-04 12:28:43.950510022 -0400 +++ b/src/genqamdemod.cc 2018-04-09 13:33:02.852412423 -0400 @@ -36,7 +36,7 @@ int nr1 (args(0).rows ()); int nc1 (args(0).columns ()); - int arg_is_empty1 = empty_arg ("genqamdemod", nr1, nc1); + int arg_is_empty1 = args(0).isempty (); Matrix y (nr1,nc1); int nr2 (args(1).rows ()); @@ -48,7 +48,7 @@ if (arg_is_empty1 > 0) return octave_value (Matrix ()); - if (args(0).is_real_type () && args(1).is_real_type ()) + if (args(0).isreal () && args(1).isreal ()) { // Real-valued signal & constellation Matrix x (args(0).matrix_value ()); ColumnVector constellation (args(1).vector_value ()); @@ -70,7 +70,7 @@ } } } - else if (args(0).is_complex_type () || args(1).is_complex_type ()) + else if (args(0).iscomplex () || args(1).iscomplex ()) { // Complex-valued input & constellation ComplexMatrix x (args(0).complex_matrix_value ()); ComplexColumnVector constellation (args(1).complex_vector_value ()); diff -uNr a/src/gf.cc b/src/gf.cc --- a/src/gf.cc 2018-04-09 13:25:42.880981256 -0400 +++ b/src/gf.cc 2018-04-09 14:16:54.029455163 -0400 @@ -30,7 +30,8 @@ */ #include <octave/defun-dld.h> -#include <octave/gripes.h> +#include <octave/errwarn.h> +#include <octave/interpreter.h> #include <octave/oct-locbuf.h> #include <octave/ov.h> #include <octave/utils.h> @@ -72,7 +73,7 @@ // functions, this can't be done at the point. So if more default primitive // polynomials are added to galoisfield.cc, need to update the "16" here // as well!! -DEFUN_DLD (gf, args, nargout, +DEFMETHOD_DLD (gf, interp, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {@var{y} =} gf (@var{x})\n\ @deftypefnx {Loadable Function} {@var{y} =} gf (@var{x}, @var{m})\n\ @@ -121,7 +122,7 @@ install_s_gm_ops (); install_gm_s_ops (); galois_type_loaded = true; - mlock (); + interp.mlock (); } retval = new octave_galois (data, m, primpoly); @@ -141,7 +142,7 @@ if ((!galois_type_loaded) || (a.type_id () != octave_galois::static_type_id ())) - gripe_wrong_type_arg ("gdiag", a); + err_wrong_type_arg ("gdiag", a); else { galois m = ((const octave_galois&) a.get_rep ()).galois_value (); @@ -190,12 +191,12 @@ else { galois r = m.diag (k); - if (r.capacity () > 0) + if (r.numel () > 0) retval = new octave_galois (r); } } else - gripe_wrong_type_arg ("gdiag", a); + err_wrong_type_arg ("gdiag", a); } return retval; } @@ -301,7 +302,7 @@ if ((!galois_type_loaded) || (args(0).type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("greshape", args(0)); + err_wrong_type_arg ("greshape", args(0)); return retval; } galois a = ((const octave_galois&) args(0).get_rep ()).galois_value (); @@ -371,7 +372,7 @@ } \ else \ { \ - gripe_wrong_type_arg (#FCN, arg); \ + err_wrong_type_arg (#FCN, arg); \ return retval; \ } \ } \ @@ -468,7 +469,7 @@ if (!galois_type_loaded || (args(0).type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("gsqrt", args(0)); + err_wrong_type_arg ("gsqrt", args(0)); return retval; } @@ -507,7 +508,7 @@ if (!galois_type_loaded || (args(0).type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("glog", args(0)); + err_wrong_type_arg ("glog", args(0)); return retval; } @@ -546,7 +547,7 @@ if (!galois_type_loaded || (args(0).type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("gexp", args(0)); + err_wrong_type_arg ("gexp", args(0)); return retval; } @@ -577,9 +578,9 @@ galois filter (galois& b, galois& a, galois& x, galois& si) { - int ab_len = (a.length () > b.length () ? a.length () : b.length ()); + int ab_len = (a.numel () > b.numel () ? a.numel () : b.numel ()); b.resize (dim_vector (ab_len, 1), 0); - galois retval (x.length (), 1, 0, b.m (), b.primpoly ()); + galois retval (x.numel (), 1, 0, b.m (), b.primpoly ()); int norm = a(0, 0); if (norm == 0) @@ -587,43 +588,43 @@ error ("gfilter: the first element of a must be non-zero"); return galois (); } - if (si.length () != ab_len - 1) + if (si.numel () != ab_len - 1) { - error ("gfilter: si must be a vector of length max(length(a), length(b)) - 1"); + error ("gfilter: si must be a vector of length max(numel(a), numel(b)) - 1"); return galois (); } if (norm != 1) { int idx_norm = b.index_of (norm); - for (int i = 0; i < b.length (); i++) + for (int i = 0; i < b.numel (); i++) { if (b(i, 0) != 0) b(i, 0) = b.alpha_to (modn (b.index_of (b(i, 0))-idx_norm+b.n (), b.m (), b.n ())); } } - if (a.length () > 1) + if (a.numel () > 1) { a.resize (dim_vector (ab_len, 1), 0); if (norm != 1) { int idx_norm = a.index_of (norm); - for (int i = 0; i < a.length (); i++) + for (int i = 0; i < a.numel (); i++) if (a(i, 0) != 0) a(i, 0) = a.alpha_to (modn (a.index_of (a(i, 0))-idx_norm+a.n (), a.m (), a.n ())); } - for (int i = 0; i < x.length (); i++) + for (int i = 0; i < x.numel (); i++) { retval(i, 0) = si(0, 0); if ((b(0, 0) != 0) && (x(i, 0) != 0)) retval(i, 0) ^= b.alpha_to (modn (b.index_of (b(0, 0)) + b.index_of (x(i, 0)), b.m (), b.n ())); - if (si.length () > 1) + if (si.numel () > 1) { - for (int j = 0; j < si.length () - 1; j++) + for (int j = 0; j < si.numel () - 1; j++) { si(j, 0) = si(j+1, 0); if ((a(j+1, 0) != 0) && (retval(i, 0) != 0)) @@ -633,13 +634,13 @@ si(j, 0) ^= b.alpha_to (modn (b.index_of (b(j+1, 0)) + b.index_of (x(i, 0)), b.m (), b.n ())); } - si(si.length ()-1, 0) = 0; - if ((a(si.length (), 0) != 0) && (retval(i, 0) != 0)) - si(si.length ()-1, 0) ^= a.alpha_to (modn (a.index_of (a(si.length (), 0)) + si(si.numel ()-1, 0) = 0; + if ((a(si.numel (), 0) != 0) && (retval(i, 0) != 0)) + si(si.numel ()-1, 0) ^= a.alpha_to (modn (a.index_of (a(si.numel (), 0)) + a.index_of (retval(i, 0)), a.m (), a.n ())); - if ((b(si.length (), 0) != 0) && (x(i, 0) != 0)) - si(si.length ()-1, 0) ^= b.alpha_to (modn (b.index_of (b(si.length (), 0)) + if ((b(si.numel (), 0) != 0) && (x(i, 0) != 0)) + si(si.numel ()-1, 0) ^= b.alpha_to (modn (b.index_of (b(si.numel (), 0)) + b.index_of (x(i, 0)), b.m (), b.n ())); } @@ -655,26 +656,26 @@ } } } - else if (si.length () > 0) + else if (si.numel () > 0) { - for (int i = 0; i < x.length (); i++) + for (int i = 0; i < x.numel (); i++) { retval(i, 0) = si(0, 0); if ((b(0, 0) != 0) && (x(i, 0) != 0)) retval(i, 0) ^= b.alpha_to (modn (b.index_of (b(0, 0)) + b.index_of (x(i, 0)), b.m (), b.n ())); - if (si.length () > 1) + if (si.numel () > 1) { - for (int j = 0; j < si.length () - 1; j++) + for (int j = 0; j < si.numel () - 1; j++) { si(j, 0) = si(j+1, 0); if ((b(j+1, 0) != 0) && (x(i, 0) != 0)) si(j, 0) ^= b.alpha_to (modn (b.index_of (b(j+1, 0)) + b.index_of (x(i, 0)), b.m (), b.n ())); } - si(si.length ()-1, 0) = 0; - if ((b(si.length (), 0) != 0) && (x(i, 0) != 0)) - si(si.length ()-1, 0) ^= b.alpha_to (modn (b.index_of (b(si.length (), 0)) + si(si.numel ()-1, 0) = 0; + if ((b(si.numel (), 0) != 0) && (x(i, 0) != 0)) + si(si.numel ()-1, 0) ^= b.alpha_to (modn (b.index_of (b(si.numel (), 0)) + b.index_of (x(i, 0)), b.m (), b.n ())); } @@ -688,7 +689,7 @@ } } else - for (int i = 0; i < x.length (); i++) + for (int i = 0; i < x.numel (); i++) if ((b(0, 0) != 0) && (x(i, 0) != 0)) retval(i, 0) = b.alpha_to (modn (b.index_of (b(0, 0)) + b.index_of (x(i, 0)), b.m (), b.n ())); @@ -717,7 +718,7 @@ @smallexample\n\ @group\n\ N M\n\ - SUM a(k+1) y(n-k) = SUM b(k+1) x(n-k) for 1<=n<=length(x)\n\ + SUM a(k+1) y(n-k) = SUM b(k+1) x(n-k) for 1<=n<=numel(x)\n\ k=0 k=0\n\ @end group\n\ @end smallexample\n\ @@ -729,7 +730,7 @@ $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\ @end tex\n\ @ifnottex\n\ - N=length(a)-1 and M=length(b)-1.\n\ + N=numel(a)-1 and M=numel(b)-1.\n\ @end ifnottex\n\ An equivalent form of this equation is:\n\ @tex\n\ @@ -743,7 +744,7 @@ @smallexample\n\ @group\n\ N M\n\ - y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k) for 1<=n<=length(x)\n\ + y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k) for 1<=n<=numel(x)\n\ k=1 k=0\n\ @end group\n\ @end smallexample\n\ @@ -838,8 +839,8 @@ } else { - int a_len = a.length (); - int b_len = b.length (); + int a_len = a.numel (); + int b_len = b.numel (); int si_len = (a_len > b_len ? a_len : b_len) - 1; @@ -961,7 +962,7 @@ if (!galois_type_loaded || (arg.type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("glu", arg); + err_wrong_type_arg ("glu", arg); return retval; } @@ -970,7 +971,7 @@ int nr = arg.rows (); int nc = arg.columns (); - int arg_is_empty = empty_arg ("glu", nr, nc); + int arg_is_empty = arg.isempty (); if (arg_is_empty < 0) return retval; @@ -1048,13 +1049,13 @@ if (!galois_type_loaded || (arg.type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("ginverse", arg); + err_wrong_type_arg ("ginverse", arg); return retval; } galois m = ((const octave_galois&) arg.get_rep ()).galois_value (); - int arg_is_empty = empty_arg ("ginverse", nr, nc); + int arg_is_empty = arg.isempty (); if (arg_is_empty < 0) return retval; @@ -1065,7 +1066,7 @@ } if (nr != nc) { - gripe_square_matrix_required ("ginverse"); + err_square_matrix_required ("ginverse", "X"); return retval; } @@ -1142,7 +1143,7 @@ if (!galois_type_loaded || (arg.type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("gdet", arg); + err_wrong_type_arg ("gdet", arg); return retval; } @@ -1151,7 +1152,7 @@ galois m = ((const octave_galois&) arg.get_rep ()).galois_value (); - int arg_is_empty = empty_arg ("gdet", nr, nc); + int arg_is_empty = arg.isempty (); if (arg_is_empty < 0) return retval; @@ -1163,7 +1164,7 @@ if (nr != nc) { - gripe_square_matrix_required ("det"); + err_square_matrix_required ("det", "A"); return retval; } @@ -1202,7 +1203,7 @@ if (!galois_type_loaded || (arg.type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("grank", arg); + err_wrong_type_arg ("grank", arg); return retval; } @@ -1211,7 +1212,7 @@ galois m = ((const octave_galois&) arg.get_rep ()).galois_value (); - int arg_is_empty = empty_arg ("grank", nr, nc); + int arg_is_empty = arg.isempty (); if (arg_is_empty > 0) retval = 0.0; @@ -1332,7 +1333,7 @@ if (!galois_type_loaded || (args(0).type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("rsenc", args(0)); + err_wrong_type_arg ("rsenc", args(0)); return retval; } @@ -1873,7 +1874,7 @@ if (!galois_type_loaded || (args(0).type_id () != octave_galois::static_type_id ())) { - gripe_wrong_type_arg ("rsdec", args(0)); + err_wrong_type_arg ("rsdec", args(0)); return retval; } diff -uNr a/src/gf.cc~ b/src/gf.cc~ --- a/src/gf.cc~ 1969-12-31 19:00:00.000000000 -0500 +++ b/src/gf.cc~ 2018-04-09 14:12:50.168734516 -0400 @@ -0,0 +1,2809 @@ +// Copyright (C) 1994-1997 Robert Morelos-Zaragoza <owner@eccpage.com> +// Copyright (C) 2002 Phil Karn <karn@ka9q.net> +// Copyright (C) 2003 David Bateman +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, see +// <http://www.gnu.org/licenses/>. +// +// In addition to the terms of the GPL, you are permitted to link this +// program with any Open Source program, as defined by the Open Source +// Initiative (www.opensource.org) + +/* +Part of the function rsenc and the function decode_rs are from Phil Karn. See +the website http://www.ka9q.net/code/fec for more details. + +Parts of the function bchenco and bchdeco are from Robert Morelos-Zaragoza. See +the website http://www.eccpage.com for more details. Permission has been granted +for a GPL release of his code +*/ + +#include <octave/defun-dld.h> +#include <octave/errwarn.h> +#include <octave/interpreter.h> +#include <octave/oct-locbuf.h> +#include <octave/ov.h> +#include <octave/utils.h> +#include <octave/variables.h> + +#include "galois.h" +#include "ov-galois.h" + +static bool galois_type_loaded = false; + +// PKG_ADD: autoload ("isgalois", "gf.oct"); +// PKG_DEL: autoload ("isgalois", "gf.oct", "remove"); +DEFUN_DLD (isgalois, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} isgalois (@var{expr})\n\ +Return 1 if the value of the expression @var{expr} is a Galois Field.\n\ +@end deftypefn") +{ + if (args.length () != 1) + print_usage (); + else if (!galois_type_loaded) + // Can be of Galois type if the type isn't load :-/ + return octave_value (0.); + else + return octave_value (args(0).type_id () == + octave_galois::static_type_id ()); + return octave_value (); +} + +/* +%% Test input validation +%!error isgalois () +%!error isgalois (1, 2) +*/ + +// FIXME: +// I want to replace the "16" below with __OCTAVE_GALOIS_MAX_M_AS_STRING, +// but as I don't run the preprocessor when getting the help from the +// functions, this can't be done at the point. So if more default primitive +// polynomials are added to galoisfield.cc, need to update the "16" here +// as well!! +DEFMETHOD_DLD (interp, gf, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} gf (@var{x})\n\ +@deftypefnx {Loadable Function} {@var{y} =} gf (@var{x}, @var{m})\n\ +@deftypefnx {Loadable Function} {@var{y} =} gf (@var{x}, @var{m}, @var{primpoly})\n\ +Creates a Galois field array GF(2^@var{m}) from the matrix @var{x}. The\n\ +Galois field has 2^@var{m} elements, where @var{m} must be between 1 and 16.\n\ +The elements of @var{x} must be between 0 and 2^@var{m} - 1. If @var{m} is\n\ +undefined it defaults to the value 1.\n\ +\n\ +The primitive polynomial to use in the creation of Galois field can be\n\ +specified with the @var{primpoly} variable. If this is undefined a default\n\ +primitive polynomial is used. It should be noted that the primitive\n\ +polynomial must be of the degree @var{m} and it must be irreducible.\n\ +\n\ +The output of this function is recognized as a Galois field by Octave and\n\ +other matrices will be converted to the same Galois field when used in an\n\ +arithmetic operation with a Galois field.\n\ +\n\ +@seealso{isprimitive, primpoly}\n\ +@end deftypefn") +{ + Matrix data; + octave_value retval; + int nargin = args.length (); + int m = 1; + int primpoly = 0; + + if (nargin < 1 || nargin > 3) + { + print_usage (); + return retval; + } + + data = args(0).matrix_value (); + if (nargin > 1) + m = args(1).int_value (); + if (nargin > 2) + primpoly = args(2).int_value (); + + if (!galois_type_loaded) + { + octave_galois::register_type (); + install_gm_gm_ops (); + install_m_gm_ops (); + install_gm_m_ops (); + install_s_gm_ops (); + install_gm_s_ops (); + galois_type_loaded = true; + interp.mlock (); + } + + retval = new octave_galois (data, m, primpoly); + return retval; +} + +/* +%% Test input validation +%!error gf () +%!error gf (1, 2, 3, 4) +*/ + +static octave_value +make_gdiag (const octave_value& a, const octave_value& b) +{ + octave_value retval; + + if ((!galois_type_loaded) || (a.type_id () != + octave_galois::static_type_id ())) + err_wrong_type_arg ("gdiag", a); + else + { + galois m = ((const octave_galois&) a.get_rep ()).galois_value (); + int k = b.nint_value (); + + if (! error_state) + { + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0) + retval = new octave_galois (m); + else if (nr == 1 || nc == 1) + { + int roff = 0; + int coff = 0; + if (k > 0) + { + roff = 0; + coff = k; + } + else if (k < 0) + { + k = -k; + roff = k; + coff = 0; + } + + if (nr == 1) + { + int n = nc + k; + galois r (n, n, 0, m.m (), m.primpoly ()); + for (int i = 0; i < nc; i++) + r (i+roff, i+coff) = m (0, i); + retval = new octave_galois (r); + } + else + { + int n = nr + k; + galois r (n, n, 0, m.m (), m.primpoly ()); + for (int i = 0; i < nr; i++) + r (i+roff, i+coff) = m (i, 0); + retval = new octave_galois (r); + } + } + else + { + galois r = m.diag (k); + if (r.numel () > 0) + retval = new octave_galois (r); + } + } + else + err_wrong_type_arg ("gdiag", a); + } + return retval; +} + +// PKG_ADD: autoload ("gdiag", "gf.oct"); +// PKG_DEL: autoload ("gdiag", "gf.oct", "remove"); +DEFUN_DLD (gdiag, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} gdiag (@var{v}, @var{k})\n\ +Return a diagonal matrix with Galois vector @var{v} on diagonal @var{k}.\n\ +The second argument is optional. If it is positive, the vector is placed on\n\ +the @var{k}-th super-diagonal. If it is negative, it is placed on the\n\ +@var{-k}-th sub-diagonal. The default value of @var{k} is 0, and the\n\ +vector is placed on the main diagonal. For example,\n\ +\n\ +@example\n\ +gdiag (gf ([1, 2, 3], 2), 1)\n\ +ans =\n\ +GF(2^2) array. Primitive Polynomial = D^2+D+1 (decimal 7)\n\ +\n\ +Array elements =\n\ +\n\ + 0 1 0 0\n\ + 0 0 2 0\n\ + 0 0 0 3\n\ + 0 0 0 0\n\ +\n\ +@end example\n\ +@seealso{diag}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 1 && args(0).is_defined ()) + retval = make_gdiag (args(0), octave_value (0.)); + else if (nargin == 2 && args(0).is_defined () && args(1).is_defined ()) + retval = make_gdiag (args(0), args(1)); + else + print_usage (); + + return retval; +} + +/* +%% Test input validation +%!error gdiag () +%!error gdiag (1, 2, 3) +*/ + +// PKG_ADD: autoload ("greshape", "gf.oct"); +// PKG_DEL: autoload ("greshape", "gf.oct", "remove"); +DEFUN_DLD (greshape, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} greshape (@var{a}, @var{m}, @var{n})\n\ +Return a matrix with @var{m} rows and @var{n} columns whose elements are\n\ +taken from the Galois array @var{a}. To decide how to order the elements,\n\ +Octave pretends that the elements of a matrix are stored in column-major\n\ +order (like Fortran arrays are stored).\n\ +\n\ +For example,\n\ +\n\ +@example\n\ +greshape (gf ([1, 2, 3, 4], 3), 2, 2)\n\ +ans =\n\ +GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)\n\ +\n\ +Array elements =\n\ +\n\ + 1 3\n\ + 2 4\n\ +\n\ +@end example\n\ +\n\ +The @code{greshape} function is equivalent to\n\ +\n\ +@example\n\ +@group\n\ +retval = gf (zeros (m, n), a.m, a.prim_poly);\n\ +retval(:) = a;\n\ +@end group\n\ +@end example\n\ +\n\ +@noindent\n\ +but it is somewhat less cryptic to use @code{reshape} instead of the\n\ +colon operator. Note that the total number of elements in the original\n\ +matrix must match the total number of elements in the new matrix.\n\ +@seealso{reshape, :}\n\ +@end deftypefn") +{ + octave_value retval; + int nargin = args.length (); + + if (nargin != 2 && nargin != 3) + { + print_usage (); + } + else + { + int mr = 0, mc = 0; + + if ((!galois_type_loaded) || (args(0).type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("greshape", args(0)); + return retval; + } + galois a = ((const octave_galois&) args(0).get_rep ()).galois_value (); + + if (nargin == 2) + { + RowVector tmp = args(1).row_vector_value (); + mr = (int)tmp(0); + mc = (int)tmp(1); + } + else if (nargin == 3) + { + mr = args(1).nint_value (); + mc = args(2).nint_value (); + } + + int nr = a.rows (); + int nc = a.cols (); + if ((nr * nc) != (mr * mc)) + error ("greshape: sizes must match"); + else + { + RowVector tmp1 (mr*mc); + for (int i = 0; i < nr; i++) + for (int j = 0; j < nc; j++) + tmp1(i+j*nr) = (double)a(i, j); + galois tmp2 (mr, mc, 0, a.m (), a.primpoly ()); + for (int i = 0; i < mr; i++) + for (int j = 0; j < mc; j++) + tmp2(i, j) = (int)tmp1(i+j*mr); + retval = new octave_galois (tmp2); + } + } + return retval; +} + +/* +%% Test input validation +%!error greshape () +%!error greshape (1) +%!error greshape (1, 2, 3, 4) +*/ + +#define DATA_REDUCTION(FCN) \ + \ + octave_value_list retval; \ + \ + int nargin = args.length (); \ + \ + if (nargin == 1 || nargin == 2) \ + { \ + octave_value arg = args(0); \ + \ + int dim = (nargin == 1 ? -1 : args(1).int_value (true) - 1); \ + \ + if (! error_state) \ + { \ + if (dim <= 1 && dim >= -1) \ + { \ + if (galois_type_loaded && (arg.type_id () == \ + octave_galois::static_type_id ())) \ + { \ + galois tmp = ((const octave_galois&)arg.get_rep ()).galois_value (); \ + \ + if (! error_state) \ + retval(0) = new octave_galois (tmp.FCN (dim)); \ + } \ + else \ + { \ + err_wrong_type_arg (#FCN, arg); \ + return retval; \ + } \ + } \ + else \ + error (#FCN ": invalid dimension argument = %d", dim + 1); \ + } \ + } \ + else \ + print_usage (); \ + \ + return retval + +// PKG_ADD: autoload ("gprod", "gf.oct"); +// PKG_DEL: autoload ("gprod", "gf.oct", "remove"); +DEFUN_DLD (gprod, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} gprod (@var{x}, @var{dim})\n\ +Product of elements along dimension @var{dim} of Galois array. If\n\ +@var{dim} is omitted, it defaults to 1 (column-wise products).\n\ +@seealso{prod}\n\ +@end deftypefn") +{ + DATA_REDUCTION (prod); +} + +/* +%% Test input validation +%!error gprod () +%!error gprod (1, 2, 3) +*/ + +// PKG_ADD: autoload ("gsum", "gf.oct"); +// PKG_DEL: autoload ("gsum", "gf.oct", "remove"); +DEFUN_DLD (gsum, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} gsum (@var{x}, @var{dim})\n\ +Sum of elements along dimension @var{dim} of Galois array. If @var{dim}\n\ +is omitted, it defaults to 1 (column-wise sum).\n\ +@seealso{sum}\n\ +@end deftypefn") +{ + DATA_REDUCTION (sum); +} + +/* +%% Test input validation +%!error gsum () +%!error gsum (1, 2, 3) +*/ + +// PKG_ADD: autoload ("gsumsq", "gf.oct"); +// PKG_DEL: autoload ("gsumsq", "gf.oct", "remove"); +DEFUN_DLD (gsumsq, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} gsumsq (@var{x}, @var{dim})\n\ +Sum of squares of elements along dimension @var{dim} of Galois array.\n\ +If @var{dim} is omitted, it defaults to 1 (column-wise sum of squares).\n\ +\n\ +This function is equivalent to computing\n\ +@example\n\ +gsum (x .* conj (x), dim)\n\ +@end example\n\ +but it uses less memory.\n\ +@seealso{sumsq}\n\ +@end deftypefn") +{ + DATA_REDUCTION (sumsq); +} + +/* +%% Test input validation +%!error gsumsq () +%!error gsumsq (1, 2, 3) +*/ + +// PKG_ADD: autoload ("gsqrt", "gf.oct"); +// PKG_DEL: autoload ("gsqrt", "gf.oct", "remove"); +DEFUN_DLD (gsqrt, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} gsqrt (@var{x})\n\ +Compute the square root of @var{x}, element by element, in a Galois Field.\n\ +@seealso{exp}\n\ +@end deftypefn") +{ + octave_value retval; + int nargin = args.length (); + + if (nargin != 1) + { + print_usage (); + return retval; + } + + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("gsqrt", args(0)); + return retval; + } + + galois a = ((const octave_galois&) args(0).get_rep ()).galois_value (); + + retval = new octave_galois (a.sqrt ()); + + return retval; +} + +/* +%% Test input validation +%!error gsqrt () +%!error gsqrt (1, 2) +*/ + +// PKG_ADD: autoload ("glog", "gf.oct"); +// PKG_DEL: autoload ("glog", "gf.oct", "remove"); +DEFUN_DLD (glog, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} glog (@var{x})\n\ +Compute the natural logarithm for each element of @var{x} for a Galois\n\ +array.\n\ +@seealso{log}\n\ +@end deftypefn") +{ + octave_value retval; + int nargin = args.length (); + + if (nargin != 1) + { + print_usage (); + return retval; + } + + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("glog", args(0)); + return retval; + } + + galois a = ((const octave_galois&) args(0).get_rep ()).galois_value (); + + retval = new octave_galois (a.log ()); + + return retval; +} + +/* +%% Test input validation +%!error glog () +%!error glog (1, 2) +*/ + +// PKG_ADD: autoload ("gexp", "gf.oct"); +// PKG_DEL: autoload ("gexp", "gf.oct", "remove"); +DEFUN_DLD (gexp, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} gexp (@var{x})\n\ +Compute the anti-logarithm for each element of @var{x} for a Galois\n\ +array.\n\ +@seealso{exp}\n\ +@end deftypefn") +{ + octave_value retval; + int nargin = args.length (); + + if (nargin != 1) + { + print_usage (); + return retval; + } + + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("gexp", args(0)); + return retval; + } + + galois a = ((const octave_galois&) args(0).get_rep ()).galois_value (); + + retval = new octave_galois (a.exp ()); + + return retval; +} + +/* +%% Test input validation +%!error gexp () +%!error gexp (1, 2) +*/ + +static inline int +modn (int x, int m, int n) +{ + while (x >= n) + { + x -= n; + x = (x >> m) + (x & n); + } + return x; +} + +galois +filter (galois& b, galois& a, galois& x, galois& si) +{ + int ab_len = (a.numel () > b.numel () ? a.numel () : b.numel ()); + b.resize (dim_vector (ab_len, 1), 0); + galois retval (x.numel (), 1, 0, b.m (), b.primpoly ()); + int norm = a(0, 0); + + if (norm == 0) + { + error ("gfilter: the first element of a must be non-zero"); + return galois (); + } + if (si.numel () != ab_len - 1) + { + error ("gfilter: si must be a vector of length max(numel(a), numel(b)) - 1"); + return galois (); + } + if (norm != 1) + { + int idx_norm = b.index_of (norm); + for (int i = 0; i < b.numel (); i++) + { + if (b(i, 0) != 0) + b(i, 0) = b.alpha_to (modn (b.index_of (b(i, 0))-idx_norm+b.n (), + b.m (), b.n ())); + } + } + if (a.numel () > 1) + { + a.resize (dim_vector (ab_len, 1), 0); + + if (norm != 1) + { + int idx_norm = a.index_of (norm); + for (int i = 0; i < a.numel (); i++) + if (a(i, 0) != 0) + a(i, 0) = a.alpha_to (modn (a.index_of (a(i, 0))-idx_norm+a.n (), + a.m (), a.n ())); + } + + for (int i = 0; i < x.numel (); i++) + { + retval(i, 0) = si(0, 0); + if ((b(0, 0) != 0) && (x(i, 0) != 0)) + retval(i, 0) ^= b.alpha_to (modn (b.index_of (b(0, 0)) + + b.index_of (x(i, 0)), b.m (), b.n ())); + if (si.numel () > 1) + { + for (int j = 0; j < si.numel () - 1; j++) + { + si(j, 0) = si(j+1, 0); + if ((a(j+1, 0) != 0) && (retval(i, 0) != 0)) + si(j, 0) ^= a.alpha_to (modn (a.index_of (a(j+1, 0)) + + a.index_of (retval(i, 0)), a.m (), a.n ())); + if ((b(j+1, 0) != 0) && (x(i, 0) != 0)) + si(j, 0) ^= b.alpha_to (modn (b.index_of (b(j+1, 0)) + + b.index_of (x(i, 0)), b.m (), b.n ())); + } + si(si.numel ()-1, 0) = 0; + if ((a(si.numel (), 0) != 0) && (retval(i, 0) != 0)) + si(si.numel ()-1, 0) ^= a.alpha_to (modn (a.index_of (a(si.numel (), 0)) + + a.index_of (retval(i, 0)), + a.m (), a.n ())); + if ((b(si.numel (), 0) != 0) && (x(i, 0) != 0)) + si(si.numel ()-1, 0) ^= b.alpha_to (modn (b.index_of (b(si.numel (), 0)) + + b.index_of (x(i, 0)), + b.m (), b.n ())); + } + else + { + si(0, 0) = 0; + if ((a(1, 0) != 0) && (retval(i, 0) != 0)) + si(0, 0) ^= a.alpha_to (modn (a.index_of (a(1, 0))+ + a.index_of (retval(i, 0)), a.m (), a.n ())); + if ((b(1, 0) != 0) && (x(i, 0) != 0)) + si(0, 0) ^= b.alpha_to (modn (b.index_of (b(1, 0))+ + b.index_of (x(i, 0)), b.m (), b.n ())); + } + } + } + else if (si.numel () > 0) + { + for (int i = 0; i < x.numel (); i++) + { + retval(i, 0) = si(0, 0); + if ((b(0, 0) != 0) && (x(i, 0) != 0)) + retval(i, 0) ^= b.alpha_to (modn (b.index_of (b(0, 0)) + + b.index_of (x(i, 0)), b.m (), b.n ())); + if (si.numel () > 1) + { + for (int j = 0; j < si.numel () - 1; j++) + { + si(j, 0) = si(j+1, 0); + if ((b(j+1, 0) != 0) && (x(i, 0) != 0)) + si(j, 0) ^= b.alpha_to (modn (b.index_of (b(j+1, 0)) + + b.index_of (x(i, 0)), b.m (), b.n ())); + } + si(si.numel ()-1, 0) = 0; + if ((b(si.numel (), 0) != 0) && (x(i, 0) != 0)) + si(si.numel ()-1, 0) ^= b.alpha_to (modn (b.index_of (b(si.numel (), 0)) + + b.index_of (x(i, 0)), + b.m (), b.n ())); + } + else + { + si(0, 0) = 0; + if ((b(1, 0) != 0) && (x(i, 0) != 0)) + si(0, 0) ^= b.alpha_to (modn (b.index_of (b(1, 0)) + + b.index_of (x(i, 0)), b.m (), b.n ())); + } + } + } + else + for (int i = 0; i < x.numel (); i++) + if ((b(0, 0) != 0) && (x(i, 0) != 0)) + retval(i, 0) = b.alpha_to (modn (b.index_of (b(0, 0)) + + b.index_of (x(i, 0)), b.m (), b.n ())); + + return retval; +} + + +// PKG_ADD: autoload ("gfilter", "gf.oct"); +// PKG_DEL: autoload ("gfilter", "gf.oct", "remove"); +DEFUN_DLD (gfilter, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {y =} gfilter (@var{b}, @var{a}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} gfilter (@var{b}, @var{a}, @var{x}, @var{si})\n\ +Digital filtering of vectors in a Galois Field. Returns the solution to\n\ +the following linear, time-invariant difference equation over a Galois\n\ +Field:\n\ +@tex\n\ +$$\n\ +\\sum_{k=0}^N a_{k+1} y_{n-k} = \\sum_{k=0}^M b_{k+1} x_{n-k}, \\qquad\n\ + 1 \\le n \\le P\n\ +$$\n\ +@end tex\n\ +@ifnottex\n\ +\n\ +@smallexample\n\ +@group\n\ + N M\n\ + SUM a(k+1) y(n-k) = SUM b(k+1) x(n-k) for 1<=n<=numel(x)\n\ + k=0 k=0\n\ +@end group\n\ +@end smallexample\n\ +@end ifnottex\n\ +\n\ +@noindent\n\ +where\n\ +@tex\n\ + $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\ +@end tex\n\ +@ifnottex\n\ + N=numel(a)-1 and M=numel(b)-1.\n\ +@end ifnottex\n\ +An equivalent form of this equation is:\n\ +@tex\n\ +$$\n\ +y_n = -\\sum_{k=1}^N c_{k+1} y_{n-k} + \\sum_{k=0}^M d_{k+1} x_{n-k}, \\qquad\n\ + 1 \\le n \\le P\n\ +$$\n\ +@end tex\n\ +@ifnottex\n\ +\n\ +@smallexample\n\ +@group\n\ + N M\n\ + y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k) for 1<=n<=numel(x)\n\ + k=1 k=0\n\ +@end group\n\ +@end smallexample\n\ +@end ifnottex\n\ +\n\ +@noindent\n\ +where\n\ +@tex\n\ +$c = a/a_1$ and $d = b/a_1$.\n\ +@end tex\n\ +@ifnottex\n\ + c = a/a(1) and d = b/a(1).\n\ +@end ifnottex\n\ +\n\ +If the fourth argument @var{si} is provided, it is taken as the initial\n\ +state of the system and the final state is returned as @var{sf}. The\n\ +state vector is a column vector whose length is equal to the length of\n\ +the longest coefficient vector minus one. If @var{si} is not supplied,\n\ +the initial state vector is set to all zeros.\n\ +@seealso{filter}\n\ +@end deftypefn") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 3 || nargin > 4) + { + print_usage (); + return retval; + } + + if (!galois_type_loaded) + { + error ("gfilter: wrong argument types"); + return retval; + } + + bool x_is_row_vector = (args(2).rows () == 1); + bool si_is_row_vector = (nargin == 4 && args(3).rows () == 1); + galois b, a, x, si; + bool ib=false, ia=false, ix = false, isi=false; + + if (args(0).type_id () == octave_galois::static_type_id ()) + { + b = ((const octave_galois&) args(0).get_rep ()).galois_value (); + ib = true; + } + if (args(1).type_id () == octave_galois::static_type_id ()) + { + a = ((const octave_galois&) args(1).get_rep ()).galois_value (); + ia = true; + } + if (args(2).type_id () == octave_galois::static_type_id ()) + { + x = ((const octave_galois&) args(2).get_rep ()).galois_value (); + ix = true; + } + if (nargin == 4) + { + if (args(3).type_id () == octave_galois::static_type_id ()) + { + si = ((const octave_galois&) args(3).get_rep ()).galois_value (); + isi = true; + } + } + + if (!ib && !ia && !ix && !isi) + { + error ("gfilter: wrong argument types"); + return retval; + } + + if (!ib) + { + if (ia) + b = galois (args(0).matrix_value (), a.m (), a.primpoly ()); + else if (ix) + b = galois (args(0).matrix_value (), x.m (), x.primpoly ()); + else if (isi) + b = galois (args(0).matrix_value (), si.m (), si.primpoly ()); + } + if (!ia) + a = galois (args(1).matrix_value (), b.m (), b.primpoly ()); + if (!ix) + x = galois (args(2).matrix_value (), b.m (), b.primpoly ()); + + if (nargin == 4) + { + if (!isi) + si = galois (args(3).matrix_value (), b.m (), b.primpoly ()); + } + else + { + int a_len = a.numel (); + int b_len = b.numel (); + + int si_len = (a_len > b_len ? a_len : b_len) - 1; + + si = galois (si_len, 1, 0, b.m (), b.primpoly ()); + } + + if ((b.m () != a.m ()) || (b.m () != x.m ()) || (b.m () != si.m ()) || + (b.primpoly () != a.primpoly ()) || (b.primpoly () != x.primpoly ()) || + (b.primpoly () != si.primpoly ())) + { + error ("gfilter: arguments must be in same galois field"); + return retval; + } + + if (b.cols () > 1) + b = b.transpose (); + if (a.cols () > 1) + a = a.transpose (); + if (x.cols () > 1) + x = x.transpose (); + if (si.cols () > 1) + si = si.transpose (); + + if (b.cols () > 1 || a.cols () > 1 || x.cols () > 1 || si.cols () > 1) + { + error ("gfilter: arguments must be vectors"); + return retval; + } + + galois y (filter (b, a, x, si)); + if (nargout == 2) + { + if (si_is_row_vector) + retval(1) = new octave_galois (si.transpose ()); + else + retval(1) = new octave_galois (si); + } + + if (x_is_row_vector) + retval(0) = new octave_galois (y.transpose ()); + else + retval(0) = new octave_galois (y); + + return retval; +} + +/* +%% Test input validation +%!error gfilter () +%!error gfilter (1) +%!error gfilter (1, 2) +%!error gfilter (1, 2, 3, 4, 5) +*/ + +// PKG_ADD: autoload ("glu", "gf.oct"); +// PKG_DEL: autoload ("glu", "gf.oct", "remove"); +DEFUN_DLD (glu, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{l}, @var{u}, @var{p}] =} glu (@var{a})\n\ +@cindex LU decomposition of Galois matrix\n\ +Compute the LU decomposition of @var{a} in a Galois Field. The result is\n\ +returned in a permuted form, according to the optional return value\n\ +@var{p}. For example, given the matrix\n\ +@code{a = gf ([1, 2; 3, 4], 3)},\n\ +\n\ +@example\n\ +[l, u, p] = glu (a)\n\ +@end example\n\ +\n\ +@noindent\n\ +returns\n\ +\n\ +@example\n\ +l =\n\ +GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)\n\ +\n\ +Array elements =\n\ +\n\ + 1 0\n\ + 6 1\n\ +\n\ +u =\n\ +GF(2^3) array. Primitive Polynomial = D^3+D+1 (decimal 11)\n\ +\n\ +Array elements =\n\ +\n\ + 3 4\n\ + 0 7\n\ +\n\ +p =\n\ +\n\ +Permutation Matrix\n\ +\n\ + 0 1\n\ + 1 0\n\ +\n\ +@end example\n\ +\n\ +Such that @code{@var{p} * @var{a} = @var{l} * @var{u}}. If the argument\n\ +@var{p} is not included then the permutations are applied to @var{l}\n\ +so that @code{@var{a} = @var{l} * @var{u}}. @var{l} is then a pseudo-\n\ +lower triangular matrix. The matrix @var{a} can be rectangular.\n\ +@seealso{lu}\n\ +@end deftypefn") +{ + octave_value_list retval; + + + int nargin = args.length (); + + if (nargin != 1 || nargout > 3) + { + print_usage (); + return retval; + } + + octave_value arg = args(0); + + if (!galois_type_loaded || (arg.type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("glu", arg); + return retval; + } + + galois m = ((const octave_galois&) arg.get_rep ()).galois_value (); + + int nr = arg.rows (); + int nc = arg.columns (); + + int arg_is_empty = arg.isempty (); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + { + retval(0) = new octave_galois (galois (0, 0, 0, m.m (), m.primpoly ())); + retval(1) = new octave_galois (galois (0, 0, 0, m.m (), m.primpoly ())); + retval(2) = new octave_galois (galois (0, 0, 0, m.m (), m.primpoly ())); + return retval; + } + + if (! error_state) + { + galoisLU fact (m); + + switch (nargout) + { + case 0: + case 1: + case 2: + { + // While we don't have sparse galois matrices converting the + // permutation matrix to a full matrix is the best we can do. + Matrix P = Matrix (fact.P ()); + galois L = P.transpose () * fact.L (); + retval(1) = new octave_galois (fact.U ()); + retval(0) = new octave_galois (L); + } + break; + + case 3: + default: + retval(2) = fact.P (); + retval(1) = new octave_galois (fact.U ()); + retval(0) = new octave_galois (fact.L ()); + break; + } + } + + return retval; +} + +/* +%% Test input validation +%!error glu () +%!error glu (1, 2) +*/ + +// PKG_ADD: autoload ("ginv", "gf.oct"); +// PKG_DEL: autoload ("ginv", "gf.oct", "remove"); +DEFUN_DLD (ginv, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{x}, @var{rcond}] =} ginv (@var{a})\n\ +Compute the inverse of the square matrix @var{a}. Return an estimate\n\ +of the reciprocal condition number if requested, otherwise warn of an\n\ +ill-conditioned matrix if the reciprocal condition number is small.\n\ +@seealso{inv}\n\ +@end deftypefn") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage (); + return retval; + } + + octave_value arg = args(0); + + int nr = arg.rows (); + int nc = arg.columns (); + + if (!galois_type_loaded || (arg.type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("ginverse", arg); + return retval; + } + + galois m = ((const octave_galois&) arg.get_rep ()).galois_value (); + + int arg_is_empty = arg.isempty (); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + { + retval(0) = new octave_galois (galois (0, 0, 0, m.m (), m.primpoly ())); + return retval; + } + if (nr != nc) + { + err_square_matrix_required ("ginverse", "X"); + return retval; + } + + if (! error_state) + { + int info; + double rcond = 0.0; + + galois result = m.inverse (info, 1); + + if (nargout > 1) + retval(1) = rcond; + + retval(0) = new octave_galois (result); + + if (nargout < 2 && info == -1) + warning ("inverse: matrix singular to machine precision, rcond = %g", rcond); + } + + return retval; +} + +/* +%% Test input validation +%!error ginv () +%!error ginv (1, 2) +*/ + +// FIXME: this should really be done with an alias, but +// alias_builtin() won't do the right thing if we are actually using +// dynamic linking. + +// PKG_ADD: autoload ("ginverse", "gf.oct"); +// PKG_DEL: autoload ("ginverse", "gf.oct", "remove"); +DEFUN_DLD (ginverse, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} ginverse (@var{a})\n\ +Compute the inverse of the square matrix @var{a}. Return an estimate\n\ +of the reciprocal condition number if requested, otherwise warn of an\n\ +ill-conditioned matrix if the reciprocal condition number is small.\n\ +@seealso{ginv}\n\ +@end deftypefn") +{ + return Fginv (args, nargout); +} + +/* +%% Test input validation +%!error ginverse () +%!error ginverse (1, 2) +*/ + +// PKG_ADD: autoload ("gdet", "gf.oct"); +// PKG_DEL: autoload ("gdet", "gf.oct", "remove"); +DEFUN_DLD (gdet, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{d} =} gdet (@var{a})\n\ +Compute the determinant of the Galois array @var{a}.\n\ +@seealso{det}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage (); + return retval; + } + + octave_value arg = args(0); + + if (!galois_type_loaded || (arg.type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("gdet", arg); + return retval; + } + + int nr = arg.rows (); + int nc = arg.columns (); + + galois m = ((const octave_galois&) arg.get_rep ()).galois_value (); + + int arg_is_empty = arg.isempty (); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + { + retval = new octave_galois (galois (1, 1, 1, m.m (), m.primpoly ())); + return retval; + } + + if (nr != nc) + { + err_square_matrix_required ("det", "A"); + return retval; + } + + retval = new octave_galois (m.determinant ()); + return retval; +} + +/* +%% Test input validation +%!error gdet () +%!error gdet (1, 2) +*/ + +// PKG_ADD: autoload ("grank", "gf.oct"); +// PKG_DEL: autoload ("grank", "gf.oct", "remove"); +DEFUN_DLD (grank, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{d} =} grank (@var{a})\n\ +Compute the rank of the Galois array @var{a} by counting the independent\n\ +rows and columns.\n\ +@seealso{rank}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage (); + return retval; + } + + octave_value arg = args(0); + + if (!galois_type_loaded || (arg.type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("grank", arg); + return retval; + } + + int nr = arg.rows (); + int nc = arg.columns (); + + galois m = ((const octave_galois&) arg.get_rep ()).galois_value (); + + int arg_is_empty = arg.isempty (); + + if (arg_is_empty > 0) + retval = 0.0; + else if (arg_is_empty == 0) + { + int d = 0; + int mm = m.m (); + int mn = m.n (); + OCTAVE_LOCAL_BUFFER (int, ci, nr); + + for (int i = 0; i < nc; i++) + { + int idx = -1; + int iel = 0; + for (int j = 0; j < nr; j++) + { + ci[j] = m.elem (j, i); + if (ci[j] != 0 && idx == -1) + { + iel = ci[j]; + idx = j; + } + } + + if (idx != -1) + { + d++; + int indx = m.index_of (iel); + for (int j = 0; j < nr; j++) + if (ci[j] != 0) + ci[j] = m.alpha_to (modn (m.index_of (ci[j]) - indx + mn, mm, mn)); + + for (int j = i+1; j < nc; j++) + { + if (m.elem (idx, j) != 0) + { + indx = m.index_of (m.elem (idx, j)); + for (int k = 0; k < nr; k++) + if (ci[k] != 0) + m.elem (k, j) ^= m.alpha_to (modn (m.index_of (ci[k]) + indx + + mn, mm, mn)); + } + } + } + } + retval = (double)d; + } + return retval; +} + +/* +%% Test input validation +%!error grank () +%!error grank (1, 2) +*/ + +// PKG_ADD: autoload ("rsenc", "gf.oct"); +// PKG_DEL: autoload ("rsenc", "gf.oct", "remove"); +DEFUN_DLD (rsenc, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{code} =} rsenc (@var{msg}, @var{n}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{code} =} rsenc (@var{msg}, @var{n}, @var{k}, @var{g})\n\ +@deftypefnx {Loadable Function} {@var{code} =} rsenc (@var{msg}, @var{n}, @var{k}, @var{fcr}, @var{prim})\n\ +@deftypefnx {Loadable Function} {@var{code} =} rsenc (@dots{}, @var{parpos})\n\ +Encodes the message @var{msg} using a [@var{n},@var{k}] Reed-Solomon coding.\n\ +The variable @var{msg} is a Galois array with @var{k} columns and an arbitrary\n\ +number of rows. Each row of @var{msg} represents a single block to be coded\n\ +by the Reed-Solomon coder. The coded message is returned in the Galois\n\ +array @var{code} containing @var{n} columns and the same number of rows as\n\ +@var{msg}.\n\ +\n\ +The use of @code{rsenc} can be seen in the following short example.\n\ +\n\ +@example\n\ +m = 3; n = 2^m -1; k = 3;\n\ +msg = gf ([1 2 3; 4 5 6], m);\n\ +code = rsenc (msg, n, k);\n\ +@end example\n\ +\n\ +If @var{n} does not equal @code{2^@var{m}-1}, where m is an integer, then a\n\ +shorten Reed-Solomon coding is used where zeros are added to the start of\n\ +each row to obtain an allowable codeword length. The returned @var{code}\n\ +has these prepending zeros stripped.\n\ +\n\ +By default the generator polynomial used in the Reed-Solomon coding is based\n\ +on the properties of the Galois Field in which @var{msg} is given. This\n\ +default generator polynomial can be overridden by a polynomial in @var{g}.\n\ +Suitable generator polynomials can be constructed with @code{rsgenpoly}.\n\ +@var{fcr} is an integer value, and it is taken to be the first consecutive\n\ +root of the generator polynomial. The variable @var{prim} is then the\n\ +primitive element used to construct the generator polynomial, such that\n\ +@tex\n\ +$g = (x - A^b) (x - A^{b+p}) \\cdots (x - A ^{b+2tp-1})$.\n\ +@end tex\n\ +@ifnottex\n\ +\n\ +@var{g} = (@var{x} - A^@var{b}) * (@var{x} - A^(@var{b}+@var{prim})) * ... * (@var{x} - A^(@var{b}+2*@var{t}*@var{prim}-1)).\n\ +@end ifnottex\n\ +\n\ +where @var{b} is equal to @code{@var{fcr} * @var{prim}}. By default @var{fcr}\n\ +and @var{prim} are both 1.\n\ +\n\ +By default the parity symbols are placed at the end of the coded message.\n\ +The variable @var{parpos} controls this positioning and can take the values\n\ +@code{\"beginning\"} or @code{\"end\"}.\n\ +@seealso{gf, rsdec, rsgenpoly}\n\ +@end deftypefn") +{ + octave_value retval; + int nargin = args.length (); + + if (nargin < 3 || nargin > 5) + { + print_usage (); + return retval; + } + + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("rsenc", args(0)); + return retval; + } + + galois msg = ((const octave_galois&) args(0).get_rep ()).galois_value (); + int nsym = msg.rows (); + int primpoly = msg.primpoly (); + int n = args(1).nint_value (); + int k = args(2).nint_value (); + + int m = 1; + while (n > (1<<m)) + m++; + int nn = (1<<m) - 1; + + if (msg.cols () != k) + { + error ("rsenc: message contains incorrect number of symbols"); + return retval; + } + + if (msg.m () != m) + { + error ("rsenc: message in incorrect galois field for codeword length"); + return retval; + } + + if ((n < 3) || (n < k) || (m > __OCTAVE_GALOIS_MAX_M)) + { + error ("rsenc: invalid values of message and codeword length"); + return retval; + } + + if ((n-k) & 1) + { + error ("rsenc: difference of message and codeword length must be even"); + return retval; + } + + int nroots = n-k; + galois genpoly; + bool have_genpoly = false; + bool parity_at_end = true; + int fcr = 0; + int prim = 0; + + for (int i = 3; i < nargin; i++) + { + if (args(i).is_string ()) + { + std::string parstr = args(i).string_value (); + for (int j = 0; j < (int)parstr.length (); j++) + parstr[j] = toupper (parstr[j]); + + if (!parstr.compare("END")) + { + parity_at_end = true; + } + else if (!parstr.compare("BEGINNING")) + { + parity_at_end = false; + } + else + { + error ("rsenc: unrecoginized parity position"); + return retval; + } + } + else + { + if (args(i).type_id () == octave_galois::static_type_id ()) + { + if (have_genpoly) + { + print_usage (); + return retval; + } + genpoly = ((const octave_galois&) args(i).get_rep ()).galois_value (); + + if (genpoly.cols () > genpoly.rows ()) + genpoly = genpoly.transpose (); + } + else + { + if (have_genpoly) + { + if (prim != 0) + { + print_usage (); + return retval; + } + prim = args(i).nint_value (); + } + else + fcr = args(i).nint_value (); + } + have_genpoly = true; + } + } + + if ((genpoly.rows () == 0) || (genpoly.cols () == 0)) + { + if (fcr == 0) + fcr = 1; + if (prim == 0) + prim = 1; + + // Create polynomial of right length. + genpoly = galois (nroots+1, 1, 0, m, primpoly); + + genpoly(nroots, 0) = 1; + int i, root; + for (i = 0, root=fcr*prim; i < nroots; i++, root += prim) + { + genpoly(nroots-i-1, 0) = 1; + + // Multiply genpoly by @**(root + x) + for (int j = i; j > 0; j--) + { + int k = nroots - j; + if (genpoly(k, 0) != 0) + genpoly(k, 0) = genpoly(k+1, 0) + ^ genpoly.alpha_to (modn (genpoly.index_of (genpoly(k, 0)) + + root, m, n)); + else + genpoly(k, 0) = genpoly(k+1, 0); + } + // genpoly(nroots,0) can never be zero + genpoly(nroots, 0) = genpoly.alpha_to (modn (genpoly.index_of (genpoly(nroots, 0)) + + root, m, n)); + } + + } + else + { + if (genpoly.cols () != 1) + { + error ("rsenc: the generator polynomial must be a vector"); + return retval; + } + + if (genpoly.primpoly () != primpoly) + { + error ("rsenc: the generator polynomial must be same galois field " + "as the message"); + return retval; + } + + if (genpoly.rows () != nroots+1) + { + error ("rsenc: generator polynomial has incorrect order"); + return retval; + } + } + + int norm = genpoly(0, 0); + + // Take logarithm of generator polynomial, for faster coding + for (int i = 0; i < nroots+1; i++) + genpoly(i, 0) = genpoly.index_of (genpoly(i, 0)); + + // Add space for parity block + msg.resize (dim_vector (nsym, n), 0); + + // The code below basically finds the parity bits by treating the + // message as a polynomial and dividing it by the generator polynomial. + // The parity bits are then the remainder of this division. If the parity + // is at the end the polynomial is treat MSB first, otherwise it is + // treated LSB first + // + // This code could just as easily be written as + // [ignore par] = gdeconv(msg, genpoly); + // But the code below has the advantage of being 20 times faster :-) + + if (parity_at_end) + { + for (int l = 0; l < nsym; l++) + { + galois par (nroots, 1, 0, m, primpoly); + for (int i = 0; i < k; i++) + { + int feedback = par.index_of (par(0, 0) ^ msg(l, i)); + if (feedback != nn) + { + if (norm != 1) + feedback = modn (nn-genpoly(0, 0)+feedback, m, nn); + for (int j = 1; j < nroots; j++) + par(j, 0) ^= par.alpha_to (modn (feedback + + genpoly(j, 0), m, nn)); + } + for (int j = 1; j < nroots; j++) + par(j-1, 0) = par(j, 0); + if (feedback != nn) + par(nroots-1, 0) = par.alpha_to (modn (feedback+ + genpoly(nroots, 0), m, nn)); + else + par(nroots-1, 0) = 0; + } + for (int j = 0; j < nroots; j++) + msg(l, k+j) = par(j, 0); + } + } + else + { + for (int l = 0; l < nsym; l++) + { + for (int i=k; i > 0; i--) + msg(l, i+nroots-1) = msg(l, i-1); + for (int i = 0; i<nroots; i++) + msg(l, i) = 0; + } + for (int l = 0; l < nsym; l++) + { + galois par (nroots, 1, 0, m, primpoly); + for (int i = n; i > nroots; i--) + { + int feedback = par.index_of (par(0, 0) ^ msg(l, i-1)); + if (feedback != nn) + { + if (norm != 1) + feedback = modn (nn-genpoly(0, 0)+feedback, m, nn); + for (int j = 1; j < nroots; j++) + par(j, 0) ^= par.alpha_to (modn (feedback + + genpoly(j, 0), m, nn)); + } + for (int j = 1; j < nroots; j++) + par(j-1, 0) = par(j, 0); + if (feedback != nn) + par(nroots-1, 0) = par.alpha_to (modn (feedback+ + genpoly(nroots, 0), m, nn)); + else + par(nroots-1, 0) = 0; + } + for (int j = 0; j < nroots; j++) + msg(l, j) = par(nroots-j-1, 0); + } + } + + retval = new octave_galois (msg); + + return retval; +} + +/* +%% Test input validation +%!error rsenc () +%!error rsenc (1) +%!error rsenc (1, 2) +%!error rsenc (1, 2, 3, 4, 5, 6) +*/ + +int +decode_rs(galois& data, const int prim, const int iprim, const int nroots, + const int fcr, const int drow, const bool msb_first) +{ + int deg_lambda, el, deg_omega; + int i, j, r, k; + int q, tmp, num1, num2, den, discr_r; + int syn_error, count; + int m = data.m (); + int n = data.n (); + int A0 = n; + + /* Err Locator and syndrome poly */ + OCTAVE_LOCAL_BUFFER (int, lambda, nroots+1); + OCTAVE_LOCAL_BUFFER (int, s, nroots); + + OCTAVE_LOCAL_BUFFER (int, b, nroots+1); + OCTAVE_LOCAL_BUFFER (int, t, nroots+1); + OCTAVE_LOCAL_BUFFER (int, omega, nroots+1); + + OCTAVE_LOCAL_BUFFER (int, root, nroots); + OCTAVE_LOCAL_BUFFER (int, reg, nroots+1); + OCTAVE_LOCAL_BUFFER (int, loc, nroots); + + /* form the syndromes; i.e., evaluate data(x) at roots of g(x) */ + if (msb_first) + { + for (i = 0; i < nroots; i++) + s[i] = data(drow, 0); + + for (j = 1; j < n; j++) + for (i = 0; i<nroots; i++) + if(s[i] == 0) + s[i] = data(drow, j); + else + s[i] = data(drow, j) ^ data.alpha_to (modn (data.index_of (s[i]) + + (fcr+i)*prim, m, n)); + } + else + { + for (i = 0; i<nroots; i++) + s[i] = data(drow, n-1); + + for (j = n-1; j>0; j--) + for (i = 0; i < nroots; i++) + if(s[i] == 0) + s[i] = data(drow, j-1); + else + s[i] = data(drow, j-1) ^ data.alpha_to (modn (data.index_of (s[i]) + + (fcr+i)*prim, m, n)); + } + + /* Convert syndromes to index form, checking for nonzero condition */ + syn_error = 0; + for (i = 0; i < nroots; i++) + { + syn_error |= s[i]; + s[i] = data.index_of (s[i]); + } + + if (!syn_error) + /* if syndrome is zero, data(drow,:) is a codeword and there are no + * errors to correct. So return data(drow,:) unmodified + */ + return 0; + + memset(&lambda[1], 0, nroots*sizeof (lambda[0])); + lambda[0] = 1; + + for (i = 0; i < nroots+1; i++) + b[i] = data.index_of (lambda[i]); + + /* + * Begin Berlekamp-Massey algorithm to determine error locator polynomial + */ + r = 0; + el = 0; + while (++r <= nroots) + {/* r is the step number */ + /* Compute discrepancy at the r-th step in poly-form */ + discr_r = 0; + for (i = 0; i < r; i++) + { + if ((lambda[i] != 0) && (s[r-i-1] != A0)) + { + discr_r ^= data.alpha_to (modn (data.index_of (lambda[i]) + + s[r-i-1], m, n)); + } + } + discr_r = data.index_of (discr_r); /* Index form */ + if (discr_r == A0) + { + /* 2 lines below: B(x) <-- x*B(x) */ + memmove(&b[1], b, nroots*sizeof (b[0])); + b[0] = A0; + } + else + { + /* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */ + t[0] = lambda[0]; + for (i = 0 ; i < nroots; i++) + { + if(b[i] != A0) + t[i+1] = lambda[i+1] ^ data.alpha_to (modn (discr_r + b[i], m, n)); + else + t[i+1] = lambda[i+1]; + } + if (2 * el <= r - 1) + { + el = r - el; + /* + * 2 lines below: B(x) <-- inv(discr_r) * + * lambda(x) + */ + for (i = 0; i <= nroots; i++) + b[i] = (lambda[i] == 0) ? A0 : modn (data.index_of (lambda[i]) - + discr_r + n, m, n); + } + else + { + /* 2 lines below: B(x) <-- x*B(x) */ + memmove(&b[1], b, nroots*sizeof (b[0])); + b[0] = A0; + } + memcpy(lambda, t, (nroots+1)*sizeof (t[0])); + } + } + + /* Convert lambda to index form and compute deg(lambda(x)) */ + deg_lambda = 0; + for (i = 0; i < nroots+1; i++) + { + lambda[i] = data.index_of (lambda[i]); + if(lambda[i] != A0) + deg_lambda = i; + } + + /* Find roots of the error locator polynomial by Chien search */ + memcpy(®[1], &lambda[1], nroots*sizeof (reg[0])); + count = 0; /* Number of roots of lambda(x) */ + for (i = 1, k = iprim-1; i <= n; i++, k = modn (k+iprim, m, n)) + { + q = 1; /* lambda[0] is always 0 */ + for (j = deg_lambda; j > 0; j--) + { + if (reg[j] != A0) + { + reg[j] = modn (reg[j] + j, m, n); + q ^= data.alpha_to (reg[j]); + } + } + if (q != 0) + continue; /* Not a root */ + /* store root (index-form) and error location number */ + root[count] = i; + loc[count] = k; + /* If we've already found max possible roots, + * abort the search to save time + */ + if(++count == deg_lambda) + break; + } + if (deg_lambda != count) + { + /* + * deg(lambda) unequal to number of roots => uncorrectable + * error detected + */ + return -1; + } + /* + * Compute err evaluator poly omega(x) = s(x)*lambda(x) (modulo + * x**nroots). in index form. Also find deg(omega). + */ + deg_omega = 0; + for (i = 0; i < nroots; i++) + { + tmp = 0; + j = (deg_lambda < i) ? deg_lambda : i; + for (; j >= 0; j--) + { + if ((s[i - j] != A0) && (lambda[j] != A0)) + tmp ^= data.alpha_to (modn (s[i - j] + lambda[j], m, n)); + } + if(tmp != 0) + deg_omega = i; + omega[i] = data.index_of (tmp); + } + omega[nroots] = A0; + + /* + * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 = + * inv(X(l))**(fcr-1) and den = lambda_pr(inv(X(l))) all in poly-form + */ + for (j = count-1; j >= 0; j--) + { + num1 = 0; + for (i = deg_omega; i >= 0; i--) + { + if (omega[i] != A0) + num1 ^= data.alpha_to (modn (omega[i] + i * root[j], m, n)); + } + num2 = data.alpha_to (modn (root[j] * (fcr - 1) + n, m, n)); + den = 0; + + /* lambda[i+1] for i even is the formal deriv lambda_pr of lambda[i] */ + for (i = (deg_lambda < nroots-1 ? deg_lambda : nroots-1) & ~1; i >= 0; + i -=2) + { + if(lambda[i+1] != A0) + den ^= data.alpha_to (modn (lambda[i+1] + i * root[j], m, n)); + } + if (den == 0) + { + count = -1; + break; + } + /* Apply error to data */ + if (num1 != 0) + { + if (msb_first) + data(drow, loc[j]) ^= data.alpha_to (modn (data.index_of (num1) + + data.index_of (num2) + + n - data.index_of (den), + m, n)); + else + data(drow, n-loc[j]-1) ^= data.alpha_to (modn (data.index_of (num1) + + data.index_of (num2) + + n - data.index_of (den), + m, n)); + } + } + + return count; +} + +// PKG_ADD: autoload ("rsdec", "gf.oct"); +// PKG_DEL: autoload ("rsdec", "gf.oct", "remove"); +DEFUN_DLD (rsdec, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{msg} =} rsdec (@var{code}, @var{n}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{msg} =} rsdec (@var{code}, @var{n}, @var{k}, @var{g})\n\ +@deftypefnx {Loadable Function} {@var{msg} =} rsdec (@var{code}, @var{n}, @var{k}, @var{fcr}, @var{prim})\n\ +@deftypefnx {Loadable Function} {@var{msg} =} rsdec (@dots{}, @var{parpos})\n\ +@deftypefnx {Loadable Function} {[@var{msg}, @var{nerr}] =} rsdec (@dots{})\n\ +@deftypefnx {Loadable Function} {[@var{msg}, @var{nerr}, @var{ccode}] =} rsdec (@dots{})\n\ +Decodes the message contained in @var{code} using a [@var{n},@var{k}]\n\ +Reed-Solomon code. The variable @var{code} must be a Galois array with\n\ +@var{n} columns and an arbitrary number of rows. Each row of @var{code}\n\ +represents a single block to be decoded by the Reed-Solomon coder. The\n\ +decoded message is returned in the variable @var{msg} containing @var{k}\n\ +columns and the same number of rows as @var{code}.\n\ +\n\ +If @var{n} does not equal @code{2^@var{m}-1}, where m is an integer, then a\n\ +shorten Reed-Solomon decoding is used where zeros are added to the start of\n\ +each row to obtain an allowable codeword length. The returned @var{msg}\n\ +has these prepending zeros stripped.\n\ +\n\ +By default the generator polynomial used in the Reed-Solomon coding is based\n\ +on the properties of the Galois Field in which @var{msg} is given. This\n\ +default generator polynomial can be overridden by a polynomial in @var{g}.\n\ +Suitable generator polynomials can be constructed with @code{rsgenpoly}.\n\ +@var{fcr} is an integer value, and it is taken to be the first consecutive\n\ +root of the generator polynomial. The variable @var{prim} is then the\n\ +primitive element used to construct the generator polynomial. By default\n\ +@var{fcr} and @var{prim} are both 1. It is significantly faster to specify\n\ +the generator polynomial in terms of @var{fcr} and @var{prim}, since @var{g}\n\ +is converted to this form in any case.\n\ +\n\ +By default the parity symbols are placed at the end of the coded message.\n\ +The variable @var{parpos} controls this positioning and can take the values\n\ +@code{\"beginning\"} or @code{\"end\"}. If the parity symbols are at the end, the message is\n\ +treated with the most-significant symbol first, otherwise the message is\n\ +treated with the least-significant symbol first.\n\ +@seealso{gf, rsenc, rsgenpoly}\n\ +@end deftypefn") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 3 || nargin > 5) + { + print_usage (); + return retval; + } + + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) + { + err_wrong_type_arg ("rsdec", args(0)); + return retval; + } + + galois code = ((const octave_galois&) args(0).get_rep ()).galois_value (); + int nsym = code.rows (); + int primpoly = code.primpoly (); + int n = args(1).nint_value (); + int k = args(2).nint_value (); + + int m = 1; + while (n > (1<<m)) + m++; + int nn = (1<<m) - 1; + + if (code.cols () != n) + { + error ("rsdec: coded message contains incorrect number of symbols"); + return retval; + } + + if (code.m () != m) + { + error ("rsdec: coded message in incorrect galois field for " + "codeword length"); + return retval; + } + + if ((n < 3) || (n < k) || (m > __OCTAVE_GALOIS_MAX_M)) + { + error ("rsdec: invalid values of message and codeword length"); + return retval; + } + + if ((n-k) & 1) + { + error ("rsdec: difference of message and codeword length must be even"); + return retval; + } + + int nroots = n-k; + galois genpoly; + bool have_genpoly = false; + bool parity_at_end = true; + int fcr = 0; + int prim = 0; + int iprim; + + for (int i = 3; i < 6; i++) + { + if (nargin > i) + { + if (args(i).is_string ()) + { + std::string parstr = args(i).string_value (); + for (int j = 0; j < (int)parstr.length (); j++) + parstr[j] = toupper (parstr[j]); + + if (!parstr.compare("END")) + { + parity_at_end = true; + } + else if (!parstr.compare("BEGINNING")) + { + parity_at_end = false; + } + else + { + error ("rsdec: unrecoginized parrity position"); + return retval; + } + } + else + { + if (args(i).type_id () == octave_galois::static_type_id ()) + { + if (have_genpoly) + { + print_usage (); + return retval; + } + genpoly = ((const octave_galois&) args(i).get_rep ()).galois_value (); + } + else + { + if (have_genpoly) + { + if (prim != 0) + { + print_usage (); + return retval; + } + prim = args(i).nint_value (); + } + else + fcr = args(i).nint_value (); + } + have_genpoly = true; + } + } + } + + if (have_genpoly) + { + if (fcr != 0) + { + if ((fcr < 1) || (fcr > nn)) + { + error ("rsdec: invalid first consecutive root of generator polynomial"); + return retval; + } + if ((prim < 1) || (prim > nn)) + { + error ("rsdec: invalid primitive element of generator polynomial"); + return retval; + } + } + else + { + if (genpoly.cols () > genpoly.rows ()) + genpoly = genpoly.transpose (); + + if (genpoly.cols () != 1) + { + error ("rsdec: the generator polynomial must be a vector"); + return retval; + } + + if (genpoly.primpoly () != primpoly) + { + error ("rsdec: the generator polynomial must be same galois " + "field as the message"); + return retval; + } + + if (genpoly.rows () != nroots+1) + { + error ("rsdec: generator polynomial has incorrect order"); + return retval; + } + + // Find the roots of the generator polynomial + int count = 0; + OCTAVE_LOCAL_BUFFER (int, roots, nroots); + for (int j = 0; j <= nn; j++) + { + // Evaluate generator polynomial at j + int val = genpoly(0, 0); + int indx = genpoly.index_of (j); + for (int i = 0; i<nroots; i++) + { + if (val == 0) + val = genpoly(i+1, 0); + else + val = genpoly(i+1, 0) ^ genpoly.alpha_to (modn (indx + + genpoly.index_of (val), + m, nn)); + } + if (val == 0) + { + roots[count] = j; + count++; + if (count == nroots) + break; + } + } + + if (count != nroots) + { + error ("rsdec: generator polynomial can not have repeated roots"); + return retval; + } + + // Logarithm of roots wrt primitive element + for (int i = 0; i < count; i++) + roots[i] = genpoly.index_of (roots[i]); + + // Find a corresponding fcr and prim that coincide with the roots. + // FIXME: This is a naive algorithm and should be improved !!! + bool found = true; + for (fcr = 1; fcr < n+1; fcr++) + { + for (prim = 1; prim < n+1; prim++) + { + found = true; + for (int i = 0; i<nroots; i++) + { + int tmp = modn ((fcr + i)*prim, m, n); + for (int j = 0; j<count; j++) + { + if (tmp == roots[j]) + { + tmp = -1; + break; + } + } + if (tmp != -1) + { + found = false; + break; + } + } + if (found) + break; + } + if (found) + break; + } + } + } + else + { + fcr = 1; + prim = 1; + } + + /* Find prim-th root of 1, used in decoding */ + for (iprim = 1; (iprim % prim) != 0; iprim += n) + ; + iprim = iprim / prim; + + galois msg (nsym, k, 0, m, primpoly); + ColumnVector nerr (nsym, 0); + + if (nn != n) + { + code.resize (dim_vector (nsym, nn), 0); + if (parity_at_end) + for (int l = 0; l < nsym; l++) + for (int i=n; i > 0; i--) + code(l, i+nn-n-1) = code(l, i-1); + } + + for (int l = 0; l < nsym; l++) + nerr(l) = decode_rs (code, prim, iprim, nroots, fcr, l, parity_at_end); + + if (nn != n) + { + if (parity_at_end) + for (int l = 0; l < nsym; l++) + for (int i = 0; i > n; i--) + code(l, i) = code(l, i+nn-n); + code.resize (dim_vector (nsym, n), 0); + } + + if (parity_at_end) + { + for (int l = 0; l < nsym; l++) + for (int i = 0; i < k; i++) + msg(l, i) = code(l, i); + } + else + { + for (int l = 0; l < nsym; l++) + for (int i = 0; i < k; i++) + msg(l, i) = code(l, nroots+i); + } + + retval(0) = new octave_galois (msg); + retval(1) = octave_value (nerr); + retval(2) = new octave_galois (code); + + return retval; +} + +/* +%% Test input validation +%!error rsdec () +%!error rsdec (1) +%!error rsdec (1, 2) +%!error rsdec (1, 2, 3, 4, 5, 6) +*/ + +// PKG_ADD: autoload ("bchenco", "gf.oct"); +// PKG_DEL: autoload ("bchenco", "gf.oct", "remove"); +DEFUN_DLD (bchenco, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{code} =} bchenco (@var{msg}, @var{n}, @var{k})\n\ +@deftypefnx {Loadable Function} {@var{code} =} bchenco (@var{msg}, @var{n}, @var{k}, @var{g})\n\ +@deftypefnx {Loadable Function} {@var{code} =} bchenco (@dots{}, @var{parpos})\n\ +Encodes the message @var{msg} using a [@var{n},@var{k}] BCH coding.\n\ +The variable @var{msg} is a binary array with @var{k} columns and an\n\ +arbitrary number of rows. Each row of @var{msg} represents a single symbol\n\ +to be coded by the BCH coder. The coded message is returned in the binary\n\ +array @var{code} containing @var{n} columns and the same number of rows as\n\ +@var{msg}.\n\ +\n\ +The use of @code{bchenco} can be seen in the following short example.\n\ +\n\ +@example\n\ +m = 3; n = 2^m -1; k = 4;\n\ +msg = randint (10,k);\n\ +code = bchenco (msg, n, k);\n\ +@end example\n\ +\n\ +Valid codes can be found using @code{bchpoly}. In general the codeword\n\ +length @var{n} should be of the form @code{2^@var{m}-1}, where m is an\n\ +integer. However, shortened BCH codes can be used such that if\n\ +@code{[2^@var{m}-1,@var{k}]} is a valid code\n\ +@code{[2^@var{m}-1-@var{x},@var{k}-@var{x}]}\n is also a valid code using\n\ +the same generator polynomial.\n\ +\n\ +By default the generator polynomial used in the BCH coding is\n\ +based on the properties of the Galois Field GF(2^@var{m}). This\n\ +default generator polynomial can be overridden by a polynomial in @var{g}.\n\ +Suitable generator polynomials can be constructed with @code{bchpoly}.\n\ +\n\ +By default the parity symbols are placed at the beginning of the coded\n\ +message. The variable @var{parpos} controls this positioning and can take\n\ +the values @code{\"beginning\"} or @code{\"end\"}.\n\ +@seealso{bchpoly, bchdeco, encode}\n\ +@end deftypefn") +{ + octave_value retval; + int nargin = args.length (); + + if (nargin < 3 || nargin > 5) + { + print_usage (); + return retval; + } + + Matrix msg = args(0).matrix_value (); + int nsym = msg.rows (); + int nn = args(1).nint_value (); + int k = args(2).nint_value (); + + int m = 1; + while (nn > (1<<m)) + m++; + + int n = (1<<m) - 1; + + if (msg.cols () != k) + { + error ("bchenco: message contains incorrect number of symbols"); + return retval; + } + + if ((n < 3) || (nn < k) || (m > __OCTAVE_GALOIS_MAX_M)) + { + error ("bchenco: invalid values of message or codeword length"); + return retval; + } + + galois genpoly; + bool have_genpoly = false; + bool parity_at_end = false; + + for (int i = 3; i < nargin; i++) + { + if (args(i).is_string ()) + { + std::string parstr = args(i).string_value (); + for (int j = 0; j < (int)parstr.length (); j++) + parstr[j] = toupper (parstr[j]); + + if (!parstr.compare("END")) + { + parity_at_end = true; + } + else if (!parstr.compare("BEGINNING")) + { + parity_at_end = false; + } + else + { + error ("bchenco: unrecoginized parity position"); + return retval; + } + } + else + { + have_genpoly = true; + genpoly = galois (args(i).matrix_value (), m); + if (genpoly.cols () > genpoly.rows ()) + genpoly = genpoly.transpose (); + + if (genpoly.cols () != 1) + { + error ("bchenco: the generator polynomial must be a vector"); + return retval; + } + + if (genpoly.rows () != nn-k+1) + { + error ("bchenco: generator polynomial has incorrect order"); + return retval; + } + } + } + + if (!have_genpoly) + { + // The code below is basically bchpoly.m in C++, so if there is a need + // it can be used to rewrite bchpoly as an oct-file... + + RowVector found (n, 0); + found(0) = 1; + galois c (1, m, 0, m); + c(0, 0) = c.index_of (1); + Array<int> cs (dim_vector (1, 1), 1); + + int nc = 1; + + // Find the cyclotomic cosets of GF(2^m) + while (found.min () == 0) + { + int idx = n; + for (int i = 0; i<n; i++) + if ((found(i) == 0) && (c.index_of (i+1) < idx)) + idx = c.index_of (i+1); + + c.resize (dim_vector (nc+1, m)); + cs.resize (dim_vector (nc+1, 1)); + c(nc, 0) = idx; + found(c.alpha_to (idx)-1) = 1; + cs(nc) = 1; + int r = idx; + while ((r = modn (r<<1, m, n)) > idx) + { + c(nc, cs(nc)) = r; + found(c.alpha_to (r)-1) = 1; + cs(nc) += 1; + } + nc++; + } + + // Re-use the found vector with 1==not-found !!! + found.resize (nc); + + galois f (1, 0, 0, m); + int t = 0; + int nf = 0; + do + { + t++; + for (int i = 0; i < nc; i++) + { + if (found(i) == 1) + { + for (int j = 2*(t-1); j<2*t; j++) + { + int flag = 0; + for (int l = 0; l < cs(i); l++) + { + if (c(i, l) == j+1) + { + f.resize (dim_vector (1, nf+cs(i))); + for (int ll = 0; ll < cs(i); ll++) + f(0, nf+ll) = c(i, ll); + found(i) = 0; + nf += cs(i); + flag = 1; + break; + } + } + if (flag) break; + } + } + } + } + while (nf < nn - k); + + if (nf != nn - k) + { + error ("bchenco: can not find valid generator polynomial for parameters"); + return retval; + } + + // Create polynomial of right length. + genpoly = galois (nf+1, 1, 0, m); + + genpoly(0, 0) = 1; + for (int i = 0; i < nf; i++) + { + genpoly(i+1, 0) = 1; + + // Multiply genpoly by @**(root + x) + for (int l = i; l > 0; l--) + { + if (genpoly(l, 0) != 0) + genpoly(l, 0) = genpoly(l-1, 0) + ^ genpoly.alpha_to (modn (genpoly.index_of (genpoly(l, 0)) + f(0, i), + m, n)); + else + genpoly(l, 0) = genpoly(l-1, 0); + } + // genpoly(0,0) can never be zero + genpoly(0, 0) = genpoly.alpha_to (modn (genpoly.index_of (genpoly(0, 0)) + + f(0, i), + m, n)); + } + } + + // Add space for parity block + msg.resize (nsym, nn, 0); + + // The code below basically finds the parity bits by treating the + // message as a polynomial and dividing it by the generator polynomial. + // The parity bits are then the remainder of this division. + // + // This code could just as easily be written as + // [ignore par] = gdeconv(gf(msg), gf(genpoly)); + // But the code below has the advantage of being 20 times faster :-) + + if (parity_at_end) + { + for (int l = 0; l < nsym; l++) + { + for (int i = 0; i < k; i++) + { + int feedback = (int)msg(l, i) ^ (int)msg(l, k); + if (feedback != 0) + { + for (int j = 0; j < nn-k-1; j++) + if (genpoly(nn-k-j-1, 0) != 0) + msg(l, k+j) = (int)msg(l, k+j+1) ^ feedback; + else + msg(l, k+j) = msg(l, k+j+1); + msg(l, nn-1) = genpoly(0, 0) & feedback; + } + else + { + for (int j = k; j < nn-1; j++) + msg(l, j) = msg(l, j+1); + msg(l, nn-1) = 0; + } + } + } + } + else + { + for (int l = 0; l < nsym; l++) + { + for (int i=k; i > 0; i--) + msg(l, i+nn-k-1) = msg(l, i-1); + for (int i = 0; i<nn-k; i++) + msg(l, i) = 0; + } + + for (int l = 0; l < nsym; l++) + { + for (int i = k-1; i >= 0; i--) + { + int feedback = (int)msg(l, nn-k+i) ^ (int)msg(l, nn-k-1); + if (feedback != 0) + { + for (int j = nn - k -1; j > 0; j--) + if (genpoly(j, 0) != 0) + msg(l, j) = (int)msg(l, j-1) ^ feedback; + else + msg(l, j) = msg(l, j-1); + msg(l, 0) = genpoly(0, 0) & feedback; + } + else + { + for (int j = nn - k - 1; j > 0; j--) + msg(l, j) = msg(l, j-1); + msg(l, 0) = 0; + } + } + } + } + + retval = msg; + return retval; +} + +/* +%% Test input validation +%!error bchenco () +%!error bchenco (1) +%!error bchenco (1, 2) +%!error bchenco (1, 2, 3, 4, 5, 6) +*/ + +// PKG_ADD: autoload ("bchdeco", "gf.oct"); +// PKG_DEL: autoload ("bchdeco", "gf.oct", "remove"); +DEFUN_DLD (bchdeco, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{msg} =} bchdeco (@var{code}, @var{k}, @var{t})\n\ +@deftypefnx {Loadable Function} {@var{msg} =} bchdeco (@var{code}, @var{k}, @var{t}, @var{prim})\n\ +@deftypefnx {Loadable Function} {@var{msg} =} bchdeco (@dots{}, @var{parpos})\n\ +@deftypefnx {Loadable Function} {[@var{msg}, @var{err}] =} bchdeco (@dots{})\n\ +@deftypefnx {Loadable Function} {[@var{msg}, @var{err}, @var{ccode}] =} bchdeco (@dots{})\n\ +Decodes the coded message @var{code} using a BCH coder. The message length\n\ +of the coder is defined in variable @var{k}, and the error correction\n\ +capability of the code is defined in @var{t}.\n\ +\n\ +The variable @var{code} is a binary array with @var{n} columns and an\n\ +arbitrary number of rows. Each row of @var{code} represents a single symbol\n\ +to be decoded by the BCH coder. The decoded message is returned in the\n\ +binary array @var{msg} containing @var{k} columns and the same number of\n\ +rows as @var{code}.\n\ +\n\ +The use of @code{bchdeco} can be seen in the following short example.\n\ +\n\ +@example\n\ +m = 3; n = 2^m -1; k = 4; t = 1;\n\ +msg = randint (10, k);\n\ +code = bchenco (msg, n, k);\n\ +noisy = mod (randerr (10,n) + code, 2);\n\ +[dec, err] = bchdeco (msg, k, t);\n\ +@end example\n\ +\n\ +Valid codes can be found using @code{bchpoly}. In general the codeword\n\ +length @var{n} should be of the form @code{2^@var{m}-1}, where m is an\n\ +integer. However, shortened BCH codes can be used such that if\n\ +@code{[2^@var{m}-1,@var{k}]} is a valid code\n\ +@code{[2^@var{m}-1-@var{x},@var{k}-@var{x}]}\n is also a valid code using\n\ +the same generator polynomial.\n\ +\n\ +By default the BCH coding is based on the properties of the Galois\n\ +Field GF(2^@var{m}). The primitive polynomial used in the Galois\n\ +can be overridden by a primitive polynomial in @var{prim}. Suitable\n\ +primitive polynomials can be constructed with @code{primpoly}. The form\n\ +of @var{prim} maybe be either a integer representation of the primitive\n\ +polynomial as given by @code{primpoly}, or a binary representation that\n\ +might be constructed like\n\ +\n\ +@example\n\ +m = 3;\n\ +prim = de2bi (primpoly (m));\n\ +@end example\n\ +\n\ +By default the parity symbols are assumed to be placed at the beginning of\n\ +the coded message. The variable @var{parpos} controls this positioning and\n\ +can take the values @code{\"beginning\"} or @code{\"end\"}.\n\ +@seealso{bchpoly, bchenco, decode, primpoly}\n\ +@end deftypefn") +{ + octave_value_list retval; + int nargin = args.length (); + + if (nargin < 3 || nargin > 5) + { + print_usage (); + return retval; + } + + Matrix code = args(0).matrix_value (); + int nsym = code.rows (); + int nn = code.cols (); + int k = args(1).nint_value (); + int t = args(2).nint_value (); + int t2 = t << 1; + + int m = 1; + while (nn > (1<<m)) + m++; + + int n = (1<<m) - 1; + + if ((n < 3) || (n < k) || (m > __OCTAVE_GALOIS_MAX_M)) + { + error ("bchdeco: invalid values of message or codeword length"); + return retval; + } + + int prim = 0; // primitve polynomial of zero flags default + bool parity_at_end = false; + + for (int i = 3; i < nargin; i++) + { + if (args(i).is_string ()) + { + std::string parstr = args(i).string_value (); + for (int j = 0; j < (int)parstr.length (); j++) + parstr[j] = toupper (parstr[j]); + + if (!parstr.compare("END")) + { + parity_at_end = true; + } + else if (!parstr.compare("BEGINNING")) + { + parity_at_end = false; + } + else + { + error ("bchdeco: unrecoginized parity position"); + return retval; + } + } + else + { + if (args(i).is_real_scalar ()) + prim = args(i).int_value (); + else + { + Matrix tmp = args(i).matrix_value (); + + if (tmp.cols () > tmp.rows ()) + tmp = tmp.transpose (); + + if (tmp.cols () != 1) + { + error ("bchdeco: the primitve polynomial must be a scalar " + "or a vector"); + return retval; + } + + prim = 0; + for (int i = 0; i < tmp.rows (); i++) + if ((int)tmp(i, 0) & 1) + prim |= (1<<i); + } + } + } + + // Create a variable in the require Galois Field to have access to the + // lookup tables alpha_to and index_of. + galois tables (1, 1, 0, m, prim); + ColumnVector nerr (nsym, 0); + + for (int lsym = 0; lsym < nsym; lsym++) + { + /* first form the syndromes */ + Array<int> s (dim_vector(t2+1, 1), 0); + bool syn_error = false; + + for (int i = 1; i <= t2; i++) + { + for (int j = 0; j < nn; j++) + { + if (parity_at_end) + { + if (code(lsym, nn-j-1) != 0) + s(i) ^= tables.alpha_to (modn (i*j, m, n)); + } + else + { + if (code(lsym, j) != 0) + s(i) ^= tables.alpha_to (modn (i*j, m, n)); + } + } + if (s(i) != 0) + syn_error = true; /* set error flag if non-zero syndrome */ + + } + + if (syn_error) + { /* if there are errors, try to correct them */ + int q, u; + Array<int> d (dim_vector (t2+2, 1)), l(dim_vector (t2+2, 1)), + u_lu(dim_vector (t2+2, 1)), reg(dim_vector (t2+2, 1)), + elp(dim_vector (t2+2, t2+2)); + + /* convert syndrome from polynomial form to index form */ + for (int i = 1; i <= t2; i++) + s(i) = tables.index_of (s(i)); + + /* + * Compute the error location polynomial via the Berlekamp + * iterative algorithm. Following the terminology of Lin and + * Costello's book : d(u) is the 'mu'th discrepancy, where + * u='mu'+1 and 'mu' (the Greek letter!) is the step number + * ranging from -1 to 2*t (see L&C), l(u) is the degree of + * the elp at that step, and u_l(u) is the difference between + * the step number and the degree of the elp. + */ + /* initialise table entries */ + d(0) = 0; /* index form */ + d(1) = s(1); /* index form */ + elp(0, 0) = 0; /* index form */ + elp(1, 0) = 1; /* polynomial form */ + for (int i = 1; i < t2; i++) + { + elp(0, i) = n; /* index form */ + elp(1, i) = 0; /* polynomial form */ + } + l(0) = 0; + l(1) = 0; + u_lu(0) = -1; + u_lu(1) = 0; + u = 0; + + do + { + u++; + if (d(u) == n) + { + l(u + 1) = l(u); + for (int i = 0; i <= l(u); i++) + { + elp(u + 1, i) = elp(u, i); + elp(u, i) = tables.index_of (elp(u, i)); + } + } + else + /* + * search for words with greatest u_lu(q) for + * which d(q)!=0 + */ + { + q = u - 1; + while ((d(q) == n) && (q > 0)) + q--; + /* have found first non-zero d(q) */ + if (q > 0) + { + int j = q; + do + { + j--; + if ((d(j) != n) && (u_lu(q) < u_lu(j))) + q = j; + } + while (j > 0); + } + + /* + * have now found q such that d(u)!=0 and + * u_lu(q) is maximum + */ + /* store degree of new elp polynomial */ + if (l(u) > l(q) + u - q) + l(u + 1) = l(u); + else + l(u + 1) = l(q) + u - q; + + /* form new elp(x) */ + for (int i = 0; i < t2; i++) + elp(u + 1, i) = 0; + for (int i = 0; i <= l(q); i++) + if (elp(q, i) != n) + elp(u + 1, i + u - q) = + tables.alpha_to (modn ((d(u) + n - d(q) + elp(q, i)), m, n)); + for (int i = 0; i <= l(u); i++) + { + elp(u + 1, i) ^= elp(u, i); + elp(u, i) = tables.index_of (elp(u, i)); + } + } + u_lu(u + 1) = u - l(u + 1); + + /* form (u+1)th discrepancy */ + if (u < t2) + { + /* no discrepancy computed on last iteration */ + d(u + 1) = tables.alpha_to (s(u + 1)); + + for (int i = 1; i <= l(u + 1); i++) + if ((s(u + 1 - i) != n) && (elp(u + 1, i) != 0)) + d(u + 1) ^= tables.alpha_to (modn (s(u + 1 - i) + + tables.index_of (elp(u + 1, i)), + m, n)); + /* put d(u+1) into index form */ + d(u + 1) = tables.index_of (d(u + 1)); + } + } + while ((u < t2) && (l(u + 1) <= t)); + + u++; + if (l(u) <= t) + {/* Can correct errors */ + int count; + Array<int> loc (dim_vector (t+2, 1)); + + /* put elp into index form */ + for (int i = 0; i <= l(u); i++) + elp(u, i) = tables.index_of (elp(u, i)); + + /* Chien search: find roots of the error location polynomial */ + for (int i = 1; i <= l(u); i++) + reg(i) = elp(u, i); + count = 0; + for (int i = 1; i <= n; i++) + { + q = 1; + for (int j = 1; j <= l(u); j++) + if (reg(j) != n) + { + reg(j) = modn ((reg(j) + j), m, n); + q ^= tables.alpha_to (reg(j)); + } + if (!q) + { /* store root and error + * location number indices */ + loc(count) = n - i; + count++; + if (count > l(u)) + break; + } + } + + if (count == l(u)) + { + /* no. roots = degree of elp hence <= t errors */ + nerr(lsym) = l(u); + for (int i = 0; i < l(u); i++) + if (parity_at_end) + code(lsym, nn-loc(i)-1) = + (int)code(lsym, nn-loc(i)-1) ^ 1; + else + code(lsym, loc(i)) = (int)code(lsym, loc(i)) ^ 1; + } + else /* elp has degree >t hence cannot solve */ + nerr(lsym) = -1; + } + else + nerr(lsym) = -1; + } + } + + Matrix msg (nsym, k); + if (parity_at_end) + { + for (int l = 0; l < nsym; l++) + for (int i = 0; i < k; i++) + msg(l, i) = code(l, i); + } + else + { + for (int l = 0; l < nsym; l++) + for (int i = 0; i < k; i++) + msg(l, i) = code(l, nn-k+i); + } + + retval(0) = octave_value (msg); + retval(1) = octave_value (nerr); + retval(2) = octave_value (code); + return retval; +} + +/* +%% Test input validation +%!error bchdeco () +%!error bchdeco (1) +%!error bchdeco (1, 2) +%!error bchdeco (1, 2, 3, 4, 5, 6) +*/ diff -uNr a/src/__gfweight__.cc b/src/__gfweight__.cc --- a/src/__gfweight__.cc 2015-04-04 12:28:43.938510295 -0400 +++ b/src/__gfweight__.cc 2018-04-09 13:51:01.973935827 -0400 @@ -68,7 +68,7 @@ if (k > 128) { octave_stdout << "__gfweight__: this is likely to take a very long time!!\n"; - flush_octave_stdout (); + octave::flush_stdout (); } Array<char> codeword (dim_vector (n, 1), 0); diff -uNr a/src/__gfweight__.cc~ b/src/__gfweight__.cc~ --- a/src/__gfweight__.cc~ 1969-12-31 19:00:00.000000000 -0500 +++ b/src/__gfweight__.cc~ 2015-04-04 12:28:43.938510295 -0400 @@ -0,0 +1,89 @@ +//Copyright (C) 2003 David Bateman +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, see +// <http://www.gnu.org/licenses/>. +// +// In addition to the terms of the GPL, you are permitted to link this +// program with any Open Source program, as defined by the Open Source +// Initiative (www.opensource.org) + +#include <octave/oct.h> + +static int +get_weight (const Array<char>& codeword, const Matrix& gen, + int weight, int depth, int start, int n, int k) +{ + int retval = weight; + + for (int i = start; i < k ; i++) + { + OCTAVE_QUIT; + + Array<char> new_codeword (codeword); + int tmp = 0; + for (int j = 0; j < n; j++) + if (new_codeword (j) ^= (char)gen(i,j)) + tmp++; + if (tmp < retval) + retval = tmp; + if (depth < retval) + retval = get_weight (new_codeword, gen, retval, depth+1, i+1, n, k); + } + return retval; +} + +DEFUN_DLD (__gfweight__, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{w} =} __gfweight__ (@var{gen})\n\ +Returns the minimum distance @var{w} of the generator matrix @var{gen}.\n\ +The codeword length is @var{k}.\n\ +\n\ +This is an internal function of @code{gfweight}. You should use\n\ +@code{gfweight} rather than use this function directly.\n\ +@seealso{gfweight}\n\ +@end deftypefn") +{ + + if (args.length () != 1) + { + print_usage (); + return octave_value (); + } + + Matrix gen = args(0).matrix_value (); + int k = gen.rows (); + int n = gen.columns (); + + if (k > 128) + { + octave_stdout << "__gfweight__: this is likely to take a very long time!!\n"; + flush_octave_stdout (); + } + + Array<char> codeword (dim_vector (n, 1), 0); + return octave_value ((double)get_weight (codeword, gen, n - k + 1, 1, + 0, n, k)); +} + +/* +%% Test input validation +%!error __gfweight__ () +%!error __gfweight__ (1, 2) +*/ + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ diff -uNr a/src/ov-galois.cc b/src/ov-galois.cc --- a/src/ov-galois.cc 2018-04-09 13:25:42.884981069 -0400 +++ b/src/ov-galois.cc 2018-04-09 13:27:40.103502409 -0400 @@ -21,7 +21,7 @@ #include <iostream> #include <octave/byte-swap.h> -#include <octave/gripes.h> +#include <octave/errwarn.h> #include <octave/lo-ieee.h> #include <octave/oct-locbuf.h> #include <octave/oct-obj.h> @@ -328,13 +328,13 @@ if (rows () > 0 && columns () > 0) { - gripe_implicit_conversion ("Octave:array-as-scalar", + warn_implicit_conversion ("Octave:array-as-scalar", "real matrix", "real scalar"); retval = (double) gval (0, 0); } else - gripe_invalid_conversion ("galois", "real scalar"); + err_invalid_conversion ("galois", "real scalar"); return retval; } @@ -348,13 +348,13 @@ if (rows () > 0 && columns () > 0) { - gripe_implicit_conversion ("Octave:array-as-scalar", + warn_implicit_conversion ("Octave:array-as-scalar", "real matrix", "real scalar"); retval = (double) gval (0, 0); } else - gripe_invalid_conversion ("galois", "complex scalar"); + err_invalid_conversion ("galois", "complex scalar"); return retval; } diff -uNr a/src/ov-galois.h b/src/ov-galois.h --- a/src/ov-galois.h 2018-04-09 13:25:42.872981630 -0400 +++ b/src/ov-galois.h 2018-04-09 14:09:33.961913785 -0400 @@ -49,7 +49,6 @@ #endif class octave_value_list; -class tree_walker; // Data structures. @@ -100,7 +99,7 @@ bool is_defined (void) const { return true; } - bool is_numeric_type (void) const { return true; } + bool isnumeric (void) const { return true; } bool is_constant (void) const { return true; } @@ -124,7 +123,7 @@ bool is_real_matrix (void) const { return false; } - bool is_real_type (void) const { return false; } + bool isreal (void) const { return false; } // FIXME bool valid_as_scalar_index (void) const { return false; } diff -uNr a/src/ov-galois.h~ b/src/ov-galois.h~ --- a/src/ov-galois.h~ 1969-12-31 19:00:00.000000000 -0500 +++ b/src/ov-galois.h~ 2018-04-09 13:35:21.137945379 -0400 @@ -0,0 +1,183 @@ +//Copyright (C) 2003 David Bateman +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, see +// <http://www.gnu.org/licenses/>. +// +// In addition to the terms of the GPL, you are permitted to link this +// program with any Open Source program, as defined by the Open Source +// Initiative (www.opensource.org) + +#if !defined (octave_galois_h) +#define octave_galois_h 1 + +#include <octave/ov.h> +#include <octave/ov-typeinfo.h> + +#include "galois.h" + +// The keys of the values in the octave map +#define __GALOIS_PRIMPOLY_STR "prim_poly" +#define __GALOIS_ORDER_STR "m" +#define __GALOIS_DATA_STR "x" +#ifdef GALOIS_DISP_PRIVATES +#define __GALOIS_LENGTH_STR "n" +#define __GALOIS_ALPHA_TO_STR "alpha_to" +#define __GALOIS_INDEX_OF_STR "index_of" +#endif + +#if !defined (HAVE_OCTAVE_HDF5_ID_TYPE) +#if defined (HAVE_HDF5) +typedef hid_t octave_hdf5_id; +#else +typedef int octave_hdf5_id; +#endif +#endif + +#if ! defined (OV_REP_TYPE) +# define OV_REP_TYPE octave_base_value +#endif + +class octave_value_list; +class tree_walker; + +// Data structures. + +class +octave_galois : public octave_base_value +{ +public: + + octave_galois (const Matrix& data = Matrix (0, 0), const int _m = 1, + const int _primpoly = 0) + { gval = galois (data, _m, _primpoly); } + + octave_galois (const galois& gm) + : octave_base_value (), gval (gm) { } + octave_galois (const octave_galois& s) + : octave_base_value (), gval (s.gval) { } + + ~octave_galois (void) { }; + + OV_REP_TYPE *clone (void) const { return new octave_galois (*this); } + OV_REP_TYPE *empty_clone (void) const { return new octave_galois (); } + + octave_value subsref (const std::string &type, + const std::list<octave_value_list>& idx); + + octave_value_list subsref (const std::string& type, + const std::list<octave_value_list>& idx, int) + { return subsref (type, idx); } + + octave_value do_index_op (const octave_value_list& idx, + bool resize_ok); + + octave_value do_index_op (const octave_value_list& idx) + { return do_index_op (idx, 0); } + + void assign (const octave_value_list& idx, const galois& rhs); + + dim_vector dims (void) const { return gval.dims (); } + + octave_value resize (const dim_vector& dv, bool) const; + + size_t byte_size (void) const { return gval.byte_size (); } + + octave_value all (int dim = 0) const { return gval.all (dim); } + octave_value any (int dim = 0) const { return gval.any (dim); } + + bool is_matrix_type (void) const { return true; } + + bool is_defined (void) const { return true; } + + bool isnumeric (void) const { return true; } + + bool is_constant (void) const { return true; } + + bool is_true (void) const; + + bool is_galois_type (void) const { return true; } + + bool print_as_scalar (void) const; + +#if defined (HAVE_OCTAVE_BASE_VALUE_PRINT_CONST) + void print (std::ostream& os, bool pr_as_read_syntax = false) const; +#else + void print (std::ostream& os, bool pr_as_read_syntax = false); +#endif + + void print_raw (std::ostream& os, bool pr_as_read_syntax = false) const; + + bool print_name_tag (std::ostream& os, const std::string& name) const; + + void print_info (std::ostream& os, const std::string& prefix) const; + + bool is_real_matrix (void) const { return false; } + + bool isreal (void) const { return false; } + + // FIXME + bool valid_as_scalar_index (void) const { return false; } + + double double_value (bool = false) const; + + double scalar_value (bool frc_str_conv = false) const + { return double_value (frc_str_conv); } + + Matrix matrix_value (bool = false) const; + + NDArray array_value (bool = false) const; + + Complex complex_value (bool = false) const; + + ComplexMatrix complex_matrix_value (bool = false) const + { return ComplexMatrix ( matrix_value ()); } + + galois galois_value (void) const { return gval; } + + octave_value_list dotref (const octave_value_list& idx); + + int m (void) const { return gval.m (); } + int primpoly (void) const { return gval.primpoly (); } + + bool save_ascii (std::ostream& os); + + bool load_ascii (std::istream& is); + + bool save_binary (std::ostream& os, bool& save_as_floats); + + bool load_binary (std::istream& is, bool swap, + oct_mach_info::float_format fmt); + + bool save_hdf5 (octave_hdf5_id loc_id, const char *name, bool save_as_floats); + + bool load_hdf5 (octave_hdf5_id loc_id, const char *name); + +private: + // The array used to managed the Galios Field data + galois gval; + +#if defined (DECLARE_OCTAVE_ALLOCATOR) + DECLARE_OCTAVE_ALLOCATOR +#endif + + DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/