Changeset 9c0326 in git

Timestamp:

Sep 10, 2014, 5:43:31 PM (9 years ago)

Author:

Yue Ren <ren@…>

Branches:

(u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')

Children:

b71400abefd0c9c66c0009063f40ddd5b1f1c4c7

Parents:

eacb78146504ca284c1888c739cb1aade67c8520

git-author:

Yue Ren <ren@mathematik.uni-kl.de>2014-09-10 17:43:31+02:00

git-committer:

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

Message:

chg: status update 12.09.

Files:

• Singular/dyn_modules/gfanlib/bbcone.cc (modified) (2 diffs)
• Singular/dyn_modules/gfanlib/flip.cc (modified) (2 diffs)
• Singular/dyn_modules/gfanlib/flip.h (modified) (1 diff)
• Singular/dyn_modules/gfanlib/groebnerCone.cc (modified) (19 diffs)
• Singular/dyn_modules/gfanlib/groebnerCone.h (modified) (3 diffs)
• Singular/dyn_modules/gfanlib/initial.cc (modified) (1 diff)
• Singular/dyn_modules/gfanlib/lift.cc (modified) (1 diff)
• Singular/dyn_modules/gfanlib/startingCone.cc (modified) (1 diff)
• Singular/dyn_modules/gfanlib/tropicalStrategy.cc (modified) (5 diffs)
• Singular/dyn_modules/gfanlib/tropicalStrategy.h (modified) (4 diffs)
• Singular/dyn_modules/gfanlib/witness.cc (modified) (3 diffs)
• Singular/dyn_modules/gfanlib/witness.h (modified) (1 diff)
• gfanlib/Makefile.in.save (deleted)
• gfanlib/gfanlib_matrix.h (modified) (2 diffs)
• gfanlib/gfanlib_vector.h (modified) (1 diff)

Singular/dyn_modules/gfanlib/bbcone.cc

@@ -1607 +1607 @@
   WerrorS("facetContaining: unexpected parameters");
   return TRUE;
+}
+
+
+BOOLEAN memleaktest(leftv res, leftv args)
+{
+  bigintmat* ineq = NULL;
+  if (args->Typ() == INTMAT_CMD)
+  {
+    intvec* ineq0 = (intvec*) args->Data();
+    ineq = iv2bim(ineq0,coeffs_BIGINT);
+  }
+  else
+    ineq = (bigintmat*) args->Data();
+  gfan::ZMatrix* zm = bigintmatToZMatrix(ineq);
+  gfan::ZCone* zc = new gfan::ZCone(*zm, gfan::ZMatrix(0, zm->getWidth()));
+  delete zm;
+  if (args->Typ() == INTMAT_CMD)
+    delete ineq;
+  res->rtyp = coneID;
+  res->data = (void*) zc;
+  return FALSE;
 }
 
@@ -1664 +1685 @@
   p->iiAddCproc("","listOfFacets",FALSE,listOfFacets);
   p->iiAddCproc("","facetContaining",FALSE,facetContaining);
+  p->iiAddCproc("","memleaktest",FALSE,memleaktest);
   coneID=setBlackboxStuff(b,"cone");
 }

Singular/dyn_modules/gfanlib/flip.cc

@@ -38 +38 @@
  *  (2) Is is a Groebner basis of the same ideal with respect to the ordering on s
  **/
-std::pair<ideal,ring> flip0(const ideal I, const ring r,
-                            const gfan::ZVector interiorPoint,
-                            const gfan::ZVector facetNormal,
-                            const gfan::ZVector adjustedInteriorPoint,
-                            const gfan::ZVector adjustedFacetNormal)
+std::pair<ideal,ring> flip(const ideal I, const ring r,
+                           const gfan::ZVector interiorPoint,
+                           const gfan::ZVector facetNormal,
+                           const gfan::ZVector adjustedInteriorPoint,
+                           const gfan::ZVector adjustedFacetNormal)
 {
   /* create a ring with weighted ordering  */
@@ -78 +78 @@
   rChangeCurrRing(sAdjusted);
   ideal inIsAdjustedGB = kStd(inIsAdjusted,currRing->qideal,testHomog,&nullVector);
+  ideal IsAdjustedGB = lift(I,r,inIsAdjustedGB,sAdjusted);
   id_Delete(&inIsAdjusted,sAdjusted);
-  ideal IsAdjustedGB = lift(I,r,inIsAdjustedGB,sAdjusted);
   id_Delete(&inIsAdjustedGB,sAdjusted);
 

Singular/dyn_modules/gfanlib/flip.h

@@ -6 +6 @@
 #include <gfanlib/gfanlib_vector.h>
 
-std::pair<ideal,ring> flip0(const ideal I, const ring r,
-                            const gfan::ZVector interiorPoint,
-                            const gfan::ZVector facetNormal,
-                            const gfan::ZVector adjustedInteriorPoint,
-                            const gfan::ZVector adjustedFacetNormal);
+std::pair<ideal,ring> flip(const ideal I, const ring r,
+                           const gfan::ZVector interiorPoint,
+                           const gfan::ZVector facetNormal,
+                           const gfan::ZVector adjustedInteriorPoint,
+                           const gfan::ZVector adjustedFacetNormal);
 
 #endif

Singular/dyn_modules/gfanlib/groebnerCone.cc

@@ -12 +12 @@
 #include <gfanlib/gfanlib_matrix.h>
 
+#include <initial.h>
 #include <tropicalStrategy.h>
 #include <groebnerCone.h>
@@ -57 +58 @@
 groebnerCone::groebnerCone():
   polynomialIdeal(NULL),
+  reducedPolynomialIdeal(NULL),
+  initialPolynomialIdeal(NULL),
   polynomialRing(NULL),
   polyhedralCone(gfan::ZCone(0)),
@@ -66 +69 @@
 groebnerCone::groebnerCone(const ideal I, const ring r, const tropicalStrategy& currentCase):
   polynomialIdeal(NULL),
+  reducedPolynomialIdeal(NULL),
+  initialPolynomialIdeal(NULL),
   polynomialRing(NULL),
   currentStrategy(&currentCase)
 {
   assume(checkPolynomialInput(I,r));
-  if (I) polynomialIdeal = id_Copy(I,r);
+  if (I)
+  {
+    polynomialIdeal = id_Copy(I,r);
+    reducedPolynomialIdeal = id_Copy(I,r);
+  }
   if (r) polynomialRing = rCopy(r);
-  currentCase.reduce(I,r);
+  currentCase.reduce(reducedPolynomialIdeal,polynomialRing);
 
   int n = rVar(polynomialRing);
@@ -81 +90 @@
   gfan::ZVector tailexpw = gfan::ZVector(n);
   gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
-  for (int i=0; i<idSize(polynomialIdeal); i++)
-  {
-    g = I->m[i];
+  for (int i=0; i<idSize(reducedPolynomialIdeal); i++)
+  {
+    g = reducedPolynomialIdeal->m[i];
     if (g)
     {
@@ -110 +119 @@
   polyhedralCone.canonicalize();
   interiorPoint = polyhedralCone.getRelativeInteriorPoint();
+  initialPolynomialIdeal = sloppyInitial(reducedPolynomialIdeal,polynomialRing,interiorPoint);
   assume(checkOrderingAndCone(polynomialRing,polyhedralCone));
 }
@@ -122 +132 @@
 groebnerCone::groebnerCone(const ideal I, const ring r, const gfan::ZVector& w, const tropicalStrategy& currentCase):
   polynomialIdeal(NULL),
+  reducedPolynomialIdeal(NULL),
+  initialPolynomialIdeal(NULL),
   polynomialRing(NULL),
   currentStrategy(&currentCase)
 {
   assume(checkPolynomialInput(I,r));
-  if (I) polynomialIdeal = id_Copy(I,r);
+  if (I)
+  {
+    polynomialIdeal = id_Copy(I,r);
+    reducedPolynomialIdeal = id_Copy(I,r);
+  }
   if (r) polynomialRing = rCopy(r);
-  currentCase.reduce(I,r);
-
-  int n = rVar(r);
+  currentCase.reduce(reducedPolynomialIdeal,polynomialRing);
+
+  int n = rVar(polynomialRing);
   gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
   gfan::ZMatrix equations = gfan::ZMatrix(0,n);
   int* expv = (int*) omAlloc((n+1)*sizeof(int));
-  for (int i=0; i<idSize(polynomialIdeal); i++)
-  {
-    poly g = polynomialIdeal->m[i];
+  for (int i=0; i<idSize(reducedPolynomialIdeal); i++)
+  {
+    poly g = reducedPolynomialIdeal->m[i];
     if (g)
     {
@@ -167 +183 @@
   polyhedralCone.canonicalize();
   interiorPoint = polyhedralCone.getRelativeInteriorPoint();
+  initialPolynomialIdeal = sloppyInitial(reducedPolynomialIdeal,polynomialRing,interiorPoint);
   assume(checkOrderingAndCone(polynomialRing,polyhedralCone));
 }
@@ -176 +193 @@
 groebnerCone::groebnerCone(const ideal I, const ring r, const gfan::ZVector& u, const gfan::ZVector& w, const tropicalStrategy& currentCase):
   polynomialIdeal(NULL),
+  reducedPolynomialIdeal(NULL),
+  initialPolynomialIdeal(NULL),
   polynomialRing(NULL),
   currentStrategy(&currentCase)
 {
   assume(checkPolynomialInput(I,r));
-  if (I) polynomialIdeal = id_Copy(I,r);
+  if (I)
+  {
+    polynomialIdeal = id_Copy(I,r);
+    reducedPolynomialIdeal = id_Copy(I,r);
+  }
   if (r) polynomialRing = rCopy(r);
-  currentCase.reduce(I,r);
-
-  int n = rVar(r);
+  currentCase.reduce(reducedPolynomialIdeal,polynomialRing);
+
+  int n = rVar(polynomialRing);
   gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
   gfan::ZMatrix equations = gfan::ZMatrix(0,n);
   int* expv = (int*) omAlloc((n+1)*sizeof(int));
-  for (int i=0; i<idSize(polynomialIdeal); i++)
-  {
-    poly g = polynomialIdeal->m[i];
+  for (int i=0; i<idSize(reducedPolynomialIdeal); i++)
+  {
+    poly g = reducedPolynomialIdeal->m[i];
     if (g)
     {
@@ -222 +245 @@
   polyhedralCone.canonicalize();
   interiorPoint = polyhedralCone.getRelativeInteriorPoint();
+  initialPolynomialIdeal = sloppyInitial(reducedPolynomialIdeal,polynomialRing,interiorPoint);
   assume(checkOrderingAndCone(polynomialRing,polyhedralCone));
 }
@@ -228 +252 @@
 groebnerCone::groebnerCone(const ideal I, const ideal inI, const ring r, const tropicalStrategy& currentCase):
   polynomialIdeal(id_Copy(I,r)),
+  reducedPolynomialIdeal(id_Copy(I,r)),
+  initialPolynomialIdeal(id_Copy(inI,r)),
   polynomialRing(rCopy(r)),
   currentStrategy(&currentCase)
@@ -233 +259 @@
   assume(checkPolynomialInput(I,r));
   assume(checkPolynomialInput(inI,r));
+  currentStrategy->reduce(reducedPolynomialIdeal,polynomialRing);
+  // currentStrategy->reduce(initialPolynomialIdeal,polynomialRing);
+
   int n = rVar(r);
   gfan::ZMatrix equations = gfan::ZMatrix(0,n);
   int* expv = (int*) omAlloc((n+1)*sizeof(int));
-  for (int i=0; i<idSize(inI); i++)
-  {
-    poly g = inI->m[i];
+  for (int i=0; i<idSize(initialPolynomialIdeal); i++)
+  {
+    poly g = initialPolynomialIdeal->m[i];
     if (g)
     {
@@ -252 +281 @@
   }
   gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
-  for (int i=0; i<idSize(I); i++)
-  {
-    poly g = I->m[i];
+  for (int i=0; i<idSize(reducedPolynomialIdeal); i++)
+  {
+    poly g = reducedPolynomialIdeal->m[i];
     if (g)
     {
@@ -283 +312 @@
 groebnerCone::groebnerCone(const groebnerCone &sigma):
   polynomialIdeal(NULL),
+  reducedPolynomialIdeal(NULL),
+  initialPolynomialIdeal(NULL),
   polynomialRing(NULL),
   polyhedralCone(gfan::ZCone(sigma.getPolyhedralCone())),
@@ -292 +323 @@
   assume(checkPolyhedralInput(sigma.getPolyhedralCone(),sigma.getInteriorPoint()));
   if (sigma.getPolynomialIdeal()) polynomialIdeal = id_Copy(sigma.getPolynomialIdeal(),sigma.getPolynomialRing());
+  if (sigma.getReducedPolynomialIdeal()) reducedPolynomialIdeal = id_Copy(sigma.getReducedPolynomialIdeal(),sigma.getPolynomialRing());
+  if (sigma.getInitialPolynomialIdeal()) initialPolynomialIdeal = id_Copy(sigma.getInitialPolynomialIdeal(),sigma.getPolynomialRing());
   if (sigma.getPolynomialRing()) polynomialRing = rCopy(sigma.getPolynomialRing());
 }
@@ -333 +366 @@
   {
     assume(!currentStrategy->restrictToLowerHalfSpace() || R[i][0].sign()<=0);
+    if (currentStrategy->restrictToLowerHalfSpace() && R[i][0].sign()==0)
+      continue;
     poly s = currentStrategy->checkInitialIdealForMonomial(I,r,R[i]);
     if (s==NULL)
@@ -355 +390 @@
     return false;
   }
+  if (polyhedralCone.containsRelatively(interiorPoint))
+  {
+    std::cout << "ERROR: interiorPoint is contained in the interior of the maximal Groebner cone!" << std::endl
+              << "cone: " << std::endl
+              << toString(&polyhedralCone)
+              << "interiorPoint:" << std::endl
+              << interiorPoint << std::endl;
+    return false;
+  }
   gfan::ZCone hopefullyAFacet = polyhedralCone.faceContaining(interiorPoint);
   if (hopefullyAFacet.dimension()!=(polyhedralCone.dimension()-1))
@@ -379 +423 @@
 }
 
-/***
+/**
  * Given an interior point on the facet and the outer normal factor on the facet,
  * returns the adjacent groebnerCone sharing that facet
- **/
+ */
 groebnerCone groebnerCone::flipCone(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const
 {
@@ -391 +435 @@
    *   to obtain an initial form with respect to interiorPoint+e*facetNormal,
    *   for e>0 sufficiently small */
-  std::pair<ideal,ring> flipped = currentStrategy->flip(polynomialIdeal,polynomialRing,interiorPoint,facetNormal);
-  std::cout << "finished flip" << std::endl;
+  std::pair<ideal,ring> flipped = currentStrategy->flip(reducedPolynomialIdeal,polynomialRing,interiorPoint,facetNormal);
   assume(checkPolynomialInput(flipped.first,flipped.second));
   groebnerCone flippedCone(flipped.first, flipped.second, interiorPoint, facetNormal, *currentStrategy);
@@ -464 +507 @@
   for (int i=0; i<interiorPoints.getHeight(); i++)
   {
-    std::cout << "performing flip " << i << " out of " << interiorPoints.getHeight() << std::endl;
+    gfan::ZVector w = interiorPoints[i];
+    gfan::ZVector v = facetNormals[i];
+    if (currentStrategy->restrictToLowerHalfSpace())
+    {
+      assume(w[0].sign()<=0);
+      if (w[0].sign()==0 && v[0].sign()>0)
+        continue;
+    }
     neighbours.insert(this->flipCone(interiorPoints[i],facetNormals[i]));
   }

Singular/dyn_modules/gfanlib/groebnerCone.h

@@ -35 +35 @@
   ideal polynomialIdeal;
   /**
+   * reduced form of the ideal from which the equations and inequalities can be read
+   */
+  ideal reducedPolynomialIdeal;
+  /**
+   * generators of the initial ideal
+   */
+  ideal initialPolynomialIdeal;
+  /**
    * ring in which the ideal exists
    */
@@ -40 +48 @@
   gfan::ZCone polyhedralCone;
   gfan::ZVector interiorPoint;
-  gfan::ZVector negativePointInLineality;
   const tropicalStrategy* currentStrategy;
 
@@ -67 +74 @@
 
   ideal getPolynomialIdeal() const { return polynomialIdeal; };
+  ideal getReducedPolynomialIdeal() const { return reducedPolynomialIdeal; };
+  ideal getInitialPolynomialIdeal() const { return initialPolynomialIdeal; };
   ring getPolynomialRing() const { return polynomialRing; };
   gfan::ZCone getPolyhedralCone() const { return polyhedralCone; };
   gfan::ZVector getInteriorPoint() const { return interiorPoint; };
-  gfan::ZVector getNegativePointInLineality() const { return negativePointInLineality; };
   const tropicalStrategy* getTropicalStrategy() const {return currentStrategy; };
 

Singular/dyn_modules/gfanlib/initial.cc

@@ -39 +39 @@
 }
 
-/***
+/**
  * Checks if p is sorted with respect to w.
- **/
+ */
 static bool checkSloppyInput(const poly p, const ring r, const gfan::ZVector w)
 {

Singular/dyn_modules/gfanlib/lift.cc

@@ -9 +9 @@
  * the ordering of s.
  **/
-ideal lift(const ideal I, const ring r, const ideal inI, const ring s)
+// ideal lift(const ideal I, const ring r, const ideal inI, const ring s)
+// {
+//   nMapFunc identity = n_SetMap(r->cf,s->cf);
+//   int k = idSize(I); ideal Is = idInit(k);
+//   for (int i=0; i<k; i++)
+//     Is->m[i] = p_PermPoly(I->m[i],NULL,r,s,identity,NULL,0);
+//   ideal J = idInit(k);
+//   for (int i=0; i<k; i++)
+//     J->m[i] = witness(inI->m[i],Is,inI,s);
+//   id_Delete(&Is,s);
+//   return J;
+// }
+
+ideal lift(const ideal J, const ring r, const ideal inI, const ring s)
 {
-  nMapFunc identity = n_SetMap(r->cf,s->cf);
-  int k = idSize(I); ideal Is = idInit(k);
+  nMapFunc identity = n_SetMap(s->cf,r->cf);
+  int k = idSize(inI); ideal inIr = idInit(k);
   for (int i=0; i<k; i++)
-    Is->m[i] = p_PermPoly(I->m[i],NULL,r,s,identity,NULL,0);
-  ideal J = idInit(k);
+    inIr->m[i] = p_PermPoly(inI->m[i],NULL,s,r,identity,NULL,0);
+  ideal Ir = witness(inIr,J,r);
+
+  identity = n_SetMap(r->cf,s->cf);
+  ideal Is = idInit(k);
   for (int i=0; i<k; i++)
-    J->m[i] = witness(inI->m[i],Is,inI,s);
-  id_Delete(&Is,s);
-  return J;
+    Is->m[i] = p_PermPoly(Ir->m[i],NULL,r,s,identity,NULL,0);
+
+  id_Delete(&inIr,r);
+  id_Delete(&Ir,r);
+  return Is;
 }

Singular/dyn_modules/gfanlib/startingCone.cc

@@ -457 +457 @@
 
     ideal J = lift(I,r,inI,s);
+    // currentStrategy.reduce(J,s);
     groebnerCone startingCone(J,inI,s,currentStrategy);
     id_Delete(&inI,s);

Singular/dyn_modules/gfanlib/tropicalStrategy.cc

@@ -4 +4 @@
 #include <ttinitialReduction.h>
 #include <tropical.h>
+#include <std_wrapper.h>
 
 // for various commands in dim(ideal I, ring r):
@@ -23 +24 @@
   {
     int i = idPosConstant(I);
-    if ((i != -1) && (n_IsUnit(pGetCoeff(I->m[i]),currRing->cf)))
+    if ((i != -1) && (n_IsUnit(p_GetCoeff(I->m[i],currRing->cf),currRing->cf)))
       return -1;
     ideal vv = id_Head(I,currRing);
@@ -66 +67 @@
 {
   assume(rField_is_Q(r) || rField_is_Zp(r));
+  startingIdeal = gfanlib_kStd_wrapper(startingIdeal,startingRing);
   if (!completelyHomogeneous)
   {
@@ -174 +176 @@
     startingIdeal->m[i+1] = p_PermPoly(J->m[i],shiftByOne,originalRing,startingRing,nMap,NULL,0);
   omFreeSize(shiftByOne,(n+1)*sizeof(int));
+  startingIdeal = gfanlib_kStd_wrapper(startingIdeal,startingRing);
 
   linealitySpace = homogeneitySpace(startingIdeal,startingRing);
@@ -247 +250 @@
   return *this;
 }
+
+bool tropicalStrategy::reduce(ideal I, const ring r) const
+{
+  rTest(r);
+  id_Test(I,r);
+  nMapFunc nMap = n_SetMap(startingRing->cf,r->cf);
+  number p = nMap(uniformizingParameter,startingRing->cf,r->cf);
+  bool b = this->extraReductionAlgorithm(I,r,p);
+  n_Delete(&p,r->cf);
+  return b;
+}
+
+ring tropicalStrategy::getShortcutRingPrependingWeight(const ring r, const gfan::ZVector v) const
+{
+  ring rShortcut = rCopy0(r);
+
+  // save old ordering
+  int* order = rShortcut->order;
+  int* block0 = rShortcut->block0;
+  int* block1 = rShortcut->block1;
+  int** wvhdl = rShortcut->wvhdl;
+
+  // adjust weight and create new ordering
+  gfan::ZVector w = adjustWeightForHomogeneity(v);
+  int h = rBlocks(r); int n = rVar(r);
+  rShortcut->order = (int*) omAlloc0((h+1)*sizeof(int));
+  rShortcut->block0 = (int*) omAlloc0((h+1)*sizeof(int));
+  rShortcut->block1 = (int*) omAlloc0((h+1)*sizeof(int));
+  rShortcut->wvhdl = (int**) omAlloc0((h+1)*sizeof(int*));
+  rShortcut->order[0] = ringorder_a;
+  rShortcut->block0[0] = 1;
+  rShortcut->block1[0] = n;
+  bool overflow;
+  rShortcut->wvhdl[0] = ZVectorToIntStar(w,overflow);
+  for (int i=1; i<=h; i++)
+  {
+    rShortcut->order[i] = order[i-1];
+    rShortcut->block0[i] = block0[i-1];
+    rShortcut->block1[i] = block1[i-1];
+    rShortcut->wvhdl[i] = wvhdl[i-1];
+  }
+
+  // if valuation non-trivial, change coefficient ring to residue field
+  if (valuationIsNonTrivial())
+  {
+    nKillChar(rShortcut->cf);
+    rShortcut->cf = nCopyCoeff(shortcutRing->cf);
+  }
+  rComplete(rShortcut);
+  rTest(rShortcut);
+
+  // delete old ordering
+  omFree(order);
+  omFree(block0);
+  omFree(block1);
+  omFree(wvhdl);
+
+  return rShortcut;
+}
+
+poly tropicalStrategy::checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector w) const
+{
+  // prepend extra weight for homogeneity
+  // switch to residue field if valuation is non trivial
+  ring rShortcut = getShortcutRingPrependingWeight(r,w);
+
+  // compute the initial ideal and map it into the constructed ring
+  // if switched to residue field, remove possibly 0 elements
+  ideal inI = initial(I,r,w);
+  int k = idSize(inI);
+  ideal inIShortcut = idInit(k);
+  nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
+  for (int i=0; i<k; i++)
+    inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
+  if (isValuationNonTrivial())
+    idSkipZeroes(inIShortcut);
+
+  // check initial ideal for monomial and
+  // if it exsists, return a copy of the monomial in the input ring
+  poly p = checkForMonomialViaSuddenSaturation(inIShortcut,rShortcut);
+  poly monomial = NULL;
+  if (p!=NULL)
+  {
+    monomial=p_One(r);
+    for (int i=1; i<n; i++)
+      p_SetExp(monomial,i,p_GetExp(p,i,rShortcut),r);
+    p_Delete(&p,rShortcut);
+  }
+  id_Delete(&inI,r);
+  id_Delete(&inIShortcut,rShortcut);
+  rDelete(rShortcut);
+  return monomial;
+}
+
+ring tropicalStrategy::copyAndChangeCoefficientRing(const ring r) const
+{
+  ring rShortcut = rCopy0(r);
+  nKillChar(rShortcut->cf);
+  rShortcut->cf = nCopyCoeff(shortcutRing->cf);
+  rComplete(rShortcut);
+  rTest(rShortcut);
+  return rShortcut;
+}
+
+ideal getWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
+{
+  // if the valuation is trivial and the ring and ideal have not been extended,
+  // then it is sufficient to return the difference between the elements of inJ
+  // and their normal forms with respect to I and r
+  if (isValuationTrivial())
+    return witness(inJ,I,r);
+  // if the valuation is non-trivial and the ring and ideal have been extended,
+  // then we can make a shortcut through the residue field
+  else
+  {
+    checkFirstGenerator(I,r);
+    assume(idSize(inI)==idSize(I));
+    /**
+     * change ground ring into finite field
+     * and map the data into it
+     */
+    ring rShortcut = copyAndChangeCoefficientRing(r);
+
+    int k = idSize(inI);
+    int l = idSize(I);
+    ideal inJShortcut = idInit(k);
+    ideal inIShortcut = idInit(l);
+    nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
+    for (int i=0; i<k; i++)
+      inJShortcut->m[i] = p_PermPoly(inJ->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
+    for (int j=0; j<l; j++)
+      inIShortcut->m[j] = p_PermPoly(inI->m[j],NULL,r,rShortcut,takingResidues,NULL,0);
+    id_Test(inJShortcut,rShortcut);
+    id_Test(inIShortcut,rShortcut);
+
+    /**
+     * Compute a division with remainder over the finite field
+     * and map the result back to r
+     */
+    matrix QShortcut = divisionDiscardingRemainder(inJShortcut,inIShortcut,rShortcut);
+    matrix Q = mpNew(k,l);
+    nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
+    for (int ij=k*l-1; ij>=0; ij--)
+      Q->m[ij] = p_PermPoly(Q->m[ij],NULL,rShortcut,r,takingRepresentatives,NULL,0);
+
+    /**
+     * Compute the normal forms
+     */
+    ideal J = idInit(k);
+    for (int j=0; j<l; j++)
+    {
+      poly q0 = p_Copy(inJ->m[j],r);
+      for (int i=1; i<k; i++)
+        q0 = p_Add_q(q0,p_Neg(p_Mult_q(p_Copy(MATELEM(Q,i,j),r),p_Copy(inI->m[i],r),r),r),r);
+      J->m[j] = p_Mult_q(q0,p_Copy(J->m[0],r),r);
+      q0 = NULL;
+      for (int i=1; i<k; i++)
+        J->m[j] = p_Add_q(I->m[j],p_Mult_q(p_Copy(MATELEM(Q,i,j),r),p_Copy(I->m[i],r),r),r);
+    }
+
+    return J;
+  }
+}
+
+ideal tropicalStrategy::getStdOfInitialIdeal(const ideal inI, const ring r) const
+{
+  ideal inJ;
+  if (isValuationTrivial())
+    // if valuation trivial, then compute std as usual
+    inJ = gfanlib_kStd_wrapper(inI,r);
+  else
+  {
+    // if valuation non-trivial, then uniformizing parameter is in ideal
+    // so switch to residue field first and compute standard basis over the residue field
+    ring rShortcut = copyAndChangeCoefficientRing(r);
+    nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
+    int k = idSize(inI);
+    ideal inIShortcut = idInit(k);
+    for (int i=0; i<k; i++)
+      inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
+    ideal inJShortcut = gfanlib_kStd_wrapper(inIShortcut,rShortcut);
+
+    // and lift the result back to the ring with valuation
+    nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
+    k = idSize(inJShortcut);
+    inJ = idInit(k+1);
+    inJ->m[0] = p_One(r);
+    p_SetCoeff(inJ->m[0],nMap(uniformizingParameter,startingRing->cf,r->cf),r);
+    for (int i=0; i<k; i++)
+      inJ->m[i+1] = p_PermPoly(inJShortcut->m[i],NULL,rShortcut,r,takingRepresentatives,NULL,0);
+  }
+  return inJ;
+}
+
+ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector w, const gfan::ZVector v)
+{
+  // copy shortcutRing and change to desired ordering
+  bool ok;
+  ring s = rCopy0(r);
+  int n = rVar(s);
+  deleteOrdering(s);
+  gfan::ZVector wAdjusted = adjustWeightForHomogeneity(w);
+  gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(wAdjusted,v);
+  s->order = (int*) omAlloc0(4*sizeof(int));
+  s->block0 = (int*) omAlloc0(4*sizeof(int));
+  s->block1 = (int*) omAlloc0(4*sizeof(int));
+  s->wvhdl = (int**) omAlloc0(4*sizeof(int**));
+  s->order[0] = ringorder_a;
+  s->block0[0] = 1;
+  s->block1[0] = n;
+  s->wvhdl[0] = ZVectorToIntStar(adjustedInteriorPoint,ok);
+  s->order[1] = ringorder_wp;
+  s->block0[1] = 1;
+  s->block1[1] = n;
+  s->wvhdl[1] = ZVectorToIntStar(adjustedFacetNormal,ok);
+  s->order[2] = ringorder_C;
+  rComplete(s);
+  rTest(s);
+
+  return s;
+}
+
+ring tropicalStrategy::copyAndChangeOrderingDP(const ring r, const gfan::ZVector w, const gfan::ZVector v)
+{
+  // copy shortcutRing and change to desired ordering
+  bool ok;
+  ring s = rCopy0(r);
+  int n = rVar(s);
+  deleteOrdering(s);
+  gfan::ZVector wAdjusted = adjustWeightForHomogeneity(w);
+  gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(wAdjusted,v);
+  s->order = (int*) omAlloc0(5*sizeof(int));
+  s->block0 = (int*) omAlloc0(5*sizeof(int));
+  s->block1 = (int*) omAlloc0(5*sizeof(int));
+  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
+  s->order[0] = ringorder_a;
+  s->block0[0] = 1;
+  s->block1[0] = n;
+  s->wvhdl[0] = ZVectorToIntStar(adjustedInteriorPoint,ok);
+  s->order[1] = ringorder_a;
+  s->block0[1] = 1;
+  s->block1[1] = n;
+  s->wvhdl[1] = ZVectorToIntStar(adjustedFacetNormal,ok);
+  s->order[2] = ringorder_dp;
+  s->block0[2] = 1;
+  s->block1[2] = n;
+  s->order[3] = ringorder_C;
+  rComplete(s);
+  rTest(s);
+
+  return s;
+}
+
+
+std::pair<ideal,ring> tropicalStrategy::getFlip(const ideal I, const ring r,
+                                                const gfan::ZVector interiorPoint,
+                                                const gfan::ZVector facetNormal) const
+{
+  assume(isValuationTrivial() || interiorPoint[0].sign()<0);
+  assume(checkFirstGenerators(I,r));
+
+  // get a generating system of the initial ideal
+  // and compute a standard basis with respect to adjacent ordering
+  ideal inIr = initial(I,r,interiorPoint);
+  ring sAdjusted = copyAndChangeOrderingWP(r,interiorPoint,facetNormal);
+  nMapFunc identity = n_SetMap(r->cf,sAdjusted->cf);
+  int k = idSize(I); ideal inIs = idInit(k);
+  for (int i=0; i<k; i++)
+    inIsAdjusted->m[i] = p_PermPoly(inIr->m[i],NULL,r,sAdjusted,identity,NULL,0);
+  ideal inJsAdjusted = getStdOfInitialIdeal(inIsAdjusted,sShortcut);
+
+  // find witnesses of the new standard basis elements of the initial ideal
+  // with the help of the old standard basis of the ideal
+  k = idSize(inJsAdjusted);
+  ideal inJr = idInit(k);
+  identity = n_SetMap(sAdjusted->cf,r->cf);
+  for (int i=0; i<k; i++)
+    inJr->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,sAdjusted,r,identity,NULL,0);
+
+  ideal Jr = getWitness(inJr,inIr,I,r);
+  ring s = copyAndChangeOrderingDP(r,interiorPoint,facetNormal);
+  identity = n_SetMap(r->cf,s->cf);
+  for (int i=0; i<k; i++)
+    Js->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,r,s,identity,NULL,0);
+
+  return std::make_pair(Js,s);
+}

Singular/dyn_modules/gfanlib/tropicalStrategy.h

@@ -11 +11 @@
 #include <flip.h>
 #include <initial.h>
+#include <witness.h>
 
 /** \file
@@ -119 +120 @@
     return b;
   }
+  bool isValuationTrivial() const
+  {
+    bool b = (uniformizingParameter==NULL);
+    return b;
+  }
+  bool isValuationNonTrivial() const
+  {
+    bool b = (uniformizingParameter!=NULL);
+    return b;
+  }
 
   /**
@@ -224 +235 @@
 
   /**
+   * If valuation trivial, returns a copy of r with a positive weight prepended,
+   * such that any ideal homogeneous with respect to w is homogeneous with respect to that weight.
+   * If valuation non-trivial, changes the coefficient ring to the residue field.
+   */
+  ring getShortcutRingPrependingWeight(const ring r, const gfan::ZVector w) const;
+
+  /**
+   * Returns a copy of the shortcutRing which is weighted with respect to w first and v second.
+   * If valuation trivial, the coefficient ring will be the field with valuation.
+   * If valuation non-trivial, the coefficient ring will be the residue field.
+   */
+  ring getShortcutRingDoublyWeighted(const gfan::ZVector w, const gfan::ZVector v) const;
+
+  /**
+   * copies r and prepends extra weight w
+   */
+  ring copyAndPrependWeight(const ring r, const gfan::ZVector w) const;
+
+  /**
+   * changes the coefficient ring of r to be the residue field
+   */
+  void changeCoefficientToResidueField(ring r) const;
+
+
+
+  /**
    * reduces the generators of an ideal I so that
    * the inequalities and equations of the Groebner cone can be read off.
    */
-  bool reduce(ideal I, const ring r) const
-  {
-    rTest(r);  id_Test(I,r);
-    nMapFunc nMap = n_SetMap(startingRing->cf,r->cf);
-    number p = nMap(uniformizingParameter,startingRing->cf,r->cf);
-    bool b = this->extraReductionAlgorithm(I,r,p);
-    n_Delete(&p,r->cf);
-    return b;
-  }
+  bool reduce(ideal I, const ring r) const;
 
   /**
@@ -241 +270 @@
    * returns false otherwise.
    **/
-  poly checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector w) const
-  {
-    gfan::ZVector v = adjustWeightForHomogeneity(w);
-    if (isConstantCoefficientCase())
-    {
-      ring rShortcut = rCopy0(r);
-      bool overflow;
-      /**
-       * prepend extra weight vector for homogeneity
-       */
-      int* order = rShortcut->order;
-      int* block0 = rShortcut->block0;
-      int* block1 = rShortcut->block1;
-      int** wvhdl = rShortcut->wvhdl;
-      int h = rBlocks(r); int n = rVar(r);
-      rShortcut->order = (int*) omAlloc0((h+1)*sizeof(int));
-      rShortcut->block0 = (int*) omAlloc0((h+1)*sizeof(int));
-      rShortcut->block1 = (int*) omAlloc0((h+1)*sizeof(int));
-      rShortcut->wvhdl = (int**) omAlloc0((h+1)*sizeof(int*));
-      rShortcut->order[0] = ringorder_a;
-      rShortcut->block0[0] = 1;
-      rShortcut->block1[0] = n;
-      rShortcut->wvhdl[0] = ZVectorToIntStar(v,overflow);
-      for (int i=1; i<=h; i++)
+  poly checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector w) const;
+
+  /**
+   * suppose w a weight in maximal groebner cone of >
+   * suppose I (initially) reduced standard basis w.r.t. > and inI initial forms of its elements w.r.t. w
+   * suppose inJ elements of initial ideal that are homogeneous w.r.t w
+   * returns J elements of ideal whose initial form w.r.t. w are inI
+   * in particular, if w lies also inthe maximal groebner cone of another ordering >'
+   * and inJ is a standard basis of the initial ideal w.r.t. >'
+   * then the returned J will be a standard baiss of the ideal w.r.t. >'
+   */
+  ideal getWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const;
+
+  /**
+   * given generators of the initial ideal, computes its standard basis
+   */
+  ideal getStdOfInitialIdeal(const ideal inI, const ring r) const;
+
+  /**
+   * given an interior point of a groebner cone
+   * computes the groebner cone adjacent to it
+   */
+  std::pair<ideal,ring> getFlip(const ideal I, const ring r, const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
+
+  /**
+   * checks whether the first elements of I is p-t
+   */
+  static bool checkFirstGenerator(ideal I, ring r)
+  {
+    // check for correct size of ideal and first generator
+    if (idSize(I)==0)
+    {
+      WerrorS("checkFirstElement: ideal empty");
+      return false;
+    }
+    poly p = I->m[0];
+    if (size(p)!=2)
+    {
+      WerrorS("checkFirstElement: no binomial");
+      return false;
+    }
+
+    // check first term
+    int n = rVar(r);
+    if (!n_IsUnit(p_GetCoeff(p,r->cf),r->cf))
+    {
+      WerrorS("checkFirstElement: wrong coefficient");
+      return false;
+    }
+    for (int i=1; i<=n; i++)
+    {
+      if (p_GetExp(p,i,r)!=0)
       {
-        rShortcut->order[i] = order[i-1];
-        rShortcut->block0[i] = block0[i-1];
-        rShortcut->block1[i] = block1[i-1];
-        rShortcut->wvhdl[i] = wvhdl[i-1];
+        WerrorS("checkFirstElement: wrong exponent vector");
+        return false;
       }
-      rComplete(rShortcut);
-      rTest(rShortcut);
-      omFree(order);
-      omFree(block0);
-      omFree(block1);
-      omFree(wvhdl);
-
-      ideal inI = initial(I,r,w);
-      int k = idSize(inI);
-      ideal inIShortcut = idInit(k);
-      nMapFunc identity = n_SetMap(r->cf,rShortcut->cf);
-      for (int i=0; i<k; i++)
-        inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,identity,NULL,0);
-
-      poly p = checkForMonomialViaSuddenSaturation(inIShortcut,rShortcut);
-      poly monomial = NULL;
-      if (p!=NULL)
+    }
+
+    // check second term
+    pIter(p);
+    if (!n_IsOne(p_GetCoeff(p,r->cf),r->cf))
+    {
+      WerrorS("checkFirstElement: wrong coefficient");
+      return false;
+    }
+    if (p_GetExp(p,1,r)!=1)
+    {
+      WerrorS("checkFirstElement: wrong exponent vector");
+      return false;
+    }
+    for (int i=2; i<=n; i++)
+    {
+      if (p_GetExp(p,i,r)!=0)
       {
-        monomial=p_One(r);
-        for (int i=1; i<n; i++)
-          p_SetExp(monomial,i,p_GetExp(p,i,rShortcut),r);
-        p_Delete(&p,rShortcut);
+        WerrorS("checkFirstElement: wrong exponent vector");
+        return false;
       }
-      id_Delete(&inI,r);
-      id_Delete(&inIShortcut,rShortcut);
-      rDelete(rShortcut);
-      return monomial;
-    }
-    else
-    {
-      ring rShortcut = rCopy0(r);
-      bool overflow;
-      /**
-       * prepend extra weight vector for homogeneity
-       */
-      int* order = rShortcut->order;
-      int* block0 = rShortcut->block0;
-      int* block1 = rShortcut->block1;
-      int** wvhdl = rShortcut->wvhdl;
-      int h = rBlocks(r); int n = rVar(r);
-      rShortcut->order = (int*) omAlloc0((h+1)*sizeof(int));
-      rShortcut->block0 = (int*) omAlloc0((h+1)*sizeof(int));
-      rShortcut->block1 = (int*) omAlloc0((h+1)*sizeof(int));
-      rShortcut->wvhdl = (int**) omAlloc0((h+1)*sizeof(int*));
-      rShortcut->order[0] = ringorder_a;
-      rShortcut->block0[0] = 1;
-      rShortcut->block1[0] = n;
-      rShortcut->wvhdl[0] = ZVectorToIntStar(v,overflow);
-      for (int i=1; i<=h; i++)
-      {
-        rShortcut->order[i] = order[i-1];
-        rShortcut->block0[i] = block0[i-1];
-        rShortcut->block1[i] = block1[i-1];
-        rShortcut->wvhdl[i] = wvhdl[i-1];
-      }
-      omFree(order);
-      omFree(block0);
-      omFree(block1);
-      omFree(wvhdl);
-      /**
-       * change ground domain into finite field
-       */
-      nKillChar(rShortcut->cf);
-      rShortcut->cf = nCopyCoeff(shortcutRing->cf);
-      rComplete(rShortcut);
-      rTest(rShortcut);
-
-      ideal inI = initial(I,r,w);
-      int k = idSize(inI);
-      ideal inIShortcut = idInit(k);
-      nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
-      for (int i=0; i<k; i++)
-        inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
-
-      idSkipZeroes(inIShortcut);
-      poly p = checkForMonomialViaSuddenSaturation(inIShortcut,rShortcut);
-      poly monomial = NULL;
-      if (p!=NULL)
-      {
-        monomial=p_One(r);
-        for (int i=1; i<n; i++)
-          p_SetExp(monomial,i,p_GetExp(p,i,rShortcut),r);
-        p_Delete(&p,rShortcut);
-      }
-      id_Delete(&inI,r);
-      id_Delete(&inIShortcut,rShortcut);
-      rDelete(rShortcut);
-      return monomial;
-    }
-  }
-  std::pair<ideal,ring> flip(const ideal I, const ring r,
-                             const gfan::ZVector interiorPoint,
-                             const gfan::ZVector facetNormal) const
-  {
-    gfan::ZVector adjustedInteriorPoint = adjustWeightForHomogeneity(interiorPoint);
-    gfan::ZVector adjustedFacetNormal = adjustWeightUnderHomogeneity(facetNormal,adjustedInteriorPoint);
-    return flip0(I,r,interiorPoint,facetNormal,adjustedInteriorPoint,adjustedFacetNormal);
+    }
+
+    return true;
   }
 };

Singular/dyn_modules/gfanlib/witness.cc

@@ -14 +14 @@
  * Returns an error if this is not possible.
  **/
-static matrix divisionDiscardingRemainder(const poly f, const ideal G, const ring r)
+matrix divisionDiscardingRemainder(const poly f, const ideal G, const ring r)
 {
   ring origin = currRing;
@@ -20 +20 @@
   ideal F = idInit(1); F->m[0]=f;
   ideal m = idLift(G,F);
-  // ideal m = idLift(G,F,NULL,FALSE,TRUE);
   F->m[0]=NULL; id_Delete(&F, currRing);
+  matrix Q = id_Module2formatedMatrix(m,IDELEMS(G),1,currRing);
+  if (origin != r) rChangeCurrRing(origin);
+  return Q;
+}
+
+/**
+ * Given F[0],...,F[l] and G[0],...,G[k],
+ * computes Q[i,j] for i=0,..,k and j=0,...,l such that
+ * F[j] = Q[0,j]*G[0]+...+Q[k,j]*G[k].
+ */
+matrix divisionDiscardingRemainder(const ideal F, const ideal G, const ring r)
+{
+  ring origin = currRing;
+  if (origin != r) rChangeCurrRing(r);
+  ideal m = idLift(G,F);
   matrix Q = id_Module2formatedMatrix(m,IDELEMS(G),1,currRing);
   if (origin != r) rChangeCurrRing(origin);
@@ -89 +103 @@
 }
 #endif
+
+/***
+ * Let w be the uppermost weight vector in the matrix defining the ordering on r.
+ * Let I be a Groebner basis of an ideal in r, inI its initial form with respect w.
+ * Given an w-homogeneous element m of inI, computes a witness g of m in I,
+ * i.e. g in I such that in_w(g)=m.
+ **/
+ideal witness(const ideal inI, const ideal J, const ring r)
+{
+  ring origin = currRing;
+  if (origin!=r)
+    rChangeCurrRing(r);
+  ideal NFinI = kNF(J,r->qideal,inI);
+  if (origin!=r)
+    rChangeCurrRing(origin);
+
+  int k = idSize(inI);
+  ideal I = idInit(k);
+  for (int i=0; i<k; i++)
+  {
+    I->m[i] = p_Add_q(p_Copy(inI->m[i],r),p_Neg(NFinI->m[i],r),r);
+    NFinI->m[i] = NULL;
+  }
+
+  return I;
+}

Singular/dyn_modules/gfanlib/witness.h

@@ -5 +5 @@
 #include <libpolys/polys/simpleideals.h>
 
+matrix divisionDiscardingRemainder(const poly f, const ideal G, const ring r);
+matrix divisionDiscardingRemainder(const ideal F, const ideal G, const ring r);
 poly witness(const poly m, const ideal I, const ideal inI, const ring r);
+ideal witness(const ideal inI, const ideal J, const ring r);
 
 #ifndef NDEBUG

gfanlib/gfanlib_matrix.h

@@ -39 +39 @@
   {
     Matrix ret(1,v.size());
-    for(int i=0;i<v.size();i++)ret[0][i]=v[i];
+    for(unsigned i=0;i<v.size();i++)ret[0][i]=v[i];
     return ret;
   }
@@ -331 +331 @@
     int retSwaps=0;
     int currentRow=0;
-
-    if (height<2)
-      return 0;
 
     for(int i=0;i<width;i++)

gfanlib/gfanlib_vector.h

@@ -85 +85 @@
   unsigned int size()const{return v.size();};
   void resize(int n){v.resize(n,typ());};
-  void grow(int i){if(size()<i)resize(i);}
+  void grow(int i){if((int)size()<i)resize(i);}
   void push_back(typ a)
   {