Mercurial > forge
changeset 884:1eb91fa57b97 octave-forge
Add sqrt over galois field, code cleanup, some optimization and html
version of docs
| author | adb014 |
|---|---|
| date | Tue, 01 Apr 2003 19:37:15 +0000 |
| parents | 125a742f9169 |
| children | 1ea43ea20fcb |
| files | main/comm/INDEX main/comm/Makefile main/comm/awgn.png main/comm/comms.m main/comm/comms.txi main/comm/eyediagram.m main/comm/eyediagram.png main/comm/galois-def.h main/comm/galois.cc main/comm/galois.h main/comm/gf.cc main/comm/gftable.m main/comm/ov-galois.cc main/comm/rsdecof.m main/comm/rsencof.m main/comm/scatterplot.m main/comm/scatterplot.png |
| diffstat | 17 files changed, 448 insertions(+), 349 deletions(-) [+] |
line wrap: on
line diff
--- a/main/comm/INDEX Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/INDEX Tue Apr 01 19:37:15 2003 +0000 @@ -57,7 +57,7 @@ rcosiir rcosine rcosfir -Galois Fields +Galois Fields of Even Characateristic + - = Addition and subtraction in a Galois Field. * / \ = Matrix multiplication and division of Galois arrays. .* ./ .\ = Element by element multiplication and division of Galois arrays. @@ -86,6 +86,7 @@ glog glu gprod + gsqrt grank greshape groots
--- a/main/comm/Makefile Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/Makefile Tue Apr 01 19:37:15 2003 +0000 @@ -2,6 +2,7 @@ MAKEINFO = makeinfo --no-split TEXI2DVI = texi2dvi --clean +TEXI2HTML = texi2html -split_chapter DVIPS = dvips GALOISTARGET = gf.oct @@ -10,9 +11,9 @@ ov-galois.cc GALOISOBJECTS = $(patsubst %.cc,%.o,$(GALOISSOURCES)) GALOISLINKTARGETS = isgalois.oct gdiag.oct greshape.oct gprod.oct gsum.oct \ - gsumsq.oct glog.oct gexp.oct gfilter.oct glu.oct ginv.oct \ - ginverse.oct gdet.oct grank.oct rsenc.oct rsdec.oct \ - bchenco.oct bchdeco.oct + gsumsq.oct gsqrt.oct glog.oct gexp.oct gfilter.oct \ + glu.oct ginv.oct ginverse.oct gdet.oct grank.oct \ + rsenc.oct rsdec.oct bchenco.oct bchdeco.oct GALOISDEPENDS = $(patsubst %.cc,%.d,$(GALOISSOURCES)) OTHERSOURCES = primpoly.cc isprimitive.cc _errcore.cc cyclpoly.cc cyclgen.cc \ @@ -86,6 +87,10 @@ @echo "Making info $@"; \ $(MAKEINFO) $< +%.html : %.texi + @echo "Making html $@"; \ + $(TEXI2HTML) $< ; \ + %.texi : %.txi @echo "Making texinfo $@"; \ $(RM) -f $@; \
--- a/main/comm/comms.m Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/comms.m Tue Apr 01 19:37:15 2003 +0000 @@ -443,7 +443,7 @@ gmat = gf(reshape(mod([0:matlen*matlen-1],2^m),matlen,matlen),m); fprintf("PASSED\n"); fprintf(" Access Galois structures: "); - if (gcol.m ~= m || gcol.prim_poly ~= primpoly(m,"min", ... + if (gcol.m != m || gcol.prim_poly != primpoly(m,"min", ... "nodisplay")) error("FAILED"); endif @@ -462,39 +462,47 @@ error("FAILED"); endif tmp = gdiag(grow); - if (size(tmp,1) ~= 2^m || size(tmp,2) ~= 2^m) + if (size(tmp,1) != 2^m || size(tmp,2) != 2^m) error("FAILED"); endif for i=1:2^m for j=1:2^m if ((i == j) && (tmp(i,j) != grow(i))) error("FAILED"); - elseif ((i ~= j) && (tmp(i,j) ~= 0)) + elseif ((i != j) && (tmp(i,j) != 0)) error("FAILED"); endif end end tmp = gdiag(gmat); - if (length(tmp) ~= matlen) + if (length(tmp) != matlen) error("FAILED"); endif for i=1:matlen - if (gmat(i,i) ~= tmp(i)) + if (gmat(i,i) != tmp(i)) error("FAILED"); endif end tmp = greshape(gmat,prod(size(gmat)),1); - if (length(tmp) ~= prod(size(gmat))) + if (length(tmp) != prod(size(gmat))) + error("FAILED"); + endif + if (gexp(glog(gf([1:n],m))) != [1:n]) error("FAILED"); endif + tmp = gsqrt(gmat); + if (tmp .* tmp != gmat) + error("FAILED"); + endif + fprintf("PASSED\n"); fprintf(" Unary operators: "); tmp = - grow; - if (tmp ~= grow) + if (tmp != grow) error("FAILED"); endif tmp = !grow; - if (tmp(1) ~= 1) + if (tmp(1) != 1) error("FAILED"); endif if (any(tmp(2:length(tmp)))) @@ -503,7 +511,7 @@ tmp = gmat'; for i=1:size(gmat,1) for j=1:size(gmat,2) - if (gmat(i,j) ~= tmp(j,i)) + if (gmat(i,j) != tmp(j,i)) error("FAILED"); endif end @@ -517,25 +525,25 @@ elsqr = gcol .* gcol; elsqr2 = gcol .^ gf(2,m); for i=1:length(elsqr) - if (elsqr(i) ~= multbl(i,i)) + if (elsqr(i) != multbl(i,i)) error("FAILED"); endif end for i=1:length(elsqr) - if (elsqr(i) ~= elsqr2(i)) + if (elsqr(i) != elsqr2(i)) error("FAILED"); endif end tmp = grow(2:length(grow)) ./ gcol(2:length(gcol))'; - if (length(tmp) ~= n || any(tmp ~= ones(1,length(grow)-1))) + if (length(tmp) != n || any(tmp != ones(1,length(grow)-1))) error("FAILED"); endif fprintf("PASSED\n"); fprintf(" Logical operators: "); - if (grow(1) ~= gzero || grow(2) ~= gone || grow(2^m) ~= n) + if (grow(1) != gzero || grow(2) != gone || grow(2^m) != n) error("FAILED"); endif - if (!(grow(1) == gzero) || any(grow ~= gcol')) + if (!(grow(1) == gzero) || any(grow != gcol')) error("FAILED"); endif fprintf("PASSED\n");
--- a/main/comm/comms.txi Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/comms.txi Tue Apr 01 19:37:15 2003 +0000 @@ -298,7 +298,8 @@ @iftex which produces a sine-wave with noise added as seen in Figure 1. -@center @image{awgn} +@center @image{awgn, , } + @center Figure 1: Sine-wave with 10dB signal-to-noise ratio @end iftex @@ -426,7 +427,8 @@ which produces a eye-diagram of a noisy signal as seen in Figure 2. Similarly an example of the use of the function @dfn{scatterplot} is -@center @image{eyediagram} +@center @image{eyediagram, , } + @center Figure 2: Eye-diagram of a QPSK like signal with 15dB signal-to-noise ratio @end iftex @ifinfo @@ -450,13 +452,24 @@ @iftex which produces a scatterplot of a noisy signal as seen in Figure 3. -@center @image{scatterplot} +@center @image{scatterplot, , } + @center Figure 3: Scatterplot of a QPSK like signal with 15dB signal-to-noise ratio @end iftex @node Source Coding, Block Coding, Random Signals, Top @chapter Source Coding +@menu +* Quantization:: +* PCM Coding:: +* Arithmetic Coding:: +* Dynamic Range Compression:: +@end menu + +@node Quantization, PCM Coding, , Source Coding +@section Quantization + An important aspect of converting an analog signal to the digital domain is quantization. This is the process of mapping a continuous signal to a set of defined values. Octave contains two functions to perform quantization, @@ -510,11 +523,23 @@ signal to be quantized, then the parititioning suggested by @dfn{lloyds} will be sub-optimal. +@node PCM Coding, Arithmetic Coding, Quantization, Source Coding +@section PCM Coding + The DPCM function @dfn{dpcmenco}, @dfn{dpcmdeco} and @dfn{dpcmopt} implement a form of preditive quantization, where the predictability of the signal is used to further compress it. These functions are not yet implemented. +@node Arithmetic Coding, Dynamic Range Compression, PCM Coding, Source Coding +@section Arithmetic Coding + +The arithmetic coding functions @dfn{arithenco} and @dfn{arithdeco} are +not yet implemented. + +@node Dynamic Range Compression, , Arithmetic Coding, Source Coding +@section Dynamic Range Compression + The final source coding function is @dfn{compand} which is used to compress and expand the dynamic range of a signal. For instance consider a logarithm quantized by a linear partitioning. Such a
--- a/main/comm/eyediagram.m Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/eyediagram.m Tue Apr 01 19:37:15 2003 +0000 @@ -15,11 +15,11 @@ ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ## -*- texinfo -*- -## @deftypefn {Function File} eyediagram (@var{x},@var{n}) -## @deftypefnx {Function File} eyediagram (@var{x},@var{n},@var{per}) -## @deftypefnx {Function File} eyediagram (@var{x},@var{n},@var{per},@var{off}) -## @deftypefnx {Function File} eyediagram (@var{x},@var{n},@var{per},@var{off},@var{str}) -## @deftypefnx {Function File} eyediagram (@var{x},@var{n},@var{per},@var{off},@var{str},@var{h}) +## @deftypefn {Function File} {} eyediagram (@var{x},@var{n}) +## @deftypefnx {Function File} {} eyediagram (@var{x},@var{n},@var{per}) +## @deftypefnx {Function File} {} eyediagram (@var{x},@var{n},@var{per},@var{off}) +## @deftypefnx {Function File} {} eyediagram (@var{x},@var{n},@var{per},@var{off},@var{str}) +## @deftypefnx {Function File} {} eyediagram (@var{x},@var{n},@var{per},@var{off},@var{str},@var{h}) ## @deftypefnx {Function File} {@var{h} =} eyediagram (@var{...}) ## ## Plot the eye-diagram of a signal. The signal @var{x} can be either in one
--- a/main/comm/galois-def.h Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/galois-def.h Tue Apr 01 19:37:15 2003 +0000 @@ -79,12 +79,12 @@ for (int i=0; i<nr; i++) \ for (int j=0; j<nc; j++) \ { \ - if ((M1.elem(i,j) < 0) || (M1. elem(i,j) > NN)) \ + if ((M1(i,j) < 0) || (M1(i,j) > NN)) \ { \ gripe_nonconformant_galois (OP, M1.m()); \ return RET (); \ } \ - if (((double)M1.elem(i,j) - (double)((int)M1.elem(i,j))) != 0.) \ + if (((double)M1(i,j) - (double)((int)M1(i,j))) != 0.) \ { \ gripe_nonconformant_galois (OP, M1.m()); \ return RET (); \ @@ -100,7 +100,7 @@ for (int i=0; i<nr; i++) \ for (int j=0; j<nc; j++) \ { \ - if (M.elem(i,j) == 0) \ + if (M(i,j) == 0) \ { \ gripe_divzero_galois(OP); \ return RET (); \ @@ -183,6 +183,7 @@ } \ 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++) \ { \ @@ -190,7 +191,7 @@ r(i,j) = 0; \ else \ { \ - r(i,j) = r.index_of((int)m1(0,0)) OP r.index_of((int)m2(i,j)) + NN; \ + 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)); \ } \ @@ -208,6 +209,7 @@ } \ 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++) \ { \ @@ -215,7 +217,7 @@ r(i,j) = 0; \ else \ { \ - r(i,j) = r.index_of((int)m1(i,j)) OP r.index_of((int)m2(0,0)) + NN; \ + 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)); \ } \ @@ -276,7 +278,7 @@ \ for (int j = 0; j < m1_nc; j++) \ for (int i = 0; i < m1_nr; i++) \ - r.elem(i, j) = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \ + r(i, j) = C1 (m1(i, j)) OP C2 (m2(i, j)); \ } \ } \ else \ @@ -288,14 +290,14 @@ r.resize(m2_nr,m2_nc); \ for (int i=0; i<m2_nr; i++) \ for (int j=0; j<m2_nc; j++) \ - r.elem(i, j) = C1 (m1.elem(0, 0)) OP C2 (m2.elem(i, 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.elem(i, j) = C1 (m1.elem(i, j)) OP C2 (m2.elem(0, 0)); \ + r(i, j) = C1 (m1(i, j)) OP C2 (m2(0, 0)); \ } \ else \ gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ @@ -335,8 +337,8 @@ for (int j = 0; j < m1_nc; j++) \ for (int i = 0; i < m1_nr; i++) \ { \ - r.elem(i, j) = (m1.elem(i, j) != ZERO) \ - OP (m2.elem(i, j) != ZERO); \ + r(i, j) = (m1(i, j) != ZERO) \ + OP (m2(i, j) != ZERO); \ } \ } \ } \ @@ -347,16 +349,16 @@ r.resize(m2_nr,m2_nc); \ for (int i=0; i<m2_nr; i++) \ for (int j=0; j<m2_nc; j++) \ - r.elem(i, j) = (m1.elem(0, 0) != ZERO) \ - OP (m2.elem(i, j) != ZERO); \ + 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.elem(i, j) = (m1.elem(i, j) != ZERO) \ - OP (m2.elem(0, 0) != ZERO); \ + 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); \ @@ -382,7 +384,7 @@ RET.resize (nr, 1); \ for (int i = 0; i < nr; i++) \ { \ - RET.elem (i, 0) = INIT_VAL; \ + RET (i, 0) = INIT_VAL; \ for (int j = 0; j < nc; j++) \ { \ ROW_EXPR; \ @@ -394,7 +396,7 @@ RET.resize (1, nc); \ for (int j = 0; j < nc; j++) \ { \ - RET.elem (0, j) = INIT_VAL; \ + RET (0, j) = INIT_VAL; \ for (int i = 0; i < nr; i++) \ { \ COL_EXPR; \
--- a/main/comm/galois.cc Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/galois.cc Tue Apr 01 19:37:15 2003 +0000 @@ -50,7 +50,7 @@ gripe_integer_galois(); return; } - elem(i,j) = (int)a(i,j); + xelem(i,j) = (int)a(i,j); } } @@ -71,7 +71,7 @@ gripe_integer_galois(); return; } - elem(i,j) = (int)a(i,j); + xelem(i,j) = (int)a(i,j); } } @@ -89,7 +89,7 @@ for (int i=0; i<rows(); i++) for (int j=0; j<columns(); j++) - elem(i,j) = val; + xelem(i,j) = val; field = stored_galois_fields.create_galois_field(_m, _primpoly); } @@ -110,7 +110,7 @@ for (int i=0; i<rows(); i++) for (int j=0; j<columns(); j++) - elem(i,j) = (int)val; + xelem(i,j) = (int)val; field = stored_galois_fields.create_galois_field(_m, _primpoly); } @@ -184,8 +184,9 @@ return *this; } - for (int i = 0; i < a.length (); i++) - elem (i, i) ^= a.elem (i, i); + for (int i=0; i<rows(); i++) + for (int j=0; j<columns(); j++) + xelem(i,j) ^= a (i, j); return *this; } @@ -215,8 +216,9 @@ return *this; } - for (int i = 0; i < a.length (); i++) - elem (i, i) ^= a.elem (i, i); + for (int i=0; i<rows(); i++) + for (int j=0; j<columns(); j++) + xelem(i,j) ^= a (i, j); return *this; } @@ -260,17 +262,17 @@ if (k > 0) { for (int i = 0; i < ndiag; i++) - retval.elem (i,0) = elem (i, i+k); + retval (i,0) = xelem (i, i+k); } else if ( k < 0) { for (int i = 0; i < ndiag; i++) - retval.elem (i,0) = elem (i-k, i); + retval (i,0) = xelem (i-k, i); } else { for (int i = 0; i < ndiag; i++) - retval.elem (i,0) = elem (i, i); + retval (i,0) = xelem (i, i); } } else @@ -291,7 +293,7 @@ for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) - b.elem (i, j) = ! elem (i, j); + b (i, j) = ! xelem (i, j); return b; } @@ -304,7 +306,7 @@ a.resize(d2,d1); for (int j = 0; j < d2; j++) for (int i = 0; i < d1; i++) - a.xelem (j, i) = xelem (i, j); + a (j, i) = xelem (i, j); return a; } @@ -320,9 +322,9 @@ galois elem_pow (const galois& a, const galois& b) { - galois result(a); int a_nr = a.rows (); int a_nc = a.cols (); + galois result(a_nr, a_nc, 0, a.m(), a.primpoly()); int b_nr = b.rows (); int b_nc = b.cols (); @@ -339,28 +341,24 @@ if (a_nr == 1 && a_nc == 1) { - result.resize(b_nr,b_nc); + result.resize(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.elem(i,j) == 0) + if (b(i,j) == 0) result (i,j) = 1; - else if (a.elem(0,0) == 0) - result (i,j) = 0; - else - result (i,j) = a.alpha_to(modn(a.index_of(a.elem(0,0)) * - b.elem(i,j), a.m(),a.n())); + 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.elem(0,0) == 0) + if (b(0,0) == 0) result (i,j) = 1; - else if (a.elem(i,j) == 0) - result (i,j) = 0; - else - result (i,j) = a.alpha_to(modn(a.index_of(a.elem(i,j)) * - b.elem(0,0),a.m(),a.n())); + 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 { @@ -372,13 +370,11 @@ for (int j = 0; j < a_nc; j++) for (int i = 0; i < a_nr; i++) - if (b.elem(i,j) == 0) + if (b(i,j) == 0) result (i,j) = 1; - else if (a.elem(i,j) == 0) - result (i,j) = 0; - else - result (i,j) = a.alpha_to(modn(a.index_of(a.elem(i,j)) * - b.elem(i,j),a.m(),a.n())); + 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; @@ -386,15 +382,15 @@ galois elem_pow (const galois& a, const Matrix& b) { - galois result(a); int a_nr = a.rows (); int a_nc = a.cols (); + galois result(a_nr, a_nc, 0, a.m(), a.primpoly()); int b_nr = b.rows (); int b_nc = b.cols (); if (b_nr == 1 && b_nc == 1) - return elem_pow(a, b.elem(0,0)); + return elem_pow(a, b(0,0)); if (a_nr != b_nr || a_nc != b_nc) { @@ -405,15 +401,13 @@ for (int j = 0; j < a_nc; j++) for (int i = 0; i < a_nr; i++) { - int tmp = (int)b.elem(i,j); + int tmp = (int)b(i,j); while (tmp < 0) tmp += a.n(); if (tmp == 0) result (i,j) = 1; - else if (a.elem(i,j) == 0) - result (i,j) = 0; - else - result (i,j) = a.alpha_to(modn(a.index_of(a.elem(i,j)) * tmp, + 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; @@ -421,9 +415,9 @@ galois elem_pow (const galois& a, double b) { - galois result(a); int a_nr = a.rows (); int a_nc = a.cols (); + galois result(a_nr, a_nc, 0, a.m(), a.primpoly()); int bi = (int) b; if ((double)bi != b) { @@ -439,10 +433,8 @@ { if (bi == 0) result (i,j) = 1; - else if (a.elem(i,j) == 0) - result (i,j) = 0; - else - result (i,j) = a.alpha_to(modn(a.index_of(a.elem(i,j)) * + 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; @@ -450,9 +442,9 @@ galois elem_pow (const galois &a, int b) { - galois result(a); int a_nr = a.rows (); int a_nc = a.cols (); + galois result(a_nr, a_nc, 0, a.m(), a.primpoly()); while (b < 0) b += a.n(); @@ -462,10 +454,8 @@ { if (b == 0) result (i,j) = 1; - else if (a.elem(i,j) == 0) - result (i,j) = 0; - else - result (i,j) = a.alpha_to(modn(a.index_of(a.elem(i,j)) * b, + 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; @@ -501,7 +491,7 @@ gripe_square_galois(); return galois (); } else - return pow (a, b.elem(0,0)); + return pow (a, b(0,0)); } galois pow (const galois& a, int b) @@ -521,7 +511,7 @@ { retval = galois(nr, nc, 0, a.m(), a.primpoly()); for (int i =0; i<nr; i++) - retval.elem(i,i) = 1; + retval(i,i) = 1; } else { @@ -600,17 +590,23 @@ galois retval(a_nr, b_nc, 0, a.m(), a.primpoly()); if (a_nr != 0 && a_nc != 0 && b_nc != 0) { - for (int i=0; i<b_nc; i++) - for (int j=0; j<b_nr; j++) - if (b.elem(j,i) != 0) { - int temp = b.elem(j, i); + // 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.elem(k,j) != 0) + if (a(k,j) != 0) retval(k,i) = retval(k,i) ^ a.alpha_to( - modn(a.index_of(temp) + - a.index_of(a.elem(k,j)),a.m(),a.n())); + modn(tmp + c(k,0),a.m(),a.n())); } } + } } return retval; } @@ -640,17 +636,17 @@ galois retval (0, 0, 0, m(), primpoly()); #define ROW_EXPR \ - if ((retval.elem (i, 0) == 0) || (elem(i,j) == 0)) \ - retval.elem (i, 0) = 0; \ + if ((retval (i, 0) == 0) || (elem(i,j) == 0)) \ + retval (i, 0) = 0; \ else \ - retval.elem (i, 0) = alpha_to(modn(index_of(retval.elem(i,0)) + \ + retval (i, 0) = alpha_to(modn(index_of(retval(i,0)) + \ index_of(elem(i,j)),m(),n())); #define COL_EXPR \ - if ((retval.elem (0, j) == 0) || (elem(i,j) == 0)) \ - retval.elem (0, j) = 0; \ + if ((retval (0, j) == 0) || (elem(i,j) == 0)) \ + retval (0, j) = 0; \ else \ - retval.elem (0, j) = alpha_to(modn(index_of(retval.elem(0,j)) + \ + retval (0, j) = alpha_to(modn(index_of(retval(0,j)) + \ index_of(elem(i,j)),m(),n())); GALOIS_REDUCTION_OP(retval, ROW_EXPR, COL_EXPR, 1, 1); @@ -672,10 +668,10 @@ #define ROW_EXPR \ - retval.elem (i, 0) ^= elem(i,j); + retval (i, 0) ^= elem(i,j); #define COL_EXPR \ - retval.elem (0, j) ^= elem(i,j); + retval (0, j) ^= elem(i,j); GALOIS_REDUCTION_OP (retval, ROW_EXPR, COL_EXPR, 0, 0); return retval; @@ -696,11 +692,11 @@ #define ROW_EXPR \ if (elem(i,j) != 0) \ - retval.elem (i, 0) ^= alpha_to(modn(2*index_of(elem(i,j)),m(),n())); + retval (i, 0) ^= alpha_to(modn(2*index_of(elem(i,j)),m(),n())); #define COL_EXPR \ if (elem(i,j) != 0) \ - retval.elem (0, j) ^= alpha_to(modn(2*index_of(elem(i,j)),m(),n())); + retval (0, j) ^= alpha_to(modn(2*index_of(elem(i,j)),m(),n())); GALOIS_REDUCTION_OP (retval, ROW_EXPR, COL_EXPR, 0, 0); return retval; @@ -710,6 +706,25 @@ } galois +galois::sqrt (void) const +{ + galois retval (*this); + int nr = rows (); + int nc = cols (); + + for (int j=0; j<nc; j++) { + for (int i=0; i<nr; i++) + if (retval.index_of(retval(i,j)) & 1) + retval(i,j) = retval.alpha_to((retval.index_of(retval(i,j)) + + retval.n()) / 2); + else + retval(i,j) = retval.alpha_to(retval.index_of(retval(i,j)) + / 2); + } + return retval; +} + +galois galois::log (void) const { bool warned = false; @@ -722,9 +737,9 @@ int nr = rows (); int nc = cols (); - for (int i=0; i<nr; i++) - for (int j=0; j<nc; j++) { - if (retval.elem(i,j) == 0) { + 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; @@ -732,9 +747,9 @@ // 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.elem(i,j) = retval.index_of(retval.elem(i,j)); + retval(i,j) = retval.index_of(retval(i,j)); } else - retval.elem(i,j) = retval.index_of(retval.elem(i,j)); + retval(i,j) = retval.index_of(retval(i,j)); } return retval; } @@ -752,9 +767,9 @@ int nr = rows (); int nc = cols (); - for (int i=0; i<nr; i++) - for (int j=0; j<nc; j++) { - if (retval.elem(i,j) == n()) { + 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; @@ -762,9 +777,9 @@ // 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.elem(i,j) = retval.alpha_to(retval.elem(i,j)); + retval(i,j) = retval.alpha_to(retval(i,j)); } else - retval.elem(i,j) = retval.alpha_to(retval.elem(i,j)); + retval(i,j) = retval.alpha_to(retval(i,j)); } return retval; } @@ -787,48 +802,47 @@ for (int j = 0; j < mn; j++) { int jp = j; - // Find the pivot and test for singularity if (ptype == LU::ROW) { for (int i = j+1; i < a_nr; i++) - if (a_fact.elem(i,j) > a_fact.elem(jp,j)) + if (a_fact(i,j) > a_fact(jp,j)) jp = i; } else { for (int i = j+1; i < a_nc; i++) - if (a_fact.elem(j,i) > a_fact.elem(j,jp)) + if (a_fact(j,i) > a_fact(j,jp)) jp = i; } - ipvt.elem(j) = jp; + ipvt(j) = jp; - if (a_fact.elem(jp,j) != 0) { + if (a_fact(jp,j) != 0) { if (ptype == LU::ROW) { // Apply the interchange to columns 1:NC. if (jp != j) for (int i = 0; i < a_nc; i++) { - int tmp = a_fact.elem(j,i); - a_fact.elem(j,i) = a_fact.elem(jp,i); - a_fact.elem(jp,i) = tmp; + 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.elem(i,j); - a_fact.elem(i,j) = a_fact.elem(i,jp); - a_fact.elem(i,jp) = tmp; + 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.elem(j,j)); + int idxj = a_fact.index_of(a_fact(j,j)); for (int i = j+1; i < a_nr; i++) { - if (a_fact.elem(i,j) == 0) - a_fact.elem(i,j) = 0; + if (a_fact(i,j) == 0) + a_fact(i,j) = 0; else - a_fact.elem(i,j) = a_fact.alpha_to(modn(a_fact.index_of( - a_fact.elem(i,j)) - idxj + a_fact.n(), a_fact.m(), + 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())); } } @@ -839,12 +853,12 @@ if (j < mn-1) { // Update trailing submatrix. for (int i=j+1; i < a_nr; i++) { - if (a_fact.elem(i,j) != 0) { - int idxi = a_fact.index_of(a_fact.elem(i,j)); + 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.elem(j,k) != 0) - a_fact.elem(i,k) ^= a_fact.alpha_to(modn(a_fact.index_of( - a_fact.elem(j,k)) + idxi, a_fact.m(), a_fact.n())); + 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())); } } } @@ -861,12 +875,12 @@ galois l (a_nr, mn, 0, a_fact.m(), a_fact.primpoly()); - for (int i = 0; i < a_nr; i++) - { - l.xelem (i, i) = 1; - for (int j = 0; j < (i > a_nc ? a_nc : i); j++) - l.xelem (i, j) = a_fact.xelem (i, j); - } + for (int i = 0; i < mn; i++) + l (i, i) = 1; + + for (int j = 0; j < mn; j++) + for (int i = j+1; i < a_nr; i++) + l (i, j) = a_fact (i, j); return l; } @@ -880,12 +894,9 @@ galois u (mn, a_nc, 0, a_fact.m(), a_fact.primpoly()); - for (int i = 0; i < mn; i++) - { - for (int j = i; j < a_nc; j++) - u.xelem (i, j) = a_fact.xelem (i, j); - } - + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < (j+1 > mn ? mn : j+1); i++) + u (i, j) = a_fact (i, j); return u; } @@ -899,24 +910,24 @@ Array<int> pvt (n); for (int i = 0; i < n; i++) - pvt.xelem (i) = i; + pvt (i) = i; for (int i = 0; i < ipvt.length(); i++) { - int k = ipvt.xelem (i); + int k = ipvt (i); if (k != i) { - int tmp = pvt.xelem (k); - pvt.xelem (k) = pvt.xelem (i); - pvt.xelem (i) = tmp; + int tmp = pvt (k); + pvt (k) = pvt (i); + pvt (i) = tmp; } } Matrix p(n, n, 0.0); for (int i = 0; i < n; i++) - p.xelem (i, pvt.xelem (i)) = 1.0; + p (i, pvt (i)) = 1.0; return p; } @@ -956,7 +967,7 @@ // Solve with identity matrix to find the inverse. galois btmp(nr, nr, 0, m(), primpoly()); for (int i=0; i < nr; i++) - btmp.elem(i,i) = 1; + btmp(i,i) = 1; galois retval = solve(btmp, info, 0); @@ -985,7 +996,7 @@ if (nr == 0 || nc == 0) { info = 0; - retval.elem(0,0) = 1; + retval(0,0) = 1; } else { LU fact (*this); @@ -993,14 +1004,14 @@ galois A (fact.A()); info = 0; - retval.elem(0,0) = A.elem(0,0); + retval(0,0) = A(0,0); for (int i=1; i<nr; i++) { - if ((retval.elem (0, 0) == 0) || (A.elem(i,i) == 0)) { - retval.elem (0,0) = 0; + if ((retval (0, 0) == 0) || (A(i,i) == 0)) { + retval (0,0) = 0; error("What the hell are we doing here!!!"); } else - retval.elem (0,0) = alpha_to(modn(index_of(retval.elem(0,0)) + \ - index_of(A.elem(i,i)),m(),n())); + retval (0,0) = alpha_to(modn(index_of(retval(0,0)) + \ + index_of(A(i,i)),m(),n())); } } } @@ -1039,6 +1050,7 @@ int nc = cols (); int b_nr = b.rows (); int b_nc = b.cols (); + galois c(nr,1,0,m(),primpoly()); // if (nr == 0 || nc == 0 || nr != nc || nr != b_nr) { if (nr == 0 || nc == 0 || nr != b_nr) { @@ -1066,37 +1078,45 @@ retval.resize(b_nr+nc-nr,b_nc,0); //Solve L*X = B, overwriting B with X. - for (int j=0; j<b_nc; j++) - for (int k=0; k<(nc < nr ? nc : nr); k++) - if (retval.elem(k,j) != 0) { - int idx = index_of(retval.elem(k,j)); + 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.elem(i,k) != 0) - retval.elem(i,j) ^= alpha_to(modn(index_of(A.elem(i,k)) + idx, - m(), n())); + 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 j=0; j<b_nc; j++) - for (int k=(nc < nr ? nc-1 : nr-1); k>=0; k--) - if (retval.elem(k,j) != 0) { - retval.elem(k,j) = alpha_to(modn(index_of(retval.elem(k,j)) - - index_of(A.elem(k,k)) + n(), m(), n())); - int idx = index_of(retval.elem(k,j)); - for (int i=0; i<(k < nr ? k : nr); i++) - if (A.elem(i,k) != 0) - retval.elem(i,j) ^= alpha_to(modn(index_of(A.elem(i,k)) + idx, - m(), n())); + for (int k=(nc < nr ? nc-1 : nr-1); k>=0; k--) { + int mn = k+1 < nr ? k+1 : nr; + for (int i=0; i<mn; i++) + c(i,0) = index_of(A(i,k)); + mn = k < nr ? k : nr; + for (int j=0; j<b_nc; j++) + if (retval(k,j) != 0) { + retval(k,j) = alpha_to(modn(index_of(retval(k,j)) - + c(k,0) + n(), m(), n())); + int idx = index_of(retval(k,j)); + for (int i=0; i<mn; i++) + if (A(i,k) != 0) + retval(i,j) ^= alpha_to(modn(c(i,0) + idx, m(), n())); } + } // Apply row interchanges to the right hand sides. //for (int j=0; j<IP.length(); j++) { for (int j=IP.length()-1; j>=0; j--) { int piv = IP(j); for (int i=0; i<b_nc; i++) { - int tmp = retval.elem(j,i); - retval.elem(j,i) = retval.elem(piv,i); - retval.elem(piv,i) = tmp; + int tmp = retval(j,i); + retval(j,i) = retval(piv,i); + retval(piv,i) = tmp; } } } @@ -1120,35 +1140,42 @@ for (int j=0; j<IP.length(); j++) { int piv = IP(j); for (int i=0; i<b_nc; i++) { - int tmp = retval.elem(j,i); - retval.elem(j,i) = retval.elem(piv,i); - retval.elem(piv,i) = tmp; + int tmp = retval(j,i); + retval(j,i) = retval(piv,i); + retval(piv,i) = tmp; } } //Solve L*X = B, overwriting B with X. - for (int j=0; j<b_nc; j++) - for (int k=0; k<(nc < nr ? nc : nr); k++) - if (retval.elem(k,j) != 0) { - int idx = index_of(retval.elem(k,j)); + 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.elem(i,k) != 0) - retval.elem(i,j) ^= alpha_to(modn(index_of(A.elem(i,k)) + idx, - m(), n())); + 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 j=0; j<b_nc; j++) - for (int k=(nc < nr ? nc-1 : nr-1); k>=0; k--) - if (retval.elem(k,j) != 0) { - retval.elem(k,j) = alpha_to(modn(index_of(retval.elem(k,j)) - - index_of(A.elem(k,k)) + n(), m(), n())); - int idx = index_of(retval.elem(k,j)); - for (int i=0; i<(k < nr ? k : nr); i++) - if (A.elem(i,k) != 0) - retval.elem(i,j) ^= alpha_to(modn(index_of(A.elem(i,k)) + idx, - m(), n())); + for (int k=(nc < nr ? nc-1 : nr-1); k>=0; k--) { + int mn = k+1 < nr ? k+1 : nr; + for (int i=0; i<mn; i++) + c(i,0) = index_of(A(i,k)); + mn = k < nr ? k : nr; + for (int j=0; j<b_nc; j++) + if (retval(k,j) != 0) { + retval(k,j) = alpha_to(modn(index_of(retval(k,j)) - + c(k,0) + n(), m(), n())); + int idx = index_of(retval(k,j)); + for (int i=0; i<mn; i++) + if (A(i,k) != 0) + retval(i,j) ^= alpha_to(modn(c(i,0) + idx, m(), n())); } + } // Resize the number of solution rows if needed if (nc < nr)
--- a/main/comm/galois.h Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/galois.h Tue Apr 01 19:37:15 2003 +0000 @@ -99,6 +99,7 @@ galois prod (int dim) const; galois sum (int dim) const; galois sumsq (int dim) const; + galois sqrt (void) const; galois log (void) const; galois exp (void) const; galois inverse (void) const;
--- a/main/comm/gf.cc Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/gf.cc Tue Apr 01 19:37:15 2003 +0000 @@ -166,13 +166,13 @@ int n = nc + k; galois r (n, n, 0, m.m(), m.primpoly()); for (int i = 0; i < nc; i++) - r.elem (i+roff, i+coff) = m.elem (0, 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.elem (i+roff, i+coff) = m.elem (i, 0); + r (i+roff, i+coff) = m (i, 0); retval = new octave_galois (r); } } else { @@ -297,11 +297,11 @@ RowVector tmp1(mr*mc); for (int i=0;i<nr;i++) for (int j=0;j<nc;j++) - tmp1(i+j*nr) = (double)a.elem(i,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.elem(i,j) = (int)tmp1(i+j*mr); + tmp2(i,j) = (int)tmp1(i+j*mr); retval = new octave_galois(tmp2); } } @@ -389,6 +389,35 @@ DATA_REDUCTION (sumsq); } +// PKG_ADD: dispatch ("sqrt", "gsqrt", "galois"); +DEFUN_DLD (gsqrt, args, , + "-*- texinfo -*-\n" +"@deftypefn {Loadable Function} {} gsqrt (@var{x})\n" +"Compute the square root of @var{x}, element by element, in a Galois Field.\n" +"@end deftypefn\n" +"@seealso{exp}") +{ + octave_value retval; + int nargin = args.length (); + + if (!galois_type_loaded || (args(0).type_id () != + octave_galois::static_type_id ())) { + gripe_wrong_type_arg ("gsqrt", args(0)); + return retval; + } + + if (nargin != 1) { + print_usage("gsqrt"); + return retval; + } + + galois a = ((const octave_galois&) args(0).get_rep()).galois_value (); + + retval = new octave_galois(a.sqrt()); + + return retval; +} + // PKG_ADD: dispatch ("log", "glog", "galois"); DEFUN_DLD (glog, args, , "-*- texinfo -*-\n" @@ -419,6 +448,7 @@ return retval; } + // PKG_ADD: dispatch ("exp", "gexp", "galois"); DEFUN_DLD (gexp, args, , "-*- texinfo -*-\n" @@ -438,7 +468,7 @@ } if (nargin != 1) { - print_usage("exp"); + print_usage("gexp"); return retval; } @@ -462,7 +492,7 @@ int ab_len = (a.length() > b.length() ? a.length() : b.length()); b.resize(ab_len, 1, 0); galois retval(x.length(), 1, 0, b.m(), b.primpoly()); - int norm = a.elem(0,0); + int norm = a(0,0); if (norm == 0) { error("gfilter: the first element of a must be non-zero"); @@ -475,8 +505,8 @@ if (norm != 1) { int idx_norm = b.index_of(norm); for (int i=0; i < b.length(); i++) { - if (b.elem(i,0) != 0) - b.elem(i,0) = b.alpha_to(modn(b.index_of(b.elem(i,0))-idx_norm+b.n(), + 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())); } } @@ -486,74 +516,74 @@ if (norm != 1) { int idx_norm = a.index_of(norm); for (int i=0; i < a.length(); i++) - if (a.elem(i,0) != 0) - a.elem(i,0) = a.alpha_to(modn(a.index_of(a.elem(i,0))-idx_norm+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.elem(i,0) = si.elem(0,0); - if ((b.elem(0,0) != 0) && (x.elem(i,0) != 0)) - retval.elem(i,0) ^= b.alpha_to(modn(b.index_of(b.elem(0,0)) + - b.index_of(x.elem(i,0)),b.m(),b.n())); + retval(i,0) = si(0,0); + if ((b(0,0) != 0) && (x(i,0) != 0)) + retval(i,0) ^= b.alpha_to(modn(b.index_of(b(0,0)) + + b.index_of(x(i,0)),b.m(),b.n())); if (si.length() > 1) { for (int j = 0; j < si.length() - 1; j++) { - si.elem(j,0) = si.elem(j+1,0); - if ((a.elem(j+1,0) != 0) && (retval.elem(i,0) != 0)) - si.elem(j,0) ^= a.alpha_to(modn(a.index_of(a.elem(j+1,0)) + - a.index_of(retval.elem(i,0)),a.m(),a.n())); - if ((b.elem(j+1,0) != 0) && (x.elem(i,0) != 0)) - si.elem(j,0) ^= b.alpha_to(modn(b.index_of(b.elem(j+1,0)) + - b.index_of(x.elem(i,0)),b.m(),b.n())); + 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.elem(si.length()-1,0) = 0; - if ((a.elem(si.length(),0) != 0) && (retval.elem(si.length(),0) != 0)) - si.elem(si.length()-1,0) ^= a.alpha_to(modn(a.index_of( - a.elem(si.length(),0))+a.index_of(retval.elem(i,0)), + si(si.length()-1,0) = 0; + if ((a(si.length(),0) != 0) && (retval(si.length(),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.elem(si.length(),0) != 0) && (x.elem(i,0) != 0)) - si.elem(si.length()-1,0) ^= b.alpha_to(modn(b.index_of( - b.elem(si.length(),0))+ b.index_of(x.elem(i,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.elem(0,0) = 0; - if ((a.elem(1,0) != 0) && (retval.elem(i,0) != 0)) - si.elem(0,0) ^= a.alpha_to(modn(a.index_of(a.elem(1,0))+ - a.index_of(retval.elem(i,0)),a.m(),a.n())); - if ((b.elem(1,0) != 0) && (x.elem(i,0) != 0)) - si.elem(0,0) ^= b.alpha_to(modn(b.index_of(b.elem(1,0))+ - b.index_of(x.elem(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.elem(i,0) = si.elem(0,0); - if ((b.elem(0,0) != 0) && (x.elem(i,0) != 0)) - retval.elem(i,0) ^= b.alpha_to(modn(b.index_of(b.elem(0,0)) + - b.index_of(x.elem(i,0)),b.m(),b.n())); + retval(i,0) = si(0,0); + if ((b(0,0) != 0) && (x(i,0) != 0)) + retval(i,0) ^= b.alpha_to(modn(b.index_of(b(0,0)) + + b.index_of(x(i,0)),b.m(),b.n())); if (si.length() > 1) { for (int j = 0; j < si.length() - 1; j++) { - si.elem(j,0) = si.elem(j+1,0); - if ((b.elem(j+1,0) != 0) && (x.elem(i,0) != 0)) - si.elem(j,0) ^= b.alpha_to(modn(b.index_of(b.elem(j+1,0)) + - b.index_of(x.elem(i,0)),b.m(),b.n())); + 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.elem(si.length()-1,0) = 0; - if ((b.elem(si.length(),0) != 0) && (x.elem(i,0) != 0)) - si.elem(si.length()-1,0) ^= b.alpha_to(modn(b.index_of( - b.elem(si.length(),0)) + b.index_of(x.elem(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.elem(0,0) = 0; - if ((b.elem(1,0) != 0) && (x.elem(i,0) != 0)) - si.elem(0,0) ^= b.alpha_to(modn(b.index_of(b.elem(1,0)) + - b.index_of(x.elem(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.elem(0,0) != 0) && (x.elem(i,0) != 0)) - retval.elem(i,0) = b.alpha_to(modn(b.index_of(b.elem(0,0)) + - b.index_of(x.elem(i,0)), b.m(),b.n())); + if ((b(0,0) != 0) && (x(i,0) != 0)) + retval(i,0) = b.alpha_to(modn(b.index_of(b(0,0)) + + b.index_of(x(i,0)), b.m(),b.n())); return retval; } @@ -1175,23 +1205,23 @@ // Create polynomial of right length. genpoly = galois(nroots+1,1,0,m,primpoly); - genpoly.elem(nroots,0) = 1; + genpoly(nroots,0) = 1; int i,root; for (i = 0,root=fcr*prim; i < nroots; i++,root += prim) { - genpoly.elem(nroots-i-1,0) = 1; + genpoly(nroots-i-1,0) = 1; // Multiply genpoly by @**(root + x) for (int j = i; j > 0; j--){ int k = nroots - j; - if (genpoly.elem(k,0) != 0) - genpoly.elem(k,0) = genpoly.elem(k+1,0) ^ genpoly.alpha_to( - modn(genpoly.index_of(genpoly.elem(k,0)) + root, m, n)); + 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.elem(k,0) = genpoly.elem(k+1,0); + genpoly(k,0) = genpoly(k+1,0); } // genpoly(nroots,0) can never be zero - genpoly.elem(nroots,0) = genpoly.alpha_to(modn(genpoly.index_of( - genpoly.elem(nroots,0)) + root, m, n)); + genpoly(nroots,0) = genpoly.alpha_to(modn(genpoly.index_of( + genpoly(nroots,0)) + root, m, n)); } } else { @@ -1211,11 +1241,11 @@ } } - int norm = genpoly.elem(0,0); + int norm = genpoly(0,0); // Take logarithm of generator polynomial, for faster coding for (int i = 0; i < nroots+1; i++) - genpoly.elem(i,0) = genpoly.index_of(genpoly.elem(i,0)); + genpoly(i,0) = genpoly.index_of(genpoly(i,0)); // Add space for parity block msg.resize(nsym,n,0); @@ -1234,24 +1264,24 @@ 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.elem(0,0) ^ msg.elem(l,i)); + int feedback = par.index_of(par(0,0) ^ msg(l,i)); if (feedback != nn) { if (norm != 1) - feedback = modn(nn-genpoly.elem(0,0)+feedback, m, nn); + feedback = modn(nn-genpoly(0,0)+feedback, m, nn); for (int j = 1; j < nroots; j++) - par.elem(j,0) ^= par.alpha_to(modn(feedback + - genpoly.elem(j,0), m, nn)); + par(j,0) ^= par.alpha_to(modn(feedback + + genpoly(j,0), m, nn)); } for (int j = 1; j < nroots; j++) - par.elem(j-1,0) = par.elem(j,0); + par(j-1,0) = par(j,0); if (feedback != nn) - par.elem(nroots-1,0) = par.alpha_to(modn(feedback+ - genpoly.elem(nroots,0), m, nn)); + par(nroots-1,0) = par.alpha_to(modn(feedback+ + genpoly(nroots,0), m, nn)); else - par.elem(nroots-1,0) = 0; + par(nroots-1,0) = 0; } for (int j = 0; j < nroots; j++) - msg.elem(l,k+j) = par.elem(j,0); + msg(l,k+j) = par(j,0); } } else { for (int l = 0; l < nsym; l++) { @@ -1263,24 +1293,24 @@ 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.elem(0,0) ^ msg.elem(l,i-1)); + int feedback = par.index_of(par(0,0) ^ msg(l,i-1)); if (feedback != nn) { if (norm != 1) - feedback = modn(nn-genpoly.elem(0,0)+feedback, m, nn); + feedback = modn(nn-genpoly(0,0)+feedback, m, nn); for (int j = 1; j < nroots; j++) - par.elem(j,0) ^= par.alpha_to(modn(feedback + - genpoly.elem(j,0), m, nn)); + par(j,0) ^= par.alpha_to(modn(feedback + + genpoly(j,0), m, nn)); } for (int j = 1; j < nroots; j++) - par.elem(j-1,0) = par.elem(j,0); + par(j-1,0) = par(j,0); if (feedback != nn) - par.elem(nroots-1,0) = par.alpha_to(modn(feedback+ - genpoly.elem(nroots,0), m, nn)); + par(nroots-1,0) = par.alpha_to(modn(feedback+ + genpoly(nroots,0), m, nn)); else - par.elem(nroots-1,0) = 0; + par(nroots-1,0) = 0; } for (int j = 0; j < nroots; j++) - msg.elem(l,j) = par.elem(nroots-j-1,0); + msg(l,j) = par(nroots-j-1,0); } } @@ -1323,25 +1353,25 @@ /* form the syndromes; i.e., evaluate data(x) at roots of g(x) */ if (msb_first) { for(i=0;i<nroots;i++) - s[i] = data.elem(drow,0); + s[i] = data(drow,0); for(j=1;j<n;j++) for(i=0;i<nroots;i++) if(s[i] == 0) - s[i] = data.elem(drow,j); + s[i] = data(drow,j); else - s[i] = data.elem(drow,j) ^ data.alpha_to(modn(data.index_of(s[i]) + + 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.elem(drow,n-1); + 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.elem(drow,j-1); + s[i] = data(drow,j-1); else - s[i] = data.elem(drow,j-1) ^ data.alpha_to(modn(data.index_of(s[i]) + + s[i] = data(drow,j-1) ^ data.alpha_to(modn(data.index_of(s[i]) + (fcr+i)*prim, m, n)); } @@ -1667,13 +1697,13 @@ OCTAVE_LOCAL_BUFFER(int, roots, nroots); for (int j=0; j <=nn; j++) { // Evaluate generator polynomial at j - int val = genpoly.elem(0,0); + int val = genpoly(0,0); int indx = genpoly.index_of(j); for (int i=0; i<nroots; i++) { if (val == 0) - val = genpoly.elem(i+1,0); + val = genpoly(i+1,0); else - val = genpoly.elem(i+1,0) ^ genpoly.alpha_to(modn(indx + + val = genpoly(i+1,0) ^ genpoly.alpha_to(modn(indx + genpoly.index_of(val), m, nn)); } if (val == 0) { @@ -1737,7 +1767,7 @@ if (parity_at_end) for (int l = 0; l < nsym; l++) for (int i=n; i > 0; i--) - code.elem(l,i+nn-n-1) = code.elem(l,i-1); + code(l,i+nn-n-1) = code(l,i-1); } for (int l = 0; l < nsym; l++) @@ -1747,18 +1777,18 @@ if (parity_at_end) for (int l = 0; l < nsym; l++) for (int i=0; i > n; i--) - code.elem(l,i) = code.elem(l,i+nn-n); + code(l,i) = code(l,i+nn-n); code.resize(nsym,n,0); } if (parity_at_end) { for (int l = 0; l < nsym; l++) for (int i=0; i < k; i++) - msg.elem(l,i) = code.elem(l,i); + msg(l,i) = code(l,i); } else { for (int l = 0; l < nsym; l++) for (int i=0; i < k; i++) - msg.elem(l,i) = code.elem(l,nroots+i); + msg(l,i) = code(l,nroots+i); } retval(0) = new octave_galois (msg); @@ -1892,12 +1922,12 @@ c.resize(nc+1,m); cs.resize(nc+1); - c.elem(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.elem(nc,cs(nc)) = r; + c(nc,cs(nc)) = r; found(c.alpha_to(r)-1) = 1; cs(nc) += 1; } @@ -1917,10 +1947,10 @@ for (int j = 2*(t-1); j<2*t; j++) { int flag = 0; for (int l=0; l<cs(i); l++) { - if (c.elem(i,l) == j+1) { + if (c(i,l) == j+1) { f.resize(1,nf+cs(i)); for (int ll=0; ll<cs(i); ll++) - f.elem(0,nf+ll) = c.elem(i,ll); + f(0,nf+ll) = c(i,ll); found(i) = 0; nf += cs(i); flag = 1; @@ -1941,21 +1971,21 @@ // Create polynomial of right length. genpoly = galois(nf+1,1,0,m); - genpoly.elem(0,0) = 1; + genpoly(0,0) = 1; for (int i = 0; i < nf; i++) { - genpoly.elem(i+1,0) = 1; + genpoly(i+1,0) = 1; // Multiply genpoly by @**(root + x) for (int l = i; l > 0; l--){ - if (genpoly.elem(l,0) != 0) - genpoly.elem(l,0) = genpoly.elem(l-1,0) ^ genpoly.alpha_to( - modn(genpoly.index_of(genpoly.elem(l,0)) + f.elem(0,i), m, n)); + 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.elem(l,0) = genpoly.elem(l-1,0); + genpoly(l,0) = genpoly(l-1,0); } // genpoly(0,0) can never be zero - genpoly.elem(0,0) = genpoly.alpha_to(modn(genpoly.index_of( - genpoly.elem(0,0)) + f.elem(0,i), m, n)); + genpoly(0,0) = genpoly.alpha_to(modn(genpoly.index_of( + genpoly(0,0)) + f(0,i), m, n)); } } @@ -1973,18 +2003,18 @@ if (parity_at_end) { for (int l = 0; l < nsym; l++) { for (int i = 0; i < k; i++) { - int feedback = (int)msg.elem(l,i) ^ (int)msg.elem(l,k); + int feedback = (int)msg(l,i) ^ (int)msg(l,k); if (feedback != 0) { for (int j = 0; j < nn-k-1; j++) - if (genpoly.elem(nn-k-j-1,0) != 0) - msg.elem(l,k+j) = (int)msg.elem(l,k+j+1) ^ feedback; + if (genpoly(nn-k-j-1,0) != 0) + msg(l,k+j) = (int)msg(l,k+j+1) ^ feedback; else - msg.elem(l,k+j) = msg.elem(l,k+j+1); - msg.elem(l,nn-1) = genpoly.elem(0,0) & feedback; + 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.elem(l,j) = msg.elem(l,j+1); - msg.elem(l,nn-1) = 0; + msg(l,j) = msg(l,j+1); + msg(l,nn-1) = 0; } } } @@ -1998,18 +2028,18 @@ for (int l = 0; l < nsym; l++) { for (int i = k-1; i >= 0; i--) { - int feedback = (int)msg.elem(l,nn-k+i) ^ (int)msg.elem(l,nn-k-1); + 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.elem(j,0) != 0) - msg.elem(l,j) = (int)msg.elem(l,j-1) ^ feedback; + if (genpoly(j,0) != 0) + msg(l,j) = (int)msg(l,j-1) ^ feedback; else - msg.elem(l,j) = msg.elem(l,j-1); - msg.elem(l,0) = genpoly.elem(0,0) & feedback; + 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.elem(l,j) = msg.elem(l,j-1); - msg.elem(l,0) = 0; + msg(l,j) = msg(l,j-1); + msg(l,0) = 0; } } } @@ -2151,10 +2181,10 @@ for (int i = 1; i <= t2; i++) { for (int j = 0; j < nn; j++) { if (parity_at_end) { - if (code.elem(lsym,nn-j-1) != 0) + if (code(lsym,nn-j-1) != 0) s(i) ^= tables.alpha_to(modn(i*j,m,n)); } else { - if (code.elem(lsym,j) != 0) + if (code(lsym,j) != 0) s(i) ^= tables.alpha_to(modn(i*j,m,n)); } } @@ -2295,10 +2325,10 @@ nerr(lsym) = l(u); for (int i = 0; i < l(u); i++) if (parity_at_end) - code.elem(lsym,nn-loc(i)-1) = - (int)code.elem(lsym,nn-loc(i)-1) ^ 1; + code(lsym,nn-loc(i)-1) = + (int)code(lsym,nn-loc(i)-1) ^ 1; else - code.elem(lsym,loc(i)) = (int)code.elem(lsym,loc(i)) ^ 1; + code(lsym,loc(i)) = (int)code(lsym,loc(i)) ^ 1; } else /* elp has degree >t hence cannot solve */ nerr(lsym) = -1; } else @@ -2310,11 +2340,11 @@ if (parity_at_end) { for (int l = 0; l < nsym; l++) for (int i = 0; i < k; i++) - msg.elem(l,i) = code.elem(l,i); + msg(l,i) = code(l,i); } else { for (int l = 0; l < nsym; l++) for (int i=0; i < k; i++) - msg.elem(l,i) = code.elem(l,nn-k+i); + msg(l,i) = code(l,nn-k+i); } retval(0) = octave_value(msg);
--- a/main/comm/gftable.m Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/gftable.m Tue Apr 01 19:37:15 2003 +0000 @@ -15,7 +15,7 @@ ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ## -*- texinfo -*- -## @deftypefn {Function File} gftable (@var{m}, @var{primpoly}) +## @deftypefn {Function File} {} gftable (@var{m}, @var{primpoly}) ## ## This function exists for compatiability with matlab. As the octave galois ## fields store a copy of the lookup tables for every field in use internally,
--- a/main/comm/ov-galois.cc Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/ov-galois.cc Tue Apr 01 19:37:15 2003 +0000 @@ -58,7 +58,7 @@ Matrix data(r,c); for (int i=0; i< r; i++) for (int j=0; j< c; j++) - data(i,j) = (double)gval.elem(i,j); + data(i,j) = (double)gval(i,j); retval(0) = octave_value (data); } #ifdef GALOIS_DISP_PRIVATES @@ -249,7 +249,7 @@ for (int i = 0; i < gval.rows(); i++) for (int j = 0; j < gval.columns(); j++) - data(i,j) = (double)gval.elem(i,j); + data(i,j) = (double)gval(i,j); octave_print_internal (os, data, false, current_print_indent_level ()); newline (os); @@ -321,7 +321,7 @@ retval.resize(rows(),columns()); for (int i=0; i<rows(); i++) for (int j=0; j<columns(); j++) - retval(i,j) = gval.elem(i,j); + retval(i,j) = gval(i,j); return retval; }
--- a/main/comm/rsdecof.m Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/rsdecof.m Tue Apr 01 19:37:15 2003 +0000 @@ -15,8 +15,8 @@ ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ## -*- texinfo -*- -## @deftypefn {Function File} rsdecof (@var{in},@var{out}) -## @deftypefnx {Function File} rsdecof (@var{in},@var{out},@var{t}) +## @deftypefn {Function File} {} rsdecof (@var{in},@var{out}) +## @deftypefnx {Function File} {} rsdecof (@var{in},@var{out},@var{t}) ## ## Decodes an ascii file using a Reed-Solomon coder. The input file is ## defined by @var{in} and the result is written to the output file @var{out}.
--- a/main/comm/rsencof.m Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/rsencof.m Tue Apr 01 19:37:15 2003 +0000 @@ -15,9 +15,9 @@ ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ## -*- texinfo -*- -## @deftypefn {Function File} rsencof (@var{in},@var{out}) -## @deftypefnx {Function File} rsencof (@var{in},@var{out},@var{t}) -## @deftypefnx {Function File} rsencof (@var{...},@var{pad}) +## @deftypefn {Function File} {} rsencof (@var{in},@var{out}) +## @deftypefnx {Function File} {} rsencof (@var{in},@var{out},@var{t}) +## @deftypefnx {Function File} {} rsencof (@var{...},@var{pad}) ## ## Encodes an ascii file using a Reed-Solomon coder. The input file is ## defined by @var{in} and the result is written to the output file @var{out}.
--- a/main/comm/scatterplot.m Mon Mar 31 20:28:10 2003 +0000 +++ b/main/comm/scatterplot.m Tue Apr 01 19:37:15 2003 +0000 @@ -15,11 +15,11 @@ ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ## -*- texinfo -*- -## @deftypefn {Function File} scatterplot (@var{x}) -## @deftypefnx {Function File} scatterplot (@var{x},@var{n}) -## @deftypefnx {Function File} scatterplot (@var{x},@var{n},@var{off}) -## @deftypefnx {Function File} scatterplot (@var{x},@var{n},@var{off},@var{str}) -## @deftypefnx {Function File} scatterplot (@var{x},@var{n},@var{off},@var{str},@var{h}) +## @deftypefn {Function File} {} scatterplot (@var{x}) +## @deftypefnx {Function File} {} scatterplot (@var{x},@var{n}) +## @deftypefnx {Function File} {} scatterplot (@var{x},@var{n},@var{off}) +## @deftypefnx {Function File} {} scatterplot (@var{x},@var{n},@var{off},@var{str}) +## @deftypefnx {Function File} {} scatterplot (@var{x},@var{n},@var{off},@var{str},@var{h}) ## @deftypefnx {Function File} {@var{h} =} scatterplot (@var{...}) ## ## Display the scatter plot of a signal. The signal @var{x} can be either in