Changeset b71400a in git

Timestamp:

Sep 15, 2014, 11:08:21 AM (9 years ago)

Author:

Yue Ren <ren@…>

Branches:

(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')

Children:

dffd1541cc74564f509bd5d284948852484430b0

Parents:

9c03260ac0404405e4994f73b30ef46f079535d8

git-author:

Yue Ren <ren@mathematik.uni-kl.de>2014-09-15 11:08:21+02:00

git-committer:

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

Message:

chg: status update 15.09.

Files:

• Singular/dyn_modules/gfanlib/adjustWeights.cc (modified) (3 diffs)
• Singular/dyn_modules/gfanlib/callgfanlib_conversion.cc (modified) (1 diff)
• Singular/dyn_modules/gfanlib/callgfanlib_conversion.h (modified) (1 diff)
• Singular/dyn_modules/gfanlib/groebnerCone.cc (modified) (10 diffs)
• Singular/dyn_modules/gfanlib/groebnerCone.h (modified) (1 diff)
• Singular/dyn_modules/gfanlib/initial.cc (modified) (3 diffs)
• Singular/dyn_modules/gfanlib/initial.h (modified) (1 diff)
• Singular/dyn_modules/gfanlib/ppinitialReduction.cc (modified) (7 diffs)
• Singular/dyn_modules/gfanlib/ppinitialReduction.h (modified) (1 diff)
• Singular/dyn_modules/gfanlib/singularWishlist.h (modified) (1 diff)
• Singular/dyn_modules/gfanlib/startingCone.cc (modified) (4 diffs)
• Singular/dyn_modules/gfanlib/tropicalCurves.cc (modified) (4 diffs)
• Singular/dyn_modules/gfanlib/tropicalCurves.h (modified) (1 diff)
• Singular/dyn_modules/gfanlib/tropicalStrategy.cc (modified) (15 diffs)
• Singular/dyn_modules/gfanlib/tropicalStrategy.h (modified) (6 diffs)
• Singular/dyn_modules/gfanlib/tropicalTraversal.cc (modified) (2 diffs)
• Singular/dyn_modules/gfanlib/witness.cc (modified) (1 diff)
• gfanlib/gfanlib_z.h (modified) (3 diffs)
• libpolys/polys/matpol.cc (modified) (2 diffs)
• libpolys/polys/matpol.h (modified) (2 diffs)

Singular/dyn_modules/gfanlib/adjustWeights.cc

@@ -18 +18 @@
 static bool checkForNonPositiveLaterEntries(const gfan::ZVector w)
 {
-  if (w[0].sign()<0)
-  {
-    std::cout << "ERROR: negative weight in weight vector first entry" << std::endl
-              << "weight: " << w << std::endl;
-    return false;
-  }
-  for (unsigned i=1; i<w.size(); i++)
+  // if (w[0].sign()<0)
+  // {
+  //   std::cout << "ERROR: negative weight in weight vector first entry" << std::endl
+  //             << "weight: " << w << std::endl;
+  //   return false;
+  // }
+  for (unsigned i=0; i<w.size(); i++)
   {
     if (w[i].sign()<=0)
@@ -71 +71 @@
   for (unsigned i=1; i<w.size(); i++)
     v[i]=-w[i]+max+1;
-  assume(checkForNonPositiveLaterEntries(v));
+  // assume(checkForNonPositiveLaterEntries(v));
   return v;
 }
@@ -108 +108 @@
 gfan::ZVector valued_adjustWeightUnderHomogeneity(const gfan::ZVector e, const gfan::ZVector w)
 {
-  assume(checkForNonPositiveLaterEntries(w));
+  // assume(checkForNonPositiveLaterEntries(w));
   /* find k such that e+k*w is strictly positive,
    * i.e. k such that e[i]+k*w[i] is strictly positive

Singular/dyn_modules/gfanlib/callgfanlib_conversion.cc

@@ -96 +96 @@
 }
 
+gfan::ZVector expvToZVector(const int n, const int* expv)
+{
+  gfan::ZVector zv(n);
+  for(int i=0; i<n; i++)
+    zv[i]=gfan::Integer(expv[i+1]);
+  return zv;
+}
+
 gfan::ZVector wvhdlEntryToZVector(const int n, const int* wvhdl0)
 {

Singular/dyn_modules/gfanlib/callgfanlib_conversion.h

@@ -19 +19 @@
 char* toString(gfan::ZMatrix const &m);
 
+gfan::ZVector expvToZVector(const int n, const int* expv);
+
 gfan::ZFan* toFanStar(std::set<gfan::ZCone> setOfCones);
 std::set<gfan::ZVector> rays(std::set<gfan::ZCone> setOfCones);

Singular/dyn_modules/gfanlib/groebnerCone.cc

@@ -81 +81 @@
   }
   if (r) polynomialRing = rCopy(r);
+  currentCase.pReduce(polynomialIdeal,r);
   currentCase.reduce(reducedPolynomialIdeal,polynomialRing);
 
@@ -96 +97 @@
     {
       p_GetExpV(g,leadexpv,r);
-      leadexpw = intStar2ZVector(n, leadexpv);
+      leadexpw = expvToZVector(n, leadexpv);
       pIter(g);
       while (g)
       {
         p_GetExpV(g,tailexpv,r);
-        tailexpw = intStar2ZVector(n, tailexpv);
+        tailexpw = expvToZVector(n, tailexpv);
         inequalities.appendRow(leadexpw-tailexpw);
         pIter(g);
@@ -109 +110 @@
   omFreeSize(leadexpv,(n+1)*sizeof(int));
   omFreeSize(tailexpv,(n+1)*sizeof(int));
-  if (currentStrategy->restrictToLowerHalfSpace())
-  {
-    gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
-    lowerHalfSpaceCondition[0] = -1;
-    inequalities.appendRow(lowerHalfSpaceCondition);
-  }
+  // if (currentStrategy->restrictToLowerHalfSpace())
+  // {
+  //   gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
+  //   lowerHalfSpaceCondition[0] = -1;
+  //   inequalities.appendRow(lowerHalfSpaceCondition);
+  // }
 
   polyhedralCone = gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
@@ -144 +145 @@
   }
   if (r) polynomialRing = rCopy(r);
+  currentCase.pReduce(polynomialIdeal,r);
   currentCase.reduce(reducedPolynomialIdeal,polynomialRing);
 
@@ -205 +207 @@
   }
   if (r) polynomialRing = rCopy(r);
+  currentCase.pReduce(polynomialIdeal,r);
   currentCase.reduce(reducedPolynomialIdeal,polynomialRing);
 
@@ -259 +262 @@
   assume(checkPolynomialInput(I,r));
   assume(checkPolynomialInput(inI,r));
+  currentCase.pReduce(polynomialIdeal,r);
   currentStrategy->reduce(reducedPolynomialIdeal,polynomialRing);
   // currentStrategy->reduce(initialPolynomialIdeal,polynomialRing);
@@ -435 +439 @@
    *   to obtain an initial form with respect to interiorPoint+e*facetNormal,
    *   for e>0 sufficiently small */
-  std::pair<ideal,ring> flipped = currentStrategy->flip(reducedPolynomialIdeal,polynomialRing,interiorPoint,facetNormal);
+  std::pair<ideal,ring> flipped = currentStrategy->getFlip(polynomialIdeal,reducedPolynomialIdeal,polynomialRing,interiorPoint,facetNormal);
   assume(checkPolynomialInput(flipped.first,flipped.second));
   groebnerCone flippedCone(flipped.first, flipped.second, interiorPoint, facetNormal, *currentStrategy);
@@ -547 +551 @@
   for (int i=1; i<r-1; i++)
   {
-    newInequalities = inequalities.submatrix(0,0,i-1,c);
+    newInequalities = inequalities.submatrix(0,0,i,c);
     newInequalities.append(inequalities.submatrix(i+1,0,r,c));
     newEquations = equations;
@@ -566 +570 @@
 
 
+bool groebnerCone::pointsOutwards(const gfan::ZVector w) const
+{
+  gfan::ZCone dual = polyhedralCone.dualCone();
+  return (!dual.contains(w));
+}
+
+
 /***
  * Returns a complete list of neighboring Groebner cones in the tropical variety.
@@ -578 +589 @@
     {
       ideal initialIdeal = initial(polynomialIdeal,polynomialRing,interiorPoints[i]);
-      std::set<gfan::ZVector> rays = raysOfTropicalStar(initialIdeal,polynomialRing,interiorPoints[i],*currentStrategy);
-      for (std::set<gfan::ZVector>::iterator ray = rays.begin(); ray!=rays.end(); ray++)
-        neighbours.insert(this->flipCone(interiorPoints[i],*ray));
+      gfan::ZMatrix ray = raysOfTropicalStar(initialIdeal,polynomialRing,interiorPoints[i],*currentStrategy);
+      for (int j=0; j<ray.getHeight(); j++)
+        if (pointsOutwards(ray[j]))
+        {
+          groebnerCone neighbour = flipCone(interiorPoints[i],ray[j]);
+          neighbours.insert(neighbour);
+        }
     }
   }

Singular/dyn_modules/gfanlib/groebnerCone.h

@@ -94 +94 @@
    **/
   bool checkFlipConeInput(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
+  bool pointsOutwards(const gfan::ZVector) const;
   groebnerCone flipCone(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
 

Singular/dyn_modules/gfanlib/initial.cc

@@ -36 +36 @@
   for (int i=0; i<W.getHeight(); i++)
     d[i] = wDeg(p,r,W[i]);
+  return d;
+}
+
+gfan::ZVector WDeg(const poly p, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
+{
+  gfan::ZVector d = gfan::ZVector(W.getHeight()+1);
+  d[0] = wDeg(p,r,w);
+  for (int i=0; i<W.getHeight(); i++)
+    d[i+1] = wDeg(p,r,W[i]);
   return d;
 }
@@ -166 +175 @@
 }
 
+poly initial(const poly p, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
+{
+  int n = rVar(r);
+  poly q0 = p_Head(p,r);
+  poly q1 = q0;
+  gfan::ZVector d = WDeg(p,r,w,W);
+  for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
+  {
+    gfan::ZVector e = WDeg(currentTerm,r,w,W);
+    if (d<e)
+    {
+      p_Delete(&q0,r);
+      q0 = p_Head(p,r);
+      q1 = q0;
+      d = e;
+    }
+    else
+      if (d==e)
+      {
+        pNext(q1) = p_Head(currentTerm,r);
+        pIter(q1);
+      }
+  }
+  return q0;
+}
+
+
 /***
  * Runs the above procedure over all generators of an ideal.
@@ -179 +215 @@
 }
 
+ideal initial(const ideal I, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
+{
+  int k = idSize(I); ideal inI = idInit(k);
+  for (int i=0; i<k; i++)
+    inI->m[i] = initial(I->m[i],r,w,W);
+  return inI;
+}
 
 #ifndef NDEBUG

Singular/dyn_modules/gfanlib/initial.h

@@ -40 +40 @@
 ideal initial(const ideal I, const ring r, const gfan::ZVector w);
 poly initial(const poly p, const ring r, const gfan::ZMatrix W);
-poly initial(const ideal I, const ring r, const gfan::ZMatrix W);
+ideal initial(const ideal I, const ring r, const gfan::ZMatrix W);
+poly initial(const poly p, const ring r, const gfan::ZVector w, const gfan::ZMatrix W);
+ideal initial(const ideal I, const ring r, const gfan::ZVector w, const gfan::ZMatrix W);
 
 

Singular/dyn_modules/gfanlib/ppinitialReduction.cc

@@ -8 +8 @@
 #include <map>
 #include <set>
+
+bool isOrderingLocalInT(const ring r)
+{
+  poly one = p_One(r);
+  poly t = p_One(r);
+  p_SetExp(t,1,1,r);
+  p_Setm(t,r);
+  int s = p_LmCmp(one,t,r);
+  p_Delete(&one,r);
+  p_Delete(&t,r);
+  return (s==1);
+}
 
 /***
@@ -80 +92 @@
 }
 
-
 #ifndef NDEBUG
 BOOLEAN pppReduce(leftv res, leftv args)
@@ -105 +116 @@
 }
 #endif //NDEBUG
+
+bool pReduce0(ideal &I, const number p, const ring r)
+{
+  int k = idSize(I);
+  for (int i=0; i<k; i++)
+  {
+    if (I->m[i]!=NULL)
+    {
+      number c = p_GetCoeff(I->m[i],r);
+      if (!n_Equal(p,c,r->cf))
+        if (pReduce(I->m[i],p,r))
+          return true;
+    }
+  }
+  return false;
+}
 
 
@@ -494 +521 @@
 
 
+static std::vector<int> synchronize(const ideal I, const ideal Hi)
+{
+  int k = idSize(Hi);
+  int l = idSize(I);
+  std::vector<int> synch(k);
+  int j;
+  for (int i=0; i<k; i++)
+  {
+    for (j=0; j<l; j++)
+    {
+      if (I->m[i]==Hi->m[j])
+        synch[i] = j;
+    }
+    if (j==l)
+      synch[i] = -1;
+  }
+  return synch;
+}
+
+static void synchronize(ideal I, const ideal Hi, const std::vector<int> synch)
+{
+  for (unsigned i=0; i<synch.size(); i++)
+    if (synch[i]>=0)
+      I->m[i] = Hi->m[synch[i]];
+}
+
 /**
  * reduces I initially with respect to itself.
@@ -525 +578 @@
 
   std::map<long,ideal>::iterator it=H.begin();
+  it++;
   ideal Hi = it->second;
-  assume(idSize(Hi)==1);
-  assume(pLength(Hi->m[0])==2);
-  it++;
-  Hi = it->second;
 
   /***
@@ -535 +585 @@
    *  and all lower components
    **/
+  std::vector<int> synch = synchronize(I,Hi);
   if (ppreduceInitially(Hi,p,r)) return true;
+  synchronize(I,Hi,synch);
+  id_Test(Hi,r);
+  id_Test(I,r);
 
   ideal G = idInit(n); int m=0;
@@ -577 +631 @@
     }
     m += l; IDELEMS(GG) = m; GG->m = &G->m[n-m];
+    std::vector<int> synch = synchronize(I,it->second);
     if (ppreduceInitially(it->second,p,GG,r)) return true;
+    synchronize(I,it->second,synch);
     idShallowDelete(&Hi); Hi = it->second;
   }

Singular/dyn_modules/gfanlib/ppinitialReduction.h

@@ -14 +14 @@
 #endif
 
-bool ppreduceInitially(ideal I, ring r, number p);
+bool isOrderingLocalInT(const ring r);
+bool pReduce0(ideal &I, const number p, const ring r);
+bool ppreduceInitially(ideal I, const ring r, const number p);
 BOOLEAN ppreduceInitially(leftv res, leftv args);
 

Singular/dyn_modules/gfanlib/singularWishlist.h

@@ -59 +59 @@
 inline void idShallowDelete (ideal *h)
 {
-  int j,elems;
-  if (*h == NULL)
-    return;
-  elems=j=(*h)->nrows*(*h)->ncols;
-  if (j>0)
-    omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*elems);
-  omFreeBin((ADDRESS)*h, sip_sideal_bin);
-  *h=NULL;
-}
+  if (*h != NULL)
+  {
+    int k;
+    k=(*h)->nrows*(*h)->ncols;
+    if (k>0)
+      omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*k);
+    omFreeBin((ADDRESS)*h, sip_sideal_bin);
+    *h=NULL;
+  }
+  return;
+ }
 
 #endif

Singular/dyn_modules/gfanlib/startingCone.cc

@@ -21 +21 @@
 static bool checkContainmentInTropicalVariety(const groebnerCone sigma)
 {
-  gfan::ZVector w = sigma.getInteriorPoint();
   ideal I = sigma.getPolynomialIdeal();
   ring r = sigma.getPolynomialRing();
   const tropicalStrategy* currentStrategy = sigma.getTropicalStrategy();
-  poly s = currentStrategy->checkInitialIdealForMonomial(I,r,w);
-  if (s)
-  {
-    p_Delete(&s,r);
-    return false;
-  }
+
+  gfan::ZCone zc = sigma.getPolyhedralCone();
+  gfan::ZMatrix zm = zc.extremeRays();
+  for (int i=0; i<zm.getHeight(); i++)
+  {
+    gfan::ZVector w = zm[i];
+    if (currentStrategy->isValuationNonTrivial() && w[0].sign()==0)
+      continue;
+    poly s = currentStrategy->checkInitialIdealForMonomial(I,r,w);
+    if (s)
+    {
+      p_Delete(&s,r);
+      return false;
+    }
+  }
+
+  zm = zc.generatorsOfLinealitySpace();
+  for (int i=0; i<zm.getHeight(); i++)
+  {
+    gfan::ZVector w = zm[i];
+    if (currentStrategy->isValuationNonTrivial() && w[0].sign()==0)
+      continue;
+    poly s = currentStrategy->checkInitialIdealForMonomial(I,r,w);
+    if (s)
+    {
+      p_Delete(&s,r);
+      return false;
+    }
+  }
+
   return true;
 }
@@ -313 +336 @@
 {
   ring r = currentStrategy.getStartingRing();
-  ideal I = currentStrategy.getStartingIdeal();
+  const ideal I = currentStrategy.getStartingIdeal();
   if (currentStrategy.isConstantCoefficientCase())
   {
@@ -359 +382 @@
     id_Delete(&J,s);
 
-    assume(checkContainmentInTropicalVariety(startingCone));
+    // assume(checkContainmentInTropicalVariety(startingCone));
     return startingCone;
   }
@@ -463 +486 @@
     rDelete(s);
 
+    assume(checkContainmentInTropicalVariety(startingCone));
     return startingCone;
   }

Singular/dyn_modules/gfanlib/tropicalCurves.cc

@@ -142 +142 @@
 
     inIs = gfanlib_kStd_wrapper(inIs,s,isHomog);
-    ideal ininIs = initial(inIs,s,W[W.getHeight()-1]);
-
-    poly mons = checkForMonomialViaSuddenSaturation(inIs,s);
+    ideal ininIs = initial(inIs,s,w,W);
+
+    poly mons = checkForMonomialViaSuddenSaturation(ininIs,s);
     if (mons)
     {
@@ -165 +165 @@
 
 
-std::set<gfan::ZVector> raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector u, const tropicalStrategy& currentStrategy)
+gfan::ZMatrix raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector u, const tropicalStrategy& currentStrategy)
 {
   std::set<gfan::ZCone> C = tropicalStar(I,r,u,currentStrategy);
-  std::set<gfan::ZVector> raysOfC;
+  gfan::ZMatrix raysOfC(0,u.size());
   if (!currentStrategy.restrictToLowerHalfSpace())
   {
@@ -174 +174 @@
     {
       assume(zc->dimensionOfLinealitySpace()+1 == zc->dimension());
-      raysOfC.insert(zc->semiGroupGeneratorOfRay());
+      raysOfC.appendRow(zc->semiGroupGeneratorOfRay());
     }
   }
@@ -182 +182 @@
     {
       assume(zc->dimensionOfLinealitySpace()+2 == zc->dimension());
-      raysOfC.insert(zc->getRelativeInteriorPoint());
+      raysOfC.appendRow(zc->getRelativeInteriorPoint());
     }
   }

Singular/dyn_modules/gfanlib/tropicalCurves.h

@@ -9 +9 @@
 std::set<gfan::ZCone> tropicalStar(const ideal I, const ring r, const gfan::ZVector u,
                                    const tropicalStrategy& currentStrategy);
-std::set<gfan::ZVector> raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector u,
-                                           const tropicalStrategy& currentStrategy);
+gfan::ZMatrix raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector u,
+                                 const tropicalStrategy& currentStrategy);
 
 #ifndef NDEBUG

Singular/dyn_modules/gfanlib/tropicalStrategy.cc

@@ -195 +195 @@
   startingRing(rCopy(currentStrategy.getStartingRing())),
   startingIdeal(id_Copy(currentStrategy.getStartingIdeal(),currentStrategy.getStartingRing())),
-  uniformizingParameter(n_Copy(currentStrategy.getUniformizingParameter(),startingRing->cf)),
-  shortcutRing(rCopy(currentStrategy.getShortcutRing())),
+  uniformizingParameter(NULL),
+  shortcutRing(NULL),
   onlyLowerHalfSpace(currentStrategy.restrictToLowerHalfSpace()),
   weightAdjustingAlgorithm1(currentStrategy.weightAdjustingAlgorithm1),
   weightAdjustingAlgorithm2(currentStrategy.weightAdjustingAlgorithm2),
   extraReductionAlgorithm(currentStrategy.extraReductionAlgorithm)
+{
+  if (originalRing) rTest(originalRing);
+  if (originalIdeal) id_Test(originalIdeal,originalRing);
+  if (startingRing) rTest(startingRing);
+  if (startingIdeal) id_Test(startingIdeal,startingRing);
+  if (currentStrategy.getUniformizingParameter())
+  {
+    uniformizingParameter = n_Copy(currentStrategy.getUniformizingParameter(),startingRing->cf);
+    n_Test(uniformizingParameter,startingRing->cf);
+  }
+  if (currentStrategy.getShortcutRing())
+  {
+    shortcutRing = rCopy(currentStrategy.getShortcutRing());
+    rTest(shortcutRing);
+  }
+}
+
+tropicalStrategy::~tropicalStrategy()
 {
   if (originalRing) rTest(originalRing);
@@ -208 +226 @@
   if (uniformizingParameter) n_Test(uniformizingParameter,startingRing->cf);
   if (shortcutRing) rTest(shortcutRing);
-}
-
-tropicalStrategy::~tropicalStrategy()
-{
-  if (originalRing) rTest(originalRing);
-  if (originalIdeal) id_Test(originalIdeal,originalRing);
-  if (startingRing) rTest(startingRing);
-  if (startingIdeal) id_Test(startingIdeal,startingRing);
-  if (uniformizingParameter) n_Test(uniformizingParameter,startingRing->cf);
-  if (shortcutRing) rTest(shortcutRing);
 
   id_Delete(&originalIdeal,originalRing);
   rDelete(originalRing);
-  id_Delete(&startingIdeal,startingRing);
-  n_Delete(&uniformizingParameter,startingRing->cf);
-  rDelete(startingRing);
-  rDelete(shortcutRing);
+  if (startingIdeal) id_Delete(&startingIdeal,startingRing);
+  if (uniformizingParameter) n_Delete(&uniformizingParameter,startingRing->cf);
+  if (startingRing) rDelete(startingRing);
+  if (shortcutRing) rDelete(shortcutRing);
 }
 
@@ -251 +259 @@
 }
 
+void tropicalStrategy::putUniformizingBinomialInFront(ideal I, const ring r, const number q) const
+{
+  poly p = p_One(r);
+  p_SetCoeff(p,q,r);
+  poly t = p_One(r);
+  p_SetExp(t,1,1,r);
+  p_Setm(t,r);
+  poly pt = p_Add_q(p,p_Neg(t,r),r);
+
+  int k = idSize(I);
+  int l;
+  for (l=0; l<k; l++)
+  {
+    if (p_EqualPolys(I->m[l],pt,r))
+      break;
+  }
+  p_Delete(&pt,r);
+
+  if (l>1)
+  {
+    pt = I->m[l];
+    for (int i=l; i>0; i--)
+      I->m[l] = I->m[l-1];
+    I->m[0] = pt;
+    pt = NULL;
+  }
+  return;
+}
+
 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);
+
+  if (isValuationTrivial())
+    return false;
+
+  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
+  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
+  bool b = extraReductionAlgorithm(I,r,p);
+  // putUniformizingBinomialInFront(I,r,p);
   n_Delete(&p,r->cf);
+
+  return b;
+}
+
+bool tropicalStrategy::pReduce(ideal I, const ring r) const
+{
+  rTest(r);
+  id_Test(I,r);
+
+  if (isValuationTrivial())
+    return false;
+
+  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
+  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
+  bool b = pReduce0(I,p,r);
+  n_Delete(&p,r->cf);
+
   return b;
 }
@@ -293 +352 @@
 
   // if valuation non-trivial, change coefficient ring to residue field
-  if (valuationIsNonTrivial())
+  if (isValuationNonTrivial())
   {
     nKillChar(rShortcut->cf);
@@ -330 +389 @@
   // if it exsists, return a copy of the monomial in the input ring
   poly p = checkForMonomialViaSuddenSaturation(inIShortcut,rShortcut);
-  poly monomial = NULL;
+  poly monomial = NULL; int n = rVar(r);
   if (p!=NULL)
   {
@@ -354 +413 @@
 }
 
-ideal getWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
+ideal tropicalStrategy::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,
@@ -365 +424 @@
   else
   {
-    checkFirstGenerator(I,r);
     assume(idSize(inI)==idSize(I));
+    int uni = findPositionOfUniformizingBinomial(I,r);
+    assume(uni>=0);
     /**
      * change ground ring into finite field
@@ -373 +433 @@
     ring rShortcut = copyAndChangeCoefficientRing(r);
 
-    int k = idSize(inI);
+    int k = idSize(inJ);
     int l = idSize(I);
     ideal inJShortcut = idInit(k);
@@ -390 +450 @@
      */
     matrix QShortcut = divisionDiscardingRemainder(inJShortcut,inIShortcut,rShortcut);
-    matrix Q = mpNew(k,l);
+    matrix Q = mpNew(l,k);
     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);
+      Q->m[ij] = p_PermPoly(QShortcut->m[ij],NULL,rShortcut,r,takingRepresentatives,NULL,0);
+
+    nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
+    number p = identity(uniformizingParameter,startingRing->cf,r->cf);
 
     /**
@@ -399 +462 @@
      */
     ideal J = idInit(k);
-    for (int j=0; j<l; j++)
+    for (int j=0; j<k; 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);
+      for (int i=0; i<l; i++)
+      {
+        poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
+        poly inIi = p_Copy(inI->m[i],r);
+        q0 = p_Add_q(q0,p_Neg(p_Mult_q(qij,inIi,r),r),r);
+      }
+      q0 = p_Div_nn(q0,p,r);
+      poly q0g0 = p_Mult_q(q0,p_Copy(I->m[uni],r),r);
+      // q0 = NULL;
+      poly qigi = NULL;
+      for (int i=0; i<l; i++)
+      {
+        poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
+        // poly inIi = p_Copy(I->m[i],r);
+        poly Ii = p_Copy(I->m[i],r);
+        qigi = p_Add_q(qigi,p_Mult_q(qij,Ii,r),r);
+      }
+      J->m[j] = p_Add_q(q0g0,qigi,r);
     }
 
+    id_Delete(&inIShortcut,rShortcut);
+    id_Delete(&inJShortcut,rShortcut);
+    mp_Delete(&QShortcut,rShortcut);
+    rDelete(rShortcut);
+    mp_Delete(&Q,r);
+    n_Delete(&p,r->cf);
     return J;
   }
@@ -437 +518 @@
     inJ = idInit(k+1);
     inJ->m[0] = p_One(r);
-    p_SetCoeff(inJ->m[0],nMap(uniformizingParameter,startingRing->cf,r->cf),r);
+    nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
+    p_SetCoeff(inJ->m[0],identity(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);
+
+    id_Delete(&inIShortcut,rShortcut);
+    rDelete(rShortcut);
   }
   return inJ;
 }
 
-ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector w, const gfan::ZVector v)
+static void deleteOrdering(ring r)
+{
+  if (r->order != NULL)
+  {
+    int i=rBlocks(r);
+    assume(r->block0 != NULL && r->block1 != NULL && r->wvhdl != NULL);
+    /* delete order */
+    omFreeSize((ADDRESS)r->order,i*sizeof(int));
+    omFreeSize((ADDRESS)r->block0,i*sizeof(int));
+    omFreeSize((ADDRESS)r->block1,i*sizeof(int));
+    /* delete weights */
+    for (int j=0; j<i; j++)
+      if (r->wvhdl[j]!=NULL)
+        omFree(r->wvhdl[j]);
+    omFreeSize((ADDRESS)r->wvhdl,i*sizeof(int *));
+  }
+  else
+    assume(r->block0 == NULL && r->block1 == NULL && r->wvhdl == NULL);
+  return;
+}
+
+ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector w, const gfan::ZVector v) const
 {
   // copy shortcutRing and change to desired ordering
@@ -452 +558 @@
   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);
+  gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(v,wAdjusted);
   s->order = (int*) omAlloc0(5*sizeof(int));
   s->block0 = (int*) omAlloc0(5*sizeof(int));
@@ -488 +566 @@
   s->block0[0] = 1;
   s->block1[0] = n;
-  s->wvhdl[0] = ZVectorToIntStar(adjustedInteriorPoint,ok);
+  s->wvhdl[0] = ZVectorToIntStar(wAdjusted,ok);
   s->order[1] = ringorder_a;
   s->block0[1] = 1;
   s->block1[1] = n;
-  s->wvhdl[1] = ZVectorToIntStar(adjustedFacetNormal,ok);
+  s->wvhdl[1] = ZVectorToIntStar(vAdjusted,ok);
   s->order[2] = ringorder_dp;
   s->block0[2] = 1;
@@ -503 +581 @@
 }
 
-
-std::pair<ideal,ring> tropicalStrategy::getFlip(const ideal I, const ring r,
+ring tropicalStrategy::copyAndChangeOrderingLS(const ring r, const gfan::ZVector w, const gfan::ZVector v) const
+{
+  // copy shortcutRing and change to desired ordering
+  bool ok;
+  ring s = rCopy0(r);
+  int n = rVar(s);
+  deleteOrdering(s);
+  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(w,ok);
+  s->order[1] = ringorder_a;
+  s->block0[1] = 1;
+  s->block1[1] = n;
+  s->wvhdl[1] = ZVectorToIntStar(v,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 ideal redI, const ring r,
                                                 const gfan::ZVector interiorPoint,
                                                 const gfan::ZVector facetNormal) const
 {
   assume(isValuationTrivial() || interiorPoint[0].sign()<0);
-  assume(checkFirstGenerators(I,r));
+  assume(checkForUniformizingBinomial(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);
+  ideal inIr = initial(redI,r,interiorPoint);
   ring sAdjusted = copyAndChangeOrderingWP(r,interiorPoint,facetNormal);
   nMapFunc identity = n_SetMap(r->cf,sAdjusted->cf);
-  int k = idSize(I); ideal inIs = idInit(k);
+  int k = idSize(I); ideal inIsAdjusted = 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);
+  ideal inJsAdjusted = getStdOfInitialIdeal(inIsAdjusted,sAdjusted);
 
   // find witnesses of the new standard basis elements of the initial ideal
@@ -529 +636 @@
     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);
+  ideal Jr = getWitness(inJr,inIr,redI,r);
+  ring s = copyAndChangeOrderingLS(r,interiorPoint,facetNormal);
   identity = n_SetMap(r->cf,s->cf);
+  ideal Js = idInit(k);
   for (int i=0; i<k; i++)
-    Js->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,r,s,identity,NULL,0);
-
+    Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identity,NULL,0);
+
+  // this->reduce(Jr,r);
+  // cleanup
+  id_Delete(&inJsAdjusted,sAdjusted);
+  rDelete(sAdjusted);
+  id_Delete(&Jr,r);
+  id_Delete(&inJr,r);
+
+  assume(isValuationTrivial() || isOrderingLocalInT(s));
   return std::make_pair(Js,s);
 }
+
+
+bool tropicalStrategy::checkForUniformizingBinomial(const ideal I, const ring r) const
+{
+  // if the valuation is trivial,
+  // then there is no special condition the first generator has to fullfill
+  if (isValuationTrivial())
+    return true;
+
+  // if the valuation is non-trivial then checks if the first generator is p-t
+  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
+  poly p = p_One(r);
+  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
+  poly t = p_One(r);
+  p_SetExp(t,1,1,r);
+  p_Setm(t,r);
+  poly pt = p_Add_q(p,p_Neg(t,r),r);
+
+  for (int i=0; i<idSize(I); i++)
+  {
+    if (p_EqualPolys(I->m[i],pt,r))
+    {
+      p_Delete(&pt,r);
+      return true;
+    }
+  }
+  p_Delete(&pt,r);
+  return false;
+}
+
+int tropicalStrategy::findPositionOfUniformizingBinomial(const ideal I, const ring r) const
+{
+  assume(isValuationNonTrivial());
+
+  // if the valuation is non-trivial then checks if the first generator is p-t
+  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
+  poly p = p_One(r);
+  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
+  poly t = p_One(r);
+  p_SetExp(t,1,1,r);
+  p_Setm(t,r);
+  poly pt = p_Add_q(p,p_Neg(t,r),r);
+
+  for (int i=0; i<idSize(I); i++)
+  {
+    if (p_EqualPolys(I->m[i],pt,r))
+    {
+      p_Delete(&pt,r);
+      return i;
+    }
+  }
+  p_Delete(&pt,r);
+  return -1;
+}
+
+bool tropicalStrategy::checkForUniformizingParameter(const ideal inI, const ring r) const
+{
+  // if the valuation is trivial,
+  // then there is no special condition the first generator has to fullfill
+  if (isValuationTrivial())
+    return true;
+
+  // if the valuation is non-trivial then checks if the first generator is p
+  if (inI->m[0]==NULL)
+    return false;
+  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
+  poly p = p_One(r);
+  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
+
+  for (int i=0; i<idSize(inI); i++)
+  {
+    if (p_EqualPolys(inI->m[i],p,r))
+    {
+      p_Delete(&p,r);
+      return true;
+    }
+  }
+  p_Delete(&p,r);
+  return false;
+}

Singular/dyn_modules/gfanlib/tropicalStrategy.h

@@ -91 +91 @@
   bool (*extraReductionAlgorithm) (ideal I, ring r, number p);
 
+  ring copyAndChangeCoefficientRing(const ring r) const;
+  ring copyAndChangeOrderingWP(const ring r, const gfan::ZVector w, const gfan::ZVector v) const;
+  ring copyAndChangeOrderingLS(const ring r, const gfan::ZVector w, const gfan::ZVector v) const;
+
+  /**
+   * if valuation non-trivial, checks whether the generating system contains p-t
+   * otherwise returns true
+   */
+  bool checkForUniformizingBinomial(const ideal I, const ring r) const;
+
+  /**
+   * if valuation non-trivial, checks whether the genearting system contains p
+   * otherwise returns true
+   */
+  bool checkForUniformizingParameter(const ideal inI, const ring r) const;
+  int findPositionOfUniformizingBinomial(const ideal I, const ring r) const;
+  void putUniformizingBinomialInFront(ideal I, const ring r, const number q) const;
+
+
 public:
 
@@ -154 +173 @@
   ring getStartingRing() const
   {
-    rTest(startingRing);
+    if (startingRing) rTest(startingRing);
     return startingRing;
   }
@@ -186 +205 @@
   ring getShortcutRing() const
   {
-    rTest(shortcutRing);
+    if (shortcutRing) rTest(shortcutRing);
     return shortcutRing;
   }
@@ -258 +277 @@
   void changeCoefficientToResidueField(ring r) const;
 
-
-
   /**
    * reduces the generators of an ideal I so that
@@ -265 +282 @@
    */
   bool reduce(ideal I, const ring r) const;
+
+  bool pReduce(ideal I, const ring r) const;
 
   /**
@@ -292 +311 @@
    * 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)
-      {
-        WerrorS("checkFirstElement: wrong exponent vector");
-        return false;
-      }
-    }
-
-    // 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)
-      {
-        WerrorS("checkFirstElement: wrong exponent vector");
-        return false;
-      }
-    }
-
-    return true;
-  }
+  std::pair<ideal,ring> getFlip(const ideal I, const ideal redI, const ring r, const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
 };
 

Singular/dyn_modules/gfanlib/tropicalTraversal.cc

@@ -8 +8 @@
   while(!workingList.empty())
   {
-    std::cout << "starting traversal" << std::endl;
+    // std::cout << "starting traversal" << std::endl;
     const groebnerCone sigma=*(workingList.begin());
     const groebnerCones neighbours = sigma.tropicalNeighbours();
@@ -18 +18 @@
     tropicalVariety.insert(sigma);
     workingList.erase(sigma);
-    std::cout << "tropicalVariety.size():" << tropicalVariety.size() << std::endl;
-    std::cout << "workingList.size():" << workingList.size() << std::endl;
+    // std::cout << "tropicalVariety.size():" << tropicalVariety.size() << std::endl;
+    // std::cout << "workingList.size():" << workingList.size() << std::endl;
   }
   return tropicalVariety;

Singular/dyn_modules/gfanlib/witness.cc

@@ -35 +35 @@
   ring origin = currRing;
   if (origin != r) rChangeCurrRing(r);
-  ideal m = idLift(G,F);
-  matrix Q = id_Module2formatedMatrix(m,IDELEMS(G),1,currRing);
+  ideal R; matrix U;
+  ideal m = idLift(G,F,&R,FALSE,FALSE,TRUE,&U);
+  matrix Q = id_Module2formatedMatrix(m,IDELEMS(G),IDELEMS(F),currRing);
+  id_Delete(&R,r);
+  mp_Delete(&U,r);
   if (origin != r) rChangeCurrRing(origin);
   return Q;

gfanlib/gfanlib_z.h

@@ -10 +10 @@
 
 #include <string.h>
+#include <iostream>
 #include <ostream>
 
@@ -68 +69 @@
     freefunc(str,strlen(str)+1);
     return f;
+  }
+  void debugPrint() const
+  {
+    void (*freefunc)(void *, size_t);
+    mp_get_memory_functions(0,0,&freefunc);
+    char *str=mpz_get_str(0,10,value);
+    std::cout << str;
+    freefunc(str,strlen(str)+1);
+    return;
   }
   Integer& operator+=(const Integer& a)
@@ -330 +340 @@
       return f;
     }
+    friend void debugPrint(IntegerTemplate const &a)
+    {
+      std::cout << a << std::endl;
+      return;
+    }
     LimbWord signExtension(LimbWord a)
     {

libpolys/polys/matpol.cc

@@ -8 +8 @@
 
 #include <stdio.h>
+#include <string.h>
+#include <sstream>
 #include <math.h>
 
@@ -756 +758 @@
   }
 }
+
+void mp_Write(matrix im, const ring r)
+{
+  int i,ii = MATROWS(im)-1;
+  int j,jj = MATCOLS(im)-1;
+  poly *pp = im->m;
+  for (i=0; i<=ii; i++)
+  {
+    for (j=0; j<=jj; j++)
+    {
+      Print("matrix[%u,%u]=",i+1,j+1);
+      if ((i<ii)||(j<jj)) p_Write(*pp++, r);
+      else                p_Write0(*pp, r);
+    }
+  }
+  return;
+}
+
 
 char * iiStringMatrix(matrix im, int dim, const ring r, char ch)

libpolys/polys/matpol.h

@@ -10 +10 @@
 // #include <kernel/structs.h>
 #include <polys/monomials/ring.h>
+#include <string.h>
 
 // THIS IS REALLY DIRTY: ip_smatrix HAS TO BE IDENTICAL TO ip_sideal
@@ -91 +92 @@
 char * iiStringMatrix(matrix im, int dim, const ring r, char ch=',');
 
+void mp_Write(matrix im, const ring r);
+
 extern omBin ip_smatrix_bin;