Index: Singular/dyn_modules/gfanlib/tropical.cc
===================================================================
--- Singular/dyn_modules/gfanlib/tropical.cc	(revision d770e6a687e6451b6e42d8b59cf74426e947eaf0)
+++ Singular/dyn_modules/gfanlib/tropical.cc	(revision fa167e9991d970432d089c5d6b7094886d1ffbb7)
@@ -5,4 +5,5 @@
 #include <bbcone.h>
 
+#include <initialReduction.h>
 
 poly initial(poly p)
@@ -236,509 +237,509 @@
 
 
-/***
- * 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;
-}
+// /***
+//  * 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;
+// }
 
 
@@ -969,8 +970,8 @@
   p->iiAddCproc("","initial",FALSE,initial);
 #ifndef NDEBUG
-  p->iiAddCproc("","divideByGcd",FALSE,divideByGcd);
-  p->iiAddCproc("","preduce",FALSE,preduce);
+  // p->iiAddCproc("","divideByGcd",FALSE,divideByGcd);
+  p->iiAddCproc("","pReduce",FALSE,pReduce);
 #endif //NDEBUG
-  p->iiAddCproc("","initialReduction",FALSE,initialReduction);
+  // p->iiAddCproc("","initialReduction",FALSE,initialReduction);
   p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
 }
