matpol.cc File Reference

#include "misc/auxiliary.h"
#include "misc/mylimits.h"
#include "misc/intvec.h"
#include "coeffs/numbers.h"
#include "reporter/reporter.h"
#include "monomials/ring.h"
#include "monomials/p_polys.h"
#include "simpleideals.h"
#include "matpol.h"
#include "prCopy.h"
#include "clapsing.h"
#include "sparsmat.h"

Go to the source code of this file.

Data Structures
class	row_col_weight
class	mp_permmatrix

Functions
static poly	mp_Exdiv (poly m, poly d, poly vars, const ring)
static poly	mp_Select (poly fro, poly what, const ring)
static poly	mp_SelectId (ideal I, poly what, const ring R)
matrix	mpNew (int r, int c)
	create a r x c zero-matrix More...
matrix	mp_Copy (matrix a, const ring r)
	copies matrix a (from ring r to r) More...
matrix	mp_Copy (const matrix a, const ring rSrc, const ring rDst)
	copies matrix a from rSrc into rDst More...
matrix	mp_InitP (int r, int c, poly p, const ring R)
	make it a p * unit matrix More...
matrix	mp_InitI (int r, int c, int v, const ring R)
	make it a v * unit matrix More...
matrix	mp_MultI (matrix a, long f, const ring R)
	c = f*a More...
matrix	mp_MultP (matrix a, poly p, const ring R)
	multiply a matrix 'a' by a poly 'p', destroy the args More...
matrix	pMultMp (poly p, matrix a, const ring R)
matrix	mp_Add (matrix a, matrix b, const ring R)
matrix	mp_Sub (matrix a, matrix b, const ring R)
matrix	mp_Mult (matrix a, matrix b, const ring R)
matrix	mp_Transp (matrix a, const ring R)
poly	mp_Trace (matrix a, const ring R)
poly	TraceOfProd (matrix a, matrix b, int n, const ring R)
matrix	mp_Coeffs (ideal I, int var, const ring R)
	corresponds to Maple's coeffs: var has to be the number of a variable More...
void	mp_Monomials (matrix c, int r, int var, matrix m, const ring R)
matrix	mp_CoeffProc (poly f, poly vars, const ring R)
matrix	mp_CoeffProcId (ideal I, poly vars, const ring R)
void	mp_Coef2 (poly v, poly mon, matrix *c, matrix *m, const ring R)
	corresponds to Macauley's coef: the exponent vector of vars has to contain the variables, eg 'xy'; then the poly f is searched for monomials in x and y, these monimials are written to the first row of the matrix co. the second row of co contains the respective factors in f. Thus f = sum co[1,i]*co[2,i], i = 1..cols, rows equals 2. More...
int	mp_Compare (matrix a, matrix b, const ring R)
BOOLEAN	mp_Equal (matrix a, matrix b, const ring R)
static poly	p_Insert (poly p1, poly p2, const ring R)
static void	mp_PartClean (matrix a, int lr, int lc, const ring R)
BOOLEAN	mp_IsDiagUnit (matrix U, const ring R)
void	iiWriteMatrix (matrix im, const char *n, int dim, const ring r, int spaces)
	set spaces to zero by default More...
char *	iiStringMatrix (matrix im, int dim, const ring r, char ch)
void	mp_Delete (matrix *a, const ring r)
static float	mp_PolyWeight (poly p, const ring r)
static void	mpReplace (int j, int n, int &sign, int *perm)
static int	mp_PivBar (matrix a, int lr, int lc, const ring r)
static void	mpSwapRow (matrix a, int pos, int lr, int lc)
static int	mp_PrepareRow (matrix a, int lr, int lc, const ring R)
static int	mp_PivRow (matrix a, int lr, int lc, const ring r)
static void	mpSwapCol (matrix a, int pos, int lr, int lc)
static int	mp_PreparePiv (matrix a, int lr, int lc, const ring r)
static void	mp_ElimBar (matrix a0, matrix re, poly div, int lr, int lc, const ring R)
void	mp_MinorToResult (ideal result, int &elems, matrix a, int r, int c, ideal R, const ring)
	entries of a are minors and go to result (only if not in R) More...
static void	mpFinalClean (matrix a)
void	mp_RecMin (int ar, ideal result, int &elems, matrix a, int lr, int lc, poly barDiv, ideal R, const ring r)
	produces recursively the ideal of all arxar-minors of a More...
poly	mp_DetBareiss (matrix a, const ring r)
	returns the determinant of the matrix m; uses Bareiss algorithm More...
matrix	mp_Wedge (matrix a, int ar, const ring R)
static ideal	sm_MultAndShift (poly f, ideal B, int s, const ring r)
static void	sm_AddSubMat (ideal res, ideal sm, int col, const ring r)
ideal	sm_Tensor (ideal A, ideal B, const ring r)
ideal	sm_Add (ideal a, ideal b, const ring R)
ideal	sm_Sub (ideal a, ideal b, const ring R)
ideal	sm_Mult (ideal a, ideal b, const ring R)
ideal	sm_Flatten (ideal a, const ring R)
ideal	sm_UnFlatten (ideal a, int col, const ring R)
poly	sm_Trace (ideal a, const ring R)
int	sm_Compare (ideal a, ideal b, const ring R)
BOOLEAN	sm_Equal (ideal a, ideal b, const ring R)
static matrix	mu (matrix A, const ring R)
poly	mp_DetMu (matrix A, const ring R)
DetVariant	mp_GetAlgorithmDet (matrix m, const ring r)
DetVariant	mp_GetAlgorithmDet (const char *s)
poly	mp_Det (matrix a, const ring r, DetVariant d)
poly	sm_Det (ideal a, const ring r, DetVariant d)

Function Documentation

iiStringMatrix()

char * iiStringMatrix	(	matrix	im,
		int	dim,
		const ring	r,
		char	ch
	)

Definition at line 848 of file matpol.cc.

{
  int i,ii = MATROWS(im);
  int j,jj = MATCOLS(im);
  poly *pp = im->m;
  char ch_s[2];
  ch_s[0]=ch;
  ch_s[1]='\0';
 
  StringSetS("");
 
  for (i=0; i<ii; i++)
  {
    for (j=0; j<jj; j++)
    {
      p_String0(*pp++, r);
      StringAppendS(ch_s);
      if (dim > 1) StringAppendS("\n");
    }
  }
  char *s=StringEndS();
  s[strlen(s)- (dim > 1 ? 2 : 1)]='\0';
  return s;
}

iiWriteMatrix()

void iiWriteMatrix	(	matrix	im,
		const char *	n,
		int	dim,
		const ring	r,
		int	spaces
	)

set spaces to zero by default

Definition at line 827 of file matpol.cc.

{
  int i,ii = MATROWS(im)-1;
  int j,jj = MATCOLS(im)-1;
  poly *pp = im->m;
 
  for (i=0; i<=ii; i++)
  {
    for (j=0; j<=jj; j++)
    {
      if (spaces>0)
        Print("%-*.*s",spaces,spaces," ");
      if (dim == 2) Print("%s[%u,%u]=",n,i+1,j+1);
      else if (dim == 1) Print("%s[%u]=",n,j+1);
      else if (dim == 0) Print("%s=",n);
      if ((i<ii)||(j<jj)) p_Write(*pp++, r);
      else                p_Write0(*pp, r);
    }
  }
}

mp_Add()

matrix mp_Add	(	matrix	a,
		matrix	b,
		const ring	R
	)

Definition at line 172 of file matpol.cc.

{
  int k, n = a->nrows, m = a->ncols;
  if ((n != b->nrows) || (m != b->ncols))
  {
/*
*    Werror("cannot add %dx%d matrix and %dx%d matrix",
*      m,n,b->cols(),b->rows());
*/
    return NULL;
  }
  matrix c = mpNew(n,m);
  for (k=m*n-1; k>=0; k--)
    c->m[k] = p_Add_q(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
  return c;
}

mp_Coef2()

void mp_Coef2	(	poly	v,
		poly	mon,
		matrix *	c,
		matrix *	m,
		const ring	R
	)

corresponds to Macauley's coef: the exponent vector of vars has to contain the variables, eg 'xy'; then the poly f is searched for monomials in x and y, these monimials are written to the first row of the matrix co. the second row of co contains the respective factors in f. Thus f = sum co[1,i]*co[2,i], i = 1..cols, rows equals 2.

Definition at line 574 of file matpol.cc.

{
  poly* s;
  poly p;
  int sl,i,j;
  int l=0;
  poly sel=mp_Select(v,mon, R);
 
  p_Vec2Polys(sel,&s,&sl, R);
  for (i=0; i<sl; i++)
    l=si_max(l,pLength(s[i]));
  *c=mpNew(sl,l);
  *m=mpNew(sl,l);
  poly h;
  int isConst;
  for (j=1; j<=sl;j++)
  {
    p=s[j-1];
    if (p_IsConstant(p, R)) /*p != NULL */
    {
      isConst=-1;
      i=l;
    }
    else
    {
      isConst=1;
      i=1;
    }
    while(p!=NULL)
    {
      h = p_Head(p, R);
      MATELEM(*m,j,i) = h;
      i+=isConst;
      p = p->next;
    }
  }
  while (v!=NULL)
  {
    i = 1;
    j = __p_GetComp(v, R);
    loop
    {
      poly mp=MATELEM(*m,j,i);
      if (mp!=NULL)
      {
        h = mp_Exdiv(v, mp /*MATELEM(*m,j,i)*/, mp, R);
        if (h!=NULL)
        {
          p_SetComp(h,0, R);
          MATELEM(*c,j,i) = p_Add_q(MATELEM(*c,j,i), h, R);
          break;
        }
      }
      if (i < l)
        i++;
      else
        break;
    }
    v = v->next;
  }
}

mp_CoeffProc()

matrix mp_CoeffProc	(	poly	f,
		poly	vars,
		const ring	R
	)

Definition at line 392 of file matpol.cc.

{
  assume(vars!=NULL);
  poly sel, h;
  int l, i;
  int pos_of_1 = -1;
  matrix co;
 
  if (f==NULL)
  {
    co = mpNew(2, 1);
    MATELEM(co,1,1) = p_One(R);
    //MATELEM(co,2,1) = NULL;
    return co;
  }
  sel = mp_Select(f, vars, R);
  l = pLength(sel);
  co = mpNew(2, l);
 
  if (rHasLocalOrMixedOrdering(R))
  {
    for (i=l; i>=1; i--)
    {
      h = sel;
      pIter(sel);
      pNext(h)=NULL;
      MATELEM(co,1,i) = h;
      //MATELEM(co,2,i) = NULL;
      if (p_IsConstant(h, R)) pos_of_1 = i;
    }
  }
  else
  {
    for (i=1; i<=l; i++)
    {
      h = sel;
      pIter(sel);
      pNext(h)=NULL;
      MATELEM(co,1,i) = h;
      //MATELEM(co,2,i) = NULL;
      if (p_IsConstant(h, R)) pos_of_1 = i;
    }
  }
  while (f!=NULL)
  {
    i = 1;
    loop
    {
      if (i!=pos_of_1)
      {
        h = mp_Exdiv(f, MATELEM(co,1,i),vars, R);
        if (h!=NULL)
        {
          MATELEM(co,2,i) = p_Add_q(MATELEM(co,2,i), h, R);
          break;
        }
      }
      if (i == l)
      {
        // check monom 1 last:
        if (pos_of_1 != -1)
        {
          h = mp_Exdiv(f, MATELEM(co,1,pos_of_1),vars, R);
          if (h!=NULL)
          {
            MATELEM(co,2,pos_of_1) = p_Add_q(MATELEM(co,2,pos_of_1), h, R);
          }
        }
        break;
      }
      i ++;
    }
    pIter(f);
  }
  return co;
}

mp_CoeffProcId()

matrix mp_CoeffProcId	(	ideal	I,
		poly	vars,
		const ring	R
	)

Definition at line 469 of file matpol.cc.

{
  assume(vars!=NULL);
  poly sel, h;
  int l, i;
  int pos_of_1 = -1;
  matrix co;
 
  if (idIs0(I))
  {
    co = mpNew(IDELEMS(I)+1,1);
    MATELEM(co,1,1) = p_One(R);
    return co;
  }
  sel = mp_SelectId(I, vars, R);
  l = pLength(sel);
  co = mpNew(IDELEMS(I)+1, l);
 
  if (rHasLocalOrMixedOrdering(R))
  {
    for (i=l; i>=1; i--)
    {
      h = sel;
      pIter(sel);
      pNext(h)=NULL;
      MATELEM(co,1,i) = h;
      //MATELEM(co,2,i) = NULL;
      if (p_IsConstant(h, R)) pos_of_1 = i;
    }
  }
  else
  {
    for (i=1; i<=l; i++)
    {
      h = sel;
      pIter(sel);
      pNext(h)=NULL;
      MATELEM(co,1,i) = h;
      //MATELEM(co,2,i) = NULL;
      if (p_IsConstant(h, R)) pos_of_1 = i;
    }
  }
  for(int j=0;j<IDELEMS(I);j++)
  {
    poly f=I->m[j];
    while (f!=NULL)
    {
      i = 1;
      loop
      {
        if (i!=pos_of_1)
        {
          h = mp_Exdiv(f, MATELEM(co,1,i),vars, R);
          if (h!=NULL)
          {
            MATELEM(co,j+2,i) = p_Add_q(MATELEM(co,j+2,i), h, R);
            break;
          }
        }
        if (i == l)
        {
          // check monom 1 last:
          if (pos_of_1 != -1)
          {
            h = mp_Exdiv(f, MATELEM(co,1,pos_of_1),vars, R);
            if (h!=NULL)
            {
              MATELEM(co,j+2,pos_of_1) = p_Add_q(MATELEM(co,j+2,pos_of_1), h, R);
            }
          }
          break;
        }
        i ++;
      }
      pIter(f);
    }
  }
  return co;
}

mp_Coeffs()

matrix mp_Coeffs	(	ideal	I,
		int	var,
		const ring	R
	)

corresponds to Maple's coeffs: var has to be the number of a variable

Definition at line 306 of file matpol.cc.

{
  poly h,f;
  int l, i, c, m=0;
  /* look for maximal power m of x_var in I */
  for (i=IDELEMS(I)-1; i>=0; i--)
  {
    f=I->m[i];
    while (f!=NULL)
    {
      l=p_GetExp(f,var, R);
      if (l>m) m=l;
      pIter(f);
    }
  }
  matrix co=mpNew((m+1)*I->rank,IDELEMS(I));
  /* divide each monomial by a power of x_var,
  * remember the power in l and the component in c*/
  for (i=IDELEMS(I)-1; i>=0; i--)
  {
    f=I->m[i];
    I->m[i]=NULL;
    while (f!=NULL)
    {
      l=p_GetExp(f,var, R);
      p_SetExp(f,var,0, R);
      c=si_max((int)p_GetComp(f, R),1);
      p_SetComp(f,0, R);
      p_Setm(f, R);
      /* now add the resulting monomial to co*/
      h=pNext(f);
      pNext(f)=NULL;
      //MATELEM(co,c*(m+1)-l,i+1)
      //  =p_Add_q(MATELEM(co,c*(m+1)-l,i+1),f, R);
      MATELEM(co,(c-1)*(m+1)+l+1,i+1)
        =p_Add_q(MATELEM(co,(c-1)*(m+1)+l+1,i+1),f, R);
      /* iterate f*/
      f=h;
    }
  }
  id_Delete(&I, R);
  return co;
}

mp_Compare()

int mp_Compare	(	matrix	a,
		matrix	b,
		const ring	R
	)

Definition at line 636 of file matpol.cc.

{
  if (MATCOLS(a)<MATCOLS(b)) return -1;
  else if (MATCOLS(a)>MATCOLS(b)) return 1;
  if (MATROWS(a)<MATROWS(b)) return -1;
  else if (MATROWS(a)<MATROWS(b)) return 1;
 
  unsigned ii=MATCOLS(a)*MATROWS(a)-1;
  unsigned j=0;
  int r=0;
  while (j<=ii)
  {
    r=p_Compare(a->m[j],b->m[j],R);
    if (r!=0) return r;
    j++;
  }
  return r;
}

mp_Copy() [1/2]

matrix mp_Copy	(	const matrix	a,
		const ring	rSrc,
		const ring	rDst
	)

copies matrix a from rSrc into rDst

Definition at line 78 of file matpol.cc.

{
  id_Test((ideal)a, rSrc);
 
  poly t;
  int i, m=MATROWS(a), n=MATCOLS(a);
 
  matrix b = mpNew(m, n);
 
  for (i=m*n-1; i>=0; i--)
  {
    t = a->m[i];
    if (t!=NULL)
    {
      b->m[i] = prCopyR_NoSort(t, rSrc, rDst);
      p_Normalize(b->m[i], rDst);
    }
  }
  b->rank=a->rank;
 
  id_Test((ideal)b, rDst);
 
  return b;
}

mp_Copy() [2/2]

matrix mp_Copy	(	matrix	a,
		const ring	r
	)

copies matrix a (from ring r to r)

Definition at line 57 of file matpol.cc.

{
  id_Test((ideal)a, r);
  poly t;
  int m=MATROWS(a), n=MATCOLS(a);
  matrix b = mpNew(m, n);
 
  for (long i=(long)m*(long)n-1; i>=0; i--)
  {
    t = a->m[i];
    if (t!=NULL)
    {
      p_Normalize(t, r);
      b->m[i] = p_Copy(t, r);
    }
  }
  b->rank=a->rank;
  return b;
}

mp_Delete()

void mp_Delete	(	matrix *	a,
		const ring	r
	)

Definition at line 873 of file matpol.cc.

{
  id_Delete((ideal *) a, r);
}

mp_Det()

poly mp_Det	(	matrix	a,
		const ring	r,
		DetVariant	d
	)

Definition at line 2136 of file matpol.cc.

{
  if ((MATCOLS(a)==0)
  && (MATROWS(a)==0))
    return p_One(r);
  if (d==DetDefault) d=mp_GetAlgorithmDet(a,r);
  switch (d)
  {
    case DetBareiss: return mp_DetBareiss(a,r);
    case DetMu: return mp_DetMu(a,r);
    case DetFactory: return singclap_det(a,r);
    case DetSBareiss:
    {
      ideal I=id_Matrix2Module(mp_Copy(a, r),r);
      poly p=sm_CallDet(I, r);
      id_Delete(&I, r);
      return p;
    }
    default:
      WerrorS("unknown algorithm for det");
      return NULL;
  }
}

mp_DetBareiss()

poly mp_DetBareiss	(	matrix	a,
		const ring	r
	)

returns the determinant of the matrix m; uses Bareiss algorithm

Definition at line 1669 of file matpol.cc.

{
  int s;
  poly div, res;
  if (MATROWS(a) != MATCOLS(a))
  {
    Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
    return NULL;
  }
  matrix c = mp_Copy(a,r);
  mp_permmatrix *Bareiss = new mp_permmatrix(c,r);
  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
 
  /* Bareiss */
  div = NULL;
  while(Bareiss->mpPivotBareiss(&w))
  {
    Bareiss->mpElimBareiss(div);
    div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
  }
  Bareiss->mpRowReorder();
  Bareiss->mpColReorder();
  Bareiss->mpSaveArray();
  s = Bareiss->mpGetSign();
  delete Bareiss;
 
  /* result */
  res = MATELEM(c,1,1);
  MATELEM(c,1,1) = NULL;
  id_Delete((ideal *)&c,r);
  if (s < 0)
    res = p_Neg(res,r);
  return res;
}

mp_DetMu()

poly mp_DetMu	(	matrix	A,
		const ring	R
	)

Definition at line 2063 of file matpol.cc.

{
  int n=MATROWS(A);
  assume(MATCOLS(A)==n);
    /*
    *   Intput:
    *   int n: Dimension der Matrix
    *   int A: n*n Matrix
    *
    *   Berechnet n-1 mal: X = mu(X)*A
    *
    *   Output: det(A)
    */
 
    //speichere A ab:
    matrix workA=mp_Copy(A,R);
 
    // berechen X = mu(X)*A
    matrix X;
    for (int i = n-1; i >0; i--)
    {
        X=mu(workA,R);
        id_Delete((ideal*)&workA,R);
        workA=mp_Mult(X,A,R);
        id_Delete((ideal*)&X,R);
    }
 
    // berrechne det(A)
    poly res;
    if (n%2 == 0)
    {
        res=p_Neg(MATELEM0(workA,0,0),R);
    }
    else
    {
        res=MATELEM0(workA,0,0);
    }
    MATELEM0(workA,0,0)=NULL;
    id_Delete((ideal*)&workA,R);
    return res;
}

mp_ElimBar()

static void mp_ElimBar ( matrix  a0, matrix  re, poly  div, int  lr, int  lc, const ring  R  )

static

Definition at line 1442 of file matpol.cc.

{
  int r=lr-1, c=lc-1;
  poly *b = a0->m, *x = re->m;
  poly piv, elim, q1, *ap, *a, *q;
  int i, j;
 
  ap = &b[r*a0->ncols];
  piv = ap[c];
  for(j=c-1; j>=0; j--)
    if (ap[j] != NULL) ap[j] = p_Neg(ap[j],R);
  for(i=r-1; i>=0; i--)
  {
    a = &b[i*a0->ncols];
    q = &x[i*re->ncols];
    if (a[c] != NULL)
    {
      elim = a[c];
      for (j=c-1; j>=0; j--)
      {
        q1 = NULL;
        if (a[j] != NULL)
        {
          q1 = sm_MultDiv(a[j], piv, div,R);
          if (ap[j] != NULL)
          {
            poly q2 = sm_MultDiv(ap[j], elim, div, R);
            q1 = p_Add_q(q1,q2,R);
          }
        }
        else if (ap[j] != NULL)
          q1 = sm_MultDiv(ap[j], elim, div, R);
        if (q1 != NULL)
        {
          if (div)
            sm_SpecialPolyDiv(q1, div,R);
          q[j] = q1;
        }
      }
    }
    else
    {
      for (j=c-1; j>=0; j--)
      {
        if (a[j] != NULL)
        {
          q1 = sm_MultDiv(a[j], piv, div, R);
          if (div)
            sm_SpecialPolyDiv(q1, div, R);
          q[j] = q1;
        }
      }
    }
  }
}

mp_Equal()

BOOLEAN mp_Equal	(	matrix	a,
		matrix	b,
		const ring	R
	)

Definition at line 655 of file matpol.cc.

{
  if ((MATCOLS(a)!=MATCOLS(b)) || (MATROWS(a)!=MATROWS(b)))
    return FALSE;
  int i=MATCOLS(a)*MATROWS(a)-1;
  while (i>=0)
  {
    if (a->m[i]==NULL)
    {
      if (b->m[i]!=NULL) return FALSE;
    }
    else if (b->m[i]==NULL) return FALSE;
    else if (p_Cmp(a->m[i],b->m[i], R)!=0) return FALSE;
    i--;
  }
  i=MATCOLS(a)*MATROWS(a)-1;
  while (i>=0)
  {
    if(!p_EqualPolys(a->m[i],b->m[i], R)) return FALSE;
    i--;
  }
  return TRUE;
}

mp_Exdiv()

static poly mp_Exdiv ( poly  m, poly  d, poly  vars, const ring  R  )

static

Definition at line 554 of file matpol.cc.

{
  int i;
  poly h = p_Head(m, R);
  for (i=1; i<=rVar(R); i++)
  {
    if (p_GetExp(vars,i, R) > 0)
    {
      if (p_GetExp(d,i, R) != p_GetExp(h,i, R))
      {
        p_Delete(&h, R);
        return NULL;
      }
      p_SetExp(h,i,0, R);
    }
  }
  p_Setm(h, R);
  return h;
}

mp_GetAlgorithmDet() [1/2]

DetVariant mp_GetAlgorithmDet

(

const char *

s

)

Definition at line 2125 of file matpol.cc.

{
  if (strcmp(s,"Bareiss")==0) return DetBareiss;
  if (strcmp(s,"SBareiss")==0) return DetSBareiss;
  if (strcmp(s,"Mu")==0) return DetMu;
  if (strcmp(s,"Factory")==0) return DetFactory;
  WarnS("unknown method for det");
  return DetDefault;
}

mp_GetAlgorithmDet() [2/2]

DetVariant mp_GetAlgorithmDet	(	matrix	m,
		const ring	r
	)

Definition at line 2105 of file matpol.cc.

{
  if (MATROWS(m)+2*r->N>20+5*rField_is_Zp(r)) return DetMu;
  if (MATROWS(m)<10+5*rField_is_Zp(r)) return DetSBareiss;
  BOOLEAN isConst=TRUE;
  int s=0;
  for(int i=MATCOLS(m)*MATROWS(m)-1;i>=0;i--)
  {
    poly p=m->m[i];
    if (p!=NULL)
    {
      if(!p_IsConstant(p,r)) isConst=FALSE;
      s++;
    }
  }
  if (isConst && rField_is_Q(r)) return DetFactory;
  if (s*2<MATCOLS(m)*MATROWS(m)) // few entries
    return DetSBareiss;
  return DetMu;
}

mp_InitI()

matrix mp_InitI	(	int	r,
		int	c,
		int	v,
		const ring	R
	)

make it a v * unit matrix

Definition at line 122 of file matpol.cc.

{
  return mp_InitP(r, c, p_ISet(v, R), R);
}

mp_InitP()

matrix mp_InitP	(	int	r,
		int	c,
		poly	p,
		const ring	R
	)

make it a p * unit matrix

Definition at line 106 of file matpol.cc.

{
  matrix rc = mpNew(r,c);
  int i=si_min(r,c), n = c*(i-1)+i-1, inc = c+1;
 
  p_Normalize(p, R);
  while (n>0)
  {
    rc->m[n] = p_Copy(p, R);
    n -= inc;
  }
  rc->m[0]=p;
  return rc;
}

mp_IsDiagUnit()

BOOLEAN mp_IsDiagUnit	(	matrix	U,
		const ring	R
	)

Definition at line 809 of file matpol.cc.

{
  if(MATROWS(U)!=MATCOLS(U))
    return FALSE;
  for(int i=MATCOLS(U);i>=1;i--)
  {
    for(int j=MATCOLS(U); j>=1; j--)
    {
      if (i==j)
      {
        if (!p_IsUnit(MATELEM(U,i,i), R)) return FALSE;
      }
      else if (MATELEM(U,i,j)!=NULL) return FALSE;
    }
  }
  return TRUE;
}

mp_MinorToResult()

void mp_MinorToResult	(	ideal	result,
		int &	elems,
		matrix	a,
		int	r,
		int	c,
		ideal	R,
		const	ring
	)

entries of a are minors and go to result (only if not in R)

Definition at line 1500 of file matpol.cc.

{
  poly *q1;
  int e=IDELEMS(result);
  int i,j;
 
  if (R != NULL)
  {
    for (i=r-1;i>=0;i--)
    {
      q1 = &(a->m)[i*a->ncols];
      //for (j=c-1;j>=0;j--)
      //{
      //  if (q1[j]!=NULL) q1[j] = kNF(R,currRing->qideal,q1[j]);
      //}
    }
  }
  for (i=r-1;i>=0;i--)
  {
    q1 = &(a->m)[i*a->ncols];
    for (j=c-1;j>=0;j--)
    {
      if (q1[j]!=NULL)
      {
        if (elems>=e)
        {
          pEnlargeSet(&(result->m),e,e);
          e += e;
          IDELEMS(result) =e;
        }
        result->m[elems] = q1[j];
        q1[j] = NULL;
        elems++;
      }
    }
  }
}

mp_Monomials()

void mp_Monomials	(	matrix	c,
		int	r,
		int	var,
		matrix	m,
		const ring	R
	)

Definition at line 355 of file matpol.cc.

{
  /* clear contents of m*/
  int k,l;
  for (k=MATROWS(m);k>0;k--)
  {
    for(l=MATCOLS(m);l>0;l--)
    {
      p_Delete(&MATELEM(m,k,l), R);
    }
  }
  omfreeSize((ADDRESS)m->m,MATROWS(m)*MATCOLS(m)*sizeof(poly));
  /* allocate monoms in the right size r x MATROWS(c)*/
  m->m=(poly*)omAlloc0(r*MATROWS(c)*sizeof(poly));
  MATROWS(m)=r;
  MATCOLS(m)=MATROWS(c);
  m->rank=r;
  /* the maximal power p of x_var: MATCOLS(m)=r*(p+1) */
  int p=MATCOLS(m)/r-1;
  /* fill in the powers of x_var=h*/
  poly h=p_One(R);
  for(k=r;k>0; k--)
  {
    MATELEM(m,k,k*(p+1))=p_One(R);
  }
  for(l=p;l>=0; l--)
  {
    p_SetExp(h,var,p-l, R);
    p_Setm(h, R);
    for(k=r;k>0; k--)
    {
      MATELEM(m,k,k*(p+1)-l)=p_Copy(h, R);
    }
  }
  p_Delete(&h, R);
}

mp_Mult()

matrix mp_Mult	(	matrix	a,
		matrix	b,
		const ring	R
	)

Definition at line 206 of file matpol.cc.

{
  int i, j, k;
  int m = MATROWS(a);
  int p = MATCOLS(a);
  int q = MATCOLS(b);
 
  if (p!=MATROWS(b))
  {
/*
*   Werror("cannot multiply %dx%d matrix and %dx%d matrix",
*     m,p,b->rows(),q);
*/
    return NULL;
  }
  matrix c = mpNew(m,q);
 
  for (i=0; i<m; i++)
  {
    for (k=0; k<p; k++)
    {
      poly aik;
      if ((aik=MATELEM0(a,i,k))!=NULL)
      {
        for (j=0; j<q; j++)
        {
          poly bkj;
          if ((bkj=MATELEM0(b,k,j))!=NULL)
          {
            poly *cij=&(MATELEM0(c,i,j));
            poly s = pp_Mult_qq(aik /*MATELEM0(a,i,k)*/, bkj/*MATELEM0(b,k,j)*/, R);
            (*cij)/*MATELEM0(c,i,j)*/ = p_Add_q((*cij) /*MATELEM0(c,i,j)*/ ,s, R);
          }
        }
      }
    }
  }
  for(i=m*q-1;i>=0;i--) p_Normalize(c->m[i], R);
  return c;
}

mp_MultI()

matrix mp_MultI	(	matrix	a,
		long	f,
		const ring	R
	)

c = f*a

Definition at line 128 of file matpol.cc.

{
  int k, n = a->nrows, m = a->ncols;
  poly p = p_ISet(f, R);
  matrix c = mpNew(n,m);
 
  for (k=m*n-1; k>0; k--)
    c->m[k] = pp_Mult_qq(a->m[k], p, R);
  c->m[0] = p_Mult_q(p_Copy(a->m[0], R), p, R);
  return c;
}

mp_MultP()

matrix mp_MultP	(	matrix	a,
		poly	p,
		const ring	R
	)

multiply a matrix 'a' by a poly 'p', destroy the args

Definition at line 141 of file matpol.cc.

{
  int k, n = a->nrows, m = a->ncols;
 
  p_Normalize(p, R);
  for (k=m*n-1; k>0; k--)
  {
    if (a->m[k]!=NULL)
      a->m[k] = p_Mult_q(a->m[k], p_Copy(p, R), R);
  }
  a->m[0] = p_Mult_q(a->m[0], p, R);
  return a;
}

mp_PartClean()

static void mp_PartClean ( matrix  a, int  lr, int  lc, const ring  R  )

static

Definition at line 797 of file matpol.cc.

{
  poly *q1;
  int i,j;
 
  for (i=lr-1;i>=0;i--)
  {
    q1 = &(a->m)[i*a->ncols];
    for (j=lc-1;j>=0;j--) if(q1[j]) p_Delete(&q1[j], R);
  }
}

mp_PivBar()

static int mp_PivBar ( matrix  a, int  lr, int  lc, const ring  r  )

static

Definition at line 1327 of file matpol.cc.

{
  float f1, f2;
  poly *q1;
  int i,j,io;
 
  io = -1;
  f1 = 1.0e30;
  for (i=lr-1;i>=0;i--)
  {
    q1 = &(a->m)[i*a->ncols];
    f2 = 0.0;
    for (j=lc-1;j>=0;j--)
    {
      if (q1[j]!=NULL)
        f2 += mp_PolyWeight(q1[j],r);
    }
    if ((f2!=0.0) && (f2<f1))
    {
      f1 = f2;
      io = i;
    }
  }
  if (io<0) return 0;
  else return io+1;
}

mp_PivRow()

static int mp_PivRow ( matrix  a, int  lr, int  lc, const ring  r  )

static

Definition at line 1387 of file matpol.cc.

{
  float f1, f2;
  poly *q1;
  int j,jo;
 
  jo = -1;
  f1 = 1.0e30;
  q1 = &(a->m)[(lr-1)*a->ncols];
  for (j=lc-1;j>=0;j--)
  {
    if (q1[j]!=NULL)
    {
      f2 = mp_PolyWeight(q1[j],r);
      if (f2<f1)
      {
        f1 = f2;
        jo = j;
      }
    }
  }
  if (jo<0) return 0;
  else return jo+1;
}

mp_PolyWeight()

static float mp_PolyWeight ( poly  p, const ring  r  )

static

Definition at line 1295 of file matpol.cc.

{
  int i;
  float res;
 
  if (pNext(p) == NULL)
  {
    res = (float)n_Size(pGetCoeff(p),r->cf);
    for (i=r->N;i>0;i--)
    {
      if(p_GetExp(p,i,r)!=0)
      {
        res += 2.0;
        break;
      }
    }
  }
  else
  {
    res = 0.0;
    do
    {
      res += (float)n_Size(pGetCoeff(p),r->cf)+2.0;
      pIter(p);
    }
    while (p);
  }
  return res;
}

mp_PreparePiv()

static int mp_PreparePiv ( matrix  a, int  lr, int  lc, const ring  r  )

static

Definition at line 1432 of file matpol.cc.

{
  int c;
 
  c = mp_PivRow(a, lr, lc,r);
  if(c==0) return 0;
  if(c<lc) mpSwapCol(a, c, lr, lc);
  return 1;
}

mp_PrepareRow()

static int mp_PrepareRow ( matrix  a, int  lr, int  lc, const ring  R  )

static

Definition at line 1374 of file matpol.cc.

{
  int r;
 
  r = mp_PivBar(a,lr,lc,R);
  if(r==0) return 0;
  if(r<lr) mpSwapRow(a, r, lr, lc);
  return 1;
}

mp_RecMin()

void mp_RecMin	(	int	ar,
		ideal	result,
		int &	elems,
		matrix	a,
		int	lr,
		int	lc,
		poly	barDiv,
		ideal	R,
		const ring	r
	)

produces recursively the ideal of all arxar-minors of a

for minors with Bareiss

Definition at line 1596 of file matpol.cc.

{
  int k;
  int kr=lr-1,kc=lc-1;
  matrix nextLevel=mpNew(kr,kc);
 
  loop
  {
/*--- look for an optimal row and bring it to last position ------------*/
    if(mp_PrepareRow(a,lr,lc,r)==0) break;
/*--- now take all pivots from the last row ------------*/
    k = lc;
    loop
    {
      if(mp_PreparePiv(a,lr,k,r)==0) break;
      mp_ElimBar(a,nextLevel,barDiv,lr,k,r);
      k--;
      if (ar>1)
      {
        mp_RecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R,r);
        mp_PartClean(nextLevel,kr,k, r);
      }
      else mp_MinorToResult(result,elems,nextLevel,kr,k,R,r);
      if (ar>k-1) break;
    }
    if (ar>=kr) break;
/*--- now we have to take out the last row...------------*/
    lr = kr;
    kr--;
  }
  mpFinalClean(nextLevel);
}

mp_Select()

static poly mp_Select ( poly  fro, poly  what, const ring  R  )

static

Definition at line 741 of file matpol.cc.

{
  int i;
  poly h, res;
  res = NULL;
  while (fro!=NULL)
  {
    h = p_One(R);
    for (i=1; i<=rVar(R); i++)
      p_SetExp(h,i, p_GetExp(fro,i, R) * p_GetExp(what, i, R), R);
    p_SetComp(h, p_GetComp(fro, R), R);
    p_Setm(h, R);
    res = p_Insert(h, res, R);
    fro = fro->next;
  }
  return res;
}

mp_SelectId()

static poly mp_SelectId ( ideal  I, poly  what, const ring  R  )

static

Definition at line 759 of file matpol.cc.

{
  int i;
  poly h, res;
  res = NULL;
  for(int j=0;j<IDELEMS(I);j++)
  {
    poly fro=I->m[j];
    while (fro!=NULL)
    {
      h = p_One(R);
      for (i=1; i<=rVar(R); i++)
        p_SetExp(h,i, p_GetExp(fro,i, R) * p_GetExp(what, i, R), R);
      p_SetComp(h, p_GetComp(fro, R), R);
      p_Setm(h, R);
      res = p_Insert(h, res, R);
      fro = fro->next;
    }
  }
  return res;
}

mp_Sub()

matrix mp_Sub	(	matrix	a,
		matrix	b,
		const ring	R
	)

Definition at line 189 of file matpol.cc.

{
  int k, n = a->nrows, m = a->ncols;
  if ((n != b->nrows) || (m != b->ncols))
  {
/*
*    Werror("cannot sub %dx%d matrix and %dx%d matrix",
*      m,n,b->cols(),b->rows());
*/
    return NULL;
  }
  matrix c = mpNew(n,m);
  for (k=m*n-1; k>=0; k--)
    c->m[k] = p_Sub(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
  return c;
}

mp_Trace()

poly mp_Trace	(	matrix	a,
		const ring	R
	)

Definition at line 268 of file matpol.cc.

{
  int i;
  int n = (MATCOLS(a)<MATROWS(a)) ? MATCOLS(a) : MATROWS(a);
  poly  t = NULL;
 
  for (i=1; i<=n; i++)
    t = p_Add_q(t, p_Copy(MATELEM(a,i,i), R), R);
  return t;
}

mp_Transp()

matrix mp_Transp	(	matrix	a,
		const ring	R
	)

Definition at line 247 of file matpol.cc.

{
  int    i, j, r = MATROWS(a), c = MATCOLS(a);
  poly *p;
  matrix b =  mpNew(c,r);
 
  p = b->m;
  for (i=0; i<c; i++)
  {
    for (j=0; j<r; j++)
    {
      if (a->m[j*c+i]!=NULL) *p = p_Copy(a->m[j*c+i], R);
      p++;
    }
  }
  return b;
}

mp_Wedge()

matrix mp_Wedge	(	matrix	a,
		int	ar,
		const ring	R
	)

Definition at line 1744 of file matpol.cc.

{
  int     i,j,k,l;
  int *rowchoise,*colchoise;
  BOOLEAN rowch,colch;
  matrix result;
  matrix tmp;
  poly p;
 
  i = binom(a->nrows,ar);
  j = binom(a->ncols,ar);
 
  rowchoise=(int *)omAlloc(ar*sizeof(int));
  colchoise=(int *)omAlloc(ar*sizeof(int));
  result = mpNew(i,j);
  tmp    = mpNew(ar,ar);
  l = 1; /* k,l:the index in result*/
  idInitChoise(ar,1,a->nrows,&rowch,rowchoise);
  while (!rowch)
  {
    k=1;
    idInitChoise(ar,1,a->ncols,&colch,colchoise);
    while (!colch)
    {
      for (i=1; i<=ar; i++)
      {
        for (j=1; j<=ar; j++)
        {
          MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
        }
      }
      p = mp_DetBareiss(tmp, R);
      if ((k+l) & 1) p=p_Neg(p, R);
      MATELEM(result,l,k) = p;
      k++;
      idGetNextChoise(ar,a->ncols,&colch,colchoise);
    }
    idGetNextChoise(ar,a->nrows,&rowch,rowchoise);
    l++;
  }
 
  /*delete the matrix tmp*/
  for (i=1; i<=ar; i++)
  {
    for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
  }
  id_Delete((ideal *) &tmp, R);
  return (result);
}

mpFinalClean()

static void mpFinalClean ( matrix  a )

static

Definition at line 1588 of file matpol.cc.

{
  omFreeSize((ADDRESS)a->m,a->nrows*a->ncols*sizeof(poly));
  omFreeBin((ADDRESS)a, sip_sideal_bin);
}

mpNew()

matrix mpNew	(	int	r,
		int	c
	)

create a r x c zero-matrix

Definition at line 37 of file matpol.cc.

{
  matrix rc = (matrix)omAllocBin(sip_sideal_bin);
  rc->nrows = r;
  rc->ncols = c;
  rc->rank = r;
  if ((c != 0)&&(r!=0))
  {
    size_t s=((size_t)r)*((size_t)c)*sizeof(poly);
    rc->m = (poly*)omAlloc0(s);
    //if (rc->m==NULL)
    //{
    //  Werror("internal error: creating matrix[%d][%d]",r,c);
    //  return NULL;
    //}
  }
  return rc;
}

mpReplace()

static void mpReplace ( int  j, int  n, int &  sign, int *  perm  )

static

Definition at line 1137 of file matpol.cc.

{
  int k;
 
  if (j != n)
  {
    k = perm[n];
    perm[n] = perm[j];
    perm[j] = k;
    sign = -sign;
  }
}

mpSwapCol()

static void mpSwapCol ( matrix  a, int  pos, int  lr, int  lc  )

static

Definition at line 1412 of file matpol.cc.

{
  poly sw;
  int j;
  poly* a2 = a->m;
  poly* a1 = &a2[pos-1];
 
  a2 = &a2[lc-1];
  for (j=a->ncols*(lr-1); j>=0; j-=a->ncols)
  {
    sw = a1[j];
    a1[j] = a2[j];
    a2[j] = sw;
  }
}

mpSwapRow()

static void mpSwapRow ( matrix  a, int  pos, int  lr, int  lc  )

static

Definition at line 1354 of file matpol.cc.

{
  poly sw;
  int j;
  poly* a2 = a->m;
  poly* a1 = &a2[a->ncols*(pos-1)];
 
  a2 = &a2[a->ncols*(lr-1)];
  for (j=lc-1; j>=0; j--)
  {
    sw = a1[j];
    a1[j] = a2[j];
    a2[j] = sw;
  }
}

mu()

static matrix mu ( matrix  A, const ring  R  )

static

Definition at line 2025 of file matpol.cc.

{
  int n=MATROWS(A);
  assume(MATCOLS(A)==n);
    /*  Die Funktion erstellt die Matrix mu
    *
    *   Input:
    *   int n: Dimension der Matrix
    *   int A: Matrix der Groesse n*n
    *   int X: Speicherplatz fuer Output
    *
    *   In der Matrix X speichert man die Matrix mu
    */
 
    // X als n*n Null-Matrix initalisieren
    matrix X=mpNew(n,n);
 
    //  Diagonaleintraege von X berrechnen
    poly sum = NULL;
    for (int i = n-1; i >= 0; i--)
    {
        MATELEM0(X,i,i) = p_Copy(sum,R);
        sum=p_Sub(sum,p_Copy(MATELEM0(A,i,i),R),R);
    }
    p_Delete(&sum,R);
 
    //  Eintraege aus dem oberen Dreieck von A nach X uebertragen
    for (int i = n-1; i >=0; i--)
    {
        for (int j = i+1; j < n; j++)
        {
            MATELEM0(X,i,j)=p_Copy(MATELEM0(A,i,j),R);
        }
    }
    return X;
}

p_Insert()

static poly p_Insert ( poly  p1, poly  p2, const ring  R  )

static

Definition at line 683 of file matpol.cc.

{
  poly a1, p, a2, a;
  int c;
 
  if (p1==NULL) return p2;
  if (p2==NULL) return p1;
  a1 = p1;
  a2 = p2;
  a = p  = p_One(R);
  loop
  {
    c = p_Cmp(a1, a2, R);
    if (c == 1)
    {
      a = pNext(a) = a1;
      pIter(a1);
      if (a1==NULL)
      {
        pNext(a) = a2;
        break;
      }
    }
    else if (c == -1)
    {
      a = pNext(a) = a2;
      pIter(a2);
      if (a2==NULL)
      {
        pNext(a) = a1;
        break;
      }
    }
    else
    {
      p_LmDelete(&a2, R);
      a = pNext(a) = a1;
      pIter(a1);
      if (a1==NULL)
      {
        pNext(a) = a2;
        break;
      }
      else if (a2==NULL)
      {
        pNext(a) = a1;
        break;
      }
    }
  }
  p_LmDelete(&p, R);
  return p;
}

pMultMp()

matrix pMultMp	(	poly	p,
		matrix	a,
		const ring	R
	)

Definition at line 158 of file matpol.cc.

{
  int k, n = a->nrows, m = a->ncols;
 
  p_Normalize(p, R);
  for (k=m*n-1; k>0; k--)
  {
    if (a->m[k]!=NULL)
      a->m[k] = p_Mult_q(p_Copy(p, R), a->m[k], R);
  }
  a->m[0] = p_Mult_q(p, a->m[0], R);
  return a;
}

sm_Add()

ideal sm_Add	(	ideal	a,
		ideal	b,
		const ring	R
	)

Definition at line 1864 of file matpol.cc.

{
  assume(IDELEMS(a)==IDELEMS(b));
  assume(a->rank==b->rank);
  ideal c=idInit(IDELEMS(a),a->rank);
  for (int k=IDELEMS(a)-1; k>=0; k--)
    c->m[k] = p_Add_q(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
  return c;
}

sm_AddSubMat()

static void sm_AddSubMat ( ideal  res, ideal  sm, int  col, const ring  r  )

static

Definition at line 1815 of file matpol.cc.

{
  for(int i=0;i<IDELEMS(sm);i++)
  {
    res->m[col+i]=p_Add_q(res->m[col+i],sm->m[i],r);
    sm->m[i]=NULL;
  }
}

sm_Compare()

int sm_Compare	(	ideal	a,
		ideal	b,
		const ring	R
	)

Definition at line 1977 of file matpol.cc.

{
  if (IDELEMS(a)<IDELEMS(b)) return -1;
  else if (IDELEMS(a)>IDELEMS(b)) return 1;
  if ((a->rank)<(b->rank)) return -1;
  else if ((a->rank)<(b->rank)) return 1;
 
  unsigned ii=IDELEMS(a)-1;
  unsigned j=0;
  int r=0;
  while (j<=ii)
  {
    r=p_Compare(a->m[j],b->m[j],R);
    if (r!=0) return r;
    j++;
  }
  return r;
}

sm_Det()

poly sm_Det	(	ideal	a,
		const ring	r,
		DetVariant	d
	)

Definition at line 2160 of file matpol.cc.

{
  if ((MATCOLS(a)==0)
  && (MATROWS(a)==0))
    return p_One(r);
  if (d==DetDefault) d=mp_GetAlgorithmDet((matrix)a,r);
  if (d==DetSBareiss) return sm_CallDet(a,r);
  matrix m=id_Module2Matrix(id_Copy(a,r),r);
  poly p=mp_Det(m,r,d);
  id_Delete((ideal *)&m,r);
  return p;
}

sm_Equal()

BOOLEAN sm_Equal	(	ideal	a,
		ideal	b,
		const ring	R
	)

Definition at line 1996 of file matpol.cc.

{
  if ((a->rank!=b->rank) || (IDELEMS(a)!=IDELEMS(b)))
    return FALSE;
  int i=IDELEMS(a)-1;
  while (i>=0)
  {
    if (a->m[i]==NULL)
    {
      if (b->m[i]!=NULL) return FALSE;
    }
    else if (b->m[i]==NULL) return FALSE;
    else if (p_Cmp(a->m[i],b->m[i], R)!=0) return FALSE;
    i--;
  }
  i=IDELEMS(a)-1;
  while (i>=0)
  {
    if(!p_EqualPolys(a->m[i],b->m[i], R)) return FALSE;
    i--;
  }
  return TRUE;
}

sm_Flatten()

ideal sm_Flatten	(	ideal	a,
		const ring	R
	)

Definition at line 1919 of file matpol.cc.

{
  if (IDELEMS(a)==0) return id_Copy(a,R);
  ideal res=idInit(1,IDELEMS(a)*a->rank);
  for(int i=0;i<IDELEMS(a);i++)
  {
    if(a->m[i]!=NULL)
    {
      poly p=p_Copy(a->m[i],R);
      if (i==0) res->m[0]=p;
      else
      {
        p_Shift(&p,i*a->rank,R);
        res->m[0]=p_Add_q(res->m[0],p,R);
      }
    }
  }
  return res;
}

sm_Mult()

ideal sm_Mult	(	ideal	a,
		ideal	b,
		const ring	R
	)

Definition at line 1884 of file matpol.cc.

{
  int i, j, k;
  int m = a->rank;
  int p = IDELEMS(a);
  int q = IDELEMS(b);
 
  assume (IDELEMS(a)==b->rank);
  ideal c = idInit(q,m);
 
  for (i=0; i<m; i++)
  {
    for (k=0; k<p; k++)
    {
      poly aik;
      if ((aik=SMATELEM(a,i,k,R))!=NULL)
      {
        for (j=0; j<q; j++)
        {
          poly bkj=SMATELEM(b,k,j,R);
          if (bkj!=NULL)
          {
            poly s = p_Mult_q(p_Copy(aik,R) /*SMATELEM(a,i,k)*/, bkj/*SMATELEM(b,k,j)*/, R);
            if (s!=NULL) p_SetComp(s,i+1,R);
            c->m[j]=p_Add_q(c->m[j],s, R);
          }
        }
        p_Delete(&aik,R);
      }
    }
  }
  for(i=q-1;i>=0;i--) p_Normalize(c->m[i], R);
  return c;
}

sm_MultAndShift()

static ideal sm_MultAndShift ( poly  f, ideal  B, int  s, const ring  r  )

static

Definition at line 1796 of file matpol.cc.

{
  assume(f!=NULL);
  ideal res=idInit(IDELEMS(B),B->rank+s);
  int q=IDELEMS(B); // p x q
  for(int j=0;j<q;j++)
  {
    poly h=pp_Mult_qq(f,B->m[j],r);
    if (h!=NULL)
    {
      if (s>0) p_Shift(&h,s,r);
      res->m[j]=h;
    }
  }
  p_Delete(&f,r);
  return res;
}

sm_Sub()

ideal sm_Sub	(	ideal	a,
		ideal	b,
		const ring	R
	)

Definition at line 1874 of file matpol.cc.

{
  assume(IDELEMS(a)==IDELEMS(b));
  assume(a->rank==b->rank);
  ideal c=idInit(IDELEMS(a),a->rank);
  for (int k=IDELEMS(a)-1; k>=0; k--)
    c->m[k] = p_Sub(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
  return c;
}

sm_Tensor()

ideal sm_Tensor	(	ideal	A,
		ideal	B,
		const ring	r
	)

Definition at line 1824 of file matpol.cc.

{
  // size of the result m*p x n*q
  int n=IDELEMS(A); // m x n
  int m=A->rank;
  int q=IDELEMS(B); // p x q
  int p=B->rank;
  ideal res=idInit(n*q,m*p);
  poly *a=(poly*)omAlloc(m*sizeof(poly));
  for(int i=0; i<n; i++)
  {
    memset(a,0,m*sizeof(poly));
    p_Vec2Array(A->m[i],a,m,r);
    for(int j=0;j<m;j++)
    {
      if (a[j]!=NULL)
      {
        ideal sm=sm_MultAndShift(a[j], // A_i_j
                                 B,
                                 j*p, // shift j*p down
                                 r);
        sm_AddSubMat(res,sm,i*q,r); // add this columns to col i*q ff
        id_Delete(&sm,r); // delete the now empty ideal
      }
    }
  }
  omFreeSize(a,m*sizeof(poly));
  return res;
}

sm_Trace()

poly sm_Trace	(	ideal	a,
		const ring	R
	)

Definition at line 1966 of file matpol.cc.

{
  int i;
  int n = (IDELEMS(a)<a->rank) ? IDELEMS(a) : a->rank;
  poly  t = NULL;
 
  for (i=0; i<=n; i++)
    t = p_Add_q(t, p_Copy(SMATELEM(a,i,i,R), R), R);
  return t;
}

sm_UnFlatten()

ideal sm_UnFlatten	(	ideal	a,
		int	col,
		const ring	R
	)

Definition at line 1939 of file matpol.cc.

{
  if ((IDELEMS(a)!=1)
  ||((a->rank % col)!=0))
  {
    Werror("wrong format: %d x %d for unflatten",(int)a->rank,IDELEMS(a));
    return NULL;
  }
  int row=a->rank/col;
  ideal res=idInit(col,row);
  poly p=a->m[0];
  while(p!=NULL)
  {
    poly h=p_Head(p,R);
    int comp=p_GetComp(h,R);
    int c=(comp-1)/row;
    int r=comp%row; if (r==0) r=row;
    p_SetComp(h,r,R); p_SetmComp(h,R);
    res->m[c]=p_Add_q(res->m[c],h,R);
    pIter(p);
  }
  return res;
}

TraceOfProd()

poly TraceOfProd	(	matrix	a,
		matrix	b,
		int	n,
		const ring	R
	)

Definition at line 282 of file matpol.cc.

{
  int i, j;
  poly  p, t = NULL;
 
  for (i=1; i<=n; i++)
  {
    for (j=1; j<=n; j++)
    {
      p = pp_Mult_qq(MATELEM(a,i,j), MATELEM(b,j,i), R);
      t = p_Add_q(t, p, R);
    }
  }
  return t;
}