Changeset 5d9aa6 in git

Timestamp:

Apr 13, 2011, 8:49:42 PM (13 years ago)

Author:

Oleksandr Motsak <motsak@…>

Branches:

(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')

Children:

503a31399a9b66dd3306a8251a9a5042a4509f5c

Parents:

0a3a629f65533731dfccd53a13103d906918010b

git-author:

Oleksandr Motsak <motsak@mathematik.uni-kl.de>2011-04-13 20:49:42+02:00

git-committer:

Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 12:12:33+01:00

Message:

FIX: further fixes in matpol
ADD: sparse matrices (sparsmat.*) - needed for matpol (Berreiss?)

Location:

libpolys/polys

Files:

• Makefile.am (modified) (3 diffs)
• matpol.cc (modified) (90 diffs)
• simpleideals.h (modified) (1 diff)
• sparsmat.cc (moved) (moved from kernel/sparsmat.cc)
• sparsmat.h (moved) (moved from kernel/sparsmat.h)

libpolys/polys/Makefile.am

@@ -36 +36 @@
         pDebug.cc pInline0.cc polys0.cc prCopy.cc \
         kbuckets.cc sbuckets.cc ${USE_P_PROCS_STATIC_CC} ${USE_P_PROCS_DYNAMIC_CC} weight.cc \
-  simpleideals.cc matpol.cc
+  simpleideals.cc matpol.cc sparsmat.cc
 
 BUILT_SOURCES = templates/p_Procs.inc
@@ -44 +44 @@
 include_HEADERS = \
         monomials/ring.h nc/nc.h \
+  nc/sca.h nc/summator.h nc/ncSAFormula.h nc/ncSACache.h nc/ncSAMult.h \
         pInline0.h operations/pShallowCopyDelete.h \
         templates/p_MemAdd.h templates/p_MemCmp.h templates/p_MemCopy.h operations/p_Mult_q.h \
@@ -49 +50 @@
         templates/p_Procs_Dynamic.h templates/p_Procs_Impl.h templates/p_Procs_Set.h templates/p_Procs_Static.h \
         monomials/p_polys.h monomials/polys-impl.h monomials/maps.h polys.h prCopy.h prCopyMacros.h \
-        kbuckets.h sbuckets.h simpleideals.h weight.h
+        kbuckets.h sbuckets.h simpleideals.h weight.h \
+  matpol.h sparsmat.h
 
 P_PROCS_CPPFLAGS_COMMON = -DHAVE_CONFIG_H -DDYNAMIC_VERSION 

libpolys/polys/matpol.cc

@@ -37 +37 @@
 #include "prCopy.h"
 
+#include "sparsmat.h"
+   
 //omBin ip_smatrix_bin = omGetSpecBin(sizeof(ip_smatrix));
 #define ip_smatrix_bin sip_sideal_bin
@@ -45 +47 @@
 static void mpReplace(int j, int n, int &sign, int *perm);
 static int mpNextperm(perm * z, int max);
-static poly mpLeibnitz(matrix a);
-static poly minuscopy (poly p);
-static poly pInsert(poly p1, poly p2);
-static poly mpExdiv ( poly m, poly d, poly vars);
-static poly mpSelect (poly fro, poly what);
-
-static void mpPartClean(matrix, int, int);
-static void mpFinalClean(matrix);
-static int mpPrepareRow (matrix, int, int);
-static int mpPreparePiv (matrix, int, int);
-static int mpPivBar(matrix, int, int);
-static int mpPivRow(matrix, int, int);
-static float mpPolyWeight(poly);
-static void mpSwapRow(matrix, int, int, int);
-static void mpSwapCol(matrix, int, int, int);
-static void mpElimBar(matrix, matrix, poly, int, int);
-
-/*2
-* create a r x c zero-matrix
-*/
+static poly mp_Leibnitz(matrix a, const ring);
+static poly minuscopy (poly p, const ring);
+static poly p_Insert(poly p1, poly p2, const ring);
+static poly mp_Exdiv ( poly m, poly d, poly vars, const ring);
+static poly mp_Select (poly fro, poly what, const ring);
+
+static void mp_PartClean(matrix, int, int, const ring);
+static void mp_FinalClean(matrix, const ring);
+static int mp_PrepareRow (matrix, int, int, const ring);
+static int mp_PreparePiv (matrix, int, int, const ring);
+static int mp_PivBar(matrix, int, int, const ring);
+static int mp_PivRow(matrix, int, int, const ring);
+static float mp_PolyWeight(poly, const ring);
+static void mp_SwapRow(matrix, int, int, int, const ring);
+static void mp_SwapCol(matrix, int, int, int, const ring);
+static void mp_ElimBar(matrix, matrix, poly, int, int, const ring);
+
+/// create a r x c zero-matrix
 matrix mpNew(int r, int c)
 {
@@ -80 +80 @@
   {
     int s=r*c*sizeof(poly);
-    rc->m = (polyset)omAlloc0(s);
+    rc->m = (poly*)omAlloc0(s);
     //if (rc->m==NULL)
     //{
@@ -90 +90 @@
 }
 
-/// copies matrix a to b
-matrix mpCopy (matrix a, const ring r)
-{
-  idTest((ideal)a);
+/// copies matrix a (from ring r to r)
+matrix mp_Copy (matrix a, const ring r)
+{
+  id_Test((ideal)a, r);
   poly t;
   int i, m=MATROWS(a), n=MATCOLS(a);
@@ -103 +103 @@
     if (t!=NULL)
     {
-      pNormalize(t);
-      b->m[i] = pCopy(t);
+      p_Normalize(t, r);
+      b->m[i] = p_Copy(t, r);
     }
   }
@@ -111 +111 @@
 }
 
-/*2
-*copies matrix a from rSrc into rDst
-*/
-matrix mpCopy(const matrix a, const ring rSrc, const ring rDst)
-{
-  const ring save = currRing;
-
-#ifndef NDEBUG
-  if( currRing != rSrc )
-    rChangeCurrRing(rSrc);
-  idTest((ideal)a);
-#endif
+/// copies matrix a from rSrc into rDst
+matrix mp_Copy(const matrix a, const ring rSrc, const ring rDst)
+{
+  id_Test((ideal)a, rSrc);
 
   poly t;
@@ -140 +132 @@
   b->rank=a->rank;
 
-#ifndef NDEBUG
-  if( currRing != rDst )
-    rChangeCurrRing(rDst);
-  idTest((ideal)b);
-#endif
-
-  if( save != currRing )
-    rChangeCurrRing(save);
+  id_Test((ideal)b, rDst);
 
   return b;
@@ -154 +139 @@
 
 
-/*2
-* make it a p * unit matrix
-*/
-matrix mpInitP(int r, int c, poly p)
+/// make it a p * unit matrix
+matrix mp_InitP(int r, int c, poly p, const ring _R)
 {
   matrix rc = mpNew(r,c);
   int i=si_min(r,c), n = c*(i-1)+i-1, inc = c+1;
 
-  pNormalize(p);
+  p_Normalize(p, _R);
   while (n>0)
   {
-    rc->m[n] = pCopy(p);
+    rc->m[n] = p_Copy(p, _R);
     n -= inc;
   }
@@ -172 +155 @@
 }
 
-/*2
-* make it a v * unit matrix
-*/
-matrix mpInitI(int r, int c, int v)
-{
-  return mpInitP(r,c,pISet(v));
-}
-
-/*2
-* c = f*a
-*/
-matrix mpMultI(matrix a, int f)
+/// make it a v * unit matrix
+matrix mp_InitI(int r, int c, int v, const ring R)
+{
+  return mp_InitP(r, c, p_ISet(v, R), R);
+}
+
+/// c = f*a
+matrix mp_MultI(matrix a, int f, const ring R)
 {
   int k, n = a->nrows, m = a->ncols;
-  poly p = pISet(f);
+  poly p = p_ISet(f, R);
   matrix c = mpNew(n,m);
 
   for (k=m*n-1; k>0; k--)
-    c->m[k] = ppMult_qq(a->m[k], p);
-  c->m[0] = pMult(pCopy(a->m[0]), p);
+    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;
 }
 
-/*2
-* multiply a matrix 'a' by a poly 'p', destroy the args
-*/
-matrix mpMultP(matrix a, poly p)
+/// multiply a matrix 'a' by a poly 'p', destroy the args
+matrix mp_MultP(matrix a, poly p, const ring R)
 {
   int k, n = a->nrows, m = a->ncols;
 
-  pNormalize(p);
+  p_Normalize(p, R);
   for (k=m*n-1; k>0; k--)
   {
     if (a->m[k]!=NULL)
-      a->m[k] = pMult(a->m[k], pCopy(p));
-  }
-  a->m[0] = pMult(a->m[0], p);
+      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;
 }
@@ -215 +192 @@
 * multiply a poly 'p' by a matrix 'a', destroy the args
 */
-matrix pMultMp(poly p, matrix a)
+matrix pMultMp(poly p, matrix a, const ring R)
 {
   int k, n = a->nrows, m = a->ncols;
 
-  pNormalize(p);
+  p_Normalize(p, R);
   for (k=m*n-1; k>0; k--)
   {
     if (a->m[k]!=NULL)
-      a->m[k] = pMult(pCopy(p), a->m[k]);
-  }
-  a->m[0] = pMult(p, a->m[0]);
+      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;
 }
 
-matrix mpAdd(matrix a, matrix b)
+matrix mp_Add(matrix a, matrix b, const ring R)
 {
   int k, n = a->nrows, m = a->ncols;
@@ -242 +219 @@
   matrix c = mpNew(n,m);
   for (k=m*n-1; k>=0; k--)
-    c->m[k] = pAdd(pCopy(a->m[k]), pCopy(b->m[k]));
+    c->m[k] = p_Add_q(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
   return c;
 }
 
-matrix mpSub(matrix a, matrix b)
+matrix mp_Sub(matrix a, matrix b, const ring R)
 {
   int k, n = a->nrows, m = a->ncols;
@@ -259 +236 @@
   matrix c = mpNew(n,m);
   for (k=m*n-1; k>=0; k--)
-    c->m[k] = pSub(pCopy(a->m[k]), pCopy(b->m[k]));
+    c->m[k] = p_Sub(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
   return c;
 }
 
-matrix mpMult(matrix a, matrix b)
+matrix mp_Mult(matrix a, matrix b, const ring R)
 {
   int i, j, k;
@@ -293 +270 @@
           {
             poly *cij=&(MATELEM(c,i,j));
-            poly s = ppMult_qq(aik /*MATELEM(a,i,k)*/, bkj/*MATELEM(b,k,j)*/);
+            poly s = pp_Mult_qq(aik /*MATELEM(a,i,k)*/, bkj/*MATELEM(b,k,j)*/, R);
             if (/*MATELEM(c,i,j)*/ (*cij)==NULL) (*cij)=s;
-            else (*cij) = pAdd((*cij) /*MATELEM(c,i,j)*/ ,s);
+            else (*cij) = p_Add_q((*cij) /*MATELEM(c,i,j)*/ ,s, R);
           }
         }
@@ -303 +280 @@
     }
   }
-  for(i=m*q-1;i>=0;i--) pNormalize(c->m[i]);
+  for(i=m*q-1;i>=0;i--) p_Normalize(c->m[i], R);
   return c;
 }
 
-matrix mpTransp(matrix a)
+matrix mp_Transp(matrix a, const ring R)
 {
   int    i, j, r = MATROWS(a), c = MATCOLS(a);
@@ -318 +295 @@
     for (j=0; j<r; j++)
     {
-      if (a->m[j*c+i]!=NULL) *p = pCopy(a->m[j*c+i]);
+      if (a->m[j*c+i]!=NULL) *p = p_Copy(a->m[j*c+i], R);
       p++;
     }
@@ -328 +305 @@
 *returns the trace of matrix a
 */
-poly mpTrace ( matrix a)
+poly mp_Trace ( matrix a, const ring R)
 {
   int i;
@@ -335 +312 @@
 
   for (i=1; i<=n; i++)
-    t = pAdd(t, pCopy(MATELEM(a,i,i)));
+    t = p_Add_q(t, p_Copy(MATELEM(a,i,i), R), R);
   return t;
 }
@@ -342 +319 @@
 *returns the trace of the product of a and b
 */
-poly TraceOfProd ( matrix a, matrix b, int n)
+poly TraceOfProd ( matrix a, matrix b, int n, const ring _R)
 {
   int i, j;
@@ -351 +328 @@
     for (j=1; j<=n; j++)
     {
-      p = ppMult_qq(MATELEM(a,i,j), MATELEM(b,j,i));
-      t = pAdd(t, p);
+      p = pp_Mult_qq(MATELEM(a,i,j), MATELEM(b,j,i), _R);
+      t = p_Add_q(t, p, _R);
     }
   }
@@ -388 +365 @@
   int       *qrow, *qcol;
   poly      *Xarray;
+  ring      _R;
+  
   void mpInitMat();
   poly * mpRowAdr(int);
@@ -396 +375 @@
   void mpColSwap(int, int);
   public:
-  mp_permmatrix() : a_m(0) {}
-  mp_permmatrix(matrix);
+  mp_permmatrix() : a_m(0), _R(NULL) {}
+  mp_permmatrix(matrix, const ring);
   mp_permmatrix(mp_permmatrix *);
   ~mp_permmatrix();
@@ -421 +400 @@
 #define SIZE_OF_SYSTEM_PAGE 4096
 #endif
-/*2
-* entries of a are minors and go to result (only if not in R)
-*/
-void mpMinorToResult(ideal result, int &elems, matrix a, int r, int c,
-                     ideal R)
+
+/*
+/// entries of a are minors and go to result (only if not in R)
+void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
+                     ideal R, const ring _R)
 {
   poly *q1;
@@ -471 +450 @@
 }
 
-/*2
-* produces recursively the ideal of all arxar-minors of a
-*/
-void mpRecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
-              poly barDiv, ideal R)
+/// produces recursively the ideal of all arxar-minors of a
+void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
+              poly barDiv, ideal R, const ring _R)
 {
   int k;
@@ -483 +460 @@
   loop
   {
-/*--- look for an optimal row and bring it to last position ------------*/
+// --- look for an optimal row and bring it to last position ------------
     if(mpPrepareRow(a,lr,lc)==0) break;
-/*--- now take all pivots from the last row ------------*/
+// --- now take all pivots from the last row ------------
     k = lc;
     loop
@@ -501 +478 @@
     }
     if (ar>=kr) break;
-/*--- now we have to take out the last row...------------*/
+// --- now we have to take out the last row...------------
     lr = kr;
     kr--;
@@ -507 +484 @@
   mpFinalClean(nextLevel);
 }
+*/
+
 
 /*2
@@ -512 +491 @@
 *uses Bareiss algorithm
 */
-poly mpDetBareiss (matrix a)
+poly mp_DetBareiss (matrix a, const ring R)
 {
   int s;
@@ -521 +500 @@
     return NULL;
   }
-  matrix c = mpCopy(a);
-  mp_permmatrix *Bareiss = new mp_permmatrix(c);
+  matrix c = mp_Copy(a, R);
+  mp_permmatrix *Bareiss = new mp_permmatrix(c, R);
   row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
 
@@ -541 +520 @@
   res = MATELEM(c,1,1);
   MATELEM(c,1,1) = NULL;
-  idDelete((ideal *)&c);
+  id_Delete((ideal *)&c, R);
   if (s < 0)
-    res = pNeg(res);
+    res = p_Neg(res, R);
   return res;
 }
@@ -551 +530 @@
 *uses Newtons formulea for symmetric functions
 */
-poly mpDet (matrix m)
+poly mp_Det (matrix m, const ring R)
 {
   int i,j,k,n;
@@ -565 +544 @@
     return NULL;
   }
-  k=rChar();
+  k=rChar(R);
   if ((k > 0) && (k <= n))
-    return mpLeibnitz(m);
-  ONE = nInit(1);
-  ma[1]=mpCopy(m);
+    return mp_Leibnitz(m, R);
+  ONE = n_Init(1, R);
+  ma[1]=mp_Copy(m, R);
   k = (n+1) / 2;
   s = mpNew(1, n);
-  MATELEM(s,1,1) = mpTrace(m);
+  MATELEM(s,1,1) = mp_Trace(m, R);
   for (i=2; i<=k; i++)
   {
     //ma[i] = mpNew(n,n);
-    ma[i]=mpMult(ma[i-1], ma[1]);
-    MATELEM(s,1,i) = mpTrace(ma[i]);
-    pTest(MATELEM(s,1,i));
+    ma[i]=mp_Mult(ma[i-1], ma[1], R);
+    MATELEM(s,1,i) = mp_Trace(ma[i], R);
+    p_Test(MATELEM(s,1,i), R);
   }
   for (i=k+1; i<=n; i++)
   {
-    MATELEM(s,1,i) = TraceOfProd(ma[i / 2], ma[(i+1) / 2], n);
-    pTest(MATELEM(s,1,i));
+    MATELEM(s,1,i) = TraceOfProd(ma[i / 2], ma[(i+1) / 2], n, R);
+    p_Test(MATELEM(s,1,i), R);
   }
   for (i=1; i<=k; i++)
-    idDelete((ideal *)&(ma[i]));
+    id_Delete((ideal *)&(ma[i]), R);
 /* the array s contains the traces of the powers of the matrix m,
 *  these are the power sums of the eigenvalues of m */
   a = mpNew(1,n);
-  MATELEM(a,1,1) = minuscopy(MATELEM(s,1,1));
+  MATELEM(a,1,1) = minuscopy(MATELEM(s,1,1), R);
   for (i=2; i<=n; i++)
   {
-    p = pCopy(MATELEM(s,1,i));
+    p = p_Copy(MATELEM(s,1,i), R);
     for (j=i-1; j>=1; j--)
     {
-      q = ppMult_qq(MATELEM(s,1,j), MATELEM(a,1,i-j));
-      pTest(q);
-      p = pAdd(p,q);
+      q = pp_Mult_qq(MATELEM(s,1,j), MATELEM(a,1,i-j), R);
+      p_Test(q, R);
+      p = p_Add_q(p,q, R);
     }
     // c= -1/i
-    d = nInit(-(int)i);
-    c = nDiv(ONE, d);
-    nDelete(&d);
-
-    pMult_nn(p, c);
-    pTest(p);
+    d = n_Init(-(int)i, R);
+    c = n_Div(ONE, d, R);
+    n_Delete(&d, R);
+
+    p_Mult_nn(p, c, R);
+    p_Test(p, R);
     MATELEM(a,1,i) = p;
-    nDelete(&c);
+    n_Delete(&c, R);
   }
 /* the array a contains the elementary symmetric functions of the
@@ -614 +593 @@
   for (i=1; i<=n-1; i++)
   {
-    //pDelete(&(MATELEM(a,1,i)));
-    pDelete(&(MATELEM(s,1,i)));
-  }
-  pDelete(&(MATELEM(s,1,n)));
+    //p_Delete(&(MATELEM(a,1,i)), R);
+    p_Delete(&(MATELEM(s,1,i)), R);
+  }
+  p_Delete(&(MATELEM(s,1,n)), R);
 /* up to a sign, the determinant is the n-th elementary symmetric function */
   if ((n/2)*2 < n)
   {
-    d = nInit(-1);
-    pMult_nn(MATELEM(a,1,n), d);
-    nDelete(&d);
-  }
-  nDelete(&ONE);
-  idDelete((ideal *)&s);
+    d = n_Init(-1, R);
+    p_Mult_nn(MATELEM(a,1,n), d, R);
+    n_Delete(&d, R);
+  }
+  n_Delete(&ONE, R);
+  id_Delete((ideal *)&s, R);
   poly result=MATELEM(a,1,n);
   MATELEM(a,1,n)=NULL;
-  idDelete((ideal *)&a);
+  id_Delete((ideal *)&a, R);
   return result;
 }
@@ -636 +615 @@
 * compute all ar-minors of the matrix a
 */
-matrix mpWedge(matrix a, int ar)
+matrix mp_Wedge(matrix a, int ar, const ring R)
 {
   int     i,j,k,l;
@@ -699 +678 @@
 //    {
 //      p=pHomogen(MATELEM(a,i,j),v);
-//      pDelete(&(MATELEM(a,i,j)));
+//      p_Delete(&(MATELEM(a,i,j)), ?);
 //      MATELEM(a,i,j)=p;
 //    }
@@ -710 +689 @@
 * var has to be the number of a variable
 */
-matrix mpCoeffs (ideal I, int var)
+matrix mp_Coeffs (ideal I, int var, const ring R)
 {
   poly h,f;
@@ -721 +700 @@
     while (f!=NULL)
     {
-      l=pGetExp(f,var);
+      l=p_GetExp(f,var, R);
       if (l>m) m=l;
       pIter(f);
@@ -735 +714 @@
     while (f!=NULL)
     {
-      l=pGetExp(f,var);
-      pSetExp(f,var,0);
-      c=si_max((int)pGetComp(f),1);
-      pSetComp(f,0);
-      pSetm(f);
+      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)
-      //  =pAdd(MATELEM(co,c*(m+1)-l,i+1),f);
+      //  =p_Add_q(MATELEM(co,c*(m+1)-l,i+1),f, R);
       MATELEM(co,(c-1)*(m+1)+l+1,i+1)
-        =pAdd(MATELEM(co,(c-1)*(m+1)+l+1,i+1),f);
+        =p_Add_q(MATELEM(co,(c-1)*(m+1)+l+1,i+1),f, R);
       /* iterate f*/
       f=h;
     }
   }
-  idDelete(&I);
+  id_Delete(&I, R);
   return co;
 }
@@ -760 +739 @@
 * build the matrix of the corresponding monomials in m
 */
-void   mpMonomials(matrix c, int r, int var, matrix m)
+void   mp_Monomials(matrix c, int r, int var, matrix m, const ring R)
 {
   /* clear contents of m*/
@@ -768 +747 @@
     for(l=MATCOLS(m);l>0;l--)
     {
-      pDelete(&MATELEM(m,k,l));
+      p_Delete(&MATELEM(m,k,l), R);
     }
   }
@@ -780 +759 @@
   int p=MATCOLS(m)/r-1;
   /* fill in the powers of x_var=h*/
-  poly h=pOne();
+  poly h=p_One(R);
   for(k=r;k>0; k--)
   {
-    MATELEM(m,k,k*(p+1))=pOne();
+    MATELEM(m,k,k*(p+1))=p_One(R);
   }
   for(l=p;l>=0; l--)
@@ -791 +770 @@
     for(k=r;k>0; k--)
     {
-      MATELEM(m,k,k*(p+1)-l)=pCopy(h);
-    }
-  }
-  pDelete(&h);
-}
-
-matrix mpCoeffProc (poly f, poly vars)
+      MATELEM(m,k,k*(p+1)-l)=p_Copy(h, R);
+    }
+  }
+  p_Delete(&h), R;
+}
+
+matrix mp_CoeffProc (poly f, poly vars, const ring R)
 {
   assume(vars!=NULL);
@@ -808 +787 @@
   {
     co = mpNew(2, 1);
-    MATELEM(co,1,1) = pOne();
+    MATELEM(co,1,1) = p_One(R);
     MATELEM(co,2,1) = NULL;
     return co;
@@ -849 +828 @@
         if (h!=NULL)
         {
-          MATELEM(co,2,i) = pAdd(MATELEM(co,2,i), h);
+          MATELEM(co,2,i) = p_Add_q(MATELEM(co,2,i), h, R);
           break;
         }
@@ -861 +840 @@
           if (h!=NULL)
           {
-            MATELEM(co,2,pos_of_1) = pAdd(MATELEM(co,2,pos_of_1), h);
+            MATELEM(co,2,pos_of_1) = p_Add_q(MATELEM(co,2,pos_of_1), h, R);
           }
         }
@@ -878 +857 @@
 * consider all variables in vars
 */
-static poly mpExdiv ( poly m, poly d, poly vars)
+static poly mp_Exdiv ( poly m, poly d, poly vars, const ring _R)
 {
   int i;
@@ -888 +867 @@
       if (pGetExp(d,i) != pGetExp(h,i))
       {
-        pDelete(&h);
+        p_Delete(&h, _R);
         return NULL;
       }
@@ -898 +877 @@
 }
 
-void mpCoef2(poly v, poly mon, matrix *c, matrix *m)
+void mp_Coef2(poly v, poly mon, matrix *c, matrix *m, const ring _R)
 {
   polyset s;
@@ -947 +926 @@
         {
           pSetComp(h,0);
-          MATELEM(*c,j,i) = pAdd(MATELEM(*c,j,i), h);
+          MATELEM(*c,j,i) = p_Add_q(MATELEM(*c,j,i), h, _R);
           break;
         }
@@ -961 +940 @@
 
 
-BOOLEAN mpEqual(matrix a, matrix b)
+BOOLEAN mp_Equal(matrix a, matrix b, const ring _R)
 {
   if ((MATCOLS(a)!=MATCOLS(b)) || (MATROWS(a)!=MATROWS(b)))
@@ -974 +953 @@
     else
       if (b->m[i]==NULL) return FALSE;
-      else if (pCmp(a->m[i],b->m[i])!=0) return FALSE;
+      else if (p_Cmp(a->m[i],b->m[i], _R)!=0) return FALSE;
     i--;
   }
@@ -981 +960 @@
   {
 #if 0
-    poly tt=pSub(pCopy(a->m[i]),pCopy(b->m[i]));
+    poly tt=p_Sub(p_Copy(a->m[i], _R),p_Copy(b->m[i], _R), _R);
     if (tt!=NULL)
     {
-      pDelete(&tt);
+      p_Delete(&tt, _R);
       return FALSE;
     }
 #else
-    if(!pEqualPolys(a->m[i],b->m[i])) return FALSE;
+    if(!p_EqualPolys(a->m[i],b->m[i], _R)) return FALSE;
 #endif
     i--;
@@ -1014 +993 @@
 }
 
-mp_permmatrix::mp_permmatrix(matrix A) : sign(1)
+mp_permmatrix::mp_permmatrix(matrix A, const ring r) : sign(1),  _R(r)
 {
   a_m = A->nrows;
@@ -1030 +1009 @@
   a_n = M->s_n;
   sign = M->sign;
+  _R = M->_R;
+  
   this->mpInitMat();
   Xarray = (poly *)omAlloc0(a_m*a_n*sizeof(poly));
@@ -1058 +1039 @@
     {
       for (k=a_m*a_n-1; k>=0; k--)
-        pDelete(&Xarray[k]);
+        p_Delete(&Xarray[k], _R);
       omFreeSize((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
     }
@@ -1086 +1067 @@
 void mp_permmatrix::mpDelElem(int r, int c)
 {
-  pDelete(&Xarray[a_n*qrow[r]+qcol[c]]);
+  p_Delete(&Xarray[a_n*qrow[r]+qcol[c]], _R);
 }
 
@@ -1116 +1097 @@
           {
             q1 = SM_MULT(a[jj], piv, div);
-            pDelete(&a[jj]);
-            q2 = pAdd(q2, q1);
+            p_Delete(&a[jj]), _R;
+            q2 = p_Add_q(q2, q1, _R);
           }
         }
@@ -1128 +1109 @@
         a[jj] = q2;
       }
-      pDelete(&a[qcol[s_n]]);
+      p_Delete(&a[qcol[s_n]], _R);
     }
     else
@@ -1138 +1119 @@
         {
           q2 = SM_MULT(a[jj], piv, div);
-          pDelete(&a[jj]);
+          p_Delete(&a[jj], _R);
           if (div)
             SM_DIV(q2, div);
@@ -1178 +1159 @@
           fo = f1;
           if (iopt >= 0)
-            pDelete(&(this->mpRowAdr(iopt)[qcol[0]]));
+            p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), _R);
           iopt = i;
         }
         else
-          pDelete(&(this->mpRowAdr(i)[qcol[0]]));
+          p_Delete(&(this->mpRowAdr(i)[qcol[0]]), _R);
       }
     }
@@ -1257 +1238 @@
         fo = f1;
         if (iopt >= 0)
-        pDelete(&(this->mpRowAdr(iopt)[qcol[0]]));
+        p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), _R);
         iopt = row;
       }
       else
-        pDelete(&(this->mpRowAdr(row)[qcol[0]]));
+        p_Delete(&(this->mpRowAdr(row)[qcol[0]]), _R);
     }
     if (iopt >= 0)
@@ -1458 +1439 @@
 }
 
-/*
-* perform replacement for pivot strategy in Bareiss algorithm
-* change sign of determinant
-*/
+/// perform replacement for pivot strategy in Bareiss algorithm
+/// change sign of determinant
 static void mpReplace(int j, int n, int &sign, int *perm)
 {
@@ -1527 +1506 @@
 }
 
-static poly mpLeibnitz(matrix a)
+static poly mp_Leibnitz(matrix a, const ring _R)
 {
   int i, e, n;
@@ -1535 +1514 @@
   n = MATROWS(a);
   memset(&z,0,(n+2)*sizeof(int));
-  p = pOne();
+  p = p_One(_R);
   for (i=1; i <= n; i++)
-    p = pMult(p, pCopy(MATELEM(a, i, i)));
+    p = p_Mult_q(p, p_Copy(MATELEM(a, i, i), _R), _R);
   d = p;
   for (i=1; i<= n; i++)
@@ -1548 +1527 @@
     {
       e = mpNextperm((perm *)&z, n);
-      p = pOne();
+      p = p_One(_R);
       for (i = 1; i <= n; i++)
-        p = pMult(p, pCopy(MATELEM(a, i, z[i])));
+        p = p_Mult_q(p, p_Copy(MATELEM(a, i, z[i]), _R), _R);
       if (z[0] > 0)
-        d = pAdd(d, p);
+        d = p_Add_q(d, p, _R);
       else
-        d = pSub(d, p);
+        d = p_Sub_q(d, p, _R);
     }
   }
@@ -1560 +1539 @@
 }
 
-static poly minuscopy (poly p)
+static poly minuscopy (poly p, const ring R)
 {
   poly w;
   number  e;
-  e = nInit(-1);
-  w = pCopy(p);
-  pMult_nn(w, e);
-  nDelete(&e);
+  e = n_Init(-1, R);
+  w = p_Copy(p, R);
+  p_Mult_nn(w, e, R);
+  n_Delete(&e, R);
   return w;
 }
@@ -1575 +1554 @@
 * arguments are destroyed
 */
-static poly pInsert(poly p1, poly p2)
+static poly p_Insert(poly p1, poly p2, const ring R)
 {
   poly a1, p, a2, a;
@@ -1584 +1563 @@
   a1 = p1;
   a2 = p2;
-  a = p  = pOne();
+  a = p  = p_One(R);
   loop
   {
@@ -1633 +1612 @@
 *    x^i*y^j (i, j >= 0) that appear in fro
 */
-static poly mpSelect (poly fro, poly what)
+static poly mp_Select (poly fro, poly what, const ring _R)
 {
   int i;
@@ -1640 +1619 @@
   while (fro!=NULL)
   {
-    h = pOne();
+    h = p_One(_R);
     for (i=1; i<=pVariables; i++)
-      pSetExp(h,i, pGetExp(fro,i) * pGetExp(what, i));
-    pSetComp(h, pGetComp(fro));
-    pSetm(h);
-    res = pInsert(h, res);
+      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;
   }
@@ -1667 +1646 @@
 */
 
-static void mpPartClean(matrix a, int lr, int lc)
+static void mp_PartClean(matrix a, int lr, int lc, const ring _R)
 {
   poly *q1;
@@ -1675 +1654 @@
   {
     q1 = &(a->m)[i*a->ncols];
-    for (j=lc-1;j>=0;j--) if(q1[j]) pDelete(&q1[j]);
-  }
-}
-
-static void mpFinalClean(matrix a)
+    for (j=lc-1;j>=0;j--) if(q1[j]) p_Delete(&q1[j], _R);
+  }
+}
+
+static void mp_FinalClean(matrix a, const ring _R)
 {
   omFreeSize((ADDRESS)a->m,a->nrows*a->ncols*sizeof(poly));
@@ -1689 +1668 @@
 *  for application in minor
 */
-static int mpPrepareRow (matrix a, int lr, int lc)
+static int mp_PrepareRow (matrix a, int lr, int lc, const ring _R)
 {
   int r;
 
-  r = mpPivBar(a,lr,lc);
+  r = mp_PivBar(a,lr,lc, _R);
   if(r==0) return 0;
-  if(r<lr) mpSwapRow(a, r, lr, lc);
+  if(r<lr) mp_SwapRow(a, r, lr, lc, _R);
   return 1;
 }
@@ -1703 +1682 @@
 *  for application in minor
 */
-static int mpPreparePiv (matrix a, int lr, int lc)
+static int mp_PreparePiv (matrix a, int lr, int lc, const ring _R)
 {
   int c;
 
-  c = mpPivRow(a, lr, lc);
+  c = mp_PivRow(a, lr, lc, _R);
   if(c==0) return 0;
-  if(c<lc) mpSwapCol(a, c, lr, lc);
+  if(c<lc) mp_SwapCol(a, c, lr, lc, _R);
   return 1;
 }
@@ -1716 +1695 @@
 * find best row
 */
-static int mpPivBar(matrix a, int lr, int lc)
+static int mp_PivBar(matrix a, int lr, int lc, const ring _R)
 {
   float f1, f2;
@@ -1731 +1710 @@
     {
       if (q1[j]!=NULL)
-        f2 += mpPolyWeight(q1[j]);
+        f2 += mp_PolyWeight(q1[j], _R);
     }
     if ((f2!=0.0) && (f2<f1))
@@ -1746 +1725 @@
 * find pivot in the last row
 */
-static int mpPivRow(matrix a, int lr, int lc)
+static int mp_PivRow(matrix a, int lr, int lc, const ring _R)
 {
   float f1, f2;
@@ -1759 +1738 @@
     if (q1[j]!=NULL)
     {
-      f2 = mpPolyWeight(q1[j]);
+      f2 = mp_PolyWeight(q1[j], _R);
       if (f2<f1)
       {
@@ -1774 +1753 @@
 * weigth of a polynomial, for pivot strategy
 */
-static float mpPolyWeight(poly p)
+static float mp_PolyWeight(poly p, const ring _R)
 {
   int i;
@@ -1781 +1760 @@
   if (pNext(p) == NULL)
   {
-    res = (float)nSize(pGetCoeff(p));
+    res = (float)n_Size(p_GetCoeff(p, _R), _R);
     for (i=pVariables;i>0;i--)
     {
-      if(pGetExp(p,i)!=0)
+      if(p_GetExp(p,i, _R)!=0)
       {
         res += 2.0;
@@ -1796 +1775 @@
     do
     {
-      res += (float)nSize(pGetCoeff(p))+2.0;
+      res += (float)n_Size(p_GetCoeff(p, _R), _R) + 2.0;
       pIter(p);
     }
@@ -1808 +1787 @@
   poly sw;
   int j;
-  polyset a2 = a->m, a1 = &a2[a->ncols*(pos-1)];
+  poly* a2 = a->m, a1 = &a2[a->ncols*(pos-1)];
 
   a2 = &a2[a->ncols*(lr-1)];
@@ -1823 +1802 @@
   poly sw;
   int j;
-  polyset a2 = a->m, a1 = &a2[pos-1];
+  poly* a2 = a->m, a1 = &a2[pos-1];
 
   a2 = &a2[lc-1];
@@ -1834 +1813 @@
 }
 
-static void mpElimBar(matrix a0, matrix re, poly div, int lr, int lc)
+static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring _R)
 {
   int r=lr-1, c=lc-1;
@@ -1861 +1840 @@
           {
             q2 = SM_MULT(ap[j], elim, div);
-            q1 = pAdd(q1,q2);
+            q1 = p_Add_q(q1,q2, _R);
           }
         }
@@ -1890 +1869 @@
 }
 
-BOOLEAN mpIsDiagUnit(matrix U)
+BOOLEAN mp_IsDiagUnit(matrix U, const ring _R)
 {
   if(MATROWS(U)!=MATCOLS(U))
@@ -1900 +1879 @@
       if (i==j)
       {
-        if (!pIsUnit(MATELEM(U,i,i))) return FALSE;
+        if (!p_IsUnit(MATELEM(U,i,i), _R)) return FALSE;
       }
       else if (MATELEM(U,i,j)!=NULL) return FALSE;
@@ -1908 +1887 @@
 }
 
-void iiWriteMatrix(matrix im, const char *n, int dim,int spaces)
+void iiWriteMatrix(matrix im, const char *n, int dim, const ring r, int spaces)
 {
   int i,ii = MATROWS(im)-1;
@@ -1923 +1902 @@
       else if (dim == 1) Print("%s[%u]=",n,j+1);
       else if (dim == 0) Print("%s=",n);
-      if ((i<ii)||(j<jj)) pWrite(*pp++);
-      else                pWrite0(*pp);
-    }
-  }
-}
-
-char * iiStringMatrix(matrix im, int dim,char ch)
+      if ((i<ii)||(j<jj)) p_Write(*pp++, r);
+      else                p_Write0(*pp, r);
+    }
+  }
+}
+
+char * iiStringMatrix(matrix im, int dim, , const ring r, char ch)
 {
   int i,ii = MATROWS(im);
@@ -1940 +1919 @@
     for (j=0; j<jj; j++)
     {
-      pString0(*pp++);
+      p_String0(*pp++, r);
       s=StringAppend("%c",ch);
       if (dim > 1) s = StringAppendS("\n");
@@ -1949 +1928 @@
 }
 
-void   mpDelete(matrix* a, const ring r)
+void   mp_Delete(matrix* a, const ring r)
 {
   id_Delete((ideal *) a, r);

libpolys/polys/simpleideals.h

@@ -102 +102 @@
 intvec *id_Sort(ideal id,BOOLEAN nolex, const ring r);
 
+
+int     binom (int n,int r);
+
 #endif