Mercurial > forge
changeset 10226:44545240fca5 octave-forge
Fix formatting (whitespace changes only)
| author | jordigh |
|---|---|
| date | Fri, 11 May 2012 18:08:42 +0000 |
| parents | 12bc4b455216 |
| children | fe6da9839797 |
| files | main/comm/src/__errcore__.cc main/comm/src/__gfweight__.cc main/comm/src/cyclgen.cc main/comm/src/cyclpoly.cc main/comm/src/galois-def.cc main/comm/src/galois-def.h main/comm/src/galois-ops.h main/comm/src/galois.cc main/comm/src/galois.h main/comm/src/galoisfield.cc main/comm/src/galoisfield.h main/comm/src/genqamdemod.cc main/comm/src/gf.cc main/comm/src/isprimitive.cc main/comm/src/op-gm-gm.cc main/comm/src/op-gm-m.cc main/comm/src/op-gm-s.cc main/comm/src/op-m-gm.cc main/comm/src/op-s-gm.cc main/comm/src/ov-galois.cc main/comm/src/ov-galois.h main/comm/src/primpoly.cc main/comm/src/syndtable.cc |
| diffstat | 23 files changed, 2248 insertions(+), 2191 deletions(-) [+] |
line wrap: on
line diff
--- a/main/comm/src/__errcore__.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/__errcore__.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 <iomanip> @@ -44,18 +45,18 @@ error("_errcore: Matrix mismatch"); return retval; } - + unsigned int sz = (sizeof(unsigned int) << 3); Matrix c(a.rows(),a.cols(),0); for (int i=0; i<a.rows(); i++) for (int j=0; j<a.cols(); j++) { unsigned int tmp = (unsigned int)a(i,j) ^ (unsigned int)b(i,j); if (tmp != 0) - for (unsigned int k=0; k < sz; k++) - if (tmp & (1<<k)) - c(i,j)++; + for (unsigned int k=0; k < sz; k++) + if (tmp & (1<<k)) + c(i,j)++; } - + retval = octave_value(c); return retval; }
--- a/main/comm/src/__gfweight__.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/__gfweight__.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 <iomanip> @@ -25,25 +26,25 @@ #include <octave/quit.h> static int -get_weight (const Array<char>& codeword, const Matrix& gen, - int weight, int depth, int start, int n, int k) +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; + { + 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); - } + 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; } @@ -57,7 +58,8 @@ "This is an internal function of @dfn{gfweight}. You should use gfweight\n" "rather than use this function directly.\n" "@end deftypefn\n" -"@seealso{gfweight}") { +"@seealso{gfweight}") +{ if (args.length() != 1) { error("_gfweight: wrong number of arguments"); @@ -70,13 +72,14 @@ unsigned long long nsym = ((unsigned long long)1 << k); if (k > 128) - { - octave_stdout << "_gfweight: this is likely to take a very long time!!\n"; - flush_octave_stdout (); - } + { + 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)); + return octave_value((double)get_weight (codeword, gen, n - k + 1, 1, + 0, n, k)); } /*
--- a/main/comm/src/cyclgen.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/cyclgen.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 <iomanip> @@ -25,8 +26,8 @@ // 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) { +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); @@ -36,13 +37,13 @@ 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; + si(j) ^= 1; if (b(j+1) && x(i)) - si(j) ^= 1; + si(j) ^= 1; } si(n-1) = 0; if (a(n) && y(i)) @@ -54,9 +55,9 @@ 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). +// 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) { @@ -118,7 +119,7 @@ "matrix @var{g}." "\n" "@end deftypefn\n" -"@seealso{hammgen,gen2par,cyclpoly}") +"@seealso{hammgen,gen2par,cyclpoly}") { octave_value_list retval; int nargin = args.length (); @@ -155,26 +156,26 @@ 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; - } + 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; - } + 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; @@ -186,12 +187,12 @@ std::string s_arg = args(2).string_value (); if (s_arg == "system") - system = true; + system = true; else if (s_arg == "nosys") - system = false; + system = false; else { - error ("cyclgen: illegal argument"); - return retval; + error ("cyclgen: illegal argument"); + return retval; } } else { error ("cyclgen: illegal argument"); @@ -210,7 +211,7 @@ } unsigned long long mask = 1; - unsigned long long *alpha_to = + 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; @@ -220,10 +221,10 @@ } Matrix parity(mm,n,0); - for (int i = 0; i < n; i++) - for (int j = 0; j < mm; j++) + 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; + parity(j,i) = 1; free(alpha_to); retval(0) = octave_value (parity); @@ -231,10 +232,10 @@ 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++) + 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);
--- a/main/comm/src/cyclpoly.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/cyclpoly.cc Fri May 11 18:08:42 2012 +0000 @@ -1,17 +1,18 @@ //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 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. +// 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/>. +// 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 @@ -31,10 +32,10 @@ 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) { +// 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); @@ -44,13 +45,13 @@ 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; + si(j) ^= 1; if (b(j+1) && x(i)) - si(j) ^= 1; + si(j) ^= 1; } si(n-1) = 0; if (a(n) && y(i)) @@ -62,11 +63,12 @@ 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). +// 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) +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); @@ -153,18 +155,18 @@ std::string s_arg = args(i).string_value (); if (s_arg == "integer") - polyrep = false; + polyrep = false; else if (s_arg == "polynomial") - polyrep = true; + polyrep = true; else if (s_arg == "min") - type = CYCLIC_POLY_MIN; + type = CYCLIC_POLY_MIN; else if (s_arg == "max") - type = CYCLIC_POLY_MAX; + type = CYCLIC_POLY_MAX; else if (s_arg == "all") - type = CYCLIC_POLY_ALL; + type = CYCLIC_POLY_ALL; else { - error ("cyclpoly: invalid argument"); - return retval; + error ("cyclpoly: invalid argument"); + return retval; } } else { error ("cyclpoly: incorrect argument type"); @@ -176,7 +178,7 @@ // 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. + // should be reversed by replacing "i+=2" with "i-=2" and visa-versa. // Thats not going to happen!!! switch (type) { @@ -184,36 +186,37 @@ 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; + 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; + 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; + 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) { + 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 (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; - } + if (do_is_cyclic_polynomial(i, n, m)) { + cyclic_polys.resize(cyclic_polys.length()+1); + cyclic_polys(cyclic_polys.length()-1) = (double)i; + } } } break; @@ -223,16 +226,16 @@ } if (cyclic_polys.length() == 0) { - octave_stdout << + 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; + 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);
--- a/main/comm/src/galois-def.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/galois-def.cc Fri May 11 18:08:42 2012 +0000 @@ -22,10 +22,11 @@ void gripe_nonconformant_galois (const char *op, int op1_m, int op1_primpoly, - int op2_m, int op2_primpoly) + int op2_m, int op2_primpoly) { (*current_liboctave_error_handler) - ("%s: nonconformant arguments. op1 is GF(2^%d) (primitive polynomial %d), op2 is GF(2^%d) (primitive polynomial %d)", + ("%s: nonconformant arguments. op1 is GF(2^%d) " + "(primitive polynomial %d), op2 is GF(2^%d) (primitive polynomial %d)", op, op1_m, op1_primpoly, op2_m, op2_primpoly); }
--- a/main/comm/src/galois-def.h Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/galois-def.h Fri May 11 18:08:42 2012 +0000 @@ -1,27 +1,28 @@ //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 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. +// 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/>. +// 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) +// 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_defs_h) #define octave_galois_defs_h 1 void gripe_nonconformant_galois (const char *op, int op1_m, int op1_primpoly, - int op2_m, int op2_primpoly); + int op2_m, int op2_primpoly); void gripe_nonconformant_galois (const char *op, int m); void gripe_divzero_galois (const char *op); void gripe_invalid_galois (void); @@ -38,263 +39,263 @@ void gripe_init_galois (void); // Compute X % N, where N is 2^M - 1, without a slow divide -#define MODN(X, M, N) \ - { \ - while (X >= N) \ - { \ - X -= N; \ - X = (X >> M) + (X & N); \ - } \ +#define MODN(X, M, N) \ + { \ + while (X >= N) \ + { \ + X -= N; \ + X = (X >> M) + (X & N); \ + } \ } -#define CHECK_GALOIS(OP, RET, M1, M2, NN) \ - { \ - if (!M1.have_field() || !M2.have_field()) \ - { \ - gripe_invalid_galois(); \ - return RET (); \ - } \ - if ((M1.primpoly() != M2.primpoly()) || (M1.m() != M2.m())) \ - { \ - gripe_nonconformant_galois (OP, M1.m(), M1.primpoly(), M2.m(), M2.primpoly()); \ - return RET (); \ - } \ +#define CHECK_GALOIS(OP, RET, M1, M2, NN) \ + { \ + if (!M1.have_field() || !M2.have_field()) \ + { \ + gripe_invalid_galois(); \ + return RET (); \ + } \ + if ((M1.primpoly() != M2.primpoly()) || (M1.m() != M2.m())) \ + { \ + gripe_nonconformant_galois (OP, M1.m(), M1.primpoly(), M2.m(), M2.primpoly()); \ + return RET (); \ + } \ } -#define CHECK_MATRIX(OP, RET, M1, M2, NN) \ - { \ - int nr = M1.rows(); \ - int nc = M1.cols(); \ - \ - if (!M1.have_field()) \ - { \ - gripe_invalid_galois(); \ - return RET (); \ - } \ - for (int i=0; i<nr; i++) \ - for (int j=0; j<nc; j++) \ - { \ - if ((M1(i,j) < 0) || (M1(i,j) > NN)) \ - { \ - gripe_nonconformant_galois (OP, M1.m()); \ - return RET (); \ - } \ - if (((double)M1(i,j) - (double)((int)M1(i,j))) != 0.) \ - { \ - gripe_nonconformant_galois (OP, M1.m()); \ - return RET (); \ - } \ - } \ - } +#define CHECK_MATRIX(OP, RET, M1, M2, NN) \ + { \ + int nr = M1.rows(); \ + int nc = M1.cols(); \ + \ + if (!M1.have_field()) \ + { \ + gripe_invalid_galois(); \ + return RET (); \ + } \ + for (int i=0; i<nr; i++) \ + for (int j=0; j<nc; j++) \ + { \ + if ((M1(i,j) < 0) || (M1(i,j) > NN)) \ + { \ + gripe_nonconformant_galois (OP, M1.m()); \ + return RET (); \ + } \ + if (((double)M1(i,j) - (double)((int)M1(i,j))) != 0.) \ + { \ + gripe_nonconformant_galois (OP, M1.m()); \ + return RET (); \ + } \ + } \ + } -#define CHECK_DIV_ZERO(OP, RET, M) \ - { \ - int nr = M.rows(); \ - int nc = M.cols(); \ - \ - for (int i=0; i<nr; i++) \ - for (int j=0; j<nc; j++) \ - { \ - if (M(i,j) == 0) \ - { \ - gripe_divzero_galois(OP); \ - return RET (); \ - } \ - } \ - } +#define CHECK_DIV_ZERO(OP, RET, M) \ + { \ + int nr = M.rows(); \ + int nc = M.cols(); \ + \ + for (int i=0; i<nr; i++) \ + for (int j=0; j<nc; j++) \ + { \ + if (M(i,j) == 0) \ + { \ + gripe_divzero_galois(OP); \ + return RET (); \ + } \ + } \ + } #define CHECK_NODIV_ZERO(OP, RET, M) -#define MM_BIN_OP1(R, OP, M1, M2, GR1, GR2, CHECKTYPE) \ - R \ - OP (const M1& m1, const M2& m2) \ - { \ - R r (m ## GR1); \ - \ - int m1_nr = m1.rows (); \ - int m1_nc = m1.cols (); \ - \ - int m2_nr = m2.rows (); \ - int m2_nc = m2.cols (); \ - \ - CHECK_ ## CHECKTYPE (#OP, R, m ## GR1, m ## GR2, r.n()); \ - \ - if (m1_nr != m2_nr || m1_nc != m2_nc) \ - { \ - if ((m1_nr == 1 && m1_nc == 1) && (m2_nr > 0 && m2_nc > 0)) \ - { \ - r.resize(dim_vector(m2_nr,m2_nc)); \ - for (int i=0; i<m2_nr; i++) \ - for (int j=0; j<m2_nc; j++) \ - r(i,j) = (int)m1(0,0) ^ (int)m2(i,j); \ - } \ - else if ((m2_nr == 1 && m2_nc == 1) && (m1_nr > 0 && m1_nc > 0)) \ - { \ - r.resize(dim_vector(m1_nr,m1_nc)); \ - for (int i=0; i<m1_nr; i++) \ - for (int j=0; j<m1_nc; j++) \ - r(i,j) = (int)m1(i,j) ^ (int)m2(0,0); \ - } \ - else \ - gripe_nonconformant (#OP, m1_nr, m1_nc, m2_nr, m2_nc); \ - } \ - else \ - { \ - if (m1_nr > 0 && m1_nc > 0) \ - for (int i=0; i<m1_nr; i++) \ - for (int j=0; j<m1_nc; j++) \ - r(i,j) ^= (int) m ## GR2 (i,j); \ - } \ - \ - return r; \ +#define MM_BIN_OP1(R, OP, M1, M2, GR1, GR2, CHECKTYPE) \ + R \ + OP (const M1& m1, const M2& m2) \ + { \ + R r (m ## GR1); \ + \ + int m1_nr = m1.rows (); \ + int m1_nc = m1.cols (); \ + \ + int m2_nr = m2.rows (); \ + int m2_nc = m2.cols (); \ + \ + CHECK_ ## CHECKTYPE (#OP, R, m ## GR1, m ## GR2, r.n()); \ + \ + if (m1_nr != m2_nr || m1_nc != m2_nc) \ + { \ + if ((m1_nr == 1 && m1_nc == 1) && (m2_nr > 0 && m2_nc > 0)) \ + { \ + r.resize(dim_vector(m2_nr,m2_nc)); \ + for (int i=0; i<m2_nr; i++) \ + for (int j=0; j<m2_nc; j++) \ + r(i,j) = (int)m1(0,0) ^ (int)m2(i,j); \ + } \ + else if ((m2_nr == 1 && m2_nc == 1) && (m1_nr > 0 && m1_nc > 0)) \ + { \ + r.resize(dim_vector(m1_nr,m1_nc)); \ + for (int i=0; i<m1_nr; i++) \ + for (int j=0; j<m1_nc; j++) \ + r(i,j) = (int)m1(i,j) ^ (int)m2(0,0); \ + } \ + else \ + gripe_nonconformant (#OP, m1_nr, m1_nc, m2_nr, m2_nc); \ + } \ + else \ + { \ + if (m1_nr > 0 && m1_nc > 0) \ + for (int i=0; i<m1_nr; i++) \ + for (int j=0; j<m1_nc; j++) \ + r(i,j) ^= (int) m ## GR2 (i,j); \ + } \ + \ + return r; \ } #define MM_BIN_OP2(R, F, OP, M1, M2, NN, GR1, GR2, CHECKTYPE, ZEROCHECK) \ - R \ - F (const M1& m1, const M2& m2) \ - { \ - R r(m ## GR1); \ - \ - int m1_nr = m1.rows (); \ - int m1_nc = m1.cols (); \ - \ - int m2_nr = m2.rows (); \ - int m2_nc = m2.cols (); \ - \ - CHECK_ ## CHECKTYPE (#F, R, m ## GR1, m ## GR2, r.n()); \ - \ - CHECK_ ## ZEROCHECK ## DIV_ZERO (#F, R, m2); \ - \ - if (m1_nr != m2_nr || m1_nc != m2_nc) \ - { \ - if ((m1_nr == 1 && m1_nc == 1) && (m2_nr > 0 && m2_nc > 0)) \ - { \ - r.resize(dim_vector(m2_nr,m2_nc)); \ - if (m1(0,0) == 0) \ - { \ - for (int i=0; i<m2_nr; i++) \ - for (int j=0; j<m2_nc; j++) \ - r(i,j) = 0; \ - } \ - else \ - { \ - int indxm1 = r.index_of((int)m1(0,0)); \ - for (int i=0; i<m2_nr; i++) \ - for (int j=0; j<m2_nc; j++) \ - { \ - if (m2(i,j) == 0) \ - r(i,j) = 0; \ - else \ - { \ - r(i,j) = indxm1 OP r.index_of((int)m2(i,j)) + NN; \ - MODN(r(i,j), r.m(), r.n()); \ - r(i,j) = r.alpha_to(r(i,j)); \ - } \ - } \ - } \ - } \ - else if ((m2_nr == 1 && m2_nc == 1) && (m1_nr > 0 && m1_nc > 0)) \ - { \ - r.resize(dim_vector(m1_nr,m1_nc)); \ - if (m2(0,0) == 0) \ - { \ - for (int i=0; i<m1_nr; i++) \ - for (int j=0; j<m1_nc; j++) \ - r(i,j) = 0; \ - } \ - else \ - { \ - int indxm2 = r.index_of((int)m2(0,0)); \ - for (int i=0; i<m1_nr; i++) \ - for (int j=0; j<m1_nc; j++) \ - { \ - if (m1(i,j) == 0) \ - r(i,j) = 0; \ - else \ - { \ - r(i,j) = r.index_of((int)m1(i,j)) OP indxm2 + NN; \ - MODN(r(i,j), r.m(), r.n()); \ - r(i,j) = r.alpha_to(r(i,j)); \ - } \ - } \ - } \ - } \ - else \ - gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ - } \ - else \ - if (m1_nr > 0 && m1_nc > 0) \ - for (int i=0; i<m1_nr; i++) \ - for (int j=0; j<m1_nc; j++) \ - { \ - if ((m1(i,j) == 0) || (m2(i,j) == 0)) \ - r(i,j) = 0; \ - else \ - { \ - r(i,j) = r.index_of((int)m1(i,j)) OP r.index_of((int)m2(i,j)) + NN; \ - MODN(r(i,j), r.m(), r.n()); \ - r(i,j) = r.alpha_to(r(i,j)); \ - } \ - } \ - \ - return r; \ + R \ + F (const M1& m1, const M2& m2) \ + { \ + R r(m ## GR1); \ + \ + int m1_nr = m1.rows (); \ + int m1_nc = m1.cols (); \ + \ + int m2_nr = m2.rows (); \ + int m2_nc = m2.cols (); \ + \ + CHECK_ ## CHECKTYPE (#F, R, m ## GR1, m ## GR2, r.n()); \ + \ + CHECK_ ## ZEROCHECK ## DIV_ZERO (#F, R, m2); \ + \ + if (m1_nr != m2_nr || m1_nc != m2_nc) \ + { \ + if ((m1_nr == 1 && m1_nc == 1) && (m2_nr > 0 && m2_nc > 0)) \ + { \ + r.resize(dim_vector(m2_nr,m2_nc)); \ + if (m1(0,0) == 0) \ + { \ + for (int i=0; i<m2_nr; i++) \ + for (int j=0; j<m2_nc; j++) \ + r(i,j) = 0; \ + } \ + else \ + { \ + int indxm1 = r.index_of((int)m1(0,0)); \ + for (int i=0; i<m2_nr; i++) \ + for (int j=0; j<m2_nc; j++) \ + { \ + if (m2(i,j) == 0) \ + r(i,j) = 0; \ + else \ + { \ + r(i,j) = indxm1 OP r.index_of((int)m2(i,j)) + NN; \ + MODN(r(i,j), r.m(), r.n()); \ + r(i,j) = r.alpha_to(r(i,j)); \ + } \ + } \ + } \ + } \ + else if ((m2_nr == 1 && m2_nc == 1) && (m1_nr > 0 && m1_nc > 0)) \ + { \ + r.resize(dim_vector(m1_nr,m1_nc)); \ + if (m2(0,0) == 0) \ + { \ + for (int i=0; i<m1_nr; i++) \ + for (int j=0; j<m1_nc; j++) \ + r(i,j) = 0; \ + } \ + else \ + { \ + int indxm2 = r.index_of((int)m2(0,0)); \ + for (int i=0; i<m1_nr; i++) \ + for (int j=0; j<m1_nc; j++) \ + { \ + if (m1(i,j) == 0) \ + r(i,j) = 0; \ + else \ + { \ + r(i,j) = r.index_of((int)m1(i,j)) OP indxm2 + NN; \ + MODN(r(i,j), r.m(), r.n()); \ + r(i,j) = r.alpha_to(r(i,j)); \ + } \ + } \ + } \ + } \ + else \ + gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ + } \ + else \ + if (m1_nr > 0 && m1_nc > 0) \ + for (int i=0; i<m1_nr; i++) \ + for (int j=0; j<m1_nc; j++) \ + { \ + if ((m1(i,j) == 0) || (m2(i,j) == 0)) \ + r(i,j) = 0; \ + else \ + { \ + r(i,j) = r.index_of((int)m1(i,j)) OP r.index_of((int)m2(i,j)) + NN; \ + MODN(r(i,j), r.m(), r.n()); \ + r(i,j) = r.alpha_to(r(i,j)); \ + } \ + } \ + \ + return r; \ } -#define MM_BIN_OPS1(R, M1, M2, GR1, GR2, CHECK) \ +#define MM_BIN_OPS1(R, M1, M2, GR1, GR2, CHECK) \ MM_BIN_OP1(R, operator +, M1, M2, GR1, GR2, CHECK) \ - MM_BIN_OP1(R, operator -, M1, M2, GR1, GR2, CHECK) \ - MM_BIN_OP2(R, product, +, M1, M2, 0, GR1, GR2, CHECK, NO) \ + MM_BIN_OP1(R, operator -, M1, M2, GR1, GR2, CHECK) \ + MM_BIN_OP2(R, product, +, M1, M2, 0, GR1, GR2, CHECK, NO) \ MM_BIN_OP2(R, quotient, -, M1, M2, r.n(), GR1, GR2, CHECK, ) -#define MM_BIN_OPS2(R, M1, M2, GR1, GR2, CHECK) \ +#define MM_BIN_OPS2(R, M1, M2, GR1, GR2, CHECK) \ MM_BIN_OP1(R, operator +, M1, M2, GR1, GR2, CHECK) -#define MM_CMP_OP1(F, OP, M1, C1, M2, C2, GR1, GR2, CHECKTYPE) \ - boolMatrix \ - F (const M1& m1, const M2& m2) \ - { \ - boolMatrix r; \ - \ - int m1_nr = m1.rows (); \ - int m1_nc = m1.cols (); \ - \ - int m2_nr = m2.rows (); \ - int m2_nc = m2.cols (); \ - \ +#define MM_CMP_OP1(F, OP, M1, C1, M2, C2, GR1, GR2, CHECKTYPE) \ + boolMatrix \ + F (const M1& m1, const M2& m2) \ + { \ + boolMatrix r; \ + \ + int m1_nr = m1.rows (); \ + int m1_nc = m1.cols (); \ + \ + int m2_nr = m2.rows (); \ + int m2_nc = m2.cols (); \ + \ CHECK_ ## CHECKTYPE (#F, boolMatrix, m ## GR1, m ## GR2, m ## GR1.n()); \ - \ - if (m1_nr == m2_nr && m1_nc == m2_nc) \ - { \ - r.resize (m1_nr, m1_nc); \ - \ - for (int j = 0; j < m1_nc; j++) \ - for (int i = 0; i < m1_nr; i++) \ - r(i, j) = C1 (m1(i, j)) OP C2 (m2(i, j)); \ - } \ - else \ - { \ - if ((m1_nr == 1 && m1_nc == 1) && (m2_nr > 0 && m2_nc > 0)) \ - { \ - r.resize(m2_nr,m2_nc); \ - for (int i=0; i<m2_nr; i++) \ - for (int j=0; j<m2_nc; j++) \ - r(i, j) = C1 (m1(0, 0)) OP C2 (m2(i, j)); \ - } \ - else if ((m2_nr == 1 && m2_nc == 1) && (m1_nr > 0 && m1_nc > 0)) \ - { \ - r.resize(m1_nr,m1_nc); \ - for (int i=0; i<m1_nr; i++) \ - for (int j=0; j<m1_nc; j++) \ - r(i, j) = C1 (m1(i, j)) OP C2 (m2(0, 0)); \ - } \ - else \ - gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ - } \ - \ - return r; \ + \ + if (m1_nr == m2_nr && m1_nc == m2_nc) \ + { \ + r.resize (m1_nr, m1_nc); \ + \ + for (int j = 0; j < m1_nc; j++) \ + for (int i = 0; i < m1_nr; i++) \ + r(i, j) = C1 (m1(i, j)) OP C2 (m2(i, j)); \ + } \ + else \ + { \ + if ((m1_nr == 1 && m1_nc == 1) && (m2_nr > 0 && m2_nc > 0)) \ + { \ + r.resize(m2_nr,m2_nc); \ + for (int i=0; i<m2_nr; i++) \ + for (int j=0; j<m2_nc; j++) \ + r(i, j) = C1 (m1(0, 0)) OP C2 (m2(i, j)); \ + } \ + else if ((m2_nr == 1 && m2_nc == 1) && (m1_nr > 0 && m1_nc > 0)) \ + { \ + r.resize(m1_nr,m1_nc); \ + for (int i=0; i<m1_nr; i++) \ + for (int j=0; j<m1_nc; j++) \ + r(i, j) = C1 (m1(i, j)) OP C2 (m2(0, 0)); \ + } \ + else \ + gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ + } \ + \ + return r; \ } -#define MM_CMP_OPS1(M1, C1, M2, C2, GR1, GR2, CHECK) \ +#define MM_CMP_OPS1(M1, C1, M2, C2, GR1, GR2, CHECK) \ MM_CMP_OP1(mx_el_lt, <, M1, C1, M2, C2, GR1, GR2, CHECK) \ MM_CMP_OP1(mx_el_le, <=, M1, C1, M2, C2, GR1, GR2, CHECK) \ MM_CMP_OP1(mx_el_ge, >=, M1, C1, M2, C2, GR1, GR2, CHECK) \ @@ -303,102 +304,102 @@ MM_CMP_OP1(mx_el_ne, !=, M1, , M2, , GR1, GR2, CHECK) \ #define MM_BOOL_OP1(F, OP, M1, M2, ZERO, GR1, GR2, CHECKTYPE) \ - boolMatrix \ - F (const M1& m1, const M2& m2) \ - { \ - boolMatrix r; \ - \ - int m1_nr = m1.rows (); \ - int m1_nc = m1.cols (); \ - \ - int m2_nr = m2.rows (); \ - int m2_nc = m2.cols (); \ - \ + boolMatrix \ + F (const M1& m1, const M2& m2) \ + { \ + boolMatrix r; \ + \ + int m1_nr = m1.rows (); \ + int m1_nc = m1.cols (); \ + \ + int m2_nr = m2.rows (); \ + int m2_nc = m2.cols (); \ + \ CHECK_ ## CHECKTYPE (#F, boolMatrix, m ## GR1, m ## GR2, m ## GR1.n()); \ - \ - if (m1_nr == m2_nr && m1_nc == m2_nc) \ - { \ - if (m1_nr != 0 || m1_nc != 0) \ - { \ - r.resize (m1_nr,m1_nc); \ - \ - for (int j = 0; j < m1_nc; j++) \ - for (int i = 0; i < m1_nr; i++) \ - { \ - r(i, j) = (m1(i, j) != ZERO) \ - OP (m2(i, j) != ZERO); \ - } \ - } \ - } \ - else \ - { \ - if ((m1_nr == 1 && m1_nc == 1) && (m2_nr > 0 && m2_nc > 0)) \ - { \ - r.resize(m2_nr,m2_nc); \ - for (int i=0; i<m2_nr; i++) \ - for (int j=0; j<m2_nc; j++) \ - r(i, j) = (m1(0, 0) != ZERO) \ - OP (m2(i, j) != ZERO); \ - } \ - else if ((m2_nr == 1 && m2_nc == 1) && (m1_nr > 0 && m1_nc > 0)) \ - { \ - r.resize(m1_nr,m1_nc); \ - for (int i=0; i<m1_nr; i++) \ - for (int j=0; j<m1_nc; j++) \ - r(i, j) = (m1(i, j) != ZERO) \ - 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); \ - } \ - \ - return r; \ + \ + if (m1_nr == m2_nr && m1_nc == m2_nc) \ + { \ + if (m1_nr != 0 || m1_nc != 0) \ + { \ + r.resize (m1_nr,m1_nc); \ + \ + for (int j = 0; j < m1_nc; j++) \ + for (int i = 0; i < m1_nr; i++) \ + { \ + r(i, j) = (m1(i, j) != ZERO) \ + OP (m2(i, j) != ZERO); \ + } \ + } \ + } \ + else \ + { \ + if ((m1_nr == 1 && m1_nc == 1) && (m2_nr > 0 && m2_nc > 0)) \ + { \ + r.resize(m2_nr,m2_nc); \ + for (int i=0; i<m2_nr; i++) \ + for (int j=0; j<m2_nc; j++) \ + r(i, j) = (m1(0, 0) != ZERO) \ + OP (m2(i, j) != ZERO); \ + } \ + else if ((m2_nr == 1 && m2_nc == 1) && (m1_nr > 0 && m1_nc > 0)) \ + { \ + r.resize(m1_nr,m1_nc); \ + for (int i=0; i<m1_nr; i++) \ + for (int j=0; j<m1_nc; j++) \ + r(i, j) = (m1(i, j) != ZERO) \ + 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); \ + } \ + \ + return r; \ } -#define MM_BOOL_OPS1(M1, M2, ZERO, GR1, GR2, CHECK) \ +#define MM_BOOL_OPS1(M1, M2, ZERO, GR1, GR2, CHECK) \ MM_BOOL_OP1 (mx_el_and, &&, M1, M2, ZERO, GR1, GR2, CHECK) \ MM_BOOL_OP1 (mx_el_or, ||, M1, M2, ZERO, GR1, GR2, CHECK) -#define GALOIS_REDUCTION_OP(RET, ROW_EXPR, COL_EXPR, INIT_VAL, \ - MT_RESULT) \ - \ - int nr = rows (); \ - int nc = cols (); \ - \ - if (nr > 0 && nc > 0) \ - { \ - if ((nr == 1 && dim == -1) || dim == 1) \ - { \ - RET.resize (dim_vector(nr, 1)); \ - for (int i = 0; i < nr; i++) \ - { \ - RET (i, 0) = INIT_VAL; \ - for (int j = 0; j < nc; j++) \ - { \ - ROW_EXPR; \ - } \ - } \ - } \ - else \ - { \ - RET.resize (dim_vector(1, nc)); \ - for (int j = 0; j < nc; j++) \ - { \ - RET (0, j) = INIT_VAL; \ - for (int i = 0; i < nr; i++) \ - { \ - COL_EXPR; \ - } \ - } \ - } \ - } \ - else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \ - RET.resize (dim_vector(1, 1), MT_RESULT); \ - else if (nr == 0 && (dim == 0 || dim == -1)) \ - RET.resize (dim_vector(1, nc), MT_RESULT); \ - else if (nc == 0 && dim == 1) \ - RET.resize (dim_vector(nr, 1), MT_RESULT); \ - else \ +#define GALOIS_REDUCTION_OP(RET, ROW_EXPR, COL_EXPR, INIT_VAL, \ + MT_RESULT) \ + \ + int nr = rows (); \ + int nc = cols (); \ + \ + if (nr > 0 && nc > 0) \ + { \ + if ((nr == 1 && dim == -1) || dim == 1) \ + { \ + RET.resize (dim_vector(nr, 1)); \ + for (int i = 0; i < nr; i++) \ + { \ + RET (i, 0) = INIT_VAL; \ + for (int j = 0; j < nc; j++) \ + { \ + ROW_EXPR; \ + } \ + } \ + } \ + else \ + { \ + RET.resize (dim_vector(1, nc)); \ + for (int j = 0; j < nc; j++) \ + { \ + RET (0, j) = INIT_VAL; \ + for (int i = 0; i < nr; i++) \ + { \ + COL_EXPR; \ + } \ + } \ + } \ + } \ + else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \ + RET.resize (dim_vector(1, 1), MT_RESULT); \ + else if (nr == 0 && (dim == 0 || dim == -1)) \ + RET.resize (dim_vector(1, nc), MT_RESULT); \ + else if (nc == 0 && dim == 1) \ + RET.resize (dim_vector(nr, 1), MT_RESULT); \ + else \ RET.resize (dim_vector(nr > 0, nc > 0)); #endif
--- a/main/comm/src/galois-ops.h Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/galois-ops.h Fri May 11 18:08:42 2012 +0000 @@ -1,117 +1,120 @@ //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 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. +// 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/>. +// 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) +// 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 (galois_octave_ops_h) #define galois_octave_ops_h 1 // Override the operator and function definition defines from Octave -#define DEFBINOP_OP_G(name, t1, t2, op) \ - BINOPDECL (name, a1, a2) \ - { \ +#define DEFBINOP_OP_G(name, t1, t2, op) \ + BINOPDECL (name, a1, a2) \ + { \ CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ - return new octave_galois \ - (v1.t1 ## _value () op v2.t2 ## _value ()); \ + return new octave_galois \ + (v1.t1 ## _value () op v2.t2 ## _value ()); \ } -#define DEFBINOP_FN_G(name, t1, t2, f) \ - BINOPDECL (name, a1, a2) \ - { \ - CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ +#define DEFBINOP_FN_G(name, t1, t2, f) \ + BINOPDECL (name, a1, a2) \ + { \ + CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ return new octave_galois (f (v1.t1 ## _value (), v2.t2 ## _value ())); \ } -#define DEFBINOP_OP_B_S1(name, t1, t2, op) \ - BINOPDECL (name, a1, a2) \ - { \ +#define DEFBINOP_OP_B_S1(name, t1, t2, op) \ + BINOPDECL (name, a1, a2) \ + { \ CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ - return octave_value \ - (v1.matrix_value () op v2.t2 ##_value ()); \ + return octave_value \ + (v1.matrix_value () op v2.t2 ##_value ()); \ } -#define DEFBINOP_FN_B_S1(name, t1, t2, f) \ - BINOPDECL (name, a1, a2) \ - { \ - CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ +#define DEFBINOP_FN_B_S1(name, t1, t2, f) \ + BINOPDECL (name, a1, a2) \ + { \ + CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ return octave_value (f (v1.matrix_value (), v2.t2 ## _value ())); \ } -#define DEFBINOP_OP_G_S1(name, t1, t2, op) \ - BINOPDECL (name, a1, a2) \ - { \ +#define DEFBINOP_OP_G_S1(name, t1, t2, op) \ + BINOPDECL (name, a1, a2) \ + { \ CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ - return new octave_galois \ - (v1.matrix_value () op v2.t2 ## _value ()); \ + return new octave_galois \ + (v1.matrix_value () op v2.t2 ## _value ()); \ } -#define DEFBINOP_FN_G_S1(name, t1, t2, f) \ - BINOPDECL (name, a1, a2) \ - { \ - CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ +#define DEFBINOP_FN_G_S1(name, t1, t2, f) \ + BINOPDECL (name, a1, a2) \ + { \ + CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ return new octave_galois (f (v1.matrix_value (), v2.t2 ## _value ())); \ } -#define DEFBINOP_OP_B_S2(name, t1, t2, op) \ - BINOPDECL (name, a1, a2) \ - { \ +#define DEFBINOP_OP_B_S2(name, t1, t2, op) \ + BINOPDECL (name, a1, a2) \ + { \ CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ - return octave_value \ - (v1.t1 ## _value () op v2.matrix_value ()); \ + return octave_value \ + (v1.t1 ## _value () op v2.matrix_value ()); \ } -#define DEFBINOP_FN_B_S2(name, t1, t2, f) \ - BINOPDECL (name, a1, a2) \ - { \ - CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ +#define DEFBINOP_FN_B_S2(name, t1, t2, f) \ + BINOPDECL (name, a1, a2) \ + { \ + CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ return octave_value (f (v1.t1 ## _value (), v2.matrix_value ())); \ } -#define DEFBINOP_OP_G_S2(name, t1, t2, op) \ - BINOPDECL (name, a1, a2) \ - { \ +#define DEFBINOP_OP_G_S2(name, t1, t2, op) \ + BINOPDECL (name, a1, a2) \ + { \ CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ - return new octave_galois \ - (v1.t1 ## _value () op v2.matrix_value ()); \ + return new octave_galois \ + (v1.t1 ## _value () op v2.matrix_value ()); \ } -#define DEFBINOP_FN_G_S2(name, t1, t2, f) \ - BINOPDECL (name, a1, a2) \ - { \ - CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ +#define DEFBINOP_FN_G_S2(name, t1, t2, f) \ + BINOPDECL (name, a1, a2) \ + { \ + CAST_BINOP_ARGS (const octave_ ## t1&, const octave_ ## t2&); \ return new octave_galois (f (v1.t1 ## _value (), v2.matrix_value ())); \ } -#define DEFCATOP_G_FN(name, t1, t2, f) \ - CATOPDECL (name, a1, a2) \ - { \ - CAST_BINOP_ARGS (octave_ ## t1&, const octave_ ## t2&); \ - return new octave_galois (f (v1.t1 ## _value (), v2.t2 ## _value (), ra_idx)); \ +#define DEFCATOP_G_FN(name, t1, t2, f) \ + CATOPDECL (name, a1, a2) \ + { \ + CAST_BINOP_ARGS (octave_ ## t1&, const octave_ ## t2&); \ + return new octave_galois (f (v1.t1 ## _value (), v2.t2 ## _value (), \ + ra_idx)); \ } -#define DEFCATOP_G_METHOD(name, t1, t2, f) \ - CATOPDECL (name, a1, a2) \ - { \ - CAST_BINOP_ARGS (octave_ ## t1&, const octave_ ## t2&); \ - return new octave_galois (v1.t1 ## _value (). f (v2.t2 ## _value (), ra_idx)); \ +#define DEFCATOP_G_METHOD(name, t1, t2, f) \ + CATOPDECL (name, a1, a2) \ + { \ + CAST_BINOP_ARGS (octave_ ## t1&, const octave_ ## t2&); \ + return new octave_galois (v1.t1 ## _value (). f (v2.t2 ## _value (),\ + ra_idx)); \ } -#define INSTALL_G_CATOP(t1, t2, f) INSTALL_CATOP(t1, t2, f) +#define INSTALL_G_CATOP(t1, t2, f) INSTALL_CATOP(t1, t2, f) #endif
--- a/main/comm/src/galois.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/galois.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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> @@ -29,15 +30,17 @@ // galois class -galois::galois (const Array<int>& a, const int& _m, const int& _primpoly) : MArray<int> (a.dims ()), field (NULL) { +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; + gripe_range_galois(_m); + return; } xelem(i,j) = (int)a(i,j); } @@ -46,15 +49,16 @@ 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) { +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; + gripe_range_galois(_m); + return; } xelem(i,j) = (int)a(i,j); } @@ -63,19 +67,21 @@ 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) { +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; + gripe_range_galois(_m); + return; } if ((a(i,j) - (double)((int)a(i,j))) != 0.) { - gripe_integer_galois(); - return; + gripe_integer_galois(); + return; } xelem(i,j) = (int)a(i,j); } @@ -84,7 +90,10 @@ 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) { +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 @@ -96,7 +105,10 @@ 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) { +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 @@ -133,7 +145,7 @@ } galois & galois::operator = (const galois& t) -{ +{ if (!t.have_field()) { gripe_copy_invalid_galois(); if (have_field()) @@ -228,14 +240,14 @@ } galois galois::index (idx_vector& i, idx_vector& j, int resize_ok, - const int& rfv) const + const int& rfv) const { galois retval(MArray<int>::index(i, j, resize_ok, rfv), m(), primpoly()); return retval; } -galois +galois galois::concat (const galois& rb, const Array<int>& ra_idx) { if (rb.numel() > 0) @@ -243,7 +255,7 @@ return *this; } -galois +galois galois::concat (const Matrix& rb, const Array<int>& ra_idx) { if (numel() == 1) @@ -259,12 +271,12 @@ 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; + gripe_range_galois(m()); + return *this; } if ((rb(i,j) - (double)((int)rb(i,j))) != 0.) { - gripe_integer_galois(); - return *this; + gripe_integer_galois(); + return *this; } tmp(i,j) = (int)rb(i,j); } @@ -292,22 +304,22 @@ 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; + gripe_range_galois(rb.m()); + return retval; } if ((ra(i,j) - (double)((int)ra(i,j))) != 0.) { - gripe_integer_galois(); - return retval; + 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; + retval(i,j) = 0; else if (tmp > _n) - retval(i,j) = _n; + retval(i,j) = _n; else - retval(i,j) = tmp; + retval(i,j) = tmp; #endif } } @@ -319,7 +331,8 @@ 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"); + (*current_liboctave_error_handler) ("inserted galois variable must " + "be in the same field"); else Array<int>::insert (t, r, c); return *this; @@ -342,27 +355,27 @@ nnr += k; if (nnr > 0 && nnc > 0) - { - int ndiag = (nnr < nnc) ? nnr : nnc; - retval.resize(dim_vector(ndiag, 1)); + { + 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); - } + if (k > 0) + { + for (int i = 0; i < ndiag; i++) + retval (i,0) = xelem (i, i+k); } - else + 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; @@ -396,7 +409,7 @@ for (int j = 0; j < d2; j++) for (int i = 0; i < d1; i++) a (j, i) = xelem (i, j); - + return a; } @@ -429,42 +442,42 @@ } 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())); - } + { + 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) { - 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())); + gripe_nonconformant ("operator .^", a_nr, a_nc, a_nr, a_nc); + return galois (); } - else - { - if (a_nr != b_nr || a_nc != b_nc) - { - gripe_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())); - } + 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; } @@ -480,25 +493,25 @@ if (b_nr == 1 && b_nc == 1) return elem_pow(a, b(0,0)); - + if (a_nr != b_nr || a_nc != b_nc) - { - gripe_nonconformant ("operator .^", a_nr, a_nc, b_nr, b_nc); - return galois (); - } + { + gripe_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())); - } + { + 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; } @@ -519,13 +532,13 @@ 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())); - } + { + 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; } @@ -540,13 +553,13 @@ 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())); - } + { + 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; } @@ -595,36 +608,36 @@ } 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; - } + 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; + { + galois atmp; - if (b < 0 ) { - atmp = a.inverse(); - b = abs(b); - } else - atmp = a; + if (b < 0 ) { + atmp = a.inverse(); + b = abs(b); + } else + atmp = a; - retval = atmp; - b--; - while (b > 0) - { - if (b & 1) - retval = retval * atmp; + retval = atmp; + b--; + while (b > 0) + { + if (b & 1) + retval = retval * atmp; - b >>= 1; + b >>= 1; - if (b > 0) - atmp = atmp * atmp; - } + if (b > 0) + atmp = atmp * atmp; } + } return retval; } @@ -640,7 +653,7 @@ } galois -operator * (const galois& a, const Matrix& b) +operator * (const galois& a, const Matrix& b) { galois tmp (b, a.m(), a.primpoly()); @@ -671,38 +684,39 @@ if ((a_nr == 1 && a_nc == 1) || (b_nr == 1 && b_nc == 1)) return product (a, b); else if (a_nc != b_nr) - { - gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc); - return galois(); - } + { + gripe_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) { - 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)); + // 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; + 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 +// Other operators boolMatrix galois::all (int dim) const { @@ -725,19 +739,19 @@ galois retval (0, 0, 0, m(), primpoly()); -#define ROW_EXPR \ +#define ROW_EXPR \ if ((retval (i, 0) == 0) || (elem(i,j) == 0)) \ - retval (i, 0) = 0; \ - else \ + retval (i, 0) = 0; \ + else \ retval (i, 0) = alpha_to(modn(index_of(retval(i,0)) + \ - index_of(elem(i,j)),m(),n())); + index_of(elem(i,j)),m(),n())); -#define COL_EXPR \ +#define COL_EXPR \ if ((retval (0, j) == 0) || (elem(i,j) == 0)) \ - retval (0, j) = 0; \ - else \ + retval (0, j) = 0; \ + else \ retval (0, j) = alpha_to(modn(index_of(retval(0,j)) + \ - index_of(elem(i,j)),m(),n())); + index_of(elem(i,j)),m(),n())); GALOIS_REDUCTION_OP(retval, ROW_EXPR, COL_EXPR, 1, 1); return retval; @@ -757,10 +771,10 @@ galois retval (0, 0, 0, m(), primpoly()); -#define ROW_EXPR \ +#define ROW_EXPR \ retval (i, 0) ^= elem(i,j); -#define COL_EXPR \ +#define COL_EXPR \ retval (0, j) ^= elem(i,j); GALOIS_REDUCTION_OP (retval, ROW_EXPR, COL_EXPR, 0, 0); @@ -780,12 +794,12 @@ galois retval (0, 0, 0, m(), primpoly()); -#define ROW_EXPR \ - if (elem(i,j) != 0) \ +#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) \ +#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); @@ -805,19 +819,19 @@ 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); + 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); - } + retval(i,j) = retval.alpha_to(retval.index_of(retval(i,j)) + / 2); + } return retval; } galois galois::log (void) const { - bool warned = false; + bool warned = false; if (!have_field()) { gripe_invalid_galois (); return galois(); @@ -829,17 +843,17 @@ 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)); + 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)); + retval(i,j) = retval.index_of(retval(i,j)); } return retval; } @@ -847,7 +861,7 @@ galois galois::exp (void) const { - bool warned = false; + bool warned = false; if (!have_field()) { gripe_invalid_galois (); return galois(); @@ -860,16 +874,16 @@ 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)); + 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)); + retval(i,j) = retval.alpha_to(retval(i,j)); } return retval; } @@ -893,64 +907,65 @@ int jp = j; // Find the pivot and test for singularity - if (ptype == galoisLU::ROW) { + if (ptype == galoisLU::ROW) { for (int i = j+1; i < a_nr; i++) - if (a_fact(i,j) > a_fact(jp,j)) - jp = 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; + 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; - } + 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; - } + // 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())); - } + 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())); - } - } + 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())); + } + } } } } @@ -1007,17 +1022,17 @@ if (nr != nc || nr == 0 || nc == 0) { (*current_liboctave_error_handler) ("inverse requires square matrix"); return galois (); - } else { + } else { int info = 0; - // Solve with identity matrix to find the inverse. + // 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) + if (info == 0) return retval; else return galois(); @@ -1052,14 +1067,14 @@ 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())); + 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; @@ -1080,7 +1095,7 @@ galois galois::solve (const galois& b, int& info, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler) const { galois retval (b); @@ -1110,9 +1125,9 @@ if (fact.singular()) { info = -1; if (sing_handler) - sing_handler (0.0); + sing_handler (0.0); else - (*current_liboctave_error_handler)("galois matrix singular"); + (*current_liboctave_error_handler)("galois matrix singular"); return galois(); } else { @@ -1121,61 +1136,61 @@ // Resize the number of solution rows if needed if (nc > nr) - retval.resize(dim_vector(b_nr+nc-nr, b_nc), 0); + 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 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())); - } + 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())); - } + 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; - } + 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); + sing_handler (0.0); else - (*current_liboctave_error_handler)("galois matrix singular"); + (*current_liboctave_error_handler)("galois matrix singular"); return galois(); } else { @@ -1184,48 +1199,48 @@ // 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; - } + 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())); - } + 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())); - } + 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)); + retval.resize(dim_vector(b_nr+nc-nr, b_nc)); } } @@ -1258,19 +1273,19 @@ // if ((a_nc != b_nc) || (b.rows () != b.cols ())) if (a_nc != b_nc) - { - int a_nr = a.rows (); - int b_nr = b.rows (); + { + int a_nr = a.rows (); + int b_nr = b.rows (); - gripe_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc); - return galois(); - } + gripe_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) + if (info == 0) return galois (result.transpose ()); else return galois(); @@ -1302,13 +1317,13 @@ // if ((a_nr != b_nr) || (a.rows () != a.columns ())) if (a_nr != b_nr) - { - int a_nc = a.cols (); - int b_nc = b.cols (); + { + int a_nc = a.cols (); + int b_nc = b.cols (); - gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc); - return galois(); - } + gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc); + return galois(); + } galois result = a.solve (b, info, 0);
--- a/main/comm/src/galois.h Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/galois.h Fri May 11 18:08:42 2012 +0000 @@ -58,15 +58,17 @@ galois (const Array<int>& a, const int& m=1, const int& primpoly=0); galois (const MArray<int>& a, const int& m=1, const int& primpoly=0); galois (const Matrix& a, const int& m=1, const int& primpoly=0); - galois (int nr, int nc, const int& val=0, const int& _m=1, const int& _primpoly=0); - galois (int nr, int nc, double val=0., const int& _m=1, const int& _primpoly=0); + galois (int nr, int nc, const int& val=0, const int& _m=1, + const int& _primpoly=0); + galois (int nr, int nc, double val=0., const int& _m=1, + const int& _primpoly=0); galois (const galois& a); ~galois (void); - + galois index (idx_vector& i, int resize_ok = 0, const int& rfv = 0) const; galois index (idx_vector& i, idx_vector& j, int resize_ok = 0, - const int& rfv = 0) const; + const int& rfv = 0) const; // unary operations @@ -80,8 +82,8 @@ galois concat (const galois& rb, const Array<int>& ra_idx); galois concat (const Matrix& rb, const Array<int>& ra_idx); - friend galois concat (const Matrix& ra, const galois& rb, - const Array<int>& ra_idx); + friend galois concat (const Matrix& ra, const galois& rb, + const Array<int>& ra_idx); galois& insert (const galois& a, int r, int c); @@ -99,20 +101,20 @@ galois solve (const galois& b) const; galois solve (const galois& b, int& info) const; galois solve (const galois& b, int& info, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler) const; galois determinant (void) const; galois determinant (int& info) const; - galois &operator = (const galois& t); - galois &operator += (const galois& a); - galois &operator -= (const galois& a); - - private: + galois &operator = (const galois& t); + galois &operator += (const galois& a); + galois &operator -= (const galois& a); + +private: // Pointer to the Galois field structure used galois_field_node *field; - - public: + +public: // Is the variable initialized?? bool have_field (void) const { return (field ? true : false); }; @@ -145,12 +147,12 @@ galoisLU (const galoisLU& a) : base_lu <galois> (a) { } galoisLU& operator = (const galoisLU& a) - { - if (this != &a) - base_lu <galois> :: operator = (a); + { + if (this != &a) + base_lu <galois> :: operator = (a); - return *this; - } + return *this; + } ~galoisLU (void) { } @@ -158,7 +160,7 @@ galois U (void) const; - bool singular (void) const { return info != 0; } + bool singular (void) const { return info != 0; } pivot_type type (void) const { return ptype; } private:
--- a/main/comm/src/galoisfield.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/galoisfield.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 "galois.h"
--- a/main/comm/src/galoisfield.h Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/galoisfield.h Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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_field_int_h) #define octave_galois_field_int_h 1
--- a/main/comm/src/genqamdemod.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/genqamdemod.cc Fri May 11 18:08:42 2012 +0000 @@ -1,17 +1,18 @@ //Copyright (C) 2006 Charalampos C. Tsimenidis // -// 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 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. +// 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/>. +// 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/oct.h> #include <octave/ov-fcn.h> @@ -25,83 +26,83 @@ @end deftypefn") { -octave_value retval; -int i, j, m; -double tmp1, tmp2; + octave_value retval; + int i, j, m; + double tmp1, tmp2; -if (args.length()<2) -{ - print_usage (); - return retval; -} + if (args.length()<2) + { + print_usage (); + return retval; + } -int nr1 (args(0).rows()); -int nc1 (args(0).columns()); -int arg_is_empty1=empty_arg("genqamdemod", nr1, nc1); -Matrix y(nr1,nc1); + int nr1 (args(0).rows()); + int nc1 (args(0).columns()); + int arg_is_empty1=empty_arg("genqamdemod", nr1, nc1); + Matrix y(nr1,nc1); -int nr2 (args(1).rows()); -int nc2 (args(1).columns()); -int M=(nr2>nc2)?nr2:nc2; + int nr2 (args(1).rows()); + int nc2 (args(1).columns()); + int M=(nr2>nc2)?nr2:nc2; -if (arg_is_empty1 < 0) + if (arg_is_empty1 < 0) return retval; -if (arg_is_empty1 > 0) + if (arg_is_empty1 > 0) return octave_value (Matrix ()); -if (args(0).is_real_type() && args(1).is_real_type()) -{ // Real-valued signal & constellation - Matrix x(args(0).matrix_value()); - ColumnVector constellation(args(1).vector_value()); - for (i=0;i<nr1;i++) - { - for (j=0;j<nc1;j++) - { - tmp1=fabs(x(i,j)-constellation(0)); - y(i,j)=0; - for (m=1;m<M;m++) - { - tmp2=fabs(x(i,j)-constellation(m)); - if (tmp2<tmp1) - { - y(i,j)=m; - tmp1=tmp2; - } - } - } - } -} -else if (args(0).is_complex_type() || args(1).is_complex_type()) -{ // Complex-valued input & constellation - ComplexMatrix x (args(0).complex_matrix_value()); - ComplexColumnVector constellation(args(1).complex_vector_value()); - if (!error_state) - { - for (i=0;i<nr1;i++) - { - for (j=0;j<nc1;j++) - { - tmp1=abs(x(i,j)-constellation(0)); - y(i,j)=0; - for (m=1;m<M;m++) - { - tmp2=abs(x(i,j)-constellation(m)); - if (tmp2<tmp1) - { - y(i,j)=m; - tmp1=tmp2; - } - - } - } - } - } - else - print_usage (); + if (args(0).is_real_type() && args(1).is_real_type()) + { // Real-valued signal & constellation + Matrix x(args(0).matrix_value()); + ColumnVector constellation(args(1).vector_value()); + for (i=0;i<nr1;i++) + { + for (j=0;j<nc1;j++) + { + tmp1=fabs(x(i,j)-constellation(0)); + y(i,j)=0; + for (m=1;m<M;m++) + { + tmp2=fabs(x(i,j)-constellation(m)); + if (tmp2<tmp1) + { + y(i,j)=m; + tmp1=tmp2; + } + } + } + } + } + else if (args(0).is_complex_type() || args(1).is_complex_type()) + { // Complex-valued input & constellation + ComplexMatrix x (args(0).complex_matrix_value()); + ComplexColumnVector constellation(args(1).complex_vector_value()); + if (!error_state) + { + for (i=0;i<nr1;i++) + { + for (j=0;j<nc1;j++) + { + tmp1=abs(x(i,j)-constellation(0)); + y(i,j)=0; + for (m=1;m<M;m++) + { + tmp2=abs(x(i,j)-constellation(m)); + if (tmp2<tmp1) + { + y(i,j)=m; + tmp1=tmp2; + } + + } + } + } + } + else + print_usage (); + } + else + { + print_usage (); + } + return retval=y; } -else -{ - print_usage (); -} -return retval=y; -}
--- a/main/comm/src/gf.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/gf.cc Fri May 11 18:08:42 2012 +0000 @@ -2,22 +2,23 @@ // 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 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. +// 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/>. +// 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) +// 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 @@ -41,17 +42,17 @@ "-*- 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") +"@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(); + 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(); } // XXX FIXME XXX @@ -104,14 +105,14 @@ } 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; - mlock (); + 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; + mlock (); } retval = new octave_galois(data, m, primpoly); @@ -123,51 +124,51 @@ { octave_value retval; - if ((!galois_type_loaded) || (a.type_id () != - octave_galois::static_type_id ())) + if ((!galois_type_loaded) || (a.type_id () != + octave_galois::static_type_id ())) gripe_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 (); + { + 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 == 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.capacity () > 0) - retval = new octave_galois (r); - } + 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.capacity () > 0) + retval = new octave_galois (r); } + } else gripe_wrong_type_arg ("gdiag", a); } @@ -260,8 +261,8 @@ } else { int mr = 0, mc = 0; - if ((!galois_type_loaded) || (args(0).type_id () != - octave_galois::static_type_id ())) { + if ((!galois_type_loaded) || (args(0).type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("greshape", args(0)); return retval; } @@ -274,7 +275,7 @@ } else if (nargin == 3) { mr = args(1).nint_value (); mc = args(2).nint_value (); - } + } int nr = a.rows(); int nc = a.cols(); @@ -283,55 +284,55 @@ 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); + 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); + for (int j=0;j<mc;j++) + tmp2(i,j) = (int)tmp1(i+j*mr); retval = new octave_galois(tmp2); } } return retval; } -#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 \ - { \ - gripe_wrong_type_arg (#FCN, arg); \ - return retval; \ - } \ - } \ - else \ - error (#FCN ": invalid dimension argument = %d", dim + 1); \ - } \ - } \ - else \ - print_usage (); \ - \ +#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 \ + { \ + gripe_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"); @@ -392,8 +393,8 @@ return retval; } - if (!galois_type_loaded || (args(0).type_id () != - octave_galois::static_type_id ())) { + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("gsqrt", args(0)); return retval; } @@ -422,8 +423,8 @@ return retval; } - if (!galois_type_loaded || (args(0).type_id () != - octave_galois::static_type_id ())) { + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("glog", args(0)); return retval; } @@ -453,8 +454,8 @@ return retval; } - if (!galois_type_loaded || (args(0).type_id () != - octave_galois::static_type_id ())) { + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("gexp", args(0)); return retval; } @@ -493,8 +494,8 @@ int idx_norm = b.index_of(norm); for (int i=0; i < b.length(); 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())); + 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) { @@ -503,74 +504,74 @@ if (norm != 1) { int idx_norm = a.index_of(norm); for (int i=0; i < a.length(); 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())); + 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++) { 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())); + 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) { - for (int j = 0; j < si.length() - 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.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))+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))+ b.index_of(x(i,0)), - b.m(),b.n())); + for (int j = 0; j < si.length() - 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.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))+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))+ 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())); + 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.length() > 0) { for (int i = 0; i < x.length(); 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())); + 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) { - for (int j = 0; j < si.length() - 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)) + b.index_of(x(i,0)),b.m(),b.n())); + for (int j = 0; j < si.length() - 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)) + 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())); + 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.length(); 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())); + 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; } @@ -696,7 +697,7 @@ if (!ib) { if (ia) b = galois(args(0).matrix_value(), a.m(), a.primpoly()); - else if (ix) + 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()); @@ -705,7 +706,7 @@ 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()); @@ -714,12 +715,12 @@ int b_len = b.length (); 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() != a.primpoly()) || (b.primpoly() != x.primpoly()) || (b.primpoly() != si.primpoly())) { error("gfilter: arguments must be in same galois field"); return retval; @@ -746,7 +747,7 @@ else retval(1) = new octave_galois (si); } - + if (x_is_row_vector) retval(0) = new octave_galois (y.transpose ()); else @@ -803,20 +804,20 @@ "@seealso{lu}") { octave_value_list retval; - + int nargin = args.length (); if (nargin != 1 || nargout > 3) - { - print_usage (); - return retval; - } + { + print_usage (); + return retval; + } octave_value arg = args(0); - if (!galois_type_loaded || (arg.type_id () != - octave_galois::static_type_id ())) { + if (!galois_type_loaded || (arg.type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("glu", arg); return retval; } @@ -845,15 +846,15 @@ 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); + // 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 (); @@ -862,7 +863,7 @@ break; } } - + return retval; } @@ -881,18 +882,18 @@ int nargin = args.length (); if (nargin != 1) - { - print_usage (); - return retval; - } + { + 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 ())) { + if (!galois_type_loaded || (arg.type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("ginverse", arg); return retval; } @@ -904,30 +905,30 @@ 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(0) = new octave_galois (galois(0, 0, 0, m.m(), m.primpoly())); return retval; } if (nr != nc) - { - gripe_square_matrix_required ("ginverse"); - return retval; - } + { + gripe_square_matrix_required ("ginverse"); + return retval; + } if (! error_state) - { - int info; - double rcond = 0.0; + { + int info; + double rcond = 0.0; - galois result = m.inverse (info, 1); + galois result = m.inverse (info, 1); - if (nargout > 1) - retval(1) = rcond; + if (nargout > 1) + retval(1) = rcond; - retval(0) = new octave_galois (result); + retval(0) = new octave_galois (result); - if (nargout < 2 && info == -1) - warning ("inverse: matrix singular to machine precision, rcond = %g", rcond); - } + if (nargout < 2 && info == -1) + warning ("inverse: matrix singular to machine precision, rcond = %g", rcond); + } return retval; } @@ -965,8 +966,8 @@ octave_value arg = args(0); - if (!galois_type_loaded || (arg.type_id () != - octave_galois::static_type_id ())) { + if (!galois_type_loaded || (arg.type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("gdet", arg); return retval; } @@ -981,7 +982,7 @@ 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())); + retval = new octave_galois (galois(1, 1, 1, m.m(), m.primpoly())); return retval; } @@ -1014,8 +1015,8 @@ octave_value arg = args(0); - if (!galois_type_loaded || (arg.type_id () != - octave_galois::static_type_id ())) { + if (!galois_type_loaded || (arg.type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("grank", arg); return retval; } @@ -1039,31 +1040,31 @@ 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; - } - } + { + 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)); + 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)); - } - } + 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; @@ -1127,14 +1128,14 @@ { 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 ())) { + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("rsenc", args(0)); return retval; } @@ -1181,35 +1182,35 @@ 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]); - + parstr[j] = toupper(parstr[j]); + if (!parstr.compare("END")) { - parity_at_end = true; + parity_at_end = true; } else if (!parstr.compare("BEGINNING")) { - parity_at_end = false; + parity_at_end = false; } else { - error ("rsenc: unrecoginized parity position"); - return retval; + 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 (have_genpoly) { + print_usage (); + return retval; + } + genpoly = ((const octave_galois&) args(i).get_rep()).galois_value (); - if (genpoly.cols() > genpoly.rows()) - genpoly = genpoly.transpose(); + 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(); + if (have_genpoly) { + if (prim != 0) { + print_usage (); + return retval; + } + prim = args(i).nint_value(); + } else + fcr = args(i).nint_value(); } have_genpoly = true; } @@ -1231,16 +1232,17 @@ // 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); + 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)); + genpoly(nroots,0) = genpoly + .alpha_to(modn(genpoly.index_of(genpoly(nroots,0)) + root, m, n)); } } else { @@ -1250,7 +1252,8 @@ } if (genpoly.primpoly() != primpoly) { - error ("rsenc: the generator polynomial must be same galois field as the message"); + error ("rsenc: the generator polynomial must be same galois field " + "as the message"); return retval; } @@ -1269,67 +1272,67 @@ // Add space for parity block msg.resize(dim_vector (nsym, n), 0); - // The code below basically finds the parity bits by treating the + // 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 + // 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 + // 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 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); + 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); + msg(l,i+nroots-1) = msg(l,i-1); for (int i=0; i<nroots; i++) - msg(l,i) = 0; + 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 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); + msg(l,j) = par(nroots-j-1,0); } } @@ -1339,7 +1342,7 @@ } int decode_rs(galois& data, const int prim, const int iprim, const int nroots, - const int fcr, const int drow, const bool msb_first) + const int fcr, const int drow, const bool msb_first) { int deg_lambda, el, deg_omega; int i, j, r, k; @@ -1352,7 +1355,7 @@ /* 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); @@ -1368,22 +1371,22 @@ 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)); + 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)); + 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 */ @@ -1404,22 +1407,22 @@ 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 */ + 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.alpha_to(modn(data.index_of(lambda[i]) + + s[r-i-1], m, n)); } } - discr_r = data.index_of(discr_r); /* Index form */ + 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])); @@ -1428,24 +1431,24 @@ /* 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(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); + 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; + /* 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])); } @@ -1461,13 +1464,13 @@ /* 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) */ + 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]); + reg[j] = modn(reg[j] + j, m, n); + q ^= data.alpha_to(reg[j]); } } if (q != 0) @@ -1498,14 +1501,14 @@ 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)); + 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 @@ -1514,16 +1517,16 @@ 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)); + 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) { + 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)); + den ^= data.alpha_to(modn(lambda[i+1] + i * root[j], m, n)); } if (den == 0) { count = -1; @@ -1532,11 +1535,15 @@ /* 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)); + 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)); + 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)); } } @@ -1587,14 +1594,14 @@ 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 ())) { + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) { gripe_wrong_type_arg ("rsdec", args(0)); return retval; } @@ -1616,7 +1623,8 @@ } if (code.m() != m) { - error ("rsdec: coded message in incorrect galois field for codeword length"); + error ("rsdec: coded message in incorrect galois field for " + "codeword length"); return retval; } @@ -1641,36 +1649,36 @@ 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; + 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; - } + 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 (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; } } } @@ -1678,129 +1686,131 @@ if (have_genpoly) { if (fcr != 0) { if ((fcr < 1) || (fcr > nn)) { - error("rsdec: invalid first consecutive root of generator polynomial"); - return retval; + 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; + error("rsdec: invalid primitive element of generator polynomial"); + return retval; } } else { if (genpoly.cols() > genpoly.rows()) - genpoly = genpoly.transpose(); + genpoly = genpoly.transpose(); if (genpoly.cols() != 1) { - error ("rsdec: the generator polynomial must be a vector"); - return retval; + 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; + 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; + 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; - } + // 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; + 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]); + roots[i] = genpoly.index_of(roots[i]); // Find a corresponding fcr and prim that coincide with the roots. // XXX FIXME XXX. 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; + 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) + 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 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) + 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); + 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); + 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); + msg(l,i) = code(l,nroots+i); } retval(0) = new octave_galois (msg); @@ -1852,7 +1862,7 @@ { octave_value retval; int nargin = args.length (); - + if ((nargin < 3) || (nargin > 5)) { print_usage (); return retval; @@ -1887,30 +1897,30 @@ 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]); - + parstr[j] = toupper(parstr[j]); + if (!parstr.compare("END")) { - parity_at_end = true; + parity_at_end = true; } else if (!parstr.compare("BEGINNING")) { - parity_at_end = false; + parity_at_end = false; } else { - error ("bchenco: unrecoginized parity position"); - return retval; + 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(); + genpoly = genpoly.transpose(); if (genpoly.cols() != 1) { - error ("bchenco: the generator polynomial must be a vector"); - return retval; + 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; + error ("bchenco: generator polynomial has incorrect order"); + return retval; } } } @@ -1931,19 +1941,19 @@ 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); + 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; + 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; + c(nc,cs(nc)) = r; + found(c.alpha_to(r)-1) = 1; + cs(nc) += 1; } nc++; } @@ -1957,23 +1967,23 @@ 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; - } - } + 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); @@ -1981,7 +1991,7 @@ error("bchenco: can not find valid generator polynomial for parameters"); return retval; } - + // Create polynomial of right length. genpoly = galois(nf+1,1,0,m); @@ -1991,70 +2001,72 @@ // 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); + 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)); + 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 + // 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. + // The parity bits are then the remainder of this division. // - // This code could just as easily be written as + // 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; - } + 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); + msg(l,i+nn-k-1) = msg(l,i-1); for (int i=0; i<nn-k; i++) - msg(l,i) = 0; + 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; - } + 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; + } } } } @@ -2120,7 +2132,7 @@ { octave_value_list retval; int nargin = args.length (); - + if ((nargin < 3) || (nargin > 5)) { print_usage (); return retval; @@ -2144,41 +2156,42 @@ return retval; } - int prim = 0; // primitve polynomial of zero flags default + 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]); - + parstr[j] = toupper(parstr[j]); + if (!parstr.compare("END")) { - parity_at_end = true; + parity_at_end = true; } else if (!parstr.compare("BEGINNING")) { - parity_at_end = false; + parity_at_end = false; } else { - error ("bchdeco: unrecoginized parity position"); - return retval; + 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 (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() > tmp.rows()) + tmp = tmp.transpose(); - if (tmp.cols() != 1) { - error ("bchdeco: the primitve polynomial must be a scalar or a vector"); - return retval; - } + 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); + prim = 0; + for (int i=0; i < tmp.rows(); i++) + if ((int)tmp(i,0) & 1) + prim |= (1<<i); } } } @@ -2195,28 +2208,28 @@ 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 (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 */ + syn_error = true; /* set error flag if non-zero syndrome */ } - if (syn_error) { /* if there are errors, try to correct them */ + 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)); + 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)); + s(i) = tables.index_of(s(i)); /* * Compute the error location polynomial via the Berlekamp @@ -2225,130 +2238,131 @@ * 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. + * 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 */ + 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 */ + 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)); + 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); + } - 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)); - } + /* + * 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)); + 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)); + /* 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; - } - } + /* 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; + 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; + nerr(lsym) = -1; } } @@ -2356,13 +2370,13 @@ if (parity_at_end) { for (int l = 0; l < nsym; l++) for (int i = 0; i < k; i++) - msg(l,i) = code(l,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); + msg(l,i) = code(l,nn-k+i); } - + retval(0) = octave_value(msg); retval(1) = octave_value(nerr); retval(2) = octave_value(code);
--- a/main/comm/src/isprimitive.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/isprimitive.cc Fri May 11 18:08:42 2012 +0000 @@ -1,14 +1,14 @@ //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 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. +// 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/>. @@ -52,7 +52,7 @@ } DEFUN_DLD (isprimitive, args, nargout, - "-*- texinfo -*-\n" +"-*- texinfo -*-\n" "@deftypefn {Loadable Function} {@var{y} =} isprimitive (@var{a})\n" "\n" "Returns 1 is the polynomial represented by @var{a} is a primitive\n" @@ -71,13 +71,13 @@ } // If only 0/1 in a, assume that each row is a polynomial representation - + bool poly = true; for (int i=0; i<a.rows(); i++) { for (int j=0; j<a.columns(); j++) { if (((int)a(i,j) != 0) && ((int)a(i,j) != 1)) { - poly = false; - break; + poly = false; + break; } } if (!poly) break; @@ -85,41 +85,43 @@ if (poly) { if (a.columns() > 24) { - error("isprimitive: order of the primitive polynomial must be less than 22"); + error("isprimitive: order of the primitive polynomial must " + "be less than 22"); return retval; } Matrix b(a.rows(),1); - + for (int i=0; i<a.rows(); i++) { int tmp = (int)a(i,0); for (int j=1; j<a.columns(); j++) - tmp = (tmp << 1) | (int)a(i,j); - + tmp = (tmp << 1) | (int)a(i,j); + int m = 1; while (tmp > (1<<(m+1))) - m++; - + m++; + b(i,0) = do_isprimitive(tmp, m); } retval = octave_value (b); } else { for (int i=0; i<a.rows(); i++) for (int j=0; j<a.columns(); j++) - if (a(i,j) > (1<<23)) { - error("isprimitive: order of the primitive polynomial must be less than 22"); - return retval; - } + if (a(i,j) > (1<<23)) { + error("isprimitive: order of the primitive polynomial must " + "be less than 22"); + return retval; + } Matrix b(a.rows(),a.columns()); for (int i=0; i<a.rows(); i++) { for (int j=0; j<a.columns(); j++) { - int m = 1; - while (a(i,j) > (1<<(m+1))) - m++; - - b(i,j) = do_isprimitive((int)a(i,j), m); + int m = 1; + while (a(i,j) > (1<<(m+1))) + m++; + + b(i,j) = do_isprimitive((int)a(i,j), m); } } retval = octave_value (b);
--- a/main/comm/src/op-gm-gm.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/op-gm-gm.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 "galois.h" @@ -51,7 +52,7 @@ // galois by galois ops. DEFBINOP_OP_G (add, galois, galois, +) -DEFBINOP_OP_G (sub, galois, galois, -) +DEFBINOP_OP_G (sub, galois, galois, -) DEFBINOP_OP_G (mul, galois, galois, *) DEFBINOP_FN_G (div, galois, galois, xdiv) DEFBINOP_FN_G (pow, galois, galois, pow)
--- a/main/comm/src/op-gm-m.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/op-gm-m.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 "galois.h" @@ -62,13 +63,13 @@ DEFCATOP_G_METHOD (gm_m, galois, matrix, concat) // Need to create temporary Galois array so that matrix values are checked -DEFASSIGNOP (assign, galois, matrix) +DEFASSIGNOP (assign, galois, matrix) { - CAST_BINOP_ARGS (octave_galois&, const octave_matrix&); + CAST_BINOP_ARGS (octave_galois&, const octave_matrix&); - v1.assign (idx, galois(v2.matrix_value (), v1.galois_value().m(), - v1.galois_value().primpoly())); - return octave_value (); + v1.assign (idx, galois(v2.matrix_value (), v1.galois_value().m(), + v1.galois_value().primpoly())); + return octave_value (); } void
--- a/main/comm/src/op-gm-s.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/op-gm-s.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 "galois.h" @@ -62,17 +63,17 @@ DEFCATOP (gm_s, galois, scalar) { CAST_BINOP_ARGS (octave_galois&, const octave_scalar&); - return new octave_galois (v1.galois_value (). concat (v2.matrix_value (), - ra_idx)); + return new octave_galois (v1.galois_value (). concat (v2.matrix_value (), + ra_idx)); } DEFASSIGNOP(assign, galois, scalar) { - CAST_BINOP_ARGS (octave_galois&, const octave_scalar&); + CAST_BINOP_ARGS (octave_galois&, const octave_scalar&); - v1.assign (idx, galois(1,1,v2.scalar_value (), v1.galois_value().m(), - v1.galois_value().primpoly())); - return octave_value (); + v1.assign (idx, galois(1,1,v2.scalar_value (), v1.galois_value().m(), + v1.galois_value().primpoly())); + return octave_value (); } void
--- a/main/comm/src/op-m-gm.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/op-m-gm.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 "galois.h"
--- a/main/comm/src/op-s-gm.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/op-s-gm.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 "galois.h" @@ -75,8 +76,8 @@ DEFCATOP (s_gm, scalar, galois) { CAST_BINOP_ARGS (octave_scalar&, const octave_galois&); - return new octave_galois (concat (v1.matrix_value (), v2.galois_value (), - ra_idx)); + return new octave_galois (concat (v1.matrix_value (), v2.galois_value (), + ra_idx)); } void
--- a/main/comm/src/ov-galois.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/ov-galois.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 <octave/config.h> @@ -31,14 +32,14 @@ DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA(octave_galois, "galois", "galois"); octave_value octave_galois::resize (const dim_vector& dv, bool) const -{ +{ if (dv.length() > 2) - { - error ("Can not resize galois structure to NDArray"); - return octave_value (); - } - galois retval (gval); - retval.resize (dv); + { + error ("Can not resize galois structure to NDArray"); + return octave_value (); + } + galois retval (gval); + retval.resize (dv); return new octave_galois (retval); } @@ -61,9 +62,9 @@ Matrix data(r,c); for (int i=0; i< r; i++) for (int j=0; j< c; j++) - data(i,j) = (double)gval(i,j); + data(i,j) = (double)gval(i,j); retval(0) = octave_value (data); - } + } #ifdef GALOIS_DISP_PRIVATES else if (nm == __GALOIS_LENGTH_STR) retval(0) = octave_value((double)gval.n()); @@ -82,80 +83,80 @@ } #endif else - error ("galois structure has no member `%s'", nm.c_str ()); + error ("galois structure has no member `%s'", nm.c_str ()); return retval; } octave_value octave_galois::do_index_op (const octave_value_list& idx, - bool resize_ok) + bool resize_ok) { octave_value retval; int len = idx.length (); switch (len) + { + case 2: { - case 2: - { - idx_vector i = idx (0).index_vector (); - idx_vector j = idx (1).index_vector (); - - retval = new octave_galois (gval.index (i, j, resize_ok)); - } - break; + idx_vector i = idx (0).index_vector (); + idx_vector j = idx (1).index_vector (); - case 1: - { - idx_vector i = idx (0).index_vector (); + retval = new octave_galois (gval.index (i, j, resize_ok)); + } + break; + + case 1: + { + idx_vector i = idx (0).index_vector (); - retval = new octave_galois (gval.index (i, resize_ok)); - } - break; + retval = new octave_galois (gval.index (i, resize_ok)); + } + break; - default: - { - std::string n = type_name (); + default: + { + std::string n = type_name (); - error ("invalid number of indices (%d) for %s value", - len, n.c_str ()); - } - break; + error ("invalid number of indices (%d) for %s value", + len, n.c_str ()); } + break; + } return retval; } octave_value octave_galois::subsref (const std::string &type, - const std::list<octave_value_list>& idx) + const std::list<octave_value_list>& idx) { octave_value retval; int skip = 1; switch (type[0]) + { + case '(': + retval = do_index_op (idx.front (), true); + break; + + case '.': { - case '(': - retval = do_index_op (idx.front (), true); - break; - - case '.': - { - octave_value_list t = dotref (idx.front ()); + octave_value_list t = dotref (idx.front ()); - retval = (t.length () == 1) ? t(0) : octave_value (t); - } - break; + retval = (t.length () == 1) ? t(0) : octave_value (t); + } + break; - case '{': - error ("%s cannot be indexed with %c", type_name().c_str(), type[0]); - break; + case '{': + error ("%s cannot be indexed with %c", type_name().c_str(), type[0]); + break; - default: - panic_impossible (); - } + default: + panic_impossible (); + } if (! error_state) retval = retval.next_subsref (type, idx, skip); @@ -169,11 +170,11 @@ bool retval = false; if (rows () > 0 && columns () > 0) - { - boolMatrix m = (gval.all () . all ()); + { + boolMatrix m = (gval.all () . all ()); - retval = (m.rows () == 1 && m.columns () == 1 && m(0,0)); - } + retval = (m.rows () == 1 && m.columns () == 1 && m(0,0)); + } return retval; } @@ -210,21 +211,21 @@ for (int i=m; i>=0; i--) { if (primpoly & (1<<i)) { - if (i > 0) { - if (first) { - first = false; - os << "D"; - } else - os << "+D"; - if (i != 1) - os << "^" << i; - } else { - if (first) { - first = false; - os << "1"; - } else - os << "+1"; - } + if (i > 0) { + if (first) { + first = false; + os << "D"; + } else + os << "+D"; + if (i != 1) + os << "^" << i; + } else { + if (first) { + first = false; + os << "1"; + } else + os << "+1"; + } } } os << " (decimal " << primpoly << ")"; @@ -257,10 +258,10 @@ #if 0 if (Vstruct_levels_to_print < 0) #else - if (false) + if (false) #endif - os << name << " = "; - else + os << name << " = "; + else { os << name << " ="; newline (os); @@ -282,12 +283,12 @@ double retval = octave_NaN; if (rows () > 0 && columns () > 0) - { - gripe_implicit_conversion ("Octave:array-as-scalar", - "real matrix", "real scalar"); + { + gripe_implicit_conversion ("Octave:array-as-scalar", + "real matrix", "real scalar"); - retval = (double) gval (0, 0); - } + retval = (double) gval (0, 0); + } else gripe_invalid_conversion ("galois", "real scalar"); @@ -302,12 +303,12 @@ Complex retval (tmp, tmp); if (rows () > 0 && columns () > 0) - { - gripe_implicit_conversion ("Octave:array-as-scalar", - "real matrix", "real scalar"); + { + gripe_implicit_conversion ("Octave:array-as-scalar", + "real matrix", "real scalar"); - retval = (double) gval (0, 0); - } + retval = (double) gval (0, 0); + } else gripe_invalid_conversion ("galois", "complex scalar"); @@ -344,7 +345,7 @@ void octave_galois::assign (const octave_value_list& idx, - const galois& rhs) + const galois& rhs) { int len = idx.length (); @@ -356,38 +357,38 @@ } switch (len) + { + case 2: { - case 2: + idx_vector i = idx (0).index_vector (); + + if (! error_state) { - idx_vector i = idx (0).index_vector (); + idx_vector j = idx (1).index_vector (); if (! error_state) - { - idx_vector j = idx (1).index_vector (); + gval.assign (i, j, rhs); + } + } + break; - if (! error_state) - gval.assign (i, j, rhs); - } - } - break; + case 1: + { + idx_vector i = idx (0).index_vector (); - case 1: - { - idx_vector i = idx (0).index_vector (); + if (! error_state) + gval.assign (i, rhs); + } + break; - if (! error_state) - gval.assign (i, rhs); - } - break; - - default: - error ("invalid number of indices (%d) for galois assignment", - len); - break; - } + default: + error ("invalid number of indices (%d) for galois assignment", + len); + break; + } } -bool +bool octave_galois::save_ascii (std::ostream& os) { dim_vector d = dims (); @@ -406,7 +407,7 @@ return true; } -bool +bool octave_galois::load_ascii (std::istream& is) { int mord, prim, mdims; @@ -414,42 +415,42 @@ if (extract_keyword (is, "m", mord) && extract_keyword (is, "prim", prim) && extract_keyword (is, "ndims", mdims)) - { - dim_vector dv; - dv.resize (mdims); - - for (int i = 0; i < mdims; i++) - is >> dv(i); - - if (dv.length() != 2) - { - error ("load: N-D galois matrix not supported"); - success = false; - } - else - { + { + dim_vector dv; + dv.resize (mdims); - Matrix tmp (dv(0), dv(1)); - is >> tmp; + for (int i = 0; i < mdims; i++) + is >> dv(i); - if (!is) - { - error ("load: failed to load matrix constant"); - success = false; - } - gval = galois (tmp, mord, prim); - } - } - else + if (dv.length() != 2) { - error ("load: failed to extract galois field order, primitive and/or imension"); + error ("load: N-D galois matrix not supported"); success = false; } + else + { + + Matrix tmp (dv(0), dv(1)); + is >> tmp; + + if (!is) + { + error ("load: failed to load matrix constant"); + success = false; + } + gval = galois (tmp, mord, prim); + } + } + else + { + error ("load: failed to extract galois field order, primitive and/or imension"); + success = false; + } return success;; } -bool +bool octave_galois::save_binary (std::ostream& os, bool& save_as_floats) { char tmp = m (); @@ -467,10 +468,10 @@ itmp = - d.length(); os.write (X_CAST (char *, &itmp), 4); for (int i=0; i < d.length (); i++) - { - itmp = d(i); - os.write (X_CAST (char *, &itmp), 4); - } + { + itmp = d(i); + os.write (X_CAST (char *, &itmp), 4); + } Matrix m = matrix_value (); save_type st; @@ -486,9 +487,9 @@ return true; } -bool +bool octave_galois::load_binary (std::istream& is, bool swap, - oct_mach_info::float_format fmt) + oct_mach_info::float_format fmt) { char mord; int32_t prim, mdims; @@ -508,33 +509,33 @@ // Don't treat N-D arrays yet if (mdims == -2) - { - mdims = - mdims; - int32_t di; - dim_vector dv; - dv.resize (mdims); + { + mdims = - mdims; + int32_t di; + dim_vector dv; + dv.resize (mdims); - for (int i = 0; i < mdims; i++) - { - if (! is.read (X_CAST (char *, &di), 4)) - return false; - if (swap) - swap_bytes <4> (X_CAST (char *, &di)); - dv(i) = di; - } + for (int i = 0; i < mdims; i++) + { + if (! is.read (X_CAST (char *, &di), 4)) + return false; + if (swap) + swap_bytes <4> (X_CAST (char *, &di)); + dv(i) = di; + } - char tmp; - if (! is.read (X_CAST (char *, &tmp), 1)) - return false; + char tmp; + if (! is.read (X_CAST (char *, &tmp), 1)) + return false; - Matrix m (dv(0), dv(1)); - double *re = m.fortran_vec (); - read_doubles (is, re, X_CAST (save_type, tmp), dv.numel (), swap, fmt); - if (error_state || ! is) - return false; + Matrix m (dv(0), dv(1)); + double *re = m.fortran_vec (); + read_doubles (is, re, X_CAST (save_type, tmp), dv.numel (), swap, fmt); + if (error_state || ! is) + return false; - gval = galois (m, mord, prim); - } + gval = galois (m, mord, prim); + } return true; } @@ -560,61 +561,61 @@ int32_t itmp; space_hid = H5Screate_simple (0, hdims, (hsize_t*) 0); - if (space_hid < 0) - { - H5Gclose (group_hid); - return false; - } + if (space_hid < 0) + { + H5Gclose (group_hid); + return false; + } #if HAVE_HDF5_18 - data_hid = H5Dcreate (group_hid, "m", H5T_NATIVE_UCHAR, space_hid, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + data_hid = H5Dcreate (group_hid, "m", H5T_NATIVE_UCHAR, space_hid, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); #else - data_hid = H5Dcreate (group_hid, "m", H5T_NATIVE_UCHAR, space_hid, - H5P_DEFAULT); + data_hid = H5Dcreate (group_hid, "m", H5T_NATIVE_UCHAR, space_hid, + H5P_DEFAULT); #endif - if (data_hid < 0) - { - H5Sclose (space_hid); - H5Gclose (group_hid); - return false; - } - + if (data_hid < 0) + { + H5Sclose (space_hid); + H5Gclose (group_hid); + return false; + } + tmp = m (); retval = H5Dwrite (data_hid, H5T_NATIVE_UCHAR, H5S_ALL, H5S_ALL, H5P_DEFAULT, - (void*) &tmp) >= 0; + (void*) &tmp) >= 0; H5Dclose (data_hid); if (!retval) - { - H5Sclose (space_hid); - H5Gclose (group_hid); - return false; - } + { + H5Sclose (space_hid); + H5Gclose (group_hid); + return false; + } #if HAVE_HDF5_18 - data_hid = H5Dcreate (group_hid, "prim", H5T_NATIVE_UINT, space_hid, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + data_hid = H5Dcreate (group_hid, "prim", H5T_NATIVE_UINT, space_hid, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); #else - data_hid = H5Dcreate (group_hid, "prim", H5T_NATIVE_UINT, space_hid, - H5P_DEFAULT); + data_hid = H5Dcreate (group_hid, "prim", H5T_NATIVE_UINT, space_hid, + H5P_DEFAULT); #endif - if (data_hid < 0) - { - H5Sclose (space_hid); - H5Gclose (group_hid); - return false; - } - + if (data_hid < 0) + { + H5Sclose (space_hid); + H5Gclose (group_hid); + return false; + } + itmp = primpoly (); retval = H5Dwrite (data_hid, H5T_NATIVE_UINT, H5S_ALL, H5S_ALL, H5P_DEFAULT, - (void*) &itmp) >= 0; + (void*) &itmp) >= 0; H5Dclose (data_hid); if (!retval) - { - H5Sclose (space_hid); - H5Gclose (group_hid); - return false; - } + { + H5Sclose (space_hid); + H5Gclose (group_hid); + return false; + } H5Sclose (space_hid); // Octave uses column-major, while HDF5 uses row-major ordering @@ -629,41 +630,41 @@ hid_t save_type_hid = H5T_NATIVE_UINT; if (tmp <= 8) - { - save_type_hid = H5T_NATIVE_UCHAR; - char *wtmp = (char *)vtmp; - for (int i = 0; i < d.numel (); i++) - wtmp[i] = (char) mtmp[i]; - } + { + save_type_hid = H5T_NATIVE_UCHAR; + char *wtmp = (char *)vtmp; + for (int i = 0; i < d.numel (); i++) + wtmp[i] = (char) mtmp[i]; + } else if (tmp <= 16) - { - save_type_hid = H5T_NATIVE_USHORT; - uint16_t *wtmp = (uint16_t *)vtmp; - for (int i = 0; i < d.numel (); i++) - wtmp[i] = (uint16_t) mtmp[i]; - } + { + save_type_hid = H5T_NATIVE_USHORT; + uint16_t *wtmp = (uint16_t *)vtmp; + for (int i = 0; i < d.numel (); i++) + wtmp[i] = (uint16_t) mtmp[i]; + } else - { - for (int i = 0; i < d.numel (); i++) - vtmp[i] = (uint32_t) mtmp[i]; - } + { + for (int i = 0; i < d.numel (); i++) + vtmp[i] = (uint32_t) mtmp[i]; + } #if HAVE_HDF5_18 - data_hid = H5Dcreate (group_hid, "val", save_type_hid, space_hid, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + data_hid = H5Dcreate (group_hid, "val", save_type_hid, space_hid, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); #else - data_hid = H5Dcreate (group_hid, "val", save_type_hid, space_hid, - H5P_DEFAULT); + data_hid = H5Dcreate (group_hid, "val", save_type_hid, space_hid, + H5P_DEFAULT); #endif if (data_hid < 0) - { - H5Sclose (space_hid); - H5Gclose (group_hid); - return false; - } + { + H5Sclose (space_hid); + H5Gclose (group_hid); + return false; + } retval = H5Dwrite (data_hid, save_type_hid, H5S_ALL, H5S_ALL, - H5P_DEFAULT, (void*) vtmp) >= 0; + H5P_DEFAULT, (void*) vtmp) >= 0; H5Dclose (data_hid); H5Sclose (space_hid); @@ -673,7 +674,7 @@ bool octave_galois::load_hdf5 (hid_t loc_id, const char *name, - bool have_h5giterate_bug) + bool have_h5giterate_bug) { bool retval = false; char mord; @@ -697,19 +698,19 @@ rank = H5Sget_simple_extent_ndims (space_id); if (rank != 0) - { - H5Dclose (data_hid); - H5Gclose (group_hid); - return false; - } + { + H5Dclose (data_hid); + H5Gclose (group_hid); + return false; + } - if (H5Dread (data_hid, H5T_NATIVE_UCHAR, H5S_ALL, H5S_ALL, - H5P_DEFAULT, (void *) &mord) < 0) - { - H5Dclose (data_hid); - H5Gclose (group_hid); - return false; - } + if (H5Dread (data_hid, H5T_NATIVE_UCHAR, H5S_ALL, H5S_ALL, + H5P_DEFAULT, (void *) &mord) < 0) + { + H5Dclose (data_hid); + H5Gclose (group_hid); + return false; + } H5Dclose (data_hid); #if HAVE_HDF5_18 @@ -721,19 +722,19 @@ rank = H5Sget_simple_extent_ndims (space_id); if (rank != 0) - { - H5Dclose (data_hid); - H5Gclose (group_hid); - return false; - } + { + H5Dclose (data_hid); + H5Gclose (group_hid); + return false; + } - if (H5Dread (data_hid, H5T_NATIVE_UINT, H5S_ALL, H5S_ALL, - H5P_DEFAULT, (void *) &prim) < 0) - { - H5Dclose (data_hid); - H5Gclose (group_hid); - return false; - } + if (H5Dread (data_hid, H5T_NATIVE_UINT, H5S_ALL, H5S_ALL, + H5P_DEFAULT, (void *) &prim) < 0) + { + H5Dclose (data_hid); + H5Gclose (group_hid); + return false; + } H5Dclose (data_hid); #if HAVE_HDF5_18 @@ -746,12 +747,12 @@ // Only handle matrices for now if (rank != 2) - { - H5Sclose (space_id); - H5Dclose (data_hid); - H5Gclose (group_hid); - return false; - } + { + H5Sclose (space_id); + H5Dclose (data_hid); + H5Gclose (group_hid); + return false; + } OCTAVE_LOCAL_BUFFER (hsize_t, hdims, rank); OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank); @@ -760,12 +761,12 @@ MArray<int> m (dim_vector (hdims[1], hdims[0])); int *re = m.fortran_vec (); - if (H5Dread (data_hid, H5T_NATIVE_UINT, H5S_ALL, H5S_ALL, - H5P_DEFAULT, (void *) re) >= 0) - { - retval = true; - gval = galois (m, mord, prim); - } + if (H5Dread (data_hid, H5T_NATIVE_UINT, H5S_ALL, H5S_ALL, + H5P_DEFAULT, (void *) re) >= 0) + { + retval = true; + gval = galois (m, mord, prim); + } H5Sclose (space_id); H5Dclose (data_hid); @@ -774,7 +775,7 @@ #endif /* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** + ;;; Local Variables: *** + ;;; mode: C++ *** + ;;; End: *** */
--- a/main/comm/src/ov-galois.h Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/ov-galois.h Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 @@ -49,10 +50,10 @@ 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) + 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) { } @@ -63,17 +64,17 @@ 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); + 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); } + 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); + bool resize_ok); octave_value do_index_op (const octave_value_list& idx) - { return do_index_op (idx, 0); } + { return do_index_op (idx, 0); } void assign (const octave_value_list& idx, const galois& rhs); @@ -109,7 +110,7 @@ void print_info (std::ostream& os, const std::string& prefix) const; bool is_real_matrix (void) const { return false; } - + bool is_real_type (void) const { return false; } // XXX FIXME XXX @@ -118,7 +119,7 @@ double double_value (bool = false) const; double scalar_value (bool frc_str_conv = false) const - { return double_value (frc_str_conv); } + { return double_value (frc_str_conv); } Matrix matrix_value (bool = false) const; @@ -127,9 +128,9 @@ Complex complex_value (bool = false) const; ComplexMatrix complex_matrix_value (bool = false) const - { return ComplexMatrix ( matrix_value()); } + { return ComplexMatrix ( matrix_value()); } - galois galois_value (void) const { return gval; } + galois galois_value (void) const { return gval; } octave_value_list dotref (const octave_value_list& idx); @@ -141,9 +142,9 @@ 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 load_binary (std::istream& is, bool swap, + oct_mach_info::float_format fmt); #if defined (HAVE_HDF5) bool save_hdf5 (hid_t loc_id, const char *name, bool save_as_floats);
--- a/main/comm/src/primpoly.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/primpoly.cc Fri May 11 18:08:42 2012 +0000 @@ -1,21 +1,22 @@ //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 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. +// 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/>. +// 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) +// 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 <iomanip> @@ -24,12 +25,12 @@ #include <octave/pager.h> enum primpoly_type -{ - PRIMPOLY_MIN=0, - PRIMPOLY_MAX, - PRIMPOLY_ALL, - PRIMPOLY_K -}; + { + PRIMPOLY_MIN=0, + PRIMPOLY_MAX, + PRIMPOLY_ALL, + PRIMPOLY_K + }; static bool do_isprimitive (const int& a, const int& m) @@ -123,16 +124,16 @@ std::string s_arg = args(1).string_value (); if (s_arg == "nodisplay") - disp = false; + disp = false; else if (s_arg == "min") - type = PRIMPOLY_MIN; + type = PRIMPOLY_MIN; else if (s_arg == "max") - type = PRIMPOLY_MAX; + type = PRIMPOLY_MAX; else if (s_arg == "all") - type = PRIMPOLY_ALL; + type = PRIMPOLY_ALL; else { - error ("primpoly: invalid argument"); - return retval; + error ("primpoly: invalid argument"); + return retval; } } else { error ("primpoly: incorrect argument type"); @@ -143,8 +144,8 @@ if (nargin > 2) { if (args(2).is_scalar_type ()) { if (type == PRIMPOLY_K) { - error ("primpoly: invalid arguments"); - return retval; + error ("primpoly: invalid arguments"); + return retval; } k = args(2).int_value(); type = PRIMPOLY_K; @@ -152,25 +153,25 @@ std::string s_arg = args(2).string_value (); if (s_arg == "nodisplay") { - if (!disp) { - error ("primpoly: invalid arguments"); - return retval; - } - disp = false; + if (!disp) { + error ("primpoly: invalid arguments"); + return retval; + } + disp = false; } else if (!disp) { - if (s_arg == "min") - type = PRIMPOLY_MIN; - else if (s_arg == "max") - type = PRIMPOLY_MAX; - else if (s_arg == "all") - type = PRIMPOLY_ALL; - else { - error ("primpoly: invalid argument"); - return retval; - } + if (s_arg == "min") + type = PRIMPOLY_MIN; + else if (s_arg == "max") + type = PRIMPOLY_MAX; + else if (s_arg == "all") + type = PRIMPOLY_ALL; + else { + error ("primpoly: invalid argument"); + return retval; + } } else { - error ("primpoly: invalid arguments"); - return retval; + error ("primpoly: invalid arguments"); + return retval; } } else { error ("primpoly: incorrect argument type"); @@ -183,36 +184,36 @@ primpolys.resize(1); for (int i = (1<<m)+1; i < (1<<(1+m)); i+=2) if (do_isprimitive(i,m)) { - primpolys(0) = (double)i; - break; + primpolys(0) = (double)i; + break; } break; case PRIMPOLY_MAX: primpolys.resize(1); for (int i = (1<<(m+1))-1; i > (1<<m); i-=2) if (do_isprimitive(i,m)) { - primpolys(0) = (double)i; - break; + primpolys(0) = (double)i; + break; } break; case PRIMPOLY_ALL: for (int i = (1<<m)+1; i < (1<<(1+m)); i+=2) if (do_isprimitive(i,m)) { - primpolys.resize(primpolys.length()+1); - primpolys(primpolys.length()-1) = (double)i; + primpolys.resize(primpolys.length()+1); + primpolys(primpolys.length()-1) = (double)i; } break; case PRIMPOLY_K: for (int i = (1<<m)+1; i < (1<<(1+m)); i+=2) { int ki = 0; for (int j=0; j < m+1; j++) - if (i & (1 << j)) - ki++; + if (i & (1 << j)) + ki++; if (ki == k) { - if (do_isprimitive(i,m)) { - primpolys.resize(primpolys.length()+1); - primpolys(primpolys.length()-1) = (double)i; - } + if (do_isprimitive(i,m)) { + primpolys.resize(primpolys.length()+1); + primpolys(primpolys.length()-1) = (double)i; + } } } break; @@ -227,27 +228,27 @@ else { octave_stdout << std::endl << "Primitive polynomial(s) =" << std::endl << std::endl; for (int i=0; i < primpolys.length(); i++) { - bool first = true; - for (int j=m; j>=0; j--) { - if ((int)primpolys(i) & (1<<j)) { - if (j > 0) { - if (first) { - first = false; - octave_stdout << "D"; - } else - octave_stdout << "+D"; - if (j != 1) - octave_stdout << "^" << j; - } else { - if (first) { - first = false; - octave_stdout << "1"; - } else - octave_stdout << "+1"; - } - } - } - octave_stdout << std::endl; + bool first = true; + for (int j=m; j>=0; j--) { + if ((int)primpolys(i) & (1<<j)) { + if (j > 0) { + if (first) { + first = false; + octave_stdout << "D"; + } else + octave_stdout << "+D"; + if (j != 1) + octave_stdout << "^" << j; + } else { + if (first) { + first = false; + octave_stdout << "1"; + } else + octave_stdout << "+1"; + } + } + } + octave_stdout << std::endl; } } octave_stdout << std::endl;
--- a/main/comm/src/syndtable.cc Fri May 11 17:59:49 2012 +0000 +++ b/main/comm/src/syndtable.cc Fri May 11 18:08:42 2012 +0000 @@ -43,9 +43,9 @@ int l = pos.rows(); pos.resize(dim_vector (l+new_pos.rows(), cols), 0); for (int j=0; j<new_pos.rows(); j++) { - for (int k=0; k<cols; k++) - pos(l+j,k) = new_pos(j,k); - pos(l+j,COL_MAJ(i)) += (1<<COL_MIN(i)); + for (int k=0; k<cols; k++) + pos(l+j,k) = new_pos(j,k); + pos(l+j,COL_MAJ(i)) += (1<<COL_MIN(i)); } } } @@ -127,11 +127,11 @@ // Now use the syndrome as the rows indices to put the error vectors // in place if (((unsigned int)syndrome < nrows) && !filled(syndrome)) { - filled(syndrome) = 1; - nfilled--; - for (int i = 0; i < n; i++) - table(syndrome,i) = ((errpos(j,COL_MAJ(i)) & - ((unsigned int)1 << COL_MIN(i))) != 0); + filled(syndrome) = 1; + nfilled--; + for (int i = 0; i < n; i++) + table(syndrome,i) = ((errpos(j,COL_MAJ(i)) & + ((unsigned int)1 << COL_MIN(i))) != 0); } }