Changeset 340ab8 in git for kernel

Timestamp:

Oct 5, 2016, 7:44:03 PM (8 years ago)

Author:

Sharwan Tiwari <stiwari@…>

Branches:

(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')

Children:

7e48fd30f913149068233c11121d6395b3cf6fa6

Parents:

637aa3957f20a1c1863a9ca0de302c346df62fd7

git-author:

Sharwan Tiwari <stiwari@gmail.com>2016-10-05 19:44:03+02:00

git-committer:

Sharwan Tiwari <stiwari@gmail.com>2016-10-11 17:40:28+02:00

Message:

nchilb

File:

• kernel/combinatorics/hilb.cc (modified) (19 diffs)

kernel/combinatorics/hilb.cc

@@ -1401 +1401 @@
 
 
-static void idInsertMonPos_nchilb(ideal I, poly p, int pos)
-{  
+static void idInsertMonomials(ideal I, poly p)
+{ 
   /*
-   * poly p != null, and there is no null poly in I
-   * adds p at the position 'pos' in I 
-   * without making  copy of it
-   */
-     
-  pEnlargeSet(&(I->m), IDELEMS(I), IDELEMS(I)+1);
-  IDELEMS(I) += 1;
-  for(int i = IDELEMS(I) - 1; i > pos; i--)
-  {
-    I->m[i] = I->m[i-1];
-  }
-  I->m[pos] = p;
-}
-
-static void idInsertMonAtEnd_nchilb(ideal I, poly p)
-{  
-  /* 
-   * poly p !=null, and there is no null poly in I
+   * adds monomial in I and if required,
+   * enlarges the size of poly-set by 16
+   * does not make copy of  p
    */
     
-  pEnlargeSet(&(I->m), IDELEMS(I), IDELEMS(I)+1);
-  IDELEMS(I) += 1;
-  I->m[IDELEMS(I)-1] = p;
-}
-
-static void idInsertMonomials(ideal h1, poly h2)
-{ 
-  /*
-   * adds monomial in h1 and if required,
-   * enlarges the size of poly-set by 16
-   * does not make copy of  h2
-   */
-    
-  if(h2 == NULL)
+  if(I == NULL)
   {
     return;
   }
 
-  int j = IDELEMS(h1) - 1;
-  while ((j >= 0) && (h1->m[j] == NULL))
+  int j = IDELEMS(I) - 1;
+  while ((j >= 0) && (I->m[j] == NULL))
   {
     j--;
   }
   j++;
-  if (j == IDELEMS(h1))
-  {
-    pEnlargeSet(&(h1->m), IDELEMS(h1), 16);
-    IDELEMS(h1) +=16;
-  }
-  h1->m[j] = h2;
+  if (j == IDELEMS(I))
+  {
+    pEnlargeSet(&(I->m), IDELEMS(I), 16);
+    IDELEMS(I) +=16;
+  }
+  I->m[j] = p;
 }
 
@@ -1470 +1442 @@
 */
 
-static ideal PlacePolyInSortedMonoIdeal_TotalDegOrder(ideal I, poly p)
-{
-  /*
-   * order of the currRing should be a total degree
-   * poly p != null
-   * inserts poly p in the I, s.t.
-   * the polys of I are in ascending order
-   * does not make copy of p2
-   */
-  int i;
-  if((I == NULL) || (idIs0(I)))
-  {
-    ideal res = idInit(1,1);
-    res->m[0] = p;
-    return(res);
-  }
-  idSkipZeroes(I);
-
-  for(i = 0; (i<IDELEMS(I)) && (p_Totaldegree(I->m[i], currRing) <= p_Totaldegree(p, currRing)); i++)
-  {
-    if(p_DivisibleBy(I->m[i], p, currRing))
-    {
-      pDelete(&p);//this poly is not required anymore
-      return(I);
-    }
-  }
-  for(i = IDELEMS(I)-1; (i>=0) && (p_Totaldegree(I->m[i],currRing) >= p_Totaldegree(p,currRing)); i--)
-  {
-    if(p_DivisibleBy(p, I->m[i], currRing))
-    {
-      I->m[i] = NULL;
-    }
-  }
-  if(idIs0(I))
-  {
-    I->m[0] = p;
-    return(I);
-  }
-
-  idSkipZeroes(I);
-  //when all polynomials of I have same degree
-  if(p_Totaldegree(I->m[0], currRing) == p_Totaldegree(I->m[IDELEMS(I)-1], currRing))
-  {
-    if(p_Totaldegree(p, currRing) < p_Totaldegree(I->m[0], currRing))
-    {
-      idInsertMonPos_nchilb(I, p, 0);
-      return(I);
-    }
-    if(p_Totaldegree(p, currRing) >= p_Totaldegree(I->m[IDELEMS(I)-1], currRing))
-    {
-
-      for(i = IDELEMS(I)-1; i>=0; i--)
-      {
-        if(p_LmCmp(I->m[i], p, currRing) == 1)
-        {        
-          if(i==0)
-          {
-            idInsertMonPos_nchilb(I, p, 0);
-            return(I);
-          }
-        }
-        else
-        {
-          if(i == (IDELEMS(I)-1))
-          {
-           idInsertMonAtEnd_nchilb(I, p);
-            return(I);
-          }
-          idInsertMonPos_nchilb(I, p, i+1);
-          return(I);
-        }
-      }
-    }
-  }
-  if(p_Totaldegree(p, currRing) <= p_Totaldegree(I->m[0], currRing))
-  {
-   
-    for(i = 0; i < IDELEMS(I); i++)
-    {
-      if(p_LmCmp(I->m[i], p, currRing) == 1)
-      {
-        idInsertMonPos_nchilb(I, p, i);
-        return(I);
-      }
-    }
-    idInsertMonAtEnd_nchilb(I, p);
-    return(I);
-  }
-  if(p_Totaldegree(p, currRing) >= p_Totaldegree(I->m[IDELEMS(I)-1], currRing))
-  {
-    for(i = IDELEMS(I)-1; i >= 0; i--)
-    {
-      if(p_LmCmp(I->m[i], p, currRing) != 1)
-      {
-        if(i == IDELEMS(I)-1)
-        {
-          idInsertMonAtEnd_nchilb(I, p);
-          return(I);
-        }
-        else
-        {
-          idInsertMonPos_nchilb(I, p, i+1);
-          return(I);
-        }
-
-      }
-    }
-
-  }
-  
-  for(i = IDELEMS(I)-2; ;)
-  {
-    if(p_Totaldegree(p, currRing) == p_Totaldegree(I->m[i], currRing))
-    {
-      for(int j = i; j >= 0; j--)
-      {
-        if(p_LmCmp(I->m[j], p, currRing) != 1)
-        {
-          idInsertMonPos_nchilb(I, p, j+1);
-          return(I);
-        }
-      }
-    }
-    if(p_Totaldegree(p, currRing) > p_Totaldegree(I->m[i], currRing))
-    {
-      idInsertMonPos_nchilb(I, p, i+1);
-      return(I);
-    }
-    i--;
-  }
-}
-
-static ideal SortMonoId_TotalDegOrder(ideal I)
-{
-  /*
-   * sorts the ideal in ascending order
-   * order must be a total degree 
-   */
-  if(idIs0(I))
-  {
-    return(I);
-  }
-  int i;
-  ideal result;
-  idSkipZeroes(I);
-  result = idInit(1, 1);
-  for(i = 0; i <= IDELEMS(I)-1; i++)
-  {
-     result = PlacePolyInSortedMonoIdeal_TotalDegOrder(result, I->m[i]);
-  }
-  idSkipZeroes(result);
-  return(result);
-}
-
 static int isMonoIdBasesSame(ideal J, ideal Ob)
 {
@@ -1632 +1450 @@
    */
   int i, j, s, cnt;
-  s=1;
+  s = 1;
   int JCount = IDELEMS(J);
   int ObCount = IDELEMS(Ob);
@@ -1669 +1487 @@
   }
   
-  int count=0;
+  int count = 0;
   for(int i = 0; i < IDELEMS(I); i++)
   {
@@ -1676 +1494 @@
       return (count);
     }
-    count=count+1;
+    count = count + 1;
   }
   
@@ -1826 +1644 @@
 }
 
-static ideal  MinimalMonomialBasis(ideal I)
+static int monCompare(const void  *m, const void *n)
+{
+  /* compares monomials */
+
+  return (p_Compare(*(poly*)m, *(poly*)n, currRing));
+}
+
+static void sortMonoIdeal_totalDegOrder(ideal I)
+{
+  /*
+   * sorts the monomial ideal in ascending order
+   * order must be a total degree
+   */
+
+  qsort(I->m, IDELEMS(I), sizeof(poly), monCompare);
+
+}
+
+
+
+static ideal  minimalMonomialsGenSet(ideal I)
 {
   /*
@@ -1832 +1670 @@
    * can be generated by others in I
    */
-  I = SortMonoId_TotalDegOrder(I);
+  //first sort monomials of the ideal
+  
+  idSkipZeroes(I);
+
+  sortMonoIdeal_totalDegOrder(I);
+  
   ideal J = idInit(1, 1);
   int i, k;
   int count = 0;
   int ICount = IDELEMS(I);
-  
-  std::vector<BOOLEAN> checks(ICount, TRUE);
-  for(k = 0; k < ICount; k++)
-  {
-    if(checks[k])
-    {
-      idInsertMonomials(J, I->m[k]);
-      count++;
-      for(i = count; i < ICount; i++)
-      {
-        if(checks[i])
-        {
-          if( pLmDivisibleByNoComp(J->m[count-1], I->m[i]))
-          {
-            checks[i]=FALSE;
-            pDelete(&(I->m[i]));
-          }
-        }
+
+  for(k = ICount - 1; k >=1; k--)
+  {
+    for(i = 0; i < k; i++)
+    {
+  
+      if(p_LmDivisibleBy(I->m[i], I->m[k], currRing))
+      { 
+        pDelete(&(I->m[k]));//this is not req.
+        break;
       }
     }
   }
-  idSkipZeroes(J);
-  return(J);
+  
+  idSkipZeroes(I);
+  return(I);
 }
 
@@ -1884 +1720 @@
     }
   }
-  p_SetExpV(smon, s, r);
+  
+  p_SetExpV(smon, s, currRing);
   omFree(e);
   omFree(s);
-  p_SetComp(smon, p_GetComp(p, r), r);
+  
+  p_SetComp(smon, p_GetComp(p, currRing), currRing);
+  p_Setm(smon, currRing);
+  
   return(smon);
 }
@@ -1914 +1754 @@
     s[j] = e[j];
   }
-  p_SetExpV(dw, s, r);//new exponents
+  
+  p_SetExpV(dw, s, currRing);//new exponents
   omFree(e);
   omFree(s);
-  p_SetComp(dw, p_GetComp(w, r), r);
+  
+  p_SetComp(dw, p_GetComp(w, currRing), currRing);
+  p_Setm(dw, currRing);
+  
   return(dw);
 }
@@ -1933 +1777 @@
    */
   int i;
-  //ring r=currRing;
   poly smon, dw;
   poly qmonp = NULL;
   bool del;
+  
   for(i = 0;i <= d - 1; i++)
   {
@@ -1942 +1786 @@
     smon = shiftInMon(p, i, lV, currRing);
     del = TRUE;
-    if(pLmDivisibleByNoComp(smon, w))
+    
+    if(pLmDivisibleBy(smon, w))
     {
       flag = TRUE;
@@ -1950 +1795 @@
       pDelete(&smon);
       
-      //make Jwi =1;
+      //delete all monomials of Jwi
+      //and make Jwi =1
+      
       for(int j = 0;j < IDELEMS(Jwi); j++)
       {
         pDelete(&Jwi->m[j]);
       }
-      
-      pEnlargeSet(&(Jwi->m), IDELEMS(Jwi), -IDELEMS(Jwi)+1);
-      IDELEMS(Jwi)=1;
-      Jwi->m[0]=p_One(currRing);
+
+      idInsertMonomials(Jwi, p_One(currRing)); 
       break;
     }
 
-    if(pLmDivisibleByNoComp(dw, smon))
+    if(pLmDivisibleBy(dw, smon))
     { 
       del = FALSE; 
       qmonp = p_Divide(smon, dw, currRing);
-
-      pEnlargeSet(&(Jwi->m), IDELEMS(Jwi), 1);
-      IDELEMS(Jwi) += 1;
-      Jwi->m[IDELEMS(Jwi)-1] = shiftInMon(qmonp, -d, lV, currRing);
+      idInsertMonomials(Jwi, shiftInMon(qmonp, -d, lV, currRing));
       
-      pDelete(&qmonp);//shiftInMon(qmonp, -d, lV, currRing):returns a new poly,
-                     // qmonp remains unchanged
+      //shiftInMon(qmonp, -d, lV, currRing):returns a new poly,
+      //qmonp remains unchanged, delete it
+      pDelete(&qmonp);
       pDelete(&dw);
       pDelete(&smon);
     }
-    //in case both if were false, delete dw and smon
+    //in case both if are false, delete dw and smon
     if(del)
     {
@@ -2012 +1855 @@
     }
   }
-  Jwi = MinimalMonomialBasis(Jwi);
+  
+  Jwi = minimalMonomialsGenSet(Jwi);
   return(Jwi);
 }
+
 
 void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE )
@@ -2029 +1874 @@
 
   int trInd;
-  S = SortMonoId_TotalDegOrder(S);
-
-
+  S = minimalMonomialsGenSet(S);
+  
   int (*POS)(ideal, poly, std::vector<ideal>, std::vector<poly>, int);
   if(IG_CASE)
@@ -2056 +1900 @@
 
   std::vector<int> C;
-
-  ring r = currRing;
 
   int ds, is, ps, sz;
@@ -2094 +1936 @@
 
       wi = pCopy(w);
-      p_SetExp(wi, (ds*lV)+is, 1, r);
-      p_Setm(wi, r);
+      p_SetExp(wi, (ds*lV)+is, 1, currRing);
+      p_Setm(wi, currRing);
 
       Jwi = NULL;
@@ -2155 +1997 @@
   //printf("\nlength of the Orbit==%ld\n", C.size());
 #endif
-
+  ring r = currRing;
   char** tt=(char**)omalloc(sizeof(char*));
   tt[0] = omStrDup("t");
@@ -2179 +2021 @@
       if(mat[rowCount][colCount] != 0)
       {
-        mR->m[lO*rowCount + colCount] = p_ISet(mat[rowCount][colCount], R);
-        p_SetCoeff(mR->m[lO*rowCount+colCount], n_Mult(pGetCoeff(mR->m[lO*rowCount+colCount]),n_Param(1, R->cf), R->cf), R);
+       MATELEM(mR, rowCount + 1, colCount + 1) = p_ISet(mat[rowCount][colCount], R);
+        p_SetCoeff(MATELEM(mR, rowCount + 1, colCount + 1), n_Mult(pGetCoeff(mR->m[lO*rowCount+colCount]),n_Param(1, R->cf), R->cf), R);
       }
     }
@@ -2219 +2061 @@
   omFree(xx[0]);
   omFree(xx);
+  rKill(R);
   rChangeCurrRing(r);
 }