bigintmat.cc File Reference

#include <misc/auxiliary.h>
#include "bigintmat.h"
#include <misc/intvec.h>
#include "rmodulon.h"
#include <math.h>
#include <string.h>

Go to the source code of this file.

Macros
#define	swap(_i, _j)
#define	MIN(a, b)   (a < b ? a : b)
#define	MAX(a, b)   (a > b ? a : b)

Functions
static coeffs	numbercoeffs (number n, coeffs c)
	create Z/nA of type n_Zn More...
bool	operator== (const bigintmat &lhr, const bigintmat &rhr)
bool	operator!= (const bigintmat &lhr, const bigintmat &rhr)
bigintmat *	bimAdd (bigintmat *a, bigintmat *b)
	Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ? : NULL as a result means an error (non-compatible matrices?) More...
bigintmat *	bimAdd (bigintmat *a, int b)
bigintmat *	bimSub (bigintmat *a, bigintmat *b)
bigintmat *	bimSub (bigintmat *a, int b)
bigintmat *	bimMult (bigintmat *a, bigintmat *b)
bigintmat *	bimMult (bigintmat *a, int b)
bigintmat *	bimMult (bigintmat *a, number b, const coeffs cf)
intvec *	bim2iv (bigintmat *b)
bigintmat *	iv2bim (intvec *b, const coeffs C)
bigintmat *	bimCopy (const bigintmat *b)
	same as copy constructor - apart from it being able to accept NULL as input More...
static int	intArrSum (int *a, int length)
static int	findLongest (int *a, int length)
static int	getShorter (int *a, int l, int j, int cols, int rows)
bigintmat *	bimChangeCoeff (bigintmat *a, coeffs cnew)
	Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen. More...
void	bimMult (bigintmat *a, bigintmat *b, bigintmat *c)
	Multipliziert Matrix a und b und speichert Ergebnis in c. More...
static void	reduce_mod_howell (bigintmat *A, bigintmat *b, bigintmat *eps, bigintmat *x)
static bigintmat *	prependIdentity (bigintmat *A)
static number	bimFarey (bigintmat *A, number N, bigintmat *L)
static number	solveAx_dixon (bigintmat *A, bigintmat *B, bigintmat *x, bigintmat *kern)
static number	solveAx_howell (bigintmat *A, bigintmat *b, bigintmat *x, bigintmat *kern)
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 d is returned. Currently available for Z, Q and Z/nZ (and possibly for all fields: d=1 there) Beware that the internal functions can find the kernel as well - but the interface is lacking. More...
void	diagonalForm (bigintmat *A, bigintmat **S, bigintmat **T)
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. More...
bool	nCoeffs_are_equal (coeffs r, coeffs s)

Macro Definition Documentation

#define MAX	(		a,
			b
	)		(a > b ? a : b)

#define MIN	(		a,
			b
	)		(a < b ? a : b)

#define swap	(		_i,
			_j
	)

Value:

int __i = (_i), __j=(_j); \
  number c = v[__i];        \
  v[__i] = v[__j];          \
  v[__j] = c                \

Function Documentation

intvec* bim2iv

(

bigintmat *

b

)

Definition at line 339 of file bigintmat.cc.

{
  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* bimAdd	(	bigintmat *	a,
		bigintmat *	b
	)

Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ? : NULL as a result means an error (non-compatible matrices?)

Definition at line 178 of file bigintmat.cc.

{
  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,
		int	b
	)

Definition at line 195 of file bigintmat.cc.

{

  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_Add((*a)[i], bb, basecoeffs), basecoeffs);

  n_Delete(&bb,basecoeffs);
  return bim;
}

bigintmat* bimChangeCoeff	(	bigintmat *	a,
		coeffs	cnew
	)

Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen.

Definition at line 1688 of file bigintmat.cc.

{
  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;
}

bigintmat* bimCopy

(

const bigintmat *

b

)

same as copy constructor - apart from it being able to accept NULL as input

Definition at line 403 of file bigintmat.cc.

{
  if (b == NULL)
    return NULL;

  return new bigintmat(b);
}

static number bimFarey ( bigintmat *  A, number  N, bigintmat *  L  )

static

Definition at line 1919 of file bigintmat.cc.

                                                             {
  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
        Print("den increasing to ");
        n_Print(den, Z);
        Print("\n");
#endif
      }
      n_Delete(&dz, Z);
      L->rawset(i, j, nz);
    }
  }

  nKillChar(Q);
  Print("bimFarey worked\n");
#if 0
  L->Print();
  Print("\n * 1/");
  n_Print(den, Z);
  Print("\n");
#endif
  return den;
}

bigintmat* bimMult	(	bigintmat *	a,
		bigintmat *	b
	)

Definition at line 251 of file bigintmat.cc.

{
  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);

        number sum2 = n_Add(sum, prod, basecoeffs); // no inplace add :(

        n_Delete(&sum, basecoeffs); n_Delete(&prod, basecoeffs);

        sum = sum2;
      }
      bim->rawset(i, j, sum, basecoeffs);
    }
  return bim;
}

bigintmat* bimMult	(	bigintmat *	a,
		int	b
	)

Definition at line 299 of file bigintmat.cc.

{

  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
	)

Definition at line 318 of file bigintmat.cc.

{
  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;
}

void bimMult	(	bigintmat *	a,
		bigintmat *	b,
		bigintmat *	c
	)

Multipliziert Matrix a und b und speichert Ergebnis in c.

Definition at line 1818 of file bigintmat.cc.

{
  if (!nCoeffs_are_equal(a->basecoeffs(), b->basecoeffs())) {
    Werror("Error in bimMult. Coeffs do not agree!");
    return;
  }
  if ((a->rows() != c->rows()) || (b->cols() != c->cols()) || (a->cols() != b->rows())) {
    Werror("Error in bimMult. Dimensions do not agree!");
    return;
  }
  bigintmat *tmp = bimMult(a, b);
  c->copy(tmp);

  delete tmp;
}

bigintmat* bimSub	(	bigintmat *	a,
		bigintmat *	b
	)

Definition at line 214 of file bigintmat.cc.

{
  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,
		int	b
	)

Definition at line 232 of file bigintmat.cc.

{
  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_Sub((*a)[i], bb, basecoeffs), basecoeffs);

  n_Delete(&bb,basecoeffs);
  return bim;
}

void diagonalForm	(	bigintmat *	A,
		bigintmat **	S,
		bigintmat **	T
	)

Definition at line 2301 of file bigintmat.cc.

{
  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();
      Print("\n");
#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
    Print("Diag ? %d\n", diag);
    a->Print();
    Print("\n");
#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);
}

static int findLongest ( int *  a, int  length  )

static

Definition at line 533 of file bigintmat.cc.

{
  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  )

static

Definition at line 548 of file bigintmat.cc.

{
  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;
}

static int intArrSum ( int *  a, int  length  )

static

Definition at line 525 of file bigintmat.cc.

{
  int sum = 0;
  for (int i=0; i<length; i++)
    sum += a[i];
  return sum;
}

bigintmat* iv2bim	(	intvec *	b,
		const coeffs	C
	)

Definition at line 347 of file bigintmat.cc.

{
  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 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 at line 2405 of file bigintmat.cc.

                                                              {
#if 0
  Print("Kernel of ");
  a->Print();
  Print(" modulo ");
  n_Print(p, q);
  Print("\n");
#endif

  coeffs coe = numbercoeffs(p, q);
  bigintmat *m = bimChangeCoeff(a, coe), *U, *V;
  diagonalForm(m, &U, &V);
#if 0
  Print("\ndiag form: ");
  m->Print();
  Print("\nU:\n");
  U->Print();
  Print("\nV:\n");
  V->Print();
  Print("\n");
#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++);

#undef MAX
#define MAX(a,b) (a > b ? a : b)
  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
	)

Definition at line 2449 of file bigintmat.cc.

                                           {
  if ((r == NULL) || (s == NULL))
    return false;
  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;
}

static coeffs numbercoeffs ( number  n, coeffs  c  )

static

create Z/nA of type n_Zn

Definition at line 21 of file bigintmat.cc.

{
  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;
}

bool operator!=	(	const bigintmat &	lhr,
		const bigintmat &	rhr
	)

Definition at line 172 of file bigintmat.cc.

{
  return !(lhr==rhr);
}

bool operator==	(	const bigintmat &	lhr,
		const bigintmat &	rhr
	)

Definition at line 155 of file bigintmat.cc.

{
  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;
}

static bigintmat* prependIdentity ( bigintmat *  A )

static

Definition at line 1907 of file bigintmat.cc.

{
  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 void reduce_mod_howell ( bigintmat *  A, bigintmat *  b, bigintmat *  eps, bigintmat *  x  )

static

Definition at line 1834 of file bigintmat.cc.

                                                                                         {
  //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
  Print("reduce_mod_howell: A:\n");
  A->Print();
  Print("\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
    Print("\nx:\n");
    x->Print();
    Print("\neps:\n");
    eps->Print();
    Print("\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
  Print("\nx:\n");
  x->Print();
  Print("\neps:\n");
  eps->Print();
  Print("\n****************************************\n");
#endif
}

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 d is returned. Currently available for Z, Q and Z/nZ (and possibly for all fields: d=1 there) Beware that the internal functions can find the kernel as well - but the interface is lacking.

Definition at line 2263 of file bigintmat.cc.

                                                         {
#if 0
  Print("Solve Ax=b for A=\n");
  A->Print();
  Print("\nb = \n");
  b->Print();
  Print("\nx = \n");
  x->Print();
  Print("\n");
#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:
      Warn("have field, should use Gauss or better");
    default:
      if (R->cfXExtGcd && R->cfAnn) { //assume it's Euclidean
        return solveAx_howell(A, b, x, NULL);
      }
      Werror("have no solve algorithm");
  }
  return NULL;
}

static number solveAx_dixon ( bigintmat *  A, bigintmat *  B, bigintmat *  x, bigintmat *  kern  )

static

Definition at line 1973 of file bigintmat.cc.

                                                                                       {
  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()) {
      Print("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;
}

static number solveAx_howell ( bigintmat *  A, bigintmat *  b, bigintmat *  x, bigintmat *  kern  )

static

Definition at line 2150 of file bigintmat.cc.

                                                                                        {
  // 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
  Print("sol = 1/");
  n_Print(den, R);
  Print(" *\n");
  x->Print();
  Print("\n");
#endif
  return den;
}