Mercurial > mxe-octave
changeset 4653:00e61c4a5657
fixes for package build errors due to API changes
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 09 Apr 2018 12:12:44 -0400 |
| parents | 2f14bc0c6d0c |
| children | b9e4ebcad82f |
| files | dist-files.mk src/of-communications-5-fixes.patch src/of-database-3-fixes.patch src/of-dicom-1-fixes.patch src/of-fits-2-fixes.patch src/of-general-1-fixes.patch src/of-linear-algebra-4-fixes.patch src/of-optim-1-fixes.patch src/of-signal-1-fixes.patch src/of-strings-1-fixes.patch src/of-struct-1-fixes.patch src/of-video-1-fixes.patch |
| diffstat | 12 files changed, 6762 insertions(+), 6 deletions(-) [+] |
line wrap: on
line diff
--- a/dist-files.mk Mon Apr 09 07:11:54 2018 -0400 +++ b/dist-files.mk Mon Apr 09 12:12:44 2018 -0400 @@ -455,22 +455,26 @@ of-communications-2-fixes.patch \ of-communications-3-fixes.patch \ of-communications-4-fixes.patch \ + of-communications-5-fixes.patch \ of-communications.mk \ of-control.mk \ of-data-smoothing.mk \ of-database-1-cross-fixes.patch \ of-database-2-dev-fixes.patch \ + of-database-3-fixes.patch \ of-database.mk \ of-dataframe.mk \ of-dicom-1-fixes.patch \ of-dicom.mk \ of-financial.mk \ of-fits-1-cross-fixes.patch \ + of-fits-2-fixes.patch \ of-fits.mk \ of-fl-core-1-fixes.patch \ of-fl-core.mk \ of-fuzzy-logic-toolkit.mk \ of-ga.mk \ + of-general-1-fixes.patch \ of-general-1-symtab-fixes.patch \ of-general.mk \ of-generate_html.mk \ @@ -484,6 +488,7 @@ of-linear-algebra-1-cross-fixes.patch \ of-linear-algebra-2-dev-fixes.patch \ of-linear-algebra-3-fixes.patch \ + of-linear-algebra-4-fixes.patch \ of-linear-algebra.mk \ of-lssa.mk \ of-ltfat-1-fixes.patch \ @@ -502,6 +507,7 @@ of-odepkg-1-fixes.patch \ of-odepkg-2-fixes.patch \ of-odepkg.mk \ + of-optim-1-fixes.patch \ of-optim.mk \ of-quaternion-1-cross-fixes.patch \ of-quaternion.mk \ @@ -518,6 +524,7 @@ of-stk.mk \ of-strings-1-fixes.patch \ of-strings.mk \ + of-struct-1-fixes.patch \ of-struct.mk \ of-tisean.mk \ of-tsa.mk \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/of-communications-5-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -0,0 +1,6045 @@ +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: *** ++*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/of-database-3-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -0,0 +1,117 @@ +diff -uNr a/src/error-helpers.cc b/src/error-helpers.cc +--- a/src/error-helpers.cc 2016-09-18 14:36:40.000000000 -0400 ++++ b/src/error-helpers.cc 2018-04-09 14:52:43.380796194 -0400 +@@ -24,7 +24,7 @@ + // call verror + #ifdef HAVE_OCTAVE_VERROR_ARG_EXC + void +-c_verror (octave_execution_exception& e, const char *fmt, ...) ++c_verror (octave::execution_exception& e, const char *fmt, ...) + { + va_list args; + va_start (args, fmt); +@@ -33,7 +33,7 @@ + } + #else + void +-c_verror (const octave_execution_exception&, const char *fmt, ...) ++c_verror (const octave::execution_exception&, const char *fmt, ...) + { + va_list args; + va_start (args, fmt); +diff -uNr a/src/error-helpers.h b/src/error-helpers.h +--- a/src/error-helpers.h 2016-09-18 14:36:40.000000000 -0400 ++++ b/src/error-helpers.h 2018-04-09 14:52:56.204181295 -0400 +@@ -21,9 +21,9 @@ + + // call verror + #ifdef HAVE_OCTAVE_VERROR_ARG_EXC +-void c_verror (octave_execution_exception&, const char *, ...); ++void c_verror (octave::execution_exception&, const char *, ...); + #else +-void c_verror (const octave_execution_exception&, const char *, ...); ++void c_verror (const octave::execution_exception&, const char *, ...); + #endif + + // call verror +@@ -33,7 +33,7 @@ + // both if Octave uses exceptions for errors and if it still uses + // error_state. In the latter case return 'retval'. + #ifdef HAVE_OCTAVE_ERROR_STATE +- // can throw octave_execution_exception despite of this ++ // can throw octave::execution_exception despite of this + #define CHECK_ERROR(code, retval, ...) \ + try \ + { \ +@@ -46,7 +46,7 @@ + return retval; \ + } \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -58,7 +58,7 @@ + { \ + code ; \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -70,7 +70,7 @@ + // Octave doesn't throw exceptions for errors but still uses + // error_state. + #ifdef HAVE_OCTAVE_ERROR_STATE +- // can throw octave_execution_exception despite of this ++ // can throw octave::execution_exception despite of this + #define CHECK_ERROR_EXIT1(code, ...) \ + try \ + { \ +@@ -83,7 +83,7 @@ + exit (1); \ + } \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -95,7 +95,7 @@ + { \ + code ; \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -107,7 +107,7 @@ + // Octave uses exceptions for errors and if it still uses + // error_state. In the latter case reset error_state to 0. + #ifdef HAVE_OCTAVE_ERROR_STATE +- // can throw octave_execution_exception despite of this ++ // can throw octave::execution_exception despite of this + #define SET_ERR(code, err) \ + err = false; \ + \ +@@ -120,7 +120,7 @@ + err = true; \ + } \ + } \ +- catch (octave_execution_exception&) \ ++ catch (octave::execution_exception&) \ + { \ + err = true; \ + } +@@ -130,7 +130,7 @@ + { \ + code ; \ + } \ +- catch (octave_execution_exception&) \ ++ catch (octave::execution_exception&) \ + { \ + err = true; \ + }
--- a/src/of-dicom-1-fixes.patch Mon Apr 09 07:11:54 2018 -0400 +++ b/src/of-dicom-1-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -117,3 +117,36 @@ std::ifstream fin(resolved_filename.c_str()); if (!fin) { +diff -uNr a/src/dicomwrite.cpp b/src/dicomwrite.cpp +--- a/src/dicomwrite.cpp 2017-02-23 01:28:24.157669038 -0500 ++++ b/src/dicomwrite.cpp 2018-04-09 14:40:14.224292229 -0400 +@@ -139,7 +139,7 @@ + } + + void struct2metadata(gdcm::ImageWriter *w, gdcm::File *file, const octave_value & ov, bool trial, int sequenceDepth) { +- if(!ov.is_map()){ ++ if(!ov.isstruct()){ + error(QUOTED(OCT_FN_NAME)": 3rd arg should be struct holding metadata. it is %s",ov.type_name().c_str()); + throw std::exception() ; + } +@@ -297,7 +297,7 @@ + if (trial) octave_stdout << '[' << buf << ']' << std::endl; + de->SetByteValue( buf, gdcm::VL((uint32_t)strlen(buf)) ); + } else if ( entry.GetVR() & gdcm::VR::SQ) { // sequence +- if (!ov->is_map()) { ++ if (!ov->isstruct()) { + warning(QUOTED(OCT_FN_NAME)": dicomdict gives VR of SQ for %s, octave value is %s", keyword.c_str(), ov->class_name().c_str()); + } + octave_stdout << std::endl; +diff -uNr a/src/isdicom.cpp b/src/isdicom.cpp +--- a/src/isdicom.cpp 2017-02-23 01:28:24.157669038 -0500 ++++ b/src/isdicom.cpp 2018-04-09 14:39:33.322200138 -0400 +@@ -32,7 +32,7 @@ + reader.SetFileName (filename.c_str ()); + // gdcm::Reader.Read() will return false if the file does not exists but + // also prints to stderr so we check it first. +- return file_stat (filename).exists () && reader.Read (); ++ return octave::sys::file_stat (filename).exists () && reader.Read (); + } + + DEFUN_DLD (isdicom, args, ,
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/of-fits-2-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -0,0 +1,36 @@ +diff -uNr a/src/read_fits_image.cc b/src/read_fits_image.cc +--- a/src/read_fits_image.cc 2015-06-11 07:19:38.000000000 -0400 ++++ b/src/read_fits_image.cc 2018-04-09 14:31:49.259838178 -0400 +@@ -194,7 +194,7 @@ + return true; + } + double val = args(1).double_value(); +- if( (D_NINT( val ) != val) || (val < 0) ) ++ if( (octave::math::x_nint( val ) != val) || (val < 0) ) + { + error( "read_fits_image: second argument must be a non-negative scalar integer value" ); + return true; +diff -uNr a/src/save_fits_image.cc b/src/save_fits_image.cc +--- a/src/save_fits_image.cc 2015-06-11 07:19:38.000000000 -0400 ++++ b/src/save_fits_image.cc 2018-04-09 14:31:37.948365332 -0400 +@@ -81,7 +81,7 @@ + else if( args(2).is_scalar_type() ) + { + double val = args(2).double_value(); +- if( (D_NINT( val ) == val) ) ++ if( (octave::math::x_nint( val ) == val) ) + { + if( BYTE_IMG == val ) + bitperpixel = BYTE_IMG; +diff -uNr a/src/save_fits_image_multi_ext.cc b/src/save_fits_image_multi_ext.cc +--- a/src/save_fits_image_multi_ext.cc 2015-06-11 07:19:38.000000000 -0400 ++++ b/src/save_fits_image_multi_ext.cc 2018-04-09 14:34:46.915556809 -0400 +@@ -78,7 +78,7 @@ + else if( args(2).is_scalar_type() ) + { + double val = args(2).double_value(); +- if( (D_NINT( val ) == val) ) ++ if( (octave::math::x_nint( val ) == val) ) + { + if( BYTE_IMG == val ) + bitperpixel = BYTE_IMG;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/of-general-1-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -0,0 +1,14 @@ +diff -uNr a/src/SHA1.cc b/src/SHA1.cc +--- a/src/SHA1.cc 2015-05-27 11:43:57.000000000 -0400 ++++ b/src/SHA1.cc 2018-04-09 12:09:29.175376040 -0400 +@@ -60,8 +60,8 @@ + hash_context c; + + if (nargin >2 || nargin ==0) { +- usage("SHA1"); +- return retval; ++ print_usage (); ++ return retval; + } + else if (nargin ==2 ){ + ColumnVector init( args(1).vector_value() );
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/of-linear-algebra-4-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -0,0 +1,43 @@ +diff -uNr a/src/gsvd.cc b/src/gsvd.cc +--- a/src/gsvd.cc 2015-01-24 14:54:11.000000000 -0500 ++++ b/src/gsvd.cc 2018-04-09 14:25:44.604821251 -0400 +@@ -23,8 +23,8 @@ + + #include "defun-dld.h" + #include "error.h" +-#include "gripes.h" +-#include "oct-obj.h" ++#include "errwarn.h" ++#include "ovl.h" + #include "pr-output.h" + #include "utils.h" + +@@ -211,7 +211,7 @@ + ? GSVD::sigma_only + : (nargout > 5) ? GSVD::std : GSVD::economy ); + +- if (argA.is_real_type () && argB.is_real_type ()) ++ if (argA.isreal () && argB.isreal ()) + { + Matrix tmpA = argA.matrix_value (); + Matrix tmpB = argB.matrix_value (); +@@ -253,7 +253,7 @@ + } + } + } +- else if (argA.is_complex_type () || argB.is_complex_type ()) ++ else if (argA.iscomplex () || argB.iscomplex ()) + { + ComplexMatrix ctmpA = argA.complex_matrix_value (); + ComplexMatrix ctmpB = argB.complex_matrix_value (); +@@ -296,8 +296,8 @@ + } + else + { +- gripe_wrong_type_arg ("gsvd", argA); +- gripe_wrong_type_arg ("gsvd", argB); ++ err_wrong_type_arg ("gsvd", argA); ++ err_wrong_type_arg ("gsvd", argB); + return retval; + } + }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/of-optim-1-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -0,0 +1,129 @@ +diff -uNr a/src/__bfgsmin.cc b/src/__bfgsmin.cc +--- a/src/__bfgsmin.cc 2016-09-18 13:31:55.000000000 -0400 ++++ b/src/__bfgsmin.cc 2018-04-09 12:24:30.969197504 -0400 +@@ -260,7 +260,7 @@ + found_improvement = __bisectionstep(step, obj, f, f_args, x, dx, minarg, verbose); + } + else found_improvement = 1; +- if (xisnan(obj)) { ++ if (octave::math::isnan(obj)) { + obj = obj_0; + if (verbose) warning("__stepsize: objective function crash in Newton step, falling back to bisection"); + found_improvement = __bisectionstep(step, obj, f, f_args, x, dx, minarg, verbose); +diff -uNr a/src/error-helpers.cc b/src/error-helpers.cc +--- a/src/error-helpers.cc 2016-09-18 13:31:55.000000000 -0400 ++++ b/src/error-helpers.cc 2018-04-09 12:17:15.913535622 -0400 +@@ -24,7 +24,7 @@ + // call verror + #ifdef HAVE_OCTAVE_VERROR_ARG_EXC + void +-c_verror (octave_execution_exception& e, const char *fmt, ...) ++c_verror (octave::execution_exception& e, const char *fmt, ...) + { + va_list args; + va_start (args, fmt); +@@ -33,7 +33,7 @@ + } + #else + void +-c_verror (const octave_execution_exception&, const char *fmt, ...) ++c_verror (const octave::execution_exception&, const char *fmt, ...) + { + va_list args; + va_start (args, fmt); +diff -uNr a/src/error-helpers.h b/src/error-helpers.h +--- a/src/error-helpers.h 2016-09-18 13:31:55.000000000 -0400 ++++ b/src/error-helpers.h 2018-04-09 12:17:51.079891166 -0400 +@@ -21,9 +21,9 @@ + + // call verror + #ifdef HAVE_OCTAVE_VERROR_ARG_EXC +-void c_verror (octave_execution_exception&, const char *, ...); ++void c_verror (octave::execution_exception&, const char *, ...); + #else +-void c_verror (const octave_execution_exception&, const char *, ...); ++void c_verror (const octave::execution_exception&, const char *, ...); + #endif + + // call verror +@@ -33,7 +33,7 @@ + // both if Octave uses exceptions for errors and if it still uses + // error_state. In the latter case return 'retval'. + #ifdef HAVE_OCTAVE_ERROR_STATE +- // can throw octave_execution_exception despite of this ++ // can throw octave::execution_exception despite of this + #define CHECK_ERROR(code, retval, ...) \ + try \ + { \ +@@ -46,7 +46,7 @@ + return retval; \ + } \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -58,7 +58,7 @@ + { \ + code ; \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -70,7 +70,7 @@ + // Octave doesn't throw exceptions for errors but still uses + // error_state. + #ifdef HAVE_OCTAVE_ERROR_STATE +- // can throw octave_execution_exception despite of this ++ // can throw octave::execution_exception despite of this + #define CHECK_ERROR_EXIT1(code, ...) \ + try \ + { \ +@@ -83,7 +83,7 @@ + exit (1); \ + } \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -95,7 +95,7 @@ + { \ + code ; \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -107,7 +107,7 @@ + // Octave uses exceptions for errors and if it still uses + // error_state. In the latter case reset error_state to 0. + #ifdef HAVE_OCTAVE_ERROR_STATE +- // can throw octave_execution_exception despite of this ++ // can throw octave::execution_exception despite of this + #define SET_ERR(code, err) \ + err = false; \ + \ +@@ -120,7 +120,7 @@ + err = true; \ + } \ + } \ +- catch (octave_execution_exception&) \ ++ catch (octave::execution_exception&) \ + { \ + err = true; \ + } +@@ -130,7 +130,7 @@ + { \ + code ; \ + } \ +- catch (octave_execution_exception&) \ ++ catch (octave::execution_exception&) \ + { \ + err = true; \ + }
--- a/src/of-signal-1-fixes.patch Mon Apr 09 07:11:54 2018 -0400 +++ b/src/of-signal-1-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -1,22 +1,202 @@ +diff -uNr a/src/cl2bp.cc b/src/cl2bp.cc +--- a/src/cl2bp.cc 2015-05-25 20:22:36.842410900 -0400 ++++ b/src/cl2bp.cc 2018-04-09 12:48:37.336964071 -0400 +@@ -84,27 +84,27 @@ + + const int m = args(0).int_value(true); + if (error_state) { +- gripe_wrong_type_arg("cl2bp", args (0)); ++ err_wrong_type_arg("cl2bp", args (0)); + return retval; + } + const double w1 = args(1).double_value(); + if (error_state) { +- gripe_wrong_type_arg("cl2bp", args (1)); ++ err_wrong_type_arg("cl2bp", args (1)); + return retval; + } + const double w2 = args(2).double_value(); + if (error_state) { +- gripe_wrong_type_arg("cl2bp", args (2)); ++ err_wrong_type_arg("cl2bp", args (2)); + return retval; + } + const ColumnVector up_vector(args(3).vector_value()); + if (error_state) { +- gripe_wrong_type_arg("cl2bp", args (3)); ++ err_wrong_type_arg("cl2bp", args (3)); + return retval; + } + const ColumnVector lo_vector(args(4).vector_value()); + if (error_state) { +- gripe_wrong_type_arg("cl2bp", args (4)); ++ err_wrong_type_arg("cl2bp", args (4)); + return retval; + } + if (up_vector.length() != 3 || lo_vector.length() != 3) { +@@ -121,7 +121,7 @@ + + const int L = args(5).int_value(true); + if (error_state) { +- gripe_wrong_type_arg("cl2bp", args (5)); ++ err_wrong_type_arg("cl2bp", args (5)); + return retval; + } + if (L > 1000000) { +diff -uNr a/src/medfilt1.cc b/src/medfilt1.cc +--- a/src/medfilt1.cc 2015-05-25 20:22:36.846411059 -0400 ++++ b/src/medfilt1.cc 2018-04-09 12:48:37.336964071 -0400 +@@ -186,14 +186,14 @@ + return retval; + } + +- if (args(0).is_complex_type()) ++ if (args(0).iscomplex()) + { + error("medfilt1 cannot process complex vectors"); + return retval; + } + + int n=3; // length of the filter (default 3) +- if (nargin > 1) n = NINT(args(1).double_value()); ++ if (nargin > 1) n = octave::math::nint(args(1).double_value()); + if (n < 1) + { + error ("medfilt1 filter length must be at least 1"); +diff -uNr a/src/remez.cc b/src/remez.cc +--- a/src/remez.cc 2015-05-25 20:22:36.846411059 -0400 ++++ b/src/remez.cc 2018-04-09 12:48:37.340963884 -0400 +@@ -784,7 +784,7 @@ + return retval; + } + +- int numtaps = NINT (args(0).double_value()) + 1; // #coeff = filter order+1 ++ int numtaps = octave::math::nint (args(0).double_value()) + 1; // #coeff = filter order+1 + if (numtaps < 4) { + error("remez: number of taps must be an integer greater than 3"); + return retval; +@@ -841,7 +841,7 @@ + if (args(4).is_string() && !args(3).is_string()) + stype = args(4).string_value(); + else if (args(4).is_real_scalar()) +- density = NINT(args(4).double_value()); ++ density = octave::math::nint(args(4).double_value()); + else { + error("remez: incorrect argument list"); + return retval; +@@ -850,7 +850,7 @@ + if (nargin > 5) { + if (args(5).is_real_scalar() + && !args(4).is_real_scalar()) +- density = NINT(args(5).double_value()); ++ density = octave::math::nint(args(5).double_value()); + else { + error("remez: incorrect argument list"); + return retval; diff -uNr a/src/sosfilt.cc b/src/sosfilt.cc --- a/src/sosfilt.cc 2015-05-25 20:22:36.846411059 -0400 -+++ b/src/sosfilt.cc 2018-01-09 15:50:16.018191246 -0500 -@@ -13,7 +13,6 @@ ++++ b/src/sosfilt.cc 2018-04-09 12:52:45.933287424 -0400 +@@ -13,10 +13,9 @@ // You should have received a copy of the GNU General Public License along with // this program; if not, see <http://www.gnu.org/licenses/>. -#include <octave/config.h> #include <octave/defun-dld.h> #include <octave/error.h> - #include <octave/gripes.h> +-#include <octave/gripes.h> ++#include <octave/errwarn.h> + #include <octave/oct-obj.h> + #include <octave/pager.h> + #include <octave/quit.h> +@@ -52,7 +51,7 @@ + + if (error_state) + { +- gripe_wrong_type_arg("sosfilt",args(0)); ++ err_wrong_type_arg("sosfilt",args(0)); + return retval; + } + +@@ -66,7 +65,7 @@ + + if (error_state) + { +- gripe_wrong_type_arg("sosfilt",args(1)); ++ err_wrong_type_arg("sosfilt",args(1)); + return retval; + } + diff -uNr a/src/upfirdn.cc b/src/upfirdn.cc --- a/src/upfirdn.cc 2015-05-25 20:22:36.846411059 -0400 -+++ b/src/upfirdn.cc 2018-01-09 15:50:13.514305277 -0500 -@@ -13,7 +13,6 @@ ++++ b/src/upfirdn.cc 2018-04-09 12:52:37.249695083 -0400 +@@ -13,10 +13,9 @@ // You should have received a copy of the GNU General Public License along with // this program; if not, see <http://www.gnu.org/licenses/>. -#include <octave/config.h> #include <octave/defun-dld.h> #include <octave/error.h> - #include <octave/gripes.h> +-#include <octave/gripes.h> ++#include <octave/errwarn.h> + #include <octave/oct-obj.h> + #include <octave/pager.h> + #include <octave/quit.h> +@@ -105,7 +104,7 @@ + + if (error_state) + { +- gripe_wrong_type_arg ("upfirdn", args (1)); ++ err_wrong_type_arg ("upfirdn", args (1)); + return retval; + } + +@@ -113,7 +112,7 @@ + + if (error_state) + { +- gripe_wrong_type_arg ("upfirdn", args (2)); ++ err_wrong_type_arg ("upfirdn", args (2)); + return retval; + } + +@@ -121,7 +120,7 @@ + + if (error_state) + { +- gripe_wrong_type_arg ("upfirdn", args (3)); ++ err_wrong_type_arg ("upfirdn", args (3)); + return retval; + } + +@@ -131,19 +130,19 @@ + Matrix x = args (0).matrix_value (); + if (error_state) + { +- gripe_wrong_type_arg ("upfirdn", args (0)); ++ err_wrong_type_arg ("upfirdn", args (0)); + return retval; + } + + Matrix y = upfirdn (x, h, p, q); + retval (0) = y; + } +- else if (args (0).is_complex_type ()) ++ else if (args (0).iscomplex ()) + { + ComplexMatrix x = args (0).complex_matrix_value (); + if (error_state) + { +- gripe_wrong_type_arg ("upfirdn", args (0)); ++ err_wrong_type_arg ("upfirdn", args (0)); + return retval; + } + +@@ -152,7 +151,7 @@ + } + else + { +- gripe_wrong_type_arg ("upfirdn", args (0)); ++ err_wrong_type_arg ("upfirdn", args (0)); + return retval; + } +
--- a/src/of-strings-1-fixes.patch Mon Apr 09 07:11:54 2018 -0400 +++ b/src/of-strings-1-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -14,3 +14,23 @@ pcregexp.oct: %.oct: %.cc $(MKOCTFILE) $(PCRE_SWITCHES) -o $@ $< +diff -uNr a/src/pcregexp.cc b/src/pcregexp.cc +--- a/src/pcregexp.cc 2015-06-06 17:40:45.000000000 -0400 ++++ b/src/pcregexp.cc 2018-04-09 15:01:27.359869527 -0400 +@@ -22,6 +22,7 @@ + // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + #include <octave/oct.h> ++#include <octave/errwarn.h> + #include <pcre.h> + #include <iostream> + #include <vector> +@@ -49,7 +50,7 @@ + std::string pattern = args(0).string_value(); + std::string input = args(1).string_value(); + if (error_state) { +- gripe_wrong_type_arg("pcregexp", args(0)); ++ err_wrong_type_arg("pcregexp", args(0)); + return retval; + } +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/of-struct-1-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -0,0 +1,117 @@ +diff -uNr a/src/error-helpers.cc b/src/error-helpers.cc +--- a/src/error-helpers.cc 2016-09-18 12:09:18.000000000 -0400 ++++ b/src/error-helpers.cc 2018-04-09 13:00:47.302704176 -0400 +@@ -24,7 +24,7 @@ + // call verror + #ifdef HAVE_OCTAVE_VERROR_ARG_EXC + void +-c_verror (octave_execution_exception& e, const char *fmt, ...) ++c_verror (octave::execution_exception& e, const char *fmt, ...) + { + va_list args; + va_start (args, fmt); +@@ -33,7 +33,7 @@ + } + #else + void +-c_verror (const octave_execution_exception&, const char *fmt, ...) ++c_verror (const octave::execution_exception&, const char *fmt, ...) + { + va_list args; + va_start (args, fmt); +diff -uNr a/src/error-helpers.h b/src/error-helpers.h +--- a/src/error-helpers.h 2016-09-18 12:09:18.000000000 -0400 ++++ b/src/error-helpers.h 2018-04-09 13:01:12.285532698 -0400 +@@ -21,9 +21,9 @@ + + // call verror + #ifdef HAVE_OCTAVE_VERROR_ARG_EXC +-void c_verror (octave_execution_exception&, const char *, ...); ++void c_verror (octave::execution_exception&, const char *, ...); + #else +-void c_verror (const octave_execution_exception&, const char *, ...); ++void c_verror (const octave::execution_exception&, const char *, ...); + #endif + + // call verror +@@ -33,7 +33,7 @@ + // both if Octave uses exceptions for errors and if it still uses + // error_state. In the latter case return 'retval'. + #ifdef HAVE_OCTAVE_ERROR_STATE +- // can throw octave_execution_exception despite of this ++ // can throw octave::execution_exception despite of this + #define CHECK_ERROR(code, retval, ...) \ + try \ + { \ +@@ -46,7 +46,7 @@ + return retval; \ + } \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -58,7 +58,7 @@ + { \ + code ; \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -70,7 +70,7 @@ + // Octave doesn't throw exceptions for errors but still uses + // error_state. + #ifdef HAVE_OCTAVE_ERROR_STATE +- // can throw octave_execution_exception despite of this ++ // can throw octave::execution_exception despite of this + #define CHECK_ERROR_EXIT1(code, ...) \ + try \ + { \ +@@ -83,7 +83,7 @@ + exit (1); \ + } \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -95,7 +95,7 @@ + { \ + code ; \ + } \ +- catch (octave_execution_exception& e) \ ++ catch (octave::execution_exception& e) \ + { \ + c_verror (e, __VA_ARGS__); \ + \ +@@ -107,7 +107,7 @@ + // Octave uses exceptions for errors and if it still uses + // error_state. In the latter case reset error_state to 0. + #ifdef HAVE_OCTAVE_ERROR_STATE +- // can throw octave_execution_exception despite of this ++ // can throw octave::execution_exception despite of this + #define SET_ERR(code, err) \ + err = false; \ + \ +@@ -120,7 +120,7 @@ + err = true; \ + } \ + } \ +- catch (octave_execution_exception&) \ ++ catch (octave::execution_exception&) \ + { \ + err = true; \ + } +@@ -130,7 +130,7 @@ + { \ + code ; \ + } \ +- catch (octave_execution_exception&) \ ++ catch (octave::execution_exception&) \ + { \ + err = true; \ + }
--- a/src/of-video-1-fixes.patch Mon Apr 09 07:11:54 2018 -0400 +++ b/src/of-video-1-fixes.patch Mon Apr 09 12:12:44 2018 -0400 @@ -123,3 +123,18 @@ AC_MSG_ERROR([mkoctfile required to install $PACKAGE_NAME]) fi +diff -uNr a/src/aviinfo.cc b/src/aviinfo.cc +--- a/src/aviinfo.cc 2017-05-28 10:27:02.824900272 -0400 ++++ b/src/aviinfo.cc 2018-04-09 15:26:03.042616699 -0400 +@@ -39,9 +39,9 @@ + // remove -Wno-deprecated-delarations from src/Makefile.in; remember + // to adjust 'Depends' to Octave >= 4.2 in DESCRIPTION) + +- octave_time file_mod = file_stat (fn).mtime (); ++ octave::sys::time file_mod = octave::sys::file_stat (fn).mtime (); + +- return octave_localtime (file_mod).strftime ("%d-%b-%Y %H:%M:%S"); ++ return octave::sys::localtime (file_mod).strftime ("%d-%b-%Y %H:%M:%S"); + } + + DEFUN_DLD(aviinfo, args, ,