Changeset b71400a in git


Ignore:
Timestamp:
Sep 15, 2014, 11:08:21 AM (10 years ago)
Author:
Yue Ren <ren@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
dffd1541cc74564f509bd5d284948852484430b0
Parents:
9c03260ac0404405e4994f73b30ef46f079535d8
git-author:
Yue Ren <ren@mathematik.uni-kl.de>2014-09-15 11:08:21+02:00
git-committer:
Yue Ren <ren@mathematik.uni-kl.de>2015-02-06 13:47:05+01:00
Message:
chg: status update 15.09.
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • Singular/dyn_modules/gfanlib/adjustWeights.cc

    r9c0326 rb71400a  
    1818static bool checkForNonPositiveLaterEntries(const gfan::ZVector w)
    1919{
    20   if (w[0].sign()<0)
    21   {
    22     std::cout << "ERROR: negative weight in weight vector first entry" << std::endl
    23               << "weight: " << w << std::endl;
    24     return false;
    25   }
    26   for (unsigned i=1; i<w.size(); i++)
     20  // if (w[0].sign()<0)
     21  // {
     22  //   std::cout << "ERROR: negative weight in weight vector first entry" << std::endl
     23  //             << "weight: " << w << std::endl;
     24  //   return false;
     25  // }
     26  for (unsigned i=0; i<w.size(); i++)
    2727  {
    2828    if (w[i].sign()<=0)
     
    7171  for (unsigned i=1; i<w.size(); i++)
    7272    v[i]=-w[i]+max+1;
    73   assume(checkForNonPositiveLaterEntries(v));
     73  // assume(checkForNonPositiveLaterEntries(v));
    7474  return v;
    7575}
     
    108108gfan::ZVector valued_adjustWeightUnderHomogeneity(const gfan::ZVector e, const gfan::ZVector w)
    109109{
    110   assume(checkForNonPositiveLaterEntries(w));
     110  // assume(checkForNonPositiveLaterEntries(w));
    111111  /* find k such that e+k*w is strictly positive,
    112112   * i.e. k such that e[i]+k*w[i] is strictly positive
  • Singular/dyn_modules/gfanlib/callgfanlib_conversion.cc

    r9c0326 rb71400a  
    9696}
    9797
     98gfan::ZVector expvToZVector(const int n, const int* expv)
     99{
     100  gfan::ZVector zv(n);
     101  for(int i=0; i<n; i++)
     102    zv[i]=gfan::Integer(expv[i+1]);
     103  return zv;
     104}
     105
    98106gfan::ZVector wvhdlEntryToZVector(const int n, const int* wvhdl0)
    99107{
  • Singular/dyn_modules/gfanlib/callgfanlib_conversion.h

    r9c0326 rb71400a  
    1919char* toString(gfan::ZMatrix const &m);
    2020
     21gfan::ZVector expvToZVector(const int n, const int* expv);
     22
    2123gfan::ZFan* toFanStar(std::set<gfan::ZCone> setOfCones);
    2224std::set<gfan::ZVector> rays(std::set<gfan::ZCone> setOfCones);
  • Singular/dyn_modules/gfanlib/groebnerCone.cc

    r9c0326 rb71400a  
    8181  }
    8282  if (r) polynomialRing = rCopy(r);
     83  currentCase.pReduce(polynomialIdeal,r);
    8384  currentCase.reduce(reducedPolynomialIdeal,polynomialRing);
    8485
     
    9697    {
    9798      p_GetExpV(g,leadexpv,r);
    98       leadexpw = intStar2ZVector(n, leadexpv);
     99      leadexpw = expvToZVector(n, leadexpv);
    99100      pIter(g);
    100101      while (g)
    101102      {
    102103        p_GetExpV(g,tailexpv,r);
    103         tailexpw = intStar2ZVector(n, tailexpv);
     104        tailexpw = expvToZVector(n, tailexpv);
    104105        inequalities.appendRow(leadexpw-tailexpw);
    105106        pIter(g);
     
    109110  omFreeSize(leadexpv,(n+1)*sizeof(int));
    110111  omFreeSize(tailexpv,(n+1)*sizeof(int));
    111   if (currentStrategy->restrictToLowerHalfSpace())
    112   {
    113     gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
    114     lowerHalfSpaceCondition[0] = -1;
    115     inequalities.appendRow(lowerHalfSpaceCondition);
    116   }
     112  // if (currentStrategy->restrictToLowerHalfSpace())
     113  // {
     114  //   gfan::ZVector lowerHalfSpaceCondition = gfan::ZVector(n);
     115  //   lowerHalfSpaceCondition[0] = -1;
     116  //   inequalities.appendRow(lowerHalfSpaceCondition);
     117  // }
    117118
    118119  polyhedralCone = gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
     
    144145  }
    145146  if (r) polynomialRing = rCopy(r);
     147  currentCase.pReduce(polynomialIdeal,r);
    146148  currentCase.reduce(reducedPolynomialIdeal,polynomialRing);
    147149
     
    205207  }
    206208  if (r) polynomialRing = rCopy(r);
     209  currentCase.pReduce(polynomialIdeal,r);
    207210  currentCase.reduce(reducedPolynomialIdeal,polynomialRing);
    208211
     
    259262  assume(checkPolynomialInput(I,r));
    260263  assume(checkPolynomialInput(inI,r));
     264  currentCase.pReduce(polynomialIdeal,r);
    261265  currentStrategy->reduce(reducedPolynomialIdeal,polynomialRing);
    262266  // currentStrategy->reduce(initialPolynomialIdeal,polynomialRing);
     
    435439   *   to obtain an initial form with respect to interiorPoint+e*facetNormal,
    436440   *   for e>0 sufficiently small */
    437   std::pair<ideal,ring> flipped = currentStrategy->flip(reducedPolynomialIdeal,polynomialRing,interiorPoint,facetNormal);
     441  std::pair<ideal,ring> flipped = currentStrategy->getFlip(polynomialIdeal,reducedPolynomialIdeal,polynomialRing,interiorPoint,facetNormal);
    438442  assume(checkPolynomialInput(flipped.first,flipped.second));
    439443  groebnerCone flippedCone(flipped.first, flipped.second, interiorPoint, facetNormal, *currentStrategy);
     
    547551  for (int i=1; i<r-1; i++)
    548552  {
    549     newInequalities = inequalities.submatrix(0,0,i-1,c);
     553    newInequalities = inequalities.submatrix(0,0,i,c);
    550554    newInequalities.append(inequalities.submatrix(i+1,0,r,c));
    551555    newEquations = equations;
     
    566570
    567571
     572bool groebnerCone::pointsOutwards(const gfan::ZVector w) const
     573{
     574  gfan::ZCone dual = polyhedralCone.dualCone();
     575  return (!dual.contains(w));
     576}
     577
     578
    568579/***
    569580 * Returns a complete list of neighboring Groebner cones in the tropical variety.
     
    578589    {
    579590      ideal initialIdeal = initial(polynomialIdeal,polynomialRing,interiorPoints[i]);
    580       std::set<gfan::ZVector> rays = raysOfTropicalStar(initialIdeal,polynomialRing,interiorPoints[i],*currentStrategy);
    581       for (std::set<gfan::ZVector>::iterator ray = rays.begin(); ray!=rays.end(); ray++)
    582         neighbours.insert(this->flipCone(interiorPoints[i],*ray));
     591      gfan::ZMatrix ray = raysOfTropicalStar(initialIdeal,polynomialRing,interiorPoints[i],*currentStrategy);
     592      for (int j=0; j<ray.getHeight(); j++)
     593        if (pointsOutwards(ray[j]))
     594        {
     595          groebnerCone neighbour = flipCone(interiorPoints[i],ray[j]);
     596          neighbours.insert(neighbour);
     597        }
    583598    }
    584599  }
  • Singular/dyn_modules/gfanlib/groebnerCone.h

    r9c0326 rb71400a  
    9494   **/
    9595  bool checkFlipConeInput(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
     96  bool pointsOutwards(const gfan::ZVector) const;
    9697  groebnerCone flipCone(const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
    9798
  • Singular/dyn_modules/gfanlib/initial.cc

    r9c0326 rb71400a  
    3636  for (int i=0; i<W.getHeight(); i++)
    3737    d[i] = wDeg(p,r,W[i]);
     38  return d;
     39}
     40
     41gfan::ZVector WDeg(const poly p, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
     42{
     43  gfan::ZVector d = gfan::ZVector(W.getHeight()+1);
     44  d[0] = wDeg(p,r,w);
     45  for (int i=0; i<W.getHeight(); i++)
     46    d[i+1] = wDeg(p,r,W[i]);
    3847  return d;
    3948}
     
    166175}
    167176
     177poly initial(const poly p, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
     178{
     179  int n = rVar(r);
     180  poly q0 = p_Head(p,r);
     181  poly q1 = q0;
     182  gfan::ZVector d = WDeg(p,r,w,W);
     183  for (poly currentTerm = p->next; currentTerm; pIter(currentTerm))
     184  {
     185    gfan::ZVector e = WDeg(currentTerm,r,w,W);
     186    if (d<e)
     187    {
     188      p_Delete(&q0,r);
     189      q0 = p_Head(p,r);
     190      q1 = q0;
     191      d = e;
     192    }
     193    else
     194      if (d==e)
     195      {
     196        pNext(q1) = p_Head(currentTerm,r);
     197        pIter(q1);
     198      }
     199  }
     200  return q0;
     201}
     202
     203
    168204/***
    169205 * Runs the above procedure over all generators of an ideal.
     
    179215}
    180216
     217ideal initial(const ideal I, const ring r, const gfan::ZVector w, const gfan::ZMatrix W)
     218{
     219  int k = idSize(I); ideal inI = idInit(k);
     220  for (int i=0; i<k; i++)
     221    inI->m[i] = initial(I->m[i],r,w,W);
     222  return inI;
     223}
    181224
    182225#ifndef NDEBUG
  • Singular/dyn_modules/gfanlib/initial.h

    r9c0326 rb71400a  
    4040ideal initial(const ideal I, const ring r, const gfan::ZVector w);
    4141poly initial(const poly p, const ring r, const gfan::ZMatrix W);
    42 poly initial(const ideal I, const ring r, const gfan::ZMatrix W);
     42ideal initial(const ideal I, const ring r, const gfan::ZMatrix W);
     43poly initial(const poly p, const ring r, const gfan::ZVector w, const gfan::ZMatrix W);
     44ideal initial(const ideal I, const ring r, const gfan::ZVector w, const gfan::ZMatrix W);
    4345
    4446
  • Singular/dyn_modules/gfanlib/ppinitialReduction.cc

    r9c0326 rb71400a  
    88#include <map>
    99#include <set>
     10
     11bool isOrderingLocalInT(const ring r)
     12{
     13  poly one = p_One(r);
     14  poly t = p_One(r);
     15  p_SetExp(t,1,1,r);
     16  p_Setm(t,r);
     17  int s = p_LmCmp(one,t,r);
     18  p_Delete(&one,r);
     19  p_Delete(&t,r);
     20  return (s==1);
     21}
    1022
    1123/***
     
    8092}
    8193
    82 
    8394#ifndef NDEBUG
    8495BOOLEAN pppReduce(leftv res, leftv args)
     
    105116}
    106117#endif //NDEBUG
     118
     119bool pReduce0(ideal &I, const number p, const ring r)
     120{
     121  int k = idSize(I);
     122  for (int i=0; i<k; i++)
     123  {
     124    if (I->m[i]!=NULL)
     125    {
     126      number c = p_GetCoeff(I->m[i],r);
     127      if (!n_Equal(p,c,r->cf))
     128        if (pReduce(I->m[i],p,r))
     129          return true;
     130    }
     131  }
     132  return false;
     133}
    107134
    108135
     
    494521
    495522
     523static std::vector<int> synchronize(const ideal I, const ideal Hi)
     524{
     525  int k = idSize(Hi);
     526  int l = idSize(I);
     527  std::vector<int> synch(k);
     528  int j;
     529  for (int i=0; i<k; i++)
     530  {
     531    for (j=0; j<l; j++)
     532    {
     533      if (I->m[i]==Hi->m[j])
     534        synch[i] = j;
     535    }
     536    if (j==l)
     537      synch[i] = -1;
     538  }
     539  return synch;
     540}
     541
     542static void synchronize(ideal I, const ideal Hi, const std::vector<int> synch)
     543{
     544  for (unsigned i=0; i<synch.size(); i++)
     545    if (synch[i]>=0)
     546      I->m[i] = Hi->m[synch[i]];
     547}
     548
    496549/**
    497550 * reduces I initially with respect to itself.
     
    525578
    526579  std::map<long,ideal>::iterator it=H.begin();
     580  it++;
    527581  ideal Hi = it->second;
    528   assume(idSize(Hi)==1);
    529   assume(pLength(Hi->m[0])==2);
    530   it++;
    531   Hi = it->second;
    532582
    533583  /***
     
    535585   *  and all lower components
    536586   **/
     587  std::vector<int> synch = synchronize(I,Hi);
    537588  if (ppreduceInitially(Hi,p,r)) return true;
     589  synchronize(I,Hi,synch);
     590  id_Test(Hi,r);
     591  id_Test(I,r);
    538592
    539593  ideal G = idInit(n); int m=0;
     
    577631    }
    578632    m += l; IDELEMS(GG) = m; GG->m = &G->m[n-m];
     633    std::vector<int> synch = synchronize(I,it->second);
    579634    if (ppreduceInitially(it->second,p,GG,r)) return true;
     635    synchronize(I,it->second,synch);
    580636    idShallowDelete(&Hi); Hi = it->second;
    581637  }
  • Singular/dyn_modules/gfanlib/ppinitialReduction.h

    r9c0326 rb71400a  
    1414#endif
    1515
    16 bool ppreduceInitially(ideal I, ring r, number p);
     16bool isOrderingLocalInT(const ring r);
     17bool pReduce0(ideal &I, const number p, const ring r);
     18bool ppreduceInitially(ideal I, const ring r, const number p);
    1719BOOLEAN ppreduceInitially(leftv res, leftv args);
    1820
  • Singular/dyn_modules/gfanlib/singularWishlist.h

    r9c0326 rb71400a  
    5959inline void idShallowDelete (ideal *h)
    6060{
    61   int j,elems;
    62   if (*h == NULL)
    63     return;
    64   elems=j=(*h)->nrows*(*h)->ncols;
    65   if (j>0)
    66     omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*elems);
    67   omFreeBin((ADDRESS)*h, sip_sideal_bin);
    68   *h=NULL;
    69 }
     61  if (*h != NULL)
     62  {
     63    int k;
     64    k=(*h)->nrows*(*h)->ncols;
     65    if (k>0)
     66      omFreeSize((ADDRESS)((*h)->m),sizeof(poly)*k);
     67    omFreeBin((ADDRESS)*h, sip_sideal_bin);
     68    *h=NULL;
     69  }
     70  return;
     71 }
    7072
    7173#endif
  • Singular/dyn_modules/gfanlib/startingCone.cc

    r9c0326 rb71400a  
    2121static bool checkContainmentInTropicalVariety(const groebnerCone sigma)
    2222{
    23   gfan::ZVector w = sigma.getInteriorPoint();
    2423  ideal I = sigma.getPolynomialIdeal();
    2524  ring r = sigma.getPolynomialRing();
    2625  const tropicalStrategy* currentStrategy = sigma.getTropicalStrategy();
    27   poly s = currentStrategy->checkInitialIdealForMonomial(I,r,w);
    28   if (s)
    29   {
    30     p_Delete(&s,r);
    31     return false;
    32   }
     26
     27  gfan::ZCone zc = sigma.getPolyhedralCone();
     28  gfan::ZMatrix zm = zc.extremeRays();
     29  for (int i=0; i<zm.getHeight(); i++)
     30  {
     31    gfan::ZVector w = zm[i];
     32    if (currentStrategy->isValuationNonTrivial() && w[0].sign()==0)
     33      continue;
     34    poly s = currentStrategy->checkInitialIdealForMonomial(I,r,w);
     35    if (s)
     36    {
     37      p_Delete(&s,r);
     38      return false;
     39    }
     40  }
     41
     42  zm = zc.generatorsOfLinealitySpace();
     43  for (int i=0; i<zm.getHeight(); i++)
     44  {
     45    gfan::ZVector w = zm[i];
     46    if (currentStrategy->isValuationNonTrivial() && w[0].sign()==0)
     47      continue;
     48    poly s = currentStrategy->checkInitialIdealForMonomial(I,r,w);
     49    if (s)
     50    {
     51      p_Delete(&s,r);
     52      return false;
     53    }
     54  }
     55
    3356  return true;
    3457}
     
    313336{
    314337  ring r = currentStrategy.getStartingRing();
    315   ideal I = currentStrategy.getStartingIdeal();
     338  const ideal I = currentStrategy.getStartingIdeal();
    316339  if (currentStrategy.isConstantCoefficientCase())
    317340  {
     
    359382    id_Delete(&J,s);
    360383
    361     assume(checkContainmentInTropicalVariety(startingCone));
     384    // assume(checkContainmentInTropicalVariety(startingCone));
    362385    return startingCone;
    363386  }
     
    463486    rDelete(s);
    464487
     488    assume(checkContainmentInTropicalVariety(startingCone));
    465489    return startingCone;
    466490  }
  • Singular/dyn_modules/gfanlib/tropicalCurves.cc

    r9c0326 rb71400a  
    142142
    143143    inIs = gfanlib_kStd_wrapper(inIs,s,isHomog);
    144     ideal ininIs = initial(inIs,s,W[W.getHeight()-1]);
    145 
    146     poly mons = checkForMonomialViaSuddenSaturation(inIs,s);
     144    ideal ininIs = initial(inIs,s,w,W);
     145
     146    poly mons = checkForMonomialViaSuddenSaturation(ininIs,s);
    147147    if (mons)
    148148    {
     
    165165
    166166
    167 std::set<gfan::ZVector> raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector u, const tropicalStrategy& currentStrategy)
     167gfan::ZMatrix raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector u, const tropicalStrategy& currentStrategy)
    168168{
    169169  std::set<gfan::ZCone> C = tropicalStar(I,r,u,currentStrategy);
    170   std::set<gfan::ZVector> raysOfC;
     170  gfan::ZMatrix raysOfC(0,u.size());
    171171  if (!currentStrategy.restrictToLowerHalfSpace())
    172172  {
     
    174174    {
    175175      assume(zc->dimensionOfLinealitySpace()+1 == zc->dimension());
    176       raysOfC.insert(zc->semiGroupGeneratorOfRay());
     176      raysOfC.appendRow(zc->semiGroupGeneratorOfRay());
    177177    }
    178178  }
     
    182182    {
    183183      assume(zc->dimensionOfLinealitySpace()+2 == zc->dimension());
    184       raysOfC.insert(zc->getRelativeInteriorPoint());
     184      raysOfC.appendRow(zc->getRelativeInteriorPoint());
    185185    }
    186186  }
  • Singular/dyn_modules/gfanlib/tropicalCurves.h

    r9c0326 rb71400a  
    99std::set<gfan::ZCone> tropicalStar(const ideal I, const ring r, const gfan::ZVector u,
    1010                                   const tropicalStrategy& currentStrategy);
    11 std::set<gfan::ZVector> raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector u,
    12                                            const tropicalStrategy& currentStrategy);
     11gfan::ZMatrix raysOfTropicalStar(ideal I, const ring r, const gfan::ZVector u,
     12                                 const tropicalStrategy& currentStrategy);
    1313
    1414#ifndef NDEBUG
  • Singular/dyn_modules/gfanlib/tropicalStrategy.cc

    r9c0326 rb71400a  
    195195  startingRing(rCopy(currentStrategy.getStartingRing())),
    196196  startingIdeal(id_Copy(currentStrategy.getStartingIdeal(),currentStrategy.getStartingRing())),
    197   uniformizingParameter(n_Copy(currentStrategy.getUniformizingParameter(),startingRing->cf)),
    198   shortcutRing(rCopy(currentStrategy.getShortcutRing())),
     197  uniformizingParameter(NULL),
     198  shortcutRing(NULL),
    199199  onlyLowerHalfSpace(currentStrategy.restrictToLowerHalfSpace()),
    200200  weightAdjustingAlgorithm1(currentStrategy.weightAdjustingAlgorithm1),
    201201  weightAdjustingAlgorithm2(currentStrategy.weightAdjustingAlgorithm2),
    202202  extraReductionAlgorithm(currentStrategy.extraReductionAlgorithm)
     203{
     204  if (originalRing) rTest(originalRing);
     205  if (originalIdeal) id_Test(originalIdeal,originalRing);
     206  if (startingRing) rTest(startingRing);
     207  if (startingIdeal) id_Test(startingIdeal,startingRing);
     208  if (currentStrategy.getUniformizingParameter())
     209  {
     210    uniformizingParameter = n_Copy(currentStrategy.getUniformizingParameter(),startingRing->cf);
     211    n_Test(uniformizingParameter,startingRing->cf);
     212  }
     213  if (currentStrategy.getShortcutRing())
     214  {
     215    shortcutRing = rCopy(currentStrategy.getShortcutRing());
     216    rTest(shortcutRing);
     217  }
     218}
     219
     220tropicalStrategy::~tropicalStrategy()
    203221{
    204222  if (originalRing) rTest(originalRing);
     
    208226  if (uniformizingParameter) n_Test(uniformizingParameter,startingRing->cf);
    209227  if (shortcutRing) rTest(shortcutRing);
    210 }
    211 
    212 tropicalStrategy::~tropicalStrategy()
    213 {
    214   if (originalRing) rTest(originalRing);
    215   if (originalIdeal) id_Test(originalIdeal,originalRing);
    216   if (startingRing) rTest(startingRing);
    217   if (startingIdeal) id_Test(startingIdeal,startingRing);
    218   if (uniformizingParameter) n_Test(uniformizingParameter,startingRing->cf);
    219   if (shortcutRing) rTest(shortcutRing);
    220228
    221229  id_Delete(&originalIdeal,originalRing);
    222230  rDelete(originalRing);
    223   id_Delete(&startingIdeal,startingRing);
    224   n_Delete(&uniformizingParameter,startingRing->cf);
    225   rDelete(startingRing);
    226   rDelete(shortcutRing);
     231  if (startingIdeal) id_Delete(&startingIdeal,startingRing);
     232  if (uniformizingParameter) n_Delete(&uniformizingParameter,startingRing->cf);
     233  if (startingRing) rDelete(startingRing);
     234  if (shortcutRing) rDelete(shortcutRing);
    227235}
    228236
     
    251259}
    252260
     261void tropicalStrategy::putUniformizingBinomialInFront(ideal I, const ring r, const number q) const
     262{
     263  poly p = p_One(r);
     264  p_SetCoeff(p,q,r);
     265  poly t = p_One(r);
     266  p_SetExp(t,1,1,r);
     267  p_Setm(t,r);
     268  poly pt = p_Add_q(p,p_Neg(t,r),r);
     269
     270  int k = idSize(I);
     271  int l;
     272  for (l=0; l<k; l++)
     273  {
     274    if (p_EqualPolys(I->m[l],pt,r))
     275      break;
     276  }
     277  p_Delete(&pt,r);
     278
     279  if (l>1)
     280  {
     281    pt = I->m[l];
     282    for (int i=l; i>0; i--)
     283      I->m[l] = I->m[l-1];
     284    I->m[0] = pt;
     285    pt = NULL;
     286  }
     287  return;
     288}
     289
    253290bool tropicalStrategy::reduce(ideal I, const ring r) const
    254291{
    255292  rTest(r);
    256293  id_Test(I,r);
    257   nMapFunc nMap = n_SetMap(startingRing->cf,r->cf);
    258   number p = nMap(uniformizingParameter,startingRing->cf,r->cf);
    259   bool b = this->extraReductionAlgorithm(I,r,p);
     294
     295  if (isValuationTrivial())
     296    return false;
     297
     298  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
     299  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
     300  bool b = extraReductionAlgorithm(I,r,p);
     301  // putUniformizingBinomialInFront(I,r,p);
    260302  n_Delete(&p,r->cf);
     303
     304  return b;
     305}
     306
     307bool tropicalStrategy::pReduce(ideal I, const ring r) const
     308{
     309  rTest(r);
     310  id_Test(I,r);
     311
     312  if (isValuationTrivial())
     313    return false;
     314
     315  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
     316  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
     317  bool b = pReduce0(I,p,r);
     318  n_Delete(&p,r->cf);
     319
    261320  return b;
    262321}
     
    293352
    294353  // if valuation non-trivial, change coefficient ring to residue field
    295   if (valuationIsNonTrivial())
     354  if (isValuationNonTrivial())
    296355  {
    297356    nKillChar(rShortcut->cf);
     
    330389  // if it exsists, return a copy of the monomial in the input ring
    331390  poly p = checkForMonomialViaSuddenSaturation(inIShortcut,rShortcut);
    332   poly monomial = NULL;
     391  poly monomial = NULL; int n = rVar(r);
    333392  if (p!=NULL)
    334393  {
     
    354413}
    355414
    356 ideal getWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
     415ideal tropicalStrategy::getWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
    357416{
    358417  // if the valuation is trivial and the ring and ideal have not been extended,
     
    365424  else
    366425  {
    367     checkFirstGenerator(I,r);
    368426    assume(idSize(inI)==idSize(I));
     427    int uni = findPositionOfUniformizingBinomial(I,r);
     428    assume(uni>=0);
    369429    /**
    370430     * change ground ring into finite field
     
    373433    ring rShortcut = copyAndChangeCoefficientRing(r);
    374434
    375     int k = idSize(inI);
     435    int k = idSize(inJ);
    376436    int l = idSize(I);
    377437    ideal inJShortcut = idInit(k);
     
    390450     */
    391451    matrix QShortcut = divisionDiscardingRemainder(inJShortcut,inIShortcut,rShortcut);
    392     matrix Q = mpNew(k,l);
     452    matrix Q = mpNew(l,k);
    393453    nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
    394454    for (int ij=k*l-1; ij>=0; ij--)
    395       Q->m[ij] = p_PermPoly(Q->m[ij],NULL,rShortcut,r,takingRepresentatives,NULL,0);
     455      Q->m[ij] = p_PermPoly(QShortcut->m[ij],NULL,rShortcut,r,takingRepresentatives,NULL,0);
     456
     457    nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
     458    number p = identity(uniformizingParameter,startingRing->cf,r->cf);
    396459
    397460    /**
     
    399462     */
    400463    ideal J = idInit(k);
    401     for (int j=0; j<l; j++)
     464    for (int j=0; j<k; j++)
    402465    {
    403466      poly q0 = p_Copy(inJ->m[j],r);
    404       for (int i=1; i<k; i++)
    405         q0 = p_Add_q(q0,p_Neg(p_Mult_q(p_Copy(MATELEM(Q,i,j),r),p_Copy(inI->m[i],r),r),r),r);
    406       J->m[j] = p_Mult_q(q0,p_Copy(J->m[0],r),r);
    407       q0 = NULL;
    408       for (int i=1; i<k; i++)
    409         J->m[j] = p_Add_q(I->m[j],p_Mult_q(p_Copy(MATELEM(Q,i,j),r),p_Copy(I->m[i],r),r),r);
     467      for (int i=0; i<l; i++)
     468      {
     469        poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
     470        poly inIi = p_Copy(inI->m[i],r);
     471        q0 = p_Add_q(q0,p_Neg(p_Mult_q(qij,inIi,r),r),r);
     472      }
     473      q0 = p_Div_nn(q0,p,r);
     474      poly q0g0 = p_Mult_q(q0,p_Copy(I->m[uni],r),r);
     475      // q0 = NULL;
     476      poly qigi = NULL;
     477      for (int i=0; i<l; i++)
     478      {
     479        poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
     480        // poly inIi = p_Copy(I->m[i],r);
     481        poly Ii = p_Copy(I->m[i],r);
     482        qigi = p_Add_q(qigi,p_Mult_q(qij,Ii,r),r);
     483      }
     484      J->m[j] = p_Add_q(q0g0,qigi,r);
    410485    }
    411486
     487    id_Delete(&inIShortcut,rShortcut);
     488    id_Delete(&inJShortcut,rShortcut);
     489    mp_Delete(&QShortcut,rShortcut);
     490    rDelete(rShortcut);
     491    mp_Delete(&Q,r);
     492    n_Delete(&p,r->cf);
    412493    return J;
    413494  }
     
    437518    inJ = idInit(k+1);
    438519    inJ->m[0] = p_One(r);
    439     p_SetCoeff(inJ->m[0],nMap(uniformizingParameter,startingRing->cf,r->cf),r);
     520    nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
     521    p_SetCoeff(inJ->m[0],identity(uniformizingParameter,startingRing->cf,r->cf),r);
    440522    for (int i=0; i<k; i++)
    441523      inJ->m[i+1] = p_PermPoly(inJShortcut->m[i],NULL,rShortcut,r,takingRepresentatives,NULL,0);
     524
     525    id_Delete(&inIShortcut,rShortcut);
     526    rDelete(rShortcut);
    442527  }
    443528  return inJ;
    444529}
    445530
    446 ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector w, const gfan::ZVector v)
     531static void deleteOrdering(ring r)
     532{
     533  if (r->order != NULL)
     534  {
     535    int i=rBlocks(r);
     536    assume(r->block0 != NULL && r->block1 != NULL && r->wvhdl != NULL);
     537    /* delete order */
     538    omFreeSize((ADDRESS)r->order,i*sizeof(int));
     539    omFreeSize((ADDRESS)r->block0,i*sizeof(int));
     540    omFreeSize((ADDRESS)r->block1,i*sizeof(int));
     541    /* delete weights */
     542    for (int j=0; j<i; j++)
     543      if (r->wvhdl[j]!=NULL)
     544        omFree(r->wvhdl[j]);
     545    omFreeSize((ADDRESS)r->wvhdl,i*sizeof(int *));
     546  }
     547  else
     548    assume(r->block0 == NULL && r->block1 == NULL && r->wvhdl == NULL);
     549  return;
     550}
     551
     552ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector w, const gfan::ZVector v) const
    447553{
    448554  // copy shortcutRing and change to desired ordering
     
    452558  deleteOrdering(s);
    453559  gfan::ZVector wAdjusted = adjustWeightForHomogeneity(w);
    454   gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(wAdjusted,v);
    455   s->order = (int*) omAlloc0(4*sizeof(int));
    456   s->block0 = (int*) omAlloc0(4*sizeof(int));
    457   s->block1 = (int*) omAlloc0(4*sizeof(int));
    458   s->wvhdl = (int**) omAlloc0(4*sizeof(int**));
    459   s->order[0] = ringorder_a;
    460   s->block0[0] = 1;
    461   s->block1[0] = n;
    462   s->wvhdl[0] = ZVectorToIntStar(adjustedInteriorPoint,ok);
    463   s->order[1] = ringorder_wp;
    464   s->block0[1] = 1;
    465   s->block1[1] = n;
    466   s->wvhdl[1] = ZVectorToIntStar(adjustedFacetNormal,ok);
    467   s->order[2] = ringorder_C;
    468   rComplete(s);
    469   rTest(s);
    470 
    471   return s;
    472 }
    473 
    474 ring tropicalStrategy::copyAndChangeOrderingDP(const ring r, const gfan::ZVector w, const gfan::ZVector v)
    475 {
    476   // copy shortcutRing and change to desired ordering
    477   bool ok;
    478   ring s = rCopy0(r);
    479   int n = rVar(s);
    480   deleteOrdering(s);
    481   gfan::ZVector wAdjusted = adjustWeightForHomogeneity(w);
    482   gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(wAdjusted,v);
     560  gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(v,wAdjusted);
    483561  s->order = (int*) omAlloc0(5*sizeof(int));
    484562  s->block0 = (int*) omAlloc0(5*sizeof(int));
     
    488566  s->block0[0] = 1;
    489567  s->block1[0] = n;
    490   s->wvhdl[0] = ZVectorToIntStar(adjustedInteriorPoint,ok);
     568  s->wvhdl[0] = ZVectorToIntStar(wAdjusted,ok);
    491569  s->order[1] = ringorder_a;
    492570  s->block0[1] = 1;
    493571  s->block1[1] = n;
    494   s->wvhdl[1] = ZVectorToIntStar(adjustedFacetNormal,ok);
     572  s->wvhdl[1] = ZVectorToIntStar(vAdjusted,ok);
    495573  s->order[2] = ringorder_dp;
    496574  s->block0[2] = 1;
     
    503581}
    504582
    505 
    506 std::pair<ideal,ring> tropicalStrategy::getFlip(const ideal I, const ring r,
     583ring tropicalStrategy::copyAndChangeOrderingLS(const ring r, const gfan::ZVector w, const gfan::ZVector v) const
     584{
     585  // copy shortcutRing and change to desired ordering
     586  bool ok;
     587  ring s = rCopy0(r);
     588  int n = rVar(s);
     589  deleteOrdering(s);
     590  s->order = (int*) omAlloc0(5*sizeof(int));
     591  s->block0 = (int*) omAlloc0(5*sizeof(int));
     592  s->block1 = (int*) omAlloc0(5*sizeof(int));
     593  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
     594  s->order[0] = ringorder_a;
     595  s->block0[0] = 1;
     596  s->block1[0] = n;
     597  s->wvhdl[0] = ZVectorToIntStar(w,ok);
     598  s->order[1] = ringorder_a;
     599  s->block0[1] = 1;
     600  s->block1[1] = n;
     601  s->wvhdl[1] = ZVectorToIntStar(v,ok);
     602  s->order[2] = ringorder_dp;
     603  s->block0[2] = 1;
     604  s->block1[2] = n;
     605  s->order[3] = ringorder_C;
     606  rComplete(s);
     607  rTest(s);
     608
     609  return s;
     610}
     611
     612
     613std::pair<ideal,ring> tropicalStrategy::getFlip(const ideal I, const ideal redI, const ring r,
    507614                                                const gfan::ZVector interiorPoint,
    508615                                                const gfan::ZVector facetNormal) const
    509616{
    510617  assume(isValuationTrivial() || interiorPoint[0].sign()<0);
    511   assume(checkFirstGenerators(I,r));
     618  assume(checkForUniformizingBinomial(I,r));
    512619
    513620  // get a generating system of the initial ideal
    514621  // and compute a standard basis with respect to adjacent ordering
    515   ideal inIr = initial(I,r,interiorPoint);
     622  ideal inIr = initial(redI,r,interiorPoint);
    516623  ring sAdjusted = copyAndChangeOrderingWP(r,interiorPoint,facetNormal);
    517624  nMapFunc identity = n_SetMap(r->cf,sAdjusted->cf);
    518   int k = idSize(I); ideal inIs = idInit(k);
     625  int k = idSize(I); ideal inIsAdjusted = idInit(k);
    519626  for (int i=0; i<k; i++)
    520627    inIsAdjusted->m[i] = p_PermPoly(inIr->m[i],NULL,r,sAdjusted,identity,NULL,0);
    521   ideal inJsAdjusted = getStdOfInitialIdeal(inIsAdjusted,sShortcut);
     628  ideal inJsAdjusted = getStdOfInitialIdeal(inIsAdjusted,sAdjusted);
    522629
    523630  // find witnesses of the new standard basis elements of the initial ideal
     
    529636    inJr->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,sAdjusted,r,identity,NULL,0);
    530637
    531   ideal Jr = getWitness(inJr,inIr,I,r);
    532   ring s = copyAndChangeOrderingDP(r,interiorPoint,facetNormal);
     638  ideal Jr = getWitness(inJr,inIr,redI,r);
     639  ring s = copyAndChangeOrderingLS(r,interiorPoint,facetNormal);
    533640  identity = n_SetMap(r->cf,s->cf);
     641  ideal Js = idInit(k);
    534642  for (int i=0; i<k; i++)
    535     Js->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,r,s,identity,NULL,0);
    536 
     643    Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identity,NULL,0);
     644
     645  // this->reduce(Jr,r);
     646  // cleanup
     647  id_Delete(&inJsAdjusted,sAdjusted);
     648  rDelete(sAdjusted);
     649  id_Delete(&Jr,r);
     650  id_Delete(&inJr,r);
     651
     652  assume(isValuationTrivial() || isOrderingLocalInT(s));
    537653  return std::make_pair(Js,s);
    538654}
     655
     656
     657bool tropicalStrategy::checkForUniformizingBinomial(const ideal I, const ring r) const
     658{
     659  // if the valuation is trivial,
     660  // then there is no special condition the first generator has to fullfill
     661  if (isValuationTrivial())
     662    return true;
     663
     664  // if the valuation is non-trivial then checks if the first generator is p-t
     665  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
     666  poly p = p_One(r);
     667  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
     668  poly t = p_One(r);
     669  p_SetExp(t,1,1,r);
     670  p_Setm(t,r);
     671  poly pt = p_Add_q(p,p_Neg(t,r),r);
     672
     673  for (int i=0; i<idSize(I); i++)
     674  {
     675    if (p_EqualPolys(I->m[i],pt,r))
     676    {
     677      p_Delete(&pt,r);
     678      return true;
     679    }
     680  }
     681  p_Delete(&pt,r);
     682  return false;
     683}
     684
     685int tropicalStrategy::findPositionOfUniformizingBinomial(const ideal I, const ring r) const
     686{
     687  assume(isValuationNonTrivial());
     688
     689  // if the valuation is non-trivial then checks if the first generator is p-t
     690  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
     691  poly p = p_One(r);
     692  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
     693  poly t = p_One(r);
     694  p_SetExp(t,1,1,r);
     695  p_Setm(t,r);
     696  poly pt = p_Add_q(p,p_Neg(t,r),r);
     697
     698  for (int i=0; i<idSize(I); i++)
     699  {
     700    if (p_EqualPolys(I->m[i],pt,r))
     701    {
     702      p_Delete(&pt,r);
     703      return i;
     704    }
     705  }
     706  p_Delete(&pt,r);
     707  return -1;
     708}
     709
     710bool tropicalStrategy::checkForUniformizingParameter(const ideal inI, const ring r) const
     711{
     712  // if the valuation is trivial,
     713  // then there is no special condition the first generator has to fullfill
     714  if (isValuationTrivial())
     715    return true;
     716
     717  // if the valuation is non-trivial then checks if the first generator is p
     718  if (inI->m[0]==NULL)
     719    return false;
     720  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
     721  poly p = p_One(r);
     722  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
     723
     724  for (int i=0; i<idSize(inI); i++)
     725  {
     726    if (p_EqualPolys(inI->m[i],p,r))
     727    {
     728      p_Delete(&p,r);
     729      return true;
     730    }
     731  }
     732  p_Delete(&p,r);
     733  return false;
     734}
  • Singular/dyn_modules/gfanlib/tropicalStrategy.h

    r9c0326 rb71400a  
    9191  bool (*extraReductionAlgorithm) (ideal I, ring r, number p);
    9292
     93  ring copyAndChangeCoefficientRing(const ring r) const;
     94  ring copyAndChangeOrderingWP(const ring r, const gfan::ZVector w, const gfan::ZVector v) const;
     95  ring copyAndChangeOrderingLS(const ring r, const gfan::ZVector w, const gfan::ZVector v) const;
     96
     97  /**
     98   * if valuation non-trivial, checks whether the generating system contains p-t
     99   * otherwise returns true
     100   */
     101  bool checkForUniformizingBinomial(const ideal I, const ring r) const;
     102
     103  /**
     104   * if valuation non-trivial, checks whether the genearting system contains p
     105   * otherwise returns true
     106   */
     107  bool checkForUniformizingParameter(const ideal inI, const ring r) const;
     108  int findPositionOfUniformizingBinomial(const ideal I, const ring r) const;
     109  void putUniformizingBinomialInFront(ideal I, const ring r, const number q) const;
     110
     111
    93112public:
    94113
     
    154173  ring getStartingRing() const
    155174  {
    156     rTest(startingRing);
     175    if (startingRing) rTest(startingRing);
    157176    return startingRing;
    158177  }
     
    186205  ring getShortcutRing() const
    187206  {
    188     rTest(shortcutRing);
     207    if (shortcutRing) rTest(shortcutRing);
    189208    return shortcutRing;
    190209  }
     
    258277  void changeCoefficientToResidueField(ring r) const;
    259278
    260 
    261 
    262279  /**
    263280   * reduces the generators of an ideal I so that
     
    265282   */
    266283  bool reduce(ideal I, const ring r) const;
     284
     285  bool pReduce(ideal I, const ring r) const;
    267286
    268287  /**
     
    292311   * computes the groebner cone adjacent to it
    293312   */
    294   std::pair<ideal,ring> getFlip(const ideal I, const ring r, const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
    295 
    296   /**
    297    * checks whether the first elements of I is p-t
    298    */
    299   static bool checkFirstGenerator(ideal I, ring r)
    300   {
    301     // check for correct size of ideal and first generator
    302     if (idSize(I)==0)
    303     {
    304       WerrorS("checkFirstElement: ideal empty");
    305       return false;
    306     }
    307     poly p = I->m[0];
    308     if (size(p)!=2)
    309     {
    310       WerrorS("checkFirstElement: no binomial");
    311       return false;
    312     }
    313 
    314     // check first term
    315     int n = rVar(r);
    316     if (!n_IsUnit(p_GetCoeff(p,r->cf),r->cf))
    317     {
    318       WerrorS("checkFirstElement: wrong coefficient");
    319       return false;
    320     }
    321     for (int i=1; i<=n; i++)
    322     {
    323       if (p_GetExp(p,i,r)!=0)
    324       {
    325         WerrorS("checkFirstElement: wrong exponent vector");
    326         return false;
    327       }
    328     }
    329 
    330     // check second term
    331     pIter(p);
    332     if (!n_IsOne(p_GetCoeff(p,r->cf),r->cf))
    333     {
    334       WerrorS("checkFirstElement: wrong coefficient");
    335       return false;
    336     }
    337     if (p_GetExp(p,1,r)!=1)
    338     {
    339       WerrorS("checkFirstElement: wrong exponent vector");
    340       return false;
    341     }
    342     for (int i=2; i<=n; i++)
    343     {
    344       if (p_GetExp(p,i,r)!=0)
    345       {
    346         WerrorS("checkFirstElement: wrong exponent vector");
    347         return false;
    348       }
    349     }
    350 
    351     return true;
    352   }
     313  std::pair<ideal,ring> getFlip(const ideal I, const ideal redI, const ring r, const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal) const;
    353314};
    354315
  • Singular/dyn_modules/gfanlib/tropicalTraversal.cc

    r9c0326 rb71400a  
    88  while(!workingList.empty())
    99  {
    10     std::cout << "starting traversal" << std::endl;
     10    // std::cout << "starting traversal" << std::endl;
    1111    const groebnerCone sigma=*(workingList.begin());
    1212    const groebnerCones neighbours = sigma.tropicalNeighbours();
     
    1818    tropicalVariety.insert(sigma);
    1919    workingList.erase(sigma);
    20     std::cout << "tropicalVariety.size():" << tropicalVariety.size() << std::endl;
    21     std::cout << "workingList.size():" << workingList.size() << std::endl;
     20    // std::cout << "tropicalVariety.size():" << tropicalVariety.size() << std::endl;
     21    // std::cout << "workingList.size():" << workingList.size() << std::endl;
    2222  }
    2323  return tropicalVariety;
  • Singular/dyn_modules/gfanlib/witness.cc

    r9c0326 rb71400a  
    3535  ring origin = currRing;
    3636  if (origin != r) rChangeCurrRing(r);
    37   ideal m = idLift(G,F);
    38   matrix Q = id_Module2formatedMatrix(m,IDELEMS(G),1,currRing);
     37  ideal R; matrix U;
     38  ideal m = idLift(G,F,&R,FALSE,FALSE,TRUE,&U);
     39  matrix Q = id_Module2formatedMatrix(m,IDELEMS(G),IDELEMS(F),currRing);
     40  id_Delete(&R,r);
     41  mp_Delete(&U,r);
    3942  if (origin != r) rChangeCurrRing(origin);
    4043  return Q;
  • gfanlib/gfanlib_z.h

    r9c0326 rb71400a  
    1010
    1111#include <string.h>
     12#include <iostream>
    1213#include <ostream>
    1314
     
    6869    freefunc(str,strlen(str)+1);
    6970    return f;
     71  }
     72  void debugPrint() const
     73  {
     74    void (*freefunc)(void *, size_t);
     75    mp_get_memory_functions(0,0,&freefunc);
     76    char *str=mpz_get_str(0,10,value);
     77    std::cout << str;
     78    freefunc(str,strlen(str)+1);
     79    return;
    7080  }
    7181  Integer& operator+=(const Integer& a)
     
    330340      return f;
    331341    }
     342    friend void debugPrint(IntegerTemplate const &a)
     343    {
     344      std::cout << a << std::endl;
     345      return;
     346    }
    332347    LimbWord signExtension(LimbWord a)
    333348    {
  • libpolys/polys/matpol.cc

    r9c0326 rb71400a  
    88
    99#include <stdio.h>
     10#include <string.h>
     11#include <sstream>
    1012#include <math.h>
    1113
     
    756758  }
    757759}
     760
     761void mp_Write(matrix im, const ring r)
     762{
     763  int i,ii = MATROWS(im)-1;
     764  int j,jj = MATCOLS(im)-1;
     765  poly *pp = im->m;
     766  for (i=0; i<=ii; i++)
     767  {
     768    for (j=0; j<=jj; j++)
     769    {
     770      Print("matrix[%u,%u]=",i+1,j+1);
     771      if ((i<ii)||(j<jj)) p_Write(*pp++, r);
     772      else                p_Write0(*pp, r);
     773    }
     774  }
     775  return;
     776}
     777
    758778
    759779char * iiStringMatrix(matrix im, int dim, const ring r, char ch)
  • libpolys/polys/matpol.h

    r9c0326 rb71400a  
    1010// #include <kernel/structs.h>
    1111#include <polys/monomials/ring.h>
     12#include <string.h>
    1213
    1314// THIS IS REALLY DIRTY: ip_smatrix HAS TO BE IDENTICAL TO ip_sideal
     
    9192char * iiStringMatrix(matrix im, int dim, const ring r, char ch=',');
    9293
     94void mp_Write(matrix im, const ring r);
     95
    9396extern omBin ip_smatrix_bin;
    9497
Note: See TracChangeset for help on using the changeset viewer.