Changeset d770e6 in git for Singular

Timestamp:

Oct 29, 2013, 1:54:17 PM (11 years ago)

Author:

Yue Ren <ren@…>

Branches:

(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '4bd32dfef92ec9f5ed8dceee82d14318ae147107')

Children:

7395dd26103d1edd6fc4505722010d2bba012d96

Parents:

a73fae59ddd878f87c48e12088f6b0046549aaf7

git-author:

Yue Ren <ren@mathematik.uni-kl.de>2013-10-29 13:54:17+01:00

git-committer:

Yue Ren <ren@mathematik.uni-kl.de>2015-02-06 13:47:01+01:00

Message:

new: initialReduction

Location:

Singular/dyn_modules/gfanlib

Files:

• bbcone.cc (modified) (1 diff)
• bbcone.h (modified) (1 diff)
• tropical.cc (modified) (5 diffs)

Singular/dyn_modules/gfanlib/bbcone.cc

@@ -113 +113 @@
     zv[j]=i[j+1];
   return zv;
+}
+
+int* ZVectorToIntStar(const gfan::ZVector &v, bool &overflow)
+{
+  int* w = (int*) omAlloc(v.size()*sizeof(int));
+  for (unsigned i=0; i<v.size(); i++)
+  {
+    if (!v[i].fitsInInt())
+    {
+      omFree(w);
+      WerrorS("intoverflow converting gfan:ZVector to int*");
+      overflow = true;
+      return NULL;
+    }
+    w[i]=v[i].toInt();
+  }
+  return w;
 }
 

Singular/dyn_modules/gfanlib/bbcone.h

@@ -27 +27 @@
 
 gfan::ZVector intStar2ZVector(const int d, const int* i);
+int* ZVectorToIntStar(const gfan::ZVector &v, bool &overflow);
 char* toString(gfan::ZMatrix const &m);
 std::string toString(const gfan::ZCone* const c);

Singular/dyn_modules/gfanlib/tropical.cc

@@ -1 +1 @@
 #include <kernel/polys.h>
+#include <kernel/kstd1.h>
 #include <libpolys/coeffs/longrat.h>
+#include <libpolys/polys/clapsing.h>
 #include <bbcone.h>
 
@@ -154 +156 @@
 
 
+gfan::ZCone* maximalGroebnerCone(const ring &r, const ideal &I)
+{
+  int n = rVar(r);
+  poly g = NULL;
+  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
+  int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
+  gfan::ZVector leadexpw = gfan::ZVector(n);
+  gfan::ZVector tailexpw = gfan::ZVector(n);
+  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
+  for (int i=0; i<IDELEMS(I); i++)
+  {
+    g = (poly) I->m[i]; pGetExpV(g,leadexpv);
+    leadexpw = intStar2ZVector(n, leadexpv);
+    pIter(g);
+    while (g != NULL)
+    {
+      pGetExpV(g,tailexpv);
+      tailexpw = intStar2ZVector(n, tailexpv);
+      inequalities.appendRow(leadexpw-tailexpw);
+      pIter(g);
+    }
+  }
+  omFreeSize(leadexpv,(n+1)*sizeof(int));
+  omFreeSize(tailexpv,(n+1)*sizeof(int));
+  return new gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
+}
+
+
 BOOLEAN maximalGroebnerCone(leftv res, leftv args)
 {
@@ -162 +192 @@
     if (v == NULL)
     {
-      int n = currRing->N;
       ideal I = (ideal) u->Data();
-      poly g = NULL;
-      int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
-      int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
-      gfan::ZVector leadexpw = gfan::ZVector(n);
-      gfan::ZVector tailexpw = gfan::ZVector(n);
-      gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
-      for (int i=0; i<IDELEMS(I); i++)
-      {
-        g = (poly) I->m[i]; pGetExpV(g,leadexpv);
-        leadexpw = intStar2ZVector(n, leadexpv);
-        pIter(g);
-        while (g != NULL)
-        {
-          pGetExpV(g,tailexpv);
-          tailexpw = intStar2ZVector(n, tailexpv);
-          inequalities.appendRow(leadexpw-tailexpw);
-          pIter(g);
-        }
-      }
-      gfan::ZCone* gCone = new gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
-      omFreeSize(leadexpv,(n+1)*sizeof(int));
-      omFreeSize(tailexpv,(n+1)*sizeof(int));
-
       res->rtyp = coneID;
-      res->data = (void*) gCone;
+      res->data = (void*) maximalGroebnerCone(currRing, I);
       return FALSE;
     }
@@ -195 +201 @@
   return TRUE;
 }
+
+/***
+ * If 1, replaces all occuring t with prime p,
+ * where theoretically feasable.
+ * Also computes GCD over integers than
+ * over univariate polynomials
+ **/
+#define T_TO_P 0
+
+/***
+ * Suppose I=g_0,...,g_{n-1}.
+ * This function uses bubble sort to sorts g_1,...,g_{n-1}
+ * such that lm(g_1)>...>lm(g_{n-1})
+ **/
+static inline void sortingLaterGenerators(ideal I)
+{
+  poly cache; int i; int n=IDELEMS(I); int newn;
+  do
+  {
+    newn=0;
+    for (i=2; i<n; i++)
+    {
+      if (pLmCmp(I->m[i-1],I->m[i])<0)
+      {
+        cache=I->m[i-1];
+        I->m[i-1]=I->m[i];
+        I->m[i]=cache;
+        newn = i;
+      }
+    }
+    n=newn;
+  } while(n);
+}
+
+
+/***
+ * replaces coefficients in g of lowest degree in t
+ * divisible by p with a suitable power of t
+ **/
+static bool preduce(poly g, const number p)
+{
+  // go along g and look for relevant terms.
+  // monomials in x which have been checked are tracked in done.
+  // because we assume the leading coefficient of g is 1,
+  // the leading term does not need to be considered.
+  poly done = p_LmInit(g,currRing);
+  p_SetExp(done,1,0,currRing); p_SetCoeff(done,n_Init(1,currRing->cf),currRing);
+  p_Setm(done,currRing);
+  poly doneEnd = done;
+  poly doneCache;
+  number dnumber; long d;
+  poly subst; number ppower;
+  while(pNext(g))
+  {
+    // check whether next term needs to be reduced:
+    // first, check whether monomial in x has come up yet
+    for (doneCache=done; doneCache; pIter(doneCache))
+    {
+      if (p_LmDivisibleBy(doneCache,pNext(g),currRing))
+        break;
+    }
+    if (!doneCache) // if for loop did not terminate prematurely,
+                    // then the monomial in x is new
+    {
+      // second, check whether coefficient is divisible by p
+      if (n_DivBy(p_GetCoeff(pNext(g),currRing->cf),p,currRing->cf))
+      {
+        // reduce the term with respect to p-t:
+        // divide coefficient by p, remove old term,
+        // add t multiple of old term
+        dnumber = n_Div(p_GetCoeff(pNext(g),currRing->cf),p,currRing->cf);
+        d = n_Int(dnumber,currRing->cf);
+        n_Delete(&dnumber,currRing->cf);
+        if (!d)
+        {
+          p_Delete(&done,currRing);
+          WerrorS("initialReduction: overflow in a t-exponent");
+          return true;
+        }
+        subst=p_LmInit(pNext(g),currRing);
+        if (d>0)
+        {
+          p_AddExp(subst,1,d,currRing);
+          p_SetCoeff(subst,n_Init(1,currRing->cf),currRing);
+        }
+        else
+        {
+          p_AddExp(subst,1,-d,currRing);
+          p_SetCoeff(subst,n_Init(-1,currRing->cf),currRing);
+        }
+        p_Setm(subst,currRing); pTest(subst);
+        pNext(g)=p_LmDeleteAndNext(pNext(g),currRing);
+        pNext(g)=p_Add_q(pNext(g),subst,currRing);
+        pTest(pNext(g));
+      }
+      // either way, add monomial in x to done
+      pNext(doneEnd)=p_LmInit(pNext(g),currRing);
+      pIter(doneEnd);
+      p_SetExp(doneEnd,1,0,currRing);
+      p_SetCoeff(doneEnd,n_Init(1,currRing->cf),currRing);
+      p_Setm(doneEnd,currRing);
+    }
+    pIter(g);
+  }
+  p_Delete(&done,currRing);
+  return false;
+}
+
+
+#ifndef NDEBUG
+BOOLEAN preduce(leftv res, leftv args)
+{
+  leftv u = args;
+  if ((u != NULL) && (u->Typ() == POLY_CMD))
+  {
+    poly g; bool b;
+    number p = n_Init(3,currRing->cf);
+    omUpdateInfo();
+    Print("usedBytes=%ld\n",om_Info.UsedBytes);
+    g = (poly) u->CopyD();
+    b = preduce(g,p);
+    p_Delete(&g,currRing);
+    if (b) return TRUE;
+    omUpdateInfo();
+    Print("usedBytes=%ld\n",om_Info.UsedBytes);
+    n_Delete(&p,currRing->cf);
+    res->rtyp = NONE;
+    res->data = NULL;
+    return FALSE;
+  }
+  return TRUE;
+}
+#endif //NDEBUG
+
+
+/***
+ * Returns the highest term in g with the matching x-monomial to m
+ * or, if it does not exist, the NULL pointer
+ **/
+static poly highestMatchingX(poly g, const poly m)
+{
+  pTest(g); pTest(m);
+  poly xalpha=p_LmInit(m,currRing);
+
+  // go along g and find the first term
+  // with the same monomial in x as xalpha
+  while (g)
+  {
+    p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing);
+    p_Setm(xalpha,currRing);
+    if (p_LmEqual(g,xalpha,currRing))
+      break;
+    pIter(g);
+  }
+
+  // gCache now either points at the wanted term
+  // or is NULL
+  p_Delete(&xalpha,currRing);
+  pTest(g);
+  return g;
+}
+
+
+/***
+ * Given g with lm(g)=t^\beta x^\alpha, returns g_\alpha
+ **/
+static poly powerSeriesCoeff(poly g)
+{
+  // the first term obviously is part of our output
+  // so we copy it...
+  poly xalpha=p_LmInit(g,currRing);
+  poly coeff=p_One(currRing);
+  p_SetCoeff(coeff,n_Copy(p_GetCoeff(g,currRing),currRing->cf),currRing);
+  p_SetExp(coeff,1,p_GetExp(g,1,currRing),currRing);
+  p_Setm(coeff,currRing);
+
+  // ..and proceed with the remaining terms,
+  // appending the relevant terms to coeff via coeffCache
+  poly coeffCache=coeff;
+  pIter(g);
+  while (g)
+  {
+    p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing);
+    p_Setm(xalpha,currRing);
+    if (p_LmEqual(g,xalpha,currRing))
+    {
+      pNext(coeffCache)=p_Init(currRing);
+      pIter(coeffCache);
+      p_SetCoeff(coeffCache,n_Copy(p_GetCoeff(g,currRing),currRing->cf),currRing);
+      p_SetExp(coeffCache,1,p_GetExp(g,1,currRing),currRing);
+      p_Setm(coeffCache,currRing);
+    }
+    pIter(g);
+  }
+
+  p_Delete(&xalpha,currRing);
+  return coeff;
+}
+
+
+/***
+ * divides g by t^b knowing that each term of g
+ * is divisible by t^b, i.e. no divisibility checks
+ * needed
+ **/
+static void divideByT(poly g, const long b)
+{
+  while (g)
+  {
+    p_SetExp(g,1,p_GetExp(g,1,currRing)-b,currRing);
+    p_Setm(g,currRing);
+    pIter(g);
+  }
+}
+
+
+static void divideByGcd(poly &g)
+{
+  // first determine all g_\alpha
+  // as alpha runs over all exponent vectors
+  // of monomials in x occuring in g
+
+  // gAlpha will store all g_\alpha,
+  // the first term will, for comparison purposes for now,
+  // also keep their monomial in x.
+  // we assume that the weight on the x are positive
+  // so that the x won't make the monomial smaller
+  ideal gAlphaFront = idInit(pLength(g));
+  gAlphaFront->m[0] = p_Head(g,currRing);
+  p_SetExp(gAlphaFront->m[0],1,0,currRing);
+  p_Setm(gAlphaFront->m[0],currRing);
+  ideal gAlphaEnd = idInit(pLength(g));
+  gAlphaEnd->m[0] = gAlphaFront->m[0];
+  unsigned long gAlpha_sev[pLength(g)];
+  gAlpha_sev[0] = p_GetShortExpVector(g,currRing);
+  long beta[pLength(g)];
+  beta[0] = p_GetExp(g,1,currRing);
+  int count=0;
+
+  poly current = pNext(g); unsigned long current_not_sev;
+  int i;
+  while (current)
+  {
+    // see whether the monomial in x of current already came up
+    // since everything is homogeneous in x and the ordering is local in t,
+    // this comes down to a divisibility test in two stages
+    current_not_sev = ~p_GetShortExpVector(current,currRing);
+    for(i=0; i<=count; i++)
+    {
+      // first stage using short exponent vectors
+      // second stage a proper test
+      if (p_LmShortDivisibleBy(gAlphaFront->m[i],gAlpha_sev[i],current,current_not_sev, currRing)
+            && p_LmDivisibleBy(gAlphaFront->m[i],current,currRing))
+      {
+        // if it already occured, add the term to the respective entry
+        // without the x part
+        pNext(gAlphaEnd->m[i])=p_Init(currRing);
+        pIter(gAlphaEnd->m[i]);
+        p_SetExp(gAlphaEnd->m[i],1,p_GetExp(current,1,currRing)-beta[i],currRing);
+        p_SetCoeff(gAlphaEnd->m[i],n_Copy(p_GetCoeff(current,currRing),currRing->cf),currRing);
+        p_Setm(gAlphaEnd->m[i],currRing);
+        break;
+      }
+    }
+    if (i>count)
+    {
+      // if it is new, create a new entry for the term
+      // including the monomial in x
+      count++;
+      gAlphaFront->m[count] = p_Head(current,currRing);
+      beta[count] = p_GetExp(current,1,currRing);
+      gAlphaEnd->m[count] = gAlphaFront->m[count];
+      gAlpha_sev[count] = p_GetShortExpVector(current,currRing);
+    }
+
+    pIter(current);
+  }
+
+  if (count)
+  {
+    // now remove all the x in the leading terms
+    // so that gAlpha only contais terms in t
+    int j; long d;
+    for (i=0; i<=count; i++)
+    {
+      for (j=2; j<=currRing->N; j++)
+        p_SetExp(gAlphaFront->m[i],j,0,currRing);
+      p_Setm(gAlphaFront->m[i],currRing);
+      gAlphaEnd->m[i]=NULL;
+    }
+
+    // now compute the gcd over all g_\alpha
+    // and set the input to null as they are deleted in the process
+    poly gcd = singclap_gcd(gAlphaFront->m[0],gAlphaFront->m[1],currRing);
+    gAlphaFront->m[0] = NULL;
+    gAlphaFront->m[1] = NULL;
+    for(i=2; i<=count; i++)
+    {
+      gcd = singclap_gcd(gcd,gAlphaFront->m[i],currRing);
+      gAlphaFront->m[i] = NULL;
+    }
+    // divide g by the gcd
+    poly h = singclap_pdivide(g,gcd,currRing);
+    p_Delete(&gcd,currRing);
+    p_Delete(&g,currRing);
+    g = h;
+
+    id_Delete(&gAlphaFront,currRing);
+    id_Delete(&gAlphaEnd,currRing);
+  }
+  else
+  {
+    // if g contains only one monomial in x,
+    // then we can get rid of all the higher t
+    p_Delete(&g,currRing);
+    g = gAlphaFront->m[0];
+    pIter(gAlphaFront->m[0]);
+    pNext(g)=NULL;
+    gAlphaEnd->m[0] = NULL;
+    id_Delete(&gAlphaFront,currRing);
+    id_Delete(&gAlphaEnd,currRing);
+  }
+}
+
+
+/***
+ * 1. For each \alpha\in\NN^n, changes (c_\beta t^\beta + c_{\beta+1} t^{\beta+1} + ...)
+ *    to (c_\beta + c_{\beta+1}*p + ...) t^\beta
+ * 2. Computes the gcd over all (c_\beta + c_{\beta+1}*p + ...) and divides g by it
+ **/
+static void simplify(poly g, const number p)
+{
+  // go along g and look for relevant terms.
+  // monomials in x which have been checked are tracked in done.
+  poly done = p_LmInit(g,currRing);
+  p_SetCoeff(done,n_Init(1,currRing->cf),currRing);
+  p_Setm(done,currRing);
+  poly doneEnd = done;
+  poly doneCurrent;
+
+  poly subst; number substCoeff, substCoeffCache;
+  unsigned long d;
+
+  poly gCurrent = g;
+  while(pNext(gCurrent))
+  {
+    // check whether next term needs to be reduced:
+    // first, check whether monomial in x has come up yet
+    for (doneCurrent=done; doneCurrent; pIter(doneCurrent))
+    {
+      if (p_LmDivisibleBy(doneCurrent,pNext(gCurrent),currRing))
+        break;
+    }
+    // if the monomial in x already occured, then we want to replace
+    // as many t with p as theoretically feasable
+    if (doneCurrent)
+    {
+      // suppose pNext(gCurrent)=3*t5x and doneCurrent=t3x
+      // then we want to replace pNext(gCurrent) with 3p2*t3x
+      // subst = ?*t3x
+      subst = p_LmInit(doneCurrent,currRing);
+      // substcoeff = p2
+      n_Power(p,p_GetExp(subst,1,currRing)-p_GetExp(doneCurrent,1,currRing),&substCoeff,currRing->cf);
+      // substcoeff = 3p2
+      n_InpMult(substCoeff,p_GetCoeff(pNext(gCurrent),currRing),currRing->cf);
+      // subst = 3p2*t3x
+      p_SetCoeff(subst,substCoeff,currRing);
+      p_Setm(subst,currRing); pTest(subst);
+
+      // g = g - pNext(gCurrent) + subst
+      pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing);
+      g=p_Add_q(g,subst,currRing);
+      pTest(pNext(gbeginning));
+    }
+    else
+    {
+      // if the monomial in x is brand new,
+      // then we check whether the coefficient is divisible by p
+      if (n_DivBy(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf))
+      {
+        // reduce the term with respect to p-t:
+        // suppose pNext(gCurrent)=4p3*tx
+        // then we want to replace it with 4*t4x
+        // divide 4p3 repeatedly by p until it is not divisible anymore,
+        // keeping track on the final value 4
+        // and the number of times it has been divided 3
+        substCoeff = n_Div(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf);
+        d = 1;
+        while (n_DivBy(substCoeff,p,currRing->cf))
+        {
+          substCoeffCache = n_Div(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf);
+          n_Delete(&substCoeff,currRing->cf);
+          substCoeff = substCoeffCache;
+          d++;
+          assume(d>0); // check for overflow
+        }
+
+        // subst = ?*tx
+        subst=p_LmInit(pNext(gCurrent),currRing);
+        // subst = ?*t4x
+        p_AddExp(subst,1,d,currRing);
+        // subst = 4*t4x
+        p_SetCoeff(subst,substCoeffCache,currRing);
+        p_Setm(subst,currRing); pTest(subst);
+
+        // g = g - pNext(gCurrent) + subst
+        pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing);
+        pNext(gCurrent)=p_Add_q(pNext(gCurrent),subst,currRing);
+        pTest(pNext(gCurrent));
+      }
+
+      // either way, add monomial in x with minimal t to done
+      pNext(doneEnd)=p_LmInit(pNext(gCurrent),currRing);
+      pIter(doneEnd);
+      p_SetCoeff(doneEnd,n_Init(1,currRing->cf),currRing);
+      p_Setm(doneEnd,currRing);
+    }
+    pIter(gCurrent);
+  }
+  p_Delete(&done,currRing);
+}
+
+
+#ifndef NDEBUG
+BOOLEAN divideByGcd(leftv res, leftv args)
+{
+  leftv u = args;
+  if ((u != NULL) && (u->Typ() == POLY_CMD))
+  {
+    poly g;
+    omUpdateInfo();
+    Print("usedBytes1=%ld\n",om_Info.UsedBytes);
+    g = (poly) u->CopyD();
+    divideByGcd(g);
+    p_Delete(&g,currRing);
+    omUpdateInfo();
+    Print("usedBytes1=%ld\n",om_Info.UsedBytes);
+    res->rtyp = NONE;
+    res->data = NULL;
+    return FALSE;
+  }
+  return TRUE;
+}
+#endif //NDEBUG
+
+
+BOOLEAN initialReduction(leftv res, leftv args)
+{
+  leftv u = args;
+  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
+  {
+    leftv v = u->next;
+    if (v == NULL)
+    {
+      ideal I = (ideal) u->Data();
+
+      /***
+       * Suppose I=g_0,...,g_n and g_0=p-t.
+       * Step 1: sort elements g_1,...,g_{n-1} such that lm(g_1)>...>lm(g_{n-1})
+       * (the following algorithm is a bubble sort)
+       **/
+      sortingLaterGenerators(I);
+
+      /***
+       * Step 2: replace coefficient of terms of lowest t-degree divisible by p with t
+       **/
+      number p = p_GetCoeff(I->m[0],currRing);
+      for (int i=1; i<IDELEMS(I); i++)
+      {
+        if (preduce(I->m[i],p))
+          return TRUE;
+      }
+
+      /***
+       * Step 3: the first pass. removing terms with the same monomials in x as lt(g_i)
+       * out of g_j for i<j
+       **/
+      int i,j;
+      poly uniti, unitj;
+      for (i=1; i<IDELEMS(I)-1; i++)
+      {
+        for (j=i+1; j<IDELEMS(I); j++)
+        {
+          unitj = highestMatchingX(I->m[j],I->m[i]);
+          if (unitj)
+          {
+            unitj = powerSeriesCoeff(unitj);
+            divideByT(unitj,p_GetExp(I->m[i],1,currRing));
+            uniti = powerSeriesCoeff(I->m[i]);
+            divideByT(uniti,p_GetExp(I->m[i],1,currRing));
+            pTest(unitj); pTest(uniti); pTest(I->m[j]); pTest(I->m[i]);
+            I->m[j]=p_Add_q(p_Mult_q(uniti,I->m[j],currRing),
+                            p_Neg(p_Mult_q(unitj,p_Copy(I->m[i],currRing),currRing),currRing),
+                            currRing);
+            divideByGcd(I->m[j]);
+          }
+        }
+      }
+      for (int i=1; i<IDELEMS(I); i++)
+      {
+        if (preduce(I->m[i],p))
+          return TRUE;
+      }
+
+      /***
+       * Step 4: the second pass. removing terms divisible by lt(g_j) out of g_i for i<j
+       **/
+      for (i=1; i<IDELEMS(I)-1; i++)
+      {
+        for (j=i+1; j<IDELEMS(I); j++)
+        {
+          uniti = highestMatchingX(I->m[i],I->m[j]);
+          if (uniti && p_GetExp(uniti,1,currRing)>=p_GetExp(I->m[j],1,currRing))
+          {
+            uniti = powerSeriesCoeff(uniti);
+            divideByT(uniti,p_GetExp(I->m[j],1,currRing));
+            unitj = powerSeriesCoeff(I->m[j]);
+            divideByT(unitj,p_GetExp(I->m[j],1,currRing));
+            I->m[i] = p_Add_q(p_Mult_q(unitj,I->m[i],currRing),
+                              p_Neg(p_Mult_q(uniti,p_Copy(I->m[j],currRing),currRing),currRing),
+                              currRing);
+            divideByGcd(I->m[j]);
+          }
+        }
+      }
+      for (int i=1; i<IDELEMS(I); i++)
+      {
+        if (preduce(I->m[i],p))
+          return TRUE;
+      }
+
+      res->rtyp = NONE;
+      res->data = NULL;
+      IDDATA((idhdl)u->data) = (char*) I;
+      return FALSE;
+    }
+  }
+  WerrorS("initialReduction: unexpected parameters");
+  return TRUE;
+}
+
+
+#if 0
+/***
+ * Given a general ring r with any ordering,
+ * changes the ordering to a(v),ws(-w)
+ **/
+bool changetoAWSRing(ring r, gfan::ZVector v, gfan::ZVector w)
+{
+  omFree(r->order);
+  r->order  = (int*) omAlloc0(4*sizeof(int));
+  omFree(r->block0);
+  r->block0 = (int*) omAlloc0(4*sizeof(int));
+  omFree(r->block1);
+  r->block1 = (int*) omAlloc0(4*sizeof(int));
+  for (int i=0; r->wvhdl[i]; i++)
+  { omFree(r->wvhdl[i]); }
+  omFree(r->wvhdl);
+  r->wvhdl  = (int**) omAlloc0(4*sizeof(int*));
+
+  bool ok = false;
+  r->order[0]  = ringorder_a;
+  r->block0[0] = 1;
+  r->block1[0] = r->N;
+  r->wvhdl[0]  = ZVectorToIntStar(v,ok);
+  r->order[1]  = ringorder_ws;
+  r->block0[1] = 1;
+  r->block1[1] = r->N;
+  r->wvhdl[1]  = ZVectorToIntStar(w,ok);
+  r->order[2]=ringorder_C;
+  return ok;
+}
+
+
+/***
+ * Given a ring with ordering a(v'),ws(w'),
+ * changes the weights to v,w
+ **/
+bool changeAWSWeights(ring r, gfan::ZVector v, gfan::ZVector w)
+{
+  omFree(r->wvhdl[0]);
+  omFree(r->wvhdl[1]);
+  bool ok = false;
+  r->wvhdl[0]  = ZVectorToIntStar(v,ok);
+  r->wvhdl[1]  = ZVectorToIntStar(w,ok);
+  return ok;
+}
+
+
+// /***
+//  * Creates an int* representing the transposition of the last two variables
+//  **/
+// static inline int* createPermutationVectorForSaturation(static const ring &r)
+// {
+//   int* w = (int*) omAlloc0((rVar(r)+1)*sizeof(int));
+//   for (int i=1; i<=rVar(r)-2; i++)
+//     w[i] = i;
+//   w[rVar(r)-1] = rVar(r);
+//   w[rVar(r)] = rVar(r)-1;
+// }
+
+
+/***
+ * Creates an int* representing the permutation
+ * 1 -> 1, ..., i-1 -> i-1, i -> n, i+1 -> n-1, ... , n -> i
+ **/
+static inline int* createPermutationVectorForSaturation(const ring &r, const int i)
+{
+  int* sigma = (int*) omAlloc0((rVar(r)+1)*sizeof(int));
+  int j;
+  for (j=1; j<i; j++)
+    sigma[j] = j;
+  for (; j<=rVar(r); j++)
+    sigma[j] = rVar(r)-j+i;
+  return(sigma);
+}
+
+
+/***
+ * Changes the int* representing the permutation
+ * 1 -> 1, ..., i -> i, i+1 -> n, i+2 -> n-1, ... , n -> i+1
+ * to an int* representing the permutation
+ * 1 -> 1, ..., i-1 -> i-1, i -> n, i+1 -> n-1, ... , n -> i
+ **/
+static void changePermutationVectorForSaturation(int* sigma, const ring &r, const int i)
+{
+  for (int j=i; j<rVar(r); j++)
+    sigma[j] = rVar(r)-j+i;
+  sigma[rVar(r)] = i;
+}
+
+
+/***
+ * returns a ring in which the weights of the ring variables are permuted
+ * if handed over a poly in which the variables are permuted, this is basically
+ * as good as permuting the variables of the ring itself.
+ **/
+static ring permuteWeighstOfRingVariables(const ring &r, const int* const sigma)
+{
+  ring s = rCopy0(r);
+  for (int j=0; j<rVar(r); j++)
+  {
+    s->wvhdl[0][j] = r->wvhdl[0][sigma[j+1]];
+    s->wvhdl[1][j] = r->wvhdl[1][sigma[j+1]];
+  }
+  rComplete(s,1);
+  return s;
+}
+
+
+/***
+ * creates a ring s that is a copy of r except with ordering wp(w)
+ **/
+static inline ring createInitialRingForSaturation(const ring &r, const gfan::ZVector &w, bool &ok)
+{
+  assume(rVar(r) == (int) w.size());
+
+  ring s = rCopy0(r); int i;
+  for (i=0; s->order[i]; i++)
+    omFreeSize(s->wvhdl[i],rVar(r)*sizeof(int));
+  i++;
+  omFreeSize(s->order,i*sizeof(int));
+  s->order  = (int*) omAlloc0(3*sizeof(int));
+  omFreeSize(s->block0,i*sizeof(int));
+  s->block0 = (int*) omAlloc0(3*sizeof(int));
+  omFreeSize(s->block1,i*sizeof(int));
+  s->block1 = (int*) omAlloc0(3*sizeof(int));
+  omFreeSize(s->wvhdl,i*sizeof(int*));
+  s->wvhdl  = (int**) omAlloc0(3*sizeof(int*));
+
+  s->order[0]  = ringorder_wp;
+  s->block0[0] = 1;
+  s->block1[0] = rVar(r);
+  s->wvhdl[0]  = ZVectorToIntStar(w,ok);
+  s->order[1]=ringorder_C;
+
+  rComplete(s,1);
+  return s;
+}
+
+
+/***
+ * Given an weighted homogeneous ideal I with respect to weight w
+ * that in standard basis form with respect to the ordering ws(-w),
+ * derives the standard basis of I:<x_n>^\infty
+ * and returns a long k such that I:<x_n>^\infty=I:<x_n>^k
+ **/
+static long deriveStandardBasisOfSaturation(ideal &I, ring &r)
+{
+  long k=0, l; poly current;
+  for (int i=0; i<IDELEMS(I); i++)
+  {
+    current = I->m[i];
+    l = p_GetExp(current,rVar(r),r);
+    if (k<l) k=l;
+    while (current)
+    {
+      p_SubExp(current,rVar(r),l,r); p_Setm(current,r);
+      pIter(current);
+    }
+  }
+  return k;
+}
+
+
+/***
+ * Given a weighted homogeneous ideal I with respect to weight w
+ * with constant first element,
+ * returns NULL if I does not contain a monomial
+ * otherwise returns the monomial contained in I
+ **/
+poly containsMonomial(const ideal &I, const gfan::ZVector &w)
+{
+  assume(rField_is_Ring_Z(currRing));
+
+  // first we switch to the ground field currRing->cf / I->m[0]
+  ring r = rCopy0(currRing);
+  nKillChar(r->cf);
+  r->cf = nInitChar(n_Zp,(void*)(long)n_Int(p_GetCoeff(I->m[0],currRing),currRing->cf));
+  rComplete(r);
+
+  ideal J = id_Copy(I, currRing); poly cache; number temp;
+  for (int i=0; i<IDELEMS(I); i++)
+  {
+    cache = J->m[i];
+    while (cache)
+    {
+      // TODO: temp = npMapGMP(p_GetCoeff(cache,currRing),currRing->cf,r->cf);
+      p_SetCoeff(cache,temp,r); pIter(cache);
+    }
+  }
+
+
+  J = kStd(J,NULL,isHomog,NULL);
+
+  bool b = false;
+  ring s = createInitialRingForSaturation(currRing, w, b);
+  if (b)
+  {
+    WerrorS("containsMonomial: overflow in weight vector");
+    return NULL;
+  }
+
+  return NULL;
+}
+
+
+gfan::ZCone* startingCone(ideal I)
+{
+  I = kStd(I,NULL,isNotHomog,NULL);
+  gfan::ZCone* zc = maximalGroebnerCone(currRing,I);
+  gfan::ZMatrix rays = zc->extremeRays();
+  gfan::ZVector v;
+  for (int i=0; i<rays.getHeight(); i++)
+  {
+    v = rays[i];
+  }
+  return zc;
+}
+#endif
 
 
@@ -202 +968 @@
   p->iiAddCproc("","maximalGroebnerCone",FALSE,maximalGroebnerCone);
   p->iiAddCproc("","initial",FALSE,initial);
+#ifndef NDEBUG
+  p->iiAddCproc("","divideByGcd",FALSE,divideByGcd);
+  p->iiAddCproc("","preduce",FALSE,preduce);
+#endif //NDEBUG
+  p->iiAddCproc("","initialReduction",FALSE,initialReduction);
   p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
 }