Changeset a5ceaf in git

Timestamp:

Jun 2, 1999, 4:09:36 PM (24 years ago)

Author:

Wilfred Pohl <pohl@…>

Branches:

(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')

Children:

c801beb15dbf5f6a7997a33c6be813e2a7b471c1

Parents:

1e7da0633882bafa6663fdd72dc3ddf832a44096

Message:

det for sparse


git-svn-id: file:///usr/local/Singular/svn/trunk@3088 2c84dea3-7e68-4137-9b89-c4e89433aadc

Location:

Singular

Files:

• sparsmat.cc (modified) (13 diffs)
• sparsmat.h (modified) (2 diffs)

Singular/sparsmat.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: sparsmat.cc,v 1.5 1999-03-15 16:55:42 Singular Exp $ */
+/* $Id: sparsmat.cc,v 1.6 1999-06-02 14:09:34 pohl Exp $ */
 
 /*
@@ -115 +115 @@
   void smMultCol();
   poly smMultPoly(smpoly);
+  void smActDel();
+  void smColDel();
+  void smPivDel();
+  void smSign();
 public:
   sparse_mat(ideal);
   ~sparse_mat();
+  int smGetSign() { return sign; }
   smpoly * smGetAct() { return m_act; }
   int smGetRed() { return tored; }
   ideal smRes2Mod();
+  poly smDet();
   void smBareiss(int, int);
   void smNewBareiss(int, int);
@@ -126 +132 @@
 };
 
-/* ----------------- basics (used from 'C') ------------------ */
-
-lists smCallBareiss(ideal I, int x, int y)
-{
-  lists res=(lists)Alloc(sizeof(slists));
-  ideal II = idCopy(I);
-  sparse_mat *bareiss = new sparse_mat(II);
-  intvec *v;
-  ideal m;
-
-  idDelete(&II);
-  if (bareiss->smGetAct() == NULL)
-  {
-    Free((ADDRESS)res, sizeof(slists));
-    return NULL;
-  }
-  bareiss->smBareiss(x, y);
-  m = bareiss->smRes2Mod();
-  v = new intvec(bareiss->smGetRed());
-  bareiss->smToIntvec(v);
-  delete bareiss;
-  res->Init(2);
-  res->m[0].rtyp=MODUL_CMD;
-  res->m[0].data=(void *)m;
-  res->m[1].rtyp=INTVEC_CMD;
-  res->m[1].data=(void *)v;
-  return res;
-}
-
-lists smCallNewBareiss(ideal I, int x, int y)
-{
-  lists res=(lists)Alloc(sizeof(slists));
+/* ----------------- ops with rings ------------------ */
+ideal smRingCopy(ideal I, ring *ri, sip_sring &tmpR)
+{
   ring origR =NULL;
-  sip_sring tmpR;
   ideal II;
   if (currRing->order[0]!=ringorder_c)
@@ -186 +162 @@
     II=idCopy(I);
   }
+  *ri = origR;
+  return II;
+}
+
+void smRingClean(ring origR, ip_sring &tmpR)
+{
+  rChangeCurrRing(origR,TRUE);
+  rUnComplete(&tmpR);
+  Free((ADDRESS)tmpR.order,3*sizeof(int));
+  Free((ADDRESS)tmpR.block0,3*sizeof(int));
+  Free((ADDRESS)tmpR.block1,3*sizeof(int));
+}
+
+/* ----------------- basics (used from 'C') ------------------ */
+/*2
+*returns the determinant of the module I;
+*uses Bareiss algorithm
+*/
+poly smCallDet(ideal I)
+{
+  int r=idRankFreeModule(I);
+  if (I->ncols != r)
+  {
+    Werror("det of %d x %d module (matrix)",r,I->ncols);
+    return NULL;
+  }
+  poly res,save;
+  ring origR;
+  sip_sring tmpR;
+  ideal II=smRingCopy(I,&origR,tmpR);
+  sparse_mat *det = new sparse_mat(II);
+
+  idDelete(&II);
+  if (det->smGetAct() == NULL)
+  {
+    delete det;
+    if (origR!=NULL) smRingClean(origR,tmpR);
+    return NULL;
+  }
+  res=det->smDet();
+  if(det->smGetSign()<0) res=pNeg(res);
+  delete det;
+  if (origR!=NULL)
+  {
+    rChangeCurrRing(origR,TRUE);
+    save = res;
+    res = pFetchCopy(&tmpR, save);
+    rChangeCurrRing(&tmpR,FALSE);
+    pDelete(&save);
+    smRingClean(origR,tmpR);
+  }
+  return res;
+}
+
+lists smCallBareiss(ideal I, int x, int y)
+{
+  lists res=(lists)Alloc(sizeof(slists));
+  ideal II = idCopy(I);
+  sparse_mat *bareiss = new sparse_mat(II);
+  intvec *v;
+  ideal m;
+
+  idDelete(&II);
+  if (bareiss->smGetAct() == NULL)
+  {
+    Free((ADDRESS)res, sizeof(slists));
+    return NULL;
+  }
+  bareiss->smBareiss(x, y);
+  m = bareiss->smRes2Mod();
+  v = new intvec(bareiss->smGetRed());
+  bareiss->smToIntvec(v);
+  delete bareiss;
+  res->Init(2);
+  res->m[0].rtyp=MODUL_CMD;
+  res->m[0].data=(void *)m;
+  res->m[1].rtyp=INTVEC_CMD;
+  res->m[1].data=(void *)v;
+  return res;
+}
+
+lists smCallNewBareiss(ideal I, int x, int y)
+{
+  ring origR;
+  sip_sring tmpR;
+  lists res=(lists)Alloc(sizeof(slists));
+  ideal II=smRingCopy(I,&origR,tmpR);
   sparse_mat *bareiss = new sparse_mat(II);
   ideal mm=II;
@@ -194 +257 @@
   {
     delete bareiss;
-    if (origR!=NULL)
-    {
-      rChangeCurrRing(origR,TRUE);
-      rUnComplete(&tmpR);
-      Free((ADDRESS)tmpR.order,3*sizeof(int));
-      Free((ADDRESS)tmpR.block0,3*sizeof(int));
-      Free((ADDRESS)tmpR.block1,3*sizeof(int));
-    }
+    if (origR!=NULL) smRingClean(origR,tmpR);
     v=new intvec(1,pVariables);
   }
@@ -217 +273 @@
       mm=idInit(IDELEMS(m),m->rank);
       int k;
-      for (k=0;k<IDELEMS(m);k++) mm->m[k] = pFetchCopy(origR, m->m[k]);
+      for (k=0;k<IDELEMS(m);k++) mm->m[k] = pFetchCopy(&tmpR, m->m[k]);
       rChangeCurrRing(&tmpR,FALSE);
       idDelete(&m);
-      rChangeCurrRing(origR,TRUE);
-      rUnComplete(&tmpR);
-      Free((ADDRESS)tmpR.order,3*sizeof(int));
-      Free((ADDRESS)tmpR.block0,3*sizeof(int));
-      Free((ADDRESS)tmpR.block1,3*sizeof(int));
+      smRingClean(origR,tmpR);
     }
     else
@@ -259 +311 @@
   tored = nrows; // without border
   i = tored+1;
-  perm = (unsigned short *)Alloc(sizeof(unsigned short)*i);
+  perm = (unsigned short *)Alloc(sizeof(unsigned short)*(i+1));
+  perm[i] = 0;
   m_row = (smpoly *)Alloc0(sizeof(smpoly)*i);
   wrw = (float *)Alloc(sizeof(float)*i);
@@ -295 +348 @@
   Free((ADDRESS)wrw, sizeof(float)*i);
   Free((ADDRESS)m_row, sizeof(smpoly)*i);
-  Free((ADDRESS)perm, sizeof(unsigned short)*i);
+  Free((ADDRESS)perm, sizeof(unsigned short)*(i+1));
 }
 
@@ -322 +375 @@
 }
 
+/* ---------------- the algorithm's ------------------ */
+/*
+* the determinant (up to sign), uses new Bareiss elimination
+*/
+poly sparse_mat::smDet()
+{
+  int y;
+  poly res = NULL;
+
+  if (sign == 0)
+  {
+    this->smActDel();
+    return NULL;
+  }
+  if (act < 2)
+  {
+    if (act != 0) res = m_act[1]->m;
+    Free((void *)m_act[1], sizeof(smprec));
+    return res;
+  }
+  this->smPivot();
+  this->smSign();
+  this->smSelectPR();
+  this->sm1Elim();
+  crd++;
+  this->smColDel();
+  act--;
+  this->smZeroElim();
+  if (sign == 0)
+  {
+    this->smActDel();
+    return NULL;
+  }
+  if (act < 2)
+  {
+    if (act != 0) res = m_act[1]->m;
+    Free((void *)m_act[1], sizeof(smprec));
+    return res;
+  }
+  loop
+  {
+    m_res[crd] = piv;
+    this->smNewPivot();
+    this->smSign();
+    this->smSelectPR();
+    this->smMultCol();
+    this->smHElim();
+    crd++;
+    this->smColDel();
+    act--;
+    this->smZeroElim();
+    if (sign == 0)
+    {
+      this->smPivDel();
+      this->smActDel();
+      return NULL;
+    }
+    if (act < 2)
+    {
+      this->smPivDel();
+      if (act != 0) res = m_act[1]->m;
+      Free((void *)m_act[1], sizeof(smprec));    
+      return res;
+    }
+  }
+}
+
 /*
 * the Bareiss elimination:
@@ -339 +459 @@
     if (act != 0)
     {
-      oldpiv = NULL;
       this->smCopToRes();
     }
@@ -402 +521 @@
     if (act != 0)
     {
-      oldpiv = NULL;
       this->smCopToRes();
     }
@@ -1302 +1420 @@
 }
 
+/*
+* delete the m_act finaly
+*/
+void sparse_mat::smActDel()
+{
+  smpoly a;
+  int i;
+
+  for (i=act; i; i--)
+  {
+    a = m_act[i];
+    do
+    {
+      smElemDelete(&a);
+    } while (a != NULL);
+  }
+}
+
+/*
+* delete the pivotcol
+*/
+void sparse_mat::smColDel()
+{
+  smpoly a = m_act[act];
+
+  while (a != NULL)
+  {
+    smElemDelete(&a);
+  }
+}
+
+/*
+* delete pivot elements
+*/
+void sparse_mat::smPivDel()
+{
+  int i=crd-1;
+
+  while (i != 0)
+  {
+    smElemDelete(&m_res[i]);
+    i--;
+  }
+}
+
+/*
+* the sign of the determinant
+*/
+void sparse_mat::smSign()
+{
+  int j,i=1;
+  if (cpiv!=act) sign=-sign;
+  if ((act%2)==0) sign=-sign;
+  j=perm[1];
+  while(j<rpiv)
+  {
+    sign=-sign;
+    i++;
+    j=perm[i];
+  }
+  while(perm[i]!=0)
+  {
+    perm[i]=perm[i+1];
+    i++;
+  }
+}
+
 /* ----------------- arithmetic ------------------ */
 
@@ -1486 +1671 @@
 }
 
-/*2
+/*n
 * exact division a/b
 * a is a result of smMultDiv

Singular/sparsmat.h

@@ -8 +8 @@
  *
  *******************************************************************/
-/* $Id: sparsmat.h,v 1.1 1999-02-10 16:00:06 Singular Exp $ */
+/* $Id: sparsmat.h,v 1.2 1999-06-02 14:09:36 pohl Exp $ */
 
 
@@ -16 +16 @@
 void smSpecialPolyDiv(poly, poly);
 
+poly smCallDet(ideal I);
 lists smCallBareiss(ideal smat, int x, int y);
 lists smCallNewBareiss(ideal smat, int x, int y);