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 
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)),
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),
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())),
325  onlyLowerHalfSpace(currentStrategy.restrictToLowerHalfSpace()),
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 
378  if (originalIdeal) id_Test(originalIdeal,originalRing);
380  if (startingIdeal) id_Test(startingIdeal,startingRing);
381  if (uniformizingParameter) n_Test(uniformizingParameter,startingRing->cf);
382  if (shortcutRing) rTest(shortcutRing);
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 }
#define idPosConstant(I)
index of generator with leading term in ground ring (if any); otherwise -1
Definition: ideals.h:37
void putUniformizingBinomialInFront(ideal I, const ring r, const number q) const
bool checkForNonPositiveEntries(const gfan::ZVector &w)
implementation of the class tropicalStrategy
ring getShortcutRing() const
const CanonicalForm int s
Definition: facAbsFact.cc:55
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...
int j
Definition: facHensel.cc:105
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 kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:2822
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
gfan::ZCone homogeneitySpace(ideal I, ring r)
Definition: tropical.cc:19
gfan::ZVector nonvalued_adjustWeightUnderHomogeneity(const gfan::ZVector &e, const gfan::ZVector &)
int findPositionOfUniformizingBinomial(const ideal I, const ring r) const
ring copyAndChangeCoefficientRing(const ring r) const
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
#define FALSE
Definition: auxiliary.h:94
ring getOriginalRing() const
returns the polynomial ring over the field with valuation
bool isOrderingLocalInT(const ring r)
ideal id_Copy(ideal h1, const ring r)
copy an ideal
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff &#39;n&#39; represents the one element.
Definition: coeffs.h:468
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...
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...
#define id_Test(A, lR)
Definition: simpleideals.h:79
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
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
{p < 2^31}
Definition: coeffs.h:30
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:586
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
ring getStartingRing() const
returns the polynomial ring over the valuation ring
poly p_Div_nn(poly p, const number n, const ring r)
Definition: p_polys.cc:1487
static ideal constructStartingIdeal(ideal originalIdeal, ring originalRing, number uniformizingParameter, ring startingRing)
g
Definition: cfModGcd.cc:4031
bool onlyLowerHalfSpace
true if valuation non-trivial, false otherwise
int * ZVectorToIntStar(const gfan::ZVector &v, bool &overflow)
int k
Definition: cfEzgcd.cc:92
static BOOLEAN rField_is_Z(const ring r)
Definition: ring.h:504
#define Q
Definition: sirandom.c:25
number getUniformizingParameter() const
returns the uniformizing parameter in the valuation ring
#define omAlloc(size)
Definition: omAllocDecl.h:210
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:411
static void swapElements(ideal I, ideal J)
static ring constructStartingRing(ring r)
Given a polynomial ring r over the rational numbers and a weighted ordering, returns a polynomial rin...
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:811
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 ...
ideal startingIdeal
preimage of the input ideal under the map that sends t to the uniformizing parameter ...
~tropicalStrategy()
destructor
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
poly * m
Definition: matpol.h:18
CanonicalForm b
Definition: cfModGcd.cc:4044
static int rBlocks(ring r)
Definition: ring.h:562
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:4036
BOOLEAN linealitySpace(leftv res, leftv args)
Definition: bbcone.cc:945
ideal originalIdeal
input ideal, assumed to be a homogeneous prime ideal
static bool noExtraReduction(ideal I, ring r, number)
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
poly p_One(const ring r)
Definition: p_polys.cc:1303
bool isValuationTrivial() const
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:3394
ring shortcutRing
polynomial ring over the residue field
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent : the integer VarOffset encodes:
Definition: p_polys.h:468
ideal gfanlib_kStd_wrapper(ideal I, ring r, tHomog h=testHomog)
Definition: std_wrapper.cc:6
#define omFree(addr)
Definition: omAllocDecl.h:261
#define assume(x)
Definition: mod2.h:390
gfan::ZCone getHomogeneitySpace() const
returns the homogeneity space of the preimage ideal
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1365
ideal getStartingIdeal() const
returns the input ideal
int scDimInt(ideal S, ideal Q)
Definition: hdegree.cc:71
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:73
ring copyAndChangeOrderingLS(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
rRingOrder_t
order stuff
Definition: ring.h:67
bool areIdealsEqual(ideal I, ring r, ideal J, ring s)
#define n_Test(a, r)
BOOLEAN n_Test(number a, const coeffs r)
Definition: coeffs.h:738
#define rTest(r)
Definition: ring.h:780
return false
Definition: cfModGcd.cc:84
void pReduce(ideal I, const ring r) const
bool checkForUniformizingParameter(const ideal inI, const ring r) const
if valuation non-trivial, checks whether the genearting system contains p otherwise returns true ...
only used if HAVE_RINGS is defined
Definition: coeffs.h:44
poly searchForMonomialViaStepwiseSaturation(const ideal I, const ring r, const gfan::ZVector w0)
int dim(ideal I, ring r)
gfan::ZVector adjustWeightForHomogeneity(gfan::ZVector w) const
Given weight w, returns a strictly positive weight u such that an ideal satisfying the valuation-sepc...
bool restrictToLowerHalfSpace() const
returns true, if valuation non-trivial, false otherwise
int i
Definition: cfEzgcd.cc:125
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:501
int IsPrime(int p)
Definition: prime.cc:61
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...
#define IDELEMS(i)
Definition: simpleideals.h:23
void mp_Delete(matrix *a, const ring r)
Definition: matpol.cc:880
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:721
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
ideal computeLift(const ideal inJs, const ring s, const ideal inIr, const ideal Ir, const ring r) const
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition: p_polys.cc:4418
gfan::ZVector valued_adjustWeightUnderHomogeneity(const gfan::ZVector &e, const gfan::ZVector &w)
void rChangeCurrRing(ring r)
Definition: polys.cc:14
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:495
ring copyAndChangeOrderingWP(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
bool checkForUniformizingBinomial(const ideal I, const ring r) const
if valuation non-trivial, checks whether the generating system contains p-t otherwise returns true ...
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
static FORCE_INLINE coeffs nCopyCoeff(const coeffs r)
"copy" coeffs, i.e. increment ref
Definition: coeffs.h:429
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:856
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
number uniformizingParameter
uniformizing parameter in the valuation ring
ring rCopy(ring r)
Definition: ring.cc:1645
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent : VarOffset encodes the position in p->exp
Definition: p_polys.h:487
ring startingRing
polynomial ring over the valuation ring extended by one extra variable t
bool isValuationNonTrivial() const
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 ...
static BOOLEAN rField_is_Ring(const ring r)
Definition: ring.h:479
#define NULL
Definition: omList.c:12
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
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...
ideal getOriginalIdeal() const
returns the input ideal over the field with valuation
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of &#39;n&#39;
Definition: coeffs.h:451
void rDelete(ring r)
unconditionally deletes fields in r
Definition: ring.cc:448
gfan::ZVector nonvalued_adjustWeightForHomogeneity(const gfan::ZVector &w)
gfan::ZVector valued_adjustWeightForHomogeneity(const gfan::ZVector &w)
ideal id_Head(ideal h, const ring r)
returns the ideals of initial terms
int getExpectedDimension() const
returns the expected Dimension of the polyhedral output
const CanonicalForm & w
Definition: facAbsFact.cc:55
#define pDelete(p_ptr)
Definition: polys.h:181
#define pNext(p)
Definition: monomials.h:36
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
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...
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:232
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 checkWeightVector(const ideal I, const ring r, const gfan::ZVector &weightVector, bool checkBorder)
#define p_GetCoeff(p, r)
Definition: monomials.h:50
tropicalStrategy(const ideal I, const ring r, const bool completelyHomogeneous=true, const bool completeSpace=true)
Constructor for the trivial valuation case.
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1042
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete &#39;p&#39;
Definition: coeffs.h:455
ring originalRing
polynomial ring over a field with valuation
int p
Definition: cfModGcd.cc:4019
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:891
static Poly * h
Definition: janet.cc:971
tropicalStrategy & operator=(const tropicalStrategy &currentStrategy)
assignment operator
gfan::ZCone linealitySpace
the homogeneity space of the Grobner fan
void nKillChar(coeffs r)
undo all initialisations
Definition: numbers.cc:510
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1049
static int sign(int x)
Definition: ring.cc:3371
#define omAlloc0(size)
Definition: omAllocDecl.h:211
int l
Definition: cfEzgcd.cc:93
int expectedDimension
the expected Dimension of the polyhedral output, i.e.
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...
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:349
ideal computeStdOfInitialIdeal(const ideal inI, const ring r) const
given generators of the initial ideal, computes its standard basis
#define omStrDup(s)
Definition: omAllocDecl.h:263