Index: configure.ac
===================================================================
--- configure.ac	(revision 66b1de771706655cf8a9c3924d1dac08150c8814)
+++ configure.ac	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
@@ -204,4 +204,5 @@
 AC_CONFIG_FILES([kernel/fglm/Makefile])
 AC_CONFIG_FILES([kernel/groebner_walk/Makefile])
+AC_CONFIG_FILES([kernel/combinatorics/Makefile])
 
 AC_CONFIG_SUBDIRS([gfanlib])
Index: kernel/Makefile.am
===================================================================
--- kernel/Makefile.am	(revision 66b1de771706655cf8a9c3924d1dac08150c8814)
+++ kernel/Makefile.am	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
@@ -2,5 +2,5 @@
 # TODO: use ${top_srcdir} instead of .. in the above?
 
-SUBDIRS=sample numeric fglm groebner_walk
+SUBDIRS=sample numeric fglm groebner_walk combinatorics
 # kernelsample   
 
@@ -13,5 +13,4 @@
 
 SOURCES = polys.cc febase.cc feread.cc \
-    hdegree.cc hilb.cc hutil.cc \
     ideals.cc \
     khstd.cc kstdfac.cc \
@@ -62,4 +61,5 @@
 ${builddir}/fglm/libfglm.la \
 ${builddir}/groebner_walk/libgroebner_walk.la \
+${builddir}/combinatorics/libcombinatorics.la \
 ${USEPPROCSDYNAMICLDFLAGS} ${USEPPROCSDYNAMICLD} \
 ${abs_top_builddir}/libpolys/polys/libpolys.la
Index: kernel/combinatorics/Makefile.am
===================================================================
--- kernel/combinatorics/Makefile.am	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
+++ kernel/combinatorics/Makefile.am	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
@@ -0,0 +1,13 @@
+ACLOCAL_AMFLAGS=-I ../../m4
+AM_CPPFLAGS=-I${top_srcdir} -I${top_builddir} -I${top_srcdir}/libpolys -I${top_builddir}/libpolys \
+${FACTORY_INCLUDES} $(OMALLOC_INCLUDES) $(RESOURCES_INCLUDES) \
+${NTL_CFLAGS} $(FLINT_CFLAGS) ${GMP_CFLAGS}
+
+noinst_LTLIBRARIES=libcombinatorics.la
+libcombinatorics_la_SOURCES=hdegree.cc hilb.cc hutil.cc
+
+# libcombinatorics_la_includedir=$(includedir)/kernel/combinatorics
+# libcombinatorics_la_include_HEADERS=hutil.h
+
+
+
Index: kernel/combinatorics/hdegree.cc
===================================================================
--- kernel/combinatorics/hdegree.cc	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
+++ kernel/combinatorics/hdegree.cc	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
@@ -0,0 +1,1464 @@
+/****************************************
+*  Computer Algebra System SINGULAR     *
+****************************************/
+/*
+*  ABSTRACT -  dimension, multiplicity, HC, kbase
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "singularconfig.h"
+#endif /* HAVE_CONFIG_H */
+#include <kernel/mod2.h>
+
+
+#include <kernel/structs.h>
+#include <kernel/febase.h>
+#include <omalloc/omalloc.h>
+#include <kernel/ideals.h>
+#include <kernel/polys.h>
+#include <misc/intvec.h>
+#include <coeffs/numbers.h>
+#include <kernel/hutil.h>
+#include <kernel/stairc.h>
+
+int  hCo, hMu, hMu2;
+omBin indlist_bin = omGetSpecBin(sizeof(indlist));
+
+/*0 implementation*/
+
+// dimension
+
+void hDimSolve(scmon pure, int Npure, scfmon rad, int Nrad,
+ varset var, int Nvar)
+{
+  int  dn, iv, rad0, b, c, x;
+  scmon pn;
+  scfmon rn;
+  if (Nrad < 2)
+  {
+    dn = Npure + Nrad;
+    if (dn < hCo)
+      hCo = dn;
+    return;
+  }
+  if (Npure+1 >= hCo)
+    return;
+  iv = Nvar;
+  while(pure[var[iv]]) iv--;
+  hStepR(rad, Nrad, var, iv, &rad0);
+  if (rad0!=0)
+  {
+    iv--;
+    if (rad0 < Nrad)
+    {
+      pn = hGetpure(pure);
+      rn = hGetmem(Nrad, rad, radmem[iv]);
+      hDimSolve(pn, Npure + 1, rn, rad0, var, iv);
+      b = rad0;
+      c = Nrad;
+      hElimR(rn, &rad0, b, c, var, iv);
+      hPure(rn, b, &c, var, iv, pn, &x);
+      hLex2R(rn, rad0, b, c, var, iv, hwork);
+      rad0 += (c - b);
+      hDimSolve(pn, Npure + x, rn, rad0, var, iv);
+    }
+    else
+    {
+      hDimSolve(pure, Npure, rad, Nrad, var, iv);
+    }
+  }
+  else
+    hCo = Npure + 1;
+}
+
+int  scDimInt(ideal S, ideal Q)
+{
+  int  mc;
+  hexist = hInit(S, Q, &hNexist);
+  if (!hNexist)
+    return (currRing->N);
+  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
+  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  mc = hisModule;
+  if (!mc)
+  {
+    hrad = hexist;
+    hNrad = hNexist;
+  }
+  else
+    hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  radmem = hCreate((currRing->N) - 1);
+  hCo = (currRing->N) + 1;
+  loop
+  {
+    if (mc)
+      hComp(hexist, hNexist, mc, hrad, &hNrad);
+    if (hNrad)
+    {
+      hNvar = (currRing->N);
+      hRadical(hrad, &hNrad, hNvar);
+      hSupp(hrad, hNrad, hvar, &hNvar);
+      if (hNvar)
+      {
+        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
+        hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
+        hLexR(hrad, hNrad, hvar, hNvar);
+        hDimSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
+      }
+    }
+    else
+    {
+      hCo = 0;
+      break;
+    }
+    mc--;
+    if (mc <= 0)
+      break;
+  }
+  hKill(radmem, (currRing->N) - 1);
+  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
+  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
+  hDelete(hexist, hNexist);
+  if (hisModule)
+    omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
+  return (currRing->N) - hCo;
+}
+
+// independent set
+static scmon hInd;
+
+static void hIndSolve(scmon pure, int Npure, scfmon rad, int Nrad,
+ varset var, int Nvar)
+{
+  int  dn, iv, rad0, b, c, x;
+  scmon pn;
+  scfmon rn;
+  if (Nrad < 2)
+  {
+    dn = Npure + Nrad;
+    if (dn < hCo)
+    {
+      hCo = dn;
+      for (iv=(currRing->N); iv; iv--)
+      {
+        if (pure[iv])
+          hInd[iv] = 0;
+        else
+          hInd[iv] = 1;
+      }
+      if (Nrad)
+      {
+        pn = *rad;
+        iv = Nvar;
+        loop
+        {
+          x = var[iv];
+          if (pn[x])
+          {
+            hInd[x] = 0;
+            break;
+          }
+          iv--;
+        }
+      }
+    }
+    return;
+  }
+  if (Npure+1 >= hCo)
+    return;
+  iv = Nvar;
+  while(pure[var[iv]]) iv--;
+  hStepR(rad, Nrad, var, iv, &rad0);
+  if (rad0)
+  {
+    iv--;
+    if (rad0 < Nrad)
+    {
+      pn = hGetpure(pure);
+      rn = hGetmem(Nrad, rad, radmem[iv]);
+      pn[var[iv + 1]] = 1;
+      hIndSolve(pn, Npure + 1, rn, rad0, var, iv);
+      pn[var[iv + 1]] = 0;
+      b = rad0;
+      c = Nrad;
+      hElimR(rn, &rad0, b, c, var, iv);
+      hPure(rn, b, &c, var, iv, pn, &x);
+      hLex2R(rn, rad0, b, c, var, iv, hwork);
+      rad0 += (c - b);
+      hIndSolve(pn, Npure + x, rn, rad0, var, iv);
+    }
+    else
+    {
+      hIndSolve(pure, Npure, rad, Nrad, var, iv);
+    }
+  }
+  else
+  {
+    hCo = Npure + 1;
+    for (x=(currRing->N); x; x--)
+    {
+      if (pure[x])
+        hInd[x] = 0;
+      else
+        hInd[x] = 1;
+    }
+    hInd[var[iv]] = 0;
+  }
+}
+
+intvec * scIndIntvec(ideal S, ideal Q)
+{
+  intvec *Set=new intvec((currRing->N));
+  int  mc,i;
+  hexist = hInit(S, Q, &hNexist);
+  if (hNexist==0)
+  {
+    for(i=0; i<(currRing->N); i++)
+      (*Set)[i]=1;
+    return Set;
+  }
+  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
+  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  hInd = (scmon)omAlloc0((1 + (currRing->N)) * sizeof(int));
+  mc = hisModule;
+  if (mc==0)
+  {
+    hrad = hexist;
+    hNrad = hNexist;
+  }
+  else
+    hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  radmem = hCreate((currRing->N) - 1);
+  hCo = (currRing->N) + 1;
+  loop
+  {
+    if (mc!=0)
+      hComp(hexist, hNexist, mc, hrad, &hNrad);
+    if (hNrad!=0)
+    {
+      hNvar = (currRing->N);
+      hRadical(hrad, &hNrad, hNvar);
+      hSupp(hrad, hNrad, hvar, &hNvar);
+      if (hNvar!=0)
+      {
+        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
+        hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
+        hLexR(hrad, hNrad, hvar, hNvar);
+        hIndSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
+      }
+    }
+    else
+    {
+      hCo = 0;
+      break;
+    }
+    mc--;
+    if (mc <= 0)
+      break;
+  }
+  for(i=0; i<(currRing->N); i++)
+    (*Set)[i] = hInd[i+1];
+  hKill(radmem, (currRing->N) - 1);
+  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  omFreeSize((ADDRESS)hInd, (1 + (currRing->N)) * sizeof(int));
+  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
+  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
+  hDelete(hexist, hNexist);
+  if (hisModule)
+    omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
+  return Set;
+}
+
+indset ISet, JSet;
+
+static BOOLEAN hNotZero(scfmon rad, int Nrad, varset var, int Nvar)
+{
+  int  k1, i;
+  k1 = var[Nvar];
+  i = 0;
+  loop
+  {
+    if (rad[i][k1]==0)
+      return FALSE;
+    i++;
+    if (i == Nrad)
+      return TRUE;
+  }
+}
+
+static void hIndep(scmon pure)
+{
+  int iv;
+  intvec *Set;
+
+  Set = ISet->set = new intvec((currRing->N));
+  for (iv=(currRing->N); iv!=0 ; iv--)
+  {
+    if (pure[iv])
+      (*Set)[iv-1] = 0;
+    else
+      (*Set)[iv-1] = 1;
+  }
+  ISet = ISet->nx = (indset)omAlloc0Bin(indlist_bin);
+  hMu++;
+}
+
+void hIndMult(scmon pure, int Npure, scfmon rad, int Nrad,
+ varset var, int Nvar)
+{
+  int  dn, iv, rad0, b, c, x;
+  scmon pn;
+  scfmon rn;
+  if (Nrad < 2)
+  {
+    dn = Npure + Nrad;
+    if (dn == hCo)
+    {
+      if (Nrad==0)
+        hIndep(pure);
+      else
+      {
+        pn = *rad;
+        for (iv = Nvar; iv!=0; iv--)
+        {
+          x = var[iv];
+          if (pn[x])
+          {
+            pure[x] = 1;
+            hIndep(pure);
+            pure[x] = 0;
+          }
+        }
+      }
+    }
+    return;
+  }
+  iv = Nvar;
+  dn = Npure+1;
+  if (dn >= hCo)
+  {
+    if (dn > hCo)
+      return;
+    loop
+    {
+      if(!pure[var[iv]])
+      {
+        if(hNotZero(rad, Nrad, var, iv))
+        {
+          pure[var[iv]] = 1;
+          hIndep(pure);
+          pure[var[iv]] = 0;
+        }
+      }
+      iv--;
+      if (!iv)
+        return;
+    }
+  }
+  while(pure[var[iv]]) iv--;
+  hStepR(rad, Nrad, var, iv, &rad0);
+  iv--;
+  if (rad0 < Nrad)
+  {
+    pn = hGetpure(pure);
+    rn = hGetmem(Nrad, rad, radmem[iv]);
+    pn[var[iv + 1]] = 1;
+    hIndMult(pn, Npure + 1, rn, rad0, var, iv);
+    pn[var[iv + 1]] = 0;
+    b = rad0;
+    c = Nrad;
+    hElimR(rn, &rad0, b, c, var, iv);
+    hPure(rn, b, &c, var, iv, pn, &x);
+    hLex2R(rn, rad0, b, c, var, iv, hwork);
+    rad0 += (c - b);
+    hIndMult(pn, Npure + x, rn, rad0, var, iv);
+  }
+  else
+  {
+    hIndMult(pure, Npure, rad, Nrad, var, iv);
+  }
+}
+
+/*3
+* consider indset x := !pure
+* (for all i) (if(sm(i) > x) return FALSE)
+* else return TRUE
+*/
+static BOOLEAN hCheck1(indset sm, scmon pure)
+{
+  int iv;
+  intvec *Set;
+  while (sm->nx != NULL)
+  {
+    Set = sm->set;
+    iv=(currRing->N);
+    loop
+    {
+      if (((*Set)[iv-1] == 0) && (pure[iv] == 0))
+        break;
+      iv--;
+      if (iv == 0)
+        return FALSE;
+    }
+    sm = sm->nx;
+  }
+  return TRUE;
+}
+
+/*3
+* consider indset x := !pure
+* (for all i) if(x > sm(i)) delete sm(i))
+* return (place for x)
+*/
+static indset hCheck2(indset sm, scmon pure)
+{
+  int iv;
+  intvec *Set;
+  indset be, a1 = NULL;
+  while (sm->nx != NULL)
+  {
+    Set = sm->set;
+    iv=(currRing->N);
+    loop
+    {
+      if ((pure[iv] == 1) && ((*Set)[iv-1] == 1))
+        break;
+      iv--;
+      if (iv == 0)
+      {
+        if (a1 == NULL)
+        {
+          a1 = sm;
+        }
+        else
+        {
+          hMu2--;
+          be->nx = sm->nx;
+          delete Set;
+          omFreeBin((ADDRESS)sm, indlist_bin);
+          sm = be;
+        }
+        break;
+      }
+    }
+    be = sm;
+    sm = sm->nx;
+  }
+  if (a1 != NULL)
+  {
+    return a1;
+  }
+  else
+  {
+    hMu2++;
+    sm->set = new intvec((currRing->N));
+    sm->nx = (indset)omAlloc0Bin(indlist_bin);
+    return sm;
+  }
+}
+
+/*2
+*  definition x >= y
+*      x(i) == 0 => y(i) == 0
+*      > ex. j with x(j) == 1 and y(j) == 0
+*/
+static void hCheckIndep(scmon pure)
+{
+  intvec *Set;
+  indset res;
+  int iv;
+  if (hCheck1(ISet, pure))
+  {
+    if (hCheck1(JSet, pure))
+    {
+      res = hCheck2(JSet,pure);
+      if (res == NULL)
+        return;
+      Set = res->set;
+      for (iv=(currRing->N); iv; iv--)
+      {
+        if (pure[iv])
+          (*Set)[iv-1] = 0;
+        else
+          (*Set)[iv-1] = 1;
+      }
+    }
+  }
+}
+
+void hIndAllMult(scmon pure, int Npure, scfmon rad, int Nrad,
+ varset var, int Nvar)
+{
+  int  dn, iv, rad0, b, c, x;
+  scmon pn;
+  scfmon rn;
+  if (Nrad < 2)
+  {
+    dn = Npure + Nrad;
+    if (dn > hCo)
+    {
+      if (!Nrad)
+        hCheckIndep(pure);
+      else
+      {
+        pn = *rad;
+        for (iv = Nvar; iv; iv--)
+        {
+          x = var[iv];
+          if (pn[x])
+          {
+            pure[x] = 1;
+            hCheckIndep(pure);
+            pure[x] = 0;
+          }
+        }
+      }
+    }
+    return;
+  }
+  iv = Nvar;
+  while(pure[var[iv]]) iv--;
+  hStepR(rad, Nrad, var, iv, &rad0);
+  iv--;
+  if (rad0 < Nrad)
+  {
+    pn = hGetpure(pure);
+    rn = hGetmem(Nrad, rad, radmem[iv]);
+    pn[var[iv + 1]] = 1;
+    hIndAllMult(pn, Npure + 1, rn, rad0, var, iv);
+    pn[var[iv + 1]] = 0;
+    b = rad0;
+    c = Nrad;
+    hElimR(rn, &rad0, b, c, var, iv);
+    hPure(rn, b, &c, var, iv, pn, &x);
+    hLex2R(rn, rad0, b, c, var, iv, hwork);
+    rad0 += (c - b);
+    hIndAllMult(pn, Npure + x, rn, rad0, var, iv);
+  }
+  else
+  {
+    hIndAllMult(pure, Npure, rad, Nrad, var, iv);
+  }
+}
+
+// multiplicity
+
+static int hZeroMult(scmon pure, scfmon stc, int Nstc, varset var, int Nvar)
+{
+  int  iv = Nvar -1, sum, a, a0, a1, b, i;
+  int  x, x0;
+  scmon pn;
+  scfmon sn;
+  if (!iv)
+    return pure[var[1]];
+  else if (!Nstc)
+  {
+    sum = 1;
+    for (i = Nvar; i; i--)
+      sum *= pure[var[i]];
+    return sum;
+  }
+  x = a = 0;
+  pn = hGetpure(pure);
+  sn = hGetmem(Nstc, stc, stcmem[iv]);
+  hStepS(sn, Nstc, var, Nvar, &a, &x);
+  if (a == Nstc)
+    return pure[var[Nvar]] * hZeroMult(pn, sn, a, var, iv);
+  else
+    sum = x * hZeroMult(pn, sn, a, var, iv);
+  b = a;
+  loop
+  {
+    a0 = a;
+    x0 = x;
+    hStepS(sn, Nstc, var, Nvar, &a, &x);
+    hElimS(sn, &b, a0, a, var, iv);
+    a1 = a;
+    hPure(sn, a0, &a1, var, iv, pn, &i);
+    hLex2S(sn, b, a0, a1, var, iv, hwork);
+    b += (a1 - a0);
+    if (a < Nstc)
+    {
+      sum += (x - x0) * hZeroMult(pn, sn, b, var, iv);
+    }
+    else
+    {
+      sum += (pure[var[Nvar]] - x0) * hZeroMult(pn, sn, b, var, iv);
+      return sum;
+    }
+  }
+}
+
+static void hProject(scmon pure, varset sel)
+{
+  int  i, i0, k;
+  i0 = 0;
+  for (i = 1; i <= (currRing->N); i++)
+  {
+    if (pure[i])
+    {
+      i0++;
+      sel[i0] = i;
+    }
+  }
+  i = hNstc;
+  memcpy(hwork, hstc, i * sizeof(scmon));
+  hStaircase(hwork, &i, sel, i0);
+  if ((i0 > 2) && (i > 10))
+    hOrdSupp(hwork, i, sel, i0);
+  memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
+  hPure(hwork, 0, &i, sel, i0, hpur0, &k);
+  hLexS(hwork, i, sel, i0);
+  hMu += hZeroMult(hpur0, hwork, i, sel, i0);
+}
+
+static void hDimMult(scmon pure, int Npure, scfmon rad, int Nrad,
+ varset var, int Nvar)
+{
+  int  dn, iv, rad0, b, c, x;
+  scmon pn;
+  scfmon rn;
+  if (Nrad < 2)
+  {
+    dn = Npure + Nrad;
+    if (dn == hCo)
+    {
+      if (!Nrad)
+        hProject(pure, hsel);
+      else
+      {
+        pn = *rad;
+        for (iv = Nvar; iv; iv--)
+        {
+          x = var[iv];
+          if (pn[x])
+          {
+            pure[x] = 1;
+            hProject(pure, hsel);
+            pure[x] = 0;
+          }
+        }
+      }
+    }
+    return;
+  }
+  iv = Nvar;
+  dn = Npure+1;
+  if (dn >= hCo)
+  {
+    if (dn > hCo)
+      return;
+    loop
+    {
+      if(!pure[var[iv]])
+      {
+        if(hNotZero(rad, Nrad, var, iv))
+        {
+          pure[var[iv]] = 1;
+          hProject(pure, hsel);
+          pure[var[iv]] = 0;
+        }
+      }
+      iv--;
+      if (!iv)
+        return;
+    }
+  }
+  while(pure[var[iv]]) iv--;
+  hStepR(rad, Nrad, var, iv, &rad0);
+  iv--;
+  if (rad0 < Nrad)
+  {
+    pn = hGetpure(pure);
+    rn = hGetmem(Nrad, rad, radmem[iv]);
+    pn[var[iv + 1]] = 1;
+    hDimMult(pn, Npure + 1, rn, rad0, var, iv);
+    pn[var[iv + 1]] = 0;
+    b = rad0;
+    c = Nrad;
+    hElimR(rn, &rad0, b, c, var, iv);
+    hPure(rn, b, &c, var, iv, pn, &x);
+    hLex2R(rn, rad0, b, c, var, iv, hwork);
+    rad0 += (c - b);
+    hDimMult(pn, Npure + x, rn, rad0, var, iv);
+  }
+  else
+  {
+    hDimMult(pure, Npure, rad, Nrad, var, iv);
+  }
+}
+
+static void hDegree(ideal S, ideal Q)
+{
+  int  di;
+  int  mc;
+  hexist = hInit(S, Q, &hNexist);
+  if (!hNexist)
+  {
+    hCo = 0;
+    hMu = 1;
+    return;
+  }
+  //hWeight();
+  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
+  hsel = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
+  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  hpur0 = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  mc = hisModule;
+  hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  if (!mc)
+  {
+    memcpy(hrad, hexist, hNexist * sizeof(scmon));
+    hstc = hexist;
+    hNrad = hNstc = hNexist;
+  }
+  else
+    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  radmem = hCreate((currRing->N) - 1);
+  stcmem = hCreate((currRing->N) - 1);
+  hCo = (currRing->N) + 1;
+  di = hCo + 1;
+  loop
+  {
+    if (mc)
+    {
+      hComp(hexist, hNexist, mc, hrad, &hNrad);
+      hNstc = hNrad;
+      memcpy(hstc, hrad, hNrad * sizeof(scmon));
+    }
+    if (hNrad)
+    {
+      hNvar = (currRing->N);
+      hRadical(hrad, &hNrad, hNvar);
+      hSupp(hrad, hNrad, hvar, &hNvar);
+      if (hNvar)
+      {
+        hCo = hNvar;
+        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
+        hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
+        hLexR(hrad, hNrad, hvar, hNvar);
+        hDimSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
+      }
+    }
+    else
+    {
+      hNvar = 1;
+      hCo = 0;
+    }
+    if (hCo < di)
+    {
+      di = hCo;
+      hMu = 0;
+    }
+    if (hNvar && (hCo == di))
+    {
+      if (di && (di < (currRing->N)))
+        hDimMult(hpure, hNpure, hrad, hNrad, hvar, hNvar);
+      else if (!di)
+        hMu++;
+      else
+      {
+        hStaircase(hstc, &hNstc, hvar, hNvar);
+        if ((hNvar > 2) && (hNstc > 10))
+          hOrdSupp(hstc, hNstc, hvar, hNvar);
+        memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
+        hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
+        hLexS(hstc, hNstc, hvar, hNvar);
+        hMu += hZeroMult(hpur0, hstc, hNstc, hvar, hNvar);
+      }
+    }
+    mc--;
+    if (mc <= 0)
+      break;
+  }
+  hCo = di;
+  hKill(stcmem, (currRing->N) - 1);
+  hKill(radmem, (currRing->N) - 1);
+  omFreeSize((ADDRESS)hpur0, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  omFreeSize((ADDRESS)hsel, ((currRing->N) + 1) * sizeof(int));
+  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
+  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
+  omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
+  hDelete(hexist, hNexist);
+  if (hisModule)
+    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
+}
+
+int  scMultInt(ideal S, ideal Q)
+{
+  hDegree(S, Q);
+  return hMu;
+}
+
+void scPrintDegree(int co, int mu)
+{
+  int di = (currRing->N)-co;
+  if (currRing->OrdSgn == 1)
+  {
+    if (di>0)
+      Print("// dimension (proj.)  = %d\n// degree (proj.)   = %d\n", di-1, mu);
+    else
+      Print("// dimension (affine) = 0\n// degree (affine)  = %d\n",       mu);
+  }
+  else
+    Print("// dimension (local)   = %d\n// multiplicity = %d\n", di, mu);
+}
+
+void scDegree(ideal S, intvec *modulweight, ideal Q)
+{
+  int co, mu, l;
+  intvec *hseries2;
+  intvec *hseries1 = hFirstSeries(S, modulweight, Q);
+  l = hseries1->length()-1;
+  if (l > 1)
+    hseries2 = hSecondSeries(hseries1);
+  else
+    hseries2 = hseries1;
+  hDegreeSeries(hseries1, hseries2, &co, &mu);
+  if ((l == 1) &&(mu == 0))
+    scPrintDegree((currRing->N)+1, 0);
+  else
+    scPrintDegree(co, mu);
+  if (l>1)
+    delete hseries1;
+  delete hseries2;
+}
+
+static void hDegree0(ideal S, ideal Q)
+{
+  int  mc;
+  hexist = hInit(S, Q, &hNexist);
+  if (!hNexist)
+  {
+    hMu = -1;
+    return;
+  }
+  else
+    hMu = 0;
+  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
+  hpur0 = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  mc = hisModule;
+  if (!mc)
+  {
+    hstc = hexist;
+    hNstc = hNexist;
+  }
+  else
+    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  stcmem = hCreate((currRing->N) - 1);
+  loop
+  {
+    if (mc)
+    {
+      hComp(hexist, hNexist, mc, hstc, &hNstc);
+      if (!hNstc)
+      {
+        hMu = -1;
+        break;
+      }
+    }
+    hNvar = (currRing->N);
+    for (int i = hNvar; i; i--)
+      hvar[i] = i;
+    hStaircase(hstc, &hNstc, hvar, hNvar);
+    hSupp(hstc, hNstc, hvar, &hNvar);
+    if ((hNvar == (currRing->N)) && (hNstc >= (currRing->N)))
+    {
+      if ((hNvar > 2) && (hNstc > 10))
+        hOrdSupp(hstc, hNstc, hvar, hNvar);
+      memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
+      hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
+      if (hNpure == hNvar)
+      {
+        hLexS(hstc, hNstc, hvar, hNvar);
+        hMu += hZeroMult(hpur0, hstc, hNstc, hvar, hNvar);
+      }
+      else
+        hMu = -1;
+    }
+    else if (hNvar)
+      hMu = -1;
+    mc--;
+    if (mc <= 0 || hMu < 0)
+      break;
+  }
+  hKill(stcmem, (currRing->N) - 1);
+  omFreeSize((ADDRESS)hpur0, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
+  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
+  hDelete(hexist, hNexist);
+  if (hisModule)
+    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
+}
+
+int  scMult0Int(ideal S, ideal Q)
+{
+  hDegree0(S, Q);
+  return hMu;
+}
+
+
+// HC
+
+static poly pWork;
+
+static void hHedge(poly hEdge)
+{
+  pSetm(pWork);
+  if (pLmCmp(pWork, hEdge) == currRing->OrdSgn)
+  {
+    for (int i = hNvar; i>0; i--)
+      pSetExp(hEdge,i, pGetExp(pWork,i));
+    pSetm(hEdge);
+  }
+}
+
+
+static void hHedgeStep(scmon pure, scfmon stc,
+                       int Nstc, varset var, int Nvar,poly hEdge)
+{
+  int  iv = Nvar -1, k = var[Nvar], a, a0, a1, b, i;
+  int  x/*, x0*/;
+  scmon pn;
+  scfmon sn;
+  if (iv==0)
+  {
+    pSetExp(pWork, k, pure[k]);
+    hHedge(hEdge);
+    return;
+  }
+  else if (Nstc==0)
+  {
+    for (i = Nvar; i>0; i--)
+      pSetExp(pWork, var[i], pure[var[i]]);
+    hHedge(hEdge);
+    return;
+  }
+  x = a = 0;
+  pn = hGetpure(pure);
+  sn = hGetmem(Nstc, stc, stcmem[iv]);
+  hStepS(sn, Nstc, var, Nvar, &a, &x);
+  if (a == Nstc)
+  {
+    pSetExp(pWork, k, pure[k]);
+    hHedgeStep(pn, sn, a, var, iv,hEdge);
+    return;
+  }
+  else
+  {
+    pSetExp(pWork, k, x);
+    hHedgeStep(pn, sn, a, var, iv,hEdge);
+  }
+  b = a;
+  loop
+  {
+    a0 = a;
+    // x0 = x;
+    hStepS(sn, Nstc, var, Nvar, &a, &x);
+    hElimS(sn, &b, a0, a, var, iv);
+    a1 = a;
+    hPure(sn, a0, &a1, var, iv, pn, &i);
+    hLex2S(sn, b, a0, a1, var, iv, hwork);
+    b += (a1 - a0);
+    if (a < Nstc)
+    {
+      pSetExp(pWork, k, x);
+      hHedgeStep(pn, sn, b, var, iv,hEdge);
+    }
+    else
+    {
+      pSetExp(pWork, k, pure[k]);
+      hHedgeStep(pn, sn, b, var, iv,hEdge);
+      return;
+    }
+  }
+}
+
+void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge, ring tailRing)
+{
+  int  i;
+  int  k = ak;
+
+  #if HAVE_RINGS
+  if (rField_is_Ring(currRing) && (currRing->OrdSgn == -1))
+  {
+  //consider just monic generators (over rings with zero-divisors)
+  ideal SS=id_Copy(S,tailRing);
+  for(i=0;i<=idElem(SS);i++)
+  	{
+  	if(pIsPurePower(SS->m[i])==0)
+  		p_Delete(&SS->m[i],tailRing);
+  	}
+  	S=id_Copy(SS,tailRing);
+  }
+  #endif
+
+  hNvar = (currRing->N);
+  hexist = hInit(S, Q, &hNexist, tailRing);
+  if (k!=0)
+    hComp(hexist, hNexist, k, hexist, &hNstc);
+  else
+    hNstc = hNexist;
+  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  hvar = (varset)omAlloc((hNvar + 1) * sizeof(int));
+  hpure = (scmon)omAlloc((1 + (hNvar * hNvar)) * sizeof(int));
+  stcmem = hCreate(hNvar - 1);
+  for (i = hNvar; i>0; i--)
+    hvar[i] = i;
+  hStaircase(hexist, &hNstc, hvar, hNvar);
+  if ((hNvar > 2) && (hNstc > 10))
+    hOrdSupp(hexist, hNstc, hvar, hNvar);
+  memset(hpure, 0, (hNvar + 1) * sizeof(int));
+  hPure(hexist, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
+  hLexS(hexist, hNstc, hvar, hNvar);
+  if (hEdge!=NULL)
+    pLmFree(hEdge);
+  hEdge = pInit();
+  pWork = pInit();
+  hHedgeStep(hpure, hexist, hNstc, hvar, hNvar,hEdge);
+  pSetComp(hEdge,ak);
+  hKill(stcmem, hNvar - 1);
+  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
+  omFreeSize((ADDRESS)hvar, (hNvar + 1) * sizeof(int));
+  omFreeSize((ADDRESS)hpure, (1 + (hNvar * hNvar)) * sizeof(int));
+  hDelete(hexist, hNexist);
+  pLmFree(pWork);
+}
+
+
+
+//  kbase
+
+static poly last;
+static scmon act;
+
+static void scElKbase()
+{
+  poly q = pInit();
+  pSetCoeff0(q,nInit(1));
+  pSetExpV(q,act);
+  pNext(q) = NULL;
+  last = pNext(last) = q;
+}
+
+static int scMax( int i, scfmon stc, int Nvar)
+{
+  int x, y=stc[0][Nvar];
+  for (; i;)
+  {
+    i--;
+    x = stc[i][Nvar];
+    if (x > y) y = x;
+  }
+  return y;
+}
+
+static int scMin( int i, scfmon stc, int Nvar)
+{
+  int x, y=stc[0][Nvar];
+  for (; i;)
+  {
+    i--;
+    x = stc[i][Nvar];
+    if (x < y) y = x;
+  }
+  return y;
+}
+
+static int scRestrict( int &Nstc, scfmon stc, int Nvar)
+{
+  int x, y;
+  int i, j, Istc = Nstc;
+
+  y = MAX_INT_VAL;
+  for (i=Nstc-1; i>=0; i--)
+  {
+    j = Nvar-1;
+    loop
+    {
+      if(stc[i][j] != 0) break;
+      j--;
+      if (j == 0)
+      {
+        Istc--;
+        x = stc[i][Nvar];
+        if (x < y) y = x;
+        stc[i] = NULL;
+        break;
+      }
+    }
+  }
+  if (Istc < Nstc)
+  {
+    for (i=Nstc-1; i>=0; i--)
+    {
+      if (stc[i] && (stc[i][Nvar] >= y))
+      {
+        Istc--;
+        stc[i] = NULL;
+      }
+    }
+    j = 0;
+    while (stc[j]) j++;
+    i = j+1;
+    for(; i<Nstc; i++)
+    {
+      if (stc[i])
+      {
+        stc[j] = stc[i];
+        j++;
+      }
+    }
+    Nstc = Istc;
+    return y;
+  }
+  else
+    return -1;
+}
+
+static void scAll( int Nvar, int deg)
+{
+  int i;
+  int d = deg;
+  if (d == 0)
+  {
+    for (i=Nvar; i; i--) act[i] = 0;
+    scElKbase();
+    return;
+  }
+  if (Nvar == 1)
+  {
+    act[1] = d;
+    scElKbase();
+    return;
+  }
+  do
+  {
+    act[Nvar] = d;
+    scAll(Nvar-1, deg-d);
+    d--;
+  } while (d >= 0);
+}
+
+static void scAllKbase( int Nvar, int ideg, int deg)
+{
+  do
+  {
+    act[Nvar] = ideg;
+    scAll(Nvar-1, deg-ideg);
+    ideg--;
+  } while (ideg >= 0);
+}
+
+static void scDegKbase( scfmon stc, int Nstc, int Nvar, int deg)
+{
+  int  Ivar, Istc, i, j;
+  scfmon sn;
+  int x, ideg;
+
+  if (deg == 0)
+  {
+    for (i=Nstc-1; i>=0; i--)
+    {
+      for (j=Nvar;j;j--){ if(stc[i][j]) break; }
+      if (j==0){return;}
+    }
+    for (i=Nvar; i; i--) act[i] = 0;
+    scElKbase();
+    return;
+  }
+  if (Nvar == 1)
+  {
+    for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
+    act[1] = deg;
+    scElKbase();
+    return;
+  }
+  Ivar = Nvar-1;
+  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
+  x = scRestrict(Nstc, sn, Nvar);
+  if (x <= 0)
+  {
+    if (x == 0) return;
+    ideg = deg;
+  }
+  else
+  {
+    if (deg < x) ideg = deg;
+    else ideg = x-1;
+    if (Nstc == 0)
+    {
+      scAllKbase(Nvar, ideg, deg);
+      return;
+    }
+  }
+  loop
+  {
+    x = scMax(Nstc, sn, Nvar);
+    while (ideg >= x)
+    {
+      act[Nvar] = ideg;
+      scDegKbase(sn, Nstc, Ivar, deg-ideg);
+      ideg--;
+    }
+    if (ideg < 0) return;
+    Istc = Nstc;
+    for (i=Nstc-1; i>=0; i--)
+    {
+      if (ideg < sn[i][Nvar])
+      {
+        Istc--;
+        sn[i] = NULL;
+      }
+    }
+    if (Istc == 0)
+    {
+      scAllKbase(Nvar, ideg, deg);
+      return;
+    }
+    j = 0;
+    while (sn[j]) j++;
+    i = j+1;
+    for (; i<Nstc; i++)
+    {
+      if (sn[i])
+      {
+        sn[j] = sn[i];
+        j++;
+      }
+    }
+    Nstc = Istc;
+  }
+}
+
+static void scInKbase( scfmon stc, int Nstc, int Nvar)
+{
+  int  Ivar, Istc, i, j;
+  scfmon sn;
+  int x, ideg;
+
+  if (Nvar == 1)
+  {
+    ideg = scMin(Nstc, stc, 1);
+    while (ideg > 0)
+    {
+      ideg--;
+      act[1] = ideg;
+      scElKbase();
+    }
+    return;
+  }
+  Ivar = Nvar-1;
+  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
+  x = scRestrict(Nstc, sn, Nvar);
+  if (x == 0) return;
+  ideg = x-1;
+  loop
+  {
+    x = scMax(Nstc, sn, Nvar);
+    while (ideg >= x)
+    {
+      act[Nvar] = ideg;
+      scInKbase(sn, Nstc, Ivar);
+      ideg--;
+    }
+    if (ideg < 0) return;
+    Istc = Nstc;
+    for (i=Nstc-1; i>=0; i--)
+    {
+      if (ideg < sn[i][Nvar])
+      {
+        Istc--;
+        sn[i] = NULL;
+      }
+    }
+    j = 0;
+    while (sn[j]) j++;
+    i = j+1;
+    for (; i<Nstc; i++)
+    {
+      if (sn[i])
+      {
+        sn[j] = sn[i];
+        j++;
+      }
+    }
+    Nstc = Istc;
+  }
+}
+
+static ideal scIdKbase()
+{
+  polyset mm;
+  ideal res;
+  poly p, q = last;
+  int i = pLength(q);
+  res = idInit(i,1);
+  mm = res->m;
+  i = 0;
+  do
+  {
+    mm[i] = q;
+    i++;
+    p = pNext(q);
+    pNext(q) = NULL;
+    q = p;
+  } while (q!=NULL);
+  return res;
+}
+
+ideal scKBase(int deg, ideal s, ideal Q, intvec * mv)
+{
+  int  i, di;
+  poly p;
+
+  if (deg < 0)
+  {
+    di = scDimInt(s, Q);
+    if (di != 0)
+    {
+      //Werror("KBase not finite");
+      return idInit(1,s->rank);
+    }
+  }
+  stcmem = hCreate((currRing->N) - 1);
+  hexist = hInit(s, Q, &hNexist);
+  p = last = pInit();
+  /*pNext(p) = NULL;*/
+  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
+  *act = 0;
+  if (!hNexist)
+  {
+    scAll((currRing->N), deg);
+    goto ende;
+  }
+  if (!hisModule)
+  {
+    if (deg < 0) scInKbase(hexist, hNexist, (currRing->N));
+    else scDegKbase(hexist, hNexist, (currRing->N), deg);
+  }
+  else
+  {
+    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
+    for (i = 1; i <= hisModule; i++)
+    {
+      *act = i;
+      hComp(hexist, hNexist, i, hstc, &hNstc);
+      int deg_ei=deg;
+      if (mv!=NULL) deg_ei -= (*mv)[i-1];
+      if ((deg < 0) || (deg_ei>=0))
+      {
+        if (hNstc)
+        {
+          if (deg < 0) scInKbase(hstc, hNstc, (currRing->N));
+          else scDegKbase(hstc, hNstc, (currRing->N), deg_ei);
+        }
+        else
+          scAll((currRing->N), deg_ei);
+      }
+    }
+    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
+  }
+ende:
+  hDelete(hexist, hNexist);
+  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
+  hKill(stcmem, (currRing->N) - 1);
+  pLmDelete(&p);
+  if (p == NULL)
+    return idInit(1,s->rank);
+  else
+  {
+    last = p;
+    ideal res=scIdKbase();
+    res->rank=s->rank;
+    return res;
+  }
+}
+
+#if 0 //-- alternative implementation of scComputeHC
+void scComputeHCw(ideal ss, ideal Q, int ak, poly &hEdge, ring tailRing)
+{
+  int  i, di;
+  poly p;
+
+  if (hEdge!=NULL)
+    pLmFree(hEdge);
+
+  ideal s=idInit(IDELEMS(ss),ak);
+  for(i=IDELEMS(ss)-1;i>=0;i--)
+  {
+    if (ss->m[i]!=NULL) s->m[i]=pHead(ss->m[i]);
+  }
+  di = scDimInt(s, Q);
+  stcmem = hCreate((currRing->N) - 1);
+  hexist = hInit(s, Q, &hNexist);
+  p = last = pInit();
+  /*pNext(p) = NULL;*/
+  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
+  *act = 0;
+  if (!hNexist)
+  {
+    scAll((currRing->N), -1);
+    goto ende;
+  }
+  if (!hisModule)
+  {
+    scInKbase(hexist, hNexist, (currRing->N));
+  }
+  else
+  {
+    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
+    for (i = 1; i <= hisModule; i++)
+    {
+      *act = i;
+      hComp(hexist, hNexist, i, hstc, &hNstc);
+      if (hNstc)
+      {
+        scInKbase(hstc, hNstc, (currRing->N));
+      }
+      else
+        scAll((currRing->N), -1);
+    }
+    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
+  }
+ende:
+  hDelete(hexist, hNexist);
+  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
+  hKill(stcmem, (currRing->N) - 1);
+  pDeleteLm(&p);
+  idDelete(&s);
+  if (p == NULL)
+  {
+    return; // no HEdge
+  }
+  else
+  {
+    last = p;
+    ideal res=scIdKbase();
+    res->rank=ss->rank;
+    poly p_ind=res->m[0]; int ind=0;
+    for(i=IDELEMS(res)-1;i>0;i--)
+    {
+      if (pCmp(res->m[i],p_ind)==-1) { p_ind=res->m[i]; ind=i; }
+    }
+    assume(p_ind!=NULL);
+    assume(res->m[ind]==p_ind);
+    hEdge=p_ind;
+    res->m[ind]=NULL;
+    nDelete(&pGetCoeff(hEdge));
+    pGetCoeff(hEdge)=NULL;
+    for(i=(currRing->N);i>0;i--)
+      pIncrExp(hEdge,i);
+    pSetm(hEdge);
+
+    idDelete(&res);
+    return;
+  }
+}
+#endif
Index: kernel/combinatorics/hilb.cc
===================================================================
--- kernel/combinatorics/hilb.cc	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
+++ kernel/combinatorics/hilb.cc	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
@@ -0,0 +1,1387 @@
+/****************************************
+*  Computer Algebra System SINGULAR     *
+****************************************/
+/*
+*  ABSTRACT -  Hilbert series
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "singularconfig.h"
+#endif /* HAVE_CONFIG_H */
+#include <kernel/mod2.h>
+#include <kernel/structs.h>
+#include <kernel/febase.h>
+#include <omalloc/omalloc.h>
+#include <kernel/polys.h>
+#include <misc/intvec.h>
+#include <kernel/hutil.h>
+#include <kernel/stairc.h>
+//ADICHANGES:
+#include <kernel/ideals.h>
+#include <kernel/kstd1.h>
+#include<gmp.h>
+#include<vector>
+
+
+static int  **Qpol;
+static int  *Q0, *Ql;
+static int  hLength;
+   
+
+static int hMinModulweight(intvec *modulweight)
+{
+  int i,j,k;
+
+  if(modulweight==NULL) return 0;
+  j=(*modulweight)[0];
+  for(i=modulweight->rows()-1;i!=0;i--)
+  {
+    k=(*modulweight)[i];
+    if(k<j) j=k;
+  }
+  return j;
+}
+
+static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
+{
+  int  i, j;
+  int  x, y, z = 1;
+  int  *p;
+  for (i = Nvar; i>0; i--)
+  {
+    x = 0;
+    for (j = 0; j < Nstc; j++)
+    {
+      y = stc[j][var[i]];
+      if (y > x)
+        x = y;
+    }
+    z += x;
+    j = i - 1;
+    if (z > Ql[j])
+    {
+      if (z>(MAX_INT_VAL)/2)
+      {
+       Werror("interal arrays too big");
+       return;
+      }
+      p = (int *)omAlloc((unsigned long)z * sizeof(int));
+      if (Ql[j]!=0)
+      {
+        if (j==0)
+          memcpy(p, Qpol[j], Ql[j] * sizeof(int));
+        omFreeSize((ADDRESS)Qpol[j], Ql[j] * sizeof(int));
+      }
+      if (j==0)
+      {
+        for (x = Ql[j]; x < z; x++)
+          p[x] = 0;
+      }
+      Ql[j] = z;
+      Qpol[j] = p;
+    }
+  }
+}
+
+static int  *hAddHilb(int Nv, int x, int *pol, int *lp)
+{
+  int  l = *lp, ln, i;
+  int  *pon;
+  *lp = ln = l + x;
+  pon = Qpol[Nv];
+  memcpy(pon, pol, l * sizeof(int));
+  if (l > x)
+  {
+    for (i = x; i < l; i++)
+      pon[i] -= pol[i - x];
+    for (i = l; i < ln; i++)
+      pon[i] = -pol[i - x];
+  }
+  else
+  {
+    for (i = l; i < x; i++)
+      pon[i] = 0;
+    for (i = x; i < ln; i++)
+      pon[i] = -pol[i - x];
+  }
+  return pon;
+}
+
+static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
+{
+  int  l = lp, x, i, j;
+  int  *p, *pl;
+  p = pol;
+  for (i = Nv; i>0; i--)
+  {
+    x = pure[var[i + 1]];
+    if (x!=0)
+      p = hAddHilb(i, x, p, &l);
+  }
+  pl = *Qpol;
+  j = Q0[Nv + 1];
+  for (i = 0; i < l; i++)
+    pl[i + j] += p[i];
+  x = pure[var[1]];
+  if (x!=0)
+  {
+    j += x;
+    for (i = 0; i < l; i++)
+      pl[i + j] -= p[i];
+  }
+  j += l;
+  if (j > hLength)
+    hLength = j;
+}
+
+static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var,
+ int Nvar, int *pol, int Lpol)
+{
+  int  iv = Nvar -1, ln, a, a0, a1, b, i;
+  int  x, x0;
+  scmon pn;
+  scfmon sn;
+  int  *pon;
+  if (Nstc==0)
+  {
+    hLastHilb(pure, iv, var, pol, Lpol);
+    return;
+  }
+  x = a = 0;
+  pn = hGetpure(pure);
+  sn = hGetmem(Nstc, stc, stcmem[iv]);
+  hStepS(sn, Nstc, var, Nvar, &a, &x);
+  Q0[iv] = Q0[Nvar];
+  ln = Lpol;
+  pon = pol;
+  if (a == Nstc)
+  {
+    x = pure[var[Nvar]];
+    if (x!=0)
+      pon = hAddHilb(iv, x, pon, &ln);
+    hHilbStep(pn, sn, a, var, iv, pon, ln);
+    return;
+  }
+  else
+  {
+    pon = hAddHilb(iv, x, pon, &ln);
+    hHilbStep(pn, sn, a, var, iv, pon, ln);
+  }
+  b = a;
+  x0 = 0;
+  loop
+  {
+    Q0[iv] += (x - x0);
+    a0 = a;
+    x0 = x;
+    hStepS(sn, Nstc, var, Nvar, &a, &x);
+    hElimS(sn, &b, a0, a, var, iv);
+    a1 = a;
+    hPure(sn, a0, &a1, var, iv, pn, &i);
+    hLex2S(sn, b, a0, a1, var, iv, hwork);
+    b += (a1 - a0);
+    ln = Lpol;
+    if (a < Nstc)
+    {
+      pon = hAddHilb(iv, x - x0, pol, &ln);
+      hHilbStep(pn, sn, b, var, iv, pon, ln);
+    }
+    else
+    {
+      x = pure[var[Nvar]];
+      if (x!=0)
+        pon = hAddHilb(iv, x - x0, pol, &ln);
+      else
+        pon = pol;
+      hHilbStep(pn, sn, b, var, iv, pon, ln);
+      return;
+    }
+  }
+}
+
+/*
+*basic routines
+*/
+static void hWDegree(intvec *wdegree)
+{
+  int i, k;
+  int x;
+
+  for (i=(currRing->N); i; i--)
+  {
+    x = (*wdegree)[i-1];
+    if (x != 1)
+    {
+      for (k=hNexist-1; k>=0; k--)
+      {
+        hexist[k][i] *= x;
+      }
+    }
+  }
+}
+// ---------------------------------- ADICHANGES ---------------------------------------------
+//!!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+//returns the degree of the monomial
+static int DegMon(poly p)
+{
+    #if 1
+    int i,deg;
+    deg = 0;
+    for(i=1;i<=currRing->N;i++)
+    {
+        deg = deg + p_GetExp(p, i, currRing);
+    }
+    return(deg);
+    #else
+    return(p_Deg(p, currRing));
+    #endif
+}
+
+//Tests if the ideal is sorted by degree
+static bool idDegSortTest(ideal I)
+{
+    if((I == NULL)||(idIs0(I)))
+    {
+        return(TRUE);
+    }
+    for(int i = 0; i<IDELEMS(I)-1; i++)
+    {
+        if(DegMon(I->m[i])>DegMon(I->m[i+1]))
+        {
+            idPrint(I);
+            Werror("Ideal is not deg sorted!!");
+            return(FALSE);
+        }
+    }
+    return(TRUE);
+}
+
+//adds the new polynomial at the coresponding position
+//and simplifies the ideal
+static ideal SortByDeg_p(ideal I, poly p)
+{
+    int i,j;
+    if((I == NULL) || (idIs0(I)))
+    {
+        ideal res = idInit(1,1);
+        res->m[0] = p;
+        return(res);
+    }
+    idSkipZeroes(I);
+    #if 1
+    for(i = 0; (i<IDELEMS(I)) && (DegMon(I->m[i])<=DegMon(p)); i++)
+    {
+        if(p_DivisibleBy( I->m[i],p, currRing))
+        {
+            return(I);
+        }
+    }
+    for(i = IDELEMS(I)-1; (i>=0) && (DegMon(I->m[i])>=DegMon(p)); i--)
+    {
+        if(p_DivisibleBy(p,I->m[i], currRing))
+        {
+            I->m[i] = NULL;
+        }
+    }
+    if(idIs0(I))
+    {
+        idSkipZeroes(I);
+        I->m[0] = p;
+        return(I);
+    }
+    #endif
+    idSkipZeroes(I);
+    //First I take the case when all generators have the same degree
+    if(DegMon(I->m[0]) == DegMon(I->m[IDELEMS(I)-1]))
+    {
+        if(DegMon(p)<DegMon(I->m[0]))
+        {
+            idInsertPoly(I,p);
+            idSkipZeroes(I);
+            for(i=IDELEMS(I)-1;i>=1; i--)
+            {
+                I->m[i] = I->m[i-1];
+            }
+            I->m[0] = p;
+            return(I);
+        }
+        if(DegMon(p)>=DegMon(I->m[IDELEMS(I)-1]))
+        {
+            idInsertPoly(I,p);
+            idSkipZeroes(I);
+            return(I);
+        }
+    }
+    if(DegMon(p)<=DegMon(I->m[0]))
+    {
+        idInsertPoly(I,p);
+        idSkipZeroes(I);
+        for(i=IDELEMS(I)-1;i>=1; i--)
+        {
+            I->m[i] = I->m[i-1];
+        }
+        I->m[0] = p;
+        return(I);
+    }
+    if(DegMon(p)>=DegMon(I->m[IDELEMS(I)-1]))
+    {
+        idInsertPoly(I,p);
+        idSkipZeroes(I);
+        return(I);
+    }
+    for(i = IDELEMS(I)-2; ;)
+    {
+        if(DegMon(p)==DegMon(I->m[i]))
+        {
+            idInsertPoly(I,p);
+            idSkipZeroes(I);
+            for(j = IDELEMS(I)-1; j>=i+1;j--)
+            {
+                I->m[j] = I->m[j-1];
+            }
+            I->m[i] = p;
+            return(I);
+        }
+        if(DegMon(p)>DegMon(I->m[i]))
+        {
+            idInsertPoly(I,p);
+            idSkipZeroes(I);
+            for(j = IDELEMS(I)-1; j>=i+2;j--)
+            {
+                I->m[j] = I->m[j-1];
+            }
+            I->m[i+1] = p;
+            return(I);
+        }
+        i--;
+    }
+}
+
+//it sorts the ideal by the degrees
+static ideal SortByDeg(ideal I)
+{
+    if(idIs0(I))
+    {
+        return(I);
+    }
+    idSkipZeroes(I);
+    int i;
+    ideal res;
+    idSkipZeroes(I);
+    res = idInit(1,1);
+    res->m[0] = poly(0);
+    for(i = 0; i<=IDELEMS(I)-1;i++)
+    {
+        res = SortByDeg_p(res, I->m[i]);
+    }
+    idSkipZeroes(res);
+    //idDegSortTest(res);
+    return(res);
+}
+
+//idQuot(I,p) for I monomial ideal, p a ideal with a single monomial.
+ideal idQuotMon(ideal Iorig, ideal p)
+{
+    if(idIs0(Iorig))
+    {
+        ideal res = idInit(1,1);
+        res->m[0] = poly(0);
+        return(res);
+    }
+    if(idIs0(p))
+    {
+        ideal res = idInit(1,1);
+        res->m[0] = pOne();
+        return(res);
+    }
+    ideal I = idCopy(Iorig);
+    ideal res = idInit(IDELEMS(I),1);
+    int i,j;
+    int dummy;
+    for(i = 0; i<IDELEMS(I); i++)
+    {
+        res->m[i] = p_Copy(I->m[i], currRing);
+        for(j = 1; (j<=currRing->N) ; j++)
+        {
+            dummy = p_GetExp(p->m[0], j, currRing);
+            if(dummy > 0)
+            {
+                if(p_GetExp(I->m[i], j, currRing) < dummy)
+                {
+                    p_SetExp(res->m[i], j, 0, currRing);
+                }
+                else
+                {
+                    p_SetExp(res->m[i], j, p_GetExp(I->m[i], j, currRing) - dummy, currRing);
+                }
+            }
+        }
+        p_Setm(res->m[i], currRing);
+        if(DegMon(res->m[i]) == DegMon(I->m[i]))
+        {
+            res->m[i] = NULL;
+        }
+        else
+        {
+            I->m[i] = NULL;
+        }
+    }
+    idSkipZeroes(res);
+    idSkipZeroes(I);
+    if(!idIs0(res))
+    {
+    for(i = 0; i<=IDELEMS(res)-1; i++)
+    {
+        I = SortByDeg_p(I,res->m[i]);
+    }
+    }
+    //idDegSortTest(I);
+    return(I);
+}
+
+//id_Add for monomials
+static ideal idAddMon(ideal I, ideal p)
+{
+    #if 1
+    I = SortByDeg_p(I,p->m[0]);
+    #else
+    I = id_Add(I,p,currRing);
+    #endif
+    //idSkipZeroes(I);
+    return(I);
+}
+
+//searches for a variable that is not yet used (assumes that the ideal is sqrfree)
+static poly ChoosePVar (ideal I)
+{
+    bool flag=TRUE;
+    int i,j;
+    poly res;
+    for(i=1;i<=currRing->N;i++)
+    {
+        flag=TRUE;
+        for(j=IDELEMS(I)-1;(j>=0)&&(flag);j--)
+        {
+            if(p_GetExp(I->m[j], i, currRing)>0)
+            {
+                flag=FALSE;
+            }
+        }
+        
+        if(flag == TRUE)
+        {
+            res = p_ISet(1, currRing);
+            p_SetExp(res, i, 1, currRing);
+            p_Setm(res,currRing);
+            return(res);
+        }
+        else
+        {
+            p_Delete(&res, currRing);
+        }
+    }
+    return(NULL); //i.e. it is the maximal ideal
+}
+
+//choice XL: last entry divided by x (xy10z15 -> y9z14)
+static poly ChoosePXL(ideal I)
+{
+    int i,j,dummy=0;
+    poly m;
+    for(i = IDELEMS(I)-1; (i>=0) && (dummy == 0); i--)
+    {
+        for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
+        {
+            if(p_GetExp(I->m[i],j, currRing)>1)
+            {
+                dummy = 1;
+            }
+        }
+    }
+    m = p_Copy(I->m[i+1],currRing);
+    for(j = 1; j<=currRing->N; j++)
+    {
+        dummy = p_GetExp(m,j,currRing);
+        if(dummy >= 1)
+        {
+            p_SetExp(m, j, dummy-1, currRing);
+        }
+    }
+    if(!p_IsOne(m, currRing))
+    {
+        p_Setm(m, currRing);
+        return(m);
+    }
+    m = ChoosePVar(I);
+    return(m);
+}
+
+//choice XF: first entry divided by x (xy10z15 -> y9z14)
+static poly ChoosePXF(ideal I)
+{
+    int i,j,dummy=0;
+    poly m;
+    for(i =0 ; (i<=IDELEMS(I)-1) && (dummy == 0); i++)
+    {
+        for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
+        {
+            if(p_GetExp(I->m[i],j, currRing)>1)
+            {
+                dummy = 1;
+            }
+        }
+    }
+    m = p_Copy(I->m[i-1],currRing);
+    for(j = 1; j<=currRing->N; j++)
+    {
+        dummy = p_GetExp(m,j,currRing);
+        if(dummy >= 1)
+        {
+            p_SetExp(m, j, dummy-1, currRing);
+        }
+    }
+    if(!p_IsOne(m, currRing))
+    {
+        p_Setm(m, currRing);
+        return(m);
+    }
+    m = ChoosePVar(I);
+    return(m);
+}
+
+//choice OL: last entry the first power (xy10z15 -> xy9z15)
+static poly ChoosePOL(ideal I)
+{
+    int i,j,dummy;
+    poly m;
+    for(i = IDELEMS(I)-1;i>=0;i--)
+    {
+        m = p_Copy(I->m[i],currRing);
+        for(j=1;j<=currRing->N;j++)
+        {
+            dummy = p_GetExp(m,j,currRing);
+            if(dummy > 0)
+            {
+                p_SetExp(m,j,dummy-1,currRing);
+                p_Setm(m,currRing);
+            }       
+        }
+        if(!p_IsOne(m, currRing))
+        {
+            return(m);
+        }
+        else
+        {
+            p_Delete(&m,currRing);
+        }
+    }
+    m = ChoosePVar(I);
+    return(m);
+}
+
+//choice OF: first entry the first power (xy10z15 -> xy9z15)
+static poly ChoosePOF(ideal I)
+{
+    int i,j,dummy;
+    poly m;
+    for(i = 0 ;i<=IDELEMS(I)-1;i++)
+    {
+        m = p_Copy(I->m[i],currRing);
+        for(j=1;j<=currRing->N;j++)
+        {
+            dummy = p_GetExp(m,j,currRing);
+            if(dummy > 0)
+            {
+                p_SetExp(m,j,dummy-1,currRing);
+                p_Setm(m,currRing);
+            }       
+        }
+        if(!p_IsOne(m, currRing))
+        {
+            return(m);
+        }
+        else
+        {
+            p_Delete(&m,currRing);
+        }
+    }
+    m = ChoosePVar(I);
+    return(m);
+}
+
+//choice VL: last entry the first variable with power (xy10z15 -> y)
+static poly ChoosePVL(ideal I)
+{
+    int i,j,dummy;
+    bool flag = TRUE;
+    poly m = p_ISet(1,currRing);
+    for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
+    {
+        flag = TRUE;
+        for(j=1;(j<=currRing->N) && (flag);j++)
+        {
+            dummy = p_GetExp(I->m[i],j,currRing);
+            if(dummy >= 2)
+            {
+                p_SetExp(m,j,1,currRing);
+                p_Setm(m,currRing);
+                flag = FALSE;
+            }       
+        }
+        if(!p_IsOne(m, currRing))
+        {
+            return(m);
+        }
+    }
+    m = ChoosePVar(I);
+    return(m);
+}
+
+//choice VF: first entry the first variable with power (xy10z15 -> y)
+static poly ChoosePVF(ideal I)
+{
+    int i,j,dummy;
+    bool flag = TRUE;
+    poly m = p_ISet(1,currRing);
+    for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
+    {
+        flag = TRUE;
+        for(j=1;(j<=currRing->N) && (flag);j++)
+        {
+            dummy = p_GetExp(I->m[i],j,currRing);
+            if(dummy >= 2)
+            {
+                p_SetExp(m,j,1,currRing);
+                p_Setm(m,currRing);
+                flag = FALSE;
+            }       
+        }
+        if(!p_IsOne(m, currRing))
+        {
+            return(m);
+        }
+    }
+    m = ChoosePVar(I);
+    return(m);
+}
+
+//choice JL: last entry just variable with power (xy10z15 -> y10)
+static poly ChoosePJL(ideal I)
+{
+    int i,j,dummy;
+    bool flag = TRUE;
+    poly m = p_ISet(1,currRing);
+    for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
+    {
+        flag = TRUE;
+        for(j=1;(j<=currRing->N) && (flag);j++)
+        {
+            dummy = p_GetExp(I->m[i],j,currRing);
+            if(dummy >= 2)
+            {
+                p_SetExp(m,j,dummy-1,currRing);
+                p_Setm(m,currRing);
+                flag = FALSE;
+            }       
+        }
+        if(!p_IsOne(m, currRing))
+        {
+            return(m);
+        }
+    }
+    m = ChoosePVar(I);
+    return(m);
+}
+
+//choice JF: last entry just variable with power -1 (xy10z15 -> y9)
+static poly ChoosePJF(ideal I)
+{
+    int i,j,dummy;
+    bool flag = TRUE;
+    poly m = p_ISet(1,currRing);
+    for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
+    {
+        flag = TRUE;
+        for(j=1;(j<=currRing->N) && (flag);j++)
+        {
+            dummy = p_GetExp(I->m[i],j,currRing);
+            if(dummy >= 2)
+            {
+                p_SetExp(m,j,dummy-1,currRing);
+                p_Setm(m,currRing);
+                flag = FALSE;
+            }       
+        }
+        if(!p_IsOne(m, currRing))
+        {
+            return(m);
+        }
+    }
+    m = ChoosePVar(I);
+    return(m);
+}
+
+//chooses 1 \neq p \not\in S. This choice should be made optimal
+static poly ChooseP(ideal I)
+{
+    poly m;
+    //  TEST TO SEE WHICH ONE IS BETTER
+    //m = ChoosePXL(I);
+    //m = ChoosePXF(I);
+    //m = ChoosePOL(I);
+    //m = ChoosePOF(I);
+    //m = ChoosePVL(I);
+    //m = ChoosePVF(I);
+    m = ChoosePJL(I);
+    //m = ChoosePJF(I);
+    return(m);
+}
+
+//searches for a monomial of degree d>=2 and divides it by a variable (result monomial of deg d-1)
+static poly SearchP(ideal I)
+{
+    int i,j,exp;
+    poly res;
+    if(DegMon(I->m[IDELEMS(I)-1])<=1)
+    {
+        res = ChoosePVar(I);
+        return(res);
+    }
+    i = IDELEMS(I)-1;
+    res = p_Copy(I->m[i], currRing);
+    for(j=1;j<=currRing->N;j++)
+    {
+        exp = p_GetExp(I->m[i], j, currRing);
+        if(exp > 0)
+        {
+            p_SetExp(res, j, exp - 1, currRing);
+            p_Setm(res,currRing);
+            return(res);
+        }
+    }
+}
+
+//test if the ideal is of the form (x1, ..., xr)
+static bool JustVar(ideal I)
+{
+    #if 0
+    int i,j;
+    bool foundone;
+    for(i=0;i<=IDELEMS(I)-1;i++)
+    {
+        foundone = FALSE;
+        for(j = 1;j<=currRing->N;j++)
+        {
+            if(p_GetExp(I->m[i], j, currRing)>0)
+            {
+                if(foundone == TRUE)
+                {
+                    return(FALSE);
+                }
+                foundone = TRUE;
+            }
+        }        
+    }
+    return(TRUE);
+    #else
+    if(DegMon(I->m[IDELEMS(I)-1])>1)
+    {
+        return(FALSE);
+    }
+    return(TRUE);
+    #endif
+}
+
+//computes the Euler Characteristic of the ideal
+static void eulerchar (ideal I, int variables, mpz_ptr ec)
+{
+  loop
+  {
+    mpz_t dummy;
+    if(JustVar(I) == TRUE)
+    {
+        if(IDELEMS(I) == variables)
+        {
+            mpz_init(dummy);
+            if((variables % 2) == 0)
+                {mpz_set_si(dummy, 1);}
+            else
+                {mpz_set_si(dummy, -1);}
+            mpz_add(ec, ec, dummy);
+        }
+        //mpz_clear(dummy);
+        return;        
+    }
+    ideal p = idInit(1,1);
+    p->m[0] = SearchP(I);
+    //idPrint(I);
+    //idPrint(p);
+    //printf("\nNow get in idQuotMon\n");
+    ideal Ip = idQuotMon(I,p);
+    //idPrint(Ip);
+    //Ip = SortByDeg(Ip);
+    int i,howmanyvarinp = 0;
+    for(i = 1;i<=currRing->N;i++)
+    {
+        if(p_GetExp(p->m[0],i,currRing)>0)
+        {
+            howmanyvarinp++;
+        }
+    }    
+    eulerchar(Ip, variables-howmanyvarinp, ec);
+    id_Delete(&Ip, currRing);
+    I = idAddMon(I,p);
+  }
+}
+
+//tests if an ideal is Square Free, if no, returns the variable which appears at powers >1
+static poly SqFree (ideal I)
+{
+    int i,j;
+    bool flag=TRUE;
+    poly notsqrfree = NULL;
+    if(DegMon(I->m[IDELEMS(I)-1])<=1)
+    {
+        return(notsqrfree);
+    }
+    for(i=IDELEMS(I)-1;(i>=0)&&(flag);i--)
+    {
+        for(j=1;(j<=currRing->N)&&(flag);j++)
+        {
+            if(p_GetExp(I->m[i],j,currRing)>1)
+            {
+                flag=FALSE;
+                notsqrfree = p_ISet(1,currRing);
+                p_SetExp(notsqrfree,j,1,currRing);
+            }
+        }
+    }
+    if(notsqrfree != NULL)
+    {
+        p_Setm(notsqrfree,currRing);
+    }
+    return(notsqrfree);
+}
+
+//checks if a polynomial is in I
+static bool IsIn(poly p, ideal I)
+{
+    //assumes that I is ordered by degree
+    if(idIs0(I))
+    {
+        if(p==poly(0))
+        {
+            return(TRUE);
+        }
+        else
+        {
+            return(FALSE);
+        }
+    }
+    if(p==poly(0))
+    {
+        return(FALSE);
+    }
+    int i,j;
+    bool flag;
+    for(i = 0;i<IDELEMS(I);i++)
+    {
+        flag = TRUE;
+        for(j = 1;(j<=currRing->N) &&(flag);j++)
+        {
+            if(p_GetExp(p, j, currRing)<p_GetExp(I->m[i], j, currRing))
+            {
+                flag = FALSE;
+            }
+        }
+        if(flag)
+        {
+            return(TRUE);
+        }
+    }
+    return(FALSE);
+}
+
+//computes the lcm of min I, I monomial ideal
+static poly LCMmon(ideal I)
+{
+    if(idIs0(I))
+    {
+        return(NULL);
+    }
+    poly m;
+    int dummy,i,j;
+    m = p_ISet(1,currRing);
+    for(i=1;i<=currRing->N;i++)
+    {
+        dummy=0;
+        for(j=IDELEMS(I)-1;j>=0;j--)
+        {
+            if(p_GetExp(I->m[j],i,currRing) > dummy)
+            {
+                dummy = p_GetExp(I->m[j],i,currRing);
+            }
+        }
+        p_SetExp(m,i,dummy,currRing);
+    }
+    p_Setm(m,currRing);
+    return(m);
+}
+
+//the Roune Slice Algorithm
+void rouneslice(ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int* &hilbpower)
+{
+  loop
+  {
+    (steps)++;
+    int i,j;
+    int dummy;
+    poly m;
+    ideal p, koszsimp;
+    //----------- PRUNING OF S ---------------
+    //S SHOULD IN THIS POINT BE ORDERED BY DEGREE
+    for(i=IDELEMS(S)-1;i>=0;i--)
+    {
+        if(IsIn(S->m[i],I))
+        {
+            S->m[i]=NULL;
+            prune++;
+        }
+    }
+    idSkipZeroes(S);
+    //----------------------------------------
+    for(i=IDELEMS(I)-1;i>=0;i--)
+    {
+        m = p_Copy(I->m[i],currRing);
+        for(j=1;j<=currRing->N;j++)
+        {
+            dummy = p_GetExp(m,j,currRing);
+            if(dummy > 0)
+            {
+                p_SetExp(m,j,dummy-1,currRing);
+            }       
+        }
+        p_Setm(m, currRing);
+        if(IsIn(m,S))
+        {
+            I->m[i]=NULL;
+            //printf("\n Deleted, since pi(m) is in S\n");pWrite(m);
+        }
+    }
+    idSkipZeroes(I);
+    //----------- MORE PRUNING OF S ------------
+    m = LCMmon(I); 
+    if(m != NULL)
+    {
+        for(i=0;i<IDELEMS(S);i++)
+        {      
+            if(!(p_DivisibleBy(S->m[i], m, currRing)))  
+            {
+                S->m[i] = NULL;
+                j++;
+                moreprune++;
+            }
+            else
+            {
+                if(pLmEqual(S->m[i],m))
+                {
+                    S->m[i] = NULL;
+                    moreprune++;
+                }
+            }
+        }
+    idSkipZeroes(S);
+    }
+    /*printf("\n---------------------------\n");
+    printf("\n      I\n");idPrint(I);
+    printf("\n      S\n");idPrint(S);
+    printf("\n      q\n");pWrite(q);
+    getchar();*/
+    
+    if(idIs0(I))
+    {
+        id_Delete(&I, currRing);
+        id_Delete(&S, currRing);
+        p_Delete(&m, currRing);
+        break;
+    }
+    m = LCMmon(I);
+    if(!p_DivisibleBy(x,m, currRing))
+    {
+        //printf("\nx does not divide lcm(I)");
+        //printf("\nEmpty set");pWrite(q);
+        id_Delete(&I, currRing);
+        id_Delete(&S, currRing);
+        p_Delete(&m, currRing);
+        break;
+    }
+    m = SqFree(I);
+    if(m==NULL)
+    {
+        //printf("\n      Corner: ");
+        //pWrite(q);
+        //printf("\n      With the facets of the dual simplex:\n");
+        //idPrint(I);
+        mpz_t ec;
+        mpz_init(ec);
+        mpz_ptr ec_ptr = ec;
+        eulerchar(I, currRing->N, ec_ptr);
+        bool flag = FALSE;
+        if(NNN==0)
+            {
+                hilbertcoef = (mpz_ptr)omAlloc((NNN+1)*sizeof(mpz_t));
+                hilbpower = (int*)omAlloc((NNN+1)*sizeof(int));
+                mpz_init( &hilbertcoef[NNN]);
+                mpz_set(  &hilbertcoef[NNN], ec);
+                mpz_clear(ec);
+                hilbpower[NNN] = DegMon(q);
+                NNN++;
+            }
+        else
+        {
+            //I look if the power appears already
+            for(i = 0;(i<NNN)&&(flag == FALSE)&&(DegMon(q)>=hilbpower[i]);i++)
+            {
+                if((hilbpower[i]) == (DegMon(q)))
+                {
+                    flag = TRUE;
+                    mpz_add(&hilbertcoef[i],&hilbertcoef[i],ec_ptr);
+                }
+            }
+            if(flag == FALSE)
+            {
+                hilbertcoef = (mpz_ptr)omRealloc(hilbertcoef, (NNN+1)*sizeof(mpz_t));
+                hilbpower = (int*)omRealloc(hilbpower, (NNN+1)*sizeof(int));
+                mpz_init(&hilbertcoef[NNN]);
+                for(j = NNN; j>i; j--)
+                {
+                    mpz_set(&hilbertcoef[j],&hilbertcoef[j-1]);
+                    hilbpower[j] = hilbpower[j-1];
+                }
+                mpz_set(  &hilbertcoef[i], ec);
+                mpz_clear(ec);
+                hilbpower[i] = DegMon(q);
+                NNN++;
+            }
+        }
+        break;
+    }
+    m = ChooseP(I);
+    p = idInit(1,1);
+    p->m[0] = m;
+    ideal Ip = idQuotMon(I,p);
+    ideal Sp = idQuotMon(S,p);
+    poly pq = pp_Mult_mm(q,m,currRing);
+    rouneslice(Ip, Sp, pq, x, prune, moreprune, steps, NNN, hilbertcoef,hilbpower);
+    //id_Delete(&Ip, currRing);
+    //id_Delete(&Sp, currRing);
+    S = idAddMon(S,p);
+    p->m[0]=NULL; 
+    id_Delete(&p, currRing); // p->m[0] was also in S
+    p_Delete(&pq,currRing);
+  }
+}
+
+//it computes the first hilbert series by means of Roune Slice Algorithm
+void slicehilb(ideal I)
+{
+    //printf("Adi changes are here: \n");
+    int i, NNN = 0;
+    int steps = 0, prune = 0, moreprune = 0;
+    mpz_ptr hilbertcoef;
+    int *hilbpower;
+    ideal S = idInit(1,1);
+    poly q = p_ISet(1,currRing);
+    ideal X = idInit(1,1);
+    X->m[0]=p_One(currRing);
+    for(i=1;i<=currRing->N;i++)
+    {
+            p_SetExp(X->m[0],i,1,currRing);   
+    }
+    p_Setm(X->m[0],currRing);
+    I = id_Mult(I,X,currRing);
+    I = SortByDeg(I);
+    //printf("\n-------------RouneSlice--------------\n");
+    rouneslice(I,S,q,X->m[0],prune, moreprune, steps, NNN, hilbertcoef, hilbpower);
+    //printf("\nIn total Prune got rid of %i elements\n",prune);
+    //printf("\nIn total More Prune got rid of %i elements\n",moreprune);
+    //printf("\nSteps of rouneslice: %i\n\n", steps);
+    mpz_t coefhilb;
+    mpz_t dummy;
+    mpz_init(coefhilb);
+    mpz_init(dummy);
+    printf("\n//  %8d t^0",1);
+    for(i = 0; i<NNN; i++)
+    {
+        if(mpz_sgn(&hilbertcoef[i])!=0)
+        {
+            gmp_printf("\n//  %8Zd t^%d",&hilbertcoef[i],hilbpower[i]);
+        }
+    }
+    omFreeSize(hilbertcoef, (NNN)*sizeof(mpz_t));
+    omFreeSize(hilbpower, (NNN)*sizeof(int));
+    //printf("\n-------------------------------------\n");
+}
+
+// -------------------------------- END OF CHANGES -------------------------------------------
+static intvec * hSeries(ideal S, intvec *modulweight,
+                int /*notstc*/, intvec *wdegree, ideal Q, ring tailRing)
+{
+  intvec *work, *hseries1=NULL;
+  int  mc;
+  int  p0;
+  int  i, j, k, l, ii, mw;
+  hexist = hInit(S, Q, &hNexist, tailRing);
+  if (hNexist==0)
+  {
+    hseries1=new intvec(2);
+    (*hseries1)[0]=1;
+    (*hseries1)[1]=0;
+    return hseries1;
+  }
+
+  #if 0
+  if (wdegree == NULL)
+    hWeight();
+  else
+    hWDegree(wdegree);
+  #else
+  if (wdegree != NULL) hWDegree(wdegree);
+  #endif
+
+  p0 = 1;
+  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
+  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  stcmem = hCreate((currRing->N) - 1);
+  Qpol = (int **)omAlloc(((currRing->N) + 1) * sizeof(int *));
+  Ql = (int *)omAlloc0(((currRing->N) + 1) * sizeof(int));
+  Q0 = (int *)omAlloc(((currRing->N) + 1) * sizeof(int));
+  *Qpol = NULL;
+  hLength = k = j = 0;
+  mc = hisModule;
+  if (mc!=0)
+  {
+    mw = hMinModulweight(modulweight);
+    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
+  }
+  else
+  {
+    mw = 0;
+    hstc = hexist;
+    hNstc = hNexist;
+  }
+  loop
+  {
+    if (mc!=0)
+    {
+      hComp(hexist, hNexist, mc, hstc, &hNstc);
+      if (modulweight != NULL)
+        j = (*modulweight)[mc-1]-mw;
+    }
+    if (hNstc!=0)
+    {
+      hNvar = (currRing->N);
+      for (i = hNvar; i>=0; i--)
+        hvar[i] = i;
+      //if (notstc) // TODO: no mon divides another
+        hStaircase(hstc, &hNstc, hvar, hNvar);
+      hSupp(hstc, hNstc, hvar, &hNvar);
+      if (hNvar!=0)
+      {
+        if ((hNvar > 2) && (hNstc > 10))
+          hOrdSupp(hstc, hNstc, hvar, hNvar);
+        hHilbEst(hstc, hNstc, hvar, hNvar);
+        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
+        hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
+        hLexS(hstc, hNstc, hvar, hNvar);
+        Q0[hNvar] = 0;
+        hHilbStep(hpure, hstc, hNstc, hvar, hNvar, &p0, 1);
+      }
+    }
+    else
+    {
+      if(*Qpol!=NULL)
+        (**Qpol)++;
+      else
+      {
+        *Qpol = (int *)omAlloc(sizeof(int));
+        hLength = *Ql = **Qpol = 1;
+      }
+    }
+    if (*Qpol!=NULL)
+    {
+      i = hLength;
+      while ((i > 0) && ((*Qpol)[i - 1] == 0))
+        i--;
+      if (i > 0)
+      {
+        l = i + j;
+        if (l > k)
+        {
+          work = new intvec(l);
+          for (ii=0; ii<k; ii++)
+            (*work)[ii] = (*hseries1)[ii];
+          if (hseries1 != NULL)
+            delete hseries1;
+          hseries1 = work;
+          k = l;
+        }
+        while (i > 0)
+        {
+          (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
+          (*Qpol)[i - 1] = 0;
+          i--;
+        }
+      }
+    }
+    mc--;
+    if (mc <= 0)
+      break;
+  }
+  if (k==0)
+  {
+    hseries1=new intvec(2);
+    (*hseries1)[0]=0;
+    (*hseries1)[1]=0;
+  }
+  else
+  {
+    l = k+1;
+    while ((*hseries1)[l-2]==0) l--;
+    if (l!=k)
+    {
+      work = new intvec(l);
+      for (ii=l-2; ii>=0; ii--)
+        (*work)[ii] = (*hseries1)[ii];
+      delete hseries1;
+      hseries1 = work;
+    }
+    (*hseries1)[l-1] = mw;
+  }
+  for (i = 0; i <= (currRing->N); i++)
+  {
+    if (Ql[i]!=0)
+      omFreeSize((ADDRESS)Qpol[i], Ql[i] * sizeof(int));
+  }
+  omFreeSize((ADDRESS)Q0, ((currRing->N) + 1) * sizeof(int));
+  omFreeSize((ADDRESS)Ql, ((currRing->N) + 1) * sizeof(int));
+  omFreeSize((ADDRESS)Qpol, ((currRing->N) + 1) * sizeof(int *));
+  hKill(stcmem, (currRing->N) - 1);
+  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
+  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
+  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
+  hDelete(hexist, hNexist);
+  if (hisModule!=0)
+    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
+  return hseries1;
+}
+
+
+intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
+{
+  return hSeries(S, modulweight, 0, wdegree, Q, tailRing);
+}
+
+intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
+{
+  return hSeries(S, modulweight, 1, wdegree, Q, tailRing);
+}
+
+intvec * hSecondSeries(intvec *hseries1)
+{
+  intvec *work, *hseries2;
+  int i, j, k, s, t, l;
+  if (hseries1 == NULL)
+    return NULL;
+  work = new intvec(hseries1);
+  k = l = work->length()-1;
+  s = 0;
+  for (i = k-1; i >= 0; i--)
+    s += (*work)[i];
+  loop
+  {
+    if ((s != 0) || (k == 1))
+      break;
+    s = 0;
+    t = (*work)[k-1];
+    k--;
+    for (i = k-1; i >= 0; i--)
+    {
+      j = (*work)[i];
+      (*work)[i] = -t;
+      s += t;
+      t += j;
+    }
+  }
+  hseries2 = new intvec(k+1);
+  for (i = k-1; i >= 0; i--)
+    (*hseries2)[i] = (*work)[i];
+  (*hseries2)[k] = (*work)[l];
+  delete work;
+  return hseries2;
+}
+
+void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
+{
+  int m, i, j, k;
+  *co = *mu = 0;
+  if ((s1 == NULL) || (s2 == NULL))
+    return;
+  i = s1->length();
+  j = s2->length();
+  if (j > i)
+    return;
+  m = 0;
+  for(k=j-2; k>=0; k--)
+    m += (*s2)[k];
+  *mu = m;
+  *co = i - j;
+}
+
+static void hPrintHilb(intvec *hseries)
+{
+  int  i, j, l, k;
+  if (hseries == NULL)
+    return;
+  l = hseries->length()-1;
+  k = (*hseries)[l];
+  for (i = 0; i < l; i++)
+  {
+    j = (*hseries)[i];
+    if (j != 0)
+    {
+      Print("//  %8d t^%d\n", j, i+k);
+    }
+  }
+}
+
+/*
+*caller
+*/
+void hLookSeries(ideal S, intvec *modulweight, ideal Q)
+{
+  int co, mu, l;
+  intvec *hseries2;
+  intvec *hseries1 = hFirstSeries(S, modulweight, Q);
+  hPrintHilb(hseries1);
+  l = hseries1->length()-1;
+  if (l > 1)
+    hseries2 = hSecondSeries(hseries1);
+  else
+    hseries2 = hseries1;
+  hDegreeSeries(hseries1, hseries2, &co, &mu);
+  PrintLn();
+  hPrintHilb(hseries2);
+  if ((l == 1) &&(mu == 0))
+    scPrintDegree((currRing->N)+1, 0);
+  else
+    scPrintDegree(co, mu);
+  if (l>1)
+    delete hseries1;
+  delete hseries2;
+}
+
+
+
Index: kernel/combinatorics/hutil.cc
===================================================================
--- kernel/combinatorics/hutil.cc	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
+++ kernel/combinatorics/hutil.cc	(revision f2f5ea2fec5d8f364d1161f380e4aa9290565847)
@@ -0,0 +1,1062 @@
+/****************************************
+*  Computer Algebra System SINGULAR     *
+****************************************/
+/*
+* ABSTRACT: Utilities for staircase operations
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "singularconfig.h"
+#endif /* HAVE_CONFIG_H */
+#include <kernel/mod2.h>
+#include <kernel/structs.h>
+#include <kernel/febase.h>
+#include <omalloc/omalloc.h>
+#include <kernel/ideals.h>
+#include <kernel/polys.h>
+#include <kernel/hutil.h>
+
+scfmon hexist, hstc, hrad, hwork;
+scmon hpure, hpur0;
+varset hvar, hsel;
+int  hNexist, hNstc, hNrad, hNvar, hNpure;
+int hisModule;
+monf stcmem, radmem;
+
+// Making a global "security" copy of the allocated exponent vectors
+// is a dirty fix for correct memory disallocation: It would be
+// better, if either the fields of heist are never touched
+// (i.e. changed) except in hInit, or, if hInit would return the
+// "security" copy as well. But then, all the relevant data is held in
+// global variables, so we might do that here, as well.
+static scfmon hsecure= NULL;
+
+scfmon hInit(ideal S, ideal Q, int *Nexist, ring tailRing)
+{
+  int  sl, ql, i, k = 0;
+  polyset si, qi, ss;
+  scfmon ex, ek;
+  if (tailRing != currRing)
+    hisModule = id_RankFreeModule(S, currRing, tailRing);
+  else
+    hisModule = id_RankFreeModule(S,currRing);
+  if (hisModule < 0)
+    hisModule = 0;
+  if (S!=NULL)
+  {
+    si = S->m;
+    sl = IDELEMS(S);
+  }
+  else
+  {
+    si = NULL;
+    sl = 0;
+  }
+  if (Q!=NULL)
+  {
+    qi = Q->m;
+    ql = IDELEMS(Q);
+  }
+  else
+  {
+    qi = NULL;
+    ql = 0;
+  }
+  if ((sl + ql) == 0)
+  {
+    *Nexist = 0;
+    return NULL;
+  }
+  ss = si;
+  for (i = sl; i>0; i--)
+  {
+    if (*ss!=0)
+      k++;
+    ss++;
+  }
+  ss = qi;
+  for (i = ql; i>0; i--)
+  {
+    if (*ss!=0)
+      k++;
+    ss++;
+  }
+  *Nexist = k;
+  if (k==0)
+    return NULL;
+  ek = ex = (scfmon)omAlloc0(k * sizeof(scmon));
+  hsecure = (scfmon) omAlloc0(k * sizeof(scmon));
+  for (i = sl; i>0; i--)
+  {
+    if (*si!=NULL)
+    {
+      *ek = (scmon) omAlloc(((currRing->N)+1)*sizeof(int));
+      pGetExpV(*si, *ek);
+      ek++;
+    }
+    si++;
+  }
+  for (i = ql; i>0; i--)
+  {
+    if (*qi!=NULL)
+    {
+      *ek = (scmon) omAlloc(((currRing->N)+1)*sizeof(int));
+      pGetExpV(*qi, *ek);
+      ek++;
+    }
+    qi++;
+  }
+  memcpy(hsecure, ex, k * sizeof(scmon));
+  return ex;
+}
+
+#if 0
+void hWeight()
+{
+  int i, k;
+  int x;
+
+  i = (currRing->N);
+  loop
+  {
+    if (pWeight(i) != 1) break;
+    i--;
+    if (i == 0) return;
+  }
+  for (i=(currRing->N); i>0; i--)
+  {
+    x = pWeight(i);
+    if (x != 1)
+    {
+      for (k=hNexist-1; k>=0; k--)
+      {
+        hexist[k][i] *= x;
+      }
+    }
+  }
+}
+#endif
+
+void hDelete(scfmon ev, int ev_length)
+{
+  int i;
+
+  if (ev_length>0)
+  {
+    for (i=ev_length-1;i>=0;i--)
+      omFreeSize(hsecure[i],((currRing->N)+1)*sizeof(int));
+    omFreeSize(hsecure, ev_length*sizeof(scmon));
+    omFreeSize(ev,  ev_length*sizeof(scmon));
+  }
+}
+
+
+void hComp(scfmon exist, int Nexist, int ak, scfmon stc, int *Nstc)
+{
+  int k = 0;
+  scfmon ex = exist, co = stc;
+  int i;
+
+  for (i = Nexist; i>0; i--)
+  {
+    if (((**ex) == 0) || ((**ex) == ak))
+    {
+      *co = *ex;
+      co++;
+      k++;
+    }
+    ex++;
+  }
+  *Nstc = k;
+}
+
+
+void hSupp(scfmon stc, int Nstc, varset var, int *Nvar)
+{
+  int  nv, i0, i1, i, j;
+  nv = i0 = *Nvar;
+  i1 = 0;
+  for (i = 1; i <= nv; i++)
+  {
+    j = 0;
+    loop
+    {
+      if (stc[j][i]>0)
+      {
+        i1++;
+        var[i1] = i;
+        break;
+      }
+      j++;
+      if (j == Nstc)
+      {
+        var[i0] = i;
+        i0--;
+        break;
+      }
+    }
+  }
+  *Nvar = i1;
+}
+
+void hOrdSupp(scfmon stc, int Nstc, varset var, int Nvar)
+{
+  int  i, i1, j, jj, k, l;
+  int  x;
+  scmon temp, count;
+  float o, h, g, *v1;
+
+  v1 = (float *)omAlloc(Nvar * sizeof(float));
+  temp = (int *)omAlloc(Nstc * sizeof(int));
+  count = (int *)omAlloc(Nstc * sizeof(int));
+  for (i = 1; i <= Nvar; i++)
+  {
+    i1 = var[i];
+    *temp = stc[0][i1];
+    *count = 1;
+    jj = 1;
+    for (j = 1; j < Nstc; j++)
+    {
+      x = stc[j][i1];
+      k = 0;
+      loop
+      {
+        if (x > temp[k])
+        {
+          k++;
+          if (k == jj)
+          {
+            temp[k] = x;
+            count[k] = 1;
+            jj++;
+            break;
+          }
+        }
+        else if (x < temp[k])
+        {
+          for (l = jj; l > k; l--)
+          {
+            temp[l] = temp[l-1];
+            count[l] = count[l-1];
+          }
+          temp[k] = x;
+          count[k] = 1;
+          jj++;
+          break;
+        }
+        else
+        {
+          count[k]++;
+          break;
+        }
+      }
+    }
+    h = 0.0;
+    o = (float)Nstc/(float)jj;
+    for(j = 0; j < jj; j++)
+    {
+       g = (float)count[j];
+       if (g > o)
+         g -= o;
+       else
+         g = o - g;
+       if (g > h)
+         h = g;
+    }
+    v1[i-1] = h * (float)jj;
+  }
+  omFreeSize((ADDRESS)count, Nstc * sizeof(int));
+  omFreeSize((ADDRESS)temp, Nstc * sizeof(int));
+  for (i = 1; i < Nvar; i++)
+  {
+    i1 = var[i+1];
+    h = v1[i];
+    j = 0;
+    loop
+    {
+      if (h > v1[j])
+      {
+        for (l = i; l > j; l--)
+        {
+          v1[l] = v1[l-1];
+          var[l+1] = var[l];
+        }
+        v1[j] = h;
+        var[j+1] = i1;
+        break;
+      }
+      j++;
+      if (j == i)
+        break;
+    }
+  }
+  omFreeSize((ADDRESS)v1, Nvar * sizeof(float));
+}
+
+
+static void hShrink(scfmon co, int a, int Nco)
+{
+  while ((co[a]!=NULL) && (a<Nco)) a++;
+  int i = a;
+  int j;
+  for (j = a; j < Nco; j++)
+  {
+    if (co[j]!=NULL)
+    {
+      co[i] = co[j];
+      i++;
+    }
+  }
+}
+
+
+void hStaircase(scfmon stc, int *Nstc, varset var, int Nvar)
+{
+  int  nc = *Nstc;
+  if (nc < 2)
+    return;
+  int z = 0;
+  int i = 0;
+  int j = 1;
+  scmon n = stc[1 /*j*/];
+  scmon o = stc[0];
+  int k = Nvar;
+  loop
+  {
+    int k1 = var[k];
+    if (o[k1] > n[k1])
+    {
+      loop
+      {
+        k--;
+        if (k==0)
+        {
+          stc[i] = NULL;
+          z++;
+          break;
+        }
+        else
+        {
+          k1 = var[k];
+          if (o[k1] < n[k1])
+            break;
+        }
+      }
+      k = Nvar;
+    }
+    else if (o[k1] < n[k1])
+    {
+      loop
+      {
+        k--;
+        if (k==0)
+        {
+          stc[j] = NULL;
+          z++;
+          break;
+        }
+        else
+        {
+          k1 = var[k];
+          if (o[k1] > n[k1])
+            break;
+        }
+      }
+      k = Nvar;
+    }
+    else
+    {
+      k--;
+      if (k==0)
+      {
+        stc[j] = NULL;
+        z++;
+        k = Nvar;
+      }
+    }
+    if (k == Nvar)
+    {
+      if (stc[j]==NULL)
+        i = j - 1;
+      loop
+      {
+        i++;
+        if (i == j)
+        {
+          i = -1;
+          j++;
+          if (j < nc)
+            n = stc[j];
+          else
+          {
+            if (z!=0)
+            {
+              *Nstc -= z;
+              hShrink(stc, 0, nc);
+            }
+            return;
+          }
+        }
+        else if (stc[i]!=NULL)
+        {
+          o = stc[i];
+          break;
+        }
+      }
+    }
+  }
+}
+
+
+void hRadical(scfmon rad, int *Nrad, int Nvar)
+{
+  int  nc = *Nrad, z = 0, i, j, k;
+  scmon n, o;
+  if (nc < 2)
+    return;
+  i = 0;
+  j = 1;
+  n = rad[j];
+  o = rad[0];
+  k = Nvar;
+  loop
+  {
+    if ((o[k]!=0) && (n[k]==0))
+    {
+      loop
+      {
+        k--;
+        if (k==0)
+        {
+          rad[i] = NULL;
+          z++;
+          break;
+        }
+        else
+        {
+          if ((o[k]==0) && (n[k]!=0))
+            break;
+        }
+      }
+      k = Nvar;
+    }
+    else if (!o[k] && n[k])
+    {
+      loop
+      {
+        k--;
+        if (!k)
+        {
+          rad[j] = NULL;
+          z++;
+          break;
+        }
+        else
+        {
+          if (o[k] && !n[k])
+            break;
+        }
+      }
+      k = Nvar;
+    }
+    else
+    {
+      k--;
+      if (!k)
+      {
+        rad[j] = NULL;
+        z++;
+        k = Nvar;
+      }
+    }
+    if (k == Nvar)
+    {
+      if (!rad[j])
+        i = j - 1;
+      loop
+      {
+        i++;
+        if (i == j)
+        {
+          i = -1;
+          j++;
+          if (j < nc)
+            n = rad[j];
+          else
+          {
+            if (z)
+            {
+              *Nrad -= z;
+              hShrink(rad, 0, nc);
+            }
+            return;
+          }
+        }
+        else if (rad[i])
+        {
+          o = rad[i];
+          break;
+        }
+      }
+    }
+  }
+}
+
+
+void hLexS(scfmon stc, int Nstc, varset var, int Nvar)
+{
+  if (Nstc < 2)
+    return;
+  int  j = 1, i = 0;
+  scmon n = stc[j];
+  scmon o = stc[0];
+  int k = Nvar;
+  loop
+  {
+    int k1 = var[k];
+    if (o[k1] < n[k1])
+    {
+      i++;
+      if (i < j)
+      {
+        o = stc[i];
+        k = Nvar;
+      }
+      else
+      {
+        j++;
+        if (j < Nstc)
+        {
+          i = 0;
+          o = stc[0];
+          n = stc[j];
+          k = Nvar;
+        }
+        else
+          return;
+      }
+    }
+    else if (o[k1] > n[k1])
+    {
+      int tmp_k;
+      for (tmp_k = j; tmp_k > i; tmp_k--)
+        stc[tmp_k] = stc[tmp_k - 1];
+      stc[i] = n;
+      j++;
+      if (j < Nstc)
+      {
+        i = 0;
+        o = stc[0];
+        n = stc[j];
+        k = Nvar;
+      }
+      else
+        return;
+    }
+    else
+    {
+      k--;
+      if (k<=0) return;
+    }
+  }
+}
+
+
+void hLexR(scfmon rad, int Nrad, varset var, int Nvar)
+{
+  int  j = 1, i = 0, k, k1;
+  scmon n, o;
+  if (Nrad < 2)
+    return;
+  n = rad[j];
+  o = rad[0];
+  k = Nvar;
+  loop
+  {
+    k1 = var[k];
+    if (!o[k1] && n[k1])
+    {
+      i++;
+      if (i < j)
+      {
+        o = rad[i];
+        k = Nvar;
+      }
+      else
+      {
+        j++;
+        if (j < Nrad)
+        {
+          i = 0;
+          o = rad[0];
+          n = rad[j];
+          k = Nvar;
+        }
+        else
+          return;
+      }
+    }
+    else if (o[k1] && !n[k1])
+    {
+      for (k = j; k > i; k--)
+        rad[k] = rad[k - 1];
+      rad[i] = n;
+      j++;
+      if (j < Nrad)
+      {
+        i = 0;
+        o = rad[0];
+        n = rad[j];
+        k = Nvar;
+      }
+      else
+        return;
+    }
+    else
+      k--;
+  }
+}
+
+
+void hPure(scfmon stc, int a, int *Nstc, varset var, int Nvar,
+ scmon pure, int *Npure)
+{
+  int  nc = *Nstc, np = 0, nq = 0, j, i, i1, c, l;
+  scmon x;
+  for (j = a; j < nc; j++)
+  {
+    x = stc[j];
+    i = Nvar;
+    c = 2;
+    l = 0;
+    loop
+    {
+      i1 = var[i];
+      if (x[i1])
+      {
+        c--;
+        if (!c)
+        {
+          l = 0;
+          break;
+        }
+        else if (c == 1)
+          l = i1;
+      }
+      i--;
+      if (!i)
+        break;
+    }
+    if (l)
+    {
+      if (!pure[l])
+      {
+        np++;
+        pure[l] = x[l];
+      }
+      else if (x[l] < pure[l])
+        pure[l] = x[l];
+      stc[j] = NULL;
+      nq++;
+    }
+  }
+  *Npure = np;
+  if (nq!=0)
+  {
+    *Nstc -= nq;
+    hShrink(stc, a, nc);
+  }
+}
+
+
+void hElimS(scfmon stc, int *e1, int a2, int e2, varset var, int Nvar)
+{
+  int  nc = *e1, z = 0, i, j, k, k1;
+  scmon n, o;
+  if (!nc || (a2 == e2))
+    return;
+  j = 0;
+  i = a2;
+  o = stc[i];
+  n = stc[0];
+  k = Nvar;
+  loop
+  {
+    k1 = var[k];
+    if (o[k1] > n[k1])
+    {
+      k = Nvar;
+      i++;
+      if (i < e2)
+        o = stc[i];
+      else
+      {
+        j++;
+        if (j < nc)
+        {
+          i = a2;
+          o = stc[i];
+          n = stc[j];
+        }
+        else
+        {
+          if (z!=0)
+          {
+            *e1 -= z;
+            hShrink(stc, 0, nc);
+          }
+          return;
+        }
+      }
+    }
+    else
+    {
+      k--;
+      if (k==0)
+      {
+        stc[j] = NULL;
+        z++;
+        j++;
+        if (j < nc)
+        {
+          i = a2;
+          o = stc[i];
+          n = stc[j];
+          k = Nvar;
+        }
+        else
+        {
+          if (z!=0)
+          {
+            *e1 -= z;
+            hShrink(stc, 0, nc);
+          }
+          return;
+        }
+      }
+    }
+  }
+}
+
+
+void hElimR(scfmon rad, int *e1, int a2, int e2, varset var, int Nvar)
+{
+  int  nc = *e1, z = 0, i, j, k, k1;
+  scmon n, o;
+  if (!nc || (a2 == e2))
+    return;
+  j = 0;
+  i = a2;
+  o = rad[i];
+  n = rad[0];
+  k = Nvar;
+  loop
+  {
+    k1 = var[k];
+    if (o[k1] && !n[k1])
+    {
+      k = Nvar;
+      i++;
+      if (i < e2)
+        o = rad[i];
+      else
+      {
+        j++;
+        if (j < nc)
+        {
+          i = a2;
+          o = rad[i];
+          n = rad[j];
+        }
+        else
+        {
+          if (z!=0)
+          {
+            *e1 -= z;
+            hShrink(rad, 0, nc);
+          }
+          return;
+        }
+      }
+    }
+    else
+    {
+      k--;
+      if (!k)
+      {
+        rad[j] = NULL;
+        z++;
+        j++;
+        if (j < nc)
+        {
+          i = a2;
+          o = rad[i];
+          n = rad[j];
+          k = Nvar;
+        }
+        else
+        {
+          if (z!=0)
+          {
+            *e1 -= z;
+            hShrink(rad, 0, nc);
+          }
+          return;
+        }
+      }
+    }
+  }
+}
+
+
+void hLex2S(scfmon rad, int e1, int a2, int e2, varset var,
+ int Nvar, scfmon w)
+{
+  int  j0 = 0, j = 0, i = a2, k, k1;
+  scmon n, o;
+  if (!e1)
+  {
+    for (; i < e2; i++)
+      rad[i - a2] = rad[i];
+    return;
+  }  else if (i == e2)
+    return;
+  n = rad[j];
+  o = rad[i];
+  loop
+  {
+    k = Nvar;
+    loop
+    {
+      k1 = var[k];
+      if (o[k1] < n[k1])
+      {
+        w[j0] = o;
+        j0++;
+        i++;
+        if (i < e2)
+        {
+          o = rad[i];
+          break;
+        }
+        else
+        {
+          for (; j < e1; j++)
+          {
+            w[j0] = rad[j];
+            j0++;
+          }
+          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
+          return;
+        }
+      }
+      else if (o[k1] > n[k1])
+      {
+        w[j0] = n;
+        j0++;
+        j++;
+        if (j < e1)
+        {
+          n = rad[j];
+          break;
+        }
+        else
+        {
+          for (; i < e2; i++)
+          {
+            w[j0] = rad[i];
+            j0++;
+          }
+          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
+          return;
+        }
+      }
+      k--;
+    }
+  }
+}
+
+
+void hLex2R(scfmon rad, int e1, int a2, int e2, varset var,
+ int Nvar, scfmon w)
+{
+  int  j0 = 0, j = 0, i = a2, k, k1;
+  scmon n, o;
+  if (!e1)
+  {
+    for (; i < e2; i++)
+      rad[i - a2] = rad[i];
+    return;
+  }
+  else if (i == e2)
+    return;
+  n = rad[j];
+  o = rad[i];
+  loop
+  {
+    k = Nvar;
+    loop
+    {
+      k1 = var[k];
+      if (!o[k1] && n[k1])
+      {
+        w[j0] = o;
+        j0++;
+        i++;
+        if (i < e2)
+        {
+          o = rad[i];
+          break;
+        }
+        else
+        {
+          for (; j < e1; j++)
+          {
+            w[j0] = rad[j];
+            j0++;
+          }
+          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
+          return;
+        }
+      }
+      else if (o[k1] && !n[k1])
+      {
+        w[j0] = n;
+        j0++;
+        j++;
+        if (j < e1)
+        {
+          n = rad[j];
+          break;
+        }
+        else
+        {
+          for (; i < e2; i++)
+          {
+            w[j0] = rad[i];
+            j0++;
+          }
+          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
+          return;
+        }
+      }
+      k--;
+    }
+  }
+}
+
+
+void hStepS(scfmon stc, int Nstc, varset var, int Nvar, int *a, int *x)
+{
+  int  k1, i;
+  int  y;
+  k1 = var[Nvar];
+  y = *x;
+  i = *a;
+  loop
+  {
+    if (y < stc[i][k1])
+    {
+      *a = i;
+      *x = stc[i][k1];
+      return;
+    }
+    i++;
+    if (i == Nstc)
+    {
+      *a = i;
+      return;
+    }
+  }
+}
+
+
+void hStepR(scfmon rad, int Nrad, varset var, int Nvar, int *a)
+{
+  int  k1, i;
+  k1 = var[Nvar];
+  i = 0;
+  loop
+  {
+    if (rad[i][k1])
+    {
+      *a = i;
+      return;
+    }
+    i++;
+    if (i == Nrad)
+    {
+      *a = i;
+      return;
+    }
+  }
+}
+
+
+monf hCreate(int Nvar)
+{
+  monf xmem;
+  int  i;
+  xmem = (monf)omAlloc((Nvar + 1) * sizeof(monp));
+  for (i = Nvar; i>0; i--)
+  {
+    xmem[i] = (monp)omAlloc(LEN_MON);
+    xmem[i]->mo = NULL;
+  }
+  return xmem;
+}
+
+
+void hKill(monf xmem, int Nvar)
+{
+  int  i;
+  for (i = Nvar; i!=0; i--)
+  {
+    if (xmem[i]->mo!=NULL)
+      omFreeSize((ADDRESS)xmem[i]->mo, xmem[i]->a * sizeof(scmon));
+    omFreeSize((ADDRESS)xmem[i], LEN_MON);
+  }
+  omFreeSize((ADDRESS)xmem, (Nvar + 1) * sizeof(monp));
+}
+
+
+scfmon hGetmem(int lm, scfmon old, monp monmem)
+{
+  scfmon x = monmem->mo;
+  int  lx = monmem->a;
+  if ((x==NULL) || (lm > lx))
+  {
+    /* according to http://www.singular.uni-kl.de:8002/trac/ticket/463#comment:4
+     * we need to work around a compiler bug:
+    * if ((x!=NULL)&&(lx>0)) omFreeSize((ADDRESS)x, lx * sizeof(scmon));
+    */
+    if (x!=NULL) if (lx>0) omFreeSize((ADDRESS)x, lx * sizeof(scmon));
+    monmem->mo = x = (scfmon)omAlloc(lm * sizeof(scmon));
+    monmem->a = lm;
+  }
+  memcpy(x, old, lm * sizeof(scmon));
+  return x;
+}
+
+/*
+* a bug in Metrowerks with "lifetime analysis"
+*scmon hGetpure(scmon p)
+*{
+*  scmon p1, pn;
+*  p1 = p + 1;
+*  pn = p1 + (currRing->N);
+*  memcpy(pn, p1, (currRing->N) * sizeof(int));
+*  return pn - 1;
+*}
+*/
+scmon hGetpure(scmon p)
+{
+  scmon p1 = p;
+  scmon pn;
+  p1++;
+  pn = p1;
+  pn += (currRing->N);
+  memcpy(pn, p1, (currRing->N) * sizeof(int));
+  return pn - 1;
+}
+
Index: rnel/hdegree.cc
===================================================================
--- kernel/hdegree.cc	(revision 66b1de771706655cf8a9c3924d1dac08150c8814)
+++ 	(revision )
@@ -1,1464 +1,0 @@
-/****************************************
-*  Computer Algebra System SINGULAR     *
-****************************************/
-/*
-*  ABSTRACT -  dimension, multiplicity, HC, kbase
-*/
-
-#ifdef HAVE_CONFIG_H
-#include "singularconfig.h"
-#endif /* HAVE_CONFIG_H */
-#include <kernel/mod2.h>
-
-
-#include <kernel/structs.h>
-#include <kernel/febase.h>
-#include <omalloc/omalloc.h>
-#include <kernel/ideals.h>
-#include <kernel/polys.h>
-#include <misc/intvec.h>
-#include <coeffs/numbers.h>
-#include <kernel/hutil.h>
-#include <kernel/stairc.h>
-
-int  hCo, hMu, hMu2;
-omBin indlist_bin = omGetSpecBin(sizeof(indlist));
-
-/*0 implementation*/
-
-// dimension
-
-void hDimSolve(scmon pure, int Npure, scfmon rad, int Nrad,
- varset var, int Nvar)
-{
-  int  dn, iv, rad0, b, c, x;
-  scmon pn;
-  scfmon rn;
-  if (Nrad < 2)
-  {
-    dn = Npure + Nrad;
-    if (dn < hCo)
-      hCo = dn;
-    return;
-  }
-  if (Npure+1 >= hCo)
-    return;
-  iv = Nvar;
-  while(pure[var[iv]]) iv--;
-  hStepR(rad, Nrad, var, iv, &rad0);
-  if (rad0!=0)
-  {
-    iv--;
-    if (rad0 < Nrad)
-    {
-      pn = hGetpure(pure);
-      rn = hGetmem(Nrad, rad, radmem[iv]);
-      hDimSolve(pn, Npure + 1, rn, rad0, var, iv);
-      b = rad0;
-      c = Nrad;
-      hElimR(rn, &rad0, b, c, var, iv);
-      hPure(rn, b, &c, var, iv, pn, &x);
-      hLex2R(rn, rad0, b, c, var, iv, hwork);
-      rad0 += (c - b);
-      hDimSolve(pn, Npure + x, rn, rad0, var, iv);
-    }
-    else
-    {
-      hDimSolve(pure, Npure, rad, Nrad, var, iv);
-    }
-  }
-  else
-    hCo = Npure + 1;
-}
-
-int  scDimInt(ideal S, ideal Q)
-{
-  int  mc;
-  hexist = hInit(S, Q, &hNexist);
-  if (!hNexist)
-    return (currRing->N);
-  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
-  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  mc = hisModule;
-  if (!mc)
-  {
-    hrad = hexist;
-    hNrad = hNexist;
-  }
-  else
-    hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  radmem = hCreate((currRing->N) - 1);
-  hCo = (currRing->N) + 1;
-  loop
-  {
-    if (mc)
-      hComp(hexist, hNexist, mc, hrad, &hNrad);
-    if (hNrad)
-    {
-      hNvar = (currRing->N);
-      hRadical(hrad, &hNrad, hNvar);
-      hSupp(hrad, hNrad, hvar, &hNvar);
-      if (hNvar)
-      {
-        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
-        hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
-        hLexR(hrad, hNrad, hvar, hNvar);
-        hDimSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
-      }
-    }
-    else
-    {
-      hCo = 0;
-      break;
-    }
-    mc--;
-    if (mc <= 0)
-      break;
-  }
-  hKill(radmem, (currRing->N) - 1);
-  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
-  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
-  hDelete(hexist, hNexist);
-  if (hisModule)
-    omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
-  return (currRing->N) - hCo;
-}
-
-// independent set
-static scmon hInd;
-
-static void hIndSolve(scmon pure, int Npure, scfmon rad, int Nrad,
- varset var, int Nvar)
-{
-  int  dn, iv, rad0, b, c, x;
-  scmon pn;
-  scfmon rn;
-  if (Nrad < 2)
-  {
-    dn = Npure + Nrad;
-    if (dn < hCo)
-    {
-      hCo = dn;
-      for (iv=(currRing->N); iv; iv--)
-      {
-        if (pure[iv])
-          hInd[iv] = 0;
-        else
-          hInd[iv] = 1;
-      }
-      if (Nrad)
-      {
-        pn = *rad;
-        iv = Nvar;
-        loop
-        {
-          x = var[iv];
-          if (pn[x])
-          {
-            hInd[x] = 0;
-            break;
-          }
-          iv--;
-        }
-      }
-    }
-    return;
-  }
-  if (Npure+1 >= hCo)
-    return;
-  iv = Nvar;
-  while(pure[var[iv]]) iv--;
-  hStepR(rad, Nrad, var, iv, &rad0);
-  if (rad0)
-  {
-    iv--;
-    if (rad0 < Nrad)
-    {
-      pn = hGetpure(pure);
-      rn = hGetmem(Nrad, rad, radmem[iv]);
-      pn[var[iv + 1]] = 1;
-      hIndSolve(pn, Npure + 1, rn, rad0, var, iv);
-      pn[var[iv + 1]] = 0;
-      b = rad0;
-      c = Nrad;
-      hElimR(rn, &rad0, b, c, var, iv);
-      hPure(rn, b, &c, var, iv, pn, &x);
-      hLex2R(rn, rad0, b, c, var, iv, hwork);
-      rad0 += (c - b);
-      hIndSolve(pn, Npure + x, rn, rad0, var, iv);
-    }
-    else
-    {
-      hIndSolve(pure, Npure, rad, Nrad, var, iv);
-    }
-  }
-  else
-  {
-    hCo = Npure + 1;
-    for (x=(currRing->N); x; x--)
-    {
-      if (pure[x])
-        hInd[x] = 0;
-      else
-        hInd[x] = 1;
-    }
-    hInd[var[iv]] = 0;
-  }
-}
-
-intvec * scIndIntvec(ideal S, ideal Q)
-{
-  intvec *Set=new intvec((currRing->N));
-  int  mc,i;
-  hexist = hInit(S, Q, &hNexist);
-  if (hNexist==0)
-  {
-    for(i=0; i<(currRing->N); i++)
-      (*Set)[i]=1;
-    return Set;
-  }
-  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
-  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  hInd = (scmon)omAlloc0((1 + (currRing->N)) * sizeof(int));
-  mc = hisModule;
-  if (mc==0)
-  {
-    hrad = hexist;
-    hNrad = hNexist;
-  }
-  else
-    hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  radmem = hCreate((currRing->N) - 1);
-  hCo = (currRing->N) + 1;
-  loop
-  {
-    if (mc!=0)
-      hComp(hexist, hNexist, mc, hrad, &hNrad);
-    if (hNrad!=0)
-    {
-      hNvar = (currRing->N);
-      hRadical(hrad, &hNrad, hNvar);
-      hSupp(hrad, hNrad, hvar, &hNvar);
-      if (hNvar!=0)
-      {
-        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
-        hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
-        hLexR(hrad, hNrad, hvar, hNvar);
-        hIndSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
-      }
-    }
-    else
-    {
-      hCo = 0;
-      break;
-    }
-    mc--;
-    if (mc <= 0)
-      break;
-  }
-  for(i=0; i<(currRing->N); i++)
-    (*Set)[i] = hInd[i+1];
-  hKill(radmem, (currRing->N) - 1);
-  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  omFreeSize((ADDRESS)hInd, (1 + (currRing->N)) * sizeof(int));
-  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
-  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
-  hDelete(hexist, hNexist);
-  if (hisModule)
-    omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
-  return Set;
-}
-
-indset ISet, JSet;
-
-static BOOLEAN hNotZero(scfmon rad, int Nrad, varset var, int Nvar)
-{
-  int  k1, i;
-  k1 = var[Nvar];
-  i = 0;
-  loop
-  {
-    if (rad[i][k1]==0)
-      return FALSE;
-    i++;
-    if (i == Nrad)
-      return TRUE;
-  }
-}
-
-static void hIndep(scmon pure)
-{
-  int iv;
-  intvec *Set;
-
-  Set = ISet->set = new intvec((currRing->N));
-  for (iv=(currRing->N); iv!=0 ; iv--)
-  {
-    if (pure[iv])
-      (*Set)[iv-1] = 0;
-    else
-      (*Set)[iv-1] = 1;
-  }
-  ISet = ISet->nx = (indset)omAlloc0Bin(indlist_bin);
-  hMu++;
-}
-
-void hIndMult(scmon pure, int Npure, scfmon rad, int Nrad,
- varset var, int Nvar)
-{
-  int  dn, iv, rad0, b, c, x;
-  scmon pn;
-  scfmon rn;
-  if (Nrad < 2)
-  {
-    dn = Npure + Nrad;
-    if (dn == hCo)
-    {
-      if (Nrad==0)
-        hIndep(pure);
-      else
-      {
-        pn = *rad;
-        for (iv = Nvar; iv!=0; iv--)
-        {
-          x = var[iv];
-          if (pn[x])
-          {
-            pure[x] = 1;
-            hIndep(pure);
-            pure[x] = 0;
-          }
-        }
-      }
-    }
-    return;
-  }
-  iv = Nvar;
-  dn = Npure+1;
-  if (dn >= hCo)
-  {
-    if (dn > hCo)
-      return;
-    loop
-    {
-      if(!pure[var[iv]])
-      {
-        if(hNotZero(rad, Nrad, var, iv))
-        {
-          pure[var[iv]] = 1;
-          hIndep(pure);
-          pure[var[iv]] = 0;
-        }
-      }
-      iv--;
-      if (!iv)
-        return;
-    }
-  }
-  while(pure[var[iv]]) iv--;
-  hStepR(rad, Nrad, var, iv, &rad0);
-  iv--;
-  if (rad0 < Nrad)
-  {
-    pn = hGetpure(pure);
-    rn = hGetmem(Nrad, rad, radmem[iv]);
-    pn[var[iv + 1]] = 1;
-    hIndMult(pn, Npure + 1, rn, rad0, var, iv);
-    pn[var[iv + 1]] = 0;
-    b = rad0;
-    c = Nrad;
-    hElimR(rn, &rad0, b, c, var, iv);
-    hPure(rn, b, &c, var, iv, pn, &x);
-    hLex2R(rn, rad0, b, c, var, iv, hwork);
-    rad0 += (c - b);
-    hIndMult(pn, Npure + x, rn, rad0, var, iv);
-  }
-  else
-  {
-    hIndMult(pure, Npure, rad, Nrad, var, iv);
-  }
-}
-
-/*3
-* consider indset x := !pure
-* (for all i) (if(sm(i) > x) return FALSE)
-* else return TRUE
-*/
-static BOOLEAN hCheck1(indset sm, scmon pure)
-{
-  int iv;
-  intvec *Set;
-  while (sm->nx != NULL)
-  {
-    Set = sm->set;
-    iv=(currRing->N);
-    loop
-    {
-      if (((*Set)[iv-1] == 0) && (pure[iv] == 0))
-        break;
-      iv--;
-      if (iv == 0)
-        return FALSE;
-    }
-    sm = sm->nx;
-  }
-  return TRUE;
-}
-
-/*3
-* consider indset x := !pure
-* (for all i) if(x > sm(i)) delete sm(i))
-* return (place for x)
-*/
-static indset hCheck2(indset sm, scmon pure)
-{
-  int iv;
-  intvec *Set;
-  indset be, a1 = NULL;
-  while (sm->nx != NULL)
-  {
-    Set = sm->set;
-    iv=(currRing->N);
-    loop
-    {
-      if ((pure[iv] == 1) && ((*Set)[iv-1] == 1))
-        break;
-      iv--;
-      if (iv == 0)
-      {
-        if (a1 == NULL)
-        {
-          a1 = sm;
-        }
-        else
-        {
-          hMu2--;
-          be->nx = sm->nx;
-          delete Set;
-          omFreeBin((ADDRESS)sm, indlist_bin);
-          sm = be;
-        }
-        break;
-      }
-    }
-    be = sm;
-    sm = sm->nx;
-  }
-  if (a1 != NULL)
-  {
-    return a1;
-  }
-  else
-  {
-    hMu2++;
-    sm->set = new intvec((currRing->N));
-    sm->nx = (indset)omAlloc0Bin(indlist_bin);
-    return sm;
-  }
-}
-
-/*2
-*  definition x >= y
-*      x(i) == 0 => y(i) == 0
-*      > ex. j with x(j) == 1 and y(j) == 0
-*/
-static void hCheckIndep(scmon pure)
-{
-  intvec *Set;
-  indset res;
-  int iv;
-  if (hCheck1(ISet, pure))
-  {
-    if (hCheck1(JSet, pure))
-    {
-      res = hCheck2(JSet,pure);
-      if (res == NULL)
-        return;
-      Set = res->set;
-      for (iv=(currRing->N); iv; iv--)
-      {
-        if (pure[iv])
-          (*Set)[iv-1] = 0;
-        else
-          (*Set)[iv-1] = 1;
-      }
-    }
-  }
-}
-
-void hIndAllMult(scmon pure, int Npure, scfmon rad, int Nrad,
- varset var, int Nvar)
-{
-  int  dn, iv, rad0, b, c, x;
-  scmon pn;
-  scfmon rn;
-  if (Nrad < 2)
-  {
-    dn = Npure + Nrad;
-    if (dn > hCo)
-    {
-      if (!Nrad)
-        hCheckIndep(pure);
-      else
-      {
-        pn = *rad;
-        for (iv = Nvar; iv; iv--)
-        {
-          x = var[iv];
-          if (pn[x])
-          {
-            pure[x] = 1;
-            hCheckIndep(pure);
-            pure[x] = 0;
-          }
-        }
-      }
-    }
-    return;
-  }
-  iv = Nvar;
-  while(pure[var[iv]]) iv--;
-  hStepR(rad, Nrad, var, iv, &rad0);
-  iv--;
-  if (rad0 < Nrad)
-  {
-    pn = hGetpure(pure);
-    rn = hGetmem(Nrad, rad, radmem[iv]);
-    pn[var[iv + 1]] = 1;
-    hIndAllMult(pn, Npure + 1, rn, rad0, var, iv);
-    pn[var[iv + 1]] = 0;
-    b = rad0;
-    c = Nrad;
-    hElimR(rn, &rad0, b, c, var, iv);
-    hPure(rn, b, &c, var, iv, pn, &x);
-    hLex2R(rn, rad0, b, c, var, iv, hwork);
-    rad0 += (c - b);
-    hIndAllMult(pn, Npure + x, rn, rad0, var, iv);
-  }
-  else
-  {
-    hIndAllMult(pure, Npure, rad, Nrad, var, iv);
-  }
-}
-
-// multiplicity
-
-static int hZeroMult(scmon pure, scfmon stc, int Nstc, varset var, int Nvar)
-{
-  int  iv = Nvar -1, sum, a, a0, a1, b, i;
-  int  x, x0;
-  scmon pn;
-  scfmon sn;
-  if (!iv)
-    return pure[var[1]];
-  else if (!Nstc)
-  {
-    sum = 1;
-    for (i = Nvar; i; i--)
-      sum *= pure[var[i]];
-    return sum;
-  }
-  x = a = 0;
-  pn = hGetpure(pure);
-  sn = hGetmem(Nstc, stc, stcmem[iv]);
-  hStepS(sn, Nstc, var, Nvar, &a, &x);
-  if (a == Nstc)
-    return pure[var[Nvar]] * hZeroMult(pn, sn, a, var, iv);
-  else
-    sum = x * hZeroMult(pn, sn, a, var, iv);
-  b = a;
-  loop
-  {
-    a0 = a;
-    x0 = x;
-    hStepS(sn, Nstc, var, Nvar, &a, &x);
-    hElimS(sn, &b, a0, a, var, iv);
-    a1 = a;
-    hPure(sn, a0, &a1, var, iv, pn, &i);
-    hLex2S(sn, b, a0, a1, var, iv, hwork);
-    b += (a1 - a0);
-    if (a < Nstc)
-    {
-      sum += (x - x0) * hZeroMult(pn, sn, b, var, iv);
-    }
-    else
-    {
-      sum += (pure[var[Nvar]] - x0) * hZeroMult(pn, sn, b, var, iv);
-      return sum;
-    }
-  }
-}
-
-static void hProject(scmon pure, varset sel)
-{
-  int  i, i0, k;
-  i0 = 0;
-  for (i = 1; i <= (currRing->N); i++)
-  {
-    if (pure[i])
-    {
-      i0++;
-      sel[i0] = i;
-    }
-  }
-  i = hNstc;
-  memcpy(hwork, hstc, i * sizeof(scmon));
-  hStaircase(hwork, &i, sel, i0);
-  if ((i0 > 2) && (i > 10))
-    hOrdSupp(hwork, i, sel, i0);
-  memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
-  hPure(hwork, 0, &i, sel, i0, hpur0, &k);
-  hLexS(hwork, i, sel, i0);
-  hMu += hZeroMult(hpur0, hwork, i, sel, i0);
-}
-
-static void hDimMult(scmon pure, int Npure, scfmon rad, int Nrad,
- varset var, int Nvar)
-{
-  int  dn, iv, rad0, b, c, x;
-  scmon pn;
-  scfmon rn;
-  if (Nrad < 2)
-  {
-    dn = Npure + Nrad;
-    if (dn == hCo)
-    {
-      if (!Nrad)
-        hProject(pure, hsel);
-      else
-      {
-        pn = *rad;
-        for (iv = Nvar; iv; iv--)
-        {
-          x = var[iv];
-          if (pn[x])
-          {
-            pure[x] = 1;
-            hProject(pure, hsel);
-            pure[x] = 0;
-          }
-        }
-      }
-    }
-    return;
-  }
-  iv = Nvar;
-  dn = Npure+1;
-  if (dn >= hCo)
-  {
-    if (dn > hCo)
-      return;
-    loop
-    {
-      if(!pure[var[iv]])
-      {
-        if(hNotZero(rad, Nrad, var, iv))
-        {
-          pure[var[iv]] = 1;
-          hProject(pure, hsel);
-          pure[var[iv]] = 0;
-        }
-      }
-      iv--;
-      if (!iv)
-        return;
-    }
-  }
-  while(pure[var[iv]]) iv--;
-  hStepR(rad, Nrad, var, iv, &rad0);
-  iv--;
-  if (rad0 < Nrad)
-  {
-    pn = hGetpure(pure);
-    rn = hGetmem(Nrad, rad, radmem[iv]);
-    pn[var[iv + 1]] = 1;
-    hDimMult(pn, Npure + 1, rn, rad0, var, iv);
-    pn[var[iv + 1]] = 0;
-    b = rad0;
-    c = Nrad;
-    hElimR(rn, &rad0, b, c, var, iv);
-    hPure(rn, b, &c, var, iv, pn, &x);
-    hLex2R(rn, rad0, b, c, var, iv, hwork);
-    rad0 += (c - b);
-    hDimMult(pn, Npure + x, rn, rad0, var, iv);
-  }
-  else
-  {
-    hDimMult(pure, Npure, rad, Nrad, var, iv);
-  }
-}
-
-static void hDegree(ideal S, ideal Q)
-{
-  int  di;
-  int  mc;
-  hexist = hInit(S, Q, &hNexist);
-  if (!hNexist)
-  {
-    hCo = 0;
-    hMu = 1;
-    return;
-  }
-  //hWeight();
-  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
-  hsel = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
-  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  hpur0 = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  mc = hisModule;
-  hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  if (!mc)
-  {
-    memcpy(hrad, hexist, hNexist * sizeof(scmon));
-    hstc = hexist;
-    hNrad = hNstc = hNexist;
-  }
-  else
-    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  radmem = hCreate((currRing->N) - 1);
-  stcmem = hCreate((currRing->N) - 1);
-  hCo = (currRing->N) + 1;
-  di = hCo + 1;
-  loop
-  {
-    if (mc)
-    {
-      hComp(hexist, hNexist, mc, hrad, &hNrad);
-      hNstc = hNrad;
-      memcpy(hstc, hrad, hNrad * sizeof(scmon));
-    }
-    if (hNrad)
-    {
-      hNvar = (currRing->N);
-      hRadical(hrad, &hNrad, hNvar);
-      hSupp(hrad, hNrad, hvar, &hNvar);
-      if (hNvar)
-      {
-        hCo = hNvar;
-        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
-        hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
-        hLexR(hrad, hNrad, hvar, hNvar);
-        hDimSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
-      }
-    }
-    else
-    {
-      hNvar = 1;
-      hCo = 0;
-    }
-    if (hCo < di)
-    {
-      di = hCo;
-      hMu = 0;
-    }
-    if (hNvar && (hCo == di))
-    {
-      if (di && (di < (currRing->N)))
-        hDimMult(hpure, hNpure, hrad, hNrad, hvar, hNvar);
-      else if (!di)
-        hMu++;
-      else
-      {
-        hStaircase(hstc, &hNstc, hvar, hNvar);
-        if ((hNvar > 2) && (hNstc > 10))
-          hOrdSupp(hstc, hNstc, hvar, hNvar);
-        memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
-        hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
-        hLexS(hstc, hNstc, hvar, hNvar);
-        hMu += hZeroMult(hpur0, hstc, hNstc, hvar, hNvar);
-      }
-    }
-    mc--;
-    if (mc <= 0)
-      break;
-  }
-  hCo = di;
-  hKill(stcmem, (currRing->N) - 1);
-  hKill(radmem, (currRing->N) - 1);
-  omFreeSize((ADDRESS)hpur0, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  omFreeSize((ADDRESS)hsel, ((currRing->N) + 1) * sizeof(int));
-  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
-  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
-  omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
-  hDelete(hexist, hNexist);
-  if (hisModule)
-    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
-}
-
-int  scMultInt(ideal S, ideal Q)
-{
-  hDegree(S, Q);
-  return hMu;
-}
-
-void scPrintDegree(int co, int mu)
-{
-  int di = (currRing->N)-co;
-  if (currRing->OrdSgn == 1)
-  {
-    if (di>0)
-      Print("// dimension (proj.)  = %d\n// degree (proj.)   = %d\n", di-1, mu);
-    else
-      Print("// dimension (affine) = 0\n// degree (affine)  = %d\n",       mu);
-  }
-  else
-    Print("// dimension (local)   = %d\n// multiplicity = %d\n", di, mu);
-}
-
-void scDegree(ideal S, intvec *modulweight, ideal Q)
-{
-  int co, mu, l;
-  intvec *hseries2;
-  intvec *hseries1 = hFirstSeries(S, modulweight, Q);
-  l = hseries1->length()-1;
-  if (l > 1)
-    hseries2 = hSecondSeries(hseries1);
-  else
-    hseries2 = hseries1;
-  hDegreeSeries(hseries1, hseries2, &co, &mu);
-  if ((l == 1) &&(mu == 0))
-    scPrintDegree((currRing->N)+1, 0);
-  else
-    scPrintDegree(co, mu);
-  if (l>1)
-    delete hseries1;
-  delete hseries2;
-}
-
-static void hDegree0(ideal S, ideal Q)
-{
-  int  mc;
-  hexist = hInit(S, Q, &hNexist);
-  if (!hNexist)
-  {
-    hMu = -1;
-    return;
-  }
-  else
-    hMu = 0;
-  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
-  hpur0 = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  mc = hisModule;
-  if (!mc)
-  {
-    hstc = hexist;
-    hNstc = hNexist;
-  }
-  else
-    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  stcmem = hCreate((currRing->N) - 1);
-  loop
-  {
-    if (mc)
-    {
-      hComp(hexist, hNexist, mc, hstc, &hNstc);
-      if (!hNstc)
-      {
-        hMu = -1;
-        break;
-      }
-    }
-    hNvar = (currRing->N);
-    for (int i = hNvar; i; i--)
-      hvar[i] = i;
-    hStaircase(hstc, &hNstc, hvar, hNvar);
-    hSupp(hstc, hNstc, hvar, &hNvar);
-    if ((hNvar == (currRing->N)) && (hNstc >= (currRing->N)))
-    {
-      if ((hNvar > 2) && (hNstc > 10))
-        hOrdSupp(hstc, hNstc, hvar, hNvar);
-      memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
-      hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
-      if (hNpure == hNvar)
-      {
-        hLexS(hstc, hNstc, hvar, hNvar);
-        hMu += hZeroMult(hpur0, hstc, hNstc, hvar, hNvar);
-      }
-      else
-        hMu = -1;
-    }
-    else if (hNvar)
-      hMu = -1;
-    mc--;
-    if (mc <= 0 || hMu < 0)
-      break;
-  }
-  hKill(stcmem, (currRing->N) - 1);
-  omFreeSize((ADDRESS)hpur0, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
-  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
-  hDelete(hexist, hNexist);
-  if (hisModule)
-    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
-}
-
-int  scMult0Int(ideal S, ideal Q)
-{
-  hDegree0(S, Q);
-  return hMu;
-}
-
-
-// HC
-
-static poly pWork;
-
-static void hHedge(poly hEdge)
-{
-  pSetm(pWork);
-  if (pLmCmp(pWork, hEdge) == currRing->OrdSgn)
-  {
-    for (int i = hNvar; i>0; i--)
-      pSetExp(hEdge,i, pGetExp(pWork,i));
-    pSetm(hEdge);
-  }
-}
-
-
-static void hHedgeStep(scmon pure, scfmon stc,
-                       int Nstc, varset var, int Nvar,poly hEdge)
-{
-  int  iv = Nvar -1, k = var[Nvar], a, a0, a1, b, i;
-  int  x/*, x0*/;
-  scmon pn;
-  scfmon sn;
-  if (iv==0)
-  {
-    pSetExp(pWork, k, pure[k]);
-    hHedge(hEdge);
-    return;
-  }
-  else if (Nstc==0)
-  {
-    for (i = Nvar; i>0; i--)
-      pSetExp(pWork, var[i], pure[var[i]]);
-    hHedge(hEdge);
-    return;
-  }
-  x = a = 0;
-  pn = hGetpure(pure);
-  sn = hGetmem(Nstc, stc, stcmem[iv]);
-  hStepS(sn, Nstc, var, Nvar, &a, &x);
-  if (a == Nstc)
-  {
-    pSetExp(pWork, k, pure[k]);
-    hHedgeStep(pn, sn, a, var, iv,hEdge);
-    return;
-  }
-  else
-  {
-    pSetExp(pWork, k, x);
-    hHedgeStep(pn, sn, a, var, iv,hEdge);
-  }
-  b = a;
-  loop
-  {
-    a0 = a;
-    // x0 = x;
-    hStepS(sn, Nstc, var, Nvar, &a, &x);
-    hElimS(sn, &b, a0, a, var, iv);
-    a1 = a;
-    hPure(sn, a0, &a1, var, iv, pn, &i);
-    hLex2S(sn, b, a0, a1, var, iv, hwork);
-    b += (a1 - a0);
-    if (a < Nstc)
-    {
-      pSetExp(pWork, k, x);
-      hHedgeStep(pn, sn, b, var, iv,hEdge);
-    }
-    else
-    {
-      pSetExp(pWork, k, pure[k]);
-      hHedgeStep(pn, sn, b, var, iv,hEdge);
-      return;
-    }
-  }
-}
-
-void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge, ring tailRing)
-{
-  int  i;
-  int  k = ak;
-
-  #if HAVE_RINGS
-  if (rField_is_Ring(currRing) && (currRing->OrdSgn == -1))
-  {
-  //consider just monic generators (over rings with zero-divisors)
-  ideal SS=id_Copy(S,tailRing);
-  for(i=0;i<=idElem(SS);i++)
-  	{
-  	if(pIsPurePower(SS->m[i])==0)
-  		p_Delete(&SS->m[i],tailRing);
-  	}
-  	S=id_Copy(SS,tailRing);
-  }
-  #endif
-
-  hNvar = (currRing->N);
-  hexist = hInit(S, Q, &hNexist, tailRing);
-  if (k!=0)
-    hComp(hexist, hNexist, k, hexist, &hNstc);
-  else
-    hNstc = hNexist;
-  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  hvar = (varset)omAlloc((hNvar + 1) * sizeof(int));
-  hpure = (scmon)omAlloc((1 + (hNvar * hNvar)) * sizeof(int));
-  stcmem = hCreate(hNvar - 1);
-  for (i = hNvar; i>0; i--)
-    hvar[i] = i;
-  hStaircase(hexist, &hNstc, hvar, hNvar);
-  if ((hNvar > 2) && (hNstc > 10))
-    hOrdSupp(hexist, hNstc, hvar, hNvar);
-  memset(hpure, 0, (hNvar + 1) * sizeof(int));
-  hPure(hexist, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
-  hLexS(hexist, hNstc, hvar, hNvar);
-  if (hEdge!=NULL)
-    pLmFree(hEdge);
-  hEdge = pInit();
-  pWork = pInit();
-  hHedgeStep(hpure, hexist, hNstc, hvar, hNvar,hEdge);
-  pSetComp(hEdge,ak);
-  hKill(stcmem, hNvar - 1);
-  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
-  omFreeSize((ADDRESS)hvar, (hNvar + 1) * sizeof(int));
-  omFreeSize((ADDRESS)hpure, (1 + (hNvar * hNvar)) * sizeof(int));
-  hDelete(hexist, hNexist);
-  pLmFree(pWork);
-}
-
-
-
-//  kbase
-
-static poly last;
-static scmon act;
-
-static void scElKbase()
-{
-  poly q = pInit();
-  pSetCoeff0(q,nInit(1));
-  pSetExpV(q,act);
-  pNext(q) = NULL;
-  last = pNext(last) = q;
-}
-
-static int scMax( int i, scfmon stc, int Nvar)
-{
-  int x, y=stc[0][Nvar];
-  for (; i;)
-  {
-    i--;
-    x = stc[i][Nvar];
-    if (x > y) y = x;
-  }
-  return y;
-}
-
-static int scMin( int i, scfmon stc, int Nvar)
-{
-  int x, y=stc[0][Nvar];
-  for (; i;)
-  {
-    i--;
-    x = stc[i][Nvar];
-    if (x < y) y = x;
-  }
-  return y;
-}
-
-static int scRestrict( int &Nstc, scfmon stc, int Nvar)
-{
-  int x, y;
-  int i, j, Istc = Nstc;
-
-  y = MAX_INT_VAL;
-  for (i=Nstc-1; i>=0; i--)
-  {
-    j = Nvar-1;
-    loop
-    {
-      if(stc[i][j] != 0) break;
-      j--;
-      if (j == 0)
-      {
-        Istc--;
-        x = stc[i][Nvar];
-        if (x < y) y = x;
-        stc[i] = NULL;
-        break;
-      }
-    }
-  }
-  if (Istc < Nstc)
-  {
-    for (i=Nstc-1; i>=0; i--)
-    {
-      if (stc[i] && (stc[i][Nvar] >= y))
-      {
-        Istc--;
-        stc[i] = NULL;
-      }
-    }
-    j = 0;
-    while (stc[j]) j++;
-    i = j+1;
-    for(; i<Nstc; i++)
-    {
-      if (stc[i])
-      {
-        stc[j] = stc[i];
-        j++;
-      }
-    }
-    Nstc = Istc;
-    return y;
-  }
-  else
-    return -1;
-}
-
-static void scAll( int Nvar, int deg)
-{
-  int i;
-  int d = deg;
-  if (d == 0)
-  {
-    for (i=Nvar; i; i--) act[i] = 0;
-    scElKbase();
-    return;
-  }
-  if (Nvar == 1)
-  {
-    act[1] = d;
-    scElKbase();
-    return;
-  }
-  do
-  {
-    act[Nvar] = d;
-    scAll(Nvar-1, deg-d);
-    d--;
-  } while (d >= 0);
-}
-
-static void scAllKbase( int Nvar, int ideg, int deg)
-{
-  do
-  {
-    act[Nvar] = ideg;
-    scAll(Nvar-1, deg-ideg);
-    ideg--;
-  } while (ideg >= 0);
-}
-
-static void scDegKbase( scfmon stc, int Nstc, int Nvar, int deg)
-{
-  int  Ivar, Istc, i, j;
-  scfmon sn;
-  int x, ideg;
-
-  if (deg == 0)
-  {
-    for (i=Nstc-1; i>=0; i--)
-    {
-      for (j=Nvar;j;j--){ if(stc[i][j]) break; }
-      if (j==0){return;}
-    }
-    for (i=Nvar; i; i--) act[i] = 0;
-    scElKbase();
-    return;
-  }
-  if (Nvar == 1)
-  {
-    for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
-    act[1] = deg;
-    scElKbase();
-    return;
-  }
-  Ivar = Nvar-1;
-  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
-  x = scRestrict(Nstc, sn, Nvar);
-  if (x <= 0)
-  {
-    if (x == 0) return;
-    ideg = deg;
-  }
-  else
-  {
-    if (deg < x) ideg = deg;
-    else ideg = x-1;
-    if (Nstc == 0)
-    {
-      scAllKbase(Nvar, ideg, deg);
-      return;
-    }
-  }
-  loop
-  {
-    x = scMax(Nstc, sn, Nvar);
-    while (ideg >= x)
-    {
-      act[Nvar] = ideg;
-      scDegKbase(sn, Nstc, Ivar, deg-ideg);
-      ideg--;
-    }
-    if (ideg < 0) return;
-    Istc = Nstc;
-    for (i=Nstc-1; i>=0; i--)
-    {
-      if (ideg < sn[i][Nvar])
-      {
-        Istc--;
-        sn[i] = NULL;
-      }
-    }
-    if (Istc == 0)
-    {
-      scAllKbase(Nvar, ideg, deg);
-      return;
-    }
-    j = 0;
-    while (sn[j]) j++;
-    i = j+1;
-    for (; i<Nstc; i++)
-    {
-      if (sn[i])
-      {
-        sn[j] = sn[i];
-        j++;
-      }
-    }
-    Nstc = Istc;
-  }
-}
-
-static void scInKbase( scfmon stc, int Nstc, int Nvar)
-{
-  int  Ivar, Istc, i, j;
-  scfmon sn;
-  int x, ideg;
-
-  if (Nvar == 1)
-  {
-    ideg = scMin(Nstc, stc, 1);
-    while (ideg > 0)
-    {
-      ideg--;
-      act[1] = ideg;
-      scElKbase();
-    }
-    return;
-  }
-  Ivar = Nvar-1;
-  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
-  x = scRestrict(Nstc, sn, Nvar);
-  if (x == 0) return;
-  ideg = x-1;
-  loop
-  {
-    x = scMax(Nstc, sn, Nvar);
-    while (ideg >= x)
-    {
-      act[Nvar] = ideg;
-      scInKbase(sn, Nstc, Ivar);
-      ideg--;
-    }
-    if (ideg < 0) return;
-    Istc = Nstc;
-    for (i=Nstc-1; i>=0; i--)
-    {
-      if (ideg < sn[i][Nvar])
-      {
-        Istc--;
-        sn[i] = NULL;
-      }
-    }
-    j = 0;
-    while (sn[j]) j++;
-    i = j+1;
-    for (; i<Nstc; i++)
-    {
-      if (sn[i])
-      {
-        sn[j] = sn[i];
-        j++;
-      }
-    }
-    Nstc = Istc;
-  }
-}
-
-static ideal scIdKbase()
-{
-  polyset mm;
-  ideal res;
-  poly p, q = last;
-  int i = pLength(q);
-  res = idInit(i,1);
-  mm = res->m;
-  i = 0;
-  do
-  {
-    mm[i] = q;
-    i++;
-    p = pNext(q);
-    pNext(q) = NULL;
-    q = p;
-  } while (q!=NULL);
-  return res;
-}
-
-ideal scKBase(int deg, ideal s, ideal Q, intvec * mv)
-{
-  int  i, di;
-  poly p;
-
-  if (deg < 0)
-  {
-    di = scDimInt(s, Q);
-    if (di != 0)
-    {
-      //Werror("KBase not finite");
-      return idInit(1,s->rank);
-    }
-  }
-  stcmem = hCreate((currRing->N) - 1);
-  hexist = hInit(s, Q, &hNexist);
-  p = last = pInit();
-  /*pNext(p) = NULL;*/
-  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
-  *act = 0;
-  if (!hNexist)
-  {
-    scAll((currRing->N), deg);
-    goto ende;
-  }
-  if (!hisModule)
-  {
-    if (deg < 0) scInKbase(hexist, hNexist, (currRing->N));
-    else scDegKbase(hexist, hNexist, (currRing->N), deg);
-  }
-  else
-  {
-    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
-    for (i = 1; i <= hisModule; i++)
-    {
-      *act = i;
-      hComp(hexist, hNexist, i, hstc, &hNstc);
-      int deg_ei=deg;
-      if (mv!=NULL) deg_ei -= (*mv)[i-1];
-      if ((deg < 0) || (deg_ei>=0))
-      {
-        if (hNstc)
-        {
-          if (deg < 0) scInKbase(hstc, hNstc, (currRing->N));
-          else scDegKbase(hstc, hNstc, (currRing->N), deg_ei);
-        }
-        else
-          scAll((currRing->N), deg_ei);
-      }
-    }
-    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
-  }
-ende:
-  hDelete(hexist, hNexist);
-  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
-  hKill(stcmem, (currRing->N) - 1);
-  pLmDelete(&p);
-  if (p == NULL)
-    return idInit(1,s->rank);
-  else
-  {
-    last = p;
-    ideal res=scIdKbase();
-    res->rank=s->rank;
-    return res;
-  }
-}
-
-#if 0 //-- alternative implementation of scComputeHC
-void scComputeHCw(ideal ss, ideal Q, int ak, poly &hEdge, ring tailRing)
-{
-  int  i, di;
-  poly p;
-
-  if (hEdge!=NULL)
-    pLmFree(hEdge);
-
-  ideal s=idInit(IDELEMS(ss),ak);
-  for(i=IDELEMS(ss)-1;i>=0;i--)
-  {
-    if (ss->m[i]!=NULL) s->m[i]=pHead(ss->m[i]);
-  }
-  di = scDimInt(s, Q);
-  stcmem = hCreate((currRing->N) - 1);
-  hexist = hInit(s, Q, &hNexist);
-  p = last = pInit();
-  /*pNext(p) = NULL;*/
-  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
-  *act = 0;
-  if (!hNexist)
-  {
-    scAll((currRing->N), -1);
-    goto ende;
-  }
-  if (!hisModule)
-  {
-    scInKbase(hexist, hNexist, (currRing->N));
-  }
-  else
-  {
-    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
-    for (i = 1; i <= hisModule; i++)
-    {
-      *act = i;
-      hComp(hexist, hNexist, i, hstc, &hNstc);
-      if (hNstc)
-      {
-        scInKbase(hstc, hNstc, (currRing->N));
-      }
-      else
-        scAll((currRing->N), -1);
-    }
-    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
-  }
-ende:
-  hDelete(hexist, hNexist);
-  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
-  hKill(stcmem, (currRing->N) - 1);
-  pDeleteLm(&p);
-  idDelete(&s);
-  if (p == NULL)
-  {
-    return; // no HEdge
-  }
-  else
-  {
-    last = p;
-    ideal res=scIdKbase();
-    res->rank=ss->rank;
-    poly p_ind=res->m[0]; int ind=0;
-    for(i=IDELEMS(res)-1;i>0;i--)
-    {
-      if (pCmp(res->m[i],p_ind)==-1) { p_ind=res->m[i]; ind=i; }
-    }
-    assume(p_ind!=NULL);
-    assume(res->m[ind]==p_ind);
-    hEdge=p_ind;
-    res->m[ind]=NULL;
-    nDelete(&pGetCoeff(hEdge));
-    pGetCoeff(hEdge)=NULL;
-    for(i=(currRing->N);i>0;i--)
-      pIncrExp(hEdge,i);
-    pSetm(hEdge);
-
-    idDelete(&res);
-    return;
-  }
-}
-#endif
Index: rnel/hilb.cc
===================================================================
--- kernel/hilb.cc	(revision 66b1de771706655cf8a9c3924d1dac08150c8814)
+++ 	(revision )
@@ -1,1387 +1,0 @@
-/****************************************
-*  Computer Algebra System SINGULAR     *
-****************************************/
-/*
-*  ABSTRACT -  Hilbert series
-*/
-
-#ifdef HAVE_CONFIG_H
-#include "singularconfig.h"
-#endif /* HAVE_CONFIG_H */
-#include <kernel/mod2.h>
-#include <kernel/structs.h>
-#include <kernel/febase.h>
-#include <omalloc/omalloc.h>
-#include <kernel/polys.h>
-#include <misc/intvec.h>
-#include <kernel/hutil.h>
-#include <kernel/stairc.h>
-//ADICHANGES:
-#include <kernel/ideals.h>
-#include <kernel/kstd1.h>
-#include<gmp.h>
-#include<vector>
-
-
-static int  **Qpol;
-static int  *Q0, *Ql;
-static int  hLength;
-   
-
-static int hMinModulweight(intvec *modulweight)
-{
-  int i,j,k;
-
-  if(modulweight==NULL) return 0;
-  j=(*modulweight)[0];
-  for(i=modulweight->rows()-1;i!=0;i--)
-  {
-    k=(*modulweight)[i];
-    if(k<j) j=k;
-  }
-  return j;
-}
-
-static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
-{
-  int  i, j;
-  int  x, y, z = 1;
-  int  *p;
-  for (i = Nvar; i>0; i--)
-  {
-    x = 0;
-    for (j = 0; j < Nstc; j++)
-    {
-      y = stc[j][var[i]];
-      if (y > x)
-        x = y;
-    }
-    z += x;
-    j = i - 1;
-    if (z > Ql[j])
-    {
-      if (z>(MAX_INT_VAL)/2)
-      {
-       Werror("interal arrays too big");
-       return;
-      }
-      p = (int *)omAlloc((unsigned long)z * sizeof(int));
-      if (Ql[j]!=0)
-      {
-        if (j==0)
-          memcpy(p, Qpol[j], Ql[j] * sizeof(int));
-        omFreeSize((ADDRESS)Qpol[j], Ql[j] * sizeof(int));
-      }
-      if (j==0)
-      {
-        for (x = Ql[j]; x < z; x++)
-          p[x] = 0;
-      }
-      Ql[j] = z;
-      Qpol[j] = p;
-    }
-  }
-}
-
-static int  *hAddHilb(int Nv, int x, int *pol, int *lp)
-{
-  int  l = *lp, ln, i;
-  int  *pon;
-  *lp = ln = l + x;
-  pon = Qpol[Nv];
-  memcpy(pon, pol, l * sizeof(int));
-  if (l > x)
-  {
-    for (i = x; i < l; i++)
-      pon[i] -= pol[i - x];
-    for (i = l; i < ln; i++)
-      pon[i] = -pol[i - x];
-  }
-  else
-  {
-    for (i = l; i < x; i++)
-      pon[i] = 0;
-    for (i = x; i < ln; i++)
-      pon[i] = -pol[i - x];
-  }
-  return pon;
-}
-
-static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
-{
-  int  l = lp, x, i, j;
-  int  *p, *pl;
-  p = pol;
-  for (i = Nv; i>0; i--)
-  {
-    x = pure[var[i + 1]];
-    if (x!=0)
-      p = hAddHilb(i, x, p, &l);
-  }
-  pl = *Qpol;
-  j = Q0[Nv + 1];
-  for (i = 0; i < l; i++)
-    pl[i + j] += p[i];
-  x = pure[var[1]];
-  if (x!=0)
-  {
-    j += x;
-    for (i = 0; i < l; i++)
-      pl[i + j] -= p[i];
-  }
-  j += l;
-  if (j > hLength)
-    hLength = j;
-}
-
-static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var,
- int Nvar, int *pol, int Lpol)
-{
-  int  iv = Nvar -1, ln, a, a0, a1, b, i;
-  int  x, x0;
-  scmon pn;
-  scfmon sn;
-  int  *pon;
-  if (Nstc==0)
-  {
-    hLastHilb(pure, iv, var, pol, Lpol);
-    return;
-  }
-  x = a = 0;
-  pn = hGetpure(pure);
-  sn = hGetmem(Nstc, stc, stcmem[iv]);
-  hStepS(sn, Nstc, var, Nvar, &a, &x);
-  Q0[iv] = Q0[Nvar];
-  ln = Lpol;
-  pon = pol;
-  if (a == Nstc)
-  {
-    x = pure[var[Nvar]];
-    if (x!=0)
-      pon = hAddHilb(iv, x, pon, &ln);
-    hHilbStep(pn, sn, a, var, iv, pon, ln);
-    return;
-  }
-  else
-  {
-    pon = hAddHilb(iv, x, pon, &ln);
-    hHilbStep(pn, sn, a, var, iv, pon, ln);
-  }
-  b = a;
-  x0 = 0;
-  loop
-  {
-    Q0[iv] += (x - x0);
-    a0 = a;
-    x0 = x;
-    hStepS(sn, Nstc, var, Nvar, &a, &x);
-    hElimS(sn, &b, a0, a, var, iv);
-    a1 = a;
-    hPure(sn, a0, &a1, var, iv, pn, &i);
-    hLex2S(sn, b, a0, a1, var, iv, hwork);
-    b += (a1 - a0);
-    ln = Lpol;
-    if (a < Nstc)
-    {
-      pon = hAddHilb(iv, x - x0, pol, &ln);
-      hHilbStep(pn, sn, b, var, iv, pon, ln);
-    }
-    else
-    {
-      x = pure[var[Nvar]];
-      if (x!=0)
-        pon = hAddHilb(iv, x - x0, pol, &ln);
-      else
-        pon = pol;
-      hHilbStep(pn, sn, b, var, iv, pon, ln);
-      return;
-    }
-  }
-}
-
-/*
-*basic routines
-*/
-static void hWDegree(intvec *wdegree)
-{
-  int i, k;
-  int x;
-
-  for (i=(currRing->N); i; i--)
-  {
-    x = (*wdegree)[i-1];
-    if (x != 1)
-    {
-      for (k=hNexist-1; k>=0; k--)
-      {
-        hexist[k][i] *= x;
-      }
-    }
-  }
-}
-// ---------------------------------- ADICHANGES ---------------------------------------------
-//!!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-//returns the degree of the monomial
-static int DegMon(poly p)
-{
-    #if 1
-    int i,deg;
-    deg = 0;
-    for(i=1;i<=currRing->N;i++)
-    {
-        deg = deg + p_GetExp(p, i, currRing);
-    }
-    return(deg);
-    #else
-    return(p_Deg(p, currRing));
-    #endif
-}
-
-//Tests if the ideal is sorted by degree
-static bool idDegSortTest(ideal I)
-{
-    if((I == NULL)||(idIs0(I)))
-    {
-        return(TRUE);
-    }
-    for(int i = 0; i<IDELEMS(I)-1; i++)
-    {
-        if(DegMon(I->m[i])>DegMon(I->m[i+1]))
-        {
-            idPrint(I);
-            Werror("Ideal is not deg sorted!!");
-            return(FALSE);
-        }
-    }
-    return(TRUE);
-}
-
-//adds the new polynomial at the coresponding position
-//and simplifies the ideal
-static ideal SortByDeg_p(ideal I, poly p)
-{
-    int i,j;
-    if((I == NULL) || (idIs0(I)))
-    {
-        ideal res = idInit(1,1);
-        res->m[0] = p;
-        return(res);
-    }
-    idSkipZeroes(I);
-    #if 1
-    for(i = 0; (i<IDELEMS(I)) && (DegMon(I->m[i])<=DegMon(p)); i++)
-    {
-        if(p_DivisibleBy( I->m[i],p, currRing))
-        {
-            return(I);
-        }
-    }
-    for(i = IDELEMS(I)-1; (i>=0) && (DegMon(I->m[i])>=DegMon(p)); i--)
-    {
-        if(p_DivisibleBy(p,I->m[i], currRing))
-        {
-            I->m[i] = NULL;
-        }
-    }
-    if(idIs0(I))
-    {
-        idSkipZeroes(I);
-        I->m[0] = p;
-        return(I);
-    }
-    #endif
-    idSkipZeroes(I);
-    //First I take the case when all generators have the same degree
-    if(DegMon(I->m[0]) == DegMon(I->m[IDELEMS(I)-1]))
-    {
-        if(DegMon(p)<DegMon(I->m[0]))
-        {
-            idInsertPoly(I,p);
-            idSkipZeroes(I);
-            for(i=IDELEMS(I)-1;i>=1; i--)
-            {
-                I->m[i] = I->m[i-1];
-            }
-            I->m[0] = p;
-            return(I);
-        }
-        if(DegMon(p)>=DegMon(I->m[IDELEMS(I)-1]))
-        {
-            idInsertPoly(I,p);
-            idSkipZeroes(I);
-            return(I);
-        }
-    }
-    if(DegMon(p)<=DegMon(I->m[0]))
-    {
-        idInsertPoly(I,p);
-        idSkipZeroes(I);
-        for(i=IDELEMS(I)-1;i>=1; i--)
-        {
-            I->m[i] = I->m[i-1];
-        }
-        I->m[0] = p;
-        return(I);
-    }
-    if(DegMon(p)>=DegMon(I->m[IDELEMS(I)-1]))
-    {
-        idInsertPoly(I,p);
-        idSkipZeroes(I);
-        return(I);
-    }
-    for(i = IDELEMS(I)-2; ;)
-    {
-        if(DegMon(p)==DegMon(I->m[i]))
-        {
-            idInsertPoly(I,p);
-            idSkipZeroes(I);
-            for(j = IDELEMS(I)-1; j>=i+1;j--)
-            {
-                I->m[j] = I->m[j-1];
-            }
-            I->m[i] = p;
-            return(I);
-        }
-        if(DegMon(p)>DegMon(I->m[i]))
-        {
-            idInsertPoly(I,p);
-            idSkipZeroes(I);
-            for(j = IDELEMS(I)-1; j>=i+2;j--)
-            {
-                I->m[j] = I->m[j-1];
-            }
-            I->m[i+1] = p;
-            return(I);
-        }
-        i--;
-    }
-}
-
-//it sorts the ideal by the degrees
-static ideal SortByDeg(ideal I)
-{
-    if(idIs0(I))
-    {
-        return(I);
-    }
-    idSkipZeroes(I);
-    int i;
-    ideal res;
-    idSkipZeroes(I);
-    res = idInit(1,1);
-    res->m[0] = poly(0);
-    for(i = 0; i<=IDELEMS(I)-1;i++)
-    {
-        res = SortByDeg_p(res, I->m[i]);
-    }
-    idSkipZeroes(res);
-    //idDegSortTest(res);
-    return(res);
-}
-
-//idQuot(I,p) for I monomial ideal, p a ideal with a single monomial.
-ideal idQuotMon(ideal Iorig, ideal p)
-{
-    if(idIs0(Iorig))
-    {
-        ideal res = idInit(1,1);
-        res->m[0] = poly(0);
-        return(res);
-    }
-    if(idIs0(p))
-    {
-        ideal res = idInit(1,1);
-        res->m[0] = pOne();
-        return(res);
-    }
-    ideal I = idCopy(Iorig);
-    ideal res = idInit(IDELEMS(I),1);
-    int i,j;
-    int dummy;
-    for(i = 0; i<IDELEMS(I); i++)
-    {
-        res->m[i] = p_Copy(I->m[i], currRing);
-        for(j = 1; (j<=currRing->N) ; j++)
-        {
-            dummy = p_GetExp(p->m[0], j, currRing);
-            if(dummy > 0)
-            {
-                if(p_GetExp(I->m[i], j, currRing) < dummy)
-                {
-                    p_SetExp(res->m[i], j, 0, currRing);
-                }
-                else
-                {
-                    p_SetExp(res->m[i], j, p_GetExp(I->m[i], j, currRing) - dummy, currRing);
-                }
-            }
-        }
-        p_Setm(res->m[i], currRing);
-        if(DegMon(res->m[i]) == DegMon(I->m[i]))
-        {
-            res->m[i] = NULL;
-        }
-        else
-        {
-            I->m[i] = NULL;
-        }
-    }
-    idSkipZeroes(res);
-    idSkipZeroes(I);
-    if(!idIs0(res))
-    {
-    for(i = 0; i<=IDELEMS(res)-1; i++)
-    {
-        I = SortByDeg_p(I,res->m[i]);
-    }
-    }
-    //idDegSortTest(I);
-    return(I);
-}
-
-//id_Add for monomials
-static ideal idAddMon(ideal I, ideal p)
-{
-    #if 1
-    I = SortByDeg_p(I,p->m[0]);
-    #else
-    I = id_Add(I,p,currRing);
-    #endif
-    //idSkipZeroes(I);
-    return(I);
-}
-
-//searches for a variable that is not yet used (assumes that the ideal is sqrfree)
-static poly ChoosePVar (ideal I)
-{
-    bool flag=TRUE;
-    int i,j;
-    poly res;
-    for(i=1;i<=currRing->N;i++)
-    {
-        flag=TRUE;
-        for(j=IDELEMS(I)-1;(j>=0)&&(flag);j--)
-        {
-            if(p_GetExp(I->m[j], i, currRing)>0)
-            {
-                flag=FALSE;
-            }
-        }
-        
-        if(flag == TRUE)
-        {
-            res = p_ISet(1, currRing);
-            p_SetExp(res, i, 1, currRing);
-            p_Setm(res,currRing);
-            return(res);
-        }
-        else
-        {
-            p_Delete(&res, currRing);
-        }
-    }
-    return(NULL); //i.e. it is the maximal ideal
-}
-
-//choice XL: last entry divided by x (xy10z15 -> y9z14)
-static poly ChoosePXL(ideal I)
-{
-    int i,j,dummy=0;
-    poly m;
-    for(i = IDELEMS(I)-1; (i>=0) && (dummy == 0); i--)
-    {
-        for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
-        {
-            if(p_GetExp(I->m[i],j, currRing)>1)
-            {
-                dummy = 1;
-            }
-        }
-    }
-    m = p_Copy(I->m[i+1],currRing);
-    for(j = 1; j<=currRing->N; j++)
-    {
-        dummy = p_GetExp(m,j,currRing);
-        if(dummy >= 1)
-        {
-            p_SetExp(m, j, dummy-1, currRing);
-        }
-    }
-    if(!p_IsOne(m, currRing))
-    {
-        p_Setm(m, currRing);
-        return(m);
-    }
-    m = ChoosePVar(I);
-    return(m);
-}
-
-//choice XF: first entry divided by x (xy10z15 -> y9z14)
-static poly ChoosePXF(ideal I)
-{
-    int i,j,dummy=0;
-    poly m;
-    for(i =0 ; (i<=IDELEMS(I)-1) && (dummy == 0); i++)
-    {
-        for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
-        {
-            if(p_GetExp(I->m[i],j, currRing)>1)
-            {
-                dummy = 1;
-            }
-        }
-    }
-    m = p_Copy(I->m[i-1],currRing);
-    for(j = 1; j<=currRing->N; j++)
-    {
-        dummy = p_GetExp(m,j,currRing);
-        if(dummy >= 1)
-        {
-            p_SetExp(m, j, dummy-1, currRing);
-        }
-    }
-    if(!p_IsOne(m, currRing))
-    {
-        p_Setm(m, currRing);
-        return(m);
-    }
-    m = ChoosePVar(I);
-    return(m);
-}
-
-//choice OL: last entry the first power (xy10z15 -> xy9z15)
-static poly ChoosePOL(ideal I)
-{
-    int i,j,dummy;
-    poly m;
-    for(i = IDELEMS(I)-1;i>=0;i--)
-    {
-        m = p_Copy(I->m[i],currRing);
-        for(j=1;j<=currRing->N;j++)
-        {
-            dummy = p_GetExp(m,j,currRing);
-            if(dummy > 0)
-            {
-                p_SetExp(m,j,dummy-1,currRing);
-                p_Setm(m,currRing);
-            }       
-        }
-        if(!p_IsOne(m, currRing))
-        {
-            return(m);
-        }
-        else
-        {
-            p_Delete(&m,currRing);
-        }
-    }
-    m = ChoosePVar(I);
-    return(m);
-}
-
-//choice OF: first entry the first power (xy10z15 -> xy9z15)
-static poly ChoosePOF(ideal I)
-{
-    int i,j,dummy;
-    poly m;
-    for(i = 0 ;i<=IDELEMS(I)-1;i++)
-    {
-        m = p_Copy(I->m[i],currRing);
-        for(j=1;j<=currRing->N;j++)
-        {
-            dummy = p_GetExp(m,j,currRing);
-            if(dummy > 0)
-            {
-                p_SetExp(m,j,dummy-1,currRing);
-                p_Setm(m,currRing);
-            }       
-        }
-        if(!p_IsOne(m, currRing))
-        {
-            return(m);
-        }
-        else
-        {
-            p_Delete(&m,currRing);
-        }
-    }
-    m = ChoosePVar(I);
-    return(m);
-}
-
-//choice VL: last entry the first variable with power (xy10z15 -> y)
-static poly ChoosePVL(ideal I)
-{
-    int i,j,dummy;
-    bool flag = TRUE;
-    poly m = p_ISet(1,currRing);
-    for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
-    {
-        flag = TRUE;
-        for(j=1;(j<=currRing->N) && (flag);j++)
-        {
-            dummy = p_GetExp(I->m[i],j,currRing);
-            if(dummy >= 2)
-            {
-                p_SetExp(m,j,1,currRing);
-                p_Setm(m,currRing);
-                flag = FALSE;
-            }       
-        }
-        if(!p_IsOne(m, currRing))
-        {
-            return(m);
-        }
-    }
-    m = ChoosePVar(I);
-    return(m);
-}
-
-//choice VF: first entry the first variable with power (xy10z15 -> y)
-static poly ChoosePVF(ideal I)
-{
-    int i,j,dummy;
-    bool flag = TRUE;
-    poly m = p_ISet(1,currRing);
-    for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
-    {
-        flag = TRUE;
-        for(j=1;(j<=currRing->N) && (flag);j++)
-        {
-            dummy = p_GetExp(I->m[i],j,currRing);
-            if(dummy >= 2)
-            {
-                p_SetExp(m,j,1,currRing);
-                p_Setm(m,currRing);
-                flag = FALSE;
-            }       
-        }
-        if(!p_IsOne(m, currRing))
-        {
-            return(m);
-        }
-    }
-    m = ChoosePVar(I);
-    return(m);
-}
-
-//choice JL: last entry just variable with power (xy10z15 -> y10)
-static poly ChoosePJL(ideal I)
-{
-    int i,j,dummy;
-    bool flag = TRUE;
-    poly m = p_ISet(1,currRing);
-    for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
-    {
-        flag = TRUE;
-        for(j=1;(j<=currRing->N) && (flag);j++)
-        {
-            dummy = p_GetExp(I->m[i],j,currRing);
-            if(dummy >= 2)
-            {
-                p_SetExp(m,j,dummy-1,currRing);
-                p_Setm(m,currRing);
-                flag = FALSE;
-            }       
-        }
-        if(!p_IsOne(m, currRing))
-        {
-            return(m);
-        }
-    }
-    m = ChoosePVar(I);
-    return(m);
-}
-
-//choice JF: last entry just variable with power -1 (xy10z15 -> y9)
-static poly ChoosePJF(ideal I)
-{
-    int i,j,dummy;
-    bool flag = TRUE;
-    poly m = p_ISet(1,currRing);
-    for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
-    {
-        flag = TRUE;
-        for(j=1;(j<=currRing->N) && (flag);j++)
-        {
-            dummy = p_GetExp(I->m[i],j,currRing);
-            if(dummy >= 2)
-            {
-                p_SetExp(m,j,dummy-1,currRing);
-                p_Setm(m,currRing);
-                flag = FALSE;
-            }       
-        }
-        if(!p_IsOne(m, currRing))
-        {
-            return(m);
-        }
-    }
-    m = ChoosePVar(I);
-    return(m);
-}
-
-//chooses 1 \neq p \not\in S. This choice should be made optimal
-static poly ChooseP(ideal I)
-{
-    poly m;
-    //  TEST TO SEE WHICH ONE IS BETTER
-    //m = ChoosePXL(I);
-    //m = ChoosePXF(I);
-    //m = ChoosePOL(I);
-    //m = ChoosePOF(I);
-    //m = ChoosePVL(I);
-    //m = ChoosePVF(I);
-    m = ChoosePJL(I);
-    //m = ChoosePJF(I);
-    return(m);
-}
-
-//searches for a monomial of degree d>=2 and divides it by a variable (result monomial of deg d-1)
-static poly SearchP(ideal I)
-{
-    int i,j,exp;
-    poly res;
-    if(DegMon(I->m[IDELEMS(I)-1])<=1)
-    {
-        res = ChoosePVar(I);
-        return(res);
-    }
-    i = IDELEMS(I)-1;
-    res = p_Copy(I->m[i], currRing);
-    for(j=1;j<=currRing->N;j++)
-    {
-        exp = p_GetExp(I->m[i], j, currRing);
-        if(exp > 0)
-        {
-            p_SetExp(res, j, exp - 1, currRing);
-            p_Setm(res,currRing);
-            return(res);
-        }
-    }
-}
-
-//test if the ideal is of the form (x1, ..., xr)
-static bool JustVar(ideal I)
-{
-    #if 0
-    int i,j;
-    bool foundone;
-    for(i=0;i<=IDELEMS(I)-1;i++)
-    {
-        foundone = FALSE;
-        for(j = 1;j<=currRing->N;j++)
-        {
-            if(p_GetExp(I->m[i], j, currRing)>0)
-            {
-                if(foundone == TRUE)
-                {
-                    return(FALSE);
-                }
-                foundone = TRUE;
-            }
-        }        
-    }
-    return(TRUE);
-    #else
-    if(DegMon(I->m[IDELEMS(I)-1])>1)
-    {
-        return(FALSE);
-    }
-    return(TRUE);
-    #endif
-}
-
-//computes the Euler Characteristic of the ideal
-static void eulerchar (ideal I, int variables, mpz_ptr ec)
-{
-  loop
-  {
-    mpz_t dummy;
-    if(JustVar(I) == TRUE)
-    {
-        if(IDELEMS(I) == variables)
-        {
-            mpz_init(dummy);
-            if((variables % 2) == 0)
-                {mpz_set_si(dummy, 1);}
-            else
-                {mpz_set_si(dummy, -1);}
-            mpz_add(ec, ec, dummy);
-        }
-        //mpz_clear(dummy);
-        return;        
-    }
-    ideal p = idInit(1,1);
-    p->m[0] = SearchP(I);
-    //idPrint(I);
-    //idPrint(p);
-    //printf("\nNow get in idQuotMon\n");
-    ideal Ip = idQuotMon(I,p);
-    //idPrint(Ip);
-    //Ip = SortByDeg(Ip);
-    int i,howmanyvarinp = 0;
-    for(i = 1;i<=currRing->N;i++)
-    {
-        if(p_GetExp(p->m[0],i,currRing)>0)
-        {
-            howmanyvarinp++;
-        }
-    }    
-    eulerchar(Ip, variables-howmanyvarinp, ec);
-    id_Delete(&Ip, currRing);
-    I = idAddMon(I,p);
-  }
-}
-
-//tests if an ideal is Square Free, if no, returns the variable which appears at powers >1
-static poly SqFree (ideal I)
-{
-    int i,j;
-    bool flag=TRUE;
-    poly notsqrfree = NULL;
-    if(DegMon(I->m[IDELEMS(I)-1])<=1)
-    {
-        return(notsqrfree);
-    }
-    for(i=IDELEMS(I)-1;(i>=0)&&(flag);i--)
-    {
-        for(j=1;(j<=currRing->N)&&(flag);j++)
-        {
-            if(p_GetExp(I->m[i],j,currRing)>1)
-            {
-                flag=FALSE;
-                notsqrfree = p_ISet(1,currRing);
-                p_SetExp(notsqrfree,j,1,currRing);
-            }
-        }
-    }
-    if(notsqrfree != NULL)
-    {
-        p_Setm(notsqrfree,currRing);
-    }
-    return(notsqrfree);
-}
-
-//checks if a polynomial is in I
-static bool IsIn(poly p, ideal I)
-{
-    //assumes that I is ordered by degree
-    if(idIs0(I))
-    {
-        if(p==poly(0))
-        {
-            return(TRUE);
-        }
-        else
-        {
-            return(FALSE);
-        }
-    }
-    if(p==poly(0))
-    {
-        return(FALSE);
-    }
-    int i,j;
-    bool flag;
-    for(i = 0;i<IDELEMS(I);i++)
-    {
-        flag = TRUE;
-        for(j = 1;(j<=currRing->N) &&(flag);j++)
-        {
-            if(p_GetExp(p, j, currRing)<p_GetExp(I->m[i], j, currRing))
-            {
-                flag = FALSE;
-            }
-        }
-        if(flag)
-        {
-            return(TRUE);
-        }
-    }
-    return(FALSE);
-}
-
-//computes the lcm of min I, I monomial ideal
-static poly LCMmon(ideal I)
-{
-    if(idIs0(I))
-    {
-        return(NULL);
-    }
-    poly m;
-    int dummy,i,j;
-    m = p_ISet(1,currRing);
-    for(i=1;i<=currRing->N;i++)
-    {
-        dummy=0;
-        for(j=IDELEMS(I)-1;j>=0;j--)
-        {
-            if(p_GetExp(I->m[j],i,currRing) > dummy)
-            {
-                dummy = p_GetExp(I->m[j],i,currRing);
-            }
-        }
-        p_SetExp(m,i,dummy,currRing);
-    }
-    p_Setm(m,currRing);
-    return(m);
-}
-
-//the Roune Slice Algorithm
-void rouneslice(ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int* &hilbpower)
-{
-  loop
-  {
-    (steps)++;
-    int i,j;
-    int dummy;
-    poly m;
-    ideal p, koszsimp;
-    //----------- PRUNING OF S ---------------
-    //S SHOULD IN THIS POINT BE ORDERED BY DEGREE
-    for(i=IDELEMS(S)-1;i>=0;i--)
-    {
-        if(IsIn(S->m[i],I))
-        {
-            S->m[i]=NULL;
-            prune++;
-        }
-    }
-    idSkipZeroes(S);
-    //----------------------------------------
-    for(i=IDELEMS(I)-1;i>=0;i--)
-    {
-        m = p_Copy(I->m[i],currRing);
-        for(j=1;j<=currRing->N;j++)
-        {
-            dummy = p_GetExp(m,j,currRing);
-            if(dummy > 0)
-            {
-                p_SetExp(m,j,dummy-1,currRing);
-            }       
-        }
-        p_Setm(m, currRing);
-        if(IsIn(m,S))
-        {
-            I->m[i]=NULL;
-            //printf("\n Deleted, since pi(m) is in S\n");pWrite(m);
-        }
-    }
-    idSkipZeroes(I);
-    //----------- MORE PRUNING OF S ------------
-    m = LCMmon(I); 
-    if(m != NULL)
-    {
-        for(i=0;i<IDELEMS(S);i++)
-        {      
-            if(!(p_DivisibleBy(S->m[i], m, currRing)))  
-            {
-                S->m[i] = NULL;
-                j++;
-                moreprune++;
-            }
-            else
-            {
-                if(pLmEqual(S->m[i],m))
-                {
-                    S->m[i] = NULL;
-                    moreprune++;
-                }
-            }
-        }
-    idSkipZeroes(S);
-    }
-    /*printf("\n---------------------------\n");
-    printf("\n      I\n");idPrint(I);
-    printf("\n      S\n");idPrint(S);
-    printf("\n      q\n");pWrite(q);
-    getchar();*/
-    
-    if(idIs0(I))
-    {
-        id_Delete(&I, currRing);
-        id_Delete(&S, currRing);
-        p_Delete(&m, currRing);
-        break;
-    }
-    m = LCMmon(I);
-    if(!p_DivisibleBy(x,m, currRing))
-    {
-        //printf("\nx does not divide lcm(I)");
-        //printf("\nEmpty set");pWrite(q);
-        id_Delete(&I, currRing);
-        id_Delete(&S, currRing);
-        p_Delete(&m, currRing);
-        break;
-    }
-    m = SqFree(I);
-    if(m==NULL)
-    {
-        //printf("\n      Corner: ");
-        //pWrite(q);
-        //printf("\n      With the facets of the dual simplex:\n");
-        //idPrint(I);
-        mpz_t ec;
-        mpz_init(ec);
-        mpz_ptr ec_ptr = ec;
-        eulerchar(I, currRing->N, ec_ptr);
-        bool flag = FALSE;
-        if(NNN==0)
-            {
-                hilbertcoef = (mpz_ptr)omAlloc((NNN+1)*sizeof(mpz_t));
-                hilbpower = (int*)omAlloc((NNN+1)*sizeof(int));
-                mpz_init( &hilbertcoef[NNN]);
-                mpz_set(  &hilbertcoef[NNN], ec);
-                mpz_clear(ec);
-                hilbpower[NNN] = DegMon(q);
-                NNN++;
-            }
-        else
-        {
-            //I look if the power appears already
-            for(i = 0;(i<NNN)&&(flag == FALSE)&&(DegMon(q)>=hilbpower[i]);i++)
-            {
-                if((hilbpower[i]) == (DegMon(q)))
-                {
-                    flag = TRUE;
-                    mpz_add(&hilbertcoef[i],&hilbertcoef[i],ec_ptr);
-                }
-            }
-            if(flag == FALSE)
-            {
-                hilbertcoef = (mpz_ptr)omRealloc(hilbertcoef, (NNN+1)*sizeof(mpz_t));
-                hilbpower = (int*)omRealloc(hilbpower, (NNN+1)*sizeof(int));
-                mpz_init(&hilbertcoef[NNN]);
-                for(j = NNN; j>i; j--)
-                {
-                    mpz_set(&hilbertcoef[j],&hilbertcoef[j-1]);
-                    hilbpower[j] = hilbpower[j-1];
-                }
-                mpz_set(  &hilbertcoef[i], ec);
-                mpz_clear(ec);
-                hilbpower[i] = DegMon(q);
-                NNN++;
-            }
-        }
-        break;
-    }
-    m = ChooseP(I);
-    p = idInit(1,1);
-    p->m[0] = m;
-    ideal Ip = idQuotMon(I,p);
-    ideal Sp = idQuotMon(S,p);
-    poly pq = pp_Mult_mm(q,m,currRing);
-    rouneslice(Ip, Sp, pq, x, prune, moreprune, steps, NNN, hilbertcoef,hilbpower);
-    //id_Delete(&Ip, currRing);
-    //id_Delete(&Sp, currRing);
-    S = idAddMon(S,p);
-    p->m[0]=NULL; 
-    id_Delete(&p, currRing); // p->m[0] was also in S
-    p_Delete(&pq,currRing);
-  }
-}
-
-//it computes the first hilbert series by means of Roune Slice Algorithm
-void slicehilb(ideal I)
-{
-    //printf("Adi changes are here: \n");
-    int i, NNN = 0;
-    int steps = 0, prune = 0, moreprune = 0;
-    mpz_ptr hilbertcoef;
-    int *hilbpower;
-    ideal S = idInit(1,1);
-    poly q = p_ISet(1,currRing);
-    ideal X = idInit(1,1);
-    X->m[0]=p_One(currRing);
-    for(i=1;i<=currRing->N;i++)
-    {
-            p_SetExp(X->m[0],i,1,currRing);   
-    }
-    p_Setm(X->m[0],currRing);
-    I = id_Mult(I,X,currRing);
-    I = SortByDeg(I);
-    //printf("\n-------------RouneSlice--------------\n");
-    rouneslice(I,S,q,X->m[0],prune, moreprune, steps, NNN, hilbertcoef, hilbpower);
-    //printf("\nIn total Prune got rid of %i elements\n",prune);
-    //printf("\nIn total More Prune got rid of %i elements\n",moreprune);
-    //printf("\nSteps of rouneslice: %i\n\n", steps);
-    mpz_t coefhilb;
-    mpz_t dummy;
-    mpz_init(coefhilb);
-    mpz_init(dummy);
-    printf("\n//  %8d t^0",1);
-    for(i = 0; i<NNN; i++)
-    {
-        if(mpz_sgn(&hilbertcoef[i])!=0)
-        {
-            gmp_printf("\n//  %8Zd t^%d",&hilbertcoef[i],hilbpower[i]);
-        }
-    }
-    omFreeSize(hilbertcoef, (NNN)*sizeof(mpz_t));
-    omFreeSize(hilbpower, (NNN)*sizeof(int));
-    //printf("\n-------------------------------------\n");
-}
-
-// -------------------------------- END OF CHANGES -------------------------------------------
-static intvec * hSeries(ideal S, intvec *modulweight,
-                int /*notstc*/, intvec *wdegree, ideal Q, ring tailRing)
-{
-  intvec *work, *hseries1=NULL;
-  int  mc;
-  int  p0;
-  int  i, j, k, l, ii, mw;
-  hexist = hInit(S, Q, &hNexist, tailRing);
-  if (hNexist==0)
-  {
-    hseries1=new intvec(2);
-    (*hseries1)[0]=1;
-    (*hseries1)[1]=0;
-    return hseries1;
-  }
-
-  #if 0
-  if (wdegree == NULL)
-    hWeight();
-  else
-    hWDegree(wdegree);
-  #else
-  if (wdegree != NULL) hWDegree(wdegree);
-  #endif
-
-  p0 = 1;
-  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
-  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  stcmem = hCreate((currRing->N) - 1);
-  Qpol = (int **)omAlloc(((currRing->N) + 1) * sizeof(int *));
-  Ql = (int *)omAlloc0(((currRing->N) + 1) * sizeof(int));
-  Q0 = (int *)omAlloc(((currRing->N) + 1) * sizeof(int));
-  *Qpol = NULL;
-  hLength = k = j = 0;
-  mc = hisModule;
-  if (mc!=0)
-  {
-    mw = hMinModulweight(modulweight);
-    hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
-  }
-  else
-  {
-    mw = 0;
-    hstc = hexist;
-    hNstc = hNexist;
-  }
-  loop
-  {
-    if (mc!=0)
-    {
-      hComp(hexist, hNexist, mc, hstc, &hNstc);
-      if (modulweight != NULL)
-        j = (*modulweight)[mc-1]-mw;
-    }
-    if (hNstc!=0)
-    {
-      hNvar = (currRing->N);
-      for (i = hNvar; i>=0; i--)
-        hvar[i] = i;
-      //if (notstc) // TODO: no mon divides another
-        hStaircase(hstc, &hNstc, hvar, hNvar);
-      hSupp(hstc, hNstc, hvar, &hNvar);
-      if (hNvar!=0)
-      {
-        if ((hNvar > 2) && (hNstc > 10))
-          hOrdSupp(hstc, hNstc, hvar, hNvar);
-        hHilbEst(hstc, hNstc, hvar, hNvar);
-        memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
-        hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
-        hLexS(hstc, hNstc, hvar, hNvar);
-        Q0[hNvar] = 0;
-        hHilbStep(hpure, hstc, hNstc, hvar, hNvar, &p0, 1);
-      }
-    }
-    else
-    {
-      if(*Qpol!=NULL)
-        (**Qpol)++;
-      else
-      {
-        *Qpol = (int *)omAlloc(sizeof(int));
-        hLength = *Ql = **Qpol = 1;
-      }
-    }
-    if (*Qpol!=NULL)
-    {
-      i = hLength;
-      while ((i > 0) && ((*Qpol)[i - 1] == 0))
-        i--;
-      if (i > 0)
-      {
-        l = i + j;
-        if (l > k)
-        {
-          work = new intvec(l);
-          for (ii=0; ii<k; ii++)
-            (*work)[ii] = (*hseries1)[ii];
-          if (hseries1 != NULL)
-            delete hseries1;
-          hseries1 = work;
-          k = l;
-        }
-        while (i > 0)
-        {
-          (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
-          (*Qpol)[i - 1] = 0;
-          i--;
-        }
-      }
-    }
-    mc--;
-    if (mc <= 0)
-      break;
-  }
-  if (k==0)
-  {
-    hseries1=new intvec(2);
-    (*hseries1)[0]=0;
-    (*hseries1)[1]=0;
-  }
-  else
-  {
-    l = k+1;
-    while ((*hseries1)[l-2]==0) l--;
-    if (l!=k)
-    {
-      work = new intvec(l);
-      for (ii=l-2; ii>=0; ii--)
-        (*work)[ii] = (*hseries1)[ii];
-      delete hseries1;
-      hseries1 = work;
-    }
-    (*hseries1)[l-1] = mw;
-  }
-  for (i = 0; i <= (currRing->N); i++)
-  {
-    if (Ql[i]!=0)
-      omFreeSize((ADDRESS)Qpol[i], Ql[i] * sizeof(int));
-  }
-  omFreeSize((ADDRESS)Q0, ((currRing->N) + 1) * sizeof(int));
-  omFreeSize((ADDRESS)Ql, ((currRing->N) + 1) * sizeof(int));
-  omFreeSize((ADDRESS)Qpol, ((currRing->N) + 1) * sizeof(int *));
-  hKill(stcmem, (currRing->N) - 1);
-  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
-  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
-  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
-  hDelete(hexist, hNexist);
-  if (hisModule!=0)
-    omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
-  return hseries1;
-}
-
-
-intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
-{
-  return hSeries(S, modulweight, 0, wdegree, Q, tailRing);
-}
-
-intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
-{
-  return hSeries(S, modulweight, 1, wdegree, Q, tailRing);
-}
-
-intvec * hSecondSeries(intvec *hseries1)
-{
-  intvec *work, *hseries2;
-  int i, j, k, s, t, l;
-  if (hseries1 == NULL)
-    return NULL;
-  work = new intvec(hseries1);
-  k = l = work->length()-1;
-  s = 0;
-  for (i = k-1; i >= 0; i--)
-    s += (*work)[i];
-  loop
-  {
-    if ((s != 0) || (k == 1))
-      break;
-    s = 0;
-    t = (*work)[k-1];
-    k--;
-    for (i = k-1; i >= 0; i--)
-    {
-      j = (*work)[i];
-      (*work)[i] = -t;
-      s += t;
-      t += j;
-    }
-  }
-  hseries2 = new intvec(k+1);
-  for (i = k-1; i >= 0; i--)
-    (*hseries2)[i] = (*work)[i];
-  (*hseries2)[k] = (*work)[l];
-  delete work;
-  return hseries2;
-}
-
-void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
-{
-  int m, i, j, k;
-  *co = *mu = 0;
-  if ((s1 == NULL) || (s2 == NULL))
-    return;
-  i = s1->length();
-  j = s2->length();
-  if (j > i)
-    return;
-  m = 0;
-  for(k=j-2; k>=0; k--)
-    m += (*s2)[k];
-  *mu = m;
-  *co = i - j;
-}
-
-static void hPrintHilb(intvec *hseries)
-{
-  int  i, j, l, k;
-  if (hseries == NULL)
-    return;
-  l = hseries->length()-1;
-  k = (*hseries)[l];
-  for (i = 0; i < l; i++)
-  {
-    j = (*hseries)[i];
-    if (j != 0)
-    {
-      Print("//  %8d t^%d\n", j, i+k);
-    }
-  }
-}
-
-/*
-*caller
-*/
-void hLookSeries(ideal S, intvec *modulweight, ideal Q)
-{
-  int co, mu, l;
-  intvec *hseries2;
-  intvec *hseries1 = hFirstSeries(S, modulweight, Q);
-  hPrintHilb(hseries1);
-  l = hseries1->length()-1;
-  if (l > 1)
-    hseries2 = hSecondSeries(hseries1);
-  else
-    hseries2 = hseries1;
-  hDegreeSeries(hseries1, hseries2, &co, &mu);
-  PrintLn();
-  hPrintHilb(hseries2);
-  if ((l == 1) &&(mu == 0))
-    scPrintDegree((currRing->N)+1, 0);
-  else
-    scPrintDegree(co, mu);
-  if (l>1)
-    delete hseries1;
-  delete hseries2;
-}
-
-
-
Index: rnel/hutil.cc
===================================================================
--- kernel/hutil.cc	(revision 66b1de771706655cf8a9c3924d1dac08150c8814)
+++ 	(revision )
@@ -1,1062 +1,0 @@
-/****************************************
-*  Computer Algebra System SINGULAR     *
-****************************************/
-/*
-* ABSTRACT: Utilities for staircase operations
-*/
-
-#ifdef HAVE_CONFIG_H
-#include "singularconfig.h"
-#endif /* HAVE_CONFIG_H */
-#include <kernel/mod2.h>
-#include <kernel/structs.h>
-#include <kernel/febase.h>
-#include <omalloc/omalloc.h>
-#include <kernel/ideals.h>
-#include <kernel/polys.h>
-#include <kernel/hutil.h>
-
-scfmon hexist, hstc, hrad, hwork;
-scmon hpure, hpur0;
-varset hvar, hsel;
-int  hNexist, hNstc, hNrad, hNvar, hNpure;
-int hisModule;
-monf stcmem, radmem;
-
-// Making a global "security" copy of the allocated exponent vectors
-// is a dirty fix for correct memory disallocation: It would be
-// better, if either the fields of heist are never touched
-// (i.e. changed) except in hInit, or, if hInit would return the
-// "security" copy as well. But then, all the relevant data is held in
-// global variables, so we might do that here, as well.
-static scfmon hsecure= NULL;
-
-scfmon hInit(ideal S, ideal Q, int *Nexist, ring tailRing)
-{
-  int  sl, ql, i, k = 0;
-  polyset si, qi, ss;
-  scfmon ex, ek;
-  if (tailRing != currRing)
-    hisModule = id_RankFreeModule(S, currRing, tailRing);
-  else
-    hisModule = id_RankFreeModule(S,currRing);
-  if (hisModule < 0)
-    hisModule = 0;
-  if (S!=NULL)
-  {
-    si = S->m;
-    sl = IDELEMS(S);
-  }
-  else
-  {
-    si = NULL;
-    sl = 0;
-  }
-  if (Q!=NULL)
-  {
-    qi = Q->m;
-    ql = IDELEMS(Q);
-  }
-  else
-  {
-    qi = NULL;
-    ql = 0;
-  }
-  if ((sl + ql) == 0)
-  {
-    *Nexist = 0;
-    return NULL;
-  }
-  ss = si;
-  for (i = sl; i>0; i--)
-  {
-    if (*ss!=0)
-      k++;
-    ss++;
-  }
-  ss = qi;
-  for (i = ql; i>0; i--)
-  {
-    if (*ss!=0)
-      k++;
-    ss++;
-  }
-  *Nexist = k;
-  if (k==0)
-    return NULL;
-  ek = ex = (scfmon)omAlloc0(k * sizeof(scmon));
-  hsecure = (scfmon) omAlloc0(k * sizeof(scmon));
-  for (i = sl; i>0; i--)
-  {
-    if (*si!=NULL)
-    {
-      *ek = (scmon) omAlloc(((currRing->N)+1)*sizeof(int));
-      pGetExpV(*si, *ek);
-      ek++;
-    }
-    si++;
-  }
-  for (i = ql; i>0; i--)
-  {
-    if (*qi!=NULL)
-    {
-      *ek = (scmon) omAlloc(((currRing->N)+1)*sizeof(int));
-      pGetExpV(*qi, *ek);
-      ek++;
-    }
-    qi++;
-  }
-  memcpy(hsecure, ex, k * sizeof(scmon));
-  return ex;
-}
-
-#if 0
-void hWeight()
-{
-  int i, k;
-  int x;
-
-  i = (currRing->N);
-  loop
-  {
-    if (pWeight(i) != 1) break;
-    i--;
-    if (i == 0) return;
-  }
-  for (i=(currRing->N); i>0; i--)
-  {
-    x = pWeight(i);
-    if (x != 1)
-    {
-      for (k=hNexist-1; k>=0; k--)
-      {
-        hexist[k][i] *= x;
-      }
-    }
-  }
-}
-#endif
-
-void hDelete(scfmon ev, int ev_length)
-{
-  int i;
-
-  if (ev_length>0)
-  {
-    for (i=ev_length-1;i>=0;i--)
-      omFreeSize(hsecure[i],((currRing->N)+1)*sizeof(int));
-    omFreeSize(hsecure, ev_length*sizeof(scmon));
-    omFreeSize(ev,  ev_length*sizeof(scmon));
-  }
-}
-
-
-void hComp(scfmon exist, int Nexist, int ak, scfmon stc, int *Nstc)
-{
-  int k = 0;
-  scfmon ex = exist, co = stc;
-  int i;
-
-  for (i = Nexist; i>0; i--)
-  {
-    if (((**ex) == 0) || ((**ex) == ak))
-    {
-      *co = *ex;
-      co++;
-      k++;
-    }
-    ex++;
-  }
-  *Nstc = k;
-}
-
-
-void hSupp(scfmon stc, int Nstc, varset var, int *Nvar)
-{
-  int  nv, i0, i1, i, j;
-  nv = i0 = *Nvar;
-  i1 = 0;
-  for (i = 1; i <= nv; i++)
-  {
-    j = 0;
-    loop
-    {
-      if (stc[j][i]>0)
-      {
-        i1++;
-        var[i1] = i;
-        break;
-      }
-      j++;
-      if (j == Nstc)
-      {
-        var[i0] = i;
-        i0--;
-        break;
-      }
-    }
-  }
-  *Nvar = i1;
-}
-
-void hOrdSupp(scfmon stc, int Nstc, varset var, int Nvar)
-{
-  int  i, i1, j, jj, k, l;
-  int  x;
-  scmon temp, count;
-  float o, h, g, *v1;
-
-  v1 = (float *)omAlloc(Nvar * sizeof(float));
-  temp = (int *)omAlloc(Nstc * sizeof(int));
-  count = (int *)omAlloc(Nstc * sizeof(int));
-  for (i = 1; i <= Nvar; i++)
-  {
-    i1 = var[i];
-    *temp = stc[0][i1];
-    *count = 1;
-    jj = 1;
-    for (j = 1; j < Nstc; j++)
-    {
-      x = stc[j][i1];
-      k = 0;
-      loop
-      {
-        if (x > temp[k])
-        {
-          k++;
-          if (k == jj)
-          {
-            temp[k] = x;
-            count[k] = 1;
-            jj++;
-            break;
-          }
-        }
-        else if (x < temp[k])
-        {
-          for (l = jj; l > k; l--)
-          {
-            temp[l] = temp[l-1];
-            count[l] = count[l-1];
-          }
-          temp[k] = x;
-          count[k] = 1;
-          jj++;
-          break;
-        }
-        else
-        {
-          count[k]++;
-          break;
-        }
-      }
-    }
-    h = 0.0;
-    o = (float)Nstc/(float)jj;
-    for(j = 0; j < jj; j++)
-    {
-       g = (float)count[j];
-       if (g > o)
-         g -= o;
-       else
-         g = o - g;
-       if (g > h)
-         h = g;
-    }
-    v1[i-1] = h * (float)jj;
-  }
-  omFreeSize((ADDRESS)count, Nstc * sizeof(int));
-  omFreeSize((ADDRESS)temp, Nstc * sizeof(int));
-  for (i = 1; i < Nvar; i++)
-  {
-    i1 = var[i+1];
-    h = v1[i];
-    j = 0;
-    loop
-    {
-      if (h > v1[j])
-      {
-        for (l = i; l > j; l--)
-        {
-          v1[l] = v1[l-1];
-          var[l+1] = var[l];
-        }
-        v1[j] = h;
-        var[j+1] = i1;
-        break;
-      }
-      j++;
-      if (j == i)
-        break;
-    }
-  }
-  omFreeSize((ADDRESS)v1, Nvar * sizeof(float));
-}
-
-
-static void hShrink(scfmon co, int a, int Nco)
-{
-  while ((co[a]!=NULL) && (a<Nco)) a++;
-  int i = a;
-  int j;
-  for (j = a; j < Nco; j++)
-  {
-    if (co[j]!=NULL)
-    {
-      co[i] = co[j];
-      i++;
-    }
-  }
-}
-
-
-void hStaircase(scfmon stc, int *Nstc, varset var, int Nvar)
-{
-  int  nc = *Nstc;
-  if (nc < 2)
-    return;
-  int z = 0;
-  int i = 0;
-  int j = 1;
-  scmon n = stc[1 /*j*/];
-  scmon o = stc[0];
-  int k = Nvar;
-  loop
-  {
-    int k1 = var[k];
-    if (o[k1] > n[k1])
-    {
-      loop
-      {
-        k--;
-        if (k==0)
-        {
-          stc[i] = NULL;
-          z++;
-          break;
-        }
-        else
-        {
-          k1 = var[k];
-          if (o[k1] < n[k1])
-            break;
-        }
-      }
-      k = Nvar;
-    }
-    else if (o[k1] < n[k1])
-    {
-      loop
-      {
-        k--;
-        if (k==0)
-        {
-          stc[j] = NULL;
-          z++;
-          break;
-        }
-        else
-        {
-          k1 = var[k];
-          if (o[k1] > n[k1])
-            break;
-        }
-      }
-      k = Nvar;
-    }
-    else
-    {
-      k--;
-      if (k==0)
-      {
-        stc[j] = NULL;
-        z++;
-        k = Nvar;
-      }
-    }
-    if (k == Nvar)
-    {
-      if (stc[j]==NULL)
-        i = j - 1;
-      loop
-      {
-        i++;
-        if (i == j)
-        {
-          i = -1;
-          j++;
-          if (j < nc)
-            n = stc[j];
-          else
-          {
-            if (z!=0)
-            {
-              *Nstc -= z;
-              hShrink(stc, 0, nc);
-            }
-            return;
-          }
-        }
-        else if (stc[i]!=NULL)
-        {
-          o = stc[i];
-          break;
-        }
-      }
-    }
-  }
-}
-
-
-void hRadical(scfmon rad, int *Nrad, int Nvar)
-{
-  int  nc = *Nrad, z = 0, i, j, k;
-  scmon n, o;
-  if (nc < 2)
-    return;
-  i = 0;
-  j = 1;
-  n = rad[j];
-  o = rad[0];
-  k = Nvar;
-  loop
-  {
-    if ((o[k]!=0) && (n[k]==0))
-    {
-      loop
-      {
-        k--;
-        if (k==0)
-        {
-          rad[i] = NULL;
-          z++;
-          break;
-        }
-        else
-        {
-          if ((o[k]==0) && (n[k]!=0))
-            break;
-        }
-      }
-      k = Nvar;
-    }
-    else if (!o[k] && n[k])
-    {
-      loop
-      {
-        k--;
-        if (!k)
-        {
-          rad[j] = NULL;
-          z++;
-          break;
-        }
-        else
-        {
-          if (o[k] && !n[k])
-            break;
-        }
-      }
-      k = Nvar;
-    }
-    else
-    {
-      k--;
-      if (!k)
-      {
-        rad[j] = NULL;
-        z++;
-        k = Nvar;
-      }
-    }
-    if (k == Nvar)
-    {
-      if (!rad[j])
-        i = j - 1;
-      loop
-      {
-        i++;
-        if (i == j)
-        {
-          i = -1;
-          j++;
-          if (j < nc)
-            n = rad[j];
-          else
-          {
-            if (z)
-            {
-              *Nrad -= z;
-              hShrink(rad, 0, nc);
-            }
-            return;
-          }
-        }
-        else if (rad[i])
-        {
-          o = rad[i];
-          break;
-        }
-      }
-    }
-  }
-}
-
-
-void hLexS(scfmon stc, int Nstc, varset var, int Nvar)
-{
-  if (Nstc < 2)
-    return;
-  int  j = 1, i = 0;
-  scmon n = stc[j];
-  scmon o = stc[0];
-  int k = Nvar;
-  loop
-  {
-    int k1 = var[k];
-    if (o[k1] < n[k1])
-    {
-      i++;
-      if (i < j)
-      {
-        o = stc[i];
-        k = Nvar;
-      }
-      else
-      {
-        j++;
-        if (j < Nstc)
-        {
-          i = 0;
-          o = stc[0];
-          n = stc[j];
-          k = Nvar;
-        }
-        else
-          return;
-      }
-    }
-    else if (o[k1] > n[k1])
-    {
-      int tmp_k;
-      for (tmp_k = j; tmp_k > i; tmp_k--)
-        stc[tmp_k] = stc[tmp_k - 1];
-      stc[i] = n;
-      j++;
-      if (j < Nstc)
-      {
-        i = 0;
-        o = stc[0];
-        n = stc[j];
-        k = Nvar;
-      }
-      else
-        return;
-    }
-    else
-    {
-      k--;
-      if (k<=0) return;
-    }
-  }
-}
-
-
-void hLexR(scfmon rad, int Nrad, varset var, int Nvar)
-{
-  int  j = 1, i = 0, k, k1;
-  scmon n, o;
-  if (Nrad < 2)
-    return;
-  n = rad[j];
-  o = rad[0];
-  k = Nvar;
-  loop
-  {
-    k1 = var[k];
-    if (!o[k1] && n[k1])
-    {
-      i++;
-      if (i < j)
-      {
-        o = rad[i];
-        k = Nvar;
-      }
-      else
-      {
-        j++;
-        if (j < Nrad)
-        {
-          i = 0;
-          o = rad[0];
-          n = rad[j];
-          k = Nvar;
-        }
-        else
-          return;
-      }
-    }
-    else if (o[k1] && !n[k1])
-    {
-      for (k = j; k > i; k--)
-        rad[k] = rad[k - 1];
-      rad[i] = n;
-      j++;
-      if (j < Nrad)
-      {
-        i = 0;
-        o = rad[0];
-        n = rad[j];
-        k = Nvar;
-      }
-      else
-        return;
-    }
-    else
-      k--;
-  }
-}
-
-
-void hPure(scfmon stc, int a, int *Nstc, varset var, int Nvar,
- scmon pure, int *Npure)
-{
-  int  nc = *Nstc, np = 0, nq = 0, j, i, i1, c, l;
-  scmon x;
-  for (j = a; j < nc; j++)
-  {
-    x = stc[j];
-    i = Nvar;
-    c = 2;
-    l = 0;
-    loop
-    {
-      i1 = var[i];
-      if (x[i1])
-      {
-        c--;
-        if (!c)
-        {
-          l = 0;
-          break;
-        }
-        else if (c == 1)
-          l = i1;
-      }
-      i--;
-      if (!i)
-        break;
-    }
-    if (l)
-    {
-      if (!pure[l])
-      {
-        np++;
-        pure[l] = x[l];
-      }
-      else if (x[l] < pure[l])
-        pure[l] = x[l];
-      stc[j] = NULL;
-      nq++;
-    }
-  }
-  *Npure = np;
-  if (nq!=0)
-  {
-    *Nstc -= nq;
-    hShrink(stc, a, nc);
-  }
-}
-
-
-void hElimS(scfmon stc, int *e1, int a2, int e2, varset var, int Nvar)
-{
-  int  nc = *e1, z = 0, i, j, k, k1;
-  scmon n, o;
-  if (!nc || (a2 == e2))
-    return;
-  j = 0;
-  i = a2;
-  o = stc[i];
-  n = stc[0];
-  k = Nvar;
-  loop
-  {
-    k1 = var[k];
-    if (o[k1] > n[k1])
-    {
-      k = Nvar;
-      i++;
-      if (i < e2)
-        o = stc[i];
-      else
-      {
-        j++;
-        if (j < nc)
-        {
-          i = a2;
-          o = stc[i];
-          n = stc[j];
-        }
-        else
-        {
-          if (z!=0)
-          {
-            *e1 -= z;
-            hShrink(stc, 0, nc);
-          }
-          return;
-        }
-      }
-    }
-    else
-    {
-      k--;
-      if (k==0)
-      {
-        stc[j] = NULL;
-        z++;
-        j++;
-        if (j < nc)
-        {
-          i = a2;
-          o = stc[i];
-          n = stc[j];
-          k = Nvar;
-        }
-        else
-        {
-          if (z!=0)
-          {
-            *e1 -= z;
-            hShrink(stc, 0, nc);
-          }
-          return;
-        }
-      }
-    }
-  }
-}
-
-
-void hElimR(scfmon rad, int *e1, int a2, int e2, varset var, int Nvar)
-{
-  int  nc = *e1, z = 0, i, j, k, k1;
-  scmon n, o;
-  if (!nc || (a2 == e2))
-    return;
-  j = 0;
-  i = a2;
-  o = rad[i];
-  n = rad[0];
-  k = Nvar;
-  loop
-  {
-    k1 = var[k];
-    if (o[k1] && !n[k1])
-    {
-      k = Nvar;
-      i++;
-      if (i < e2)
-        o = rad[i];
-      else
-      {
-        j++;
-        if (j < nc)
-        {
-          i = a2;
-          o = rad[i];
-          n = rad[j];
-        }
-        else
-        {
-          if (z!=0)
-          {
-            *e1 -= z;
-            hShrink(rad, 0, nc);
-          }
-          return;
-        }
-      }
-    }
-    else
-    {
-      k--;
-      if (!k)
-      {
-        rad[j] = NULL;
-        z++;
-        j++;
-        if (j < nc)
-        {
-          i = a2;
-          o = rad[i];
-          n = rad[j];
-          k = Nvar;
-        }
-        else
-        {
-          if (z!=0)
-          {
-            *e1 -= z;
-            hShrink(rad, 0, nc);
-          }
-          return;
-        }
-      }
-    }
-  }
-}
-
-
-void hLex2S(scfmon rad, int e1, int a2, int e2, varset var,
- int Nvar, scfmon w)
-{
-  int  j0 = 0, j = 0, i = a2, k, k1;
-  scmon n, o;
-  if (!e1)
-  {
-    for (; i < e2; i++)
-      rad[i - a2] = rad[i];
-    return;
-  }  else if (i == e2)
-    return;
-  n = rad[j];
-  o = rad[i];
-  loop
-  {
-    k = Nvar;
-    loop
-    {
-      k1 = var[k];
-      if (o[k1] < n[k1])
-      {
-        w[j0] = o;
-        j0++;
-        i++;
-        if (i < e2)
-        {
-          o = rad[i];
-          break;
-        }
-        else
-        {
-          for (; j < e1; j++)
-          {
-            w[j0] = rad[j];
-            j0++;
-          }
-          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
-          return;
-        }
-      }
-      else if (o[k1] > n[k1])
-      {
-        w[j0] = n;
-        j0++;
-        j++;
-        if (j < e1)
-        {
-          n = rad[j];
-          break;
-        }
-        else
-        {
-          for (; i < e2; i++)
-          {
-            w[j0] = rad[i];
-            j0++;
-          }
-          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
-          return;
-        }
-      }
-      k--;
-    }
-  }
-}
-
-
-void hLex2R(scfmon rad, int e1, int a2, int e2, varset var,
- int Nvar, scfmon w)
-{
-  int  j0 = 0, j = 0, i = a2, k, k1;
-  scmon n, o;
-  if (!e1)
-  {
-    for (; i < e2; i++)
-      rad[i - a2] = rad[i];
-    return;
-  }
-  else if (i == e2)
-    return;
-  n = rad[j];
-  o = rad[i];
-  loop
-  {
-    k = Nvar;
-    loop
-    {
-      k1 = var[k];
-      if (!o[k1] && n[k1])
-      {
-        w[j0] = o;
-        j0++;
-        i++;
-        if (i < e2)
-        {
-          o = rad[i];
-          break;
-        }
-        else
-        {
-          for (; j < e1; j++)
-          {
-            w[j0] = rad[j];
-            j0++;
-          }
-          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
-          return;
-        }
-      }
-      else if (o[k1] && !n[k1])
-      {
-        w[j0] = n;
-        j0++;
-        j++;
-        if (j < e1)
-        {
-          n = rad[j];
-          break;
-        }
-        else
-        {
-          for (; i < e2; i++)
-          {
-            w[j0] = rad[i];
-            j0++;
-          }
-          memcpy(rad, w, (e1 + e2 - a2) * sizeof(scmon));
-          return;
-        }
-      }
-      k--;
-    }
-  }
-}
-
-
-void hStepS(scfmon stc, int Nstc, varset var, int Nvar, int *a, int *x)
-{
-  int  k1, i;
-  int  y;
-  k1 = var[Nvar];
-  y = *x;
-  i = *a;
-  loop
-  {
-    if (y < stc[i][k1])
-    {
-      *a = i;
-      *x = stc[i][k1];
-      return;
-    }
-    i++;
-    if (i == Nstc)
-    {
-      *a = i;
-      return;
-    }
-  }
-}
-
-
-void hStepR(scfmon rad, int Nrad, varset var, int Nvar, int *a)
-{
-  int  k1, i;
-  k1 = var[Nvar];
-  i = 0;
-  loop
-  {
-    if (rad[i][k1])
-    {
-      *a = i;
-      return;
-    }
-    i++;
-    if (i == Nrad)
-    {
-      *a = i;
-      return;
-    }
-  }
-}
-
-
-monf hCreate(int Nvar)
-{
-  monf xmem;
-  int  i;
-  xmem = (monf)omAlloc((Nvar + 1) * sizeof(monp));
-  for (i = Nvar; i>0; i--)
-  {
-    xmem[i] = (monp)omAlloc(LEN_MON);
-    xmem[i]->mo = NULL;
-  }
-  return xmem;
-}
-
-
-void hKill(monf xmem, int Nvar)
-{
-  int  i;
-  for (i = Nvar; i!=0; i--)
-  {
-    if (xmem[i]->mo!=NULL)
-      omFreeSize((ADDRESS)xmem[i]->mo, xmem[i]->a * sizeof(scmon));
-    omFreeSize((ADDRESS)xmem[i], LEN_MON);
-  }
-  omFreeSize((ADDRESS)xmem, (Nvar + 1) * sizeof(monp));
-}
-
-
-scfmon hGetmem(int lm, scfmon old, monp monmem)
-{
-  scfmon x = monmem->mo;
-  int  lx = monmem->a;
-  if ((x==NULL) || (lm > lx))
-  {
-    /* according to http://www.singular.uni-kl.de:8002/trac/ticket/463#comment:4
-     * we need to work around a compiler bug:
-    * if ((x!=NULL)&&(lx>0)) omFreeSize((ADDRESS)x, lx * sizeof(scmon));
-    */
-    if (x!=NULL) if (lx>0) omFreeSize((ADDRESS)x, lx * sizeof(scmon));
-    monmem->mo = x = (scfmon)omAlloc(lm * sizeof(scmon));
-    monmem->a = lm;
-  }
-  memcpy(x, old, lm * sizeof(scmon));
-  return x;
-}
-
-/*
-* a bug in Metrowerks with "lifetime analysis"
-*scmon hGetpure(scmon p)
-*{
-*  scmon p1, pn;
-*  p1 = p + 1;
-*  pn = p1 + (currRing->N);
-*  memcpy(pn, p1, (currRing->N) * sizeof(int));
-*  return pn - 1;
-*}
-*/
-scmon hGetpure(scmon p)
-{
-  scmon p1 = p;
-  scmon pn;
-  p1++;
-  pn = p1;
-  pn += (currRing->N);
-  memcpy(pn, p1, (currRing->N) * sizeof(int));
-  return pn - 1;
-}
-
