Changeset 3c0aa5 in git


Ignore:
Timestamp:
Jul 9, 2014, 10:19:05 AM (10 years ago)
Author:
Yue Ren <ren@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'c1c65551e4b3a0de0b13ddf30446f6e9569681b5')
Children:
e744d9067cafc8137765bab643ab9ebf8120f339
Parents:
7aa26b206f76709abff0bab4adaac2dfc4ddb4ef
git-author:
Yue Ren <ren@mathematik.uni-kl.de>2014-07-09 10:19:05+02:00
git-committer:
Yue Ren <ren@mathematik.uni-kl.de>2015-02-06 13:47:03+01:00
Message:
status update for natalia
Files:
6 added
31 edited

Legend:

Unmodified
Added
Removed
  • Singular/dyn_modules/gfanlib/Makefile.am

    r7aa26b2 r3c0aa5  
    11ACLOCAL_AMFLAGS = -I ../../m4
    22
    3 SOURCES = gfanlib_exceptions.h callgfanlib_conversion.cc callgfanlib_conversion.h bbcone.cc bbcone.h bbfan.cc bbfan.h bbpolytope.cc bbpolytope.h gfan.h gitfan.cc gitfan.h tropicalVarietyOfPolynomials.h tropicalVarietyOfPolynomials.cc ppinitialReduction.cc ppinitialReduction.h ttinitialReduction.cc ttinitialReduction.h containsMonomial.cc containsMonomial.h adjustWeights.cc adjustWeights.h tropicalStrategy.cc tropicalStrategy.h initial.cc initial.h witness.cc witness.h tropicalCurves.cc tropicalCurves.h groebnerCone.cc groebnerCone.h neighbours.cc neighbours.h startingCone.cc startingCone.h tropicalVarietyOfIdeals.cc tropicalVarietyOfIdeals.h tropicalVariety.cc tropicalVariety.h tropical.cc tropical.h gfanlib.cc
    4 
     3SOURCES = gfanlib_exceptions.h callgfanlib_conversion.cc callgfanlib_conversion.h bbcone.cc bbcone.h bbfan.cc bbfan.h bbpolytope.cc bbpolytope.h gfan.h gitfan.cc gitfan.h tropicalVarietyOfPolynomials.h tropicalVarietyOfPolynomials.cc ppinitialReduction.cc ppinitialReduction.h ttinitialReduction.cc ttinitialReduction.h containsMonomial.cc containsMonomial.h adjustWeights.cc adjustWeights.h tropicalStrategy.cc tropicalStrategy.h initial.cc initial.h witness.cc witness.h lift.cc lift.h flip.cc flip.h tropicalCurves.cc tropicalCurves.h groebnerCone.cc groebnerCone.h startingCone.cc startingCone.h tropicalTraversal.cc tropicalTraversal.h tropicalVarietyOfIdeals.cc tropicalVarietyOfIdeals.h tropicalVariety.cc tropicalVariety.h tropical.cc tropical.h gfanlib.cc
    54
    65MY_CPPFLAGS =  -I${top_srcdir} -I${top_builddir} \
  • Singular/dyn_modules/gfanlib/adjustWeights.cc

    r7aa26b2 r3c0aa5  
    11#include <gfanlib/gfanlib_vector.h>
     2#include <kernel/mod2.h>
     3
     4static bool checkForNonPositiveEntries(const gfan::ZVector w)
     5{
     6  for (unsigned i=1; i<w.size(); i++)
     7  {
     8    if (w[i].sign()<=0)
     9    {
     10      std::cout << "ERROR: non-positive weight in weight vector" << std::endl
     11                << "weight: " << w << std::endl;
     12      return false;
     13    }
     14  }
     15  return true;
     16}
    217
    318/***
     
    823gfan::ZVector nonvalued_adjustWeightForHomogeneity(const gfan::ZVector w)
    924{
    10   gfan::ZVector v=gfan::ZVector(w.size());
    11   gfan::Integer min=-1;
     25  /* find the smallest entry min of w */
     26  gfan::Integer min=w[0];
    1227  for (unsigned i=1; i<w.size(); i++)
    1328    if (w[i]<min) min=w[i];
     29  /* compute w+(1-min)*(1,...,1) and return it */
     30  gfan::ZVector v=gfan::ZVector(w.size());
    1431  if (min<1)
    1532  {
    16     for (unsigned i=1; i<w.size(); i++)
     33    for (unsigned i=0; i<w.size(); i++)
    1734      v[i]=w[i]-min+1;
    1835  }
     36  assume(checkForNonPositiveEntries(v));
    1937  return v;
    2038}
     
    2240gfan::ZVector valued_adjustWeightForHomogeneity(const gfan::ZVector w)
    2341{
     42  /* find the biggest entry max of w
     43   * amongst the later entries w[1],...,w[n] */
     44  gfan::Integer max=w[1];
     45  for (unsigned i=2; i<w.size(); i++)
     46    if (max<w[i]) max=w[i];
     47  /* compute -w+(1+max)*(0,1,...,1) and return it */
    2448  gfan::ZVector v=gfan::ZVector(w.size());
    2549  v[0]=-w[0];
    26   gfan::Integer max=-1;
    27   for (unsigned i=1; i<w.size(); i++)
    28     if (max<w[i]) max=w[i];
    2950  for (unsigned i=1; i<w.size(); i++)
    3051    v[i]=-w[i]+max+1;
     52  assume(checkForNonPositiveEntries(v));
    3153  return v;
    3254}
     
    3961gfan::ZVector nonvalued_adjustWeightUnderHomogeneity(const gfan::ZVector e, const gfan::ZVector w)
    4062{
    41   gfan::ZVector v=gfan::ZVector(e.size());
    42   gfan::Integer min=-1;
     63  assume(checkForNonPositiveEntries(w));
     64  /* find k such that e+k*w is strictly positive,
     65   * i.e. k such that e[i]+k*w[i] is strictly positive
     66   * for all i=0,...,n */
     67  gfan::Integer k((long)0);
     68  if (e[0].sign()<=0)
     69    k = gfan::Integer((long)1)-(gfan::Integer((long)e[0].sign())*e[0])/w[0];
    4370  for (unsigned i=1; i<e.size(); i++)
    44     if (e[i]<min) min=w[i];
    45   if (min<1)
    4671  {
    47     for (unsigned i=1; i<e.size(); i++)
    48       v[i]=e[i]-min+1;
     72    gfan::Integer kk = gfan::Integer((long)1)-(gfan::Integer((long)e[i].sign())*e[i])/w[i];
     73    if (k<kk)
     74      k = kk;
    4975  }
     76  /* compute e+k*w, check it for correctness and return it */
     77  gfan::ZVector v = e+k*w;
     78  assume(checkForNonPositiveEntries(v));
    5079  return v;
    5180}
     
    5382gfan::ZVector valued_adjustWeightUnderHomogeneity(const gfan::ZVector e, const gfan::ZVector w)
    5483{
    55   gfan::ZVector v=e-(e[0]/w[0]+1)*w;
    56   gfan::Integer min=1;
    57   for (unsigned i=1; i<v.size(); i++)
    58     if (v[i]<min) min=v[i];
    59   for (unsigned i=1; i<v.size(); i++)
    60     v[i]=v[i]-min+1;
     84  assume(checkForNonPositiveEntries(w));
     85  /* find k such that e+k*w is strictly positive,
     86   * i.e. k such that e[i]+k*w[i] is strictly positive
     87   * for all i=0,...,n */
     88  gfan::Integer k((long)0);
     89  if (e[0].sign()<=0)
     90    k = gfan::Integer((long)1)-(e[0]/w[0]);
     91  for (unsigned i=1; i<e.size(); i++)
     92  {
     93    gfan::Integer kk = gfan::Integer((long)1)-(gfan::Integer((long)e[i].sign())*e[i])/w[i];
     94    if (k<kk)
     95      k = kk;
     96  }
     97  /* compute e+k*w, check it for correctness and return it */
     98  gfan::ZVector v = e+k*w;
     99  assume(checkForNonPositiveEntries(v));
    61100  return v;
    62101}
  • Singular/dyn_modules/gfanlib/callgfanlib_conversion.cc

    r7aa26b2 r3c0aa5  
    9696}
    9797
     98gfan::ZVector wvhdlEntryToZVector(const int n, const int* wvhdl0)
     99{
     100  gfan::ZVector zv(n);
     101  for(int j=0; j<n; j++)
     102    zv[j]=wvhdl0[j];
     103  return zv;
     104}
     105
    98106int* ZVectorToIntStar(const gfan::ZVector &v, bool &overflow)
    99107{
  • Singular/dyn_modules/gfanlib/callgfanlib_conversion.h

    r7aa26b2 r3c0aa5  
    1515
    1616gfan::ZVector intStar2ZVector(const int d, const int* i);
     17gfan::ZVector wvhdlEntryToZVector(const int n, const int* wvhdl0);
    1718int* ZVectorToIntStar(const gfan::ZVector &v, bool &overflow);
    1819char* toString(gfan::ZMatrix const &m);
  • Singular/dyn_modules/gfanlib/groebnerCone.cc

    r7aa26b2 r3c0aa5  
     1#include <utility>
     2
    13#include <kernel/kstd1.h>
    24#include <kernel/ideals.h>
     
    810
    911#include <gfanlib/gfanlib.h>
    10 
     12#include <gfanlib/gfanlib_matrix.h>
     13
     14#include <tropicalStrategy.h>
     15#include <groebnerCone.h>
    1116#include <callgfanlib_conversion.h>
    12 #include <groebnerCone.h>
     17#include <containsMonomial.h>
    1318#include <initial.h>
    14 
    15 /***
    16  * Computes the Groebner cone of a polynomial g in ring r containing w relatively.
    17  * Assumes that r has a weighted ordering with weight in the said Groebner cone.
    18  **/
    19 gfan::ZCone sloppyGroebnerCone(const poly g, const ring r, const gfan::ZVector w)
    20 {
    21   int n = r->N;
    22   gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
    23   gfan::ZMatrix equations = gfan::ZMatrix(0,n);
    24 
    25   int* expv = (int*) omAlloc((n+1)*sizeof(int));
    26   p_GetExpV(g,expv,r);
    27   gfan::ZVector leadexp = intStar2ZVector(n,expv);
    28   long d = wDeg(g,r,w);
    29 
    30   poly h=g->next;
    31   for (; h && wDeg(h,r,w)==d; pIter(h))
    32   {
    33     p_GetExpV(h,expv,r);
    34     equations.appendRow(leadexp-intStar2ZVector(n,expv));
    35   }
    36 
    37   for (; h; pIter(h))
    38   {
    39     p_GetExpV(h,expv,r);
    40     inequalities.appendRow(leadexp-intStar2ZVector(n,expv));
    41   }
    42 
    43   omFreeSize(expv,(n+1)*sizeof(int));
    44   return gfan::ZCone(inequalities,equations);
    45 }
    46 
    47 /***
    48  * Computes the Groebner cone of an ideal I in ring r containing w relatively.
    49  * Assumes that r has a weighted ordering with weight in the said Groebner cone.
    50  **/
    51 gfan::ZCone sloppyGroebnerCone(const ideal I, const ring r, const gfan::ZVector w)
    52 {
    53   int k = idSize(I);
    54   gfan::ZCone zc = gfan::ZCone(r->N);
    55   for (int i=0; i<k; i++)
    56     zc = intersection(zc,sloppyGroebnerCone(I->m[i],r,w));
    57   return zc;
    58 }
    59 
    60 /***
    61  * Computes the Groebner cone of a polynomial g in ring r containing w relatively.
    62  **/
    63 gfan::ZCone groebnerCone(const poly g, const ring r, const gfan::ZVector w)
    64 {
    65   int n = r->N;
    66   gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
    67   gfan::ZMatrix equations = gfan::ZMatrix(0,n);
    68 
    69   int* expv = (int*) omAlloc((n+1)*sizeof(int));
    70   p_GetExpV(g,expv,r);
    71   gfan::ZVector leadexp = intStar2ZVector(n,expv);
    72   long d = wDeg(g,r,w);
    73 
    74   for (poly h=g->next; h; pIter(h))
    75   {
    76     p_GetExpV(h,expv,r);
    77     if (wDeg(h,r,w)<d)
    78       inequalities.appendRow(leadexp-intStar2ZVector(n,expv));
    79     else
    80       equations.appendRow(leadexp-intStar2ZVector(n,expv));
    81   }
    82 
    83   omFreeSize(expv,(n+1)*sizeof(int));
    84   return gfan::ZCone(inequalities,equations);
    85 }
    86 
    87 /***
    88  * Computes the Groebner cone of an ideal I in ring r containing w relatively.
    89  * Assumes that r has a weighted ordering with weight in the said Groebner cone.
    90  **/
    91 gfan::ZCone groebnerCone(const ideal I, const ring r, const gfan::ZVector w)
    92 {
    93   int k = idSize(I);
    94   gfan::ZCone zc = gfan::ZCone(r->N);
    95   for (int i=0; i<k; i++)
    96     zc = intersection(zc,groebnerCone(I->m[i],r,w));
    97   return zc;
    98 }
    99 
    100 gfan::ZCone fullGroebnerCone(const ideal &I, const ring &r)
    101 {
    102   int n = rVar(r);
     19#include <flip.h>
     20#include <tropicalCurves.h>
     21#include <bbcone.h>
     22
     23static bool checkPolynomialInput(const ideal I, const ring r)
     24{
     25  if (r) rTest(r);
     26  if (I && r) id_Test(I,r);
     27  return ((!I) || (I && r));
     28}
     29
     30static bool checkOrderingAndCone(const ring r, const gfan::ZCone zc)
     31{
     32  if (r)
     33  {
     34    int n = rVar(r); int* w = r->wvhdl[0];
     35    gfan::ZVector v = wvhdlEntryToZVector(n,w);
     36    if (!zc.contains(v))
     37    {
     38      std::cout << "ERROR: weight of ordering not inside Groebner cone!" << std::endl
     39                << "cone: " << std::endl
     40                << toString(&zc)
     41                << "weight: " << std::endl
     42                << v << std::endl;
     43    }
     44    return zc.contains(v);
     45  }
     46  return (zc.dimension()==0);
     47}
     48
     49static bool checkPolyhedralInput(const gfan::ZCone zc, const gfan::ZVector p)
     50{
     51  return zc.containsRelatively(p);
     52}
     53
     54groebnerCone::groebnerCone():
     55  polynomialIdeal(NULL),
     56  polynomialRing(NULL),
     57  polyhedralCone(gfan::ZCone(0)),
     58  interiorPoint(gfan::ZVector(0)),
     59  currentStrategy(NULL)
     60{
     61}
     62
     63groebnerCone::groebnerCone(const ideal I, const ring r, const tropicalStrategy& currentCase):
     64  polynomialIdeal(NULL),
     65  polynomialRing(NULL),
     66  currentStrategy(&currentCase)
     67{
     68  assume(checkPolynomialInput(I,r));
     69  if (I) polynomialIdeal = id_Copy(I,r);
     70  if (r) polynomialRing = rCopy(r);
     71
     72  int n = rVar(polynomialRing);
    10373  poly g = NULL;
    10474  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
     
    10777  gfan::ZVector tailexpw = gfan::ZVector(n);
    10878  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
    109   for (int i=0; i<IDELEMS(I); i++)
    110   {
    111     g = (poly) I->m[i]; pGetExpV(g,leadexpv);
     79  for (int i=0; i<idSize(polynomialIdeal); i++)
     80  {
     81    g = polynomialIdeal->m[i];
     82    pGetExpV(g,leadexpv);
    11283    leadexpw = intStar2ZVector(n, leadexpv);
    11384    pIter(g);
    114     while (g != NULL)
     85    while (g)
    11586    {
    11687      pGetExpV(g,tailexpv);
     
    12293  omFreeSize(leadexpv,(n+1)*sizeof(int));
    12394  omFreeSize(tailexpv,(n+1)*sizeof(int));
    124   return gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
    125 }
    126 
    127 groebnerCone::groebnerCone():
    128   I(NULL),
    129   r(NULL),
    130   c(gfan::ZCone(0)),
    131   p(gfan::ZVector(0))
    132 {
     95  polyhedralCone = gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
     96  interiorPoint = polyhedralCone.getRelativeInteriorPoint();
     97  assume(checkOrderingAndCone(polynomialRing,polyhedralCone));
     98}
     99
     100static bool checkOrderingAndWeight(const ideal I, const ring r, const gfan::ZVector w, const tropicalStrategy& currentCase)
     101{
     102  groebnerCone sigma(I,r,currentCase);
     103  gfan::ZCone zc = sigma.getPolyhedralCone();
     104  return zc.contains(w);
     105}
     106
     107groebnerCone::groebnerCone(const ideal I, const ring r, const gfan::ZVector& w, const tropicalStrategy& currentCase):
     108  polynomialIdeal(NULL),
     109  polynomialRing(NULL),
     110  currentStrategy(&currentCase)
     111{
     112  assume(checkPolynomialInput(I,r));
     113  if (I) polynomialIdeal = id_Copy(I,r);
     114  if (r) polynomialRing = rCopy(r);
     115
     116  int n = rVar(r);
     117  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
     118  gfan::ZMatrix equations = gfan::ZMatrix(0,n);
     119  int* expv = (int*) omAlloc((n+1)*sizeof(int));
     120  for (int i=0; i<idSize(polynomialIdeal); i++)
     121  {
     122    poly g = polynomialIdeal->m[i];
     123    pGetExpV(g,expv);
     124    gfan::ZVector leadexpv = intStar2ZVector(n,expv);
     125    long d = wDeg(g,polynomialRing,w);
     126    for (pIter(g); g; pIter(g))
     127    {
     128      pGetExpV(g,expv); gfan::ZVector tailexpv = intStar2ZVector(n,expv);
     129      if (wDeg(g,polynomialRing,w)==d)
     130        equations.appendRow(leadexpv-tailexpv);
     131      else
     132      {
     133        assume(wDeg(g,polynomialRing,w)<d);
     134        inequalities.appendRow(leadexpv-tailexpv);
     135      }
     136    }
     137  }
     138  omFreeSize(expv,(n+1)*sizeof(int));
     139
     140  polyhedralCone = gfan::ZCone(inequalities,equations);
     141  interiorPoint = polyhedralCone.getRelativeInteriorPoint();
     142  assume(checkOrderingAndCone(polynomialRing,polyhedralCone));
    133143}
    134144
    135145groebnerCone::groebnerCone(const groebnerCone &sigma):
    136   I(id_Copy(sigma.getIdeal(),sigma.getRing())),
    137   r(rCopy(sigma.getRing())),
    138   c(gfan::ZCone(sigma.getCone())),
    139   p(gfan::ZVector(sigma.getInteriorPoint()))
    140 {
    141 }
    142 
    143 groebnerCone::groebnerCone(const ideal &J, const ring &s, const gfan::ZCone &d, const gfan::ZVector &q):
    144   I(J),
    145   r(s),
    146   c(d),
    147   p(q)
    148 {
     146  polynomialIdeal(NULL),
     147  polynomialRing(NULL),
     148  polyhedralCone(gfan::ZCone(sigma.getPolyhedralCone())),
     149  interiorPoint(gfan::ZVector(sigma.getInteriorPoint())),
     150  currentStrategy(sigma.getTropicalStrategy())
     151{
     152  assume(checkPolynomialInput(sigma.getPolynomialIdeal(),sigma.getPolynomialRing()));
     153  assume(checkOrderingAndCone(sigma.getPolynomialRing(),sigma.getPolyhedralCone()));
     154  assume(checkPolyhedralInput(sigma.getPolyhedralCone(),sigma.getInteriorPoint()));
     155  if (sigma.getPolynomialIdeal()) polynomialIdeal = id_Copy(sigma.getPolynomialIdeal(),sigma.getPolynomialRing());
     156  if (sigma.getPolynomialRing()) polynomialRing = rCopy(sigma.getPolynomialRing());
    149157}
    150158
    151159groebnerCone::~groebnerCone()
    152160{
    153   if (I) id_Delete(&I,r);
    154   if (r) rDelete(r);
    155 }
    156 
    157 
    158 groebnerCone maximalGroebnerConeData(ideal I, const ring r)
    159 {
    160   int n = rVar(r);
    161   poly g = NULL;
    162   int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
    163   int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
    164   gfan::ZVector leadexpw = gfan::ZVector(n);
    165   gfan::ZVector tailexpw = gfan::ZVector(n);
    166   gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
    167   for (int i=0; i<IDELEMS(I); i++)
    168   {
    169     g = (poly) I->m[i]; pGetExpV(g,leadexpv);
    170     leadexpw = intStar2ZVector(n, leadexpv);
    171     pIter(g);
    172     while (g != NULL)
     161  assume(checkPolynomialInput(polynomialIdeal,polynomialRing));
     162  assume(checkOrderingAndCone(polynomialRing,polyhedralCone));
     163  assume(checkPolyhedralInput(polyhedralCone,interiorPoint));
     164  if (polynomialIdeal) id_Delete(&polynomialIdeal,polynomialRing);
     165  if (polynomialRing) rDelete(polynomialRing);
     166}
     167
     168groebnerCone& groebnerCone::operator=(const groebnerCone& sigma)
     169{
     170  assume(checkPolynomialInput(sigma.getPolynomialIdeal(),sigma.getPolynomialRing()));
     171  assume(checkOrderingAndCone(sigma.getPolynomialRing(),sigma.getPolyhedralCone()));
     172  assume(checkPolyhedralInput(sigma.getPolyhedralCone(),sigma.getInteriorPoint()));
     173  if (sigma.getPolynomialIdeal()) polynomialIdeal = id_Copy(sigma.getPolynomialIdeal(),sigma.getPolynomialRing());
     174  if (sigma.getPolynomialRing()) polynomialRing = rCopy(sigma.getPolynomialRing());
     175  polyhedralCone = sigma.getPolyhedralCone();
     176  interiorPoint = sigma.getInteriorPoint();
     177  currentStrategy = sigma.getTropicalStrategy();
     178  return *this;
     179}
     180
     181
     182/***
     183 * Returns a point in the tropical variety, if the groebnerCone contains one.
     184 * Returns an empty vector otherwise.
     185 **/
     186gfan::ZVector groebnerCone::tropicalPoint() const
     187{
     188  ideal I = polynomialIdeal;
     189  ring r = polynomialRing;
     190  gfan::ZCone zc = polyhedralCone;
     191  gfan::ZMatrix R = zc.extremeRays();
     192  assume(checkOrderingAndCone(r,zc));
     193  for (int i=0; i<R.getHeight(); i++)
     194  {
     195    ideal inI = initial(I,r,R[i]);
     196    poly s = checkForMonomialViaSuddenSaturation(inI,r);
     197    id_Delete(&inI,r);
     198    if (s == NULL)
    173199    {
    174       pGetExpV(g,tailexpv);
    175       tailexpw = intStar2ZVector(n, tailexpv);
    176       inequalities.appendRow(leadexpw-tailexpw);
    177       pIter(g);
     200      p_Delete(&s,r);
     201      return R[i];
    178202    }
    179   }
    180   omFreeSize(leadexpv,(n+1)*sizeof(int));
    181   omFreeSize(tailexpv,(n+1)*sizeof(int));
    182   gfan::ZCone zc = gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
    183   gfan::ZVector p = zc.getRelativeInteriorPoint();
    184   return groebnerCone(I,r,zc,p);
    185 }
    186 
    187 /***
    188  * Given a general ring r with any ordering, changes the ordering to wp(-w)
    189  **/
    190 void changeOrderingTo_wp(ring r, const gfan::ZVector w)
    191 {
    192   omFree(r->order);
    193   r->order  = (int*) omAlloc0(3*sizeof(int));
    194   omFree(r->block0);
    195   r->block0 = (int*) omAlloc0(3*sizeof(int));
    196   omFree(r->block1);
    197   r->block1 = (int*) omAlloc0(3*sizeof(int));
    198   for (int i=0; r->wvhdl[i]; i++) omFree(r->wvhdl[i]);
    199   omFree(r->wvhdl);
    200   r->wvhdl  = (int**) omAlloc0(3*sizeof(int*));
    201 
    202   bool ok = false;
    203   r->order[0]  = ringorder_wp;
    204   r->block0[0] = 1;
    205   r->block1[0] = r->N;
    206   r->wvhdl[0]  = ZVectorToIntStar(w,ok);
    207   r->order[1]  = ringorder_C;
    208   rComplete(r,1);
    209 }
    210 
    211 groebnerCone::groebnerCone(const ideal J, const ring s, const gfan::ZVector w)
    212 {
    213   r = rCopy(s);
    214   changeOrderingTo_wp(r,w);
    215   rChangeCurrRing(r);
    216 
    217   int k = idSize(J); I = idInit(k);
    218   nMapFunc identityMap = n_SetMap(s->cf,r->cf);
    219   for (int i=0; i<k; i++)
    220     I->m[i] = p_PermPoly(J->m[i],NULL,s,r,identityMap,NULL,0);
    221   intvec* nullVector = NULL;
    222   I = kStd(I,currQuotient,testHomog,&nullVector);
    223 
    224   c = sloppyGroebnerCone(I,r,w);
    225   p = c.getRelativeInteriorPoint();
    226 }
    227 
    228 groebnerCone::groebnerCone(const ideal J, const ring s, const gfan::ZCone d)
    229 {
    230   c = d;
    231   p = d.getRelativeInteriorPoint();
    232 
    233   r = rCopy(s);
    234   changeOrderingTo_wp(r,p);
    235   rChangeCurrRing(r);
    236 
    237   int k = idSize(J); I = idInit(k);
    238   nMapFunc identityMap = n_SetMap(s->cf,r->cf);
    239   for (int i=0; i<k; i++)
    240     I->m[i] = p_PermPoly(J->m[i],NULL,s,r,identityMap,NULL,0);
    241   intvec* nullVector = NULL;
    242   I = kStd(I,currQuotient,testHomog,&nullVector);
    243 }
    244 
    245 // /***
    246 //  * Given a general ring r with any ordering, changes the ordering to a(v),ws(-w)
    247 //  **/
    248 // bool changetoAWSRing(ring r, const gfan::ZVector v, const gfan::ZVector w)
    249 // {
    250 //   omFree(r->order);
    251 //   r->order  = (int*) omAlloc0(4*sizeof(int));
    252 //   omFree(r->block0);
    253 //   r->block0 = (int*) omAlloc0(4*sizeof(int));
    254 //   omFree(r->block1);
    255 //   r->block1 = (int*) omAlloc0(4*sizeof(int));
    256 //   for (int i=0; r->wvhdl[i]; i++)
    257 //   { omFree(r->wvhdl[i]); }
    258 //   omFree(r->wvhdl);
    259 //   r->wvhdl  = (int**) omAlloc0(4*sizeof(int*));
    260 
    261 //   bool ok = false;
    262 //   r->order[0]  = ringorder_a;
    263 //   r->block0[0] = 1;
    264 //   r->block1[0] = r->N;
    265 //   r->wvhdl[0]  = ZVectorToIntStar(v,ok);
    266 //   r->order[1]  = ringorder_ws;
    267 //   r->block0[1] = 1;
    268 //   r->block1[1] = r->N;
    269 //   r->wvhdl[1]  = ZVectorToIntStar(w,ok);
    270 //   r->order[2]  = ringorder_C;
    271 //   return ok;
    272 // }
    273 
    274 // groebnerCone::groebnerCone(const ideal J, const ring s, const gfan::ZVector w)
    275 // {
    276 //   r = rCopy(s);
    277 //   changeOrderingTo_wp(r,w);
    278 //   rChangeCurrRing(r);
    279 
    280 //   int k = idSize(J); I = idInit(k);
    281 //   nMapFunc identityMap = n_SetMap(s->cf,r->cf);
    282 //   for (int i=0; i<k; i++)
    283 //     I->m[i] = p_PermPoly(J->m[i],NULL,s,r,nMap,NULL,0);
    284 //   intvec* nullVector = NULL;
    285 //   I = kStd(I,currQuotient,testHomog,&nullVector);
    286 
    287 //   c = sloppyGroebnerCone(I,r,w);
    288 //   p = c.getRelativeInteriorPoint();
    289 // }
     203    p_Delete(&s,r);
     204  }
     205  return gfan::ZVector();
     206}
     207
     208bool groebnerCone::checkFlipConeInput(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const
     209{
     210  /* check first whether interiorPoint lies on the boundary of the cone */
     211  if (!polyhedralCone.contains(interiorPoint))
     212  {
     213    std::cout << "ERROR: interiorPoint is not contained in the Groebner cone!" << std::endl
     214              << "cone: " << std::endl
     215              << toString(&polyhedralCone)
     216              << "interiorPoint:" << std::endl
     217              << interiorPoint << std::endl;
     218    return false;
     219  }
     220  gfan::ZCone hopefullyAFacet = polyhedralCone.faceContaining(interiorPoint);
     221  if (hopefullyAFacet.dimension()!=(polyhedralCone.dimension()-1))
     222  {
     223    std::cout << "ERROR: interiorPoint is not contained in the interior of a facet!" << std::endl
     224              << "cone: " << std::endl
     225              << toString(&polyhedralCone)
     226              << "interiorPoint:" << std::endl
     227              << interiorPoint << std::endl;
     228    return false;
     229  }
     230  /* check whether facet normal points outwards */
     231  gfan::ZMatrix halfSpaceInequality(0,facetNormal.size());
     232  halfSpaceInequality.appendRow(facetNormal);
     233  gfan::ZCone halfSpace = gfan::ZCone(halfSpaceInequality,gfan::ZMatrix(0,facetNormal.size()));
     234  hopefullyAFacet = intersection(polyhedralCone, halfSpace);
     235  if (hopefullyAFacet.dimension()!=(polyhedralCone.dimension()-1))
     236  {
     237    std::cout << "ERROR: facetNormal is not pointing outwards!" << std::endl
     238              << "cone: " << std::endl
     239              << toString(&polyhedralCone)
     240              << "facetNormal:" << std::endl
     241              << facetNormal << std::endl;
     242    return false;
     243  }
     244  return true;
     245}
     246
     247/***
     248 * Given an interior point on the facet and the outer normal factor on the facet,
     249 * returns the adjacent groebnerCone sharing that facet
     250 **/
     251groebnerCone groebnerCone::flipCone(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const
     252{
     253  assume(this->checkFlipConeInput(interiorPoint,facetNormal));
     254  /* Note: the polynomial ring created will have a weighted ordering with respect to interiorPoint
     255   *   and with a weighted ordering with respect to facetNormal as tiebreaker.
     256   *   Hence it is sufficient to compute the initial form with respect to facetNormal,
     257   *   to obtain an initial form with respect to interiorPoint+e*facetNormal,
     258   *   for e>0 sufficiently small */
     259  std::pair<ideal,ring> flipped = flip(polynomialIdeal,polynomialRing,interiorPoint,facetNormal,*currentStrategy);
     260  assume(checkPolynomialInput(flipped.first,flipped.second));
     261  groebnerCone flippedCone(flipped.first, flipped.second, facetNormal, *currentStrategy);
     262  return flippedCone;
     263}
     264
     265
     266/***
     267 * Computes a relative interior point and an outer normal vector for each facet of zc
     268 **/
     269static std::pair<gfan::ZMatrix,gfan::ZMatrix> interiorPointsAndNormalsOfFacets(const gfan::ZCone zc)
     270{
     271  gfan::ZMatrix inequalities = zc.getFacets();
     272  gfan::ZMatrix equations = zc.getImpliedEquations();
     273  int r = inequalities.getHeight();
     274  int c = inequalities.getWidth();
     275
     276  /* our cone has r facets, if r==0 return empty matrices */
     277  gfan::ZMatrix relativeInteriorPoints = gfan::ZMatrix(0,c);
     278  gfan::ZMatrix outerFacetNormals = gfan::ZMatrix(0,c);
     279  if (r==0)
     280    return std::make_pair(relativeInteriorPoints,outerFacetNormals);
     281
     282  /* next we iterate over each of the r facets,
     283   * build the respective cone and add it to the list
     284   * this is the i=0 case */
     285  gfan::ZMatrix newInequalities = inequalities.submatrix(1,0,r,c);
     286  gfan::ZMatrix newEquations = equations;
     287  newEquations.appendRow(inequalities[0]);
     288  gfan::ZCone facet = gfan::ZCone(newInequalities,newEquations);
     289  relativeInteriorPoints.appendRow(facet.getRelativeInteriorPoint());
     290  assume(zc.contains(relativeInteriorPoints[0]) && !zc.containsRelatively(relativeInteriorPoints[0]));
     291  outerFacetNormals.appendRow(-inequalities[0]);
     292
     293  /* these are the cases i=1,...,r-2 */
     294  for (int i=1; i<r-1; i++)
     295  {
     296    newInequalities = inequalities.submatrix(0,0,i,c);
     297    newInequalities.append(inequalities.submatrix(i+1,0,r,c));
     298    newEquations = equations;
     299    newEquations.appendRow(inequalities[i]);
     300    facet = gfan::ZCone(newInequalities,newEquations);
     301    relativeInteriorPoints.appendRow(facet.getRelativeInteriorPoint());
     302    assume(zc.contains(relativeInteriorPoints[i]) && !zc.containsRelatively(relativeInteriorPoints[i]));
     303    outerFacetNormals.appendRow(-inequalities[i]);
     304  }
     305
     306  /* this is the i=r-1 case */
     307  newInequalities = inequalities.submatrix(0,0,r-1,c);
     308  newEquations = equations;
     309  newEquations.appendRow(inequalities[r-1]);
     310  facet = gfan::ZCone(newInequalities,newEquations);
     311  relativeInteriorPoints.appendRow(facet.getRelativeInteriorPoint());
     312  assume(zc.contains(relativeInteriorPoints[r-1]) && !zc.containsRelatively(relativeInteriorPoints[r-1]));
     313  outerFacetNormals.appendRow(-inequalities[r-1]);
     314
     315  return std::make_pair(relativeInteriorPoints,outerFacetNormals);
     316}
     317
     318
     319/***
     320 * Returns a complete list of neighboring Groebner cones.
     321 **/
     322groebnerCones groebnerCone::groebnerNeighbours() const
     323{
     324  std::pair<gfan::ZMatrix, gfan::ZMatrix> facetsData = interiorPointsAndNormalsOfFacets(polyhedralCone);
     325  gfan::ZMatrix interiorPoints = facetsData.first;
     326  gfan::ZMatrix facetNormals = facetsData.second;
     327
     328  groebnerCones neighbours;
     329  for (int i=0; i<interiorPoints.getHeight(); i++)
     330    neighbours.insert(this->flipCone(interiorPoints[i],facetNormals[i]));
     331
     332  return neighbours;
     333}
     334
     335
     336/***
     337 * Computes a relative interior point for each facet of zc
     338 **/
     339static gfan::ZMatrix interiorPointsOfFacets(const gfan::ZCone zc)
     340{
     341  gfan::ZMatrix inequalities = zc.getFacets();
     342  gfan::ZMatrix equations = zc.getImpliedEquations();
     343  int r = inequalities.getHeight();
     344  int c = inequalities.getWidth();
     345
     346  /* our cone has r facets, if r==0 return empty matrices */
     347  gfan::ZMatrix relativeInteriorPoints = gfan::ZMatrix(0,c);
     348  if (r==0) return relativeInteriorPoints;
     349
     350  /* next we iterate over each of the r facets,
     351   * build the respective cone and add it to the list
     352   * this is the i=0 case */
     353  gfan::ZMatrix newInequalities = inequalities.submatrix(1,0,r,c);
     354  gfan::ZMatrix newEquations = equations;
     355  newEquations.appendRow(inequalities[0]);
     356  gfan::ZCone facet = gfan::ZCone(newInequalities,newEquations);
     357  relativeInteriorPoints.appendRow(facet.getRelativeInteriorPoint());
     358
     359  /* these are the cases i=1,...,r-2 */
     360  for (int i=1; i<r-1; i++)
     361  {
     362    newInequalities = inequalities.submatrix(0,0,i-1,c);
     363    newInequalities.append(inequalities.submatrix(i+1,0,r,c));
     364    newEquations = equations;
     365    newEquations.appendRow(inequalities[i]);
     366    facet = gfan::ZCone(newInequalities,newEquations);
     367    relativeInteriorPoints.appendRow(facet.getRelativeInteriorPoint());
     368  }
     369
     370  /* this is the i=r-1 case */
     371  newInequalities = inequalities.submatrix(0,0,r-1,c);
     372  newEquations = equations;
     373  newEquations.appendRow(inequalities[r-1]);
     374  facet = gfan::ZCone(newInequalities,newEquations);
     375  relativeInteriorPoints.appendRow(facet.getRelativeInteriorPoint());
     376
     377  return relativeInteriorPoints;
     378}
     379
     380
     381/***
     382 * Returns a complete list of neighboring Groebner cones in the tropical variety.
     383 **/
     384groebnerCones groebnerCone::tropicalNeighbours() const
     385{
     386  gfan::ZMatrix interiorPoints = interiorPointsOfFacets(polyhedralCone);
     387  groebnerCones neighbours;
     388  for (int i=0; i<interiorPoints.getHeight(); i++)
     389  {
     390    ideal initialIdeal = initial(polynomialIdeal,polynomialRing,interiorPoints[i]);
     391    std::set<gfan::ZVector> rays = raysOfTropicalCurve(initialIdeal,polynomialRing,*currentStrategy);
     392    groebnerCones neighbours;
     393    for (std::set<gfan::ZVector>::iterator ray = rays.begin(); ray!=rays.end(); ray++)
     394      neighbours.insert(this->flipCone(interiorPoints[i],*ray));
     395  }
     396  return neighbours;
     397}
     398
    290399
    291400gfan::ZFan* toFanStar(groebnerCones setOfCones)
     
    293402  if (setOfCones.size() > 0)
    294403  {
    295     setOfGroebnerConeData::iterator cone = setOfCones.begin();
    296     gfan::ZFan* zf = new gfan::ZFan(cone->getCone().ambientDimension());
    297     for (; cone!=setOfCones.end(); cone++)
    298       zf->insert(cone->getCone());
     404    groebnerCones::iterator sigma = setOfCones.begin();
     405    gfan::ZFan* zf = new gfan::ZFan(sigma->getPolyhedralCone().ambientDimension());
     406    for (; sigma!=setOfCones.end(); sigma++)
     407      zf->insert(sigma->getPolyhedralCone());
    299408    return zf;
    300409  }
  • Singular/dyn_modules/gfanlib/groebnerCone.h

    r7aa26b2 r3c0aa5  
    1111#include <set>
    1212
    13 gfan::ZCone sloppyGroebnerCone(const poly g, const ring r, const gfan::ZVector w);
    14 gfan::ZCone sloppyGroebnerCone(const ideal I, const ring r, const gfan::ZVector w);
    15 gfan::ZCone groebnerCone(const poly g, const ring r, const gfan::ZVector w);
    16 gfan::ZCone groebnerCone(const ideal I, const ring r, const gfan::ZVector w);
    17 gfan::ZCone fullGroebnerCone(const ideal &I, const ring &r);
     13#include <tropicalStrategy.h>
     14
     15class groebnerCone;
     16struct groebnerCone_compare;
     17typedef std::set<groebnerCone,groebnerCone_compare> groebnerCones;
     18
    1819
    1920class groebnerCone
     
    2122
    2223private:
    23 
    2424  ideal polynomialIdeal;
    2525  ring polynomialRing;
    2626  gfan::ZCone polyhedralCone;
    2727  gfan::ZVector interiorPoint;
     28  const tropicalStrategy* currentStrategy;
    2829
    2930public:
     31  groebnerCone();
     32  groebnerCone(const ideal I, const ring r, const tropicalStrategy& currentCase);
     33  groebnerCone(const ideal I, const ring r, const gfan::ZVector& w, const tropicalStrategy& currentCase);
     34  groebnerCone(const groebnerCone& sigma);
     35  ~groebnerCone();
     36  groebnerCone& operator=(const groebnerCone& sigma);
    3037
    31   groebnerCone();
    32   groebnerCone(const groebnerCone &sigma);
    33   groebnerCone(const ideal I, const ring r, const gfan::ZCone c, const gfan::ZVector w);
    34   groebnerCone(const ideal I, const ring r, const gfan::ZVector w);
    35   groebnerCone(const ideal I, const ring r, const gfan::ZCone c);
    36   ~groebnerCone();
     38  void deletePolynomialIdeal()
     39  {
     40    assume ((!polynomialIdeal) || (polynomialIdeal && polynomialRing));
     41    if (polynomialIdeal) id_Delete(&polynomialIdeal,polynomialRing);
     42  }
    3743
    38   ideal getPolynomialIdeal() const { return this->polynomialIdeal; };
    39   ring getPolynomialRing() const { return this->poynomialRing; };
    40   gfan::ZCone getPolyhedralCone() const { return this->polyhedralCone; };
    41   gfan::ZVector getInteriorPoint() const { return this->interiorPoint; };
     44  void deletePolynomialRing()
     45  {
     46    assume ((!polynomialIdeal) || (polynomialIdeal && polynomialRing));
     47    if (polynomialRing) rDelete(polynomialRing);
     48    polynomialRing = NULL;
     49  }
     50
     51  ideal getPolynomialIdeal() const { return polynomialIdeal; };
     52  ring getPolynomialRing() const { return polynomialRing; };
     53  gfan::ZCone getPolyhedralCone() const { return polyhedralCone; };
     54  gfan::ZVector getInteriorPoint() const { return interiorPoint; };
     55  const tropicalStrategy* getTropicalStrategy() const {return currentStrategy; };
    4256
    4357  friend struct groebnerCone_compare;
     
    5064
    5165  /***
    52    * Given a point w in the relative interior of a facet
    53    * and a direction u away from the groebnerCone,
    54    * computes the adjacent neighbour in direction u
    55    * sharing the facet that contains w.
     66   * Given an interior point on the facet and the outer normal factor on the facet,
     67   * returns the adjacent groebnerCone sharing that facet
    5668   **/
    57   groebnerCone adjacentCone(const gfan::ZVector w, const gfan::ZVector u) const;
     69  bool checkFlipConeInput(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
     70  groebnerCone flipCone(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
     71
     72  /***
     73   * Returns a complete list of neighboring Groebner cones.
     74   **/
     75  groebnerCones groebnerNeighbours() const;
     76
     77  /***
     78   * Returns a complete list of neighboring Groebner cones in the tropical variety.
     79   **/
     80  groebnerCones tropicalNeighbours() const;
    5881};
    5982
     
    6992};
    7093
    71 groebnerCone maximalGroebnerCone(ideal I, const ring r);
    72 
    73 typedef std::set<groebnerCone,groebnerCone_compare> groebnerCones;
    74 
    7594gfan::ZFan* toFanStar(groebnerCones setOfCones);
    7695
  • Singular/dyn_modules/gfanlib/initial.cc

    r7aa26b2 r3c0aa5  
    2727
    2828/***
    29  * Returns the first terms of p of same weighted degree under w,
    30  * this is not necessarily the initial form of p with respect to w!
     29 * Checks if p is sorted with respect to w.
     30 **/
     31static bool checkSloppyInput(const poly p, const ring r, const gfan::ZVector w)
     32{
     33  long d = wDeg(p,r,w);
     34  for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
     35  {
     36    long e = wDeg(currentTerm,r,w);
     37    if (e>d)
     38      return false;
     39  }
     40  return true;
     41}
     42
     43/***
     44 * Returns the terms of p of same weighted degree under w as the leading term.
     45 * Coincides with the initial form of p with respect to w if and only if p was already
     46 * sorted with respect to w in the sense that the leading term is of highest w-degree.
    3147 **/
    3248poly sloppyInitial(const poly p, const ring r, const gfan::ZVector w)
    3349{
    34   int n = r->N;
     50  assume(checkSloppyInput(p,r,w));
     51  int n = rVar(r);
    3552  int* expv = (int*) omAlloc(n*sizeof(int));
    3653  poly q0 = p_Head(p,r);
    3754  poly q1 = q0;
    3855  long d = wDeg(p,r,w);
    39   for (poly currentTerm = p->next; wDeg(currentTerm,r,w)==d; pIter(currentTerm))
     56  for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
    4057  {
    41     pNext(q1) = p_Head(currentTerm,r);
    42     pIter(q1);
     58    if (wDeg(currentTerm,r,w) == d)
     59    {
     60      pNext(q1) = p_Head(currentTerm,r);
     61      pIter(q1);
     62    }
    4363  }
    4464  omFreeSize(expv,n*sizeof(int));
     
    4868/***
    4969 * Runs the above procedure over all generators of an ideal.
     70 * Coincides with the initial ideal of I with respect to w if and only if
     71 * the elements of I were already sorted with respect to w and
     72 * I is a standard basis form with respect to w.
    5073 **/
    5174ideal sloppyInitial(const ideal I, const ring r, const gfan::ZVector w)
  • Singular/dyn_modules/gfanlib/initial.h

    r7aa26b2 r3c0aa5  
    11#ifndef INITIAL_H
    22#define INITIAL_H
     3
     4/***
     5 * various functions to compute the initial form of polynomials and ideals
     6 **/
    37
    48#include <gfanlib/gfanlib_vector.h>
     
    610#include <Singular/ipid.h>
    711
     12/***
     13 * Computes the weighted degree of the leading monomial of p with respect to w
     14 **/
    815long wDeg(const poly p, const ring r, const gfan::ZVector w);
     16
     17/***
     18 * Returns the first terms of p of same weighted degree under w.
     19 * Coincides with the initial form of p with respect to w if and only if p was already
     20 * sorted with respect to w.
     21 **/
     22poly sloppyInitial(const poly p, const ring r, const gfan::ZVector w);
     23
     24/***
     25 * Runs the above procedure over all generators of an ideal.
     26 * Coincides with the initial ideal of I with respect to w if and only if
     27 * the elements of I were already sorted with respect to w and
     28 * I is a standard basis form with respect to w.
     29 **/
     30ideal sloppyInitial(const ideal I, const ring r, const gfan::ZVector w);
     31
    932poly initial(const poly p, const ring r, const gfan::ZVector w);
    1033ideal initial(const ideal I, const ring r, const gfan::ZVector w);
     34
     35
    1136
    1237poly initial(const poly p, const ring r);
  • Singular/dyn_modules/gfanlib/neighbours.cc

    r7aa26b2 r3c0aa5  
    1111#include <tropicalCurves.h>
    1212
     13
    1314/***
    1415 * Computes a relative interior point and outer normal vector for each facet of zc
     
    5152  return relativeInteriorPoints;
    5253}
     54
    5355
    5456/***
     
    98100}
    99101
    100 setOfGroebnerConeData groebnerNeighbours(const groebnerConeData sigma, const tropicalStrategy currentCase)
    101 {
    102   bool (*reduce)(ideal I, ring r);
     102
     103/***
     104 * Assuming that r has a ordering whose maximal Groebner cone contains w,
     105 * i.e. assuming that terms of highest w-degree can be found at the beginning of g,
     106 * computes the Groebner cone of an ideal I in r containing w relatively.
     107 **/
     108gfan::ZCone sloppyGroebnerCone(const poly g, const ring r, const gfan::ZVector w)
     109{
     110  int n = r->N;
     111  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
     112  gfan::ZMatrix equations = gfan::ZMatrix(0,n);
     113
     114  /* read off the leading exponent of g and compute its weighted degree */
     115  int* expv = (int*) omAlloc((n+1)*sizeof(int));
     116  p_GetExpV(g,expv,r);
     117  gfan::ZVector leadexp = intStar2ZVector(n,expv);
     118  long d = wDeg(g,r,w);
     119
     120  /* run through the next terms of same weighted degree and construct the equations */
     121  poly h=g->next;
     122  for (; h && wDeg(h,r,w)==d; pIter(h))
     123  {
     124    p_GetExpV(h,expv,r);
     125    equations.appendRow(leadexp-intStar2ZVector(n,expv));
     126  }
     127
     128  /* run through the remaining terms and construct the inequalities */
     129  for (; h; pIter(h))
     130  {
     131    p_GetExpV(h,expv,r);
     132    inequalities.appendRow(leadexp-intStar2ZVector(n,expv));
     133  }
     134
     135  omFreeSize(expv,(n+1)*sizeof(int));
     136  return gfan::ZCone(inequalities,equations);
     137}
     138gfan::ZCone sloppyGroebnerCone(const ideal I, const ring r, const gfan::ZVector w)
     139{
     140  int k = idSize(I);
     141  gfan::ZCone zc = gfan::ZCone(r->N);
     142  for (int i=0; i<k; i++)
     143    zc = intersection(zc,sloppyGroebnerCone(I->m[i],r,w));
     144  return zc;
     145}
     146
     147
     148groebnerCones groebnerNeighbours(const groebnerCone sigma, const tropicalStrategy currentCase)
     149{
     150  bool (*reduce)(ideal I, ring r, number p);
    103151  reduce = currentCase.reduce;
    104 
    105   gfan::ZCone zc = sigma.getCone();
     152  number p = currentCase.uniformizingParameter;
     153
     154  gfan::ZCone zc = sigma.getPolyhedralCone();
    106155  std::pair<gfan::ZMatrix, gfan::ZMatrix> facets = interiorPointsAndFacetNormals(zc);
    107156  gfan::ZMatrix interiorPoints = facets.first;
    108157  gfan::ZMatrix facetNormals = facets.second;
    109   setOfGroebnerConeData neighbours;
    110 
     158
     159  groebnerCones neighbours;
    111160  for (int i=0; i<interiorPoints.getHeight(); i++)
    112161  {
    113     std::pair<ideal,ring> flipped = flip(sigma.getIdeal(),sigma.getRing(),interiorPoints[i],facetNormals[i],currentCase);
     162    std::pair<ideal,ring> flipped = flip(sigma.getPolynomialIdeal(),sigma.getPolynomialRing(),interiorPoints[i],facetNormals[i],currentCase);
    114163    ideal I = flipped.first;
    115164    ring r = flipped.second;
    116     reduce(I,r);
    117     gfan::ZCone c = sloppyGroebnerCone(I,r,facetNormals[i]);
    118     gfan::ZVector p = c.getRelativeInteriorPoint();
    119     neighbours.insert(groebnerConeData(I,r,c,p));
     165    reduce(I,r,p);
     166    gfan::ZCone zc = sloppyGroebnerCone(I,r,facetNormals[i]);
     167    gfan::ZVector p = zc.getRelativeInteriorPoint();
     168    neighbours.insert(groebnerCone(I,r,zc,p));
    120169  }
    121170
     
    123172}
    124173
    125 setOfGroebnerConeData tropicalNeighbours(const groebnerConeData sigma, const tropicalStrategy currentCase)
    126 {
    127   gfan::ZCone zc = sigma.getCone();
     174groebnerCones tropicalNeighbours(const groebnerCone sigma, const tropicalStrategy currentCase)
     175{
     176  gfan::ZCone zc = sigma.getPolyhedralCone();
    128177  int d = zc.dimension();
    129178  gfan::ZMatrix interiorPoints = facetInteriorPoints(zc);
    130   setOfGroebnerConeData neighbours;
    131   bool (*red)(ideal I, ring r);
     179  groebnerCones neighbours;
     180  bool (*red)(ideal I, ring r, number p);
    132181  red = currentCase.reduce;
     182  number p = currentCase.uniformizingParameter;
    133183
    134184  for (int i=0; i<interiorPoints.getHeight(); i++)
    135185  {
    136186    gfan::ZVector w = interiorPoints[i];
    137     ideal I = sigma.getIdeal();
    138     ring r = sigma.getRing();
     187    ideal I = sigma.getPolynomialIdeal();
     188    ring r = sigma.getPolynomialRing();
    139189    ideal inI = initial(I,r,w);
    140     std::set<gfan::ZCone> C = tropicalCurve(inI,r,d-1,currentCase);
     190    std::set<gfan::ZCone> C = tropicalCurve(inI,r,currentCase,d-1);
    141191    std::set<gfan::ZVector> setOfRays = rays(C);
    142192
     
    146196      ideal I = flipped.first;
    147197      ring r = flipped.second;
    148       red(I,r);
     198      red(I,r,p);
    149199      gfan::ZCone c = sloppyGroebnerCone(I,r,*ray);
    150200      gfan::ZVector p = c.getRelativeInteriorPoint();
    151       neighbours.insert(groebnerConeData(I,r,c,p));
     201      neighbours.insert(groebnerCone(I,r,c,p));
    152202    }
    153203  }
  • Singular/dyn_modules/gfanlib/neighbours.h

    r7aa26b2 r3c0aa5  
    55#include <tropicalStrategy.h>
    66
    7 setOfGroebnerConeData groebnerNeighbours(const groebnerConeData sigma, const tropicalStrategy currentCase);
    8 setOfGroebnerConeData tropicalNeighbours(const groebnerConeData sigma, const tropicalStrategy currentCase);
     7gfan::ZCone sloppyGroebnerCone(const ideal I, const ring r, const gfan::ZVector w);
     8groebnerCones groebnerNeighbours(const groebnerCone sigma, const tropicalStrategy currentCase);
     9groebnerCones tropicalNeighbours(const groebnerCone sigma, const tropicalStrategy currentCase);
    910BOOLEAN tropicalNeighbours(leftv res, leftv args);
    1011
  • Singular/dyn_modules/gfanlib/ppinitialReduction.cc

    r7aa26b2 r3c0aa5  
    44#include <libpolys/polys/monomials/p_polys.h>
    55#include <singularWishlist.h>
     6#include <tropicalStrategy.h>
    67
    78#include <map>
     
    475476 * assumes that the generators of I are homogeneous in x and that p-t is in I.
    476477 **/
    477 bool ppreduceInitially(ideal I, ring r)
     478bool ppreduceInitially(ideal I, ring r, number p)
    478479{
    479480  /***
     
    500501  std::map<long,ideal>::iterator it=H.begin();
    501502  ideal Hi = it->second;
    502   assume(idSize(Hi)==1);
    503   assume(pLength(Hi->m[0])==2 && p_GetExp(Hi->m[0],1,r)==0
    504            && p_GetExp(Hi->m[0]->next,1,r)==1);
    505   number p = p_GetCoeff(Hi->m[0],r);
     503  // assume(idSize(Hi)==1);
     504  // assume(pLength(Hi->m[0])==2 && p_GetExp(Hi->m[0],1,r)==0
     505  //          && p_GetExp(Hi->m[0]->next,1,r)==1);
     506  // number p = currentCase.uniformizingParameter;
    506507  assume(!n_IsUnit(p,r->cf));
    507   idShallowDelete(&it->second);
     508  // idShallowDelete(&it->second);
    508509
    509510  /***
     
    511512   *  and all lower components
    512513   **/
    513   it++; Hi = it->second; n--;
     514  // it++; Hi = it->second; n--;
    514515  if (ppreduceInitially(Hi,p,r)) return true;
    515516
     
    560561
    561562
    562 #ifndef NDEBUG
    563 BOOLEAN ppreduceInitially4(leftv res, leftv args)
    564 {
    565   leftv u = args;
    566   if ((u != NULL) && (u->Typ() == IDEAL_CMD))
    567   {
    568     ideal I;
    569     omUpdateInfo();
    570     Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
    571     I = (ideal) u->CopyD();
    572     (void) ppreduceInitially(I,currRing);
    573     id_Delete(&I,currRing);
    574     omUpdateInfo();
    575     Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
    576     I = (ideal) u->CopyD();
    577     (void) ppreduceInitially(I,currRing);
    578     res->rtyp = IDEAL_CMD;
    579     res->data = (char*) I;
    580     return FALSE;
    581   }
    582   return TRUE;
    583 }
    584 #endif
    585 
    586 
    587 BOOLEAN ppreduceInitially(leftv res, leftv args)
    588 {
    589   leftv u = args;
    590   if ((u != NULL) && (u->Typ() == IDEAL_CMD))
    591   {
    592     ideal I = (ideal) u->CopyD();
    593     (void) ppreduceInitially(I,currRing);
    594     res->rtyp = IDEAL_CMD;
    595     res->data = (char*) I;
    596     return FALSE;
    597   }
    598   return TRUE;
    599 }
     563// #ifndef NDEBUG
     564// BOOLEAN ppreduceInitially4(leftv res, leftv args)
     565// {
     566//   leftv u = args;
     567//   if ((u != NULL) && (u->Typ() == IDEAL_CMD))
     568//   {
     569//     ideal I;
     570//     omUpdateInfo();
     571//     Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
     572//     I = (ideal) u->CopyD();
     573//     (void) ppreduceInitially(I,currRing);
     574//     id_Delete(&I,currRing);
     575//     omUpdateInfo();
     576//     Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
     577//     I = (ideal) u->CopyD();
     578//     (void) ppreduceInitially(I,currRing);
     579//     res->rtyp = IDEAL_CMD;
     580//     res->data = (char*) I;
     581//     return FALSE;
     582//   }
     583//   return TRUE;
     584// }
     585// #endif
     586
     587
     588// BOOLEAN ppreduceInitially(leftv res, leftv args)
     589// {
     590//   leftv u = args;
     591//   if ((u != NULL) && (u->Typ() == IDEAL_CMD))
     592//   {
     593//     ideal I = (ideal) u->CopyD();
     594//     (void) ppreduceInitially(I,currRing);
     595//     res->rtyp = IDEAL_CMD;
     596//     res->data = (char*) I;
     597//     return FALSE;
     598//   }
     599//   return TRUE;
     600// }
  • Singular/dyn_modules/gfanlib/ppinitialReduction.h

    r7aa26b2 r3c0aa5  
    33
    44#include <kernel/structs.h>
     5#include <tropicalStrategy.h>
    56
    67#ifndef NDEBUG
     
    1314#endif
    1415
    15 bool ppreduceInitially(ideal I, ring r);
     16bool ppreduceInitially(ideal I, ring r, number p);
    1617BOOLEAN ppreduceInitially(leftv res, leftv args);
    1718
  • Singular/dyn_modules/gfanlib/startingCone.cc

    r7aa26b2 r3c0aa5  
    33#include <tropical.h>
    44#include <initial.h>
     5#include <lift.h>
    56#include <groebnerCone.h>
    67#include <neighbours.h>
     
    1112
    1213
    13 // /***
    14 //  * Given a ring with ordering a(v'),ws(w'),
    15 //  * changes the weights to v,w
    16 //  **/
    17 // bool changeAWSWeights(ring r, gfan::ZVector v, gfan::ZVector w)
    18 // {
    19 //   omFree(r->wvhdl[0]);
    20 //   omFree(r->wvhdl[1]);
    21 //   bool ok = false;
    22 //   r->wvhdl[0]  = ZVectorToIntStar(v,ok);
    23 //   r->wvhdl[1]  = ZVectorToIntStar(w,ok);
    24 //   return ok;
    25 // }
    26 
    27 
    28 gfan::ZVector findTropicalPoint(const groebnerConeData sigma)
     14/***
     15 * checks whether sigma is contained in the tropical variety
     16 * by testing whether the initial Ideal with respect to the interior point
     17 * is monomial free.
     18 **/
     19static bool checkContainmentInTropicalVariety(const groebnerCone sigma)
    2920{
    30   ideal I = sigma.getIdeal();
    31   ring r = sigma.getRing();
    32   gfan::ZCone zc = sigma.getCone();
    33   gfan::ZMatrix R = zc.extremeRays();
    34   for (int i=0; i<R.getHeight(); i++)
     21  gfan::ZVector w = sigma.getInteriorPoint();
     22  ideal I = sigma.getPolynomialIdeal();
     23  ring r = sigma.getPolynomialRing();
     24  ideal inI = sloppyInitial(I,r,w);
     25  poly s = checkForMonomialViaSuddenSaturation(inI,r);
     26  id_Delete(&inI,r);
     27  if (s)
    3528  {
    36     // std::cout << "checking ray generated by " << R[i] << std::endl;
    37     ideal inI = initial(I,r,R[i]);
    38     poly s = checkForMonomialViaSuddenSaturation(inI,r);
    39     id_Delete(&inI,r);
    40     if (s == NULL)
    41     {
    42       p_Delete(&s,r);
    43       return R[i];
    44     }
    4529    p_Delete(&s,r);
     30    return false;
    4631  }
    47   return gfan::ZVector();
     32  return true;
    4833}
    4934
    5035
    5136/***
    52  * Computes a starting point by computing the Groebner fan,
    53  * checking each new cone whether its contains a ray in the tropical variety.
     37 * Computes a starting point by traversing the Groebner fan,
     38 * checking each cone whether it contains a ray in the tropical variety.
     39 * Returns a point in the tropical variety and a maximal Groebner cone
     40 * containing the point.
    5441 **/
    55 gfan::ZVector tropicalStartingPointViaGroebnerFan(const ideal I, const ring r, const tropicalStrategy currentCase)
     42std::pair<gfan::ZVector,groebnerCone> tropicalStartingPointViaGroebnerFan(const ideal I, const ring r, const tropicalStrategy& currentStrategy)
    5643{
    57   bool (*red)(ideal I, ring r);
    58   red = currentCase.reduce; red(I,r);
     44  currentStrategy.reduce(I,r);
    5945
    60   groebnerConeData sigma = maximalGroebnerConeData(id_Copy(I,r),rCopy(r));
    61   gfan::ZVector startingPoint = findTropicalPoint(sigma);
     46  const groebnerCone sigma(I,r,currentStrategy);
     47  gfan::ZVector startingPoint = sigma.tropicalPoint();
    6248  if (startingPoint.size() > 0)
    63     return startingPoint;
    64   setOfGroebnerConeData groebnerFan;
    65   setOfGroebnerConeData workingList;
     49    return std::make_pair(startingPoint,sigma);
     50  groebnerCones groebnerFan;
     51  groebnerCones workingList;
    6652  workingList.insert(sigma);
    6753
    6854  while (!workingList.empty())
    6955  {
    70     sigma = *(workingList.begin());
    71     setOfGroebnerConeData neighbours = groebnerNeighbours(sigma,currentCase);
    72     for (setOfGroebnerConeData::iterator tau = neighbours.begin(); tau!=neighbours.end(); tau++)
     56    const groebnerCone sigma = *(workingList.begin());
     57    groebnerCones neighbours = sigma.groebnerNeighbours();
     58    for (groebnerCones::iterator tau = neighbours.begin(); tau!=neighbours.end(); tau++)
    7359    {
    7460      if (groebnerFan.count(*tau) == 0)
     
    7662        if (workingList.count(*tau) == 0)
    7763        {
    78           startingPoint = findTropicalPoint(*tau);
     64          startingPoint = tau->tropicalPoint();
    7965          if (startingPoint.size() > 0)
    80             return startingPoint;
     66            return std::make_pair(startingPoint,*tau);
    8167        }
    8268        workingList.insert(*tau);
    8369      }
    8470    }
    85     groebnerFan.insert(groebnerConeData(sigma));
     71    groebnerFan.insert(sigma);
    8672    workingList.erase(sigma);
    8773  }
    88   return startingPoint;
     74  gfan::ZVector emptyVector = gfan::ZVector(0);
     75  groebnerCone emptyCone = groebnerCone();
     76  return std::pair<gfan::ZVector,groebnerCone>(emptyVector,emptyCone);
    8977}
    9078
    91 BOOLEAN tropicalStartingPoint0(leftv res, leftv args)
     79BOOLEAN tropicalStartingPoint(leftv res, leftv args)
    9280{
    9381  leftv u = args;
    94   ideal I = (ideal) u->CopyD();
    95   gfan::ZVector startingPoint = tropicalStartingPoint(I,currRing,nonValuedCase);
    96   id_Delete(&I, currRing);
    97   res->rtyp = BIGINTMAT_CMD;
    98   res->data = (char*) zVectorToBigintmat(startingPoint);
    99   return FALSE;
     82  if ((u!=NULL) && (u->Typ()==IDEAL_CMD))
     83  {
     84    ideal I = (ideal) u->Data();
     85    tropicalStrategy currentStrategy(I,currRing);
     86    std::pair<gfan::ZVector,groebnerCone> startingData = tropicalStartingPointViaGroebnerFan(I,currRing,currentStrategy);
     87    gfan::ZVector startingPoint = startingData.first;
     88    res->rtyp = BIGINTMAT_CMD;
     89    res->data = (void*) zVectorToBigintmat(startingPoint);
     90    return FALSE;
     91  }
     92  WerrorS("tropicalStartingPoint: unexpected parameters");
     93  return TRUE;
    10094}
    10195
    102 BOOLEAN tropicalStartingPoint1(leftv res, leftv args)
     96/***
     97 * returs the lineality space of the Groebner fan
     98 **/
     99static gfan::ZCone linealitySpaceOfGroebnerFan(const ideal I, const ring r)
     100{
     101  int n = rVar(r);
     102  gfan::ZMatrix equations = gfan::ZMatrix(0,n);
     103  int* expv = (int*) omAlloc((n+1)*sizeof(int));
     104  int k = idSize(I);
     105  for (int i=0; i<k; i++)
     106  {
     107    poly g = I->m[i];
     108    p_GetExpV(g,expv,r);
     109    gfan::ZVector leadexp = intStar2ZVector(n,expv);
     110    for (pIter(g); g; pIter(g))
     111    {
     112      p_GetExpV(g,expv,r);
     113      equations.appendRow(leadexp-intStar2ZVector(n,expv));
     114    }
     115  }
     116  omFreeSize(expv,(n+1)*sizeof(int));
     117  return gfan::ZCone(gfan::ZMatrix(0,n),equations);
     118}
     119
     120
     121/***
     122 * Computes a starting cone in the tropical variety.
     123 **/
     124groebnerCone tropicalStartingCone(const ideal I, const ring r, const tropicalStrategy& currentStrategy)
     125{
     126  ring s = rCopy(r);
     127  int k = idSize(I); ideal inI = idInit(k);
     128  nMapFunc identityMap = n_SetMap(r->cf,s->cf);
     129  for (int i=0; i<k; i++)
     130    inI->m[i] = p_PermPoly(I->m[i],NULL,r,s,identityMap,NULL,0);
     131
     132  gfan::ZCone zc = linealitySpaceOfGroebnerFan(inI,s);
     133  gfan::ZVector potentialStartingPoint; groebnerCone ambientMaximalCone;
     134  while (zc.dimension()<currentStrategy.getDimensionOfIdeal())
     135  {
     136    std::pair<gfan::ZVector,groebnerCone> startingData = tropicalStartingPointViaGroebnerFan(inI,s,currentStrategy);
     137    potentialStartingPoint = startingData.first;
     138    ambientMaximalCone = groebnerCone(startingData.second);
     139
     140    id_Delete(&inI,s); rDelete(s);
     141    inI = ambientMaximalCone.getPolynomialIdeal();
     142    s = ambientMaximalCone.getPolynomialRing();
     143
     144    inI = sloppyInitial(inI,s,potentialStartingPoint);
     145    zc = linealitySpaceOfGroebnerFan(inI,s);
     146  }
     147
     148  ideal J = lift(I,r,inI,s);
     149  id_Delete(&inI,s);
     150  groebnerCone tropicalStartingCone(J,s,potentialStartingPoint,currentStrategy);
     151  id_Delete(&J,s);
     152
     153  assume(checkContainmentInTropicalVariety(tropicalStartingCone));
     154  return tropicalStartingCone;
     155}
     156
     157
     158BOOLEAN tropicalStartingCone(leftv res, leftv args)
    103159{
    104160  leftv u = args;
    105   ideal I = (ideal) u->CopyD();
    106   gfan::ZVector startingPoint = tropicalStartingPoint(I,currRing,valuedCase);
    107   id_Delete(&I, currRing);
    108   res->rtyp = BIGINTMAT_CMD;
    109   res->data = (char*) zVectorToBigintmat(startingPoint);
    110   return FALSE;
     161  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
     162  {
     163    ideal I = (ideal) u->CopyD();
     164    leftv v = u->next;
     165    if ((v != NULL) && (v->Typ() == NUMBER_CMD))
     166    {
     167      number p = (number) v->Data();
     168      leftv w = v->next;
     169      if (w==NULL)
     170      {
     171        tropicalStrategy valuedCase(I,p,currRing);
     172        groebnerCone sigma = tropicalStartingCone(I,currRing,valuedCase);
     173        res->rtyp = coneID;
     174        res->data = (char*) new gfan::ZCone(sigma.getPolyhedralCone());
     175        return FALSE;
     176      }
     177    }
     178    else
     179    {
     180      if (v==NULL)
     181      {
     182        tropicalStrategy nonValuedCase(I,currRing);
     183        groebnerCone sigma = tropicalStartingCone(I,currRing,nonValuedCase);
     184        res->rtyp = coneID;
     185        res->data = (char*) new gfan::ZCone(sigma.getPolyhedralCone());
     186        return FALSE;
     187      }
     188    }
     189  }
     190  WerrorS("tropicalStartingCone: unexpected parameters");
     191  return TRUE;
    111192}
    112 
    113 gfan::ZCone tropicalStartingCone(const ideal I, const ring r, const tropicalStrategy currentCase)
    114 {
    115   // groebnerConeData sigma = maximalGroebnerConeData(I,r);
    116   // gfan::ZCone zc = sigma.getCone();
    117   // if (zc.dimensionOfLinealitySpace == d)
    118   //   return zc.linealitySpace();
    119   // gfan::ZVector w = tropicalStartingPoint(I,r,currentCase);
    120   // ideal inI = initial(I,w,r);
    121   // gfan::ZCone zd =
    122 
    123   gfan::ZVector w = tropicalStartingPoint(I,r,currentCase);
    124   if (w.size() == 0)
    125   {
    126     gfan::ZCone zc = fullGroebnerCone(I,r);
    127     return zc.linealitySpace();
    128   }
    129   ideal inI = initial(I,r,w);
    130   std::set<gfan::ZCone> T = tropicalCurve(inI,r,0,currentCase);
    131   return *T.begin();
    132 }
    133 
    134 BOOLEAN tropicalStartingCone0(leftv res, leftv args)
    135 {
    136   leftv u = args;
    137   ideal I = (ideal) u->CopyD();
    138   gfan::ZCone startingCone = tropicalStartingCone(I,currRing,nonValuedCase);
    139   id_Delete(&I, currRing);
    140   res->rtyp = coneID;
    141   res->data = (char*) new gfan::ZCone(startingCone);
    142   return FALSE;
    143 }
    144 
    145 BOOLEAN tropicalStartingCone1(leftv res, leftv args)
    146 {
    147   leftv u = args;
    148   ideal I = (ideal) u->CopyD();
    149   gfan::ZCone startingCone = tropicalStartingCone(I,currRing,valuedCase);
    150   id_Delete(&I, currRing);
    151   res->rtyp = coneID;
    152   res->data = (char*) new gfan::ZCone(startingCone);
    153   return FALSE;
    154 }
  • Singular/dyn_modules/gfanlib/startingCone.h

    r7aa26b2 r3c0aa5  
    55#include <tropicalStrategy.h>
    66
    7 gfan::ZVector containsTropicalPoint(const groebnerConeData sigma);
    8 gfan::ZVector tropicalStartingPoint(ideal I, ring r);
    9 gfan::ZCone tropicalStartingCone(const ideal I, const ring r, const tropicalStrategy currentCase);
     7groebnerCone tropicalStartingCone(const ideal I, const ring r, const tropicalStrategy& currentCase);
    108
    11 BOOLEAN tropicalStartingPoint0(leftv res, leftv args);
    12 BOOLEAN tropicalStartingPoint1(leftv res, leftv args);
    13 
    14 BOOLEAN tropicalStartingCone0(leftv res, leftv args);
    15 BOOLEAN tropicalStartingCone1(leftv res, leftv args);
     9BOOLEAN tropicalStartingPoint(leftv res, leftv args);
     10BOOLEAN tropicalStartingCone(leftv res, leftv args);
    1611
    1712#endif
  • Singular/dyn_modules/gfanlib/tropical.cc

    r7aa26b2 r3c0aa5  
    1111#include <tropicalCurves.h>
    1212#include <neighbours.h>
    13 #include <tropicalVarietyOfPolynomials.h>
    1413#include <tropicalStrategy.h>
    1514#include <startingCone.h>
    1615#include <tropicalVariety.h>
     16
    1717
    1818BOOLEAN homogeneitySpace(leftv res, leftv args)
     
    163163
    164164
    165 // gfan::ZCone* startingCone(ideal I)
    166 // {
    167 //   I = kStd(I,NULL,isNotHomog,NULL);
    168 //   gfan::ZCone* zc = maximalGroebnerCone(currRing,I);
    169 //   gfan::ZMatrix rays = zc->extremeRays();
    170 //   gfan::ZVector v;
    171 //   for (int i=0; i<rays.getHeight(); i++)
    172 //   {
    173 //     v = rays[i];
    174 //   }
    175 //   return zc;
    176 // }
    177 
    178 
    179165void tropical_setup(SModulFunctions* p)
    180166{
    181   valuedCase = initializeValuedCase();
    182   nonValuedCase = initializeNonValuedCase();
    183167  p->iiAddCproc("","groebnerCone",FALSE,groebnerCone);
    184168  p->iiAddCproc("","maximalGroebnerCone",FALSE,maximalGroebnerCone);
     
    186170  // p->iiAddCproc("","tropicalNeighbours",FALSE,tropicalNeighbours);
    187171#ifndef NDEBUG
    188   p->iiAddCproc("","initial0",FALSE,initial0);
    189   p->iiAddCproc("","pppReduce",FALSE,pppReduce);
    190   p->iiAddCproc("","ppreduceInitially0",FALSE,ppreduceInitially0);
    191   p->iiAddCproc("","ppreduceInitially1",FALSE,ppreduceInitially1);
    192   p->iiAddCproc("","ppreduceInitially2",FALSE,ppreduceInitially2);
    193   p->iiAddCproc("","ppreduceInitially3",FALSE,ppreduceInitially3);
    194   p->iiAddCproc("","ppreduceInitially4",FALSE,ppreduceInitially4);
    195   p->iiAddCproc("","ttpReduce",FALSE,ttpReduce);
    196   p->iiAddCproc("","ttreduceInitially0",FALSE,ttreduceInitially0);
    197   p->iiAddCproc("","ttreduceInitially1",FALSE,ttreduceInitially1);
    198   p->iiAddCproc("","ttreduceInitially2",FALSE,ttreduceInitially2);
    199   p->iiAddCproc("","ttreduceInitially3",FALSE,ttreduceInitially3);
    200   p->iiAddCproc("","ttreduceInitially4",FALSE,ttreduceInitially4);
    201   p->iiAddCproc("","checkForMonomial",FALSE,checkForMonomial);
    202   p->iiAddCproc("","dwr0",FALSE,dwr0);
    203   p->iiAddCproc("","witness0",FALSE,witness0);
    204   p->iiAddCproc("","tropicalVariety00",FALSE,tropicalVariety00);
    205   p->iiAddCproc("","tropicalVariety01",FALSE,tropicalVariety01);
    206   p->iiAddCproc("","tropicalCurve0",FALSE,tropicalCurve0);
    207   p->iiAddCproc("","tropicalCurve1",FALSE,tropicalCurve1);
    208   p->iiAddCproc("","tropicalStartingPoint0",FALSE,tropicalStartingPoint0);
    209   p->iiAddCproc("","tropicalStartingPoint1",FALSE,tropicalStartingPoint1);
    210   p->iiAddCproc("","tropicalStartingCone0",FALSE,tropicalStartingCone0);
    211   p->iiAddCproc("","tropicalStartingCone1",FALSE,tropicalStartingCone1);
     172  // p->iiAddCproc("","initial0",FALSE,initial0);
     173  // p->iiAddCproc("","pppReduce",FALSE,pppReduce);
     174  // p->iiAddCproc("","ppreduceInitially0",FALSE,ppreduceInitially0);
     175  // p->iiAddCproc("","ppreduceInitially1",FALSE,ppreduceInitially1);
     176  // p->iiAddCproc("","ppreduceInitially2",FALSE,ppreduceInitially2);
     177  // p->iiAddCproc("","ppreduceInitially3",FALSE,ppreduceInitially3);
     178  // p->iiAddCproc("","ppreduceInitially4",FALSE,ppreduceInitially4);
     179  // p->iiAddCproc("","ttpReduce",FALSE,ttpReduce);
     180  // p->iiAddCproc("","ttreduceInitially0",FALSE,ttreduceInitially0);
     181  // p->iiAddCproc("","ttreduceInitially1",FALSE,ttreduceInitially1);
     182  // p->iiAddCproc("","ttreduceInitially2",FALSE,ttreduceInitially2);
     183  // p->iiAddCproc("","ttreduceInitially3",FALSE,ttreduceInitially3);
     184  // p->iiAddCproc("","ttreduceInitially4",FALSE,ttreduceInitially4);
     185  // p->iiAddCproc("","checkForMonomial",FALSE,checkForMonomial);
     186  // p->iiAddCproc("","dwr0",FALSE,dwr0);
     187  // p->iiAddCproc("","witness0",FALSE,witness0);
     188  // p->iiAddCproc("","tropicalVariety00",FALSE,tropicalVariety00);
     189  // p->iiAddCproc("","tropicalVariety01",FALSE,tropicalVariety01);
     190  // p->iiAddCproc("","tropicalCurve0",FALSE,tropicalCurve0);
     191  // p->iiAddCproc("","tropicalCurve1",FALSE,tropicalCurve1);
     192  p->iiAddCproc("","tropicalStartingPoint",FALSE,tropicalStartingPoint);
     193  p->iiAddCproc("","tropicalStartingCone",FALSE,tropicalStartingCone);
    212194  p->iiAddCproc("","tropicalVariety",FALSE,tropicalVariety);
    213195#endif //NDEBUG
    214   p->iiAddCproc("","ppreduceInitially",FALSE,ppreduceInitially);
    215   p->iiAddCproc("","ttreduceInitially",FALSE,ttreduceInitially);
    216   p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
     196  // p->iiAddCproc("","ppreduceInitially",FALSE,ppreduceInitially);
     197  // p->iiAddCproc("","ttreduceInitially",FALSE,ttreduceInitially);
     198  // p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
    217199}
  • Singular/dyn_modules/gfanlib/tropicalCurves.cc

    r7aa26b2 r3c0aa5  
    88#include <witness.h>
    99#include <tropicalStrategy.h>
    10 #include <tropicalCurves.h>
     10#include <tropicalVarietyOfPolynomials.h>
    1111#include <set>
    1212#ifndef NDEBUG
     
    1414#endif
    1515
    16 std::set<gfan::ZCone> intersect(const std::set<gfan::ZCone> setA, const std::set<gfan::ZCone> setB, int d)
    17 {
     16/***
     17 * Given two sets of cones A,B and a dimensional bound d,
     18 * computes the intersections of all cones of A with all cones of B,
     19 * and throws away those of lower dimension than d.
     20 **/
     21static std::set<gfan::ZCone> intersect(const std::set<gfan::ZCone> setA,
     22                                       const std::set<gfan::ZCone> setB,
     23                                       int d=0)
     24{
     25  if (setA.empty())
     26    return setB;
     27  if (setB.empty())
     28    return setA;
    1829  std::set<gfan::ZCone> setAB;
    1930  for (std::set<gfan::ZCone>::iterator coneOfA=setA.begin(); coneOfA!=setA.end(); coneOfA++)
     
    3344
    3445/***
    35  * Given a relative interior point w of a cone in the weight space
    36  * and generators E of its span, returns a ring with an ordering that coincides
    37  * with a weighted ordering on the ideal with respect to a weight in the cone
    38  * which does not lie on a Groebner cone of lower dimension than that of the cone.
    39  **/
    40 static ring ringWithGenericlyWeightedOrdering(const ring r,
    41                                               const gfan::ZVector w,
    42                                               const gfan::ZMatrix E,
    43                                               const tropicalStrategy currentCase)
    44 {
    45   gfan::ZVector (*adjustWeight)(gfan::ZVector v, gfan::ZVector w);
    46   adjustWeight = currentCase.adjustWeightUnderHomogeneity;
    47 
     46 * Given a weight w and k more weights E[1], ..., E[h] stored as row vectors of a matrix E,
     47 + returns a ring whose ordering is weighted with respect to any
     48 * w+\varepsilon_1*E[1]+...+\varepsilon_k*E[h] for \varepsilon_i sufficiently small.
     49 * In particular, if E[1], ..., E[h] generate a vector space of dimension d,
     50 * this ordering can be used to compute a Groebner cone of dimension d containing w.
     51 **/
     52static ring genericlyWeightedOrdering(const ring r, const gfan::ZVector w, const gfan::ZMatrix E,
     53                                      const tropicalStrategy& currentStrategy)
     54{
    4855  int n = r->N;
    4956  int h = E.getHeight();
    5057
    51   // create a copy of r and delete its old ordering
     58  /* create a copy s of r and delete its ordering */
    5259  ring s = rCopy0(r);
    5360  omFree(s->order);
     
    6168  s->wvhdl  = (int**) omAlloc0((h+3)*sizeof(int*));
    6269
    63   // construct a new ordering and keep an eye out for weight overflows
     70  /* construct a new ordering and keep an eye out for weight overflows */
    6471  bool overflow;
    6572  s->order[0] = ringorder_a;
     
    7279    s->block0[j] = 1;
    7380    s->block1[j] = n;
    74     s->wvhdl[j] = ZVectorToIntStar(adjustWeight(E[j-1],w),overflow);
     81    s->wvhdl[j] = ZVectorToIntStar(currentStrategy.adjustWeightUnderHomogeneity(E[j-1],w),overflow);
    7582  }
    7683  s->order[h] = ringorder_wp;
    7784  s->block0[h] = 1;
    7885  s->block1[h] = n;
    79   s->wvhdl[h] = ZVectorToIntStar(adjustWeight(E[h-1],w),overflow);
     86  s->wvhdl[h] = ZVectorToIntStar(currentStrategy.adjustWeightUnderHomogeneity(E[h-1],w),overflow);
    8087  s->order[h+1] = ringorder_C;
    8188
    82   rComplete(s,1);
    83   // Return s.
    8489  if (overflow)
    8590    throw 0; //weightOverflow;
     91
     92  /* complete the ring and return it */
     93  rComplete(s,1);
    8694  return s;
    8795}
     
    8997
    9098/***
    91  * Computes the tropical curve of an x-homogeneous ideal I in ring r
    92  **/
    93 std::set<gfan::ZCone> tropicalCurve(const ideal I, const ring r, const int d,
    94                          const tropicalStrategy currentCase)
     99 * Given an ideal I which is homogeneous in the later variables of r,
     100 * whose tropicalization is, modulo its lineality space, a tropical Curve,
     101 * computes the said tropicalization.
     102 * At the end of the computation, I will contain a tropical basis.
     103 * If the dimension of the tropical variety is known beforehand,
     104 * it can be provided in order to speed up the computation.
     105 **/
     106std::set<gfan::ZCone> tropicalCurve(ideal I, const ring r, const tropicalStrategy currentStrategy)
    95107{
    96108  int k = idSize(I);
    97   std::set<gfan::ZCone> (*tropicalVariety)(const poly &g, const ring &r);
    98   tropicalVariety = currentCase.tropicalVarietyOfPolynomial;
    99   // Compute the common refinement of the tropical varieties
    100   // of the generating set
    101   std::set<gfan::ZCone> C = tropicalVariety(I->m[0],r);
     109
     110  /* Compute the common refinement over all tropical varieties
     111   * of the polynomials in the generating set */
     112  std::set<gfan::ZCone> C = tropicalVariety(I->m[0],r,currentStrategy);
    102113  for (int i=1; i<k; i++)
    103     C = intersect(C,tropicalVariety(I->m[i],r),d);
    104 
    105   // cycle through all maximal cones of the refinement
    106   // pick a monomial ordering corresponding to a generic weight vector in it
    107   // check if the initial ideal is monomial free
    108   // if the initial ideal is not monomial free, pick a witness for the monomial
    109   // and intersect with its corresponding tropical variety
    110   // if all initial ideals are monomial free, the we have our tropical curve
    111   poly mon = NULL;
    112   do
    113   {
    114     for (std::set<gfan::ZCone>::iterator zc=C.begin(); zc!=C.end(); zc++)
     114    C = intersect(C,tropicalVariety(I->m[i],r,currentStrategy),currentStrategy.getDimensionOfIdeal());
     115
     116  /* Cycle through all maximal cones of the refinement.
     117   * Pick a monomial ordering corresponding to a generic weight vector in it
     118   * and check if the initial ideal is monomial free, generic meaning
     119   * that it lies in a maximal Groebner cone in the maximal cone of the refinement.
     120   * If the initial ideal is not monomial free, compute a witness for the monomial
     121   * and compute the common refinement with its tropical variety.
     122   * If all initial ideals are monomial free, then we have our tropical curve */
     123  for (std::set<gfan::ZCone>::iterator zc=C.begin(); zc!=C.end(); zc++)
     124  {
     125    gfan::ZVector v = zc->getRelativeInteriorPoint();
     126    gfan::ZMatrix E = zc->generatorsOfSpan();
     127
     128    ring s = genericlyWeightedOrdering(r,v,E,currentStrategy);
     129    nMapFunc nMap = n_SetMap(r->cf,s->cf);
     130    ideal Is = idInit(k);
     131    for (int j=0; j<k; j++)
     132      Is->m[j] = p_PermPoly(I->m[j],NULL,r,s,nMap,NULL,0);
     133    ideal inIs = initial(Is,s,E[E.getHeight()-1]);
     134
     135    poly mons = checkForMonomialViaSuddenSaturation(inIs,s);
     136    if (mons)
    115137    {
    116       gfan::ZVector v = zc->getRelativeInteriorPoint();
    117       gfan::ZMatrix E = zc->generatorsOfSpan();
    118 
    119       // std::cout << "interiorPoint: " << v << std::endl;
    120       // std::cout << "generators of span: " << E << std::endl;
    121       // ideal inIr = initial(I,r,E[E.getHeight()-1]);
    122       ring s = ringWithGenericlyWeightedOrdering(r,v,E,currentCase);
    123       nMapFunc nMap = n_SetMap(r->cf,s->cf);
    124       ideal Is = idInit(k);
    125       for (int j=0; j<k; j++)
    126       {
    127         Is->m[j] = p_PermPoly(I->m[j],NULL,r,s,nMap,NULL,0);
    128         p_Test(Is->m[j],s);
    129       }
    130       ideal inIs = initial(Is,s,E[E.getHeight()-1]);
    131       id_Delete(&Is,s);
    132       // id_Delete(&inIr,r);
    133 
    134       mon = checkForMonomialViaSuddenSaturation(inIs,s);
    135       if (mon)
    136       {
    137         ideal Is = idInit(k);
    138         for (int j=0; j<k; j++)
    139         {
    140           Is->m[j] = p_PermPoly(I->m[j],NULL,r,s,nMap,NULL,0);
    141           p_Test(Is->m[j],s);
    142         }
    143         poly g = witness(mon,Is,inIs,s);
    144         C = intersect(C,tropicalVariety(g,s),d);
    145         p_Delete(&mon,s);
    146         p_Delete(&g,s);
    147         id_Delete(&inIs,s);
    148         id_Delete(&Is,s);
    149         rDelete(s);
    150         break;
    151       }
    152       id_Delete(&inIs,s);
    153       rDelete(s);
     138      poly gs = witness(mons,Is,inIs,s);
     139      C = intersect(C,tropicalVariety(gs,s,currentStrategy),currentStrategy.getDimensionOfIdeal());
     140      nMapFunc mMap = n_SetMap(s->cf,r->cf);
     141      poly gr = p_PermPoly(gs,NULL,s,r,mMap,NULL,0);
     142      idInsertPoly(I,gr);
     143      k++;
     144      p_Delete(&mons,s);
     145      p_Delete(&gs,s);
     146      zc = C.begin();
    154147    }
    155   } while (mon);
     148    id_Delete(&Is,s);
     149    id_Delete(&inIs,s);
     150    rDelete(s);
     151  }
    156152  return C;
    157153}
     154
     155
     156std::set<gfan::ZVector> raysOfTropicalCurve(ideal I, const ring r, const tropicalStrategy& currentStrategy)
     157{
     158  std::set<gfan::ZCone> C = tropicalCurve(I,r,currentStrategy);
     159  std::set<gfan::ZVector> raysOfC;
     160  for (std::set<gfan::ZCone>::iterator zc=C.begin(); zc!=C.end(); zc++)
     161  {
     162    assume(zc->dimensionOfLinealitySpace()+1 == zc->dimension());
     163    raysOfC.insert(zc->semiGroupGeneratorOfRay());
     164  }
     165  return raysOfC;
     166}
     167
     168
    158169/***
    159170 * Computes the tropical curve of an x-homogeneous ideal I
    160171 * which is weighted homogeneous with respect to weight w in ring r
    161172 **/
    162 
    163173#ifndef NDEBUG
    164 BOOLEAN tropicalCurve0(leftv res, leftv args)
    165 {
    166   leftv u = args;
    167   ideal I = (ideal) u->CopyD();
    168   leftv v = u->next;
    169   int d = (long)(int) v->CopyD();
    170   tropicalStrategy currentCase = nonValuedCase;
    171   std::set<gfan::ZCone> C = tropicalCurve(I,currRing,d,currentCase);
    172   id_Delete(&I,currRing);
    173   omUpdateInfo();
    174   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
    175   res->rtyp = fanID;
    176   res->data = (char*) toFanStar(C);
    177   return FALSE;
    178 }
    179 BOOLEAN tropicalCurve1(leftv res, leftv args)
    180 {
    181   leftv u = args;
    182   ideal I = (ideal) u->CopyD();
    183   leftv v = u->next;
    184   int d = (long)(int) v->CopyD();
    185   tropicalStrategy currentCase = valuedCase;
    186   std::set<gfan::ZCone> C = tropicalCurve(I,currRing,d,currentCase);
    187   id_Delete(&I,currRing);
    188   omUpdateInfo();
    189   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
    190   res->rtyp = fanID;
    191   res->data = (char*) toFanStar(C);
    192   return FALSE;
    193 }
     174// BOOLEAN tropicalCurve0(leftv res, leftv args)
     175// {
     176//   leftv u = args;
     177//   ideal I = (ideal) u->CopyD();
     178//   leftv v = u->next;
     179//   int d = (int)(long) v->CopyD();
     180//   tropicalStrategy currentCase = nonValuedCase;
     181//   std::set<gfan::ZCone> C = tropicalCurve(I,currRing,d,currentCase);
     182//   id_Delete(&I,currRing);
     183//   omUpdateInfo();
     184//   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
     185//   res->rtyp = fanID;
     186//   res->data = (char*) toFanStar(C);
     187//   return FALSE;
     188// }
     189// BOOLEAN tropicalCurve1(leftv res, leftv args)
     190// {
     191//   leftv u = args;
     192//   ideal I = (ideal) u->CopyD();
     193//   leftv v = u->next;
     194//   int d = (int)(long) v->CopyD();
     195//   tropicalStrategy currentCase = valuedCase;
     196//   std::set<gfan::ZCone> C = tropicalCurve(I,currRing,d,currentCase);
     197//   id_Delete(&I,currRing);
     198//   omUpdateInfo();
     199//   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
     200//   res->rtyp = fanID;
     201//   res->data = (char*) toFanStar(C);
     202//   return FALSE;
     203// }
    194204#endif
  • Singular/dyn_modules/gfanlib/tropicalCurves.h

    r7aa26b2 r3c0aa5  
    77#include <tropicalStrategy.h>
    88
    9 std::set<gfan::ZCone> tropicalCurve(const ideal I, const ring r, const int d, const tropicalStrategy currentCase);
     9std::set<gfan::ZCone> tropicalCurve(const ideal I, const ring r, const tropicalStrategy& currentStrategy);
     10std::set<gfan::ZVector> raysOfTropicalCurve(ideal I, const ring r, const tropicalStrategy& currentStrategy);
    1011
    1112#ifndef NDEBUG
  • Singular/dyn_modules/gfanlib/tropicalStrategy.cc

    r7aa26b2 r3c0aa5  
    11#include <tropicalStrategy.h>
    2 #include <tropicalVarietyOfPolynomials.h>
    32#include <adjustWeights.h>
    43#include <ppinitialReduction.h>
    54#include <ttinitialReduction.h>
    65
    7 tropicalStrategy valuedCase, nonValuedCase;
     6// for various commands in dim(ideal I, ring r):
     7#include <kernel/ideals.h>
     8#include <kernel/stairc.h>
    89
    9 tropicalStrategy initializeValuedCase()
     10/***
     11 * Computes the dimension of an ideal I in ring r
     12 * Copied from jjDim in iparith.cc
     13 **/
     14static int dim(ideal I, ring r)
    1015{
    11   tropicalStrategy currentCase;
    12   currentCase.tropicalVarietyOfPolynomial = lowerTropicalVarietyOfPolynomial;
    13   currentCase.adjustWeightForHomogeneity = valued_adjustWeightForHomogeneity;
    14   currentCase.adjustWeightUnderHomogeneity = valued_adjustWeightUnderHomogeneity;
    15   currentCase.reduce = ppreduceInitially;
    16   return currentCase;
     16  ring origin = currRing;
     17  if (origin != r)
     18    rChangeCurrRing(r);
     19  int d;
     20  if (rField_is_Ring(currRing))
     21  {
     22    int i = idPosConstant(I);
     23    if ((i != -1) && (n_IsUnit(pGetCoeff(I->m[i]),currRing->cf)))
     24      return -1;
     25    ideal vv = id_Head(I,currRing);
     26    if (i != -1) pDelete(&vv->m[i]);
     27    d = scDimInt(vv, currQuotient);
     28    if (rField_is_Ring_Z(currRing) && (i==-1)) d++;
     29    idDelete(&vv);
     30    return d;
     31  }
     32  else
     33    d = scDimInt(I,currQuotient);
     34  if (origin != r)
     35    rChangeCurrRing(origin);
     36  return d;
    1737}
    1838
    19 static bool doNothing(ideal I, ring r)
     39
     40/***
     41 * Initializes all relevant structures and information for the valued case,
     42 * i.e. computing a tropical variety over the rational numbers with p-adic valuation
     43 **/
     44tropicalStrategy::tropicalStrategy(ideal J, number q, ring s)
     45{
     46  /* assume that the ground field of the originalRing is Q */
     47  assume(rField_is_Q(s));
     48  originalRing = rCopy(s);
     49
     50  /* replace Q with Z for the startingRing */
     51  startingRing = rCopy0(originalRing);
     52  nKillChar(startingRing->cf);
     53  startingRing->cf = nInitChar(n_Z,NULL);
     54  rComplete(startingRing);
     55
     56  /* map p into the new coefficient domain */
     57  nMapFunc nMap = n_SetMap(originalRing->cf,startingRing->cf);
     58  uniformizingParameter = nMap(q,originalRing->cf,startingRing->cf);
     59
     60  /* map J into the new polynomial ring */
     61  int k = idSize(J);
     62  startingIdeal = idInit(k);
     63  for (int i=0; i<k; i++)
     64    startingIdeal->m[i] = p_PermPoly(J->m[i],NULL,originalRing,startingRing,nMap,NULL,0);
     65
     66  /* compute the dimension of the ideal in the original ring */
     67  dimensionOfIdeal = dim(J,s);
     68
     69  /* set the flag that convex computations only occur in the lower half space to true */
     70  onlyLowerHalfSpace = true;
     71
     72  /* set the function pointers to the appropiate functions */
     73  weightAdjustingAlgorithm1 = valued_adjustWeightForHomogeneity;
     74  weightAdjustingAlgorithm2 = valued_adjustWeightUnderHomogeneity;
     75  reductionAlgorithm = ppreduceInitially;
     76}
     77
     78static bool doNothing(ideal /*I*/, ring /*r*/, number /*p*/)
    2079{
    2180  return false;
    2281}
    2382
    24 tropicalStrategy initializeNonValuedCase()
     83/***
     84 * Initializes all relevant structures and information for the valued case,
     85 * i.e. computing a tropical variety without any valuation.
     86 **/
     87tropicalStrategy::tropicalStrategy(ideal I, ring r)
    2588{
    26   tropicalStrategy currentCase;
    27   currentCase.tropicalVarietyOfPolynomial = fullTropicalVarietyOfPolynomial;
    28   currentCase.adjustWeightForHomogeneity = nonvalued_adjustWeightForHomogeneity;
    29   currentCase.adjustWeightUnderHomogeneity = nonvalued_adjustWeightUnderHomogeneity;
    30   currentCase.reduce = doNothing;
    31   return currentCase;
     89  /* assume that the ground field of the originalRing is Q */
     90  assume(rField_is_Q(r));
     91  originalRing = rCopy(r);
     92
     93  /* the starting ring is the originaRing */
     94  startingRing = rCopy(originalRing);
     95
     96  /* the uniformizing parameter is non-existant*/
     97  uniformizingParameter = NULL;
     98
     99  /* set the startingIdeal */
     100  startingIdeal = id_Copy(I,startingRing);
     101
     102  /* compute the dimension of the ideal */
     103  dimensionOfIdeal = dim(startingIdeal,startingRing);
     104
     105  /* set the flag that convex computations only occur in the lower half space to false */
     106  onlyLowerHalfSpace = false;
     107
     108  /* set the function pointers to the appropiate functions */
     109  weightAdjustingAlgorithm1 = nonvalued_adjustWeightForHomogeneity;
     110  weightAdjustingAlgorithm2 = nonvalued_adjustWeightUnderHomogeneity;
     111  reductionAlgorithm = doNothing;
    32112}
     113
     114tropicalStrategy::tropicalStrategy(const tropicalStrategy& currentStrategy):
     115  originalRing(rCopy(currentStrategy.getOriginalRing())),
     116  startingRing(rCopy(currentStrategy.getStartingRing())),
     117  uniformizingParameter(n_Copy(currentStrategy.getUniformizingParameter(),startingRing->cf)),
     118  startingIdeal(id_Copy(currentStrategy.getStartingIdeal(),startingRing)),
     119  dimensionOfIdeal(currentStrategy.getDimensionOfIdeal()),
     120  onlyLowerHalfSpace(currentStrategy.restrictToLowerHalfSpace()),
     121  weightAdjustingAlgorithm1(currentStrategy.getWeightAdjustingAlgorithm1()),
     122  weightAdjustingAlgorithm2(currentStrategy.getWeightAdjustingAlgorithm2()),
     123  reductionAlgorithm(currentStrategy.getReductionAlgorithm())
     124{
     125  if (startingRing) rTest(startingRing);
     126  if (startingIdeal) id_Test(startingIdeal,startingRing);
     127}
     128
     129tropicalStrategy::~tropicalStrategy()
     130{
     131  if (startingRing) rTest(startingRing);
     132  if (startingIdeal) id_Test(startingIdeal,startingRing);
     133
     134  id_Delete(&startingIdeal,startingRing);
     135  n_Delete(&uniformizingParameter,startingRing->cf);
     136  rDelete(originalRing);
     137  rDelete(startingRing);
     138}
     139
     140tropicalStrategy& tropicalStrategy::operator=(const tropicalStrategy& currentStrategy)
     141{
     142  originalRing = rCopy(currentStrategy.getOriginalRing());
     143  startingRing = rCopy(currentStrategy.getStartingRing());
     144  uniformizingParameter = n_Copy(currentStrategy.getUniformizingParameter(),startingRing->cf);
     145  startingIdeal = id_Copy(currentStrategy.getStartingIdeal(),startingRing);
     146  dimensionOfIdeal = currentStrategy.getDimensionOfIdeal();
     147  onlyLowerHalfSpace = currentStrategy.restrictToLowerHalfSpace();
     148  weightAdjustingAlgorithm1 = currentStrategy.getWeightAdjustingAlgorithm1();
     149  weightAdjustingAlgorithm2 = currentStrategy.getWeightAdjustingAlgorithm2();
     150  reductionAlgorithm = currentStrategy.getReductionAlgorithm();
     151
     152  if (startingRing) rTest(startingRing);
     153  if (startingIdeal) id_Test(startingIdeal,startingRing);
     154
     155  return *this;
     156}
  • Singular/dyn_modules/gfanlib/tropicalStrategy.h

    r7aa26b2 r3c0aa5  
    77#include <set>
    88
     9typedef gfan::ZVector (*wAdjAlg1)(gfan::ZVector);
     10typedef gfan::ZVector (*wAdjAlg2)(gfan::ZVector,gfan::ZVector);
     11typedef bool (*redAlg)(ideal,ring,number);
     12
    913class tropicalStrategy
    1014{
     15private:
     16  ring originalRing;
     17  ring startingRing;
     18  number uniformizingParameter;
     19  ideal startingIdeal;
     20  int dimensionOfIdeal;
     21  bool onlyLowerHalfSpace;
     22
     23  gfan::ZVector (*weightAdjustingAlgorithm1) (gfan::ZVector w);
     24  gfan::ZVector (*weightAdjustingAlgorithm2) (gfan::ZVector v, gfan::ZVector w);
     25  bool (*reductionAlgorithm) (ideal I, ring r, number p);
     26
    1127public:
    12   std::set<gfan::ZCone> (*tropicalVarietyOfPolynomial)(const poly &g, const ring &r);
    13   gfan::ZVector (*adjustWeightForHomogeneity)(gfan::ZVector w);
    14   gfan::ZVector (*adjustWeightUnderHomogeneity)(gfan::ZVector v, gfan::ZVector w);
    15   bool (*reduce)(ideal I, ring r);
     28
     29  tropicalStrategy(const ideal I, const ring r);
     30  tropicalStrategy(const ideal J, const number p, const ring s);
     31  tropicalStrategy(const tropicalStrategy& currentStrategy);
     32  ~tropicalStrategy();
     33  tropicalStrategy& operator=(const tropicalStrategy& currentStrategy);
     34
     35
     36  ring getOriginalRing() const
     37  {
     38    rTest(originalRing);
     39    return originalRing;
     40  }
     41
     42  ring getStartingRing() const
     43  {
     44    rTest(startingRing);
     45    return startingRing;
     46  }
     47
     48  number getUniformizingParameter() const
     49  {
     50    return uniformizingParameter;
     51  }
     52
     53  ideal getStartingIdeal() const
     54  {
     55    id_Test(startingIdeal,startingRing);
     56    return startingIdeal;
     57  }
     58
     59  int getDimensionOfIdeal() const
     60  {
     61    return dimensionOfIdeal;
     62  }
     63
     64  bool restrictToLowerHalfSpace() const
     65  {
     66    return onlyLowerHalfSpace;
     67  }
     68
     69  wAdjAlg1 getWeightAdjustingAlgorithm1() const
     70  {
     71    return weightAdjustingAlgorithm1;
     72  }
     73
     74  wAdjAlg2 getWeightAdjustingAlgorithm2() const
     75  {
     76    return weightAdjustingAlgorithm2;
     77  }
     78
     79  redAlg getReductionAlgorithm() const
     80  {
     81    return reductionAlgorithm;
     82  }
     83
     84  gfan::ZVector adjustWeightForHomogeneity(gfan::ZVector w) const
     85  {
     86    return this->weightAdjustingAlgorithm1(w);
     87  }
     88
     89  gfan::ZVector adjustWeightUnderHomogeneity(gfan::ZVector v, gfan::ZVector w) const
     90  {
     91    return this->weightAdjustingAlgorithm2(v,w);
     92  }
     93
     94  bool reduce(ideal I, ring r) const
     95  {
     96    rTest(r);  id_Test(I,r);
     97    return this->reductionAlgorithm(I,r,uniformizingParameter);
     98  }
    1699};
    17100
    18 extern tropicalStrategy valuedCase;
    19 extern tropicalStrategy nonValuedCase;
    20 
    21 tropicalStrategy initializeValuedCase();
    22 tropicalStrategy initializeNonValuedCase();
    23 
    24101#endif
  • Singular/dyn_modules/gfanlib/tropicalVarietyOfPolynomials.cc

    r7aa26b2 r3c0aa5  
    33#include <libpolys/polys/monomials/p_polys.h>
    44#include <callgfanlib_conversion.h>
     5
     6#include <tropicalStrategy.h>
    57
    68#ifndef NDEBUG
     
    1012
    1113/***
    12  * Returns the tropical variety of a polynomial g in ring r
     14 * Returns the tropical variety of a polynomial g in ring r.
     15 * Depending on the current case, it will compute either the whole tropical variety
     16 * or just the part in the first lower half-space.
     17 * An empty set means that the tropical Variety is the whole space,
     18 * i.e. g is either 0 or a non-zero monomial.
    1319 **/
    14 std::set<gfan::ZCone> fullTropicalVarietyOfPolynomial(const poly &g, const ring &r)
     20std::set<gfan::ZCone> tropicalVariety(const poly g, const ring r, const tropicalStrategy &currentCase)
    1521{
    16   int n = r->N;
    17   std::set<gfan::ZCone> setOfCones;
     22  int n = rVar(r);
     23  std::set<gfan::ZCone> tropVar;
    1824  if (g && g->next)
    1925  {
     
    2935
    3036    int l = exponents.getHeight();
     37    gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
     38    lowerHalfSpaceCondition[0] = -1;
    3139    for (int i=0; i<l; i++)
    3240    {
     
    3644        equation.appendRow(exponents[i]-exponents[j]);
    3745        gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
     46        if (currentCase.restrictToLowerHalfSpace())
     47          inequalities.appendRow(lowerHalfSpaceCondition);
    3848        for (int k=0; k<l; k++)
    3949          if (k!=i && k!=j) inequalities.appendRow(exponents[i]-exponents[k]);
     
    4252        {
    4353          zc.canonicalize();
    44           setOfCones.insert(zc);
     54          tropVar.insert(zc);
    4555        }
    4656      }
    4757    }
    4858  }
    49   return setOfCones;
     59  return tropVar;
    5060}
    5161
    52 /***
    53  * Returns the tropical variety of a polynomial g in ring r
    54  * intersected with the lower halfspace
    55  **/
    56 std::set<gfan::ZCone> lowerTropicalVarietyOfPolynomial(const poly &g, const ring &r)
    57 {
    58   int n = r->N;
    59   std::set<gfan::ZCone> setOfCones;
    60   if (g && g->next)
    61   {
    62     int* expv = (int*) omAlloc((n+1)*sizeof(int));
    63     gfan::ZMatrix exponents = gfan::ZMatrix(0,n);
    64     for (poly s=g; s; pIter(s))
    65     {
    66       p_GetExpV(s,expv,r);
    67       gfan::ZVector zv = intStar2ZVector(n,expv);
    68       exponents.appendRow(intStar2ZVector(n,expv));
    69     }
    70     omFreeSize(expv,(n+1)*sizeof(int));
    7162
    72     int l = exponents.getHeight();
    73     for (int i=0; i<l; i++)
    74     {
    75       for (int j=i+1; j<l; j++)
    76       {
    77         gfan::ZMatrix equation = gfan::ZMatrix(0,n);
    78         equation.appendRow(exponents[i]-exponents[j]);
    79         gfan::ZMatrix inequalities = gfan::ZMatrix(1,n);
    80         inequalities[0][0] = -1;
    81         for (int k=0; k<l; k++)
    82           if (k!=i && k!=j) inequalities.appendRow(exponents[i]-exponents[k]);
    83         gfan::ZCone zc = gfan::ZCone(inequalities,equation);
    84         if (zc.dimension()>=n-1)
    85         {
    86           zc.canonicalize();
    87           setOfCones.insert(zc);
    88         }
    89       }
    90     }
    91   }
    92   return setOfCones;
    93 }
    94 
    95 #ifndef NDEBUG
    96 BOOLEAN tropicalVariety00(leftv res, leftv args)
    97 {
    98   leftv u = args;
    99   poly g = (poly) u->CopyD();
    100   omUpdateInfo();
    101   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
    102   std::set<gfan::ZCone> C = fullTropicalVarietyOfPolynomial(g,currRing);
    103   p_Delete(&g,currRing);
    104   res->rtyp = fanID;
    105   res->data = (char*) toFanStar(C);
    106   return FALSE;
    107 }
    108 BOOLEAN tropicalVariety01(leftv res, leftv args)
    109 {
    110   leftv u = args;
    111   poly g = (poly) u->CopyD();
    112   omUpdateInfo();
    113   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
    114   std::set<gfan::ZCone> C = lowerTropicalVarietyOfPolynomial(g,currRing);
    115   p_Delete(&g,currRing);
    116   res->rtyp = fanID;
    117   res->data = (char*) toFanStar(C);
    118   return FALSE;
    119 }
    120 #endif
     63// #ifndef NDEBUG
     64// BOOLEAN tropicalVariety00(leftv res, leftv args)
     65// {
     66//   leftv u = args;
     67//   poly g = (poly) u->CopyD();
     68//   omUpdateInfo();
     69//   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
     70//   std::set<gfan::ZCone> C = fullTropicalVarietyOfPolynomial(g,currRing);
     71//   p_Delete(&g,currRing);
     72//   res->rtyp = fanID;
     73//   res->data = (char*) toFanStar(C);
     74//   return FALSE;
     75// }
     76// BOOLEAN tropicalVariety01(leftv res, leftv args)
     77// {
     78//   leftv u = args;
     79//   poly g = (poly) u->CopyD();
     80//   omUpdateInfo();
     81//   Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
     82//   std::set<gfan::ZCone> C = lowerTropicalVarietyOfPolynomial(g,currRing);
     83//   p_Delete(&g,currRing);
     84//   res->rtyp = fanID;
     85//   res->data = (char*) toFanStar(C);
     86//   return FALSE;
     87// }
     88// #endif
  • Singular/dyn_modules/gfanlib/tropicalVarietyOfPolynomials.h

    r7aa26b2 r3c0aa5  
    66#include <set>
    77
     8#include <tropicalStrategy.h>
     9
    810#ifndef NDEBUG
    911#include <Singular/subexpr.h> // for leftv
     
    1113#endif
    1214
    13 std::set<gfan::ZCone> fullTropicalVarietyOfPolynomial(const poly &g, const ring &r);
    14 std::set<gfan::ZCone> lowerTropicalVarietyOfPolynomial(const poly &g, const ring &r);
     15std::set<gfan::ZCone> tropicalVariety(const poly g, const ring r, const tropicalStrategy &currentCase);
    1516
    1617#ifndef NDEBUG
  • Singular/dyn_modules/gfanlib/witness.cc

    r7aa26b2 r3c0aa5  
    1717{
    1818  ring origin = currRing;
    19   if (currRing != r) rChangeCurrRing(r);
     19  rComplete(r);
     20  if (origin != r) rChangeCurrRing(r);
    2021  ideal F = idInit(1); F->m[0]=f;
    21   ideal Rest; matrix U;
    22   ideal m = idLift(G,F,&Rest,FALSE,TRUE,TRUE,&U);
    23   id_Delete((ideal*)&U, r);
     22  ideal m = idLift(G,F,NULL,FALSE,FALSE);
    2423  F->m[0]=NULL; id_Delete(&F, currRing);
    2524  matrix Q = id_Module2formatedMatrix(m,IDELEMS(G),1,currRing);
    26   assume(Rest->m[0] == NULL);
    27   id_Delete((ideal*) &Rest, r);
    28   if (currRing != origin) rChangeCurrRing(origin);
     25  if (origin != r) rChangeCurrRing(origin);
    2926  return Q;
    3027}
     
    4946
    5047/***
    51  * Computes a witness g of the w-homogeneous m in ideal I in the ring r,
    52  * such that in_w(g)=m, where w is the uppermost weight vector of r.
     48 * Let w be the uppermost weight vector in the matrix defining the ordering on r.
     49 * Let I be a Groebner basis of an ideal in r, inI its initial form with respect w.
     50 * Given an w-homogeneous element m of inI, computes a witness g of m in I,
     51 * i.e. g in I such that in_w(g)=m.
    5352 **/
    5453poly witness(const poly m, const ideal I, const ideal inI, const ring r)
     
    8988}
    9089#endif
    91 
    92 
    93 /***
    94  * Suppose r and s the same ring but with two adjacent orderings,
    95  * w a weight vector in the relative interior of their common facet.
    96  * Given a standard basis of an ideal I with respect to the ordering of r,
    97  * a standard basis of its w-initial ideal inI with respect to
    98  * the ordering of s, computes a standard basis of I with respect to
    99  * the ordering of s.
    100  **/
    101 ideal lift(const ideal I, const ring r, const ideal inI, const ring s)
    102 {
    103   nMapFunc identity = n_SetMap(r->cf,s->cf);
    104   int k = idSize(I); ideal Is = idInit(k);
    105   for (int i=0; i<k; i++)
    106     Is->m[i] = p_PermPoly(I->m[i],NULL,r,s,identity,NULL,0);
    107   ideal J = idInit(k);
    108   for (int i=0; i<k; i++)
    109     J->m[i] = witness(inI->m[i],Is,inI,s);
    110   id_Delete(&Is,s);
    111   return J;
    112 }
    113 
    114 void deleteOrdering(ring r)
    115 {
    116   if (r->order != NULL)
    117   {
    118     int i=rBlocks(r);
    119     assume(r->block0 != NULL && r->block1 != NULL && r->wvhdl != NULL);
    120     // delete order
    121     omFreeSize((ADDRESS)r->order,i*sizeof(int));
    122     omFreeSize((ADDRESS)r->block0,i*sizeof(int));
    123     omFreeSize((ADDRESS)r->block1,i*sizeof(int));
    124     // delete weights
    125     for (int j=0; j<i; j++)
    126     {
    127       if (r->wvhdl[j]!=NULL)
    128         omFree(r->wvhdl[j]);
    129     }
    130     omFreeSize((ADDRESS)r->wvhdl,i*sizeof(int *));
    131   }
    132   else
    133     assume(r->block0 == NULL && r->block1 == NULL && r->wvhdl == NULL);
    134   return;
    135 }
    136 
    137 std::pair<ideal,ring> flip(const ideal I, const ring r, const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal, const tropicalStrategy currentCase)
    138 {
    139   gfan::ZVector (*adjustWeight0)(gfan::ZVector w);
    140   adjustWeight0 = currentCase.adjustWeightForHomogeneity;
    141   gfan::ZVector (*adjustWeight1)(gfan::ZVector v, gfan::ZVector w);
    142   adjustWeight1 = currentCase.adjustWeightUnderHomogeneity;
    143 
    144   gfan::ZVector adjustedInteriorPoint = adjustWeight0(interiorPoint);
    145   gfan::ZVector adjustedFacetNormal = adjustWeight1(facetNormal,adjustedInteriorPoint);
    146 
    147   ring origin = currRing;
    148   ideal inIr = initial(I,r,interiorPoint);
    149 
    150   bool ok;
    151   ring s = rCopy0(r);
    152   int n = rVar(s);
    153   deleteOrdering(s);
    154   s->order = (int*) omAlloc0(4*sizeof(int));
    155   s->block0 = (int*) omAlloc0(4*sizeof(int));
    156   s->block1 = (int*) omAlloc0(4*sizeof(int));
    157   s->wvhdl = (int**) omAlloc0(4*sizeof(int));
    158   s->order[0] = ringorder_a;
    159   s->block0[0] = 1;
    160   s->block1[0] = n;
    161   s->wvhdl[0] = ZVectorToIntStar(adjustedInteriorPoint,ok);
    162   s->order[1] = ringorder_wp;
    163   s->block0[1] = 1;
    164   s->block1[1] = n;
    165   s->wvhdl[1] = ZVectorToIntStar(adjustedFacetNormal,ok);
    166   s->order[2] = ringorder_C;
    167   rComplete(s,1);
    168   rChangeCurrRing(s);
    169   nMapFunc identity = n_SetMap(r->cf,currRing->cf);
    170 
    171   int k = idSize(I); ideal inIs = idInit(k);
    172   for (int i=0; i<k; i++)
    173     inIs->m[i] = p_PermPoly(inIr->m[i],NULL,r,s,identity,NULL,0);
    174   intvec* nullVector = NULL;
    175   ideal inIsGB = kStd(inIs,currQuotient,testHomog,&nullVector);
    176   ideal IsGB = lift(I,r,inIsGB,s);
    177 
    178   rChangeCurrRing(origin);
    179   id_Delete(&inIr,r);
    180   id_Delete(&inIs,s);
    181   id_Delete(&inIsGB,s);
    182   return std::make_pair(IsGB,s);
    183 }
  • Singular/dyn_modules/gfanlib/witness.h

    r7aa26b2 r3c0aa5  
    22#define WITNESS_H
    33
    4 #include <utility>
    5 #include <kernel/polys.h>
    6 #include <Singular/ipid.h>
    7 #include <gfanlib/gfanlib_vector.h>
    8 #include <tropicalStrategy.h>
     4#include <libpolys/polys/monomials/monomials.h>
     5#include <libpolys/polys/simpleideals.h>
    96
    107poly witness(const poly m, const ideal I, const ideal inI, const ring r);
    11 std::pair<ideal,ring> flip(const ideal I, const ring r, const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal, const tropicalStrategy currentCase);
    128
    139#ifndef NDEBUG
     10#include <Singular/ipid.h>
     11
    1412BOOLEAN dwr0(leftv res, leftv args);
    1513BOOLEAN witness0(leftv res, leftv args);
  • dyn_modules/callgfanlib/tropicalVariety.cc

    r7aa26b2 r3c0aa5  
    11#include <callgfanlib_conversion.h>
     2#include <bbfan.h>
    23#include <groebnerCone.h>
    3 #include <tropicalVarietyOfPolynomials.h>
    44#include <tropicalVarietyOfIdeals.h>
    55#include <libpolys/coeffs/numbers.h>
    66#include <kernel/structs.h>
    77
    8 static poly adjustPoly(const poly g, const number p, const ring r)
    9 {
    10   poly currentTerm = g;
    11   poly fBeginning = p_Head(currentTerm,r);
    12   poly f = fBeginning;
    13   number coeff0 = p_GetCoeff(currentTerm,r);
    14   unsigned power = 0;
    15   while (n_DivBy(coeff0,p,r->cf))
    16   {
    17     power++;
    18     number coeff1 = n_Div(coeff0,p,r->cf);
    19     n_Delete(&coeff0,r->cf);
    20     coeff0 = coeff1;
    21   }
    22   p_SetCoeff(f,coeff0,r);
    23   p_AddExp(f,1,power,r); p_Setm(f,r);
    24 
    25   for (pIter(currentTerm); currentTerm; pIter(currentTerm))
    26   {
    27     f->next = p_Head(currentTerm,r);
    28     pIter(f);
    29     coeff0 = p_GetCoeff(currentTerm,r);
    30     power = 0;
    31     while (n_DivBy(coeff0,p,r->cf))
    32     {
    33       power++;
    34       number coeff1 = n_Div(coeff0,p,r->cf);
    35       n_Delete(&coeff0,r->cf);
    36       coeff0 = coeff1;
    37     }
    38     p_SetCoeff(f,coeff0,r);
    39     p_AddExp(f,1,power,r); p_Setm(f,r);
    40   }
    41   return fBeginning;
    42 }
    43 
    44 static ideal adjustIdeal(const ideal I, const number p, const ring r)
    45 {
    46   int k = idSize(I);
    47   ideal J = idInit(k+1);
    48   for (int i = 0; i<k; i++)
    49     J->m[i+1] = adjustPoly(I->m[i],p,r);
    50   poly h0 = p_One(r);
    51   p_SetCoeff(h0,n_Copy(p,r->cf),r);
    52   poly h1 = p_One(r);
    53   p_SetExp(h1,1,1,r);
    54   p_Setm(h1,r);
    55   J->m[0] = p_Add_q(h1,p_Neg(h0,r),r);
    56   return J;
    57 }
    58 
    598BOOLEAN tropicalVariety(leftv res, leftv args)
    609{
    61   ring origin = currRing;
    62   ring r = rCopy(currRing);
    63   rChangeCurrRing(r);
    6410  leftv u = args;
    65   if ((u!=NULL) && (u->Typ()==POLY_CMD))
    66   {
    67     poly g = (poly) u->Data();
    68     leftv v = u->next;
    69     if (v==NULL)
    70     {
    71       std::set<gfan::ZCone> setOfCones = fullTropicalVarietyOfPolynomial(g,currRing);
    72       gfan::ZFan* tropG = toFanStar(setOfCones);
    73       res->rtyp = fanID;
    74       res->data = (char*) tropG;
    75       rChangeCurrRing(origin);
    76       return FALSE;
    77     }
    78     if ((v!=NULL) && (v->Typ()==NUMBER_CMD))
    79     {
    80       number p = (number) v->CopyD();
    81       poly f = adjustPoly(g,p,currRing);
    82       std::set<gfan::ZCone> setOfCones = lowerTropicalVarietyOfPolynomial(f,currRing);
    83       gfan::ZFan* tropG = toFanStar(setOfCones);
    84       res->rtyp = fanID;
    85       res->data = (char*) tropG;
    86       rChangeCurrRing(origin);
    87       return FALSE;
    88     }
    89     if ((v!=NULL) && (v->Typ()==POLY_CMD))
    90     {
    91       std::set<gfan::ZCone> setOfCones = lowerTropicalVarietyOfPolynomial(g,currRing);
    92       gfan::ZFan* tropG = toFanStar(setOfCones);
    93       res->rtyp = fanID;
    94       res->data = (char*) tropG;
    95       rChangeCurrRing(origin);
    96       return FALSE;
    97     }
    98   }
    9911  if ((u!=NULL) && (u->Typ()==IDEAL_CMD))
    10012  {
     
    10315    if (v==NULL)
    10416    {
    105       setOfGroebnerConeData TropI = fullTropicalVariety(I,currRing);
    106       gfan::ZFan* tropI = toFanStar(TropI);
     17      tropicalStrategy currentStrategy(I,currRing);
     18      gfan::ZFan* tropI = tropicalVariety(I,currRing,currentStrategy);
    10719      res->rtyp = fanID;
    10820      res->data = (char*) tropI;
    109       rChangeCurrRing(origin);
    11021      return FALSE;
    11122    }
     
    11324    {
    11425      number p = (number) v->CopyD();
    115       ideal J = adjustIdeal(I,p,currRing);
    116       setOfGroebnerConeData TropI = lowerTropicalVariety(J,currRing);
    117       gfan::ZFan* tropI = toFanStar(TropI);
     26      tropicalStrategy currentStrategy(I,p,currRing);
     27      gfan::ZFan* tropI = tropicalVariety(I,currRing,currentStrategy);
    11828      res->rtyp = fanID;
    11929      res->data = (char*) tropI;
    120       rChangeCurrRing(origin);
    12130      return FALSE;
    12231    }
    123     if ((v!=NULL) && (v->Typ()==POLY_CMD))
    124     {
    125       setOfGroebnerConeData TropI = lowerTropicalVariety(I,currRing);
    126       gfan::ZFan* tropI = toFanStar(TropI);
    127       res->rtyp = fanID;
    128       res->data = (char*) tropI;
    129       rChangeCurrRing(origin);
    130       return FALSE;
    131     }
    132     rChangeCurrRing(origin);
    13332    return FALSE;
    13433  }
    135   rChangeCurrRing(origin);
    13634  WerrorS("tropicalVariety: unexpected parameters");
    13735  return TRUE;
  • dyn_modules/callgfanlib/tropicalVarietyOfIdeals.cc

    r7aa26b2 r3c0aa5  
     1#include <groebnerCone.h>
    12#include <startingCone.h>
    2 #include <neighbours.h>
     3#include <tropicalTraversal.h>
    34
    4 setOfGroebnerConeData fullTropicalVariety(const ideal I, const ring r)
     5gfan::ZFan* tropicalVariety(const ideal I, const ring r, const tropicalStrategy currentStrategy)
    56{
    6   groebnerConeData currentCone = groebnerConeData(I,r,tropicalStartingCone(I,r,nonValuedCase));
    7   setOfGroebnerConeData tropicalVariety;
    8   setOfGroebnerConeData workingList;
    9   workingList.insert(currentCone);
    10   while (!workingList.empty())
    11   {
    12     currentCone = *workingList.begin();
    13     tropicalVariety.insert(currentCone);
    14     workingList.erase(workingList.begin());
    15 
    16     setOfGroebnerConeData neighbours = tropicalNeighbours(currentCone,nonValuedCase);
    17     for (setOfGroebnerConeData::iterator it=neighbours.begin(); it!=neighbours.end(); it++)
    18       if (!tropicalVariety.count(*it)) workingList.insert(*it);
    19   }
    20   return tropicalVariety;
     7  groebnerCone startingCone = tropicalStartingCone(I,r,currentStrategy);
     8  groebnerCones tropicalVariety = tropicalTraversal(startingCone);
     9  return toFanStar(tropicalVariety);
    2110}
    22 
    23 setOfGroebnerConeData lowerTropicalVariety(const ideal I, const ring r)
    24 {
    25   groebnerConeData currentCone = groebnerConeData(I,r,tropicalStartingCone(I,r,valuedCase));
    26   setOfGroebnerConeData tropicalVariety;
    27   setOfGroebnerConeData workingList;
    28   workingList.insert(currentCone);
    29   while (!workingList.empty())
    30   {
    31     currentCone = *workingList.begin();
    32     tropicalVariety.insert(currentCone);
    33     workingList.erase(workingList.begin());
    34 
    35     setOfGroebnerConeData neighbours = tropicalNeighbours(currentCone,valuedCase);
    36     for (setOfGroebnerConeData::iterator it=neighbours.begin(); it!=neighbours.end(); it++)
    37       if (!tropicalVariety.count(*it)) workingList.insert(*it);
    38   }
    39   return tropicalVariety;
    40 }
  • dyn_modules/callgfanlib/tropicalVarietyOfIdeals.h

    r7aa26b2 r3c0aa5  
    22#define CALLGfANLIB_TROPICALVARIETYOFIDEALS
    33
    4 #include <groebnerCone.h>
     4#include <libpolys/polys/simpleideals.h>
     5#include <tropicalStrategy.h>
    56
    6 setOfGroebnerConeData fullTropicalVariety(const ideal I, const ring r);
    7 setOfGroebnerConeData lowerTropicalVariety(const ideal I, const ring r);
     7gfan::ZFan* tropicalVariety(const ideal I, const ring r, const tropicalStrategy currentStrategy);
    88
    99#endif
  • gfanlib/gfanlib_matrix.h

    r7aa26b2 r3c0aa5  
    99#define LIB_ZMATRIX_H_
    1010
     11#include <sstream>
    1112#include <vector>
    1213#include <algorithm>
     
    3536  Matrix():width(0),height(0){
    3637  };
    37   static Matrix rowVectorMatrix(Vector<typ> const &v)
     38  Matrix rowVectorMatrix(Vector<typ> const &v)
    3839  {
    3940    Matrix ret(1,v.size());
  • gfanlib/gfanlib_vector.h

    r7aa26b2 r3c0aa5  
    310310    assert(!typ::isField());
    311311    return (*this)/gcd();
     312  }
     313
     314  void debugPrint()const
     315  {
     316    std::stringstream s;
     317    s<<"(";
     318    for(typename std::vector<typ>::const_iterator i=this->v.begin();i!=this->v.end();i++)
     319    {
     320      if(i!=this->v.begin()) s<<",";
     321      s<<*i;
     322    }
     323    s<<")"<<std::endl;
     324    std::cout << s.str();
     325    return;
    312326  }
    313327};
  • gfanlib/gfanlib_zcone.h

    r7aa26b2 r3c0aa5  
    1010
    1111#include "gfanlib_matrix.h"
     12#include <sstream>
    1213
    1314namespace gfan{
     
    360361
    361362
    362 
    363 
    364363#endif /* LIB_CONE_H_ */
  • kernel/ideals.cc

    r7aa26b2 r3c0aa5  
    26592659  omFreeSize((ADDRESS)(id_sort), idsize*sizeof(poly_sort));
    26602660}
     2661
     2662void id_Write(ideal id, ring r)
     2663{
     2664  for (int i=0; i<idSize(id); i++)
     2665    p_Write(id->m[i], r, r);
     2666}
  • kernel/ideals.h

    r7aa26b2 r3c0aa5  
    209209ideal id_TensorModuleMult(const int m, const ideal M, const ring rRing); // image of certain map for BGG
    210210
     211void id_Write(ideal id, ring r);
    211212
    212213#endif
Note: See TracChangeset for help on using the changeset viewer.