Loading...
Searching...
No Matches
bigintmat.cc
Go to the documentation of this file.
/*****************************************
* Computer Algebra System SINGULAR *
*****************************************/
/*
* ABSTRACT: class bigintmat: matrices of numbers.
* a few functions might be limited to bigint or euclidean rings.
*/
#include "misc/auxiliary.h"
#include "coeffs/bigintmat.h"
#include "misc/intvec.h"
#include "coeffs/rmodulon.h"
#include <cmath>
///create Z/nA of type n_Zn
static coeffs numbercoeffs(number n, coeffs c) // TODO: FIXME: replace with n_CoeffRingQuot1
{
mpz_t p;
number2mpz(n, c, p);
ZnmInfo *pp = new ZnmInfo;
pp->base = p;
pp->exp = 1;
coeffs nc = nInitChar(n_Zn, (void*)pp);
mpz_clear(p);
delete pp;
return nc;
}
//#define BIMATELEM(M,I,J) (M)[ (M).index(I,J) ]
bigintmat * bigintmat::transpose()
{
bigintmat * t = new bigintmat(col, row, basecoeffs());
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
t->set(j, i, BIMATELEM(*this,i,j));
}
}
return t;
}
void bigintmat::inpTranspose()
{
int n = row,
m = col,
nm = n<m?n : m; // the min, describing the square part of the matrix
//CF: this is not optimal, but so far, it seems to work
#define swap(_i, _j) \
int __i = (_i), __j=(_j); \
number c = v[__i]; \
v[__i] = v[__j]; \
v[__j] = c \
for (int i=0; i< nm; i++)
for (int j=i+1; j< nm; j++)
{
swap(i*m+j, j*n+i);
}
if (n<m)
for (int i=nm; i<m; i++)
for(int j=0; j<n; j++)
{
swap(j*n+i, i*m+j);
}
if (n>m)
for (int i=nm; i<n; i++)
for(int j=0; j<m; j++)
{
swap(i*m+j, j*n+i);
}
#undef swap
row = m;
col = n;
}
// Beginnt bei [0]
void bigintmat::set(int i, number n, const coeffs C)
{
assume (C == NULL || C == basecoeffs());
rawset(i, n_Copy(n, basecoeffs()), basecoeffs());
}
// Beginnt bei [1,1]
void bigintmat::set(int i, int j, number n, const coeffs C)
{
assume (C == NULL || C == basecoeffs());
assume (i > 0 && j > 0);
assume (i <= rows() && j <= cols());
set(index(i, j), n, C);
}
number bigintmat::get(int i) const
{
assume (i >= 0);
assume (i<rows()*cols());
return n_Copy(v[i], basecoeffs());
}
number bigintmat::view(int i) const
{
assume (i >= 0);
assume (i<rows()*cols());
return v[i];
}
number bigintmat::get(int i, int j) const
{
assume (i > 0 && j > 0);
assume (i <= rows() && j <= cols());
return get(index(i, j));
}
number bigintmat::view(int i, int j) const
{
assume (i >= 0 && j >= 0);
assume (i <= rows() && j <= cols());
return view(index(i, j));
}
// Ueberladener *=-Operator (für int und bigint)
// Frage hier: *= verwenden oder lieber = und * einzeln?
void bigintmat::operator*=(int intop)
{
number iop = n_Init(intop, basecoeffs());
inpMult(iop, basecoeffs());
n_Delete(&iop, basecoeffs());
}
void bigintmat::inpMult(number bintop, const coeffs C)
{
assume (C == NULL || C == basecoeffs());
const int l = rows() * cols();
for (int i=0; i < l; i++)
n_InpMult(v[i], bintop, basecoeffs());
}
// Stimmen Parameter?
// Welche der beiden Methoden?
// Oder lieber eine comp-Funktion?
bool operator==(const bigintmat & lhr, const bigintmat & rhr)
{
if (&lhr == &rhr) { return true; }
if (lhr.cols() != rhr.cols()) { return false; }
if (lhr.rows() != rhr.rows()) { return false; }
if (lhr.basecoeffs() != rhr.basecoeffs()) { return false; }
const int l = (lhr.rows())*(lhr.cols());
for (int i=0; i < l; i++)
{
if (!n_Equal(lhr[i], rhr[i], lhr.basecoeffs())) { return false; }
}
return true;
}
bool operator!=(const bigintmat & lhr, const bigintmat & rhr)
{
return !(lhr==rhr);
}
// Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ?
bigintmat * bimAdd(bigintmat * a, bigintmat * b)
{
if (a->cols() != b->cols()) return NULL;
if (a->rows() != b->rows()) return NULL;
if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
const coeffs basecoeffs = a->basecoeffs();
int i;
bigintmat * bim = new bigintmat(a->rows(), a->cols(), basecoeffs);
for (i=a->rows()*a->cols()-1;i>=0; i--)
bim->rawset(i, n_Add((*a)[i], (*b)[i], basecoeffs), basecoeffs);
return bim;
}
bigintmat * bimAdd(bigintmat * a, long b)
{
const int mn = si_min(a->rows(),a->cols());
const coeffs basecoeffs = a->basecoeffs();
number bb=n_Init(b,basecoeffs);
int i;
bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
for (i=1; i<=mn; i++)
BIMATELEM(*bim,i,i)=n_Add(BIMATELEM(*a,i,i), bb, basecoeffs);
n_Delete(&bb,basecoeffs);
return bim;
}
bigintmat * bimSub(bigintmat * a, bigintmat * b)
{
if (a->cols() != b->cols()) return NULL;
if (a->rows() != b->rows()) return NULL;
if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
const coeffs basecoeffs = a->basecoeffs();
int i;
bigintmat * bim = new bigintmat(a->rows(), a->cols(), basecoeffs);
for (i=a->rows()*a->cols()-1;i>=0; i--)
bim->rawset(i, n_Sub((*a)[i], (*b)[i], basecoeffs), basecoeffs);
return bim;
}
bigintmat * bimSub(bigintmat * a, long b)
{
const int mn = si_min(a->rows(),a->cols());
const coeffs basecoeffs = a->basecoeffs();
number bb=n_Init(b,basecoeffs);
int i;
bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
for (i=1; i<=mn; i++)
BIMATELEM(*bim,i,i)=n_Sub(BIMATELEM(*a,i,i), bb, basecoeffs);
n_Delete(&bb,basecoeffs);
return bim;
}
//TODO: make special versions for certain rings!
bigintmat * bimMult(bigintmat * a, bigintmat * b)
{
const int ca = a->cols();
const int cb = b->cols();
const int ra = a->rows();
const int rb = b->rows();
if (ca != rb)
{
#ifndef SING_NDEBUG
Werror("wrong bigintmat sizes at multiplication a * b: acols: %d != brows: %d\n", ca, rb);
#endif
return NULL;
}
assume (ca == rb);
if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
const coeffs basecoeffs = a->basecoeffs();
int i, j, k;
number sum;
bigintmat * bim = new bigintmat(ra, cb, basecoeffs);
for (i=1; i<=ra; i++)
for (j=1; j<=cb; j++)
{
sum = n_Init(0, basecoeffs);
for (k=1; k<=ca; k++)
{
number prod = n_Mult( BIMATELEM(*a, i, k), BIMATELEM(*b, k, j), basecoeffs);
n_InpAdd(sum, prod, basecoeffs);
n_Delete(&prod, basecoeffs);
}
bim->rawset(i, j, sum, basecoeffs);
}
return bim;
}
bigintmat * bimMult(bigintmat * a, long b)
{
const int mn = a->rows()*a->cols();
const coeffs basecoeffs = a->basecoeffs();
number bb=n_Init(b,basecoeffs);
int i;
bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
for (i=0; i<mn; i++)
bim->rawset(i, n_Mult((*a)[i], bb, basecoeffs), basecoeffs);
n_Delete(&bb,basecoeffs);
return bim;
}
bigintmat * bimMult(bigintmat * a, number b, const coeffs cf)
{
if (cf!=a->basecoeffs()) return NULL;
const int mn = a->rows()*a->cols();
const coeffs basecoeffs = a->basecoeffs();
int i;
bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
for (i=0; i<mn; i++)
bim->rawset(i, n_Mult((*a)[i], b, basecoeffs), basecoeffs);
return bim;
}
// ----------------------------------------------------------------- //
// Korrekt?
intvec * bim2iv(bigintmat * b)
{
intvec * iv = new intvec(b->rows(), b->cols(), 0);
for (int i=0; i<(b->rows())*(b->cols()); i++)
(*iv)[i] = n_Int((*b)[i], b->basecoeffs()); // Geht das so?
return iv;
}
bigintmat * iv2bim(intvec * b, const coeffs C)
{
const int l = (b->rows())*(b->cols());
bigintmat * bim = new bigintmat(b->rows(), b->cols(), C);
for (int i=0; i < l; i++)
bim->rawset(i, n_Init((*b)[i], C), C);
return bim;
}
// ----------------------------------------------------------------- //
int bigintmat::compare(const bigintmat* op) const
{
assume (basecoeffs() == op->basecoeffs() );
#ifndef SING_NDEBUG
if (basecoeffs() != op->basecoeffs() )
WerrorS("wrong bigintmat comparison: different basecoeffs!\n");
#endif
if ((col!=1) ||(op->cols()!=1))
{
if((col!=op->cols())
|| (row!=op->rows()))
return -2;
}
int i;
for (i=0; i<si_min(row*col,op->rows()*op->cols()); i++)
{
if ( n_Greater(v[i], (*op)[i], basecoeffs()) )
return 1;
else if (! n_Equal(v[i], (*op)[i], basecoeffs()))
return -1;
}
for (; i<row; i++)
{
if ( n_GreaterZero(v[i], basecoeffs()) )
return 1;
else if (! n_IsZero(v[i], basecoeffs()) )
return -1;
}
for (; i<op->rows(); i++)
{
if ( n_GreaterZero((*op)[i], basecoeffs()) )
return -1;
else if (! n_IsZero((*op)[i], basecoeffs()) )
return 1;
}
return 0;
}
bigintmat * bimCopy(const bigintmat * b)
{
if (b == NULL)
return NULL;
return new bigintmat(b);
}
void bigintmat::Write()
{
int n = cols(), m=rows();
for(int i=1; i<= m; i++)
{
for(int j=1; j< n; j++)
{
n_Write(v[(i-1)*n+j-1], basecoeffs());
StringAppendS(", ");
}
if (n) n_Write(v[i*n-1], basecoeffs());
if (i<m)
{
StringAppendS(", ");
}
}
}
char* bigintmat::String()
{
StringSetS("");
Write();
return StringEndS();
}
void bigintmat::Print()
{
char * s = String();
PrintS(s);
omFree(s);
}
char* bigintmat::StringAsPrinted()
{
if ((col==0) || (row==0))
return NULL;
else
{
int * colwid = getwid(80);
if (colwid == NULL)
{
WerrorS("not enough space to print bigintmat");
WerrorS("try string(...) for a unformatted output");
return NULL;
}
char * ps;
int slength = 0;
for (int j=0; j<col; j++)
slength += colwid[j]*row;
slength += col*row+row;
ps = (char*) omAlloc0(sizeof(char)*(slength));
int pos = 0;
for (int i=0; i<col*row; i++)
{
StringSetS("");
n_Write(v[i], basecoeffs());
char * ts = StringEndS();
const int _nl = strlen(ts);
int cj = i%col;
if (_nl > colwid[cj])
{
StringSetS("");
int ci = i/col;
StringAppend("[%d,%d]", ci+1, cj+1);
char * ph = StringEndS();
int phl = strlen(ph);
if (phl > colwid[cj])
{
for (int j=0; j<colwid[cj]-1; j++)
ps[pos+j] = ' ';
ps[pos+colwid[cj]-1] = '*';
}
else
{
for (int j=0; j<colwid[cj]-phl; j++)
ps[pos+j] = ' ';
for (int j=0; j<phl; j++)
ps[pos+colwid[cj]-phl+j] = ph[j];
}
omFree(ph);
}
else // Mit Leerzeichen auffüllen und zahl reinschreiben
{
for (int j=0; j<(colwid[cj]-_nl); j++)
ps[pos+j] = ' ';
for (int j=0; j<_nl; j++)
ps[pos+colwid[cj]-_nl+j] = ts[j];
}
// ", " und (evtl) "\n" einfügen
if ((i+1)%col == 0)
{
if (i != col*row-1)
{
ps[pos+colwid[cj]] = ',';
ps[pos+colwid[cj]+1] = '\n';
pos += colwid[cj]+2;
}
}
else
{
ps[pos+colwid[cj]] = ',';
pos += colwid[cj]+1;
}
omFree(ts); // Hier ts zerstören
}
return(ps);
// omFree(ps);
}
}
static int intArrSum(int * a, int length)
{
int sum = 0;
for (int i=0; i<length; i++)
sum += a[i];
return sum;
}
static int findLongest(int * a, int length)
{
int l = 0;
int index;
for (int i=0; i<length; i++)
{
if (a[i] > l)
{
l = a[i];
index = i;
}
}
return index;
}
static int getShorter (int * a, int l, int j, int cols, int rows)
{
int sndlong = 0;
int min;
for (int i=0; i<rows; i++)
{
int index = cols*i+j;
if ((a[index] > sndlong) && (a[index] < l))
{
min = floor(log10((double)cols))+floor(log10((double)rows))+5;
if ((a[index] < min) && (min < l))
sndlong = min;
else
sndlong = a[index];
}
}
if (sndlong == 0)
{
min = floor(log10((double)cols))+floor(log10((double)rows))+5;
if (min < l)
sndlong = min;
else
sndlong = 1;
}
return sndlong;
}
int * bigintmat::getwid(int maxwid)
{
int const c = /*2**/(col-1)+1;
int * wv = (int*)omAlloc(sizeof(int)*col*row);
int * cwv = (int*)omAlloc(sizeof(int)*col);
for (int j=0; j<col; j++)
{
cwv[j] = 0;
for (int i=0; i<row; i++)
{
StringSetS("");
n_Write(v[col*i+j], basecoeffs());
char * tmp = StringEndS();
const int _nl = strlen(tmp);
wv[col*i+j] = _nl;
if (_nl > cwv[j]) cwv[j]=_nl;
omFree(tmp);
}
}
// Groesse verkleinern, bis < maxwid
if (intArrSum(cwv, col)+c > maxwid)
{
int j = findLongest(cwv, col);
cwv[j] = getShorter(wv, cwv[j], j, col, row);
}
omFree(wv);
return cwv;
}
void bigintmat::pprint(int maxwid)
{
if ((col==0) || (row==0))
PrintS("");
else
{
int * colwid = getwid(maxwid);
char * ps;
int slength = 0;
for (int j=0; j<col; j++)
slength += colwid[j]*row;
slength += col*row+row;
ps = (char*) omAlloc0(sizeof(char)*(slength));
int pos = 0;
for (int i=0; i<col*row; i++)
{
StringSetS("");
n_Write(v[i], basecoeffs());
char * ts = StringEndS();
const int _nl = strlen(ts);
int cj = i%col;
if (_nl > colwid[cj])
{
StringSetS("");
int ci = i/col;
StringAppend("[%d,%d]", ci+1, cj+1);
char * ph = StringEndS();
int phl = strlen(ph);
if (phl > colwid[cj])
{
for (int j=0; j<colwid[cj]-1; j++)
ps[pos+j] = ' ';
ps[pos+colwid[cj]-1] = '*';
}
else
{
for (int j=0; j<colwid[cj]-phl; j++)
ps[pos+j] = ' ';
for (int j=0; j<phl; j++)
ps[pos+colwid[cj]-phl+j] = ph[j];
}
omFree(ph);
}
else // Mit Leerzeichen auffüllen und zahl reinschreiben
{
for (int j=0; j<colwid[cj]-_nl; j++)
ps[pos+j] = ' ';
for (int j=0; j<_nl; j++)
ps[pos+colwid[cj]-_nl+j] = ts[j];
}
// ", " und (evtl) "\n" einfügen
if ((i+1)%col == 0)
{
if (i != col*row-1)
{
ps[pos+colwid[cj]] = ',';
ps[pos+colwid[cj]+1] = '\n';
pos += colwid[cj]+2;
}
}
else
{
ps[pos+colwid[cj]] = ',';
pos += colwid[cj]+1;
}
omFree(ts); // Hier ts zerstören
}
PrintS(ps);
omFree(ps);
omFree(colwid);
}
}
//swaps columns i and j
void bigintmat::swap(int i, int j)
{
if ((i <= col) && (j <= col) && (i>0) && (j>0))
{
number tmp;
number t;
for (int k=1; k<=row; k++)
{
tmp = get(k, i);
t = view(k, j);
set(k, i, t);
set(k, j, tmp);
n_Delete(&tmp, basecoeffs());
}
}
else
WerrorS("Error in swap");
}
void bigintmat::swaprow(int i, int j)
{
if ((i <= row) && (j <= row) && (i>0) && (j>0))
{
number tmp;
number t;
for (int k=1; k<=col; k++)
{
tmp = get(i, k);
t = view(j, k);
set(i, k, t);
set(j, k, tmp);
n_Delete(&tmp, basecoeffs());
}
}
else
WerrorS("Error in swaprow");
}
int bigintmat::findnonzero(int i)
{
for (int j=1; j<=col; j++)
{
if (!n_IsZero(view(i,j), basecoeffs()))
{
return j;
}
}
return 0;
}
int bigintmat::findcolnonzero(int j)
{
for (int i=row; i>=1; i--)
{
if (!n_IsZero(view(i,j), basecoeffs()))
{
return i;
}
}
return 0;
}
void bigintmat::getcol(int j, bigintmat *a)
{
assume((j<=col) && (j>=1));
if (((a->rows() != row) || (a->cols() != 1)) && ((a->rows() != 1) || (a->cols() != row)))
{
assume(0);
WerrorS("Error in getcol. Dimensions must agree!");
return;
}
if (!nCoeffs_are_equal(basecoeffs(), a->basecoeffs()))
{
nMapFunc f = n_SetMap(basecoeffs(), a->basecoeffs());
number t1, t2;
for (int i=1; i<=row;i++)
{
t1 = get(i,j);
t2 = f(t1, basecoeffs(), a->basecoeffs());
a->set(i-1,t1);
n_Delete(&t1, basecoeffs());
n_Delete(&t2, a->basecoeffs());
}
return;
}
number t1;
for (int i=1; i<=row;i++)
{
t1 = view(i,j);
a->set(i-1,t1);
}
}
void bigintmat::getColRange(int j, int no, bigintmat *a)
{
number t1;
for(int ii=0; ii< no; ii++)
{
for (int i=1; i<=row;i++)
{
t1 = view(i, ii+j);
a->set(i, ii+1, t1);
}
}
}
void bigintmat::getrow(int i, bigintmat *a)
{
if ((i>row) || (i<1))
{
WerrorS("Error in getrow: Index out of range!");
return;
}
if (((a->rows() != 1) || (a->cols() != col)) && ((a->rows() != col) || (a->cols() != 1)))
{
WerrorS("Error in getrow. Dimensions must agree!");
return;
}
if (!nCoeffs_are_equal(basecoeffs(), a->basecoeffs()))
{
nMapFunc f = n_SetMap(basecoeffs(), a->basecoeffs());
number t1, t2;
for (int j=1; j<=col;j++)
{
t1 = get(i,j);
t2 = f(t1, basecoeffs(), a->basecoeffs());
a->set(j-1,t2);
n_Delete(&t1, basecoeffs());
n_Delete(&t2, a->basecoeffs());
}
return;
}
number t1;
for (int j=1; j<=col;j++)
{
t1 = get(i,j);
a->set(j-1,t1);
n_Delete(&t1, basecoeffs());
}
}
void bigintmat::setcol(int j, bigintmat *m)
{
if ((j>col) || (j<1))
{
WerrorS("Error in setcol: Index out of range!");
return;
}
if (((m->rows() != row) || (m->cols() != 1)) && ((m->rows() != 1) || (m->cols() != row)))
{
WerrorS("Error in setcol. Dimensions must agree!");
return;
}
if (!nCoeffs_are_equal(basecoeffs(), m->basecoeffs()))
{
nMapFunc f = n_SetMap(m->basecoeffs(), basecoeffs());
number t1,t2;
for (int i=1; i<=row; i++)
{
t1 = m->get(i-1);
t2 = f(t1, m->basecoeffs(),basecoeffs());
set(i, j, t2);
n_Delete(&t2, basecoeffs());
n_Delete(&t1, m->basecoeffs());
}
return;
}
number t1;
for (int i=1; i<=row; i++)
{
t1 = m->view(i-1);
set(i, j, t1);
}
}
void bigintmat::setrow(int j, bigintmat *m)
{
if ((j>row) || (j<1))
{
WerrorS("Error in setrow: Index out of range!");
return;
}
if (((m->rows() != 1) || (m->cols() != col)) && ((m->rows() != col) || (m->cols() != 1)))
{
WerrorS("Error in setrow. Dimensions must agree!");
return;
}
if (!nCoeffs_are_equal(basecoeffs(), m->basecoeffs()))
{
nMapFunc f = n_SetMap(m->basecoeffs(), basecoeffs());
number tmp1,tmp2;
for (int i=1; i<=col; i++)
{
tmp1 = m->get(i-1);
tmp2 = f(tmp1, m->basecoeffs(),basecoeffs());
set(j, i, tmp2);
n_Delete(&tmp2, basecoeffs());
n_Delete(&tmp1, m->basecoeffs());
}
return;
}
number tmp;
for (int i=1; i<=col; i++)
{
tmp = m->view(i-1);
set(j, i, tmp);
}
}
bool bigintmat::add(bigintmat *b)
{
if ((b->rows() != row) || (b->cols() != col))
{
WerrorS("Error in bigintmat::add. Dimensions do not agree!");
return false;
}
if (!nCoeffs_are_equal(basecoeffs(), b->basecoeffs()))
{
WerrorS("Error in bigintmat::add. coeffs do not agree!");
return false;
}
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
rawset(i, j, n_Add(b->view(i,j), view(i,j), basecoeffs()));
}
}
return true;
}
bool bigintmat::sub(bigintmat *b)
{
if ((b->rows() != row) || (b->cols() != col))
{
WerrorS("Error in bigintmat::sub. Dimensions do not agree!");
return false;
}
if (!nCoeffs_are_equal(basecoeffs(), b->basecoeffs()))
{
WerrorS("Error in bigintmat::sub. coeffs do not agree!");
return false;
}
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
rawset(i, j, n_Sub(view(i,j), b->view(i,j), basecoeffs()));
}
}
return true;
}
bool bigintmat::skalmult(number b, coeffs c)
{
if (!nCoeffs_are_equal(c, basecoeffs()))
{
WerrorS("Wrong coeffs\n");
return false;
}
number t1, t2;
if ( n_IsOne(b,c)) return true;
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
t1 = view(i, j);
t2 = n_Mult(t1, b, basecoeffs());
rawset(i, j, t2);
}
}
return true;
}
bool bigintmat::addcol(int i, int j, number a, coeffs c)
{
if ((i>col) || (j>col) || (i<1) || (j<1))
{
WerrorS("Error in addcol: Index out of range!");
return false;
}
if (!nCoeffs_are_equal(c, basecoeffs()))
{
WerrorS("Error in addcol: coeffs do not agree!");
return false;
}
number t1, t2, t3;
for (int k=1; k<=row; k++)
{
t1 = view(k, j);
t2 = view(k, i);
t3 = n_Mult(t1, a, basecoeffs());
n_InpAdd(t3, t2, basecoeffs());
rawset(k, i, t3);
}
return true;
}
bool bigintmat::addrow(int i, int j, number a, coeffs c)
{
if ((i>row) || (j>row) || (i<1) || (j<1))
{
WerrorS("Error in addrow: Index out of range!");
return false;
}
if (!nCoeffs_are_equal(c, basecoeffs()))
{
WerrorS("Error in addrow: coeffs do not agree!");
return false;
}
number t1, t2, t3;
for (int k=1; k<=col; k++)
{
t1 = view(j, k);
t2 = view(i, k);
t3 = n_Mult(t1, a, basecoeffs());
n_InpAdd(t3, t2, basecoeffs());
rawset(i, k, t3);
}
return true;
}
void bigintmat::colskalmult(int i, number a, coeffs c)
{
if ((i>=1) && (i<=col) && (nCoeffs_are_equal(c, basecoeffs())))
{
number t, tmult;
for (int j=1; j<=row; j++)
{
t = view(j, i);
tmult = n_Mult(a, t, basecoeffs());
rawset(j, i, tmult);
}
}
else
WerrorS("Error in colskalmult");
}
void bigintmat::rowskalmult(int i, number a, coeffs c)
{
if ((i>=1) && (i<=row) && (nCoeffs_are_equal(c, basecoeffs())))
{
number t, tmult;
for (int j=1; j<=col; j++)
{
t = view(i, j);
tmult = n_Mult(a, t, basecoeffs());
rawset(i, j, tmult);
}
}
else
WerrorS("Error in rowskalmult");
}
void bigintmat::concatrow(bigintmat *a, bigintmat *b)
{
int ay = a->cols();
int ax = a->rows();
int by = b->cols();
int bx = b->rows();
number tmp;
if (!((col == ay) && (col == by) && (ax+bx == row)))
{
WerrorS("Error in concatrow. Dimensions must agree!");
return;
}
if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs())))
{
WerrorS("Error in concatrow. coeffs do not agree!");
return;
}
for (int i=1; i<=ax; i++)
{
for (int j=1; j<=ay; j++)
{
tmp = a->get(i,j);
set(i, j, tmp);
n_Delete(&tmp, basecoeffs());
}
}
for (int i=1; i<=bx; i++)
{
for (int j=1; j<=by; j++)
{
tmp = b->get(i,j);
set(i+ax, j, tmp);
n_Delete(&tmp, basecoeffs());
}
}
}
void bigintmat::extendCols(int i)
{
bigintmat * tmp = new bigintmat(rows(), i, basecoeffs());
appendCol(tmp);
delete tmp;
}
void bigintmat::appendCol (bigintmat *a)
{
coeffs R = basecoeffs();
int ay = a->cols();
int ax = a->rows();
assume(row == ax);
assume(nCoeffs_are_equal(a->basecoeffs(), R));
bigintmat * tmp = new bigintmat(rows(), cols() + ay, R);
tmp->concatcol(this, a);
this->swapMatrix(tmp);
delete tmp;
}
void bigintmat::concatcol (bigintmat *a, bigintmat *b) {
int ay = a->cols();
int ax = a->rows();
int by = b->cols();
int bx = b->rows();
number tmp;
assume(row==ax && row == bx && ay+by ==col);
assume(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs()));
for (int i=1; i<=ax; i++)
{
for (int j=1; j<=ay; j++)
{
tmp = a->view(i,j);
set(i, j, tmp);
}
}
for (int i=1; i<=bx; i++)
{
for (int j=1; j<=by; j++)
{
tmp = b->view(i,j);
set(i, j+ay, tmp);
}
}
}
void bigintmat::splitrow(bigintmat *a, bigintmat *b)
{
int ay = a->cols();
int ax = a->rows();
int by = b->cols();
int bx = b->rows();
number tmp;
if (!(ax + bx == row))
{
WerrorS("Error in splitrow. Dimensions must agree!");
}
else if (!((col == ay) && (col == by)))
{
WerrorS("Error in splitrow. Dimensions must agree!");
}
else if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs())))
{
WerrorS("Error in splitrow. coeffs do not agree!");
}
else
{
for(int i = 1; i<=ax; i++)
{
for(int j = 1; j<=ay;j++)
{
tmp = get(i,j);
a->set(i,j,tmp);
n_Delete(&tmp, basecoeffs());
}
}
for (int i =1; i<=bx; i++)
{
for (int j=1;j<=col;j++)
{
tmp = get(i+ax, j);
b->set(i,j,tmp);
n_Delete(&tmp, basecoeffs());
}
}
}
}
void bigintmat::splitcol(bigintmat *a, bigintmat *b)
{
int ay = a->cols();
int ax = a->rows();
int by = b->cols();
int bx = b->rows();
number tmp;
if (!((row == ax) && (row == bx)))
{
WerrorS("Error in splitcol. Dimensions must agree!");
}
else if (!(ay+by == col))
{
WerrorS("Error in splitcol. Dimensions must agree!");
}
else if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs())))
{
WerrorS("Error in splitcol. coeffs do not agree!");
}
else
{
for (int i=1; i<=ax; i++)
{
for (int j=1; j<=ay; j++)
{
tmp = view(i,j);
a->set(i,j,tmp);
}
}
for (int i=1; i<=bx; i++)
{
for (int j=1; j<=by; j++)
{
tmp = view(i,j+ay);
b->set(i,j,tmp);
}
}
}
}
void bigintmat::splitcol(bigintmat *a, int i)
{
number tmp;
if ((a->rows() != row) || (a->cols()+i-1 > col) || (i<1))
{
WerrorS("Error in splitcol. Dimensions must agree!");
return;
}
if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs())))
{
WerrorS("Error in splitcol. coeffs do not agree!");
return;
}
int width = a->cols();
for (int j=1; j<=width; j++)
{
for (int k=1; k<=row; k++)
{
tmp = get(k, j+i-1);
a->set(k, j, tmp);
n_Delete(&tmp, basecoeffs());
}
}
}
void bigintmat::splitrow(bigintmat *a, int i)
{
number tmp;
if ((a->cols() != col) || (a->rows()+i-1 > row) || (i<1))
{
WerrorS("Error in Marco-splitrow");
return;
}
if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs())))
{
WerrorS("Error in splitrow. coeffs do not agree!");
return;
}
int height = a->rows();
for (int j=1; j<=height; j++)
{
for (int k=1; k<=col; k++)
{
tmp = view(j+i-1, k);
a->set(j, k, tmp);
}
}
}
bool bigintmat::copy(bigintmat *b)
{
if ((b->rows() != row) || (b->cols() != col))
{
WerrorS("Error in bigintmat::copy. Dimensions do not agree!");
return false;
}
if (!nCoeffs_are_equal(basecoeffs(), b->basecoeffs()))
{
WerrorS("Error in bigintmat::copy. coeffs do not agree!");
return false;
}
number t1;
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
t1 = b->view(i, j);
set(i, j, t1);
}
}
return true;
}
/// copy the submatrix of b, staring at (a,b) having n rows, m cols into
/// the given matrix at pos. (c,d)
/// needs c+n, d+m <= rows, cols
/// a+n, b+m <= b.rows(), b.cols()
void bigintmat::copySubmatInto(bigintmat *B, int a, int b, int n, int m, int c, int d)
{
number t1;
for (int i=1; i<=n; i++)
{
for (int j=1; j<=m; j++)
{
t1 = B->view(a+i-1, b+j-1);
set(c+i-1, d+j-1, t1);
}
}
}
int bigintmat::isOne()
{
coeffs r = basecoeffs();
if (row==col)
{
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
if (i==j)
{
if (!n_IsOne(view(i, j), r))
return 0;
}
else
{
if (!n_IsZero(view(i,j), r))
return 0;
}
}
}
}
return 1;
}
void bigintmat::one()
{
if (row==col)
{
number one = n_Init(1, basecoeffs()),
zero = n_Init(0, basecoeffs());
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
if (i==j)
{
set(i, j, one);
}
else
{
set(i, j, zero);
}
}
}
n_Delete(&one, basecoeffs());
n_Delete(&zero, basecoeffs());
}
}
void bigintmat::zero()
{
number tmp = n_Init(0, basecoeffs());
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
set(i, j, tmp);
}
}
n_Delete(&tmp,basecoeffs());
}
int bigintmat::isZero()
{
for (int i=1; i<=row; i++) {
for (int j=1; j<=col; j++) {
if (!n_IsZero(view(i,j), basecoeffs()))
return FALSE;
}
}
return TRUE;
}
//****************************************************************************
//
//****************************************************************************
//used in the det function. No idea what it does.
//looks like it return the submatrix where the i-th row
//and j-th column has been removed in the LaPlace generic
//determinant algorithm
bigintmat *bigintmat::elim(int i, int j)
{
if ((i<=0) || (i>row) || (j<=0) || (j>col))
return NULL;
int cx, cy;
cx=1;
cy=1;
number t;
bigintmat *b = new bigintmat(row-1, col-1, basecoeffs());
for (int k=1; k<=row; k++) {
if (k!=i)
{
cy=1;
for (int l=1; l<=col; l++)
{
if (l!=j)
{
t = get(k, l);
b->set(cx, cy, t);
n_Delete(&t, basecoeffs());
cy++;
}
}
cx++;
}
}
return b;
}
//returns d such that a/d is the inverse of the input
//TODO: make work for p not prime using the euc stuff.
//long term: rewrite for Z using p-adic lifting
//and Dixon. Possibly even the sparse recent Storjohann stuff
number bigintmat::pseudoinv(bigintmat *a) {
// Falls Matrix über reellen Zahlen nicht invertierbar, breche ab
assume((a->rows() == row) && (a->rows() == a->cols()) && (row == col));
number det = this->det(); //computes the HNF, so should e reused.
if ((n_IsZero(det, basecoeffs())))
return det;
// Hänge Einheitsmatrix über Matrix und wendet HNF an. An Stelle der Einheitsmatrix steht im Ergebnis die Transformationsmatrix dazu
a->one();
bigintmat *m = new bigintmat(2*row, col, basecoeffs());
m->concatrow(a,this);
m->hnf();
// Arbeite weiterhin mit der zusammengehängten Matrix
// Laufe durch die Diagonalelemente, und multipliziere jede Spalte rechts davon damit, speichere aber den alten Eintrag der Spalte, temp, der in der Zeile des Diagonalelements liegt, zwischen. Dann addiere das -temp-Fache der Diagonalspalte zur entsprechenenden Spalte rechts davon. Dadurch entsteht überall rechts der Diagonalen eine 0
number diag;
number temp, ttemp;
for (int i=1; i<=col; i++) {
diag = m->get(row+i, i);
for (int j=i+1; j<=col; j++) {
temp = m->get(row+i, j);
m->colskalmult(j, diag, basecoeffs());
temp = n_InpNeg(temp, basecoeffs());
m->addcol(j, i, temp, basecoeffs());
n_Delete(&temp, basecoeffs());
}
n_Delete(&diag, basecoeffs());
}
// Falls wir nicht modulo n arbeiten, können wir die Spalten durch den ggT teilen, um die Einträge kleiner zu bekommen
// Bei Z/n sparen wir uns das, da es hier sinnlos ist
number g;
number gcd;
for (int j=1; j<=col; j++) {
g = n_Init(0, basecoeffs());
for (int i=1; i<=2*row; i++) {
temp = m->get(i,j);
gcd = n_Gcd(g, temp, basecoeffs());
n_Delete(&g, basecoeffs());
n_Delete(&temp, basecoeffs());
g = n_Copy(gcd, basecoeffs());
n_Delete(&gcd, basecoeffs());
}
if (!(n_IsOne(g, basecoeffs())))
m->colskaldiv(j, g);
n_Delete(&g, basecoeffs());
}
// Nun müssen die Diagonalelemente durch Spaltenmultiplikation gleich gesett werden. Bei Z können wir mit dem kgV arbeiten, bei Z/n bringen wir jedes Diagonalelement auf 1 (wir arbeiten immer mit n = Primzahl. Für n != Primzahl muss noch an anderen Stellen etwas geändert werden)
g = n_Init(0, basecoeffs());
number prod = n_Init(1, basecoeffs());
for (int i=1; i<=col; i++) {
gcd = n_Gcd(g, m->get(row+i, i), basecoeffs());
n_Delete(&g, basecoeffs());
g = n_Copy(gcd, basecoeffs());
n_Delete(&gcd, basecoeffs());
ttemp = n_Copy(prod, basecoeffs());
temp = m->get(row+i, i);
n_Delete(&prod, basecoeffs());
prod = n_Mult(ttemp, temp, basecoeffs());
n_Delete(&ttemp, basecoeffs());
n_Delete(&temp, basecoeffs());
}
number lcm;
lcm = n_Div(prod, g, basecoeffs());
for (int j=1; j<=col; j++) {
ttemp = m->get(row+j,j);
temp = n_QuotRem(lcm, ttemp, NULL, basecoeffs());
m->colskalmult(j, temp, basecoeffs());
n_Delete(&ttemp, basecoeffs());
n_Delete(&temp, basecoeffs());
}
n_Delete(&lcm, basecoeffs());
n_Delete(&prod, basecoeffs());
number divisor = m->get(row+1, 1);
m->splitrow(a, 1);
delete m;
n_Delete(&det, basecoeffs());
return divisor;
}
number bigintmat::trace()
{
assume (col == row);
number t = get(1,1),
h;
coeffs r = basecoeffs();
for(int i=2; i<= col; i++) {
h = n_Add(t, view(i,i), r);
n_Delete(&t, r);
t = h;
}
return t;
}
number bigintmat::det()
{
assume (row==col);
if (col == 1)
return get(1, 1);
// should work as well in Z/pZ of type n_Zp?
// relies on XExtGcd and the other euc. functions.
if ( getCoeffType(basecoeffs())== n_Z || getCoeffType(basecoeffs() )== n_Zn) {
return hnfdet();
}
number sum = n_Init(0, basecoeffs());
number t1, t2, t3, t4;
bigintmat *b;
for (int i=1; i<=row; i++) {
b = elim(i, 1);
t1 = get(i, 1);
t2 = b->det();
t3 = n_Mult(t1, t2, basecoeffs());
t4 = n_Copy(sum, basecoeffs());
n_Delete(&sum, basecoeffs());
if ((i+1)>>1<<1==(i+1))
sum = n_Add(t4, t3, basecoeffs());
else
sum = n_Sub(t4, t3, basecoeffs());
n_Delete(&t1, basecoeffs());
n_Delete(&t2, basecoeffs());
n_Delete(&t3, basecoeffs());
n_Delete(&t4, basecoeffs());
}
return sum;
}
number bigintmat::hnfdet()
{
assume (col == row);
if (col == 1)
return get(1, 1);
bigintmat *m = new bigintmat(this);
m->hnf();
number prod = n_Init(1, basecoeffs());
number temp, temp2;
for (int i=1; i<=col; i++) {
temp = m->get(i, i);
temp2 = n_Mult(temp, prod, basecoeffs());
n_Delete(&prod, basecoeffs());
prod = temp2;
n_Delete(&temp, basecoeffs());
}
delete m;
return prod;
}
void bigintmat::swapMatrix(bigintmat *a)
{
int n = rows(), m = cols();
row = a->rows();
col = a->cols();
number * V = v;
v = a->v;
a->v = V;
a->row = n;
a->col = m;
}
int bigintmat::colIsZero(int j)
{
coeffs R = basecoeffs();
for(int i=1; i<=rows(); i++)
if (!n_IsZero(view(i, j), R)) return FALSE;
return TRUE;
}
void bigintmat::howell()
{
coeffs R = basecoeffs();
hnf(); // as a starting point...
if (getCoeffType(R)== n_Z) return; //wrong, need to prune!
int n = cols(), m = rows(), i, j, k;
//make sure, the matrix has enough space. We need no rows+1 columns.
//The resulting Howell form will be pruned to be at most square.
bigintmat * t = new bigintmat(m, m+1, R);
t->copySubmatInto(this, 1, n>m ? n-m+1 : 1, m, n>m ? m : n, 1, n>m ? 2 : m+2-n );
swapMatrix(t);
delete t;
for(i=1; i<= cols(); i++) {
if (!colIsZero(i)) break;
}
assume (i>1);
if (i>cols()) {
t = new bigintmat(rows(), 0, R);
swapMatrix(t);
delete t;
return; // zero matrix found, clearly normal.
}
int last_zero_col = i-1;
for (int c = cols(); c>0; c--) {
for(i=rows(); i>0; i--) {
if (!n_IsZero(view(i, c), R)) break;
}
if (i==0) break; // matrix SHOULD be zero from here on
number a = n_Ann(view(i, c), R);
addcol(last_zero_col, c, a, R);
n_Delete(&a, R);
for(j = c-1; j>last_zero_col; j--) {
for(k=rows(); k>0; k--) {
if (!n_IsZero(view(k, j), R)) break;
if (!n_IsZero(view(k, last_zero_col), R)) break;
}
if (k==0) break;
if (!n_IsZero(view(k, last_zero_col), R)) {
number gcd, co1, co2, co3, co4;
gcd = n_XExtGcd(view(k, last_zero_col), view(k, j), &co1, &co2, &co3, &co4, R);
if (n_Equal(gcd, view(k, j), R)) {
number q = n_Div(view(k, last_zero_col), gcd, R);
q = n_InpNeg(q, R);
addcol(last_zero_col, j, q, R);
n_Delete(&q, R);
} else if (n_Equal(gcd, view(k, last_zero_col), R)) {
swap(last_zero_col, k);
number q = n_Div(view(k, last_zero_col), gcd, R);
q = n_InpNeg(q, R);
addcol(last_zero_col, j, q, R);
n_Delete(&q, R);
} else {
coltransform(last_zero_col, j, co3, co4, co1, co2);
}
n_Delete(&gcd, R);
n_Delete(&co1, R);
n_Delete(&co2, R);
n_Delete(&co3, R);
n_Delete(&co4, R);
}
}
for(k=rows(); k>0; k--) {
if (!n_IsZero(view(k, last_zero_col), R)) break;
}
if (k) last_zero_col--;
}
t = new bigintmat(rows(), cols()-last_zero_col, R);
t->copySubmatInto(this, 1, last_zero_col+1, rows(), cols()-last_zero_col, 1, 1);
swapMatrix(t);
delete t;
}
void bigintmat::hnf()
{
// Laufen von unten nach oben und von links nach rechts
// CF: TODO: for n_Z: write a recursive version. This one will
// have exponential blow-up. Look at Michianchio
// Alternatively, do p-adic det and modular method
#if 0
char * s;
::PrintS("mat over Z is \n");
::Print("%s\n", s = nCoeffString(basecoeffs()));
omFree(s);
Print();
::Print("\n(%d x %d)\n", rows(), cols());
#endif
int i = rows();
int j = cols();
number q = n_Init(0, basecoeffs());
number one = n_Init(1, basecoeffs());
number minusone = n_Init(-1, basecoeffs());
number tmp1 = n_Init(0, basecoeffs());
number tmp2 = n_Init(0, basecoeffs());
number co1, co2, co3, co4;
number ggt = n_Init(0, basecoeffs());
while ((i>0) && (j>0))
{
// Falls erstes Nicht-Null-Element in Zeile i nicht existiert, oder hinter Spalte j vorkommt, gehe in nächste Zeile
if ((findnonzero(i)==0) || (findnonzero(i)>j))
{
i--;
}
else
{
// Laufe von links nach rechts durch die Zeile:
for (int l=1; l<=j-1; l++)
{
n_Delete(&tmp1, basecoeffs());
tmp1 = get(i, l);
// Falls Eintrag (im folgenden x genannt) gleich 0, gehe eine Spalte weiter. Ansonsten...
if (!n_IsZero(tmp1, basecoeffs()))
{
n_Delete(&tmp2, basecoeffs());
tmp2 = get(i, l+1);
// Falls Eintrag (i.f. y g.) rechts daneben gleich 0, tausche beide Spalten, sonst...
if (!n_IsZero(tmp2, basecoeffs()))
{
n_Delete(&ggt, basecoeffs());
ggt = n_XExtGcd(tmp1, tmp2, &co1, &co2, &co3, &co4, basecoeffs());
// Falls x=ggT(x, y), tausche die beiden Spalten und ziehe die (neue) rechte Spalte so häufig von der linken ab, dass an der ehemaligen Stelle von x nun eine 0 steht. Dazu:
if (n_Equal(tmp1, ggt, basecoeffs()))
{
swap(l, l+1);
n_Delete(&q, basecoeffs());
q = n_Div(tmp2, ggt, basecoeffs());
q = n_InpNeg(q, basecoeffs());
// Dann addiere das -q-fache der (neuen) rechten Spalte zur linken dazu. Damit erhalten wir die gewünschte 0
addcol(l, l+1, q, basecoeffs());
n_Delete(&q, basecoeffs());
}
else if (n_Equal(tmp1, minusone, basecoeffs()))
{
// Falls x=-1, so ist x=-ggt(x, y). Dann gehe wie oben vor, multipliziere aber zuerst die neue rechte Spalte (die mit x) mit -1
// Die Berechnung von q (=y/ggt) entfällt, da ggt=1
swap(l, l+1);
colskalmult(l+1, minusone, basecoeffs());
tmp2 = n_InpNeg(tmp2, basecoeffs());
addcol(l, l+1, tmp2, basecoeffs());
}
else
{
// CF: use the 2x2 matrix (co1, co2)(co3, co4) to
// get the gcd in position and the 0 in the other:
#ifdef CF_DEB
::PrintS("applying trafo\n");
StringSetS("");
n_Write(co1, basecoeffs()); StringAppendS("\t");
n_Write(co2, basecoeffs()); StringAppendS("\t");
n_Write(co3, basecoeffs()); StringAppendS("\t");
n_Write(co4, basecoeffs()); StringAppendS("\t");
::Print("%s\nfor l=%d\n", StringEndS(), l);
{char * s = String();
::Print("to %s\n", s);omFree(s);};
#endif
coltransform(l, l+1, co3, co4, co1, co2);
#ifdef CF_DEB
{char * s = String();
::Print("gives %s\n", s);}
#endif
}
n_Delete(&co1, basecoeffs());
n_Delete(&co2, basecoeffs());
n_Delete(&co3, basecoeffs());
n_Delete(&co4, basecoeffs());
}
else
{
swap(l, l+1);
}
// Dann betrachte die vormals rechte Spalte als neue linke, und die rechts daneben als neue rechte.
}
}
// normalize by units:
if (!n_IsZero(view(i, j), basecoeffs()))
{
number u = n_GetUnit(view(i, j), basecoeffs());
if (!n_IsOne(u, basecoeffs()))
{
colskaldiv(j, u);
}
n_Delete(&u, basecoeffs());
}
// Zum Schluss mache alle Einträge rechts vom Diagonalelement betragsmäßig kleiner als dieses
for (int l=j+1; l<=col; l++)
{
n_Delete(&q, basecoeffs());
q = n_QuotRem(view(i, l), view(i, j), NULL, basecoeffs());
q = n_InpNeg(q, basecoeffs());
addcol(l, j, q, basecoeffs());
}
i--;
j--;
// Dann betrachte die Zeile darüber und gehe dort wie vorher vor
}
}
n_Delete(&q, basecoeffs());
n_Delete(&tmp1, basecoeffs());
n_Delete(&tmp2, basecoeffs());
n_Delete(&ggt, basecoeffs());
n_Delete(&one, basecoeffs());
n_Delete(&minusone, basecoeffs());
#if 0
::PrintS("hnf over Z is \n");
Print();
::Print("\n(%d x %d)\n", rows(), cols());
#endif
}
bigintmat * bimChangeCoeff(bigintmat *a, coeffs cnew)
{
coeffs cold = a->basecoeffs();
bigintmat *b = new bigintmat(a->rows(), a->cols(), cnew);
// Erzeugt Karte von alten coeffs nach neuen
nMapFunc f = n_SetMap(cold, cnew);
number t1;
number t2;
// apply map to all entries.
for (int i=1; i<=a->rows(); i++)
{
for (int j=1; j<=a->cols(); j++)
{
t1 = a->get(i, j);
t2 = f(t1, cold, cnew);
b->set(i, j, t2);
n_Delete(&t1, cold);
n_Delete(&t2, cnew);
}
}
return b;
}
//OK: a HNF of (this | p*I)
//so the result will always have FULL rank!!!!
//(This is different form a lift of the HNF mod p: consider the matrix (p)
//to see the difference. It CAN be computed as HNF mod p^2 usually..)
bigintmat * bigintmat::modhnf(number p, coeffs R)
{
coeffs Rp = numbercoeffs(p, R); // R/pR
bigintmat *m = bimChangeCoeff(this, Rp);
m->howell();
bigintmat *a = bimChangeCoeff(m, R);
delete m;
bigintmat *C = new bigintmat(rows(), rows(), R);
int piv = rows(), i = a->cols();
while (piv)
{
if (!i || n_IsZero(a->view(piv, i), R))
{
C->set(piv, piv, p, R);
}
else
{
C->copySubmatInto(a, 1, i, rows(), 1, 1, piv);
i--;
}
piv--;
}
delete a;
return C;
}
//exactly divide matrix by b
void bigintmat::skaldiv(number b)
{
number tmp1, tmp2;
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
tmp1 = view(i, j);
tmp2 = n_Div(tmp1, b, basecoeffs());
rawset(i, j, tmp2);
}
}
}
//exactly divide col j by b
void bigintmat::colskaldiv(int j, number b)
{
number tmp1, tmp2;
for (int i=1; i<=row; i++)
{
tmp1 = view(i, j);
tmp2 = n_Div(tmp1, b, basecoeffs());
rawset(i, j, tmp2);
}
}
// col(j, k) <- col(j,k)*matrix((a, c)(b, d))
// mostly used internally in the hnf and Howell stuff
void bigintmat::coltransform(int j, int k, number a, number b, number c, number d)
{
number tmp1, tmp2, tmp3, tmp4;
for (int i=1; i<=row; i++)
{
tmp1 = get(i, j);
tmp2 = get(i, k);
tmp3 = n_Mult(tmp1, a, basecoeffs());
tmp4 = n_Mult(tmp2, b, basecoeffs());
n_InpAdd(tmp3, tmp4, basecoeffs());
n_Delete(&tmp4, basecoeffs());
n_InpMult(tmp1, c, basecoeffs());
n_InpMult(tmp2, d, basecoeffs());
n_InpAdd(tmp1, tmp2, basecoeffs());
n_Delete(&tmp2, basecoeffs());
set(i, j, tmp3);
set(i, k, tmp1);
n_Delete(&tmp1, basecoeffs());
n_Delete(&tmp3, basecoeffs());
}
}
//reduce all entries mod p. Does NOT change the coeffs type
void bigintmat::mod(number p)
{
// produce the matrix in Z/pZ
number tmp1, tmp2;
for (int i=1; i<=row; i++)
{
for (int j=1; j<=col; j++)
{
tmp1 = get(i, j);
tmp2 = n_IntMod(tmp1, p, basecoeffs());
n_Delete(&tmp1, basecoeffs());
set(i, j, tmp2);
}
}
}
void bimMult(bigintmat *a, bigintmat *b, bigintmat *c)
{
if (!nCoeffs_are_equal(a->basecoeffs(), b->basecoeffs()))
{
WerrorS("Error in bimMult. Coeffs do not agree!");
return;
}
if ((a->rows() != c->rows()) || (b->cols() != c->cols()) || (a->cols() != b->rows()))
{
WerrorS("Error in bimMult. Dimensions do not agree!");
return;
}
bigintmat *tmp = bimMult(a, b);
c->copy(tmp);
delete tmp;
}
static void reduce_mod_howell(bigintmat *A, bigintmat *b, bigintmat * eps, bigintmat *x)
{
//write b = Ax + eps where eps is "small" in the sense of bounded by the
//pivot entries in H. H does not need to be Howell (or HNF) but need
//to be triagonal in the same direction.
//b can have multiple columns.
#if 0
PrintS("reduce_mod_howell: A:\n");
A->Print();
PrintS("\nb:\n");
b->Print();
#endif
coeffs R = A->basecoeffs();
assume(x->basecoeffs() == R);
assume(b->basecoeffs() == R);
assume(eps->basecoeffs() == R);
if (!A->cols())
{
x->zero();
eps->copy(b);
#if 0
PrintS("\nx:\n");
x->Print();
PrintS("\neps:\n");
eps->Print();
PrintS("\n****************************************\n");
#endif
return;
}
bigintmat * B = new bigintmat(b->rows(), 1, R);
for(int i=1; i<= b->cols(); i++)
{
int A_col = A->cols();
b->getcol(i, B);
for(int j = B->rows(); j>0; j--)
{
number Ai = A->view(A->rows() - B->rows() + j, A_col);
if (n_IsZero(Ai, R) &&
n_IsZero(B->view(j, 1), R))
{
continue; //all is fine: 0*x = 0
}
else if (n_IsZero(B->view(j, 1), R))
{
x->rawset(x->rows() - B->rows() + j, i, n_Init(0, R));
A_col--;
}
else if (n_IsZero(Ai, R))
{
A_col--;
}
else
{
// "solve" ax=b, possibly enlarging d
number Bj = B->view(j, 1);
number q = n_Div(Bj, Ai, R);
x->rawset(x->rows() - B->rows() + j, i, q);
for(int k=j; k>B->rows() - A->rows(); k--)
{
//B[k] = B[k] - x[k]A[k][j]
number s = n_Mult(q, A->view(A->rows() - B->rows() + k, A_col), R);
B->rawset(k, 1, n_Sub(B->view(k, 1), s, R));
n_Delete(&s, R);
}
A_col--;
}
if (!A_col)
{
break;
}
}
eps->setcol(i, B);
}
delete B;
#if 0
PrintS("\nx:\n");
x->Print();
PrintS("\neps:\n");
eps->Print();
PrintS("\n****************************************\n");
#endif
}
static bigintmat * prependIdentity(bigintmat *A)
{
coeffs R = A->basecoeffs();
bigintmat *m = new bigintmat(A->rows()+A->cols(), A->cols(), R);
m->copySubmatInto(A, 1, 1, A->rows(), A->cols(), A->cols()+1, 1);
number one = n_Init(1, R);
for(int i=1; i<= A->cols(); i++)
m->set(i,i,one);
n_Delete(&one, R);
return m;
}
static number bimFarey(bigintmat *A, number N, bigintmat *L)
{
coeffs Z = A->basecoeffs(),
Q = nInitChar(n_Q, 0);
number den = n_Init(1, Z);
nMapFunc f = n_SetMap(Q, Z);
for(int i=1; i<= A->rows(); i++)
{
for(int j=1; j<= A->cols(); j++)
{
number ad = n_Mult(den, A->view(i, j), Z);
number re = n_IntMod(ad, N, Z);
n_Delete(&ad, Z);
number q = n_Farey(re, N, Z);
n_Delete(&re, Z);
if (!q)
{
n_Delete(&ad, Z);
n_Delete(&den, Z);
return NULL;
}
number d = n_GetDenom(q, Q),
n = n_GetNumerator(q, Q);
n_Delete(&q, Q);
n_Delete(&ad, Z);
number dz = f(d, Q, Z),
nz = f(n, Q, Z);
n_Delete(&d, Q);
n_Delete(&n, Q);
if (!n_IsOne(dz, Z))
{
L->skalmult(dz, Z);
n_InpMult(den, dz, Z);
#if 0
PrintS("den increasing to ");
n_Print(den, Z);
PrintLn();
#endif
}
n_Delete(&dz, Z);
L->rawset(i, j, nz);
}
}
nKillChar(Q);
PrintS("bimFarey worked\n");
#if 0
L->Print();
PrintS("\n * 1/");
n_Print(den, Z);
PrintLn();
#endif
return den;
}
static number solveAx_dixon(bigintmat *A, bigintmat *B, bigintmat *x, bigintmat *kern) {
coeffs R = A->basecoeffs();
assume(getCoeffType(R) == n_Z);
number p = n_Init(536870909, R); // PreviousPrime(2^29); not clever
coeffs Rp = numbercoeffs(p, R); // R/pR
bigintmat *Ap = bimChangeCoeff(A, Rp),
*m = prependIdentity(Ap),
*Tp, *Hp;
delete Ap;
m->howell();
Hp = new bigintmat(A->rows(), A->cols(), Rp);
Hp->copySubmatInto(m, A->cols()+1, 1, A->rows(), A->cols(), 1, 1);
Tp = new bigintmat(A->cols(), A->cols(), Rp);
Tp->copySubmatInto(m, 1, 1, A->cols(), A->cols(), 1, 1);
int i, j;
for(i=1; i<= A->cols(); i++)
{
for(j=m->rows(); j>A->cols(); j--)
{
if (!n_IsZero(m->view(j, i), Rp)) break;
}
if (j>A->cols()) break;
}
// Print("Found nullity (kern dim) of %d\n", i-1);
bigintmat * kp = new bigintmat(A->cols(), i-1, Rp);
kp->copySubmatInto(Tp, 1, 1, A->cols(), i-1, 1, 1);
kp->howell();
delete m;
//Hp is the mod-p howell form
//Tp the transformation, mod p
//kp a basis for the kernel, in howell form, mod p
bigintmat * eps_p = new bigintmat(B->rows(), B->cols(), Rp),
* x_p = new bigintmat(A->cols(), B->cols(), Rp),
* fps_p = new bigintmat(kp->cols(), B->cols(), Rp);
//initial solution
number zero = n_Init(0, R);
x->skalmult(zero, R);
n_Delete(&zero, R);
bigintmat * b = new bigintmat(B);
number pp = n_Init(1, R);
i = 1;
do
{
bigintmat * b_p = bimChangeCoeff(b, Rp), * s;
bigintmat * t1, *t2;
reduce_mod_howell(Hp, b_p, eps_p, x_p);
delete b_p;
if (!eps_p->isZero())
{
PrintS("no solution, since no modular solution\n");
delete eps_p;
delete x_p;
delete Hp;
delete kp;
delete Tp;
delete b;
n_Delete(&pp, R);
n_Delete(&p, R);
nKillChar(Rp);
return NULL;
}
t1 = bimMult(Tp, x_p);
delete x_p;
x_p = t1;
reduce_mod_howell(kp, x_p, x_p, fps_p); //we're not all interested in fps_p
s = bimChangeCoeff(x_p, R);
t1 = bimMult(A, s);
t2 = bimSub(b, t1);
t2->skaldiv(p);
delete b;
delete t1;
b = t2;
s->skalmult(pp, R);
t1 = bimAdd(x, s);
delete s;
x->swapMatrix(t1);
delete t1;
if(kern && i==1)
{
bigintmat * ker = bimChangeCoeff(kp, R);
t1 = bimMult(A, ker);
t1->skaldiv(p);
t1->skalmult(n_Init(-1, R), R);
b->appendCol(t1);
delete t1;
x->appendCol(ker);
delete ker;
x_p->extendCols(kp->cols());
eps_p->extendCols(kp->cols());
fps_p->extendCols(kp->cols());
}
n_InpMult(pp, p, R);
if (b->isZero())
{
//exact solution found, stop
delete eps_p;
delete fps_p;
delete x_p;
delete Hp;
delete kp;
delete Tp;
delete b;
n_Delete(&pp, R);
n_Delete(&p, R);
nKillChar(Rp);
return n_Init(1, R);
}
else
{
bigintmat *y = new bigintmat(x->rows(), x->cols(), R);
number d = bimFarey(x, pp, y);
if (d)
{
bigintmat *c = bimMult(A, y);
bigintmat *bd = new bigintmat(B);
bd->skalmult(d, R);
if (kern)
{
bd->extendCols(kp->cols());
}
if (*c == *bd)
{
x->swapMatrix(y);
delete y;
delete c;
if (kern)
{
y = new bigintmat(x->rows(), B->cols(), R);
c = new bigintmat(x->rows(), kp->cols(), R);
x->splitcol(y, c);
x->swapMatrix(y);
delete y;
kern->swapMatrix(c);
delete c;
}
delete bd;
delete eps_p;
delete fps_p;
delete x_p;
delete Hp;
delete kp;
delete Tp;
delete b;
n_Delete(&pp, R);
n_Delete(&p, R);
nKillChar(Rp);
return d;
}
delete c;
delete bd;
n_Delete(&d, R);
}
delete y;
}
i++;
} while (1);
delete eps_p;
delete fps_p;
delete x_p;
delete Hp;
delete kp;
delete Tp;
n_Delete(&pp, R);
n_Delete(&p, R);
nKillChar(Rp);
return NULL;
}
//TODO: re-write using reduce_mod_howell
static number solveAx_howell(bigintmat *A, bigintmat *b, bigintmat *x, bigintmat *kern)
{
// try to solve Ax=b, more precisely, find
// number d
// bigintmat x
// sth. Ax=db
// where d is small-ish (divides the determinant of A if this makes sense)
// return 0 if there is no solution.
//
// if kern is non-NULL, return a basis for the kernel
//Algo: we do row-howell (triangular matrix). The idea is
// Ax = b <=> AT T^-1x = b
// y := T^-1 x, solve AT y = b
// and return Ty.
//Howell does not compute the trafo, hence we need to cheat:
//B := (I_n | A^t)^t, then the top part of the Howell form of
//B will give a useful trafo
//Then we can find x by back-substitution and lcm/gcd to find the denominator
//The defining property of Howell makes this work.
coeffs R = A->basecoeffs();
bigintmat *m = prependIdentity(A);
m->howell(); // since m contains the identity, we'll have A->cols()
// many cols.
number den = n_Init(1, R);
bigintmat * B = new bigintmat(A->rows(), 1, R);
for(int i=1; i<= b->cols(); i++)
{
int A_col = A->cols();
b->getcol(i, B);
B->skalmult(den, R);
for(int j = B->rows(); j>0; j--)
{
number Ai = m->view(m->rows()-B->rows() + j, A_col);
if (n_IsZero(Ai, R) &&
n_IsZero(B->view(j, 1), R))
{
continue; //all is fine: 0*x = 0
}
else if (n_IsZero(B->view(j, 1), R))
{
x->rawset(x->rows() - B->rows() + j, i, n_Init(0, R));
A_col--;
}
else if (n_IsZero(Ai, R))
{
delete m;
delete B;
n_Delete(&den, R);
return 0;
}
else
{
// solve ax=db, possibly enlarging d
// so x = db/a
number Bj = B->view(j, 1);
number g = n_Gcd(Bj, Ai, R);
number xi;
if (n_Equal(Ai, g, R))
{ //good: den stable!
xi = n_Div(Bj, Ai, R);
}
else
{ //den <- den * (a/g), so old sol. needs to be adjusted
number inc_d = n_Div(Ai, g, R);
n_InpMult(den, inc_d, R);
x->skalmult(inc_d, R);
B->skalmult(inc_d, R);
xi = n_Div(Bj, g, R);
n_Delete(&inc_d, R);
} //now for the back-substitution:
x->rawset(x->rows() - B->rows() + j, i, xi);
for(int k=j; k>0; k--)
{
//B[k] = B[k] - x[k]A[k][j]
number s = n_Mult(xi, m->view(m->rows()-B->rows() + k, A_col), R);
B->rawset(k, 1, n_Sub(B->view(k, 1), s, R));
n_Delete(&s, R);
}
n_Delete(&g, R);
A_col--;
}
if (!A_col)
{
if (B->isZero()) break;
else
{
delete m;
delete B;
n_Delete(&den, R);
return 0;
}
}
}
}
delete B;
bigintmat *T = new bigintmat(A->cols(), A->cols(), R);
T->copySubmatInto(m, 1, 1, A->cols(), A->cols(), 1, 1);
if (kern)
{
int i, j;
for(i=1; i<= A->cols(); i++)
{
for(j=m->rows(); j>A->cols(); j--)
{
if (!n_IsZero(m->view(j, i), R)) break;
}
if (j>A->cols()) break;
}
Print("Found nullity (kern dim) of %d\n", i-1);
bigintmat * ker = new bigintmat(A->rows(), i-1, R);
ker->copySubmatInto(T, 1, 1, A->rows(), i-1, 1, 1);
kern->swapMatrix(ker);
delete ker;
}
delete m;
bigintmat * y = bimMult(T, x);
x->swapMatrix(y);
delete y;
x->simplifyContentDen(&den);
#if 0
PrintS("sol = 1/");
n_Print(den, R);
PrintS(" *\n");
x->Print();
PrintLn();
#endif
return den;
}
number solveAx(bigintmat *A, bigintmat *b, bigintmat *x)
{
#if 0
PrintS("Solve Ax=b for A=\n");
A->Print();
PrintS("\nb = \n");
b->Print();
PrintS("\nx = \n");
x->Print();
PrintLn();
#endif
coeffs R = A->basecoeffs();
assume (R == b->basecoeffs());
assume (R == x->basecoeffs());
assume ((x->cols() == b->cols()) && (x->rows() == A->cols()) && (A->rows() == b->rows()));
switch (getCoeffType(R))
{
case n_Z:
return solveAx_dixon(A, b, x, NULL);
case n_Zn:
case n_Znm:
case n_Z2m:
return solveAx_howell(A, b, x, NULL);
case n_Zp:
case n_Q:
case n_GF:
case n_algExt:
case n_transExt:
WarnS("have field, should use Gauss or better");
break;
default:
if (R->cfXExtGcd && R->cfAnn)
{ //assume it's Euclidean
return solveAx_howell(A, b, x, NULL);
}
WerrorS("have no solve algorithm");
break;
}
return NULL;
}
void diagonalForm(bigintmat *A, bigintmat ** S, bigintmat ** T)
{
bigintmat * t, *s, *a=A;
coeffs R = a->basecoeffs();
if (T)
{
*T = new bigintmat(a->cols(), a->cols(), R),
(*T)->one();
t = new bigintmat(*T);
}
else
{
t = *T;
}
if (S)
{
*S = new bigintmat(a->rows(), a->rows(), R);
(*S)->one();
s = new bigintmat(*S);
}
else
{
s = *S;
}
int flip=0;
do
{
bigintmat * x, *X;
if (flip)
{
x = s;
X = *S;
}
else
{
x = t;
X = *T;
}
if (x)
{
x->one();
bigintmat * r = new bigintmat(a->rows()+a->cols(), a->cols(), R);
bigintmat * rw = new bigintmat(1, a->cols(), R);
for(int i=0; i<a->cols(); i++)
{
x->getrow(i+1, rw);
r->setrow(i+1, rw);
}
for (int i=0; i<a->rows(); i++)
{
a->getrow(i+1, rw);
r->setrow(i+a->cols()+1, rw);
}
r->hnf();
for(int i=0; i<a->cols(); i++)
{
r->getrow(i+1, rw);
x->setrow(i+1, rw);
}
for(int i=0; i<a->rows(); i++)
{
r->getrow(i+a->cols()+1, rw);
a->setrow(i+1, rw);
}
delete rw;
delete r;
#if 0
Print("X: %ld\n", X);
X->Print();
Print("\nx: %ld\n", x);
x->Print();
#endif
bimMult(X, x, X);
#if 0
Print("\n2:X: %ld %ld %ld\n", X, *S, *T);
X->Print();
Print("\n2:x: %ld\n", x);
x->Print();
PrintLn();
#endif
}
else
{
a->hnf();
}
int diag = 1;
for(int i=a->rows(); diag && i>0; i--)
{
for(int j=a->cols(); j>0; j--)
{
if ((a->rows()-i)!=(a->cols()-j) && !n_IsZero(a->view(i, j), R))
{
diag = 0;
break;
}
}
}
#if 0
PrintS("Diag ? %d\n", diag);
a->Print();
PrintLn();
#endif
if (diag) break;
a = a->transpose(); // leaks - I need to write inpTranspose
flip = 1-flip;
} while (1);
if (flip)
a = a->transpose();
if (S) *S = (*S)->transpose();
if (s) delete s;
if (t) delete t;
A->copy(a);
}
//a "q-base" for the kernel of a.
//Should be re-done to use Howell rather than smith.
//can be done using solveAx now.
int kernbase (bigintmat *a, bigintmat *c, number p, coeffs q)
{
#if 0
PrintS("Kernel of ");
a->Print();
PrintS(" modulo ");
n_Print(p, q);
PrintLn();
#endif
coeffs coe = numbercoeffs(p, q);
bigintmat *m = bimChangeCoeff(a, coe), *U, *V;
diagonalForm(m, &U, &V);
#if 0
PrintS("\ndiag form: ");
m->Print();
PrintS("\nU:\n");
U->Print();
PrintS("\nV:\n");
V->Print();
PrintLn();
#endif
int rg = 0;
#undef MIN
#define MIN(a,b) (a < b ? a : b)
for(rg=0; rg<MIN(m->rows(), m->cols()) && !n_IsZero(m->view(m->rows()-rg,m->cols()-rg), coe); rg++);
bigintmat * k = new bigintmat(m->cols(), m->rows(), coe);
for(int i=0; i<rg; i++)
{
number A = n_Ann(m->view(m->rows()-i, m->cols()-i), coe);
k->set(m->cols()-i, i+1, A);
n_Delete(&A, coe);
}
for(int i=rg; i<m->cols(); i++)
{
k->set(m->cols()-i, i+1-rg, n_Init(1, coe));
}
bimMult(V, k, k);
c->copy(bimChangeCoeff(k, q));
return c->cols();
}
bool nCoeffs_are_equal(coeffs r, coeffs s)
{
if ((r == NULL) || (s == NULL))
return false;
if (r == s)
return true;
if ((getCoeffType(r)==n_Z) && (getCoeffType(s)==n_Z))
return true;
if ((getCoeffType(r)==n_Zp) && (getCoeffType(s)==n_Zp))
{
if (r->ch == s->ch)
return true;
else
return false;
}
// n_Zn stimmt wahrscheinlich noch nicht
if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
{
if (r->ch == s->ch)
return true;
else
return false;
}
if ((getCoeffType(r)==n_Q) && (getCoeffType(s)==n_Q))
return true;
// FALL n_Zn FEHLT NOCH!
//if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
return false;
}
number bigintmat::content()
{
coeffs r = basecoeffs();
number g = get(1,1), h;
int n=rows()*cols();
for(int i=1; i<n && !n_IsOne(g, r); i++)
{
h = n_Gcd(g, view(i), r);
n_Delete(&g, r);
g=h;
}
return g;
}
void bigintmat::simplifyContentDen(number *d)
{
coeffs r = basecoeffs();
number g = n_Copy(*d, r), h;
int n=rows()*cols();
for(int i=0; i<n && !n_IsOne(g, r); i++)
{
h = n_Gcd(g, view(i), r);
n_Delete(&g, r);
g=h;
}
*d = n_Div(*d, g, r);
if (!n_IsOne(g, r))
skaldiv(g);
}
auxiliary.h
All the auxiliary stuff.
TRUE
#define TRUE
Definition auxiliary.h:101
FALSE
#define FALSE
Definition auxiliary.h:97
si_min
static int si_min(const int a, const int b)
Definition auxiliary.h:126
reduce_mod_howell
static void reduce_mod_howell(bigintmat *A, bigintmat *b, bigintmat *eps, bigintmat *x)
Definition bigintmat.cc:1940
kernbase
int kernbase(bigintmat *a, bigintmat *c, number p, coeffs q)
a basis for the nullspace of a mod p: only used internally in Round2. Don't use it.
Definition bigintmat.cc:2585
swap
#define swap(_i, _j)
solveAx
number solveAx(bigintmat *A, bigintmat *b, bigintmat *x)
solve Ax=b*d. x needs to be pre-allocated to the same number of columns as b. the minimal denominator...
Definition bigintmat.cc:2418
solveAx_dixon
static number solveAx_dixon(bigintmat *A, bigintmat *B, bigintmat *x, bigintmat *kern)
Definition bigintmat.cc:2097
bimChangeCoeff
bigintmat * bimChangeCoeff(bigintmat *a, coeffs cnew)
Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen.
Definition bigintmat.cc:1797
MIN
#define MIN(a, b)
bimCopy
bigintmat * bimCopy(const bigintmat *b)
same as copy constructor - apart from it being able to accept NULL as input
Definition bigintmat.cc:403
prependIdentity
static bigintmat * prependIdentity(bigintmat *A)
Definition bigintmat.cc:2026
nCoeffs_are_equal
bool nCoeffs_are_equal(coeffs r, coeffs s)
Definition bigintmat.cc:2629
bim2iv
intvec * bim2iv(bigintmat *b)
Definition bigintmat.cc:339
bimMult
bigintmat * bimMult(bigintmat *a, bigintmat *b)
Definition bigintmat.cc:253
intArrSum
static int intArrSum(int *a, int length)
Definition bigintmat.cc:523
bimSub
bigintmat * bimSub(bigintmat *a, bigintmat *b)
Definition bigintmat.cc:216
solveAx_howell
static number solveAx_howell(bigintmat *A, bigintmat *b, bigintmat *x, bigintmat *kern)
Definition bigintmat.cc:2286
diagonalForm
void diagonalForm(bigintmat *A, bigintmat **S, bigintmat **T)
Definition bigintmat.cc:2461
iv2bim
bigintmat * iv2bim(intvec *b, const coeffs C)
Definition bigintmat.cc:347
findLongest
static int findLongest(int *a, int length)
Definition bigintmat.cc:531
getShorter
static int getShorter(int *a, int l, int j, int cols, int rows)
Definition bigintmat.cc:546
bimFarey
static number bimFarey(bigintmat *A, number N, bigintmat *L)
Definition bigintmat.cc:2038
operator!=
bool operator!=(const bigintmat &lhr, const bigintmat &rhr)
Definition bigintmat.cc:174
numbercoeffs
static coeffs numbercoeffs(number n, coeffs c)
create Z/nA of type n_Zn
Definition bigintmat.cc:20
bimAdd
bigintmat * bimAdd(bigintmat *a, bigintmat *b)
Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ? @Note: NULL as a result means an error (non-compati...
Definition bigintmat.cc:180
operator==
bool operator==(const bigintmat &lhr, const bigintmat &rhr)
Definition bigintmat.cc:157
bigintmat.h
BIMATELEM
#define BIMATELEM(M, I, J)
Definition bigintmat.h:133
bimChangeCoeff
bigintmat * bimChangeCoeff(bigintmat *a, coeffs cnew)
Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen.
Definition bigintmat.cc:1797
nCoeffs_are_equal
bool nCoeffs_are_equal(coeffs r, coeffs s)
Definition bigintmat.cc:2629
pp
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f ).
Definition cf_gcd.cc:676
den
CanonicalForm den(const CanonicalForm &f)
Definition canonicalform.h:340
N
const CanonicalForm CFMap CFMap & N
Definition cfEzgcd.cc:56
l
int l
Definition cfEzgcd.cc:100
m
int m
Definition cfEzgcd.cc:128
i
int i
Definition cfEzgcd.cc:132
k
int k
Definition cfEzgcd.cc:99
x
Variable x
Definition cfModGcd.cc:4090
p
int p
Definition cfModGcd.cc:4086
g
g
Definition cfModGcd.cc:4098
cf
CanonicalForm cf
Definition cfModGcd.cc:4091
b
CanonicalForm b
Definition cfModGcd.cc:4111
f
FILE * f
Definition checklibs.c:9
bigintmat
Matrices of numbers.
Definition bigintmat.h:51
bigintmat::Print
void Print()
IO: simply prints the matrix to the current output (screen?).
Definition bigintmat.cc:437
bigintmat::colskaldiv
void colskaldiv(int j, number b)
Macht Ganzzahldivision aller j-ten Spalteneinträge mit b.
Definition bigintmat.cc:1866
bigintmat::det
number det()
det (via LaPlace in general, hnf for euc. rings)
Definition bigintmat.cc:1507
bigintmat::splitrow
void splitrow(bigintmat *a, bigintmat *b)
Speichert in Matrix a den oberen, in b den unteren Teil der Matrix, vorausgesetzt die Dimensionen sti...
Definition bigintmat.cc:1122
bigintmat::isOne
int isOne()
is matrix is identity
Definition bigintmat.cc:1295
bigintmat::zero
void zero()
Setzt alle Einträge auf 0.
Definition bigintmat.cc:1345
bigintmat::appendCol
void appendCol(bigintmat *a)
horizontally join the matrices, m <- m|a
Definition bigintmat.cc:1078
bigintmat::trace
number trace()
the trace ....
Definition bigintmat.cc:1493
bigintmat::addrow
bool addrow(int i, int j, number a, coeffs c)
... Zeile ...
Definition bigintmat.cc:978
bigintmat::swap
void swap(int i, int j)
swap columns i and j
Definition bigintmat.cc:680
bigintmat::coltransform
void coltransform(int i, int j, number a, number b, number c, number d)
transforms cols (i,j) using the 2x2 matrix ((a,b)(c,d)) (hopefully)
Definition bigintmat.cc:1879
bigintmat::addcol
bool addcol(int i, int j, number a, coeffs c)
addiert a-faches der j-ten Spalte zur i-ten dazu
Definition bigintmat.cc:954
bigintmat::v
number * v
Definition bigintmat.h:54
bigintmat::hnf
void hnf()
transforms INPLACE to HNF
Definition bigintmat.cc:1655
bigintmat::StringAsPrinted
char * StringAsPrinted()
Returns a string as it would have been printed in the interpreter.
Definition bigintmat.cc:445
bigintmat::swapMatrix
void swapMatrix(bigintmat *a)
Definition bigintmat.cc:1561
bigintmat::cols
int cols() const
Definition bigintmat.h:144
bigintmat::isZero
int isZero()
Definition bigintmat.cc:1358
bigintmat::hnfdet
number hnfdet()
det via HNF Primzahlen als long long int, müssen noch in number umgewandelt werden?
Definition bigintmat.cc:1540
bigintmat::findnonzero
int findnonzero(int i)
find index of 1st non-zero entry in row i
Definition bigintmat.cc:718
bigintmat::modhnf
bigintmat * modhnf(number p, coeffs c)
computes HNF(this | p*I)
Definition bigintmat.cc:1824
bigintmat::setcol
void setcol(int j, bigintmat *m)
Setzt j-te Spalte gleich übergebenem Vektor (Matrix) m.
Definition bigintmat.cc:821
bigintmat::add
bool add(bigintmat *b)
Addiert zur Matrix die Matrix b dazu. Return false => an error occurred.
Definition bigintmat.cc:889
bigintmat::getwid
int * getwid(int maxwid)
Definition bigintmat.cc:574
bigintmat::inpMult
void inpMult(number bintop, const coeffs C=NULL)
inplace version of skalar mult. CHANGES input.
Definition bigintmat.cc:143
bigintmat::transpose
bigintmat * transpose()
Definition bigintmat.cc:35
bigintmat::setrow
void setrow(int i, bigintmat *m)
Setzt i-te Zeile gleich übergebenem Vektor (Matrix) m.
Definition bigintmat.cc:855
bigintmat::Write
void Write()
IO: writes the matrix into the current internal string buffer which must be started/ allocated before...
Definition bigintmat.cc:411
bigintmat::rowskalmult
void rowskalmult(int i, number a, coeffs c)
... Zeile ...
Definition bigintmat.cc:1018
bigintmat::skalmult
bool skalmult(number b, coeffs c)
Multipliziert zur Matrix den Skalar b hinzu.
Definition bigintmat.cc:933
bigintmat::simplifyContentDen
void simplifyContentDen(number *den)
ensures that Gcd(den, content)=1 enden hier wieder
Definition bigintmat.cc:2672
bigintmat::extendCols
void extendCols(int i)
append i zero-columns to the matrix
Definition bigintmat.cc:1071
bigintmat::splitcol
void splitcol(bigintmat *a, bigintmat *b)
... linken ... rechten ...
Definition bigintmat.cc:1164
bigintmat::content
number content()
the content, the gcd of all entries. Only makes sense for Euclidean rings (or possibly constructive P...
Definition bigintmat.cc:2659
bigintmat::skaldiv
void skaldiv(number b)
Macht Ganzzahldivision aller Matrixeinträge mit b.
Definition bigintmat.cc:1851
bigintmat::colskalmult
void colskalmult(int i, number a, coeffs c)
Multipliziert zur i-ten Spalte den Skalar a hinzu.
Definition bigintmat.cc:1002
bigintmat::colIsZero
int colIsZero(int i)
Definition bigintmat.cc:1572
bigintmat::copy
bool copy(bigintmat *b)
Kopiert Einträge von b auf Bigintmat.
Definition bigintmat.cc:1254
bigintmat::index
int index(int r, int c) const
helper function to map from 2-dim coordinates, starting by 1 to 1-dim coordinate, starting by 0
Definition bigintmat.h:161
bigintmat::getcol
void getcol(int j, bigintmat *a)
copies the j-th column into the matrix a - which needs to be pre-allocated with the correct size.
Definition bigintmat.cc:742
bigintmat::pseudoinv
number pseudoinv(bigintmat *a)
Speichert in Matrix a die Pseudoinverse, liefert den Nenner zurück.
Definition bigintmat.cc:1410
bigintmat::col
int col
Definition bigintmat.h:56
bigintmat::rows
int rows() const
Definition bigintmat.h:145
bigintmat::get
number get(int i, int j) const
get a copy of an entry. NOTE: starts at [1,1]
Definition bigintmat.cc:117
bigintmat::pprint
void pprint(int maxwid)
Definition bigintmat.cc:604
bigintmat::getColRange
void getColRange(int j, int no, bigintmat *a)
copies the no-columns staring by j (so j...j+no-1) into the pre-allocated a
Definition bigintmat.cc:773
bigintmat::concatcol
void concatcol(bigintmat *a, bigintmat *b)
Definition bigintmat.cc:1093
bigintmat::findcolnonzero
int findcolnonzero(int j)
find index of 1st non-zero entry in column j
Definition bigintmat.cc:730
bigintmat::copySubmatInto
void copySubmatInto(bigintmat *, int sr, int sc, int nr, int nc, int tr, int tc)
copy the submatrix of b, staring at (a,b) having n rows, m cols into the given matrix at pos....
Definition bigintmat.cc:1282
bigintmat::view
number view(int i, int j) const
view an entry an entry. NOTE: starts at [1,1]
Definition bigintmat.cc:125
bigintmat::inpTranspose
void inpTranspose()
transpose in place
Definition bigintmat.cc:48
bigintmat::one
void one()
Macht Matrix (Falls quadratisch) zu Einheitsmatrix.
Definition bigintmat.cc:1320
bigintmat::howell
void howell()
dito, but Howell form (only different for zero-divsors)
Definition bigintmat.cc:1580
bigintmat::rawset
void rawset(int i, number n, const coeffs C=NULL)
replace an entry with the given number n (only delete old). NOTE: starts at [0]. Should be named set_...
Definition bigintmat.h:196
bigintmat::operator*=
void operator*=(int intop)
UEberladener *=-Operator (fuer int und bigint) Frage hier: *= verwenden oder lieber = und * einzeln?...
Definition bigintmat.cc:134
bigintmat::basecoeffs
coeffs basecoeffs() const
Definition bigintmat.h:146
bigintmat::getrow
void getrow(int i, bigintmat *a)
Schreibt i-te Zeile in Vektor (Matrix) a.
Definition bigintmat.cc:786
bigintmat::concatrow
void concatrow(bigintmat *a, bigintmat *b)
Fügt zwei Matrixen untereinander/nebeneinander in gegebene Matrix ein, bzw spaltet gegebenen Matrix a...
Definition bigintmat.cc:1034
bigintmat::bigintmat
bigintmat()
Definition bigintmat.h:59
bigintmat::set
void set(int i, int j, number n, const coeffs C=NULL)
replace an entry with a copy (delete old + copy new!). NOTE: starts at [1,1]
Definition bigintmat.cc:93
bigintmat::elim
bigintmat * elim(int i, int j)
Liefert Streichungsmatrix (i-te Zeile und j-te Spalte gestrichen) zurück.
Definition bigintmat.cc:1376
bigintmat::compare
int compare(const bigintmat *op) const
Definition bigintmat.cc:360
bigintmat::sub
bool sub(bigintmat *b)
Subtrahiert ...
Definition bigintmat.cc:911
bigintmat::row
int row
Definition bigintmat.h:55
bigintmat::String
char * String()
IO: String returns a singular string containing the matrix, needs freeing afterwards.
Definition bigintmat.cc:430
bigintmat::swaprow
void swaprow(int i, int j)
swap rows i and j
Definition bigintmat.cc:699
bigintmat::mod
void mod(number p)
Reduziert komplette Matrix modulo p.
Definition bigintmat.cc:1906
intvec
Definition intvec.h:24
n_Mult
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of 'a' and 'b', i.e., a*b
Definition coeffs.h:639
number2mpz
static FORCE_INLINE void number2mpz(number n, coeffs c, mpz_t m)
Definition coeffs.h:990
n_Int
static FORCE_INLINE long n_Int(number &n, const coeffs r)
conversion of n to an int; 0 if not possible in Z/pZ: the representing int lying in (-p/2 ....
Definition coeffs.h:550
n_Copy
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition coeffs.h:457
n_Add
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of 'a' and 'b', i.e., a+b
Definition coeffs.h:653
n_GetDenom
static FORCE_INLINE number n_GetDenom(number &n, const coeffs r)
return the denominator of n (if elements of r are by nature not fractional, result is 1)
Definition coeffs.h:606
n_GF
@ n_GF
\GF{p^n < 2^16}
Definition coeffs.h:32
n_Q
@ n_Q
rational (GMP) numbers
Definition coeffs.h:30
n_Znm
@ n_Znm
only used if HAVE_RINGS is defined
Definition coeffs.h:45
n_algExt
@ n_algExt
used for all algebraic extensions, i.e., the top-most extension in an extension tower is algebraic
Definition coeffs.h:35
n_Zn
@ n_Zn
only used if HAVE_RINGS is defined
Definition coeffs.h:44
n_Z2m
@ n_Z2m
only used if HAVE_RINGS is defined
Definition coeffs.h:46
n_Zp
@ n_Zp
\F{p < 2^31}
Definition coeffs.h:29
n_transExt
@ n_transExt
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition coeffs.h:38
n_Z
@ n_Z
only used if HAVE_RINGS is defined
Definition coeffs.h:43
n_Gcd
static FORCE_INLINE number n_Gcd(number a, number b, const coeffs r)
in Z: return the gcd of 'a' and 'b' in Z/nZ, Z/2^kZ: computed as in the case Z in Z/pZ,...
Definition coeffs.h:667
n_QuotRem
static FORCE_INLINE number n_QuotRem(number a, number b, number *q, const coeffs r)
Definition coeffs.h:684
nCoeffString
static FORCE_INLINE char * nCoeffString(const coeffs cf)
TODO: make it a virtual method of coeffs, together with: Decompose & Compose, rParameter & rPar.
Definition coeffs.h:958
n_GreaterZero
static FORCE_INLINE BOOLEAN n_GreaterZero(number n, const coeffs r)
ordered fields: TRUE iff 'n' is positive; in Z/pZ: TRUE iff 0 < m <= roundedBelow(p/2),...
Definition coeffs.h:500
n_Ann
static FORCE_INLINE number n_Ann(number a, const coeffs r)
if r is a ring with zero divisors, return an annihilator!=0 of b otherwise return NULL
Definition coeffs.h:682
n_Print
void n_Print(number &a, const coeffs r)
print a number (BEWARE of string buffers!) mostly for debugging
Definition numbers.cc:662
n_SetMap
static FORCE_INLINE nMapFunc n_SetMap(const coeffs src, const coeffs dst)
set the mapping function pointers for translating numbers from src to dst
Definition coeffs.h:703
n_InpNeg
static FORCE_INLINE number n_InpNeg(number n, const coeffs r)
in-place negation of n MUST BE USED: n = n_InpNeg(n) (no copy is returned)
Definition coeffs.h:560
n_Farey
static FORCE_INLINE number n_Farey(number a, number b, const coeffs r)
Definition coeffs.h:762
n_Div
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition coeffs.h:618
nInitChar
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition numbers.cc:412
n_Greater
static FORCE_INLINE BOOLEAN n_Greater(number a, number b, const coeffs r)
ordered fields: TRUE iff 'a' is larger than 'b'; in Z/pZ: TRUE iff la > lb, where la and lb are the l...
Definition coeffs.h:517
n_IsZero
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
Definition coeffs.h:470
n_GetUnit
static FORCE_INLINE number n_GetUnit(number n, const coeffs r)
in Z: 1 in Z/kZ (where k is not a prime): largest divisor of n (taken in Z) that is co-prime with k i...
Definition coeffs.h:537
n_Sub
static FORCE_INLINE number n_Sub(number a, number b, const coeffs r)
return the difference of 'a' and 'b', i.e., a-b
Definition coeffs.h:658
getCoeffType
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition coeffs.h:431
n_Delete
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition coeffs.h:461
n_Write
static FORCE_INLINE void n_Write(number n, const coeffs r, const BOOLEAN bShortOut=TRUE)
Definition coeffs.h:594
n_Init
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition coeffs.h:541
n_IntMod
static FORCE_INLINE number n_IntMod(number a, number b, const coeffs r)
for r a field, return n_Init(0,r) always: n_Div(a,b,r)*b+n_IntMod(a,b,r)==a n_IntMod(a,...
Definition coeffs.h:631
n_InpMult
static FORCE_INLINE void n_InpMult(number &a, number b, const coeffs r)
multiplication of 'a' and 'b'; replacement of 'a' by the product a*b
Definition coeffs.h:644
n_Equal
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff 'a' and 'b' represent the same number; they may have different representations.
Definition coeffs.h:466
n_GetNumerator
static FORCE_INLINE number n_GetNumerator(number &n, const coeffs r)
return the numerator of n (if elements of r are by nature not fractional, result is n)
Definition coeffs.h:611
nMapFunc
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition coeffs.h:80
nKillChar
void nKillChar(coeffs r)
undo all initialisations
Definition numbers.cc:563
n_InpAdd
static FORCE_INLINE void n_InpAdd(number &a, number b, const coeffs r)
addition of 'a' and 'b'; replacement of 'a' by the sum a+b
Definition coeffs.h:649
n_IsOne
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition coeffs.h:474
n_XExtGcd
static FORCE_INLINE number n_XExtGcd(number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
Definition coeffs.h:676
Print
#define Print
Definition emacs.cc:80
WarnS
#define WarnS
Definition emacs.cc:78
StringAppend
#define StringAppend
Definition emacs.cc:79
s
const CanonicalForm int s
Definition facAbsFact.cc:51
y
const CanonicalForm int const CFList const Variable & y
Definition facAbsFact.cc:53
B
b *CanonicalForm B
Definition facBivar.cc:52
tmp1
CFList tmp1
Definition facFqBivar.cc:75
tmp2
CFList tmp2
Definition facFqBivar.cc:75
j
int j
Definition facHensel.cc:110
prod
fq_nmod_poly_t prod
Definition facHensel.cc:100
min
static int min(int a, int b)
Definition fast_mult.cc:268
WerrorS
void WerrorS(const char *s)
Definition feFopen.cc:24
flip
std::pair< ideal, ring > flip(const ideal I, const ring r, const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal, const gfan::ZVector adjustedInteriorPoint, const gfan::ZVector adjustedFacetNormal)
Definition flip.cc:17
length
static BOOLEAN length(leftv result, leftv arg)
Definition interval.cc:257
intvec.h
T
STATIC_VAR jList * T
Definition janet.cc:30
h
STATIC_VAR Poly * h
Definition janet.cc:971
lcm
int lcm(unsigned long *l, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
Definition minpoly.cc:709
assume
#define assume(x)
Definition mod2.h:389
coeffs
The main handler for Singular numbers which are suitable for Singular polynomials.
omAlloc
#define omAlloc(size)
Definition omAllocDecl.h:210
omFree
#define omFree(addr)
Definition omAllocDecl.h:261
omAlloc0
#define omAlloc0(size)
Definition omAllocDecl.h:211
NULL
#define NULL
Definition omList.c:12
index
static int index(p_Length length, p_Ord ord)
Definition p_Procs_Impl.h:592
StringSetS
void StringSetS(const char *st)
Definition reporter.cc:128
StringAppendS
void StringAppendS(const char *st)
Definition reporter.cc:107
PrintS
void PrintS(const char *s)
Definition reporter.cc:288
StringEndS
char * StringEndS()
Definition reporter.cc:151
PrintLn
void PrintLn()
Definition reporter.cc:314
Werror
void Werror(const char *fmt,...)
Definition reporter.cc:189
rmodulon.h
ZnmInfo
Definition rmodulon.h:17
R
#define R
Definition sirandom.c:27
A
#define A
Definition sirandom.c:24
Q
#define Q
Definition sirandom.c:26
gcd
int gcd(int a, int b)
Definition walkSupport.cc:836