Changeset 1aa55bf in git

Timestamp:

Oct 16, 2000, 2:06:41 PM (24 years ago)

Author:

Olaf Bachmann <obachman@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

dc078459b7830eda29e3e82e01acd64423bbddc0

Parents:

b559f83182d4517c4869319c1caa185e21439112

Message:

* sparsemat.cc : faster smSelectCopy and relatives
* ring: VarL_* stuff
* divisiblity tests based on longs


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

Location:

Singular

Files:

• Makefile.in (modified) (2 diffs)
• extra.cc (modified) (2 diffs)
• kspoly.cc (modified) (6 diffs)
• kstd2.cc (modified) (4 diffs)
• kutil.cc (modified) (2 diffs)
• kutil.h (modified) (3 diffs)
• pDebug.cc (modified) (3 diffs)
• pInline1.h (modified) (4 diffs)
• pInline2.h (modified) (2 diffs)
• p_MemAdd.h (modified) (2 diffs)
• p_Procs.cc (modified) (8 diffs)
• p_Procs.h (modified) (3 diffs)
• p_polys.h (modified) (3 diffs)
• polys.cc (modified) (2 diffs)
• polys.h (modified) (3 diffs)
• pp_Mult_nn__Template.cc (modified) (1 diff)
• ring.cc (modified) (10 diffs)
• ring.h (modified) (3 diffs)
• sparsmat.cc (modified) (6 diffs)
• structs.h (modified) (3 diffs)

Singular/Makefile.in

@@ -116 +116 @@
         pp_Mult_mm__Template.cc p_Mult_mm__Template.cc \
         p_Minus_mm_Mult_qq__Template.cc p_Add_q__Template.cc \
-        p_Neg__Template.cc kInline.cc
+        p_Neg__Template.cc pp_Mult_Coeff_mm_DivSelect__Template.cc \
+        kInline.cc utils.cc utils.h 
 
 
@@ -129 +130 @@
         fglm.h kstd1.h modulop.h sing_dbm.h weight.h \
         fglmgauss.h fglmvec.h mpsr.h sing_mp.h \
-        kstdfac.h mpsr_Get.h \
+        kstdfac.h mpsr_Get.h kmatrix.h\
         kutil.h mpsr_Put.h sing_dld.h\
         ndbm.h polys-impl.h libparse.h \

Singular/extra.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR      *
 *****************************************/
-/* $Id: extra.cc,v 1.146 2000-09-25 12:26:30 obachman Exp $ */
+/* $Id: extra.cc,v 1.147 2000-10-16 12:06:34 obachman Exp $ */
 /*
 * ABSTRACT: general interface to internals of Singular ("system" command)
@@ -741 +741 @@
 #endif
 /*==================== pDivStat =============================*/
-#ifdef PDEBUG
+#if defined(PDEBUG) || defined(PDIV_DEBUG)
     if(strcmp(sys_cmd,"pDivStat")==0)
     {

Singular/kspoly.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: kspoly.cc,v 1.11 2000-10-04 13:12:01 obachman Exp $ */
+/* $Id: kspoly.cc,v 1.12 2000-10-16 12:06:34 obachman Exp $ */
 /*
 *  ABSTRACT -  Routines for Spoly creation and reductions
@@ -22 +22 @@
 
 
-#define NEW_STUFF
-
-#ifdef NEW_STUFF
 /***************************************************************
  *
@@ -99 +96 @@
   p_LmDelete(lm, tailRing);
 }
-#endif
 
 /***************************************************************
@@ -107 +103 @@
  *
  ***************************************************************/
-#ifdef NEW_STUFF
 void ksCreateSpoly(LObject* Pair,poly spNoether, 
                    int use_buckets, ring tailRing)
@@ -189 +184 @@
   }
 }
-#endif
-
-//////////////////////////////////////////////////////////////////////////
-// Reduces PR at Current->next with PW
-// Assumes PR != NULL, Current contained in PR 
-//         Current->next != NULL, LM(PW) devides LM(Current->next)
-// Changes: PR
-// Const:   PW
-#ifdef NEW_STUFF
-void ksSpolyTail(LObject* PR, TObject* PW, poly Current, poly spNoether)
-{
-  poly Lp = PR->p;
-  number coef;
-  poly Save = PW->p;
-  ring lmRing = PR->lmRing;
-  
-  assume(Lp != NULL && Current != NULL && pNext(Current) != NULL);
-  assume(PR->bucket == NULL);
-  pAssume(pIsMonomOf(Lp, Current));
-
-  if (Lp == Save)
-    PW->p = p_Copy(Save, currRing, PR->tailRing);
-
-  
-  PR->p = pNext(Current);
-  PR->lmRing = PR->tailRing;
-  ksReducePoly(PR, PW, spNoether, &coef);
-  
-  if (! n_IsOne(coef, currRing))
-  {
-    pNext(Current) = NULL;
-    p_Mult_nn(Lp, coef, currRing, PR->tailRing);
-  }
-  n_Delete(&coef, currRing);
-  pNext(Current) = PR->p;
-  PR->p = Lp;
-  PR->lmRing = lmRing;
-  if (PW->p != Save)
-  {
-    p_Delete(&(PW->p), currRing, PR->tailRing);
-    PW->p = Save; // == Lp
-  }
-}
-
-#endif
-
-#ifndef NEW_STUFF
-
-void ksReducePoly(LObject* PR,
-                  TObject* PW,
-                  poly spNoether,
-                  number *coef)
-{
-  poly p1 = PR->p;
-  poly p2 = PW->p;
-
-  assume(p2 != NULL && p1 != NULL && pDivisibleBy(p2,  p1));
-  assume(pGetComp(p1) == pGetComp(p2) || 
-         (pMaxComp(p2) == 0));
-  
-  poly a2 = pNext(p2), lm = p1;
-
-  p1 = pNext(p1);
-
-  if (a2==NULL)
-  {
-    pDeleteLm(&lm);
-    PR->p = p1;
-    if (coef != NULL) *coef = nInit(1);
-    return;
-  }
-
-  if (! nIsOne(pGetCoeff(p2)))
-  {
-    number bn = pGetCoeff(lm);
-    number an = pGetCoeff(p2);
-    int ct = ksCheckCoeff(&an, &bn);
-    pSetCoeff(lm, bn);
-    if ((ct == 0) || (ct == 2)) 
-      p1 = pMult_nn(p1, an);
-    if (coef != NULL) *coef = an;
-    else nDelete(&an);
-  }
-  else
-  {
-    if (coef != NULL) *coef = nInit(1);
-  }
-  
-  
-  // pMonSubFrom(lm, p2);
-  pExpVectorSub(lm, p2);
-  int dummy;
-  PR->p = currRing->p_Procs->p_Minus_mm_Mult_qq(p1, lm, a2, 
-                                                dummy, spNoether, currRing);
-
-  pDeleteLm(&lm);
-}
-#endif
-
-#ifndef NEW_STUFF
-
-/***************************************************************
- *
- * Creates S-Poly of p1 and p2
- * 
- *
- ***************************************************************/
-void ksCreateSpoly(LObject* Pair, 
-                   poly spNoether,
-                   int use_buckets, ring tailRing)
-{
-  assume(kTest_L(Pair));
-  poly p1 = Pair->p1;
-  poly p2 = Pair->p2;
-
-  assume(p1 != NULL);
-  assume(p2 != NULL);
-
-  poly a1 = pNext(p1), a2 = pNext(p2);
-  number lc1 = pGetCoeff(p1), lc2 = pGetCoeff(p2);
-  poly m1, m2;
-  int co=0, ct = ksCheckCoeff(&lc1, &lc2);
-  int x, l1;
-
-  if (pGetComp(p1)!=pGetComp(p2))
-  {
-    if (pGetComp(p1)==0)
-    {
-      co=1;
-      pSetCompP(p1,pGetComp(p2));
-    }
-    else
-    {
-      co=2;
-      pSetCompP(p2,pGetComp(p1));
-    }
-  }
-
-  // get m1 = LCM(LM(p1), LM(p2))/LM(p1)
-  //     m2 = LCM(LM(p1), LM(p2))/LM(p2)
-  m1 = pInit();
-  m2 = pInit();
-  for (int i = pVariables; i; i--)
-  {
-    x = pGetExpDiff(p1, p2, i);
-    if (x > 0)
-    {
-      pSetExp(m2,i,x);
-      pSetExp(m1,i,0);
-    }
-    else
-    {
-      pSetExp(m1,i,-x);
-      pSetExp(m2,i,0);
-    }
-  }
-  pSetm(m1);
-  pSetm(m2);           // now we have m1 * LM(p1) == m2 * LM(p2)
-  pSetCoeff0(m1, lc2);
-  pSetCoeff0(m2, lc1); // and now, m1 * LT(p1) == m2 * LT(p2)
-
-  // get m2 * a2
-  a2 = currRing->p_Procs->pp_Mult_mm(a2, m2, spNoether, currRing);
-
-  // and, finally, the spoly
-  int dummy;
-  Pair->p = currRing->p_Procs->p_Minus_mm_Mult_qq(a2, m1, a1, 
-                                                  dummy, spNoether,
-                                                  currRing);
-  
-  // Clean-up time
-  pDeleteLm(&m1);
-  pDeleteLm(&m2);
-  
-  if (co != 0)
-  {
-    if (co==1)
-    {
-      pSetCompP(p1,0);
-    }
-    else
-    {
-      pSetCompP(p2,0);
-    }
-  }
-}
-#endif
-
-#ifndef NEW_STUFF
+
 //////////////////////////////////////////////////////////////////////////
 // Reduces PR at Current->next with PW
@@ -389 +196 @@
   number coef;
   poly Save = PW->p;
-  
+  ring lmRing = PR->lmRing;
+  
+  assume(Lp != NULL && Current != NULL && pNext(Current) != NULL);
+  assume(PR->bucket == NULL);
+  pAssume(pIsMonomOf(Lp, Current));
+
   if (Lp == Save)
-    PW->p = pCopy(Save);
-    
-  assume(Lp != NULL && Current != NULL && pNext(Current) != NULL);
-  assume(pIsMonomOf(Lp, Current));
+    PW->p = p_Copy(Save, currRing, PR->tailRing);
+
   
   PR->p = pNext(Current);
+  PR->lmRing = PR->tailRing;
   ksReducePoly(PR, PW, spNoether, &coef);
   
-  if (! nIsOne(coef))
+  if (! n_IsOne(coef, currRing))
   {
     pNext(Current) = NULL;
-    pMult_nn(Lp, coef);
-  }
-  nDelete(&coef);
+    p_Mult_nn(Lp, coef, currRing, PR->tailRing);
+  }
+  n_Delete(&coef, currRing);
   pNext(Current) = PR->p;
   PR->p = Lp;
+  PR->lmRing = lmRing;
   if (PW->p != Save)
   {
-    pDelete(&(PW->p));
+    p_Delete(&(PW->p), currRing, PR->tailRing);
     PW->p = Save; // == Lp
   }
 }
-
-#endif
 
 /***************************************************************

Singular/kstd2.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: kstd2.cc,v 1.49 2000-10-04 13:12:02 obachman Exp $ */
+/* $Id: kstd2.cc,v 1.50 2000-10-16 12:06:35 obachman Exp $ */
 /*
 *  ABSTRACT -  Kernel: alg. of Buchberger
@@ -188 +188 @@
 
   poly pi;
-  int i,j,at,reddeg,d,pass,ei, ii;
+  int i,j,at,reddeg,d,pass,ei, ii, h_d;
   unsigned long not_sev;
 
@@ -300 +300 @@
     h->sev = pGetShortExpVector(h->p);
     not_sev = ~ h->sev;
+    h_d = pFDeg(h->p);
     /* compute the ecart */
     if (ei <= (*h).ecart)
-      (*h).ecart = d-pFDeg((*h).p);
+      (*h).ecart = d-h_d;
     else
-      (*h).ecart = d-pFDeg((*h).p)+ei-(*h).ecart;
+      (*h).ecart = d-h_d+ei-(*h).ecart;
     /*
      * try to reduce the s-polynomial h
@@ -312 +313 @@
      */
     pass++;
-    d = pFDeg((*h).p)+(*h).ecart;
+    d = h_d +(*h).ecart;
     if ((strat->Ll >= 0) && ((d > reddeg) || (pass > strat->LazyPass)))
     {

Singular/kutil.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: kutil.cc,v 1.64 2000-10-04 13:12:03 obachman Exp $ */
+/* $Id: kutil.cc,v 1.65 2000-10-16 12:06:35 obachman Exp $ */
 /*
 * ABSTRACT: kernel: utils for kStd
@@ -27 +27 @@
 #include "kutil.h"
 
-/* Hmm ... this should be inlined or made more efficient:
-   see Long/mregular.tst */
-/*2
-*should return 1 if p divides q and p<q,
-*             -1 if q divides p and q<p
-*              0 otherwise
-*/
-static inline int     pDivComp(poly p, poly q)
-{
-  if (pGetComp(p) == pGetComp(q))
-  {
-    int i=pVariables;
-    long d;
-    BOOLEAN a=FALSE, b=FALSE;
-    for (; i>0; i--)
-    {
-      d = pGetExpDiff(p, q, i);
-      if (d)
-      {
-        if (d < 0)
-        {
-          if (b) return 0;
-          a =TRUE;
-        }
-        else
-        {
-          if (a) return 0;
-          b = TRUE;
-        }
-      }
-    }
-    if (a) return 1;
-    else if (b)  return -1;
-  }
-  return 0;
-}
 
 static poly redMora (poly h,int maxIndex,kStrategy strat);
 static poly redBba (poly h,int maxIndex,kStrategy strat);
+
+static inline int pDivComp(poly p, poly q)
+{
+  if (pGetComp(p) == pGetComp(q))
+  {
+    BOOLEAN a=FALSE, b=FALSE;
+    int i;
+    unsigned long la, lb;
+    unsigned long divmask = currRing->divmask;
+    for (i=0; i<currRing->VarL_Size; i++)
+    {
+      la = p->exp[currRing->VarL_Offset[i]];
+      lb = q->exp[currRing->VarL_Offset[i]];
+      if (la != lb)
+      {
+        if (la < lb)
+        {
+          if (b) return 0;
+          if (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask))
+            return 0;
+          a = TRUE;
+        }
+        else
+        {
+          if (a) return 0;
+          if (((la & divmask) ^ (lb & divmask)) != ((la - lb) & divmask))
+            return 0;
+          b = TRUE;
+        }
+      }
+    }
+    if (a) return 1;
+    if (b) return -1;
+  }
+  return 0;
+}
+
 
 BITSET  test=(BITSET)0;

Singular/kutil.h

@@ -4 +4 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: kutil.h,v 1.32 2000-10-04 13:12:03 obachman Exp $ */
+/* $Id: kutil.h,v 1.33 2000-10-16 12:06:36 obachman Exp $ */
 /*
 * ABSTRACT: kernel: utils for kStd
@@ -34 +34 @@
     { 
       memset((void*) this, 0, sizeof(sTObject)); 
+    }
+  sTObject(poly _p)
+    {
+      memset((void*) this, 0, sizeof(sTObject)); 
+      p = _p;
     }
 };
@@ -52 +57 @@
       memset((void*) this, 0, sizeof(sLObject));
       lmRing = tailRing = currRing;
+    }
+  sLObject(poly _p) 
+    { 
+      memset((void*) this, 0, sizeof(sLObject));
+      lmRing = tailRing = currRing;
+      p = _p;
     }
   // spoly related things

Singular/pDebug.cc

@@ -7 +7 @@
  *  Author:  obachman (Olaf Bachmann)
  *  Created: 8/00
- *  Version: $Id: pDebug.cc,v 1.5 2000-09-25 12:26:35 obachman Exp $
+ *  Version: $Id: pDebug.cc,v 1.6 2000-10-16 12:06:36 obachman Exp $
  *******************************************************************/
 
@@ -100 +100 @@
  *
  ***************************************************************/
-static int pDivisibleBy_number = 1;
-static int pDivisibleBy_FALSE = 1;
-static int pDivisibleBy_ShortFalse = 1;
-BOOLEAN pDebugLmShortDivisibleBy(poly p1, unsigned long sev_1, ring r_1,
-                               poly p2, unsigned long not_sev_2, ring r_2)
-{
-  _pPolyAssume(p_GetShortExpVector(p1, r_1) == sev_1, p1, r_1);
-  _pPolyAssume(p_GetShortExpVector(p2, r_2) == ~ not_sev_2, p2, r_2);
-
-  pDivisibleBy_number++;
-  BOOLEAN ret = p_LmDivisibleBy(p1, r_1, p2, r_2);
-  if (! ret) pDivisibleBy_FALSE++;
-  if (sev_1 & not_sev_2)
-  {
-    pDivisibleBy_ShortFalse++;
-    if (ret)
-      dReportError("p1 divides p2, but sev's are wrong");
-  }
-  return ret;
-}
-
-void pPrintDivisbleByStat()
-{
-  Print("#Tests: %d; #FALSE %d(%d); #SHORT %d(%d)\n",
-        pDivisibleBy_number,
-        pDivisibleBy_FALSE, pDivisibleBy_FALSE*100/pDivisibleBy_number,
-        pDivisibleBy_ShortFalse, pDivisibleBy_ShortFalse*100/pDivisibleBy_FALSE);
-}
+BOOLEAN p_DebugLmDivisibleByNoComp(poly a, poly b, ring r)
+{
+  int i=r->N;
+
+  do
+  {
+    if (p_GetExp(a,i,r) > p_GetExp(b,i,r))
+      return FALSE;
+    i--;
+  }
+  while (i);
+  return TRUE;
+}
+
 
 /***************************************************************
@@ -297 +283 @@
   
 #endif // PDEBUG
+
+#include "pInline1.h"
+
+#if defined(PDEBUG) || defined(PDIV_DEBUG)
+static unsigned long pDivisibleBy_number = 1;
+static unsigned long pDivisibleBy_FALSE = 1;
+static unsigned long pDivisibleBy_ShortFalse = 1;
+BOOLEAN pDebugLmShortDivisibleBy(poly p1, unsigned long sev_1, ring r_1,
+                               poly p2, unsigned long not_sev_2, ring r_2)
+{
+  _pPolyAssume(p_GetShortExpVector(p1, r_1) == sev_1, p1, r_1);
+  _pPolyAssume(p_GetShortExpVector(p2, r_2) == ~ not_sev_2, p2, r_2);
+
+  pDivisibleBy_number++;
+  BOOLEAN ret;
+  if (r_1 == r_2)
+    ret = p_LmDivisibleBy(p1, p2, r_1);
+  else
+    ret = p_LmDivisibleBy(p1, r_1, p2, r_2);
+
+  if (! ret) pDivisibleBy_FALSE++;
+  if (sev_1 & not_sev_2)
+  {
+    pDivisibleBy_ShortFalse++;
+    if (ret)
+      dReportError("p1 divides p2, but sev's are wrong");
+  }
+  return ret;
+}
+
+void pPrintDivisbleByStat()
+{
+  Print("#Tests: %d; #FALSE %d(%d); #SHORT %d(%d)\n",
+        pDivisibleBy_number,
+        pDivisibleBy_FALSE, (unsigned long) ((double)pDivisibleBy_FALSE*((double) 100)/(double)pDivisibleBy_number),
+        pDivisibleBy_ShortFalse, (unsigned long) ((double)pDivisibleBy_ShortFalse*((double)100)/(double)pDivisibleBy_FALSE));
+}
+#endif
+
 #endif // PDEBUG_CC
 

Singular/pInline1.h

@@ -7 +7 @@
  *  Author:  obachman (Olaf Bachmann)
  *  Created: 8/00
- *  Version: $Id: pInline1.h,v 1.5 2000-10-04 13:12:04 obachman Exp $
+ *  Version: $Id: pInline1.h,v 1.6 2000-10-16 12:06:36 obachman Exp $
  *******************************************************************/
 #ifndef PINLINE1_H
@@ -44 +44 @@
 #ifdef PDIV_DEBUG
 BOOLEAN pDebugLmShortDivisibleBy(poly p1, unsigned long sev_1, ring r_1,
-                               poly p2, unsigned long not_sev_2, ring r_2);
+                                 poly p2, unsigned long not_sev_2, ring r_2);
+BOOLEAN p_DebugLmDivisibleByNoComp(poly a, poly b, ring r);
+#define pDivAssume  pAssume
+#else
+#define pDivAssume(x)   ((void)0)
 #endif
 
@@ -149 +153 @@
   p_MemSub_LengthGeneral(p1->exp, p2->exp, r->ExpLSize);
 }
+// ExpVector(p1) += ExpVector(p2) - ExpVector(p3)
+PINLINE1 void p_ExpVectorAddSub(poly p1, poly p2, poly p3, ring r)
+{
+  p_CheckPolyRing1(p1, r);
+  p_CheckPolyRing1(p2, r);
+  p_CheckPolyRing1(p3, r);
+#if PDEBUG >= 1
+  for (int i=1; i<=r->N; i++)
+    pAssume1(p_GetExp(p1, i, r) + p_GetExp(p2, i, r) >= p_GetExp(p3, i, r));
+  pAssume1(p_GetComp(p1, r) == 0 || 
+           (p_GetComp(p2, r) - p_GetComp(p3, r) == 0) ||
+           (p_GetComp(p1, r) == p_GetComp(p2, r) - p_GetComp(p3, r)));
+#endif 
+ 
+  p_MemAddSub_LengthGeneral(p1->exp, p2->exp, p3->exp, r->ExpLSize);
+}
 // ExpVector(pr) = ExpVector(p1) + ExpVector(p2)
 PINLINE1 void p_ExpVectorSum(poly pr, poly p1, poly p2, ring r)
@@ -249 +269 @@
  *
  ***************************************************************/
+// return: FALSE, if there exists i, such that a->exp[i] > b->exp[i]
+//         TRUE, otherwise
+// (1) Consider long vars, instead of single exponents
+// (2) Clearly, if la > lb, then FALSE
+// (3) Suppose la <= lb, and consider first bits of single exponents in l: 
+//     if TRUE, then value of these bits is la ^ lb
+//     if FALSE, then la-lb causes an "overflow" into one of those bits, i.e.,
+//               la ^ lb != la - lb
 static inline BOOLEAN _p_LmDivisibleByNoComp(poly a, poly b, ring r)
 {
-  int i=r->N;
-
-  do
-  {
-    if (p_GetExp(a,i,r) > p_GetExp(b,i,r))
-      return FALSE;
-    i--;
-  }
-  while (i);
+  int i=r->VarL_Size - 1;
+  unsigned long divmask = r->divmask;
+  unsigned long la, lb;
+  
+  if (r->VarL_LowIndex >= 0)
+  {
+    i += r->VarL_LowIndex;
+    do
+    {
+      la = a->exp[i];
+      lb = b->exp[i];
+      if ((la > lb) || 
+          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
+      {
+        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
+        return FALSE;
+      }
+      i--;
+    }
+    while (i>=r->VarL_LowIndex);
+  }
+  else
+  {
+    do
+    {
+      la = a->exp[r->VarL_Offset[i]];
+      lb = b->exp[r->VarL_Offset[i]];
+      if ((la > lb) || 
+          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
+      {
+        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
+        return FALSE;
+      }
+      i--;
+    }
+    while (i>=0);
+  }
+  pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == TRUE);
   return TRUE;
 }
+
 static inline BOOLEAN _p_LmDivisibleByNoComp(poly a, ring r_a, poly b, ring r_b)
 {

Singular/pInline2.h

@@ -7 +7 @@
  *  Author:  obachman (Olaf Bachmann)
  *  Created: 8/00
- *  Version: $Id: pInline2.h,v 1.7 2000-10-04 13:12:04 obachman Exp $
+ *  Version: $Id: pInline2.h,v 1.8 2000-10-16 12:06:37 obachman Exp $
  *******************************************************************/
 #ifndef PINLINE2_H
@@ -387 +387 @@
 }
 
+PINLINE2 poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r)
+{
+  return r->p_Procs->pp_Mult_Coeff_mm_DivSelect(p, m, r);
+}
+
 // returns -p, destroys p
 PINLINE2 poly p_Neg(poly p, const ring r)

Singular/p_MemAdd.h

@@ -7 +7 @@
  *  Author:  obachman (Olaf Bachmann)
  *  Created: 8/00
- *  Version: $Id: p_MemAdd.h,v 1.2 2000-09-12 16:01:06 obachman Exp $
+ *  Version: $Id: p_MemAdd.h,v 1.3 2000-10-16 12:06:37 obachman Exp $
  *******************************************************************/
 #ifndef P_MEM_ADD_H
@@ -366 +366 @@
 while (0)
 
+#define _p_MemAddSub_Declare(r, s, t)                \
+  const unsigned long* _s = ((unsigned long*) s); \
+  const unsigned long* _t = ((unsigned long*) t); \
+  unsigned long* _r = ((unsigned long*) r)
+
+#define p_MemAddSub_LengthGeneral(r, s, t, length)  \
+do                                                  \
+{                                                   \
+  _p_MemAddSub_Declare(r,s, t);                     \
+  const unsigned long _l = (unsigned long) length;  \
+  unsigned long _i = 0;                             \
+                                                    \
+  do                                                \
+  {                                                 \
+    _r[_i] += _s[_i] - _t[_i];                      \
+    _i++;                                           \
+  }                                                 \
+  while (_i != _l);                                 \
+}                                                   \
+while (0)
+
 #endif P_MEM_ADD_H
+

Singular/p_Procs.cc

@@ -7 +7 @@
  *  Author:  obachman (Olaf Bachmann)
  *  Created: 8/00
- *  Version: $Id: p_Procs.cc,v 1.14 2000-10-04 15:37:54 obachman Exp $
+ *  Version: $Id: p_Procs.cc,v 1.15 2000-10-16 12:06:37 obachman Exp $
  *******************************************************************/
 #include <string.h>
@@ -46 +46 @@
 //   5 -- all Field*_Length*_Ord* procs
 #ifdef NDEBUG
-const int HAVE_FAST_P_PROCS = 2;
+const int HAVE_FAST_P_PROCS = 3;
 #else
 const int HAVE_FAST_P_PROCS = 0;
@@ -64 +64 @@
 //   3 -- special cases for length <= 4
 //   4 -- special cases for length <= 8
-const int HAVE_FAST_LENGTH = 4;
+const int HAVE_FAST_LENGTH = 3;
 
 // Set HAVE_FAST_ORD to:
@@ -180 +180 @@
   p_Minus_mm_Mult_qq_Proc,
   p_Neg_Proc,
+  pp_Mult_Coeff_mm_DivSelect_Proc,
   p_Unknown_Proc
 };
@@ -268 +269 @@
       case p_Minus_mm_Mult_qq_Proc: return "p_Minus_mm_Mult_qq_Proc";
       case p_Neg_Proc: return "p_Neg_Proc";
+      case pp_Mult_Coeff_mm_DivSelect_Proc: return "pp_Mult_Coeff_mm_DivSelect_Proc";
       case p_Unknown_Proc: return "p_Unknown_Proc";
   }
@@ -483 +485 @@
       case p_Mult_mm_Proc:
       case pp_Mult_nn_Proc:
+      case pp_Mult_Coeff_mm_DivSelect_Proc:
         return index(field, length);
 
@@ -636 +639 @@
     (p_Procs->p_Add_q != NULL) &&
     (p_Procs->p_Neg != NULL) &&
+    (p_Procs->pp_Mult_Coeff_mm_DivSelect != NULL) &&
     (p_Procs->p_Minus_mm_Mult_qq != NULL));
 }
@@ -901 +905 @@
   SetProc(p_Minus_mm_Mult_qq, field, length, ord);
   SetProc(p_Neg, field, LengthGeneral, OrdGeneral);
-}
-
-
-
+  SetProc(pp_Mult_Coeff_mm_DivSelect, field, length, OrdGeneral);
+}
+
+
+

Singular/p_Procs.h

@@ -8 +8 @@
  *  Author:  obachman (Olaf Bachmann)
  *  Created: 8/00
- *  Version: $Id: p_Procs.h,v 1.5 2000-09-18 09:19:27 obachman Exp $
+ *  Version: $Id: p_Procs.h,v 1.6 2000-10-16 12:06:38 obachman Exp $
  *******************************************************************/
 #ifndef P_PROCS_H
@@ -28 +28 @@
                                                 const ring r);
 typedef poly (*p_Neg_Proc_Ptr)(poly p, const ring r);
+typedef poly (*pp_Mult_Coeff_mm_DivSelect_Proc_Ptr)(poly p, const poly m, const ring r);
 
 typedef struct p_Procs_s
@@ -41 +42 @@
   p_Minus_mm_Mult_qq_Proc_Ptr   p_Minus_mm_Mult_qq;
   p_Neg_Proc_Ptr                p_Neg;
+  pp_Mult_Coeff_mm_DivSelect_Proc_Ptr pp_Mult_Coeff_mm_DivSelect;
 } pProcs_s;
 

Singular/p_polys.h

@@ -8 +8 @@
  *  Author:  obachman (Olaf Bachmann)
  *  Created: 9/00
- *  Version: $Id: p_polys.h,v 1.3 2000-10-04 13:12:05 obachman Exp $
+ *  Version: $Id: p_polys.h,v 1.4 2000-10-16 12:06:38 obachman Exp $
  *******************************************************************/
 #ifndef P_POLYS_H
@@ -122 +122 @@
 // ExpVector(p1) -= ExpVector(p2)
 PINLINE1 void p_ExpVectorSub(poly p1, poly p2, ring r);
+// ExpVector(p1) += ExpVector(p2) - ExpVector(p3)
+PINLINE1 void p_ExpVectorAddSub(poly p1, poly p2, poly p3, ring r);
 // ExpVector(pr) = ExpVector(p1) + ExpVector(p2)
 PINLINE1 void p_ExpVectorSum(poly pr, poly p1, poly p2, ring r);
@@ -232 +234 @@
 PINLINE2 poly pp_Mult_qq(poly p, poly q, const ring r);
 
+// returns p*Coeff(m) for such monomials pm of p, for which m is divisble by pm
+PINLINE2 poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r);
+
 /***************************************************************
  *

Singular/polys.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: polys.cc,v 1.62 2000-09-18 09:19:29 obachman Exp $ */
+/* $Id: polys.cc,v 1.63 2000-10-16 12:06:38 obachman Exp $ */
 
 /*
@@ -780 +780 @@
     }
   }
-  if (pFDeg!=pWTotaldegree) pFDeg=pTotaldegree;
+  if (pFDeg!=pWTotaldegree) 
+  {
+    if (rOrd_is_Totaldegree_Ordering(r))
+    {
+      pFDeg = pDeg;
+    }
+    else
+    {
+      pFDeg=pTotaldegree;
+    }
+  }
 }
 

Singular/polys.h

@@ -4 +4 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: polys.h,v 1.39 2000-09-20 12:56:37 obachman Exp $ */
+/* $Id: polys.h,v 1.40 2000-10-16 12:06:39 obachman Exp $ */
 /*
 * ABSTRACT - all basic methods to manipulate polynomials of the
@@ -125 +125 @@
 #define pExpVectorAdd(p1, p2)       p_ExpVectorAdd(p1, p2, currRing)
 #define pExpVectorSub(p1, p2)       p_ExpVectorSub(p1, p2, currRing)
+#define pExpVectorAddSub(p1, p2, p3)p_ExpVectorAddSub(p1, p2, p3, currRing)
 #define pExpVectorSum(pr, p1, p2)   p_ExpVectorSum(pr, p1, p2, currRing)
 #define pExpVectorDiff(pr, p1, p2)  p_ExpVectorDiff(pr, p1, p2, currRing)
@@ -206 +207 @@
 #define pMult(p, q)                 p_Mult_q(p, q, currRing)
 #define ppMult_qq(p, q)             pp_Mult_qq(p, q, currRing)
+// p*Coeff(m) for such monomials pm of p, for which m is divisble by pm
+#define ppMult_Coeff_mm_DivSelect(p, m)   pp_Mult_Coeff_mm_DivSelect(p, m, currRing)
 
 /***************************************************************

Singular/pp_Mult_nn__Template.cc

@@ -3 +3 @@
 ****************************************/
 /***************************************************************
- *  File:    pp_Mult_n__Template.cc
- *  Purpose: template for pp_Mult_n
+ *  File:    pp_Mult_nn__Template.cc
+ *  Purpose: template for pp_Mult_nn
  *  Author:  obachman (Olaf Bachmann)
  *  Created: 8/00
- *  Version: $Id: pp_Mult_nn__Template.cc,v 1.3 2000-09-18 09:19:31 obachman Exp $
+ *  Version: $Id: pp_Mult_nn__Template.cc,v 1.4 2000-10-16 12:06:39 obachman Exp $
  *******************************************************************/
 

Singular/ring.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: ring.cc,v 1.119 2000-09-20 12:33:29 obachman Exp $ */
+/* $Id: ring.cc,v 1.120 2000-10-16 12:06:39 obachman Exp $ */
 
 /*
@@ -63 +63 @@
 // unconditionally deletes fields in r
 static void rDelete(ring r);
+// set r->VarL_Size, r->VarL_Offset, r->VarL_LowIndex
+static void rSetVarL(ring r);
+// get r->divmask depending on bits per exponent
+static unsigned long rGetDivMask(int bits);
 
 /*0 implementation*/
@@ -1935 +1939 @@
 }
 
+BOOLEAN rOrder_is_DegOrdering(rRingOrder_t order)
+{
+  switch(order)
+  {
+      case ringorder_dp:
+      case ringorder_Dp:
+      case ringorder_ds:
+      case ringorder_Ds:
+        return TRUE;
+        
+      default:
+        return FALSE;
+  }
+}
+
+// return TRUE if p->exp[r->pOrdIndex] holds total degree of p */
+BOOLEAN rOrd_is_Totaldegree_Ordering(ring r =currRing)
+{
+  // Hmm.... what about Syz orderings?
+  return (r->N > 1 &&
+          rHasSimpleOrder(r) &&
+          (rOrder_is_DegOrdering((rRingOrder_t)r->order[0]) ||
+           rOrder_is_DegOrdering(( rRingOrder_t)r->order[1])));
+}
+                             
 BOOLEAN rIsPolyVar(int v)
 {
@@ -2210 +2239 @@
 }
 
-unsigned long rGetExpSize(unsigned long bitmask, int & bits)
+static unsigned long rGetExpSize(unsigned long bitmask, int & bits)
 {
   if (bitmask == 0)
@@ -2295 +2324 @@
 * optimize rGetExpSize for a block of N variables, exp <=bitmask
 */
-unsigned long rGetExpSize(unsigned long bitmask, int & bits, int N)
+static unsigned long rGetExpSize(unsigned long bitmask, int & bits, int N)
 {
 #if SIZEOF_LONG == 8
@@ -2449 +2478 @@
   int bits;
   r->bitmask=rGetExpSize(r->bitmask,bits);
+  r->BitsPerExp = bits;
+  r->divmask=rGetDivMask(bits);
+  
   // will be used for ordsgn:
   long *tmp_ordsgn=(long *)omAlloc0(2*(n+r->N)*sizeof(long));
@@ -2784 +2816 @@
   r->p_Procs = (p_Procs_s*)omAlloc(sizeof(p_Procs_s));
   p_SetProcs(r, r->p_Procs);
+
+  // ----------------------------
+  // set VarL_*
+  rSetVarL(r);
   return FALSE;
 }
@@ -2812 +2848 @@
     if (r->p_Procs != NULL)
       omFreeSize(r->p_Procs, sizeof(p_Procs_s));
-  }
+    omfreeSize(r->VarL_Offset, r->VarL_Size*sizeof(int));
+  }
+}
+
+// set r->VarL_Size, r->VarL_Offset, r->VarL_LowIndex
+static void rSetVarL(ring r)
+{
+  poly p = p_Init(r);
+  int* VarL_Offset = (int*) omAlloc0(r->ExpLSize*sizeof(int));
+  int i,j;
+  
+  for (i=1; i<=r->N; i++)
+    p_SetExp(p, i, 1, r);
+
+  r->VarL_LowIndex = 0;
+  for (i=0, j=0; i<r->ExpLSize; i++)
+  {
+    if (p->exp[i] != 0)
+    {
+      VarL_Offset[j] = i;
+      if (j > 0 && VarL_Offset[j-1] != VarL_Offset[j] - 1)
+        r->VarL_LowIndex = -1;
+      j++;
+    }
+  }
+  r->VarL_Size = j;
+  if (r->VarL_LowIndex >= 0)
+    r->VarL_LowIndex = VarL_Offset[0];
+  r->VarL_Offset = (int*) omReallocSize(VarL_Offset,r->ExpLSize*sizeof(int),
+                                        j*sizeof(int));
+  p_LmFree(p, r);
+}
+    
+// get r->divmask depending on bits per exponent
+static unsigned long rGetDivMask(int bits)
+{
+  unsigned long divmask = 1;
+  int i = bits;
+  
+  while (i < BIT_SIZEOF_LONG)
+  {
+    divmask |= (1 << (unsigned long) i);
+    i += bits;
+  }
+  return divmask;
 }
 
@@ -2828 +2908 @@
   for(j=0;j<=r->N;j++) Print("  v%d at e-pos %d, bit %d\n",
      j,r->VarOffset[j] & 0xffffff, r->VarOffset[j] >>24);
+  Print("BitsPerExp=%d\n", r->BitsPerExp);
   Print("bitmask=0x%x\n",r->bitmask);
+  Print("divmask=%p\n", r->divmask);
   PrintS("ordsgn:\n");
   for(j=0;j<r->pCompLSize;j++)
@@ -2854 +2936 @@
   Print("pVarLowIndex:%d ",r->pVarLowIndex);
   Print("pVarHighIndex:%d\n",r->pVarHighIndex);
+  Print("VarL_Size:%d\nVarL_LowIndex:%d\n", r->VarL_Size, r->VarL_LowIndex);
 #ifdef LONG_MONOMS
   Print("pDivLow:%d ",r->pDivLow);

Singular/ring.h

@@ -7 +7 @@
 * ABSTRACT - the interpreter related ring operations
 */
-/* $Id: ring.h,v 1.51 2000-09-19 15:22:25 Singular Exp $ */
+/* $Id: ring.h,v 1.52 2000-10-16 12:06:40 obachman Exp $ */
 
 /* includes */
@@ -210 +210 @@
 // returns TRUE, if simple lp or ls ordering
 BOOLEAN rHasSimpleLexOrder(ring r);
+// return TRUE if p->exp[r->pOrdIndex] holds total degree of p */
+BOOLEAN rOrd_is_Totaldegree_Ordering(ring r =currRing);
 rOrderType_t    rGetOrderType(ring r);
 BOOLEAN rIsPolyVar(int i); /* returns TRUE if var(i) belongs to p-block */
-
-void rOptimizeOrder(ring r);
 
 #ifdef RDEBUG
@@ -222 +222 @@
 #endif
 
-unsigned long rGetExpSize(unsigned long bitmask, int & bits);
-unsigned long rGetExpSize(unsigned long bitmask, int & bits, int N);
 ring rModifyRing(ring r, BOOLEAN omit_degree,
                          BOOLEAN omit_comp,

Singular/sparsmat.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: sparsmat.cc,v 1.40 2000-09-25 14:07:59 pohl Exp $ */
+/* $Id: sparsmat.cc,v 1.41 2000-10-16 12:06:41 obachman Exp $ */
 
 /*
@@ -35 +35 @@
 };
 */
+
+#if OLD > 0
+static poly smSelectCopy(poly, const poly);
+#else
+#define smSelectCopy ppMult_Coeff_mm_DivSelect
+#endif
+
 /* declare internal 'C' stuff */
 static void smExactPolyDiv(poly, poly);
 static BOOLEAN smIsNegQuot(poly, const poly, const poly);
-static BOOLEAN smCheckLead(const poly, const poly);
-static poly smSelectCopy(poly, const poly);
 static void smExpMultDiv(poly, const poly, const poly);
 static void smPolyDivN(poly, const number);
@@ -1805 +1810 @@
     if (smIsNegQuot(e, b, c))
     {
-      lead = smCheckLead(a, e);
+      lead = pLmDivisibleByNoComp(e, a);
       r = smSelectCopy(a, e);
       smExpMultDiv(r, b, c);
@@ -1843 +1848 @@
       r = smSelectCopy(a, e);
       smExpMultDiv(r, b, c);
-      if (smCheckLead(a, e))
+      if (pLmDivisibleByNoComp(e, a))
         smCombineChain(&pa, r);
       else
@@ -1902 +1907 @@
 }
 
+// obachman --> Wilfried: check the following
 static BOOLEAN smIsNegQuot(poly a, const poly b, const poly c)
 {
-  int i=pVariables;
-
-  while ((i>0) && (pGetExp(b,i) >= pGetExp(c,i))) i--;
-  if(i!=0)
-  {
+  if (pLmDivisibleByNoComp(c, b))
+  {
+    pExpVectorDiff(a, b, c);
+    // Hmm: here used to be a pSetm(a): but it is unnecessary,
+    // if b and c are correct
+    return FALSE;
+  }
+  else
+  {
+    int i;
     for (i=pVariables; i>0; i--)
     {
@@ -1916 +1927 @@
         pSetExp(a,i,0);
     }
+    // here we actually might need a pSetm, if a is to be used in
+    // comparisons
     return TRUE;
   }
-  else
-  {
-    for (i=pVariables; i>0; i--)
-    {
-      pSetExp(a,i, pGetExp(b,i)-pGetExp(c,i));
-    }
-    pSetm(a);
-    return FALSE;
-  }
-}
-
-static BOOLEAN smCheckLead(const poly t, const poly e)
+}
+
+static void smExpMultDiv(poly t, const poly b, const poly c)
 {
   int i;
-  Exponent_t w;
-  for (i=pVariables; i; i--)
-  {
-    w = pGetExp(e,i);
-    if (w&&(w>pGetExp(t,i)))
-      return FALSE;
-  }
-  return TRUE;
-}
-
-/*
-* e is a monomial with exponents e(i)
-* t is a polynom with monomials t_j and
-*   exponents t_j(i)
-* make a copy of a part of t:
-*   for the monomials t_j of t in the copy
-*   we have
-*     ((e(i)!=0)&&(e(i)<=t_j(i)) == TRUE
-*/
-static poly smSelectCopy(poly t, const poly e)
-{
-  const number y = pGetCoeff(e);
-  poly res, h;
-
-  loop
-  {
-    if(smCheckLead(t,e)) break;
-    pIter(t);
-    if(t==NULL) return NULL;
-  }
-  h = res = pNew();
-  loop
-  {
-    pExpVectorCopy(h,t);
-    pSetCoeff0(h,nMult(y,pGetCoeff(t)));
-    loop
-    {
-      pIter(t);
-      if(t==NULL)
-      {
-        pNext(h)=NULL;
-        return res;
-      }
-      if(smCheckLead(t,e)) break;
-    }
-    h=pNext(h)=pNew();
-  }
-}
-
-static void smExpMultDiv(poly t, const poly b, const poly c)
-{
-  int i;
+  pTest(t);
+  pLmTest(b);
+  pLmTest(c);
   while(t!=NULL)
   {
-    pExpVectorAdd(t,b);
-    pExpVectorSub(t,c);
+    pExpVectorAddSub(t, b, c);
     pIter(t);
   }

Singular/structs.h

@@ -4 +4 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: structs.h,v 1.37 2000-09-19 15:22:25 Singular Exp $ */
+/* $Id: structs.h,v 1.38 2000-10-16 12:06:41 obachman Exp $ */
 /*
 * ABSTRACT
@@ -252 +252 @@
   ideal      qideal; /* extension to the ring structure: qring */
 
-  unsigned long bitmask;
 
   int      *VarOffset;
@@ -283 +282 @@
   short      ExpLSize; /* size of exponent vector in long */
   short      OrdSize; /* size of ord vector (in sro_ord) */
+  short      BitsPerExp; /* number of bits per exponent */
+
+  short      ref; /* reference counter to the ring */
+
+  /* number of long vars in exp vector: 
+     long vars are those longs in the exponent vector which are 
+     occupied by variables, only */
+  short     VarL_Size;   
+  /* if >= 0, long vars in exp vector are consecutive and start there
+     if <  0, long vars in exp vector are not consecutive */
+  short     VarL_LowIndex;
+  /* array of size VarL_Size, 
+     VarL_Offset[i] gets i-th long var in exp vector */
+  int*      VarL_Offset;
+
+  /* mask for getting single exponents */
+  unsigned long bitmask;
+  /* mask used for divisiblity tests */
+  unsigned long divmask;
 
   p_Procs_s* p_Procs;
-  short      ref; /* reference counter to the ring */
 };