Changeset 50cbdc in git for Singular/walk.cc

Timestamp:

Aug 27, 2001, 4:48:02 PM (23 years ago)

Author:

Hans Schönemann <hannes@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

2567b5a6cb7109be5a83e53eb94abb1c38fb9945

Parents:

3de58c9ca0aeaafdf5cb29f986967bffa405b542

Message:

*hannes: merge-2-0-2


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

File:

• Singular/walk.cc (modified) (5 diffs)

Singular/walk.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR      *
 *****************************************/
-/* $Id: walk.cc,v 1.4 2000-09-07 13:39:45 sulandra Exp $ */
+/* $Id: walk.cc,v 1.5 2001-08-27 14:47:45 Singular Exp $ */
 /*
 * ABSTRACT: Implementation of the Groebner walk
 */
 
+/* includes */
 #include "mod2.h"
 #include "walk.h"
@@ -13 +14 @@
 #include "intvec.h"
 #include "ipid.h"
-
-// add two intvecs:
-intvec* walkAddIntVec(intvec* v1, intvec* v2)
-{
-  int n = v1->length();
+#include "tok.h"
+#include <omalloc.h>
+#include "febase.h"
+#include "numbers.h"
+#include "ipid.h"
+#include "ring.h"
+#include "kstd1.h"
+#include "matpol.h"
+#include "weight.h"
+#include "intvec.h"
+#include "syz.h"
+#include "lists.h"
+#include "prCopy.h"
+#include <string.h>
+#include "structs.h"
+#include "longalg.h"
+#ifdef HAVE_FACTORY
+#include "clapsing.h"
+#endif
+
+static void* idString(ideal L)
+{
   int i;
-  intvec *result = new intvec(n);
-  if (v2->length() > n) n = v2->length();
-  
-  for (i=0; i<n; i++)
-  {
-    (*result)[i] = (*v1)[i] + (*v2)[i];
-  }
-  
-  return result;
-}
-
-  
-  
-
-// scalar product of weights and exponent vector of p
-// assumes that weights and exponent vector have length n
-inline long walkWeightDegree(const poly p, const int* weights,
-                             const long n)
-{
-  assume(p != NULL && weights != NULL);
-
-  long i, res = 0;
-
-  for (i=0; i<n; i++) res += pGetExp(p, i+1) * weights[i];
-
-  return res;
+  printf("//ideal Itmp: ");
+  for(i=0; i<IDELEMS(L); i++)
+    printf(" %s, ", pString(L->m[i]));
+  printf("\n");
 }
 
 
 // returns gcd of integers a and b
-inline long gcd(const long a, const long b)
+static inline long gcd(const long a, const long b)
 {
   long r, p0 = a, p1 = b;
-  assume(p0 >= 0 && p1 >= 0);
-
+  //assume(p0 >= 0 && p1 >= 0);
+  if(p0 < 0)
+    p0 = -p0;
+
+  if(p1 < 0)
+    p1 = -p1;
   while(p1 != 0)
   {
@@ -64 +64 @@
 
 // cancel gcd of integers zaehler and nenner
-inline void cancel(long &zaehler, long &nenner)
+static inline void cancel(long &zaehler, long &nenner)
 {
   assume(zaehler >= 0 && nenner > 0);
@@ -75 +75 @@
 }
 
-// Returns the next Weight vector for the Groebner walk
-// Assumes monoms of polys of G are ordered decreasingly w.r.t. curr_weight
-int* walkNextWeight(const int* curr_weight,
-                    const int* target_weight,
-                    const ideal G)
-{
-  assume(currRing != NULL && target_weight != NULL && curr_weight != NULL &&
-         G != NULL);
-
-  int* diff_weight =
-    (int*)omAlloc(currRing->N*sizeof(int));
-  long j, t_zaehler = 0, t_nenner = 0;
-
-  for (j=0; j<currRing->N; j++)
-    diff_weight[j] = target_weight[j] - curr_weight[j];
-
-  for (j=0; j<IDELEMS(G); j++)
-  {
-    poly g = G->m[j];
-    if (g != 0)
-    {
-      long deg_w0_p1 = pGetOrder(g);
-      long deg_d0_p1 = walkWeightDegree(g, diff_weight, currRing->N);
-
+/********************************************************************
+ * compute a weight degree of a monomial p w.r.t. a weight_vector   *
+ ********************************************************************/
+static inline int MLmWeightedDegree(const poly p, intvec* weight)
+{
+  int i, d = 0;
+  
+  for (i=1; i<=pVariables; i++)
+    d += pGetExp(p, i) * (*weight)[i-1];
+
+  return d;
+}
+
+/********************************************************************
+ * compute a weight degree of a polynomial p w.r.t. a weight_vector *
+ ********************************************************************/
+static inline int MwalkWeightDegree(poly p, intvec* weight_vector)
+{
+  assume(weight_vector->length() >= pVariables);
+  int max = 0, maxtemp;
+  poly hp;
+
+  while(p != NULL)
+  {
+    hp = pHead(p); 
+    pIter(p);
+    maxtemp = MLmWeightedDegree(hp, weight_vector);
+
+    if (maxtemp > max) 
+      max = maxtemp;
+  }
+  return max;
+}
+
+/*****************************************************************************
+ * return an initial form of the polynom g w.r.t. a weight vector curr_weight *
+ *****************************************************************************/
+static poly MpolyInitialForm(poly g, intvec* curr_weight)
+{
+  if(g == NULL)
+    return g;
+  
+  int maxtmp, max = 0;
+  poly in_w_g = NULL, hg;    
+ 
+  while(g != NULL)
+  {
+    hg = pHead(g);
+    pIter(g);
+    maxtmp = MwalkWeightDegree(hg, curr_weight);
+
+    if(maxtmp > max)  
+    {
+      max = maxtmp;
+      in_w_g = hg;
+    } else {
+      if(maxtmp == max)
+        in_w_g = pAdd(in_w_g, hg);
+    }
+  }
+  return in_w_g;
+}
+
+
+/*****************************************************************************
+ * compute the initial form of an ideal "G" w.r.t. weight vector curr_weight *
+ ****************************************************************************/
+ideal MwalkInitialForm(ideal G, intvec* curr_weight)
+{
+  int i;
+  int nG = IDELEMS(G);
+  ideal Gomega  = idInit(nG, G->rank); 
+
+  for(i=0; i<nG; i++)
+    Gomega->m[i] = MpolyInitialForm(G->m[i], curr_weight);
+
+  //return Gomega;
+  ideal result = idCopy(Gomega);
+  idDelete(&Gomega);
+  return result;
+}
+
+/************************************************************************
+ * test that does the weight vector iv exist in the cone of the ideal G *
+ *            i.e. does in(in_w(g)) =? in(g), for all g in G            *
+ ************************************************************************/
+void* test_w_in_Cone(ideal G, intvec* iv)
+{
+  int nG = IDELEMS(G);
+  int i;
+  BOOLEAN ok = TRUE;
+  poly mi, in_mi,  gi;
+  for(i=0; i<nG; i++)
+  {
+    mi = MpolyInitialForm(G->m[i], iv);
+    in_mi = pHead(mi);
+    gi = pHead(G->m[i]);
+    if(pEqualPolys(in_mi, gi) != ok)
+    {
+      printf("//ring Test_W_in_Cone = %s ;\n", rString(currRing));
+      printf("//the computed next weight vector doesn't exist in the cone\n");
+      break;
+    }
+  }
+}
+
+//compute a least common multiple of two integers
+static inline long Mlcm(long &i1, long &i2)
+{
+  long temp = gcd(i1, i2);  
+  return ((i1*i2) / temp);
+}
+
+
+/***************************************************
+ * return  the dot product of two intvecs a and b  *
+ ***************************************************/
+static inline long  MivDotProduct(intvec* a, intvec* b)
+{
+  assume( a->length() ==  b->length());
+  int i, n = a->length();
+  long result = 0;
+  
+  for(i=0; i<n; i++)
+    result += (*a)[i] * (*b)[i];
+
+  return result;
+}
+
+
+/**21.10.00*******************************************
+ * return the "intvec" lead exponent of a polynomial *
+ *****************************************************/
+static intvec* MExpPol(poly f)
+{
+  int nR = currRing->N;
+
+  intvec* result = new intvec(nR);
+  int i;
+  
+  for(i=0; i<nR; i++)
+    (*result)[i] = pGetExp(f,i+1);    
+
+  intvec *res = ivCopy(result);
+  omFree((ADDRESS) result);
+  return res;
+}
+
+
+/***23-24.10.00******************************************
+ * compute a division of two monoms, "a" by a monom "b" *
+ *    i.e. leading term of two polynoms a and b         *
+ ********************************************************/
+static poly MpDiv(poly a, poly b)
+{
+  assume (b != NULL);
+  BOOLEAN ok = TRUE;
+  
+  if(a == NULL)
+    return a;  
+  
+  int nR = currRing->N;
+  
+  number nn = (number) omAllocBin(rnumber_bin);
+      
+  poly  ptmp, ppotenz; 
+  poly result = pISet(1); 
+
+  intvec* iva = MExpPol(a);  //head exponent of a   
+  intvec* ivb = MExpPol(b);  //head exponent of a  
+
+  int nab;
+  for(int i=0; i<nR; i++)
+  {
+    nab = (*iva)[i] - (*ivb)[i];
+    // b does not divide a
+    if(nab < 0)
+    {
+      result = NULL;
+      return result;
+    }
+    //define a polynomial which is a variable of the basering
+    ptmp = (poly) pmInit(currRing->names[i], ok);  //p:=xi
+    ppotenz = pPower(ptmp, nab); 
+    result = pMult(result, ppotenz);
+  }
+  nn = nDiv(pGetCoeff(a), pGetCoeff(b));
+  result = pMult_nn(result, nn);
+  nDelete(&nn);
+
+  return result;
+}
+
+
+/***24.10.00 *****************************************
+ *      compute a product of two monoms a and b      *
+ *      i.e. leading term of two polynoms a and b    *
+ *****************************************************/
+static poly MpMult(poly a, poly b)
+{
+  if(a == NULL || b == NULL)
+    return a;  
+  
+  int nR = currRing->N;
+  BOOLEAN ok = TRUE;  
+  
+  poly  ptmp, ppotenz; 
+  poly result = pISet(1);  // result := 1  
+  intvec* ivab = ivAdd(MExpPol(a), MExpPol(b));
+
+  for(int i=0; i<nR; i++)
+  {
+    //define a polynomial which is a variable of the basering
+    ptmp = pmInit(currRing->names[i], ok);
+    ppotenz = pPower(ptmp, (*ivab)[i]); 
+    result = pMult(result, ppotenz);
+  }
+  number nn = nMult(pGetCoeff(a), pGetCoeff(b)); 
+  result = pMult_nn(result, nn);
+
+  return result;
+}
+
+
+poly  MivSame(intvec* u , intvec* v)
+{ 
+  assume(u->length() == v->length());
+
+  int i, niv = u->length();
+  
+  for (i=0; i<niv; i++)
+    if ((*u)[i] != (*v)[i])
+      return pISet(1);
+
+  return (poly) NULL; 
+}
+
+poly  M3ivSame(intvec* temp, intvec* u , intvec* v)
+{ 
+  assume(temp->length() == u->length() && u->length() == v->length());
+
+  if(MivSame(temp, u) == NULL)
+    return (poly) NULL;
+  if(MivSame(temp, v) == NULL)
+    return pISet(1);
+  return pISet(2);  
+}
+
+
+/************************
+ *  Define a monom x^iv *
+ ************************/
+poly MPolVar(intvec* iv)
+{
+  int i, niv = pVariables;
+
+  poly ptemp = pOne();
+  poly pvar, ppotenz;
+  BOOLEAN ok = TRUE;
+  
+  for(i=0; i<niv; i++)
+  {
+    pvar = (poly) pmInit(currRing->names[i], ok);  //p:=x_i
+    ppotenz = pPower(pvar, (*iv)[i]); 
+    ptemp = pMult(ptemp, ppotenz);
+  }
+  return ptemp;
+}
+
+
+/* compute a Groebner basis of an ideal */
+ideal Mstd(ideal G)
+{
+  G = kStd(G, NULL, testHomog, NULL);
+  G = kInterRed(G, NULL);//5.02
+  idSkipZeroes(G);
+  return G;
+}
+
+/* compute a Groebner basis of a homogenoues ideal */
+ideal Mstdhom(ideal G)
+{
+  G = kStd(G, NULL, isHomog, NULL);
+  G = kInterRed(G, NULL);//21.02
+  idSkipZeroes(G);  
+  return G;
+}
+
+/* compute a reduced Groebner basis of a Groebner basis */
+ideal MkInterRed(ideal G)
+{  
+  if(G == NULL)
+    return G;
+  
+  G = kInterRed(G, NULL);
+  idSkipZeroes(G);
+  return G;
+}
+
+
+/*****************************************************
+ *              PERTURBATION WALK                    *
+ *****************************************************/
+
+/***************************************
+* create an ordering matrix as intvec  *
+****************************************/
+intvec* MivMatrixOrder(intvec* iv)
+{
+  int i,j;
+  int nR = currRing->N; 
+  intvec* ivm = new intvec(nR*nR);
+
+  for(i=0; i<nR; i++)
+    (*ivm)[i] = (*iv)[i];
+
+  for(i=1; i<nR; i++)
+    (*ivm)[i*nR+i-1] = (int) 1;
+
+  return ivm;
+}
+
+static intvec* MivMatUnit(void)
+{
+  int nR = currRing->N; 
+  intvec* ivm = new intvec(nR);
+
+  (*ivm)[0] = 1;
+
+  return ivm;
+}
+
+/* return iv = (1, ..., 1) */
+intvec* Mivdp(int nR)
+{
+  int i; 
+  intvec* ivm = new intvec(nR);
+
+  for(i=0; i<nR; i++)
+    (*ivm)[i] = 1;
+
+  return ivm;
+}
+
+/* return iv = (1,0, ..., 0) */
+intvec* Mivlp(int nR)
+{
+  int i; 
+  intvec* ivm = new intvec(nR);
+  (*ivm)[0] = 1;
+
+  return ivm;  
+}
+
+intvec* Mivdp0(int nR)
+{
+  int i; 
+  intvec* ivm = new intvec(nR);
+  (*ivm)[nR-1] = 0;
+  for(i=0; i<nR-1; i++)
+    (*ivm)[i] = 1;
+
+  return ivm;
+}
+
+/*****************************************************************************
+* If target_ord = intmat(A1, ..., An) then calculate the perturbation        *
+* vectors                                                                    *
+*   tau_p_dep = inveps^(p_deg-1)*A1 + inveps^(p_deg-2)*A2 +... + A_p_deg     *
+* where                                                                      *
+*      inveps > totaldegree(G)*(max(A2)+...+max(A_p_deg))                    *
+* intmat target_ord is an integer order matrix of the monomial ordering of   *
+* basering.                                                                  *
+* This programm computes a perturbated vector with a p_deg perturbation      *
+* degree which smaller than the numbers of varibles                          * 
+******************************************************************************/
+/* ivtarget is a matrix of a degree reverse lex. order */
+intvec* MPertVectors(ideal G, intvec* ivtarget, int pdeg)
+{
+  int nV = currRing->N;
+  //assume(pdeg <= nV && pdeg >= 0);
+
+  int i, j;
+  intvec* pert_vector =  new intvec(nV);
+  
+  //Checking that the perturbated degree is valid
+  if(pdeg > nV || pdeg <= 0)
+  {  
+    WerrorS("The perturbed degree is wrong!!");
+    return pert_vector;
+  }
+  for(i=0; i<nV; i++)
+    (*pert_vector)[i]=(*ivtarget)[i];
+
+  if(pdeg == 1)
+    return pert_vector;
+     
+  // Calculate max1 = Max(A2)+Max(A3)+...+Max(Apdeg),
+  // where the Ai are the i-te rows of the matrix target_ord.
+
+  int ntemp, maxAi, maxA=0;
+  //for(i=1; i<pdeg; i++)
+  for(i=0; i<pdeg; i++) //for "dp"
+  {
+    maxAi = (*ivtarget)[i*nV]; 
+    for(j=i*nV+1; j<(i+1)*nV; j++)
+    {
+      ntemp = (*ivtarget)[j];
+      if(ntemp > maxAi)
+        maxAi = ntemp;
+    }
+    maxA += maxAi;    
+  }
+  
+  // Calculate inveps = 1/eps, where 1/eps > totaldeg(p)*max1 for all p in G.
+  int inveps, tot_deg = 0, maxdeg;
+
+  intvec* ivUnit = Mivdp(nV);//19.02 
+  for(i=0; i<IDELEMS(G); i++)
+  {
+    //maxdeg = pTotaldegree(G->m[i], currRing); //it's wrong for ex1,2,rose
+    maxdeg = MwalkWeightDegree(G->m[i], ivUnit);
+    if (maxdeg > tot_deg ) 
+      tot_deg = maxdeg;
+  }
+  inveps = (tot_deg * maxA) + 1;
+  
+  // pert(A1) = inveps^(pdeg-1)*A1 + inveps^(pdeg-2)*A2+...+A_pdeg,
+  // pert_vector := A1
+  for ( i=1; i < pdeg; i++ )
+    for(j=0; j<nV; j++)
+       (*pert_vector)[j] = inveps*(*pert_vector)[j] + (*ivtarget)[i*nV+j];
+
+  
+  int temp = (*pert_vector)[0];
+  for(i=1; i<nV; i++)
+  {
+    temp = gcd(temp, (*pert_vector)[i]);
+    if(temp == 1)
+      break;
+  }
+  if(temp != 1)
+    for(i=0; i<nV; i++)
+      (*pert_vector)[i] = (*pert_vector)[i] / temp;
+  
+  //test_w_in_Cone(G, pert_vector);
+  return pert_vector;
+}
+
+/* ivtarget is a matrix of the lex. order */
+intvec* MPertVectorslp(ideal G, intvec* ivtarget, int pdeg)
+{
+  int nV = currRing->N;
+  //assume(pdeg <= nV && pdeg >= 0);
+
+  int i, j;
+  intvec* pert_vector =  new intvec(nV);
+  
+  //Checking that the perturbated degree is valid
+  if(pdeg > nV || pdeg <= 0)
+  {  
+    WerrorS("The perturbed degree is wrong!!");
+    return pert_vector;
+  }
+  for(i=0; i<nV; i++)
+    (*pert_vector)[i]=(*ivtarget)[i];
+
+  if(pdeg == 1)
+    return pert_vector;
+     
+  // Calculate max1 = Max(A2)+Max(A3)+...+Max(Apdeg),
+  // where the Ai are the i-te rows of the matrix target_ord.
+  int ntemp, maxAi, maxA=0;
+  for(i=1; i<pdeg; i++)
+  //for(i=0; i<pdeg; i++) //for "dp"
+  {
+    maxAi = (*ivtarget)[i*nV]; 
+    for(j=i*nV+1; j<(i+1)*nV; j++)
+    {
+      ntemp = (*ivtarget)[j];
+      if(ntemp > maxAi)
+        maxAi = ntemp;
+    }
+    maxA += maxAi;    
+  }
+  
+  // Calculate inveps := 1/eps, where 1/eps > deg(p)*max1 for all p in G.
+  int inveps, tot_deg = 0, maxdeg;
+
+  intvec* ivUnit = Mivdp(nV);//19.02 
+  for(i=0; i<IDELEMS(G); i++)
+  {
+    //maxdeg = pTotaldegree(G->m[i], currRing); //it's wrong for ex1,2,rose
+    maxdeg = MwalkWeightDegree(G->m[i], ivUnit);
+    if (maxdeg > tot_deg ) 
+      tot_deg = maxdeg;
+  }
+  inveps = (tot_deg * maxA) + 1;
+
+  // Pert(A1) = inveps^(pdeg-1)*A1 + inveps^(pdeg-2)*A2+...+A_pdeg,
+
+  for ( i=1; i < pdeg; i++ )
+    for(j=0; j<nV; j++)
+      (*pert_vector)[j] = inveps*((*pert_vector)[j]) + (*ivtarget)[i*nV+j];
+    
+  int temp = (*pert_vector)[0];
+  for(i=1; i<nV; i++)
+  {
+    temp = gcd(temp, (*pert_vector)[i]);
+    if(temp == 1)
+      break;
+  }
+  if(temp != 1)
+    for(i=0; i<nV; i++)
+      (*pert_vector)[i] = (*pert_vector)[i] / temp;
+
+  //test_w_in_Cone(G, pert_vector);
+  return pert_vector;
+}
+
+
+/***************************************************************
+ *                    FRACTAL WALK                             *
+ ***************************************************************/
+
+/*****  define a lexicographic order matrix as intvec  ********/
+intvec* MivMatrixOrderlp(int nV)
+{
+  int i;
+  intvec* ivM = new intvec(nV*nV);
+      
+  for(i=0; i<nV; i++)
+    (*ivM)[i*nV + i] = 1;
+  
+  return(ivM);
+}
+
+intvec* MivMatrixOrderdp(int nV)
+{
+  int i;
+  intvec* ivM = new intvec(nV*nV);
+      
+  for(i=0; i<nV; i++)
+    (*ivM)[i] = 1;
+
+  for(i=1; i<nV; i++)
+    (*ivM)[(i+1)*nV - i] = -1;
+  
+  return(ivM);
+}
+
+//creates an intvec of the monomial order Wp(ivstart)
+intvec* MivWeightOrderlp(intvec* ivstart)
+{
+  int i;
+  int nV = ivstart->length();
+  intvec* ivM = new intvec(nV*nV);
+      
+  for(i=0; i<nV; i++)
+    (*ivM)[i] = (*ivstart)[i];
+
+  for(i=1; i<nV; i++)
+    (*ivM)[i*nV + i-1] = 1;
+  
+  return(ivM);
+}
+
+intvec* MivWeightOrderdp(intvec* ivstart)
+{
+  int i;
+  int nV = ivstart->length();
+  intvec* ivM = new intvec(nV*nV);
+      
+  for(i=0; i<nV; i++)
+    (*ivM)[i] = (*ivstart)[i];
+
+  for(i=0; i<nV; i++)
+    (*ivM)[nV+i] = 1;
+
+  for(i=2; i<nV; i++)
+    (*ivM)[(i+1)*nV - i] = -1;
+  
+  return(ivM);
+}
+
+intvec* MivUnit(int nV)
+{
+  int i;
+  intvec* ivM = new intvec(nV);
+      
+  for(i=0; i<nV; i++)
+    (*ivM)[i] = 1;
+  
+  return(ivM);
+}
+
+/************************************************************************
+*  compute a perturbed weight vector of a matrix order w.r.t. an ideal  *
+*************************************************************************/
+intvec* Mfpertvector(ideal G, intvec* ivtarget)
+//intvec* Mfpertvector(ideal G) 
+{
+  int i, j;
+  int nV = currRing->N;
+  int niv = nV*nV;
+
+  // Calculate max1 = Max(A2) + Max(A3) + ... + Max(AnV),
+  // where the Ai are the i-te rows of the matrix 'targer_ord'.
+  int ntemp, maxAi, maxA=0;
+  for(i=1; i<nV; i++) //30.03   
+    //for(i=0; i<nV; i++) //for "dp"
+  {
+    maxAi = (*ivtarget)[i*nV]; 
+    for(j=i*nV+1; j<(i+1)*nV; j++)
+    {
+      ntemp = (*ivtarget)[j];
+      if(ntemp > maxAi)
+        maxAi = ntemp;
+    }
+    maxA = maxA + maxAi;    
+  }
+  intvec* ivUnit = Mivdp(nV);
+  
+  // Calculate inveps = 1/eps, where 1/eps > deg(p)*max1 for all p in G.
+  int inveps, tot_deg = 0, maxdeg;
+  for(i=0; i<IDELEMS(G); i++)
+  {
+    maxdeg = MwalkWeightDegree(G->m[i], ivUnit);
+    //maxdeg = pTotaldegree(G->m[i]);
+    if (maxdeg > tot_deg ) 
+      tot_deg = maxdeg;
+  }
+  inveps = (tot_deg * maxA) + 1;  
+
+  // Calculate the perturbed target orders:
+  intvec* ivtemp = new intvec(nV);  
+  intvec* pert_vector = new intvec(niv);
+  for(i=0; i<nV; i++)
+  {
+    ntemp = (*ivtarget)[i];
+    (* ivtemp)[i] = ntemp;
+    (* pert_vector)[i] = ntemp;
+  }
+  for(i=1; i<nV; i++)
+  {
+    for(j=0; j<nV; j++)
+      (* ivtemp)[j] = inveps*(*ivtemp)[j] + (*ivtarget)[i*nV+j];
+    for(j=0; j<nV; j++) 
+      (* pert_vector)[i*nV+j] = (* ivtemp)[j];   
+  }
+  omFree((ADDRESS)ivtemp);
+  return pert_vector;
+}
+
+
+/**********************************************************************
+ *  computes a transformation matrix as an ideal L
+    an i-th element of L is a representasion of an i-th element M w.r.t.
+    the generators of Gomega  
+********************************************************************/
+
+ideal MidLift(ideal Gomega, ideal M)
+{
+  //M = idLift(Gomega, M, NULL, FALSE, FALSE, TRUE, NULL); 
+  //return M;
+  //17.04.01
+  ideal Mtmp = idInit(IDELEMS(M),1); 
+  Mtmp = idLift(Gomega, M, NULL, FALSE, FALSE, TRUE, NULL);  
+  idSkipZeroes(Mtmp);
+  M = idCopy(Mtmp); 
+
+  omFree((ADDRESS)Mtmp); 
+  return M; 
+}
+
+/****************************************************************
+ * Multiplikation of two ideals by elementwise                  *
+ * i.e. Let be A := (a_i) and B := (b_i), return C := (a_i*b_i) *
+ ****************************************************************/
+ideal MidMultLift(ideal A, ideal B)
+{
+  int mA = IDELEMS(A), mB = IDELEMS(B);
+  ideal result;
+
+  if(A==NULL || B==NULL)
+    return result;
+  
+  if(mB < mA)
+  {    
+    mA = mB;
+    result = idInit(mA, 1);
+  }
+  else
+  result = idInit(mA, 1);
+  
+  int i, k=0;
+  for(i=0; i<mA; i++)
+    if(A->m[i] != NULL)
+    {  
+      result->m[k] = pMult(pCopy(A->m[i]), pCopy(B->m[i]));
+      k++;
+    }      
+
+  //idSkipZeroes(result); //walkalp_CON
+  ideal res = idCopy(result);
+  idDelete(&result);
+  return res;  
+}
+
+//computes a  multiplication of two ideals L and G, ie. L[i]*G[i] 
+ideal MLiftLmalG(ideal L, ideal G)
+{
+  int i, j;
+  ideal Gtemp = idInit(IDELEMS(L),1);
+  ideal mG = idInit(IDELEMS(G),1);
+  poly pGtmp = NULL, pmG;
+  ideal T;
+
+  for(i=0; i<IDELEMS(L); i++)
+  {
+    T = idVec2Ideal(L->m[i]);
+    mG = MidMultLift(T,G);
+    idSkipZeroes(mG);
+ 
+    for(j=0; j<IDELEMS(mG); j++)
+    {
+      pGtmp = pAdd(pGtmp, mG->m[j]);
+    }
+    Gtemp->m[i] = pCopy(pGtmp);
+  }
+  idSkipZeroes(Gtemp);
+
+  //compute a reduced Groebner basis of GF
+  //Gtemp = kInterRed(Gtemp, NULL); 
+  L = idCopy(Gtemp);
+
+  omFree((ADDRESS)mG);
+  omFree((ADDRESS)Gtemp);  
+  return L; 
+}
+
+/*********************************************************************
+ * G is a red. Groebner basis w.r.t. <_1                             *
+ * Gomega is an initial form ideal of <G> w.r.t. a weight vector w   *
+ * M is a subideal of <Gomega> and M selft is a red. Groebner basis  *
+ *    of the ideal <Gomega> w.r.t. <_w                               *
+ * Let m_i = h1.gw1 + ... + hs.gws for each m_i in M; gwi in Gomega  *
+ * return F with n(F) = n(M) and f_i = h1.g1 + ... + hs.gs for each i*
+ ********************************************************************/
+ /* MidLift + MLiftLmalG */
+ideal MLiftLmalGNew(ideal Gomega, ideal M, ideal G)
+{
+  int i, j;
+  M = idLift(Gomega, M, NULL, FALSE, FALSE, TRUE, NULL); 
+  int nM = IDELEMS(M);
+  ideal Gtemp = idInit(IDELEMS(M),1);
+  ideal mG = idInit(IDELEMS(G),1);
+  poly pmG, pGtmp = NULL;
+  ideal T;
+  
+  for(i=0; i<nM; i++)
+  {
+    T = idVec2Ideal(M->m[i]);
+    mG = MidMultLift(T, G);
+     
+    for(j=0; j<IDELEMS(mG); j++)
+      pGtmp = pAdd(pGtmp, mG->m[j]);
+    
+    Gtemp->m[i] = pCopy(pGtmp);
+  }
+  idSkipZeroes(Gtemp);
+
+  M = idCopy(Gtemp);
+
+  omFree((ADDRESS)mG); 
+  omFree((ADDRESS)Gtemp);
+  return M;
+}
+
+/******************************************************************************
+ * compute a next weight vector on the line from curr_weight to target_weight *
+ * and it still stays in the cone of the ideal G where the curr_weight too    *
+ *****************************************************************************/
+intvec* MwalkNextWeight(intvec* curr_weight, intvec* target_weight, ideal G)
+{
+  assume(currRing != NULL && curr_weight != NULL && 
+         target_weight != NULL && G != NULL);
+
+  int nRing = currRing->N;
+  int nG = IDELEMS(G);
+  intvec* ivtemp;
+
+  long t_zaehler = 0, t_nenner = 0;
+  long s_zaehler, s_nenner, temp, MwWd;
+  long deg_w0_p1, deg_d0_p1;
+  int j;
+  
+  intvec* diff_weight = ivSub(target_weight, curr_weight);
+  poly g;
+  for (j=0; j<nG; j++)
+  { 
+    g = G->m[j];
+    if (g != NULL)  
+    {
+      ivtemp = MExpPol(g);
+      deg_w0_p1 = MivDotProduct(ivtemp, curr_weight);
+      deg_d0_p1 = MivDotProduct(ivtemp, diff_weight);
+      
       pIter(g);
 
       while (g != NULL)
       {
-        // compute s = s_zahler / s_nenner
-        long s_zaehler = deg_w0_p1 - pGetOrder(g);
+        //s_zaehler = deg_w0_p1 - pGetOrder(g);
+        ivtemp = MExpPol(g);
+        MwWd = MivDotProduct(ivtemp, curr_weight);
+        s_zaehler = deg_w0_p1 - MwWd;
 
         if (s_zaehler != 0)
-        {
-          long s_nenner =
-            walkWeightDegree(g, diff_weight, currRing->N) - deg_d0_p1;
+        {          
+          //s_nenner = MwalkWeightDegree(g, diff_weight) - deg_d0_p1;
+          MwWd = MivDotProduct(ivtemp, diff_weight);
+          s_nenner = MwWd - deg_d0_p1;
+
           // check for 0 < s <= 1
           if ( (s_zaehler > 0 && s_nenner >= s_zaehler) ||
@@ -121 +909 @@
             }
 
-            // look whether s < t
-            if (t_nenner == 0 ||
-                s_zaehler*t_nenner < t_zaehler * s_nenner)
+            // compute a simply fraction of s
+            cancel(s_zaehler, s_nenner);
+
+            if(t_nenner != 0)
             {
-              cancel(s_zaehler, s_nenner);
-              t_zaehler = s_zaehler;
+              if(s_zaehler * t_nenner < s_nenner * t_zaehler)
+              {
+                t_nenner = s_nenner;
+                t_zaehler = s_zaehler;
+              }       
+            }
+            else
+            {
               t_nenner = s_nenner;
+              t_zaehler = s_zaehler; 
             }
           }
         }
-        pIter(g);
+        pIter(g); //g = g - pHead(g);
       }
     }
   }
-
-  // return if no t or if t == 1
-  if (t_nenner == 0 || t_nenner == 1)
-  {
-    omFreeSize(diff_weight, currRing->N*sizeof(int));
-    return  (int*) t_nenner;
-  }
-
-  // construct new weight vector
-  for (j=0; j<currRing->N; j++)
-    diff_weight[j] = t_nenner*curr_weight[j] + t_zaehler*diff_weight[j];
+  if(t_nenner == 0)
+  {
+    diff_weight = ivCopy(curr_weight); 
+    return diff_weight; 
+  }
+  
+  if(t_nenner == 1 && t_zaehler == 1)
+  {
+    diff_weight = ivCopy(target_weight);
+    return diff_weight;
+  }
+    
+  // construct a new weight vector
+  for (j=0; j<nRing; j++)
+  {    
+    (*diff_weight)[j] = t_nenner*(*curr_weight)[j] +
+      t_zaehler*(*diff_weight)[j];
+  }
+  
   // and take out the content
-  long temp = diff_weight[0];
-
-  for (j=1; j<currRing->N && temp != 1; j++)
-  {
-    temp = gcd(temp, diff_weight[j]);
-    if (temp == 1) goto Finish;
-  }
-
-  for (j=0; j<currRing->N; j++)
-      diff_weight[j] = diff_weight[j] / temp;
-
-  Finish:
+  temp = (*diff_weight)[0];
+  if(temp != 1)
+    for (j=1; j<nRing; j++)
+    {
+      temp = gcd(temp, (*diff_weight)[j]);
+      if (temp == 1) 
+        return diff_weight;
+    }
+
+  for (j=0; j<nRing; j++)
+      (*diff_weight)[j] = (*diff_weight)[j] / temp;
+
   return diff_weight;
 }
 
-
-// next weight vector given weights as intvecs
-intvec* walkNextWeight(intvec* curr_weight, intvec* target_weight, ideal G)
-{
-  assume(curr_weight->length() == currRing->N);
-  assume(target_weight->length() == currRing->N);
-
-  int* nw = walkNextWeight(curr_weight->ivGetVec(),
-                           target_weight->ivGetVec(),
-                           G);
-  intvec* next_weight;
-
-  if (nw != NULL && nw != (int*) 1)
-  {
-    next_weight = new intvec(currRing->N);
-    int *nw_i = next_weight->ivGetVec();
-    int i;
-
-    for (i=0; i<currRing->N; i++)
-      nw_i[i] = nw[i];
-    omFreeSize(nw, (currRing->N)*sizeof(int));
-  }
-  else
-  {
-    next_weight = (intvec*) nw;
-  }
-
-  return next_weight;
-}
-
-
-// returns ideals of initials (w.r.t. curr_weight) of ideal G
-// assumes that monoms are ordered by descending W-degree (w.r.t curr_weight)
-
-poly walkInitials(poly p)
-{
-  assume(p != NULL);
-
-  poly pi = pHead(p);
-  poly pr = pi;
-  long d_lm = pGetOrder(p);
-
-  pIter(p);
-
-  while (p != NULL && pGetOrder(p) == d_lm)
-  {
-    pNext(pi) = pHead(p);
-    pIter(pi);
-    pIter(p);
-  }
-
-  assume(p == NULL || pGetOrder(p) < d_lm);
-
-  pNext(pi) = NULL;
-  pTest(pr);
-  return pr;
-}
-
-ideal walkInitials(ideal G)
-{
-  ideal GI = idInit(IDELEMS(G),1);
+/* Condition: poly f must be divided by the ideal G */
+/* Let f = h1g1+...+hsgs, then the result is (h1, ... ,hs) */
+static ideal MNormalForm(poly f, ideal G)
+{
+  int nG = IDELEMS(G);
+  ideal lmG = idInit(nG, 1);
+  ideal result = idInit(nG, 1);
+  int i, ind=0, ncheck=0;
+
+  for(i=0; i<nG; i++)
+  {
+    lmG->m[i] = pHead(G->m[i]);
+    result->m[i] = NULL;
+  }
+
+  poly h = f;
+  poly lmh, q, pmax = pISet(1), quot, qtmp=NULL;
+  int ntest = 0;
+  while(h != NULL)
+  {
+    lmh = pHead(h);
+    for(i=0; i<nG; i++)
+    {
+      q = MpDiv(lmh, lmG->m[i]);
+
+      if(q != NULL)
+      {
+          if(ncheck == 0)
+          {
+            result->m[i] = pCopy(pAdd(result->m[i], q));
+            h = pSub(h, pMult(q, pCopy(G->m[i])));
+            break;
+          }
+          else {
+            h = pSub(f, pMult(q, pCopy(G->m[i])));
+            if(quot != NULL)
+            {
+              ntest = 1;
+              qtmp = q;
+              ind = i;
+              ncheck = 0;
+            }
+          }
+      }
+    }
+    if(i==nG)
+    {
+      f = h;
+      pIter(h);
+      ncheck = 1;
+    }
+
+    if(ntest == 1) 
+    {
+      result->m[ind] = pCopy(pAdd(result->m[ind], qtmp));
+      ntest = 0; 
+    }
+  }
+  ideal rest = idCopy(result);
+  idDelete(&result);
+  idDelete(&lmG);
+  return rest;
+}
+
+/* GW is an initial form of the ideal G w.r.t. a weight vector */
+/* polynom f is divided by the ideal GW                        */
+/* Let f := h_1.gw_1 + ... + h_s.gw_s, then the result is      */
+/*          h_1.g_1 + ... + h_s.g_s                            */
+static poly MpolyConversion(poly f, ideal GW, ideal G)
+{
+  ideal H = MNormalForm(f, GW);
+  ideal HG = MidMultLift(H, G);
+
+  poly result = NULL;
   int i;
-
-  for (i=0; i<IDELEMS(G); i++)
-    GI->m[i] = walkInitials(G->m[i]);
-
-  return GI;
-}
+  int nG = IDELEMS(G);
+  
+  for(i=0; i<nG; i++)
+   result = pCopy(pAdd(result, HG->m[i]));
+
+  return result;
+}
+
+/* GW is an initial form of the ideal G w.r.t. a weight vector */
+/* Each polynom f of the ideal M is divided by the ideal GW    */
+/* Let m_i := h_1.gw_1 + ... + h_s.gw_s, then the i-th polynom */
+/* of result is  f_i := h_1.g_1 + ... + h_s.g_s                */
+ideal MidealConversion(ideal M, ideal GW, ideal G)
+{
+  int nM = IDELEMS(M);
+  int i;
+
+  for(i=0; i<nM; i++)
+  {
+    M->m[i] = MpolyConversion(M->m[i], GW, G);
+  }
+  ideal result = idCopy(M);
+  return result;
+}
+
+/* check that the monomial f is reduced by a monomial ideal G */
+static inline int MCheckpRedId(poly f, ideal G)
+{
+  int nG = IDELEMS(G);
+  int i;
+  poly q;
+
+  for(i=0; i<nG; i++)
+  {
+    q = MpDiv(f, G->m[i]);
+    if(q != NULL)
+      return 0; 
+  }
+  return 1;
+}
+
+poly MpReduceId(poly f, ideal GO)
+{
+  ideal G = idCopy(GO);
+  int nG = IDELEMS(G);
+  int i, pcheck;
+  ideal HG = idInit(nG, 1);
+
+  for(i=0; i<nG; i++)
+    HG->m[i] = pHead(G->m[i]);
+
+  poly h = f;
+  poly lmh, q,qg, result = NULL;
+
+  while(h!=NULL)
+  {
+    f = pCopy( h);
+    lmh = pHead(h); 
+
+    if(MCheckpRedId(lmh, HG) != 0)
+    {
+      result = pCopy(pAdd(result, lmh));
+      pIter(h);
+    }
+    else 
+      for(i=0; i<nG; i++)
+      {
+        q = MpDiv(lmh, HG->m[i]);
+        if(q != NULL)
+        {
+          f = pAdd(result, f);
+          qg = pMult(q, G->m[i]);
+          h = pSub(f, qg);
+          result = NULL;
+
+          lmh = pHead(h);
+          pcheck = MCheckpRedId(lmh, HG);
+          if(pcheck != 0)
+          {
+            break;
+          }
+        }
+      }
+ }
+ idDelete(&HG);
+ return result;
+}
+
+/* return f, if a head term of f is not divided by a head ideal M */
+static poly MpMinimId(poly f, ideal M)
+{
+  int nM = IDELEMS(M);
+  ideal HM = idInit(nM, 1);
+  int i, pcheck;
+
+  for(i=0; i<nM; i++)
+    HM->m[i] = pCopy(M->m[i]);
+
+  poly result = pCopy(f);
+  poly hf=pHead(f), q, qtmp, h=f; 
+
+  if(MCheckpRedId(pHead(f), HM) != 0)
+    goto FINISH;
+
+  while(1)
+  {
+    for(i=0; i<nM; i++)
+    {
+      q = MpDiv(hf, HM->m[i]);
+      if(q != NULL)
+        {
+          qtmp = pMult(q, M->m[i]);
+          h = pSub(h, qtmp);
+
+          hf = pHead(h);
+          pcheck = MCheckpRedId(hf, HM);
+          if(pcheck != 0)
+          { 
+            result = pCopy(h);
+            goto FINISH;           
+          }
+          break;
+        }
+    }
+ }
+
+ FINISH:
+ idDelete(&HM);
+ return result;
+}
+
+/* return a minimal ideal of an ideal M */
+ideal MidMinimId(ideal M)
+{
+  int i,j=0;
+  ideal result = idInit(IDELEMS(M),1);
+  poly pmin;
+  for(i=0; i<IDELEMS(M); i++)
+  {
+    ideal Mtmp = idCopy(M);
+    Mtmp->m[j] = NULL; 
+    idSkipZeroes(Mtmp);
+    pmin = MpMinimId(pCopy(M->m[i]), Mtmp);
+    M->m[i] = pCopy(pmin);    
+    result->m[j] = pmin;
+    if(pmin == NULL)
+    {
+      i--; 
+      j--;
+      idSkipZeroes(M);
+    }
+    j++;
+    idDelete(&Mtmp);
+  }
+  idSkipZeroes(result);
+  ideal res = idCopy(result);
+  idDelete(&result);
+  return res;
+}
+
+
+ideal MidMinBase(ideal G)
+{  
+  if(G == NULL)
+    return G;
+
+    intvec * wth;
+    lists re = min_std(G,currQuotient,(tHomog)TRUE,&wth,NULL,0,3);
+    G = (ideal)re->m[1].data;
+    re->m[1].data = NULL;
+    re->m[1].rtyp = NONE;
+    re->Clean();
+
+  return G;
+}
+
+
+/* compute a Groebner basis of a homogenoues ideal */
+ideal MNWstdhomRed(ideal G, intvec* iv)
+{
+
+  ideal Gomega = MwalkInitialForm(G, iv);
+  G = kStd(Gomega, NULL, isHomog, NULL);
+  Gomega = MkInterRed(G);
+
+  return Gomega;
+}
+
+/*****************************************************************************
+* If target_ord = intmat(A1,..., An) then calculate the perturbation vectors * 
+*     tau_p_dep = inveps^(p_deg-1)*A1 + inveps^(p_deg-2)*A2 +... + A_p_deg   *
+* where                                                                      *
+*     inveps > totaldegree(G)*(max(A2)+...+max(A_p_deg))                     *
+* and                                                                        *
+*     p_deg <= the number of variables of a basering                         *
+******************************************************************************/
+intvec* Mfivpert(ideal G, intvec* target, int p_deg)
+{
+  int i, j;
+  int nV = currRing->N;
+  
+  //Checking that the perturbation degree is valid
+  if(p_deg > nV || p_deg <= 0)
+  {
+    WerrorS("The perturbed degree is wrong!!");
+    return NULL;
+  }
+
+  // Calculate max_el_rows = Max(A2)+Max(A3)+...+Max(Ap_deg),
+  //    where the Ai are the rows of the target-order matrix.
+  int nmax = 0, maxAi, ntemp;
+
+  for(i=1; i < p_deg; i++) 
+  {
+    maxAi = (*target)[i*nV];
+    for(j=1; j < nV; j++) 
+    {
+      ntemp = (*target)[i*nV + j];
+      if(ntemp > maxAi) 
+        maxAi = ntemp;
+    }
+    nmax += maxAi;
+  }
+
+  // Calculate inv_eps := 1/eps, where 1/eps > deg(p)*max_el_rows
+  //        for all p in G.
+  int inv_eps, degG, max_deg=0;
+  intvec* ivUnit = Mivdp(nV);
+  
+  for (i=0; i<IDELEMS(G); i++) 
+  {
+    degG = MwalkWeightDegree(G->m[i], ivUnit);
+    if(degG > max_deg) 
+      max_deg = degG;
+  }
+  inv_eps = (max_deg * nmax) + 1;
+
+  
+  // Calculate the perturbed target order:
+  // Since a weight vector in Singular has to be in integral, we compute
+  // tau_p_deg = inv_eps^(p_deg-1)*A1 - inv_eps^(p_deg-2)*A2+...+A_p_deg,
+
+  intvec* ivtemp = new intvec(nV);  
+  intvec* pert_vector = new intvec(nV);
+
+  for(i=0; i<nV; i++)
+  {
+    ntemp = (*target)[i];
+    (* ivtemp)[i] = ntemp;
+    (* pert_vector)[i] = ntemp;
+  }
+
+  for(i=1; i<p_deg; i++)
+  {
+    for(j=0; j<nV; j++)
+      (* ivtemp)[j] = inv_eps*(*ivtemp)[j] + (*target)[i*nV+j];
+
+    pert_vector = ivtemp;  
+  }
+  omFree((ADDRESS) ivtemp);
+  return pert_vector;
+}
+
+ideal MpHeadIdeal(ideal G)
+{
+  int i, nG = IDELEMS(G);
+  ideal result = idInit(nG,1);
+
+  for(i=0; i<nG; i++)
+  {
+    result->m[i] = pHead(G->m[i]);
+  }
+
+  ideal res = idCopy(result);
+  idDelete(&result);
+  return res;
+}
+
+void* checkideal(ideal G)
+{
+  int i;
+  printf("//** Size(G)= %d, and ", IDELEMS(G));
+
+  for(i=0; i<IDELEMS(G); i++)
+  {
+    printf("G[%d] = %d, ", i, pLength(G->m[i]));
+  }
+  printf("\n");
+}
+