My Project
tropicalStrategy.cc
Go to the documentation of this file.
1 #include "tropicalStrategy.h"
2 #include "singularWishlist.h"
3 #include "adjustWeights.h"
4 #include "ppinitialReduction.h"
5 // #include "ttinitialReduction.h"
6 #include "tropical.h"
7 #include "std_wrapper.h"
8 #include "tropicalCurves.h"
9 #include "tropicalDebug.h"
10 #include "containsMonomial.h"
11 
12 
13 // for various commands in dim(ideal I, ring r):
14 #include "kernel/ideals.h"
16 #include "kernel/GBEngine/kstd1.h"
17 #include "misc/prime.h" // for isPrime(int i)
18 
19 /***
20  * Computes the dimension of an ideal I in ring r
21  * Copied from jjDim in iparith.cc
22  **/
23 int dim(ideal I, ring r)
24 {
25  ring origin = currRing;
26  if (origin != r)
27  rChangeCurrRing(r);
28  int d;
30  {
31  int i = idPosConstant(I);
32  if ((i != -1)
33  #ifdef HAVE_RINGS
34  && (n_IsUnit(p_GetCoeff(I->m[i],currRing->cf),currRing->cf))
35  #endif
36  )
37  return -1;
38  ideal vv = id_Head(I,currRing);
39  if (i != -1) pDelete(&vv->m[i]);
40  d = scDimInt(vv, currRing->qideal);
41  if (rField_is_Z(currRing) && (i==-1)) d++;
42  idDelete(&vv);
43  return d;
44  }
45  else
46  d = scDimInt(I,currRing->qideal);
47  if (origin != r)
48  rChangeCurrRing(origin);
49  return d;
50 }
51 
52 static void swapElements(ideal I, ideal J)
53 {
54  assume(IDELEMS(I)==IDELEMS(J));
55 
56  for (int i=IDELEMS(I)-1; i>=0; i--)
57  {
58  poly cache = I->m[i];
59  I->m[i] = J->m[i];
60  J->m[i] = cache;
61  }
62 }
63 
64 static bool noExtraReduction(ideal I, ring r, number /*p*/)
65 {
66  int n = rVar(r);
67  gfan::ZVector allOnes(n);
68  for (int i=0; i<n; i++)
69  allOnes[i] = 1;
70  ring rShortcut = rCopy0(r);
71 
72  rRingOrder_t* order = rShortcut->order;
73  int* block0 = rShortcut->block0;
74  int* block1 = rShortcut->block1;
75  int** wvhdl = rShortcut->wvhdl;
76 
77  int h = rBlocks(r);
78  rShortcut->order = (rRingOrder_t*) omAlloc0((h+2)*sizeof(rRingOrder_t));
79  rShortcut->block0 = (int*) omAlloc0((h+2)*sizeof(int));
80  rShortcut->block1 = (int*) omAlloc0((h+2)*sizeof(int));
81  rShortcut->wvhdl = (int**) omAlloc0((h+2)*sizeof(int*));
82  rShortcut->order[0] = ringorder_a;
83  rShortcut->block0[0] = 1;
84  rShortcut->block1[0] = n;
85  bool overflow;
86  rShortcut->wvhdl[0] = ZVectorToIntStar(allOnes,overflow);
87  for (int i=1; i<=h; i++)
88  {
89  rShortcut->order[i] = order[i-1];
90  rShortcut->block0[i] = block0[i-1];
91  rShortcut->block1[i] = block1[i-1];
92  rShortcut->wvhdl[i] = wvhdl[i-1];
93  }
94  //rShortcut->order[h+1] = (rRingOrder_t)0; -- done by omAlloc0
95  //rShortcut->block0[h+1] = 0;
96  //rShortcut->block1[h+1] = 0;
97  //rShortcut->wvhdl[h+1] = NULL;
98 
99  rComplete(rShortcut);
100  rTest(rShortcut);
101 
102  omFree(order);
103  omFree(block0);
104  omFree(block1);
105  omFree(wvhdl);
106 
107  int k = IDELEMS(I);
108  ideal IShortcut = idInit(k);
109  nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
110  for (int i=0; i<k; i++)
111  {
112  if(I->m[i]!=NULL)
113  {
114  IShortcut->m[i] = p_PermPoly(I->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
115  }
116  }
117 
118  ideal JShortcut = gfanlib_kStd_wrapper(IShortcut,rShortcut);
119 
120  ideal J = idInit(k);
121  nMapFunc outofShortcut = n_SetMap(rShortcut->cf,r->cf);
122  for (int i=0; i<k; i++)
123  J->m[i] = p_PermPoly(JShortcut->m[i],NULL,rShortcut,r,outofShortcut,NULL,0);
124 
125  assume(areIdealsEqual(J,r,I,r));
126  swapElements(I,J);
127  id_Delete(&IShortcut,rShortcut);
128  id_Delete(&JShortcut,rShortcut);
129  rDelete(rShortcut);
130  id_Delete(&J,r);
131  return false;
132 }
133 
134 /**
135  * Initializes all relevant structures and information for the trivial valuation case,
136  * i.e. computing a tropical variety without any valuation.
137  */
138 tropicalStrategy::tropicalStrategy(const ideal I, const ring r,
139  const bool completelyHomogeneous,
140  const bool completeSpace):
141  originalRing(rCopy(r)),
142  originalIdeal(id_Copy(I,r)),
143  expectedDimension(dim(originalIdeal,originalRing)),
144  linealitySpace(homogeneitySpace(originalIdeal,originalRing)),
145  startingRing(rCopy(originalRing)),
146  startingIdeal(id_Copy(originalIdeal,originalRing)),
147  uniformizingParameter(NULL),
148  shortcutRing(NULL),
149  onlyLowerHalfSpace(false),
150  weightAdjustingAlgorithm1(nonvalued_adjustWeightForHomogeneity),
151  weightAdjustingAlgorithm2(nonvalued_adjustWeightUnderHomogeneity),
152  extraReductionAlgorithm(noExtraReduction)
153 {
154  assume(rField_is_Q(r) || rField_is_Zp(r) || rField_is_Z(r));
155  if (!completelyHomogeneous)
156  {
159  }
160  if (!completeSpace)
161  onlyLowerHalfSpace = true;
162 }
163 
164 /**
165  * Given a polynomial ring r over the rational numbers and a weighted ordering,
166  * returns a polynomial ring s over the integers with one extra variable, which is weighted -1.
167  */
168 static ring constructStartingRing(ring r)
169 {
170  assume(rField_is_Q(r));
171 
172  ring s = rCopy0(r,FALSE,FALSE);
173  nKillChar(s->cf);
174  s->cf = nInitChar(n_Z,NULL);
175 
176  int n = rVar(s)+1;
177  s->N = n;
178  char** oldNames = s->names;
179  s->names = (char**) omAlloc((n+1)*sizeof(char**));
180  s->names[0] = omStrDup("t");
181  for (int i=1; i<n; i++)
182  s->names[i] = oldNames[i-1];
183  omFree(oldNames);
184 
185  s->order = (rRingOrder_t*) omAlloc0(3*sizeof(rRingOrder_t));
186  s->block0 = (int*) omAlloc0(3*sizeof(int));
187  s->block1 = (int*) omAlloc0(3*sizeof(int));
188  s->wvhdl = (int**) omAlloc0(3*sizeof(int**));
189  s->order[0] = ringorder_ws;
190  s->block0[0] = 1;
191  s->block1[0] = n;
192  s->wvhdl[0] = (int*) omAlloc(n*sizeof(int));
193  s->wvhdl[0][0] = 1;
194  if (r->order[0] == ringorder_lp)
195  {
196  s->wvhdl[0][1] = 1;
197  }
198  else if (r->order[0] == ringorder_ls)
199  {
200  s->wvhdl[0][1] = -1;
201  }
202  else if (r->order[0] == ringorder_dp)
203  {
204  for (int i=1; i<n; i++)
205  s->wvhdl[0][i] = -1;
206  }
207  else if (r->order[0] == ringorder_ds)
208  {
209  for (int i=1; i<n; i++)
210  s->wvhdl[0][i] = 1;
211  }
212  else if (r->order[0] == ringorder_ws)
213  {
214  for (int i=1; i<n; i++)
215  s->wvhdl[0][i] = r->wvhdl[0][i-1];
216  }
217  else
218  {
219  for (int i=1; i<n; i++)
220  s->wvhdl[0][i] = -r->wvhdl[0][i-1];
221  }
222  s->order[1] = ringorder_C;
223 
224  rComplete(s);
225  rTest(s);
226  return s;
227 }
228 
229 static ideal constructStartingIdeal(ideal originalIdeal, ring originalRing, number uniformizingParameter, ring startingRing)
230 {
231  // construct p-t
232  poly g = p_One(startingRing);
233  p_SetCoeff(g,uniformizingParameter,startingRing);
234  pNext(g) = p_One(startingRing);
235  p_SetExp(pNext(g),1,1,startingRing);
236  p_SetCoeff(pNext(g),n_Init(-1,startingRing->cf),startingRing);
237  p_Setm(pNext(g),startingRing);
238  ideal pt = idInit(1);
239  pt->m[0] = g;
240 
241  // map originalIdeal from originalRing into startingRing
242  int k = IDELEMS(originalIdeal);
243  ideal J = idInit(k+1);
244  nMapFunc nMap = n_SetMap(originalRing->cf,startingRing->cf);
245  int n = rVar(originalRing);
246  int* shiftByOne = (int*) omAlloc((n+1)*sizeof(int));
247  for (int i=1; i<=n; i++)
248  shiftByOne[i]=i+1;
249  for (int i=0; i<k; i++)
250  {
251  if(originalIdeal->m[i]!=NULL)
252  {
253  J->m[i] = p_PermPoly(originalIdeal->m[i],shiftByOne,originalRing,startingRing,nMap,NULL,0);
254  }
255  }
256  omFreeSize(shiftByOne,(n+1)*sizeof(int));
257 
258  ring origin = currRing;
259  rChangeCurrRing(startingRing);
260  ideal startingIdeal = kNF(pt,startingRing->qideal,J); // mathematically redundant,
261  rChangeCurrRing(origin); // but helps with upcoming std computation
262  // ideal startingIdeal = J; J = NULL;
263  assume(startingIdeal->m[k]==NULL);
264  startingIdeal->m[k] = pt->m[0];
265  startingIdeal = gfanlib_kStd_wrapper(startingIdeal,startingRing);
266 
267  id_Delete(&J,startingRing);
268  pt->m[0] = NULL;
269  id_Delete(&pt,startingRing);
270  return startingIdeal;
271 }
272 
273 /***
274  * Initializes all relevant structures and information for the valued case,
275  * i.e. computing a tropical variety over the rational numbers with p-adic valuation
276  **/
277 tropicalStrategy::tropicalStrategy(ideal J, number q, ring s):
278  originalRing(rCopy(s)),
279  originalIdeal(id_Copy(J,s)),
280  expectedDimension(dim(originalIdeal,originalRing)+1),
281  linealitySpace(gfan::ZCone()), // to come, see below
282  startingRing(NULL), // to come, see below
283  startingIdeal(NULL), // to come, see below
284  uniformizingParameter(NULL), // to come, see below
285  shortcutRing(NULL), // to come, see below
286  onlyLowerHalfSpace(true),
287  weightAdjustingAlgorithm1(valued_adjustWeightForHomogeneity),
288  weightAdjustingAlgorithm2(valued_adjustWeightUnderHomogeneity),
289  extraReductionAlgorithm(ppreduceInitially)
290 {
291  /* assume that the ground field of the originalRing is Q */
292  assume(rField_is_Q(s));
293 
294  /* replace Q with Z for the startingRing
295  * and add an extra variable for tracking the uniformizing parameter */
297 
298  /* map the uniformizing parameter into the new coefficient domain */
299  nMapFunc nMap = n_SetMap(originalRing->cf,startingRing->cf);
301 
302  /* map the input ideal into the new polynomial ring */
305 
307 
308  /* construct the shorcut ring */
309  shortcutRing = rCopy0(startingRing,FALSE); // do not copy q-ideal
310  nKillChar(shortcutRing->cf);
314 }
315 
317  originalRing(rCopy(currentStrategy.getOriginalRing())),
318  originalIdeal(id_Copy(currentStrategy.getOriginalIdeal(),currentStrategy.getOriginalRing())),
319  expectedDimension(currentStrategy.getExpectedDimension()),
320  linealitySpace(currentStrategy.getHomogeneitySpace()),
321  startingRing(rCopy(currentStrategy.getStartingRing())),
322  startingIdeal(id_Copy(currentStrategy.getStartingIdeal(),currentStrategy.getStartingRing())),
323  uniformizingParameter(NULL),
324  shortcutRing(NULL),
325  onlyLowerHalfSpace(currentStrategy.restrictToLowerHalfSpace()),
326  weightAdjustingAlgorithm1(currentStrategy.weightAdjustingAlgorithm1),
327  weightAdjustingAlgorithm2(currentStrategy.weightAdjustingAlgorithm2),
328  extraReductionAlgorithm(currentStrategy.extraReductionAlgorithm)
329 {
334  if (currentStrategy.getUniformizingParameter())
335  {
338  }
339  if (currentStrategy.getShortcutRing())
340  {
341  shortcutRing = rCopy(currentStrategy.getShortcutRing());
343  }
344 }
345 
347 {
354 
361 }
362 
364 {
365  originalRing = rCopy(currentStrategy.getOriginalRing());
366  originalIdeal = id_Copy(currentStrategy.getOriginalIdeal(),currentStrategy.getOriginalRing());
367  expectedDimension = currentStrategy.getExpectedDimension();
368  startingRing = rCopy(currentStrategy.getStartingRing());
369  startingIdeal = id_Copy(currentStrategy.getStartingIdeal(),currentStrategy.getStartingRing());
371  shortcutRing = rCopy(currentStrategy.getShortcutRing());
372  onlyLowerHalfSpace = currentStrategy.restrictToLowerHalfSpace();
376 
383 
384  return *this;
385 }
386 
387 void tropicalStrategy::putUniformizingBinomialInFront(ideal I, const ring r, const number q) const
388 {
389  poly p = p_One(r);
390  p_SetCoeff(p,q,r);
391  poly t = p_One(r);
392  p_SetExp(t,1,1,r);
393  p_Setm(t,r);
394  poly pt = p_Add_q(p,p_Neg(t,r),r);
395 
396  int k = IDELEMS(I);
397  int l;
398  for (l=0; l<k; l++)
399  {
400  if (p_EqualPolys(I->m[l],pt,r))
401  break;
402  }
403  p_Delete(&pt,r);
404 
405  if (l>1)
406  {
407  pt = I->m[l];
408  for (int i=l; i>0; i--)
409  I->m[l] = I->m[l-1];
410  I->m[0] = pt;
411  pt = NULL;
412  }
413  return;
414 }
415 
416 bool tropicalStrategy::reduce(ideal I, const ring r) const
417 {
418  rTest(r);
419  id_Test(I,r);
420 
421  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
422  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
423  bool b = extraReductionAlgorithm(I,r,p);
424  // putUniformizingBinomialInFront(I,r,p);
425  n_Delete(&p,r->cf);
426 
427  return b;
428 }
429 
430 void tropicalStrategy::pReduce(ideal I, const ring r) const
431 {
432  rTest(r);
433  id_Test(I,r);
434 
435  if (isValuationTrivial())
436  return;
437 
438  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
439  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
440  ::pReduce(I,p,r);
441  n_Delete(&p,r->cf);
442 
443  return;
444 }
445 
446 ring tropicalStrategy::getShortcutRingPrependingWeight(const ring r, const gfan::ZVector &v) const
447 {
448  ring rShortcut = rCopy0(r,FALSE); // do not copy q-ideal
449 
450  // save old ordering
451  rRingOrder_t* order = rShortcut->order;
452  int* block0 = rShortcut->block0;
453  int* block1 = rShortcut->block1;
454  int** wvhdl = rShortcut->wvhdl;
455 
456  // adjust weight and create new ordering
457  gfan::ZVector w = adjustWeightForHomogeneity(v);
458  int h = rBlocks(r); int n = rVar(r);
459  rShortcut->order = (rRingOrder_t*) omAlloc0((h+2)*sizeof(rRingOrder_t));
460  rShortcut->block0 = (int*) omAlloc0((h+2)*sizeof(int));
461  rShortcut->block1 = (int*) omAlloc0((h+2)*sizeof(int));
462  rShortcut->wvhdl = (int**) omAlloc0((h+2)*sizeof(int*));
463  rShortcut->order[0] = ringorder_a;
464  rShortcut->block0[0] = 1;
465  rShortcut->block1[0] = n;
466  bool overflow;
467  rShortcut->wvhdl[0] = ZVectorToIntStar(w,overflow);
468  for (int i=1; i<=h; i++)
469  {
470  rShortcut->order[i] = order[i-1];
471  rShortcut->block0[i] = block0[i-1];
472  rShortcut->block1[i] = block1[i-1];
473  rShortcut->wvhdl[i] = wvhdl[i-1];
474  }
475 
476  // if valuation non-trivial, change coefficient ring to residue field
477  if (isValuationNonTrivial())
478  {
479  nKillChar(rShortcut->cf);
480  rShortcut->cf = nCopyCoeff(shortcutRing->cf);
481  }
482  rComplete(rShortcut);
483  rTest(rShortcut);
484 
485  // delete old ordering
486  omFree(order);
487  omFree(block0);
488  omFree(block1);
489  omFree(wvhdl);
490 
491  return rShortcut;
492 }
493 
494 std::pair<poly,int> tropicalStrategy::checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector &w) const
495 {
496  // quick check whether I already contains an ideal
497  int k = IDELEMS(I);
498  for (int i=0; i<k; i++)
499  {
500  poly g = I->m[i];
501  if (g!=NULL
502  && pNext(g)==NULL
503  && (isValuationTrivial() || n_IsOne(p_GetCoeff(g,r),r->cf)))
504  return std::pair<poly,int>(g,i);
505  }
506 
507  ring rShortcut;
508  ideal inIShortcut;
509  if (w.size()>0)
510  {
511  // if needed, prepend extra weight for homogeneity
512  // switch to residue field if valuation is non trivial
513  rShortcut = getShortcutRingPrependingWeight(r,w);
514 
515  // compute the initial ideal and map it into the constructed ring
516  // if switched to residue field, remove possibly 0 elements
517  ideal inI = initial(I,r,w);
518  inIShortcut = idInit(k);
519  nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
520  for (int i=0; i<k; i++)
521  inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
522  if (isValuationNonTrivial())
523  idSkipZeroes(inIShortcut);
524  id_Delete(&inI,r);
525  }
526  else
527  {
528  rShortcut = r;
529  inIShortcut = I;
530  }
531 
532  gfan::ZCone C0 = homogeneitySpace(inIShortcut,rShortcut);
533  gfan::ZCone pos = gfan::ZCone::positiveOrthant(C0.ambientDimension());
534  gfan::ZCone C0pos = intersection(C0,pos);
535  C0pos.canonicalize();
536  gfan::ZVector wpos = C0pos.getRelativeInteriorPoint();
538 
539  // check initial ideal for monomial and
540  // if it exsists, return a copy of the monomial in the input ring
541  poly p = searchForMonomialViaStepwiseSaturation(inIShortcut,rShortcut,wpos);
542  poly monomial = NULL;
543  if (p!=NULL)
544  {
545  monomial=p_One(r);
546  for (int i=1; i<=rVar(r); i++)
547  p_SetExp(monomial,i,p_GetExp(p,i,rShortcut),r);
548  p_Setm(monomial,r);
549  p_Delete(&p,rShortcut);
550  }
551 
552 
553  if (w.size()>0)
554  {
555  // if needed, cleanup
556  id_Delete(&inIShortcut,rShortcut);
557  rDelete(rShortcut);
558  }
559  return std::pair<poly,int>(monomial,-1);
560 }
561 
563 {
564  ring rShortcut = rCopy0(r,FALSE); // do not copy q-ideal
565  nKillChar(rShortcut->cf);
566  rShortcut->cf = nCopyCoeff(shortcutRing->cf);
567  rComplete(rShortcut);
568  rTest(rShortcut);
569  return rShortcut;
570 }
571 
572 ideal tropicalStrategy::computeWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
573 {
574  // if the valuation is trivial and the ring and ideal have not been extended,
575  // then it is sufficient to return the difference between the elements of inJ
576  // and their normal forms with respect to I and r
577  if (isValuationTrivial())
578  return witness(inJ,I,r);
579  // if the valuation is non-trivial and the ring and ideal have been extended,
580  // then we can make a shortcut through the residue field
581  else
582  {
583  assume(IDELEMS(inI)==IDELEMS(I));
584  int uni = findPositionOfUniformizingBinomial(I,r);
585  assume(uni>=0);
586  /**
587  * change ground ring into finite field
588  * and map the data into it
589  */
590  ring rShortcut = copyAndChangeCoefficientRing(r);
591 
592  int k = IDELEMS(inJ);
593  int l = IDELEMS(I);
594  ideal inJShortcut = idInit(k);
595  ideal inIShortcut = idInit(l);
596  nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
597  for (int i=0; i<k; i++)
598  inJShortcut->m[i] = p_PermPoly(inJ->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
599  for (int j=0; j<l; j++)
600  inIShortcut->m[j] = p_PermPoly(inI->m[j],NULL,r,rShortcut,takingResidues,NULL,0);
601  id_Test(inJShortcut,rShortcut);
602  id_Test(inIShortcut,rShortcut);
603 
604  /**
605  * Compute a division with remainder over the finite field
606  * and map the result back to r
607  */
608  matrix QShortcut = divisionDiscardingRemainder(inJShortcut,inIShortcut,rShortcut);
609  matrix Q = mpNew(l,k);
610  nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
611  for (int ij=k*l-1; ij>=0; ij--)
612  Q->m[ij] = p_PermPoly(QShortcut->m[ij],NULL,rShortcut,r,takingRepresentatives,NULL,0);
613 
614  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
615  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
616 
617  /**
618  * Compute the normal forms
619  */
620  ideal J = idInit(k);
621  for (int j=0; j<k; j++)
622  {
623  poly q0 = p_Copy(inJ->m[j],r);
624  for (int i=0; i<l; i++)
625  {
626  poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
627  poly inIi = p_Copy(inI->m[i],r);
628  q0 = p_Add_q(q0,p_Neg(p_Mult_q(qij,inIi,r),r),r);
629  }
630  q0 = p_Div_nn(q0,p,r);
631  poly q0g0 = p_Mult_q(q0,p_Copy(I->m[uni],r),r);
632  // q0 = NULL;
633  poly qigi = NULL;
634  for (int i=0; i<l; i++)
635  {
636  poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
637  // poly inIi = p_Copy(I->m[i],r);
638  poly Ii = p_Copy(I->m[i],r);
639  qigi = p_Add_q(qigi,p_Mult_q(qij,Ii,r),r);
640  }
641  J->m[j] = p_Add_q(q0g0,qigi,r);
642  }
643 
644  id_Delete(&inIShortcut,rShortcut);
645  id_Delete(&inJShortcut,rShortcut);
646  mp_Delete(&QShortcut,rShortcut);
647  rDelete(rShortcut);
648  mp_Delete(&Q,r);
649  n_Delete(&p,r->cf);
650  return J;
651  }
652 }
653 
654 ideal tropicalStrategy::computeStdOfInitialIdeal(const ideal inI, const ring r) const
655 {
656  // if valuation trivial, then compute std as usual
657  if (isValuationTrivial())
658  return gfanlib_kStd_wrapper(inI,r);
659 
660  // if valuation non-trivial, then uniformizing parameter is in ideal
661  // so switch to residue field first and compute standard basis over the residue field
662  ring rShortcut = copyAndChangeCoefficientRing(r);
663  nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
664  int k = IDELEMS(inI);
665  ideal inIShortcut = idInit(k);
666  for (int i=0; i<k; i++)
667  inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
668  ideal inJShortcut = gfanlib_kStd_wrapper(inIShortcut,rShortcut);
669 
670  // and lift the result back to the ring with valuation
671  nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
672  k = IDELEMS(inJShortcut);
673  ideal inJ = idInit(k+1);
674  inJ->m[0] = p_One(r);
675  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
676  p_SetCoeff(inJ->m[0],identity(uniformizingParameter,startingRing->cf,r->cf),r);
677  for (int i=0; i<k; i++)
678  inJ->m[i+1] = p_PermPoly(inJShortcut->m[i],NULL,rShortcut,r,takingRepresentatives,NULL,0);
679 
680  id_Delete(&inJShortcut,rShortcut);
681  id_Delete(&inIShortcut,rShortcut);
682  rDelete(rShortcut);
683  return inJ;
684 }
685 
686 ideal tropicalStrategy::computeLift(const ideal inJs, const ring s, const ideal inIr, const ideal Ir, const ring r) const
687 {
688  int k = IDELEMS(inJs);
689  ideal inJr = idInit(k);
690  nMapFunc identitysr = n_SetMap(s->cf,r->cf);
691  for (int i=0; i<k; i++)
692  inJr->m[i] = p_PermPoly(inJs->m[i],NULL,s,r,identitysr,NULL,0);
693 
694  ideal Jr = computeWitness(inJr,inIr,Ir,r);
695  nMapFunc identityrs = n_SetMap(r->cf,s->cf);
696  ideal Js = idInit(k);
697  for (int i=0; i<k; i++)
698  Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identityrs,NULL,0);
699  return Js;
700 }
701 
702 ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
703 {
704  // copy shortcutRing and change to desired ordering
705  bool ok;
706  ring s = rCopy0(r,FALSE,FALSE);
707  int n = rVar(s);
708  gfan::ZVector wAdjusted = adjustWeightForHomogeneity(w);
709  gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(v,wAdjusted);
710  s->order = (rRingOrder_t*) omAlloc0(5*sizeof(rRingOrder_t));
711  s->block0 = (int*) omAlloc0(5*sizeof(int));
712  s->block1 = (int*) omAlloc0(5*sizeof(int));
713  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
714  s->order[0] = ringorder_a;
715  s->block0[0] = 1;
716  s->block1[0] = n;
717  s->wvhdl[0] = ZVectorToIntStar(wAdjusted,ok);
718  s->order[1] = ringorder_a;
719  s->block0[1] = 1;
720  s->block1[1] = n;
721  s->wvhdl[1] = ZVectorToIntStar(vAdjusted,ok);
722  s->order[2] = ringorder_lp;
723  s->block0[2] = 1;
724  s->block1[2] = n;
725  s->order[3] = ringorder_C;
726  rComplete(s);
727  rTest(s);
728 
729  return s;
730 }
731 
732 ring tropicalStrategy::copyAndChangeOrderingLS(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
733 {
734  // copy shortcutRing and change to desired ordering
735  bool ok;
736  ring s = rCopy0(r,FALSE,FALSE);
737  int n = rVar(s);
738  s->order = (rRingOrder_t*) omAlloc0(5*sizeof(rRingOrder_t));
739  s->block0 = (int*) omAlloc0(5*sizeof(int));
740  s->block1 = (int*) omAlloc0(5*sizeof(int));
741  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
742  s->order[0] = ringorder_a;
743  s->block0[0] = 1;
744  s->block1[0] = n;
745  s->wvhdl[0] = ZVectorToIntStar(w,ok);
746  s->order[1] = ringorder_a;
747  s->block0[1] = 1;
748  s->block1[1] = n;
749  s->wvhdl[1] = ZVectorToIntStar(v,ok);
750  s->order[2] = ringorder_lp;
751  s->block0[2] = 1;
752  s->block1[2] = n;
753  s->order[3] = ringorder_C;
754  rComplete(s);
755  rTest(s);
756 
757  return s;
758 }
759 
760 std::pair<ideal,ring> tropicalStrategy::computeFlip(const ideal Ir, const ring r,
761  const gfan::ZVector &interiorPoint,
762  const gfan::ZVector &facetNormal) const
763 {
764  assume(isValuationTrivial() || interiorPoint[0].sign()<0);
766  assume(checkWeightVector(Ir,r,interiorPoint));
767 
768  // get a generating system of the initial ideal
769  // and compute a standard basis with respect to adjacent ordering
770  ideal inIr = initial(Ir,r,interiorPoint);
771  ring sAdjusted = copyAndChangeOrderingWP(r,interiorPoint,facetNormal);
772  nMapFunc identity = n_SetMap(r->cf,sAdjusted->cf);
773  int k = IDELEMS(Ir);
774  ideal inIsAdjusted = idInit(k);
775  for (int i=0; i<k; i++)
776  inIsAdjusted->m[i] = p_PermPoly(inIr->m[i],NULL,r,sAdjusted,identity,NULL,0);
777  ideal inJsAdjusted = computeStdOfInitialIdeal(inIsAdjusted,sAdjusted);
778 
779  // find witnesses of the new standard basis elements of the initial ideal
780  // with the help of the old standard basis of the ideal
781  k = IDELEMS(inJsAdjusted);
782  ideal inJr = idInit(k);
783  identity = n_SetMap(sAdjusted->cf,r->cf);
784  for (int i=0; i<k; i++)
785  inJr->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,sAdjusted,r,identity,NULL,0);
786 
787  ideal Jr = computeWitness(inJr,inIr,Ir,r);
788  ring s = copyAndChangeOrderingLS(r,interiorPoint,facetNormal);
789  identity = n_SetMap(r->cf,s->cf);
790  ideal Js = idInit(k);
791  for (int i=0; i<k; i++)
792  Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identity,NULL,0);
793 
794  reduce(Js,s);
795  assume(areIdealsEqual(Js,s,Ir,r));
797  assume(checkWeightVector(Js,s,interiorPoint));
798 
799  // cleanup
800  id_Delete(&inIsAdjusted,sAdjusted);
801  id_Delete(&inJsAdjusted,sAdjusted);
802  rDelete(sAdjusted);
803  id_Delete(&inIr,r);
804  id_Delete(&Jr,r);
805  id_Delete(&inJr,r);
806 
808  return std::make_pair(Js,s);
809 }
810 
811 
812 bool tropicalStrategy::checkForUniformizingBinomial(const ideal I, const ring r) const
813 {
814  // if the valuation is trivial,
815  // then there is no special condition the first generator has to fullfill
816  if (isValuationTrivial())
817  return true;
818 
819  // if the valuation is non-trivial then checks if the first generator is p-t
820  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
821  poly p = p_One(r);
822  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
823  poly t = p_One(r);
824  p_SetExp(t,1,1,r);
825  p_Setm(t,r);
826  poly pt = p_Add_q(p,p_Neg(t,r),r);
827 
828  for (int i=0; i<IDELEMS(I); i++)
829  {
830  if (p_EqualPolys(I->m[i],pt,r))
831  {
832  p_Delete(&pt,r);
833  return true;
834  }
835  }
836  p_Delete(&pt,r);
837  return false;
838 }
839 
840 int tropicalStrategy::findPositionOfUniformizingBinomial(const ideal I, const ring r) const
841 {
843 
844  // if the valuation is non-trivial then checks if the first generator is p-t
845  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
846  poly p = p_One(r);
847  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
848  poly t = p_One(r);
849  p_SetExp(t,1,1,r);
850  p_Setm(t,r);
851  poly pt = p_Add_q(p,p_Neg(t,r),r);
852 
853  for (int i=0; i<IDELEMS(I); i++)
854  {
855  if (p_EqualPolys(I->m[i],pt,r))
856  {
857  p_Delete(&pt,r);
858  return i;
859  }
860  }
861  p_Delete(&pt,r);
862  return -1;
863 }
864 
865 bool tropicalStrategy::checkForUniformizingParameter(const ideal inI, const ring r) const
866 {
867  // if the valuation is trivial,
868  // then there is no special condition the first generator has to fullfill
869  if (isValuationTrivial())
870  return true;
871 
872  // if the valuation is non-trivial then checks if the first generator is p
873  if (inI->m[0]==NULL)
874  return false;
875  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
876  poly p = p_One(r);
877  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
878 
879  for (int i=0; i<IDELEMS(inI); i++)
880  {
881  if (p_EqualPolys(inI->m[i],p,r))
882  {
883  p_Delete(&p,r);
884  return true;
885  }
886  }
887  p_Delete(&p,r);
888  return false;
889 }
gfan::ZVector nonvalued_adjustWeightForHomogeneity(const gfan::ZVector &w)
gfan::ZVector nonvalued_adjustWeightUnderHomogeneity(const gfan::ZVector &e, const gfan::ZVector &)
gfan::ZVector valued_adjustWeightForHomogeneity(const gfan::ZVector &w)
gfan::ZVector valued_adjustWeightUnderHomogeneity(const gfan::ZVector &e, const gfan::ZVector &w)
#define FALSE
Definition: auxiliary.h:96
BOOLEAN linealitySpace(leftv res, leftv args)
Definition: bbcone.cc:945
int * ZVectorToIntStar(const gfan::ZVector &v, bool &overflow)
int l
Definition: cfEzgcd.cc:100
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
int p
Definition: cfModGcd.cc:4078
return false
Definition: cfModGcd.cc:84
g
Definition: cfModGcd.cc:4090
CanonicalForm b
Definition: cfModGcd.cc:4103
poly * m
Definition: matpol.h:18
int expectedDimension
the expected Dimension of the polyhedral output, i.e.
bool isValuationTrivial() const
ideal getOriginalIdeal() const
returns the input ideal over the field with valuation
tropicalStrategy(const ideal I, const ring r, const bool completelyHomogeneous=true, const bool completeSpace=true)
Constructor for the trivial valuation case.
bool isValuationNonTrivial() const
std::pair< ideal, ring > computeFlip(const ideal Ir, const ring r, const gfan::ZVector &interiorPoint, const gfan::ZVector &facetNormal) const
given an interior point of a groebner cone computes the groebner cone adjacent to it
tropicalStrategy & operator=(const tropicalStrategy &currentStrategy)
assignment operator
ring copyAndChangeOrderingLS(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
void putUniformizingBinomialInFront(ideal I, const ring r, const number q) const
gfan::ZVector adjustWeightUnderHomogeneity(gfan::ZVector v, gfan::ZVector w) const
Given strictly positive weight w and weight v, returns a strictly positive weight u such that on an i...
bool reduce(ideal I, const ring r) const
reduces the generators of an ideal I so that the inequalities and equations of the Groebner cone can ...
bool onlyLowerHalfSpace
true if valuation non-trivial, false otherwise
gfan::ZCone linealitySpace
the homogeneity space of the Grobner fan
int getExpectedDimension() const
returns the expected Dimension of the polyhedral output
ring startingRing
polynomial ring over the valuation ring extended by one extra variable t
ideal originalIdeal
input ideal, assumed to be a homogeneous prime ideal
gfan::ZVector(* weightAdjustingAlgorithm1)(const gfan::ZVector &w)
A function such that: Given weight w, returns a strictly positive weight u such that an ideal satisfy...
void pReduce(ideal I, const ring r) const
~tropicalStrategy()
destructor
int findPositionOfUniformizingBinomial(const ideal I, const ring r) const
ideal computeWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
suppose w a weight in maximal groebner cone of > suppose I (initially) reduced standard basis w....
ring shortcutRing
polynomial ring over the residue field
bool(* extraReductionAlgorithm)(ideal I, ring r, number p)
A function that reduces the generators of an ideal I so that the inequalities and equations of the Gr...
ring getStartingRing() const
returns the polynomial ring over the valuation ring
gfan::ZVector adjustWeightForHomogeneity(gfan::ZVector w) const
Given weight w, returns a strictly positive weight u such that an ideal satisfying the valuation-sepc...
ring getShortcutRingPrependingWeight(const ring r, const gfan::ZVector &w) const
If valuation trivial, returns a copy of r with a positive weight prepended, such that any ideal homog...
number uniformizingParameter
uniformizing parameter in the valuation ring
ring copyAndChangeCoefficientRing(const ring r) const
ring copyAndChangeOrderingWP(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
ideal computeLift(const ideal inJs, const ring s, const ideal inIr, const ideal Ir, const ring r) const
ideal startingIdeal
preimage of the input ideal under the map that sends t to the uniformizing parameter
bool checkForUniformizingParameter(const ideal inI, const ring r) const
if valuation non-trivial, checks whether the genearting system contains p otherwise returns true
ideal getStartingIdeal() const
returns the input ideal
bool restrictToLowerHalfSpace() const
returns true, if valuation non-trivial, false otherwise
gfan::ZVector(* weightAdjustingAlgorithm2)(const gfan::ZVector &v, const gfan::ZVector &w)
A function such that: Given strictly positive weight w and weight v, returns a strictly positive weig...
ideal computeStdOfInitialIdeal(const ideal inI, const ring r) const
given generators of the initial ideal, computes its standard basis
ring getOriginalRing() const
returns the polynomial ring over the field with valuation
number getUniformizingParameter() const
returns the uniformizing parameter in the valuation ring
std::pair< poly, int > checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector &w=0) const
If given w, assuming w is in the Groebner cone of the ordering on r and I is a standard basis with re...
bool checkForUniformizingBinomial(const ideal I, const ring r) const
if valuation non-trivial, checks whether the generating system contains p-t otherwise returns true
ring getShortcutRing() const
ring originalRing
polynomial ring over a field with valuation
static FORCE_INLINE long n_Int(number &n, const coeffs r)
conversion of n to an int; 0 if not possible in Z/pZ: the representing int lying in (-p/2 ....
Definition: coeffs.h:547
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:451
#define n_Test(a, r)
BOOLEAN n_Test(number a, const coeffs r)
Definition: coeffs.h:712
@ n_Zp
\F{p < 2^31}
Definition: coeffs.h:29
@ n_Z
only used if HAVE_RINGS is defined
Definition: coeffs.h:43
static FORCE_INLINE BOOLEAN n_IsUnit(number n, const coeffs r)
TRUE iff n has a multiplicative inverse in the given coeff field/ring r.
Definition: coeffs.h:515
static FORCE_INLINE nMapFunc n_SetMap(const coeffs src, const coeffs dst)
set the mapping function pointers for translating numbers from src to dst
Definition: coeffs.h:700
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:354
static FORCE_INLINE coeffs nCopyCoeff(const coeffs r)
"copy" coeffs, i.e. increment ref
Definition: coeffs.h:429
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:455
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:538
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:73
void nKillChar(coeffs r)
undo all initialisations
Definition: numbers.cc:522
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition: coeffs.h:468
poly searchForMonomialViaStepwiseSaturation(const ideal I, const ring r, const gfan::ZVector w0)
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm & w
Definition: facAbsFact.cc:51
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
int scDimInt(ideal S, ideal Q)
ideal dimension
Definition: hdegree.cc:77
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
ideal id_Copy(ideal h1, const ring r)
copy an ideal
#define idPosConstant(I)
index of generator with leading term in ground ring (if any); otherwise -1
Definition: ideals.h:37
poly initial(const poly p, const ring r, const gfan::ZVector &w)
Returns the initial form of p with respect to w.
Definition: initial.cc:30
STATIC_VAR Poly * h
Definition: janet.cc:971
STATIC_VAR jList * Q
Definition: janet.cc:30
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:3167
void mp_Delete(matrix *a, const ring r)
Definition: matpol.cc:880
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define assume(x)
Definition: mod2.h:387
#define pNext(p)
Definition: monomials.h:36
#define p_GetCoeff(p, r)
Definition: monomials.h:50
#define omStrDup(s)
Definition: omAllocDecl.h:263
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omFree(addr)
Definition: omAllocDecl.h:261
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define NULL
Definition: omList.c:12
poly p_PermPoly(poly p, const int *perm, const ring oldRing, const ring dst, nMapFunc nMap, const int *par_perm, int OldPar, BOOLEAN use_mult)
Definition: p_polys.cc:4163
poly p_Div_nn(poly p, const number n, const ring r)
Definition: p_polys.cc:1497
poly p_One(const ring r)
Definition: p_polys.cc:1309
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition: p_polys.cc:4545
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1079
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:908
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1086
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent @Note: VarOffset encodes the position in p->exp
Definition: p_polys.h:488
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:233
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:412
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent @Note: the integer VarOffset encodes:
Definition: p_polys.h:469
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:873
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:818
void rChangeCurrRing(ring r)
Definition: polys.cc:15
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define pDelete(p_ptr)
Definition: polys.h:186
bool isOrderingLocalInT(const ring r)
bool ppreduceInitially(poly *hStar, const poly g, const ring r)
reduces h initially with respect to g, returns false if h was initially reduced in the first place,...
int IsPrime(int p)
Definition: prime.cc:61
BOOLEAN rComplete(ring r, int force)
this needs to be called whenever a new ring is created: new fields in ring are created (like VarOffse...
Definition: ring.cc:3395
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1363
void rDelete(ring r)
unconditionally deletes fields in r
Definition: ring.cc:449
static int sign(int x)
Definition: ring.cc:3372
ring rCopy(ring r)
Definition: ring.cc:1645
static int rBlocks(ring r)
Definition: ring.h:569
static BOOLEAN rField_is_Z(const ring r)
Definition: ring.h:510
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:501
rRingOrder_t
order stuff
Definition: ring.h:68
@ ringorder_lp
Definition: ring.h:77
@ ringorder_a
Definition: ring.h:70
@ ringorder_C
Definition: ring.h:73
@ ringorder_ds
Definition: ring.h:84
@ ringorder_dp
Definition: ring.h:78
@ ringorder_ws
Definition: ring.h:86
@ ringorder_ls
Definition: ring.h:83
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:507
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:593
#define rTest(r)
Definition: ring.h:786
#define rField_is_Ring(R)
Definition: ring.h:486
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
ideal id_Head(ideal h, const ring r)
returns the ideals of initial terms
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
#define IDELEMS(i)
Definition: simpleideals.h:23
#define id_Test(A, lR)
Definition: simpleideals.h:78
ideal gfanlib_kStd_wrapper(ideal I, ring r, tHomog h=testHomog)
Definition: std_wrapper.cc:6
bool checkWeightVector(const ideal I, const ring r, const gfan::ZVector &weightVector, bool checkBorder)
bool checkForNonPositiveEntries(const gfan::ZVector &w)
bool areIdealsEqual(ideal I, ring r, ideal J, ring s)
static bool noExtraReduction(ideal I, ring r, number)
int dim(ideal I, ring r)
static ideal constructStartingIdeal(ideal originalIdeal, ring originalRing, number uniformizingParameter, ring startingRing)
static ring constructStartingRing(ring r)
Given a polynomial ring r over the rational numbers and a weighted ordering, returns a polynomial rin...
static void swapElements(ideal I, ideal J)
implementation of the class tropicalStrategy
gfan::ZCone homogeneitySpace(ideal I, ring r)
Definition: tropical.cc:19
matrix divisionDiscardingRemainder(const poly f, const ideal G, const ring r)
Computes a division discarding remainder of f with respect to G.
Definition: witness.cc:9
poly witness(const poly m, const ideal I, const ideal inI, const ring r)
Let w be the uppermost weight vector in the matrix defining the ordering on r.
Definition: witness.cc:34