/*
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 2, or (at your option) any
later version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with this program; see the file COPYING. If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
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 (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
#pragma implementation
#endif
#include <iostream>
#include "galois.h"
#include "galoisfield.h"
#include "galois-def.h"
galois_field_list stored_galois_fields;
// galois class
galois::galois (const MArray2<int>& a, const int& _m, const int& _primpoly) : MArray2<int> (a.rows(), a.cols()), field (NULL) {
int _n = (1<<_m) - 1;
// Check the validity of the data in the matrix
for (int i=0; i<rows(); i++) {
for (int j=0; j<columns(); j++) {
if ((a(i,j) < 0) || (a(i,j) > _n)) {
gripe_range_galois(_m);
return;
}
xelem(i,j) = (int)a(i,j);
}
}
field = stored_galois_fields.create_galois_field(_m, _primpoly);
}
galois::galois (const Matrix& a, const int& _m, const int& _primpoly) : MArray2<int> (a.rows(), a.cols()), field (NULL) {
int _n = (1<<_m) - 1;
// Check the validity of the data in the matrix
for (int i=0; i<rows(); i++) {
for (int j=0; j<columns(); j++) {
if ((a(i,j) < 0) || (a(i,j) > _n)) {
gripe_range_galois(_m);
return;
}
if ((a(i,j) - (double)((int)a(i,j))) != 0.) {
gripe_integer_galois();
return;
}
xelem(i,j) = (int)a(i,j);
}
}
field = stored_galois_fields.create_galois_field(_m, _primpoly);
}
galois::galois (int nr, int nc, const int& val, const int& _m, const int& _primpoly) : MArray2<int> (nr, nc, val), field (NULL) {
int _n = (1<<_m) - 1;
// Check the validity of the data in the matrix
if ((val < 0) || (val > _n)) {
gripe_range_galois(_m);
return;
}
field = stored_galois_fields.create_galois_field(_m, _primpoly);
}
galois::galois (int nr, int nc, double val, const int& _m, const int& _primpoly) : MArray2<int> (nr, nc, (int)val), field (NULL) {
int _n = (1<<_m) - 1;
// Check the validity of the data in the matrix
if ((val < 0) || (val > _n)) {
gripe_range_galois(_m);
return;
}
if ((val - (double)((int)val)) != 0.) {
gripe_integer_galois();
return;
}
field = stored_galois_fields.create_galois_field(_m, _primpoly);
}
galois :: galois (const galois& a) :
MArray2<int>(a) {
if (!a.have_field()) {
gripe_copy_invalid_galois();
field = NULL;
return;
}
// This call to create_galois_field will just increment the usage counter
field = stored_galois_fields.create_galois_field(a.m(), a.primpoly());
}
galois :: ~galois (void)
{
stored_galois_fields.delete_galois_field (field);
field = NULL;
}
galois & galois::operator = (const galois& t)
{
if (!t.have_field()) {
gripe_copy_invalid_galois();
if (have_field())
stored_galois_fields.delete_galois_field (field);
field = NULL;
return *this;
}
if (have_field()) {
if ((m() != t.m()) || (primpoly() != t.primpoly())) {
stored_galois_fields.delete_galois_field (field);
field = stored_galois_fields.create_galois_field(t.m(), t.primpoly());
}
} else
field = stored_galois_fields.create_galois_field(t.m(), t.primpoly());
// Copy the data
MArray2<int>::operator = (t);
return *this;
}
galois&
galois::operator += (const galois& a)
{
int nr = rows ();
int nc = cols ();
int a_nr = a.rows ();
int a_nc = a.cols ();
if (have_field() && a.have_field()) {
if ((m() != a.m()) || (primpoly() != a.primpoly())) {
gripe_differ_galois ();
return *this;
}
} else {
gripe_invalid_galois ();
return *this;
}
if (nr != a_nr || nc != a_nc)
{
gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
return *this;
}
for (int i=0; i<rows(); i++)
for (int j=0; j<columns(); j++)
xelem(i,j) ^= a (i, j);
return *this;
}
galois&
galois::operator -= (const galois& a)
{
int nr = rows ();
int nc = cols ();
int a_nr = a.rows ();
int a_nc = a.cols ();
if (have_field() && a.have_field()) {
if ((m() != a.m()) || (primpoly() != a.primpoly())) {
gripe_differ_galois ();
return *this;
}
} else {
gripe_invalid_galois ();
return *this;
}
if (nr != a_nr || nc != a_nc)
{
gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
return *this;
}
for (int i=0; i<rows(); i++)
for (int j=0; j<columns(); j++)
xelem(i,j) ^= a (i, j);
return *this;
}
galois galois::index (idx_vector& i, int resize_ok, const int& rfv) const
{
galois retval(MArray2<int>::index(i, resize_ok, rfv), m(), primpoly());
return retval;
}
galois galois::index (idx_vector& i, idx_vector& j, int resize_ok,
const int& rfv) const
{
galois retval(MArray2<int>::index(i, j, resize_ok, rfv), m(), primpoly());
return retval;
}
#ifdef HAVE_OLD_OCTAVE_CONCAT
galois concat (const galois& ra, const galois& rb, const Array<int>& ra_idx)
{
galois retval (ra);
if (rb.numel() > 0)
retval.insert (rb, ra_idx(0), ra_idx(1));
return retval;
}
galois concat (const galois& ra, const Matrix& rb, const Array<int>& ra_idx)
{
galois retval (ra);
if (ra.numel() == 1)
return retval;
galois tmp (0, 0, 0, ra.m(), ra.primpoly());
int _n = (1<<ra.m()) - 1;
int r = rb.rows();
int c = rb.columns();
tmp.resize (r, c);
// Check the validity of the data in the matrix
for (int i=0; i<r; i++) {
for (int j=0; j<c; j++) {
if ((rb(i,j) < 0) || (rb(i,j) > _n)) {
gripe_range_galois(ra.m());
return retval;
}
if ((rb(i,j) - (double)((int)rb(i,j))) != 0.) {
gripe_integer_galois();
return retval;
}
tmp(i,j) = (int)rb(i,j);
}
}
retval.insert (tmp, ra_idx(0), ra_idx(1));
return retval;
}
#endif
#ifdef HAVE_OCTAVE_CONCAT
galois
galois::concat (const galois& rb, const Array<int>& ra_idx)
{
if (rb.numel() > 0)
insert (rb, ra_idx(0), ra_idx(1));
return *this;
}
galois
galois::concat (const Matrix& rb, const Array<int>& ra_idx)
{
if (numel() == 1)
return *this;
galois tmp (0, 0, 0, m(), primpoly());
int _n = (1<<m()) - 1;
int r = rb.rows();
int c = rb.columns();
tmp.resize (r, c);
// Check the validity of the data in the matrix
for (int i=0; i<r; i++) {
for (int j=0; j<c; j++) {
if ((rb(i,j) < 0) || (rb(i,j) > _n)) {
gripe_range_galois(m());
return *this;
}
if ((rb(i,j) - (double)((int)rb(i,j))) != 0.) {
gripe_integer_galois();
return *this;
}
tmp(i,j) = (int)rb(i,j);
}
}
insert (tmp, ra_idx(0), ra_idx(1));
return *this;
}
#endif
#if defined (HAVE_OLD_OCTAVE_CONCAT) || defined (HAVE_OCTAVE_CONCAT)
galois concat (const Matrix& ra, const galois& rb, const Array<int>& ra_idx)
{
galois retval (0, 0, 0, rb.m(), rb.primpoly());
int _n = (1<<rb.m()) - 1;
int r = ra.rows();
int c = ra.columns();
retval.resize (r, c);
if (ra.numel() < 1)
return retval;
// XXX FIXME XXX
// Check the validity of the data in the matrix. This is problematic
// as "ra" is not initialized on the initial resize and so contains
// random data that will be replaced. Humm, disable for now
for (int i=0; i<r; i++) {
for (int j=0; j<c; j++) {
#if 0
if ((ra(i,j) < 0) || (ra(i,j) > _n)) {
gripe_range_galois(rb.m());
return retval;
}
if ((ra(i,j) - (double)((int)ra(i,j))) != 0.) {
gripe_integer_galois();
return retval;
}
retval(i,j) = (int)ra(i,j);
#else
int tmp = (int)ra(i,j);
if (tmp < 0)
retval(i,j) = 0;
else if (tmp > _n)
retval(i,j) = _n;
else
retval(i,j) = tmp;
#endif
}
}
retval.insert (rb, ra_idx(0), ra_idx(1));
return retval;
}
galois& galois::insert (const galois& t, int r, int c)
{
if ((m() != t.m()) || (primpoly() != t.primpoly()))
(*current_liboctave_error_handler) ("inserted galois variable must be in the same field");
else
Array<int>::insert (t, r, c);
return *this;
}
#endif
galois galois::diag (void) const
{
return diag(0);
}
galois galois::diag (int k) const
{
int nnr = rows ();
int nnc = cols ();
galois retval(0, 0, 0, m(), primpoly());
if (k > 0)
nnc -= k;
else if (k < 0)
nnr += k;
if (nnr > 0 && nnc > 0)
{
int ndiag = (nnr < nnc) ? nnr : nnc;
retval.resize(ndiag, 1);
if (k > 0)
{
for (int i = 0; i < ndiag; i++)
retval (i,0) = xelem (i, i+k);
}
else if ( k < 0)
{
for (int i = 0; i < ndiag; i++)
retval (i,0) = xelem (i-k, i);
}
else
{
for (int i = 0; i < ndiag; i++)
retval (i,0) = xelem (i, i);
}
}
else
error("diag: requested diagonal out of range");
return retval;
}
// unary operations
boolMatrix
galois::operator ! (void) const
{
int nr = rows ();
int nc = cols ();
boolMatrix b (nr, nc);
for (int j = 0; j < nc; j++)
for (int i = 0; i < nr; i++)
b (i, j) = ! xelem (i, j);
return b;
}
galois
galois :: transpose (void) const
{
galois a(Matrix(0,0), m(), primpoly());
int d1 = rows();
int d2 = cols();
a.resize(d2,d1);
for (int j = 0; j < d2; j++)
for (int i = 0; i < d1; i++)
a (j, i) = xelem (i, j);
return a;
}
static inline int modn(int x, int m, int n)
{
while (x >= n) {
x -= n;
x = (x >> m) + (x & n);
}
return x;
}
galois elem_pow (const galois& a, const galois& b)
{
int a_nr = a.rows ();
int a_nc = a.cols ();
galois result(a_nr, a_nc, 0, a.m(), a.primpoly());
int b_nr = b.rows ();
int b_nc = b.cols ();
if (a.have_field() && b.have_field()) {
if ((a.m() != b.m()) || (a.primpoly() != b.primpoly())) {
gripe_differ_galois ();
return galois ();
}
} else {
gripe_invalid_galois ();
return galois ();
}
if (a_nr == 1 && a_nc == 1)
{
result.resize(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)
{
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()));
}
return result;
}
galois elem_pow (const galois& a, const Matrix& b)
{
int a_nr = a.rows ();
int a_nc = a.cols ();
galois result(a_nr, a_nc, 0, a.m(), a.primpoly());
int b_nr = b.rows ();
int b_nc = b.cols ();
if (b_nr == 1 && b_nc == 1)
return elem_pow(a, b(0,0));
if (a_nr != b_nr || a_nc != b_nc)
{
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()));
}
return result;
}
galois elem_pow (const galois& a, double b)
{
int a_nr = a.rows ();
int a_nc = a.cols ();
galois result(a_nr, a_nc, 0, a.m(), a.primpoly());
int bi = (int) b;
if ((double)bi != b) {
gripe_integer_galois ();
return galois ();
}
while (bi < 0)
bi += a.n();
for (int j = 0; j < a_nc; j++)
for (int i = 0; i < a_nr; i++)
{
if (bi == 0)
result (i,j) = 1;
else if (a(i,j) != 0)
result (i,j) = a.alpha_to(modn(a.index_of(a(i,j)) *
bi,a.m(),a.n()));
}
return result;
}
galois elem_pow (const galois &a, int b)
{
int a_nr = a.rows ();
int a_nc = a.cols ();
galois result(a_nr, a_nc, 0, a.m(), a.primpoly());
while (b < 0)
b += a.n();
for (int j = 0; j < a_nc; j++)
for (int i = 0; i < a_nr; i++)
{
if (b == 0)
result (i,j) = 1;
else if (a(i,j) != 0)
result (i,j) = a.alpha_to(modn(a.index_of(a(i,j)) * b,
a.m(),a.n()));
}
return result;
}
galois pow (const galois& a, double b)
{
int bi = (int)b;
if ((double)bi != b) {
gripe_integer_power_galois ();
return galois ();
}
return pow(a, bi);
}
galois pow (const galois& a, const galois& b)
{
int nr = b.rows ();
int nc = b.cols ();
if (a.have_field() && b.have_field()) {
if ((a.m() != b.m()) || (a.primpoly() != b.primpoly())) {
gripe_differ_galois ();
return galois ();
}
} else {
gripe_invalid_galois ();
return galois ();
}
if (nr != 1 || nc != 1) {
gripe_square_galois();
return galois ();
} else
return pow (a, b(0,0));
}
galois pow (const galois& a, int b)
{
galois retval;
int nr = a.rows ();
int nc = a.cols ();
if (!a.have_field()) {
gripe_invalid_galois ();
return retval;
}
if (nr == 0 || nc == 0 || nr != nc)
gripe_square_galois();
else if (b == 0)
{
retval = galois(nr, nc, 0, a.m(), a.primpoly());
for (int i =0; i<nr; i++)
retval(i,i) = 1;
}
else
{
galois atmp;
if (b < 0 ) {
atmp = a.inverse();
b = abs(b);
} else
atmp = a;
retval = atmp;
b--;
while (b > 0)
{
if (b & 1)
retval = retval * atmp;
b >>= 1;
if (b > 0)
atmp = atmp * atmp;
}
}
return retval;
}
galois
operator * (const Matrix& a, const galois& b)
{
galois tmp (a, b.m(), b.primpoly());
OCTAVE_QUIT;
return tmp * b;
}
galois
operator * (const galois& a, const Matrix& b)
{
galois tmp (b, a.m(), a.primpoly());
OCTAVE_QUIT;
return a * tmp;
}
galois
operator * (const galois& a, const galois& b)
{
if (a.have_field() && b.have_field()) {
if ((a.m() != b.m()) || (a.primpoly() != b.primpoly())) {
gripe_differ_galois ();
return galois ();
}
} else {
gripe_invalid_galois ();
return galois ();
}
int a_nr = a.rows ();
int a_nc = a.cols ();
int b_nr = b.rows ();
int b_nc = b.cols ();
if ((a_nr == 1 && a_nc == 1) || (b_nr == 1 && b_nc == 1))
return product (a, b);
else if (a_nc != b_nr)
{
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)
{
// This is not optimum for referencing b, but can use vector
// to represent index(a(k,j)). Seems to be the fastest.
galois c(a_nr, 1, 0, a.m(), a.primpoly());
for (int j=0; j<b_nr; j++) {
for (int k=0; k<a_nr; k++)
c(k,0) = a.index_of(a(k,j));
for (int i=0; i<b_nc; i++)
if (b(j,i) != 0) {
int tmp = a.index_of(b(j,i));
for (int k=0; k<a_nr; k++) {
if (a(k,j) != 0)
retval(k,i) = retval(k,i) ^ a.alpha_to(
modn(tmp + c(k,0),a.m(),a.n()));
}
}
}
}
return retval;
}
}
// Other operators
boolMatrix
galois::all (int dim) const
{
MX_ALL_OP (dim);
}
boolMatrix
galois::any (int dim) const
{
MX_ANY_OP (dim);
}
galois
galois::prod (int dim) const
{
if (!have_field()) {
gripe_invalid_galois ();
return galois();
}
galois retval (0, 0, 0, m(), primpoly());
#define ROW_EXPR \
if ((retval (i, 0) == 0) || (elem(i,j) == 0)) \
retval (i, 0) = 0; \
else \
retval (i, 0) = alpha_to(modn(index_of(retval(i,0)) + \
index_of(elem(i,j)),m(),n()));
#define COL_EXPR \
if ((retval (0, j) == 0) || (elem(i,j) == 0)) \
retval (0, j) = 0; \
else \
retval (0, j) = alpha_to(modn(index_of(retval(0,j)) + \
index_of(elem(i,j)),m(),n()));
GALOIS_REDUCTION_OP(retval, ROW_EXPR, COL_EXPR, 1, 1);
return retval;
#undef ROW_EXPR
#undef COL_EXPR
}
galois
galois::sum (int dim) const
{
if (!have_field()) {
gripe_invalid_galois ();
return galois();
}
galois retval (0, 0, 0, m(), primpoly());
#define ROW_EXPR \
retval (i, 0) ^= elem(i,j);
#define COL_EXPR \
retval (0, j) ^= elem(i,j);
GALOIS_REDUCTION_OP (retval, ROW_EXPR, COL_EXPR, 0, 0);
return retval;
#undef ROW_EXPR
#undef COL_EXPR
}
galois
galois::sumsq (int dim) const
{
if (!have_field()) {
gripe_invalid_galois ();
return galois();
}
galois retval (0, 0, 0, m(), primpoly());
#define ROW_EXPR \
if (elem(i,j) != 0) \
retval (i, 0) ^= alpha_to(modn(2*index_of(elem(i,j)),m(),n()));
#define COL_EXPR \
if (elem(i,j) != 0) \
retval (0, j) ^= alpha_to(modn(2*index_of(elem(i,j)),m(),n()));
GALOIS_REDUCTION_OP (retval, ROW_EXPR, COL_EXPR, 0, 0);
return retval;
#undef ROW_EXPR
#undef COL_EXPR
}
galois
galois::sqrt (void) const
{
galois retval (*this);
int nr = rows ();
int nc = cols ();
for (int j=0; j<nc; j++) {
for (int i=0; i<nr; i++)
if (retval.index_of(retval(i,j)) & 1)
retval(i,j) = retval.alpha_to((retval.index_of(retval(i,j))
+ retval.n()) / 2);
else
retval(i,j) = retval.alpha_to(retval.index_of(retval(i,j))
/ 2);
}
return retval;
}
galois
galois::log (void) const
{
bool warned = false;
if (!have_field()) {
gripe_invalid_galois ();
return galois();
}
galois retval (*this);
int nr = rows ();
int nc = cols ();
for (int j=0; j<nc; j++)
for (int i=0; i<nr; i++) {
if (retval(i,j) == 0) {
if (!warned) {
warning("log of zero undefined in Galois field");
warned = true;
}
// How do I flag a NaN without either
// 1) Having to check everytime that the data is valid
// 2) Causing overflow in alpha_to or index_of!!
retval(i,j) = retval.index_of(retval(i,j));
} else
retval(i,j) = retval.index_of(retval(i,j));
}
return retval;
}
galois
galois::exp (void) const
{
bool warned = false;
if (!have_field()) {
gripe_invalid_galois ();
return galois();
}
galois retval (*this);
int nr = rows ();
int nc = cols ();
for (int j=0; j<nc; j++)
for (int i=0; i<nr; i++) {
if (retval(i,j) == n()) {
if (!warned) {
warning("warning: exp of 2^m-1 undefined in Galois field");
warned = true;
}
// How do I flag a NaN without either
// 1) Having to check everytime that the data is valid
// 2) Causing overflow in alpha_to or index_of!!
retval(i,j) = retval.alpha_to(retval(i,j));
} else
retval(i,j) = retval.alpha_to(retval(i,j));
}
return retval;
}
template class base_lu <galois, int, Matrix, double>;
void
LU::factor (const galois& a, const pivot_type& typ)
{
int a_nr = a.rows ();
int a_nc = a.cols ();
int mn = (a_nr > a_nc ? a_nc : a_nr);
ptype = typ;
info = 0;
ipvt.resize (mn);
a_fact = a;
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(i,j) > a_fact(jp,j))
jp = i;
} else {
for (int i = j+1; i < a_nc; i++)
if (a_fact(j,i) > a_fact(j,jp))
jp = i;
}
ipvt(j) = jp;
if (a_fact(jp,j) != 0) {
if (ptype == LU::ROW) {
// Apply the interchange to columns 1:NC.
if (jp != j)
for (int i = 0; i < a_nc; i++) {
int tmp = a_fact(j,i);
a_fact(j,i) = a_fact(jp,i);
a_fact(jp,i) = tmp;
}
} else {
// Apply the interchange to rows 1:NR.
if (jp != j)
for (int i = 0; i < a_nr; i++) {
int tmp = a_fact(i,j);
a_fact(i,j) = a_fact(i,jp);
a_fact(i,jp) = tmp;
}
}
// Compute elements J+1:M of J-th column.
if ( j < a_nr-1) {
int idxj = a_fact.index_of(a_fact(j,j));
for (int i = j+1; i < a_nr; i++) {
if (a_fact(i,j) == 0)
a_fact(i,j) = 0;
else
a_fact(i,j) = a_fact.alpha_to(modn(a_fact.index_of(
a_fact(i,j)) - idxj + a_fact.n(), a_fact.m(),
a_fact.n()));
}
}
} else {
info = 1;
}
if (j < mn-1) {
// Update trailing submatrix.
for (int i=j+1; i < a_nr; i++) {
if (a_fact(i,j) != 0) {
int idxi = a_fact.index_of(a_fact(i,j));
for (int k=j+1; k < a_nc; k++) {
if (a_fact(j,k) != 0)
a_fact(i,k) ^= a_fact.alpha_to(modn(a_fact.index_of(
a_fact(j,k)) + idxi, a_fact.m(), a_fact.n()));
}
}
}
}
}
}
galois
LU::L (void) const
{
int a_nr = a_fact.rows ();
int a_nc = a_fact.cols ();
int mn = (a_nr < a_nc ? a_nr : a_nc);
galois l (a_nr, mn, 0, a_fact.m(), a_fact.primpoly());
for (int i = 0; i < mn; i++)
l (i, i) = 1;
for (int j = 0; j < mn; j++)
for (int i = j+1; i < a_nr; i++)
l (i, j) = a_fact (i, j);
return l;
}
galois
LU::U (void) const
{
int a_nr = a_fact.rows ();
int a_nc = a_fact.cols ();
int mn = (a_nr < a_nc ? a_nr : a_nc);
galois u (mn, a_nc, 0, a_fact.m(), a_fact.primpoly());
for (int j = 0; j < a_nc; j++)
for (int i = 0; i < (j+1 > mn ? mn : j+1); i++)
u (i, j) = a_fact (i, j);
return u;
}
Matrix
LU::P (void) const
{
int a_nr = a_fact.rows ();
int a_nc = a_fact.cols ();
int n = (ptype == LU::ROW ? a_nr : a_nc);
Array<int> pvt (n);
for (int i = 0; i < n; i++)
pvt (i) = i;
for (int i = 0; i < ipvt.length(); i++)
{
int k = ipvt (i);
if (k != i)
{
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 (i, pvt (i)) = 1.0;
return p;
}
Array<int>
LU::IP (void) const
{
return Array<int> (ipvt);
}
galois
LU::A (void) const
{
return galois (a_fact);
}
galois
galois::inverse (void) const
{
int info;
return inverse(info);
}
galois
galois::inverse (int& info, int force) const
{
int nr = rows ();
int nc = cols ();
info = 0;
if (nr != nc || nr == 0 || nc == 0) {
(*current_liboctave_error_handler) ("inverse requires square matrix");
return galois ();
} else {
int info = 0;
// Solve with identity matrix to find the inverse.
galois btmp(nr, nr, 0, m(), primpoly());
for (int i=0; i < nr; i++)
btmp(i,i) = 1;
galois retval = solve(btmp, info, 0);
if (info == 0)
return retval;
else
return galois();
}
}
galois
galois::determinant (void) const
{
int info;
return determinant (info);
}
galois
galois::determinant (int& info) const
{
galois retval(1, 1, 0, m(), primpoly());
int nr = rows ();
int nc = cols ();
info = -1;
if (nr == 0 || nc == 0) {
info = 0;
retval(0,0) = 1;
} else {
LU fact (*this);
if ( ! fact.singular()) {
galois A (fact.A());
info = 0;
retval(0,0) = A(0,0);
for (int i=1; i<nr; i++) {
if ((retval (0, 0) == 0) || (A(i,i) == 0)) {
retval (0,0) = 0;
error("What the hell are we doing here!!!");
} else
retval (0,0) = alpha_to(modn(index_of(retval(0,0)) + \
index_of(A(i,i)),m(),n()));
}
}
}
return retval;
}
galois
galois::solve (const galois& b) const
{
int info;
return solve (b, info);
}
galois
galois::solve (const galois& b, int& info) const
{
return solve (b, info, 0);
}
galois
galois::solve (const galois& b, int& info,
solve_singularity_handler sing_handler) const
{
galois retval (b);
if (!have_field() || !b.have_field()) {
gripe_invalid_galois ();
return galois();
} else if ((m() != b.m()) || (primpoly() != b.primpoly())) {
gripe_differ_galois ();
return galois();
}
int nr = rows ();
int nc = cols ();
int b_nr = b.rows ();
int b_nc = b.cols ();
galois c(nr,1,0,m(),primpoly());
// if (nr == 0 || nc == 0 || nr != nc || nr != b_nr) {
if (nr == 0 || nc == 0 || nr != b_nr) {
(*current_liboctave_error_handler)
("matrix dimension mismatch solution of linear equations");
return galois();
} else if (nc > nr) {
// Under-determined system, use column interchanges.
LU fact ((*this), LU::COL);
if (fact.singular()) {
info = -1;
if (sing_handler)
sing_handler (0.0);
else
(*current_liboctave_error_handler)("galois matrix singular");
return galois();
} else {
galois A (fact.A());
Array<int> IP (fact.IP());
// Resize the number of solution rows if needed
if (nc > nr)
retval.resize(b_nr+nc-nr,b_nc,0);
//Solve L*X = B, overwriting B with X.
int mn = (nc < nr ? nc : nr);
for (int k=0; k<mn; k++) {
for (int i=k+1; i<nr; i++)
c(i,0) = index_of(A(i,k));
for (int j=0; j<b_nc; j++)
if (retval(k,j) != 0) {
int idx = index_of(retval(k,j));
for (int i=k+1; i<nr; i++)
if (A(i,k) != 0)
retval(i,j) ^= alpha_to(modn(c(i,0) + idx, m(), n()));
}
}
// Solve U*X = B, overwriting B with X.
for (int k=(nc < nr ? nc-1 : nr-1); k>=0; k--) {
int mn = k+1 < nr ? k+1 : nr;
for (int i=0; i<mn; i++)
c(i,0) = index_of(A(i,k));
mn = k < nr ? k : nr;
for (int j=0; j<b_nc; j++)
if (retval(k,j) != 0) {
retval(k,j) = alpha_to(modn(index_of(retval(k,j)) -
c(k,0) + n(), m(), n()));
int idx = index_of(retval(k,j));
for (int i=0; i<mn; i++)
if (A(i,k) != 0)
retval(i,j) ^= alpha_to(modn(c(i,0) + idx, m(), n()));
}
}
// Apply row interchanges to the right hand sides.
//for (int j=0; j<IP.length(); j++) {
for (int j=IP.length()-1; j>=0; j--) {
int piv = IP(j);
for (int i=0; i<b_nc; i++) {
int tmp = retval(j,i);
retval(j,i) = retval(piv,i);
retval(piv,i) = tmp;
}
}
}
} else {
LU fact (*this);
if (fact.singular()) {
info = -1;
if (sing_handler)
sing_handler (0.0);
else
(*current_liboctave_error_handler)("galois matrix singular");
return galois();
} else {
galois A (fact.A());
Array<int> IP (fact.IP());
// Apply row interchanges to the right hand sides.
for (int j=0; j<IP.length(); j++) {
int piv = IP(j);
for (int i=0; i<b_nc; i++) {
int tmp = retval(j,i);
retval(j,i) = retval(piv,i);
retval(piv,i) = tmp;
}
}
//Solve L*X = B, overwriting B with X.
int mn = (nc < nr ? nc : nr);
for (int k=0; k<mn; k++) {
for (int i=k+1; i<nr; i++)
c(i,0) = index_of(A(i,k));
for (int j=0; j<b_nc; j++)
if (retval(k,j) != 0) {
int idx = index_of(retval(k,j));
for (int i=k+1; i<nr; i++)
if (A(i,k) != 0)
retval(i,j) ^= alpha_to(modn(c(i,0) + idx, m(), n()));
}
}
// Solve U*X = B, overwriting B with X.
for (int k=(nc < nr ? nc-1 : nr-1); k>=0; k--) {
int mn = k+1 < nr ? k+1 : nr;
for (int i=0; i<mn; i++)
c(i,0) = index_of(A(i,k));
mn = k < nr ? k : nr;
for (int j=0; j<b_nc; j++)
if (retval(k,j) != 0) {
retval(k,j) = alpha_to(modn(index_of(retval(k,j)) -
c(k,0) + n(), m(), n()));
int idx = index_of(retval(k,j));
for (int i=0; i<mn; i++)
if (A(i,k) != 0)
retval(i,j) ^= alpha_to(modn(c(i,0) + idx, m(), n()));
}
}
// Resize the number of solution rows if needed
if (nc < nr)
retval.resize(b_nr+nc-nr,b_nc);
}
}
return retval;
}
galois
xdiv (const galois& a, const Matrix& b)
{
galois btmp(b, a.m(), a.primpoly());
return xdiv(a, btmp);
}
galois
xdiv (const Matrix& a, const galois& b)
{
galois atmp(a, b.m(), b.primpoly());
return xdiv(atmp, b);
}
galois
xdiv (const galois& a, const galois& b)
{
int info = 0;
int a_nc = a.cols ();
int b_nc = b.cols ();
// if ((a_nc != b_nc) || (b.rows () != b.cols ()))
if (a_nc != b_nc)
{
int a_nr = a.rows ();
int b_nr = b.rows ();
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)
return galois (result.transpose ());
else
return galois();
}
galois
xleftdiv (const galois& a, const Matrix& b)
{
galois btmp(b, a.m(), a.primpoly());
return xleftdiv(a, btmp);
}
galois
xleftdiv (const Matrix& a, const galois& b)
{
galois atmp(a, b.m(), b.primpoly());
return xleftdiv(atmp, b);
}
galois
xleftdiv (const galois& a, const galois& b)
{
int info = 0;
int a_nr = a.rows ();
int b_nr = b.rows ();
// if ((a_nr != b_nr) || (a.rows () != a.columns ()))
if (a_nr != b_nr)
{
int a_nc = a.cols ();
int b_nc = b.cols ();
gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc);
return galois();
}
galois result = a.solve (b, info, 0);
if (info == 0)
return result;
else
return galois();
}
MM_BIN_OPS1(galois, galois, galois, 1, 2, GALOIS)
MM_BIN_OPS1(galois, galois, Matrix, 1, 2, MATRIX)
MM_BIN_OPS1(galois, Matrix, galois, 2, 1, MATRIX)
MM_CMP_OPS1(galois, , galois, , 1, 2, GALOIS)
MM_CMP_OPS1(galois, , Matrix, , 1, 2, MATRIX)
MM_CMP_OPS1(Matrix, , galois, , 2, 1, MATRIX)
MM_BOOL_OPS1(galois, galois, 0.0, 1, 2, GALOIS)
MM_BOOL_OPS1(galois, Matrix, 0.0, 1, 2, MATRIX)
MM_BOOL_OPS1(Matrix, galois, 0.0, 2, 1, MATRIX)
/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
syntax highlighted by Code2HTML, v. 0.9.1