Changeset 6d09f28 in git

Timestamp:

Aug 14, 2006, 7:08:37 PM (18 years ago)

Author:

Oliver Wienand <wienand@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

d39ce40a2e8c985712a746d769cc723fe80f080f

Parents:

2a5c2f4b88044ea4f445da171e3d11679bd20548

Message:

introduced def HAVE_VANGB
--> compute GB in the ring of function rather than polynomials

kstd2.cc:
* moved ind2 and ind_fact_2 to kutil.h/cc
* bba: enterpairs changes size of T in VANGB case

kutil.h/cc:
* ind2/ind_fact_2 see above
* kCreateZeroPoly: create a zero polynomial with given exponent
* createG0: creates the groebner basis of the vanishing polynomials


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

Location:

kernel

Files:

• kstd2.cc (modified) (9 diffs)
• kutil.cc (modified) (6 diffs)
• kutil.h (modified) (2 diffs)

kernel/kstd2.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: kstd2.cc,v 1.17 2006-06-07 18:44:23 wienand Exp $ */
+/* $Id: kstd2.cc,v 1.18 2006-08-14 17:08:36 wienand Exp $ */
 /*
 *  ABSTRACT -  Kernel: alg. of Buchberger
@@ -186 +186 @@
 }
 
-long ind2(long arg)
-{
-  long ind = 0;
-  if (arg <= 0) return 0;
-  while (arg%2 == 0)
-  {
-    arg = arg / 2;
-    ind++;
-  }
-  return ind;
-}
-
-long ind_fact_2(long arg)
-{
-  long ind = 0;
-  if (arg <= 0) return 0;
-  if (arg%2 == 1) { arg--; }
-  while (arg > 0)
-  {
-    ind += ind2(arg);
-    arg = arg - 2;
-  }
-  return ind;
-}
-
 poly kFindZeroPoly(poly input_p, ring leadRing, ring tailRing)
 {
@@ -225 +200 @@
 
   long k = 1;
-  // of interest is only k_ind2, special routine for improvement ... TOTO OLIVER
+  // of interest is only k_ind2, special routine for improvement ... TODO OLIVER
   for (int i = 1; i <= leadRing->N; i++)
   {
@@ -282 +257 @@
     return tmp2;
   }
-  long alpha_k = twoPow(leadRing->ch - k_ind2);
+/*  long alpha_k = twoPow(leadRing->ch - k_ind2);
   if (1 == 0 && alpha_k <= a) {  // Temporarly disabled, reducing coefficients not compatible with std TODO Oliver
     zeroPoly = p_ISet((a / alpha_k)*alpha_k, tailRing);
@@ -308 +283 @@
     pNext(tmp2) = zeroPoly;
     return tmp2;
-  }
+  } */
   return NULL;
 }
@@ -343 +318 @@
   loop
   {
-    zeroPoly = NULL; //kFindDivisibleByZeroPoly(h);
+#ifdef HAVE_VANGB
+    zeroPoly = kFindDivisibleByZeroPoly(h);
     if (zeroPoly != NULL)
     {
@@ -366 +342 @@
     }
     else
+#endif
     {
       j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, h);
@@ -994 +971 @@
       enterT(strat->P, strat);
 #ifdef HAVE_RING2TOM
+#ifdef HAVE_VANGB
+      int at_R = strat->tl;
+#endif
       if (currRing->cring == 1)
         superenterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
@@ -1000 +980 @@
         enterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
       // posInS only depends on the leading term
+#ifdef HAVE_VANGB
+      strat->enterS(strat->P, pos, strat, at_R);
+#else
       strat->enterS(strat->P, pos, strat, strat->tl);
+#endif
       if (hilb!=NULL) khCheck(Q,w,hilb,hilbeledeg,hilbcount,strat);
 //      Print("[%d]",hilbeledeg);

kernel/kutil.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: kutil.cc,v 1.29 2006-06-20 21:44:18 wienand Exp $ */
+/* $Id: kutil.cc,v 1.30 2006-08-14 17:08:36 wienand Exp $ */
 /*
 * ABSTRACT: kernel: utils for kStd
@@ -23 +23 @@
 #undef KDEBUG
 #define KDEBUG 2
+#endif
+
+#ifdef HAVE_RING2TOM
+#include "ideals.h"
 #endif
 
@@ -1018 +1022 @@
   (*Lp).length = 0;
 }
-
+         
 void initEcartPairMora (LObject* Lp,poly f,poly g,int ecartF,int ecartG)
 {
@@ -1156 +1160 @@
   else
   {
-    Lp.p = ksCreateShortSpoly(strat->S[i],p, strat->tailRing);
+    Lp.p = ksCreateShortSpoly(strat->S[i], p, strat->tailRing);
   }
   if (Lp.p == NULL)
@@ -1942 +1946 @@
 }
 
+long ind2(long arg)
+{
+  long ind = 0;
+  if (arg <= 0) return 0;
+  while (arg%2 == 0)
+  {
+    arg = arg / 2;
+    ind++;
+  }
+  return ind;
+}
+
+long ind_fact_2(long arg)
+{
+  long ind = 0;
+  if (arg <= 0) return 0;
+  if (arg%2 == 1) { arg--; }
+  while (arg > 0)
+  {
+    ind += ind2(arg);
+    arg = arg - 2;
+  }
+  return ind;
+}
+
+/*2
+* put the pair (p, f) in B and f in T
+*/
+void enterOneZeroPairRing (poly f, poly t_p, poly p, int ecart, kStrategy strat, int atR = -1)
+{
+  int      l,j,compare,compareCoeff;
+  LObject  Lp;
+
+  if (strat->interred_flag) return;
+#ifdef KDEBUG
+  Lp.ecart=0; Lp.length=0;
+#endif
+  /*- computes the lcm(s[i],p) -*/
+  Lp.lcm = pInit();
+
+  pLcm(p,f,Lp.lcm);
+  pSetm(Lp.lcm);
+  pSetCoeff(Lp.lcm, nLcm(pGetCoeff(p), pGetCoeff(f), currRing));
+  assume(!strat->sugarCrit);
+  assume(!strat->fromT);
+  /*
+  *the set B collects the pairs of type (S[j],p)
+  *suppose (r,p) is in B and (s,p) is the new pair and lcm(s,p) != lcm(r,p)
+  *if the leading term of s devides lcm(r,p) then (r,p) will be canceled
+  *if the leading term of r devides lcm(s,p) then (s,p) will not enter B
+  */
+  for(j = strat->Bl;j>=0;j--)
+  {
+    compare=pDivCompRing(strat->B[j].lcm,Lp.lcm);
+    compareCoeff = nComp((long) pGetCoeff(strat->B[j].lcm), (long) pGetCoeff(Lp.lcm));
+    if (compareCoeff == 0 || compare == compareCoeff)
+    {
+      if (compare == 1)
+      {
+        strat->c3++;
+        pLmFree(Lp.lcm);
+        return;
+      }
+      else
+      if (compare == -1)
+      {
+        deleteInL(strat->B,&strat->Bl,j,strat);
+        strat->c3++;
+      }
+    }
+    if (compare == pDivComp_EQUAL)
+    {
+      // Add hint for same LM and direction of LC (later) (TODO Oliver)
+      if (compareCoeff == 1)
+      {
+        strat->c3++;
+        pLmFree(Lp.lcm);
+        return;
+      }
+      else
+      if (compareCoeff == -1)
+      {
+        deleteInL(strat->B,&strat->Bl,j,strat);
+        strat->c3++;
+      }
+    }
+  }
+  /*
+  *the pair (S[i],p) enters B if the spoly != 0
+  */
+  /*-  compute the short s-polynomial -*/
+  if ((f==NULL) || (p==NULL)) return;
+  pNorm(p);
+  {
+    Lp.p = ksCreateShortSpoly(f, p, strat->tailRing);
+  }
+  if (Lp.p == NULL) //deactivated, as we are adding pairs with zeropoly and not from S
+  {
+    /*- the case that the s-poly is 0 -*/
+//    if (strat->pairtest==NULL) initPairtest(strat);
+//    strat->pairtest[i] = TRUE;/*- hint for spoly(S^[i],p)=0 -*/
+//    strat->pairtest[strat->sl+1] = TRUE;
+    /*hint for spoly(S[i],p) == 0 for some i,0 <= i <= sl*/
+    /*
+    *suppose we have (s,r),(r,p),(s,p) and spoly(s,p) == 0 and (r,p) is
+    *still in B (i.e. lcm(r,p) == lcm(s,p) or the leading term of s does not
+    *devide lcm(r,p)). In the last case (s,r) can be canceled if the leading
+    *term of p devides the lcm(s,r)
+    *(this canceling should be done here because
+    *the case lcm(s,p) == lcm(s,r) is not covered in chainCrit)
+    *the first case is handeled in chainCrit
+    */
+    if (Lp.lcm!=NULL) pLmFree(Lp.lcm);
+  }
+  else
+  {
+    /*- the pair (S[i],p) enters B -*/
+    Lp.p1 = f;
+    Lp.p2 = p;
+
+    pNext(Lp.p) = strat->tail;
+
+    LObject tmp_h(f, currRing, strat->tailRing);
+    tmp_h.SetShortExpVector();
+    strat->initEcart(&tmp_h);
+    tmp_h.sev = pGetShortExpVector(tmp_h.p);
+    tmp_h.SetpFDeg();
+    tmp_h.t_p = t_p;
+
+    enterT(tmp_h, strat, strat->tl + 1);
+
+    if (atR >= 0)
+    {
+      Lp.i_r2 = atR;
+      Lp.i_r1 = strat->tl;
+    }
+
+    strat->initEcartPair(&Lp,f,p,0/*strat->ecartS[i]*/,ecart);     // Attention: TODO: break ecart
+    l = strat->posInL(strat->B,strat->Bl,&Lp,strat);
+    // enterL(&strat->B, &strat->Bl, &strat->Bmax, Lp, l);
+  }
+}
+
+/* Helper for kCreateZeroPoly
+ * enumerating the exponents
+ */
+
+int nextZeroSimplexExponent (long exp[], long ind[], long cexp[], long cind[], long* cabsind, long bound, long N)
+/* gives the next exponent from the set H_1 */
+{
+  if (*cabsind < bound)
+  {
+    cexp[1] += 2;
+    long add = ind2(cexp[1]);
+    cind[1] += add;
+    *cabsind += add;
+  }
+  else
+  {
+    // cabsind >= habsind
+    if (N == 1) return 0;
+    int i = 1;
+    while (exp[i] == cexp[i] && i <= N) i++;
+    cexp[i] = exp[i];
+    *cabsind -= cind[i];
+    cind[i] = ind[i];
+    *cabsind += cind[i];
+    // Print("in: %d\n", *cabsind);
+    i += 1;
+    if (i > N) return 0;
+    cexp[i] += 2;
+    long add = ind2(cexp[i]);
+    cind[i] += add;
+    *cabsind += add;
+  }
+  return 1;
+}
+
+/*
+ * Creates the zero Polynomial on position exp
+ * long exp[] : exponent of leading term
+ * cabsind    : total 2-ind of exp (if -1 will be computed)
+ * poly* t_p  : will hold the LT in tailRing
+ * leadRing   : ring for the LT
+ * tailRing   : ring for the tail
+ */
+
+poly kCreateZeroPoly(long exp[], long cabsind, poly* t_p, ring leadRing, ring tailRing)
+{
+
+  poly zeroPoly = NULL;
+
+  number tmp1;
+  poly tmp2, tmp3;
+
+  if (cabsind == -1)
+  {
+    cabsind = 0;
+    for (int i = 1; i <= leadRing->N; i++)
+    {
+      cabsind += ind_fact_2(exp[i]);
+    }
+//    Print("cabsind: %d\n", cabsind);
+  }
+  if (cabsind < leadRing->ch)
+  {
+    zeroPoly = p_ISet(twoPow(leadRing->ch - cabsind), tailRing);
+  }
+  else
+  {
+    zeroPoly = p_ISet(1, tailRing);
+  }
+  for (int i = 1; i <= leadRing->N; i++)
+  {
+    for (long j = 1; j <= exp[i]; j++)
+    {
+      tmp1 = nInit(j);
+      tmp2 = p_ISet(1, tailRing);
+      p_SetExp(tmp2, i, 1, tailRing);
+      p_Setm(tmp2, tailRing);
+      if (nIsZero(tmp1))
+      { // should nowbe obsolet, test ! TODO OLIVER
+        zeroPoly = p_Mult_q(zeroPoly, tmp2, tailRing);
+      }
+      else
+      {
+        tmp3 = p_NSet(nCopy(tmp1), tailRing);
+        zeroPoly = p_Mult_q(zeroPoly, p_Add_q(tmp3, tmp2, tailRing), tailRing);
+      }
+    }
+  }
+  tmp2 = p_NSet(nCopy(pGetCoeff(zeroPoly)), leadRing);
+  for (int i = 1; i <= leadRing->N; i++) {
+    pSetExp(tmp2, i, p_GetExp(zeroPoly, i, tailRing));
+  }
+  p_Setm(tmp2, leadRing);
+  *t_p = zeroPoly;
+  zeroPoly = pNext(zeroPoly);
+  pNext(*t_p) = NULL;
+  pNext(tmp2) = zeroPoly;
+  return tmp2;
+}
+
+// #define OLI_DEBUG
+
+/*
+ * Generate the s-polynomial for the virtual set of zero-polynomials
+ */
+
+void initenterzeropairsRing (poly p, int ecart, kStrategy strat, int atR)
+{
+  // Initialize
+  long exp[50];            // The exponent of \hat{X} (basepoint)
+  long cexp[50];           // The current exponent for iterating over all
+  long ind[50];            // The power of 2 in the i-th component of exp
+  long cind[50];           // analog for cexp
+  long mult[50];           // How to multiply the elements of G
+  long cabsind = 0;        // The abs. index of cexp, i.e. the sum of cind
+  long habsind = 0;        // The abs. index of the coefficient of h
+  for (int i = 1; i <= currRing->N; i++)
+  {
+    exp[i] = p_GetExp(p, i, currRing);
+    if (exp[i] & 1 != 0)
+    {
+      exp[i] = exp[i] - 1;
+      mult[i] = 1;
+    }
+    cexp[i] = exp[i];
+    ind[i] = ind_fact_2(exp[i]);
+    cabsind += ind[i];
+    cind[i] = ind[i];
+  }
+  habsind = ind2((long) p_GetCoeff(p, currRing));
+  long bound = currRing->ch - habsind;
+#ifdef OLI_DEBUG
+  PrintS("-------------\npoly  :");
+  wrp(p);
+  Print("\nexp   : (%d, %d)\n", exp[1] + mult[1], exp[2] + mult[1]);
+  Print("cexp  : (%d, %d)\n", cexp[1], cexp[2]);
+  Print("cind  : (%d, %d)\n", cind[1], cind[2]);
+  Print("bound : %d\n", bound);
+  Print("cind  : %d\n", cabsind);
+#endif
+  if (cabsind == 0)
+  {
+    if (!(nextZeroSimplexExponent(exp, ind, cexp, cind, &cabsind, bound, currRing->N)))
+    {
+      return;
+    }
+  }
+  // Now the whole simplex
+  do
+  {
+    // Build s-polynomial
+    // 2**ind-def * mult * g - exp-def * h
+    poly t_p;
+    poly zeroPoly = kCreateZeroPoly(cexp, cabsind, &t_p, currRing, strat->tailRing);
+#ifdef OLI_DEBUG
+    Print("%d, (%d, %d), ind = (%d, %d)\n", cabsind, cexp[1], cexp[2], cind[1], cind[2]);
+    Print("zPoly : ");
+    wrp(zeroPoly);
+    Print("\n");
+#endif
+    enterOneZeroPairRing(zeroPoly, t_p, p, ecart, strat, atR);
+  }
+  while ( nextZeroSimplexExponent(exp, ind, cexp, cind, &cabsind, bound, currRing->N) );
+}
+
+/* 
+ * Create the Groebner basis of the vanishing polynomials.
+ */
+
+ideal createG0()
+{
+  // Initialize
+  long exp[50];            // The exponent of \hat{X} (basepoint)
+  long cexp[50];           // The current exponent for iterating over all
+  long ind[50];            // The power of 2 in the i-th component of exp
+  long cind[50];           // analog for cexp
+  long mult[50];           // How to multiply the elements of G
+  long cabsind = 0;        // The abs. index of cexp, i.e. the sum of cind
+  long habsind = 0;        // The abs. index of the coefficient of h
+  for (int i = 1; i <= currRing->N; i++)
+  {
+    exp[i] = 0;
+    cexp[i] = exp[i];
+    ind[i] = 0;
+    cind[i] = ind[i];
+  }
+  long bound = currRing->ch;
+#ifdef OLI_DEBUG
+  PrintS("-------------\npoly  :");
+  wrp(p);
+  Print("\nexp   : (%d, %d)\n", exp[1] + mult[1], exp[2] + mult[1]);
+  Print("cexp  : (%d, %d)\n", cexp[1], cexp[2]);
+  Print("cind  : (%d, %d)\n", cind[1], cind[2]);
+  Print("bound : %d\n", bound);
+  Print("cind  : %d\n", cabsind);
+#endif
+  if (cabsind == 0)
+  {
+    if (!(nextZeroSimplexExponent(exp, ind, cexp, cind, &cabsind, bound, currRing->N)))
+    {
+      return idInit(1, 1);
+    }
+  }
+  ideal G0 = idInit(1, 1);
+  // Now the whole simplex
+  do
+  {
+    // Build s-polynomial
+    // 2**ind-def * mult * g - exp-def * h
+    poly t_p;
+    poly zeroPoly = kCreateZeroPoly(cexp, cabsind, &t_p, currRing, currRing);
+#ifdef OLI_DEBUG
+    Print("%d, (%d, %d), ind = (%d, %d)\n", cabsind, cexp[1], cexp[2], cind[1], cind[2]);
+    Print("zPoly : ");
+    wrp(zeroPoly);
+    Print("\n");
+#endif
+    // Add to ideal
+    pEnlargeSet(&(G0->m), IDELEMS(G0), 1);
+    IDELEMS(G0) += 1;
+    G0->m[IDELEMS(G0) - 1] = zeroPoly;
+  }
+  while ( nextZeroSimplexExponent(exp, ind, cexp, cind, &cabsind, bound, currRing->N) );
+  idSkipZeroes(G0);
+  return G0;
+}
+
 /*2
 *(s[0],h),...,(s[k],h) will be put to the pairset L
@@ -1990 +2364 @@
       }
     }
+
+#ifdef HAVE_VANGB
+    initenterzeropairsRing(h, ecart, strat, atR);
+#endif
 
     if (new_pair) chainCritRing(h,ecart,strat);

kernel/kutil.h

@@ -4 +4 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: kutil.h,v 1.17 2006-06-12 17:40:10 wienand Exp $ */
+/* $Id: kutil.h,v 1.18 2006-08-14 17:08:37 wienand Exp $ */
 /*
 * ABSTRACT: kernel: utils for kStd
@@ -397 +397 @@
 #ifdef HAVE_RING2TOM
 int redRing2toM (LObject* h,kStrategy strat);
+long ind2(long arg);
+long ind_fact_2(long arg);
 long twoPow(long arg);
 void enterExtendedSpoly(poly h,kStrategy strat);
 void superenterpairs (poly h,int k,int ecart,int pos,kStrategy strat, int atR = -1);
+poly kCreateZeroPoly(long exp[], long cabsind, poly* t_p, ring leadRing, ring tailRing);
+ideal createG0();
 #endif
 int redLazy (LObject* h,kStrategy strat);