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_Ring_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  int* 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 = (int*) omAlloc0((h+1)*sizeof(int));
79  rShortcut->block0 = (int*) omAlloc0((h+1)*sizeof(int));
80  rShortcut->block1 = (int*) omAlloc0((h+1)*sizeof(int));
81  rShortcut->wvhdl = (int**) omAlloc0((h+1)*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 
95  rComplete(rShortcut);
96  rTest(rShortcut);
97 
98  omFree(order);
99  omFree(block0);
100  omFree(block1);
101  omFree(wvhdl);
102 
103  int k = IDELEMS(I);
104  ideal IShortcut = idInit(k);
105  nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
106  for (int i=0; i<k; i++)
107  {
108  if(I->m[i]!=NULL)
109  {
110  IShortcut->m[i] = p_PermPoly(I->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
111  }
112  }
113 
114  ideal JShortcut = gfanlib_kStd_wrapper(IShortcut,rShortcut);
115 
116  ideal J = idInit(k);
117  nMapFunc outofShortcut = n_SetMap(rShortcut->cf,r->cf);
118  for (int i=0; i<k; i++)
119  J->m[i] = p_PermPoly(JShortcut->m[i],NULL,rShortcut,r,outofShortcut,NULL,0);
120 
121  assume(areIdealsEqual(J,r,I,r));
122  swapElements(I,J);
123  id_Delete(&IShortcut,rShortcut);
124  id_Delete(&JShortcut,rShortcut);
125  rDelete(rShortcut);
126  id_Delete(&J,r);
127  return false;
128 }
129 
130 /**
131  * Initializes all relevant structures and information for the trivial valuation case,
132  * i.e. computing a tropical variety without any valuation.
133  */
134 tropicalStrategy::tropicalStrategy(const ideal I, const ring r,
135  const bool completelyHomogeneous,
136  const bool completeSpace):
137  originalRing(rCopy(r)),
138  originalIdeal(id_Copy(I,r)),
139  expectedDimension(dim(originalIdeal,originalRing)),
140  linealitySpace(homogeneitySpace(originalIdeal,originalRing)),
141  startingRing(rCopy(originalRing)),
142  startingIdeal(id_Copy(originalIdeal,originalRing)),
143  uniformizingParameter(NULL),
144  shortcutRing(NULL),
145  onlyLowerHalfSpace(false),
146  weightAdjustingAlgorithm1(nonvalued_adjustWeightForHomogeneity),
147  weightAdjustingAlgorithm2(nonvalued_adjustWeightUnderHomogeneity),
148  extraReductionAlgorithm(noExtraReduction)
149 {
151  if (!completelyHomogeneous)
152  {
155  }
156  if (!completeSpace)
157  onlyLowerHalfSpace = true;
158 }
159 
160 /**
161  * Given a polynomial ring r over the rational numbers and a weighted ordering,
162  * returns a polynomial ring s over the integers with one extra variable, which is weighted -1.
163  */
164 static ring constructStartingRing(ring r)
165 {
166  assume(rField_is_Q(r));
167 
168  ring s = rCopy0(r,FALSE,FALSE);
169  nKillChar(s->cf);
170  s->cf = nInitChar(n_Z,NULL);
171 
172  int n = rVar(s)+1;
173  s->N = n;
174  char** oldNames = s->names;
175  s->names = (char**) omAlloc((n+1)*sizeof(char**));
176  s->names[0] = omStrDup("t");
177  for (int i=1; i<n; i++)
178  s->names[i] = oldNames[i-1];
179  omFree(oldNames);
180 
181  s->order = (int*) omAlloc0(3*sizeof(int));
182  s->block0 = (int*) omAlloc0(3*sizeof(int));
183  s->block1 = (int*) omAlloc0(3*sizeof(int));
184  s->wvhdl = (int**) omAlloc0(3*sizeof(int**));
185  s->order[0] = ringorder_ws;
186  s->block0[0] = 1;
187  s->block1[0] = n;
188  s->wvhdl[0] = (int*) omAlloc(n*sizeof(int));
189  s->wvhdl[0][0] = 1;
190  if (r->order[0] == ringorder_lp)
191  {
192  s->wvhdl[0][1] = 1;
193  }
194  else if (r->order[0] == ringorder_ls)
195  {
196  s->wvhdl[0][1] = -1;
197  }
198  else if (r->order[0] == ringorder_dp)
199  {
200  for (int i=1; i<n; i++)
201  s->wvhdl[0][i] = -1;
202  }
203  else if (r->order[0] == ringorder_ds)
204  {
205  for (int i=1; i<n; i++)
206  s->wvhdl[0][i] = 1;
207  }
208  else if (r->order[0] == ringorder_ws)
209  {
210  for (int i=1; i<n; i++)
211  s->wvhdl[0][i] = r->wvhdl[0][i-1];
212  }
213  else
214  {
215  for (int i=1; i<n; i++)
216  s->wvhdl[0][i] = -r->wvhdl[0][i-1];
217  }
218  s->order[1] = ringorder_C;
219 
220  rComplete(s);
221  rTest(s);
222  return s;
223 }
224 
226 {
227  // construct p-t
228  poly g = p_One(startingRing);
229  p_SetCoeff(g,uniformizingParameter,startingRing);
230  pNext(g) = p_One(startingRing);
231  p_SetExp(pNext(g),1,1,startingRing);
232  p_SetCoeff(pNext(g),n_Init(-1,startingRing->cf),startingRing);
233  p_Setm(pNext(g),startingRing);
234  ideal pt = idInit(1);
235  pt->m[0] = g;
236 
237  // map originalIdeal from originalRing into startingRing
238  int k = IDELEMS(originalIdeal);
239  ideal J = idInit(k+1);
240  nMapFunc nMap = n_SetMap(originalRing->cf,startingRing->cf);
241  int n = rVar(originalRing);
242  int* shiftByOne = (int*) omAlloc((n+1)*sizeof(int));
243  for (int i=1; i<=n; i++)
244  shiftByOne[i]=i+1;
245  for (int i=0; i<k; i++)
246  {
247  if(originalIdeal->m[i]!=NULL)
248  {
249  J->m[i] = p_PermPoly(originalIdeal->m[i],shiftByOne,originalRing,startingRing,nMap,NULL,0);
250  }
251  }
252  omFreeSize(shiftByOne,(n+1)*sizeof(int));
253 
254  ring origin = currRing;
255  rChangeCurrRing(startingRing);
256  ideal startingIdeal = kNF(pt,startingRing->qideal,J); // mathematically redundant,
257  rChangeCurrRing(origin); // but helps with upcoming std computation
258  // ideal startingIdeal = J; J = NULL;
259  assume(startingIdeal->m[k]==NULL);
260  startingIdeal->m[k] = pt->m[0];
261  startingIdeal = gfanlib_kStd_wrapper(startingIdeal,startingRing);
262 
263  id_Delete(&J,startingRing);
264  pt->m[0] = NULL;
265  id_Delete(&pt,startingRing);
266  return startingIdeal;
267 }
268 
269 /***
270  * Initializes all relevant structures and information for the valued case,
271  * i.e. computing a tropical variety over the rational numbers with p-adic valuation
272  **/
273 tropicalStrategy::tropicalStrategy(ideal J, number q, ring s):
274  originalRing(rCopy(s)),
275  originalIdeal(id_Copy(J,s)),
277  linealitySpace(gfan::ZCone()), // to come, see below
278  startingRing(NULL), // to come, see below
279  startingIdeal(NULL), // to come, see below
280  uniformizingParameter(NULL), // to come, see below
281  shortcutRing(NULL), // to come, see below
282  onlyLowerHalfSpace(true),
286 {
287  /* assume that the ground field of the originalRing is Q */
288  assume(rField_is_Q(s));
289 
290  /* replace Q with Z for the startingRing
291  * and add an extra variable for tracking the uniformizing parameter */
293 
294  /* map the uniformizing parameter into the new coefficient domain */
295  nMapFunc nMap = n_SetMap(originalRing->cf,startingRing->cf);
297 
298  /* map the input ideal into the new polynomial ring */
301 
303 
304  /* construct the shorcut ring */
306  nKillChar(shortcutRing->cf);
310 }
311 
313  originalRing(rCopy(currentStrategy.getOriginalRing())),
314  originalIdeal(id_Copy(currentStrategy.getOriginalIdeal(),currentStrategy.getOriginalRing())),
315  expectedDimension(currentStrategy.getExpectedDimension()),
316  linealitySpace(currentStrategy.getHomogeneitySpace()),
317  startingRing(rCopy(currentStrategy.getStartingRing())),
318  startingIdeal(id_Copy(currentStrategy.getStartingIdeal(),currentStrategy.getStartingRing())),
321  onlyLowerHalfSpace(currentStrategy.restrictToLowerHalfSpace()),
325 {
330  if (currentStrategy.getUniformizingParameter())
331  {
334  }
335  if (currentStrategy.getShortcutRing())
336  {
337  shortcutRing = rCopy(currentStrategy.getShortcutRing());
339  }
340 }
341 
343 {
350 
357 }
358 
360 {
361  originalRing = rCopy(currentStrategy.getOriginalRing());
362  originalIdeal = id_Copy(currentStrategy.getOriginalIdeal(),currentStrategy.getOriginalRing());
363  expectedDimension = currentStrategy.getExpectedDimension();
364  startingRing = rCopy(currentStrategy.getStartingRing());
365  startingIdeal = id_Copy(currentStrategy.getStartingIdeal(),currentStrategy.getStartingRing());
367  shortcutRing = rCopy(currentStrategy.getShortcutRing());
368  onlyLowerHalfSpace = currentStrategy.restrictToLowerHalfSpace();
372 
374  if (originalIdeal) id_Test(originalIdeal,originalRing);
376  if (startingIdeal) id_Test(startingIdeal,startingRing);
377  if (uniformizingParameter) n_Test(uniformizingParameter,startingRing->cf);
378  if (shortcutRing) rTest(shortcutRing);
379 
380  return *this;
381 }
382 
383 void tropicalStrategy::putUniformizingBinomialInFront(ideal I, const ring r, const number q) const
384 {
385  poly p = p_One(r);
386  p_SetCoeff(p,q,r);
387  poly t = p_One(r);
388  p_SetExp(t,1,1,r);
389  p_Setm(t,r);
390  poly pt = p_Add_q(p,p_Neg(t,r),r);
391 
392  int k = IDELEMS(I);
393  int l;
394  for (l=0; l<k; l++)
395  {
396  if (p_EqualPolys(I->m[l],pt,r))
397  break;
398  }
399  p_Delete(&pt,r);
400 
401  if (l>1)
402  {
403  pt = I->m[l];
404  for (int i=l; i>0; i--)
405  I->m[l] = I->m[l-1];
406  I->m[0] = pt;
407  pt = NULL;
408  }
409  return;
410 }
411 
412 bool tropicalStrategy::reduce(ideal I, const ring r) const
413 {
414  rTest(r);
415  id_Test(I,r);
416 
417  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
418  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
419  bool b = extraReductionAlgorithm(I,r,p);
420  // putUniformizingBinomialInFront(I,r,p);
421  n_Delete(&p,r->cf);
422 
423  return b;
424 }
425 
426 void tropicalStrategy::pReduce(ideal I, const ring r) const
427 {
428  rTest(r);
429  id_Test(I,r);
430 
431  if (isValuationTrivial())
432  return;
433 
434  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
435  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
436  ::pReduce(I,p,r);
437  n_Delete(&p,r->cf);
438 
439  return;
440 }
441 
442 ring tropicalStrategy::getShortcutRingPrependingWeight(const ring r, const gfan::ZVector &v) const
443 {
444  ring rShortcut = rCopy0(r);
445 
446  // save old ordering
447  int* order = rShortcut->order;
448  int* block0 = rShortcut->block0;
449  int* block1 = rShortcut->block1;
450  int** wvhdl = rShortcut->wvhdl;
451 
452  // adjust weight and create new ordering
453  gfan::ZVector w = adjustWeightForHomogeneity(v);
454  int h = rBlocks(r); int n = rVar(r);
455  rShortcut->order = (int*) omAlloc0((h+1)*sizeof(int));
456  rShortcut->block0 = (int*) omAlloc0((h+1)*sizeof(int));
457  rShortcut->block1 = (int*) omAlloc0((h+1)*sizeof(int));
458  rShortcut->wvhdl = (int**) omAlloc0((h+1)*sizeof(int*));
459  rShortcut->order[0] = ringorder_a;
460  rShortcut->block0[0] = 1;
461  rShortcut->block1[0] = n;
462  bool overflow;
463  rShortcut->wvhdl[0] = ZVectorToIntStar(w,overflow);
464  for (int i=1; i<=h; i++)
465  {
466  rShortcut->order[i] = order[i-1];
467  rShortcut->block0[i] = block0[i-1];
468  rShortcut->block1[i] = block1[i-1];
469  rShortcut->wvhdl[i] = wvhdl[i-1];
470  }
471 
472  // if valuation non-trivial, change coefficient ring to residue field
473  if (isValuationNonTrivial())
474  {
475  nKillChar(rShortcut->cf);
476  rShortcut->cf = nCopyCoeff(shortcutRing->cf);
477  }
478  rComplete(rShortcut);
479  rTest(rShortcut);
480 
481  // delete old ordering
482  omFree(order);
483  omFree(block0);
484  omFree(block1);
485  omFree(wvhdl);
486 
487  return rShortcut;
488 }
489 
490 std::pair<poly,int> tropicalStrategy::checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector &w) const
491 {
492  // quick check whether I already contains an ideal
493  int k = IDELEMS(I);
494  for (int i=0; i<k; i++)
495  {
496  poly g = I->m[i];
497  if (g!=NULL
498  && pNext(g)==NULL
499  && (isValuationTrivial() || n_IsOne(p_GetCoeff(g,r),r->cf)))
500  return std::pair<poly,int>(g,i);
501  }
502 
503  ring rShortcut;
504  ideal inIShortcut;
505  if (w.size()>0)
506  {
507  // if needed, prepend extra weight for homogeneity
508  // switch to residue field if valuation is non trivial
509  rShortcut = getShortcutRingPrependingWeight(r,w);
510 
511  // compute the initial ideal and map it into the constructed ring
512  // if switched to residue field, remove possibly 0 elements
513  ideal inI = initial(I,r,w);
514  inIShortcut = idInit(k);
515  nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
516  for (int i=0; i<k; i++)
517  inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
518  if (isValuationNonTrivial())
519  idSkipZeroes(inIShortcut);
520  id_Delete(&inI,r);
521  }
522  else
523  {
524  rShortcut = r;
525  inIShortcut = I;
526  }
527 
528  gfan::ZCone C0 = homogeneitySpace(inIShortcut,rShortcut);
529  gfan::ZCone pos = gfan::ZCone::positiveOrthant(C0.ambientDimension());
530  gfan::ZCone C0pos = intersection(C0,pos);
531  C0pos.canonicalize();
532  gfan::ZVector wpos = C0pos.getRelativeInteriorPoint();
534 
535  // check initial ideal for monomial and
536  // if it exsists, return a copy of the monomial in the input ring
537  poly p = searchForMonomialViaStepwiseSaturation(inIShortcut,rShortcut,wpos);
538  poly monomial = NULL;
539  if (p!=NULL)
540  {
541  monomial=p_One(r);
542  for (int i=1; i<=rVar(r); i++)
543  p_SetExp(monomial,i,p_GetExp(p,i,rShortcut),r);
544  p_Setm(monomial,r);
545  p_Delete(&p,rShortcut);
546  }
547 
548 
549  if (w.size()>0)
550  {
551  // if needed, cleanup
552  id_Delete(&inIShortcut,rShortcut);
553  rDelete(rShortcut);
554  }
555  return std::pair<poly,int>(monomial,-1);
556 }
557 
559 {
560  ring rShortcut = rCopy0(r);
561  nKillChar(rShortcut->cf);
562  rShortcut->cf = nCopyCoeff(shortcutRing->cf);
563  rComplete(rShortcut);
564  rTest(rShortcut);
565  return rShortcut;
566 }
567 
568 ideal tropicalStrategy::computeWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
569 {
570  // if the valuation is trivial and the ring and ideal have not been extended,
571  // then it is sufficient to return the difference between the elements of inJ
572  // and their normal forms with respect to I and r
573  if (isValuationTrivial())
574  return witness(inJ,I,r);
575  // if the valuation is non-trivial and the ring and ideal have been extended,
576  // then we can make a shortcut through the residue field
577  else
578  {
579  assume(IDELEMS(inI)==IDELEMS(I));
580  int uni = findPositionOfUniformizingBinomial(I,r);
581  assume(uni>=0);
582  /**
583  * change ground ring into finite field
584  * and map the data into it
585  */
586  ring rShortcut = copyAndChangeCoefficientRing(r);
587 
588  int k = IDELEMS(inJ);
589  int l = IDELEMS(I);
590  ideal inJShortcut = idInit(k);
591  ideal inIShortcut = idInit(l);
592  nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
593  for (int i=0; i<k; i++)
594  inJShortcut->m[i] = p_PermPoly(inJ->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
595  for (int j=0; j<l; j++)
596  inIShortcut->m[j] = p_PermPoly(inI->m[j],NULL,r,rShortcut,takingResidues,NULL,0);
597  id_Test(inJShortcut,rShortcut);
598  id_Test(inIShortcut,rShortcut);
599 
600  /**
601  * Compute a division with remainder over the finite field
602  * and map the result back to r
603  */
604  matrix QShortcut = divisionDiscardingRemainder(inJShortcut,inIShortcut,rShortcut);
605  matrix Q = mpNew(l,k);
606  nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
607  for (int ij=k*l-1; ij>=0; ij--)
608  Q->m[ij] = p_PermPoly(QShortcut->m[ij],NULL,rShortcut,r,takingRepresentatives,NULL,0);
609 
610  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
611  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
612 
613  /**
614  * Compute the normal forms
615  */
616  ideal J = idInit(k);
617  for (int j=0; j<k; j++)
618  {
619  poly q0 = p_Copy(inJ->m[j],r);
620  for (int i=0; i<l; i++)
621  {
622  poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
623  poly inIi = p_Copy(inI->m[i],r);
624  q0 = p_Add_q(q0,p_Neg(p_Mult_q(qij,inIi,r),r),r);
625  }
626  q0 = p_Div_nn(q0,p,r);
627  poly q0g0 = p_Mult_q(q0,p_Copy(I->m[uni],r),r);
628  // q0 = NULL;
629  poly qigi = NULL;
630  for (int i=0; i<l; i++)
631  {
632  poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
633  // poly inIi = p_Copy(I->m[i],r);
634  poly Ii = p_Copy(I->m[i],r);
635  qigi = p_Add_q(qigi,p_Mult_q(qij,Ii,r),r);
636  }
637  J->m[j] = p_Add_q(q0g0,qigi,r);
638  }
639 
640  id_Delete(&inIShortcut,rShortcut);
641  id_Delete(&inJShortcut,rShortcut);
642  mp_Delete(&QShortcut,rShortcut);
643  rDelete(rShortcut);
644  mp_Delete(&Q,r);
645  n_Delete(&p,r->cf);
646  return J;
647  }
648 }
649 
650 ideal tropicalStrategy::computeStdOfInitialIdeal(const ideal inI, const ring r) const
651 {
652  // if valuation trivial, then compute std as usual
653  if (isValuationTrivial())
654  return gfanlib_kStd_wrapper(inI,r);
655 
656  // if valuation non-trivial, then uniformizing parameter is in ideal
657  // so switch to residue field first and compute standard basis over the residue field
658  ring rShortcut = copyAndChangeCoefficientRing(r);
659  nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
660  int k = IDELEMS(inI);
661  ideal inIShortcut = idInit(k);
662  for (int i=0; i<k; i++)
663  inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
664  ideal inJShortcut = gfanlib_kStd_wrapper(inIShortcut,rShortcut);
665 
666  // and lift the result back to the ring with valuation
667  nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
668  k = IDELEMS(inJShortcut);
669  ideal inJ = idInit(k+1);
670  inJ->m[0] = p_One(r);
671  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
672  p_SetCoeff(inJ->m[0],identity(uniformizingParameter,startingRing->cf,r->cf),r);
673  for (int i=0; i<k; i++)
674  inJ->m[i+1] = p_PermPoly(inJShortcut->m[i],NULL,rShortcut,r,takingRepresentatives,NULL,0);
675 
676  id_Delete(&inJShortcut,rShortcut);
677  id_Delete(&inIShortcut,rShortcut);
678  rDelete(rShortcut);
679  return inJ;
680 }
681 
682 ideal tropicalStrategy::computeLift(const ideal inJs, const ring s, const ideal inIr, const ideal Ir, const ring r) const
683 {
684  int k = IDELEMS(inJs);
685  ideal inJr = idInit(k);
686  nMapFunc identitysr = n_SetMap(s->cf,r->cf);
687  for (int i=0; i<k; i++)
688  inJr->m[i] = p_PermPoly(inJs->m[i],NULL,s,r,identitysr,NULL,0);
689 
690  ideal Jr = computeWitness(inJr,inIr,Ir,r);
691  nMapFunc identityrs = n_SetMap(r->cf,s->cf);
692  ideal Js = idInit(k);
693  for (int i=0; i<k; i++)
694  Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identityrs,NULL,0);
695  return Js;
696 }
697 
698 ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
699 {
700  // copy shortcutRing and change to desired ordering
701  bool ok;
702  ring s = rCopy0(r);
703  int n = rVar(s);
704  deleteOrdering(s);
705  gfan::ZVector wAdjusted = adjustWeightForHomogeneity(w);
706  gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(v,wAdjusted);
707  s->order = (int*) omAlloc0(5*sizeof(int));
708  s->block0 = (int*) omAlloc0(5*sizeof(int));
709  s->block1 = (int*) omAlloc0(5*sizeof(int));
710  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
711  s->order[0] = ringorder_a;
712  s->block0[0] = 1;
713  s->block1[0] = n;
714  s->wvhdl[0] = ZVectorToIntStar(wAdjusted,ok);
715  s->order[1] = ringorder_a;
716  s->block0[1] = 1;
717  s->block1[1] = n;
718  s->wvhdl[1] = ZVectorToIntStar(vAdjusted,ok);
719  s->order[2] = ringorder_lp;
720  s->block0[2] = 1;
721  s->block1[2] = n;
722  s->order[3] = ringorder_C;
723  rComplete(s);
724  rTest(s);
725 
726  return s;
727 }
728 
729 ring tropicalStrategy::copyAndChangeOrderingLS(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
730 {
731  // copy shortcutRing and change to desired ordering
732  bool ok;
733  ring s = rCopy0(r);
734  int n = rVar(s);
735  deleteOrdering(s);
736  s->order = (int*) omAlloc0(5*sizeof(int));
737  s->block0 = (int*) omAlloc0(5*sizeof(int));
738  s->block1 = (int*) omAlloc0(5*sizeof(int));
739  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
740  s->order[0] = ringorder_a;
741  s->block0[0] = 1;
742  s->block1[0] = n;
743  s->wvhdl[0] = ZVectorToIntStar(w,ok);
744  s->order[1] = ringorder_a;
745  s->block0[1] = 1;
746  s->block1[1] = n;
747  s->wvhdl[1] = ZVectorToIntStar(v,ok);
748  s->order[2] = ringorder_lp;
749  s->block0[2] = 1;
750  s->block1[2] = n;
751  s->order[3] = ringorder_C;
752  rComplete(s);
753  rTest(s);
754 
755  return s;
756 }
757 
758 std::pair<ideal,ring> tropicalStrategy::computeFlip(const ideal Ir, const ring r,
759  const gfan::ZVector &interiorPoint,
760  const gfan::ZVector &facetNormal) const
761 {
762  assume(isValuationTrivial() || interiorPoint[0].sign()<0);
764  assume(checkWeightVector(Ir,r,interiorPoint));
765 
766  // get a generating system of the initial ideal
767  // and compute a standard basis with respect to adjacent ordering
768  ideal inIr = initial(Ir,r,interiorPoint);
769  ring sAdjusted = copyAndChangeOrderingWP(r,interiorPoint,facetNormal);
770  nMapFunc identity = n_SetMap(r->cf,sAdjusted->cf);
771  int k = IDELEMS(Ir);
772  ideal inIsAdjusted = idInit(k);
773  for (int i=0; i<k; i++)
774  inIsAdjusted->m[i] = p_PermPoly(inIr->m[i],NULL,r,sAdjusted,identity,NULL,0);
775  ideal inJsAdjusted = computeStdOfInitialIdeal(inIsAdjusted,sAdjusted);
776 
777  // find witnesses of the new standard basis elements of the initial ideal
778  // with the help of the old standard basis of the ideal
779  k = IDELEMS(inJsAdjusted);
780  ideal inJr = idInit(k);
781  identity = n_SetMap(sAdjusted->cf,r->cf);
782  for (int i=0; i<k; i++)
783  inJr->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,sAdjusted,r,identity,NULL,0);
784 
785  ideal Jr = computeWitness(inJr,inIr,Ir,r);
786  ring s = copyAndChangeOrderingLS(r,interiorPoint,facetNormal);
787  identity = n_SetMap(r->cf,s->cf);
788  ideal Js = idInit(k);
789  for (int i=0; i<k; i++)
790  Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identity,NULL,0);
791 
792  reduce(Js,s);
793  assume(areIdealsEqual(Js,s,Ir,r));
795  assume(checkWeightVector(Js,s,interiorPoint));
796 
797  // cleanup
798  id_Delete(&inIsAdjusted,sAdjusted);
799  id_Delete(&inJsAdjusted,sAdjusted);
800  rDelete(sAdjusted);
801  id_Delete(&inIr,r);
802  id_Delete(&Jr,r);
803  id_Delete(&inJr,r);
804 
806  return std::make_pair(Js,s);
807 }
808 
809 
810 bool tropicalStrategy::checkForUniformizingBinomial(const ideal I, const ring r) const
811 {
812  // if the valuation is trivial,
813  // then there is no special condition the first generator has to fullfill
814  if (isValuationTrivial())
815  return true;
816 
817  // if the valuation is non-trivial then checks if the first generator is p-t
818  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
819  poly p = p_One(r);
820  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
821  poly t = p_One(r);
822  p_SetExp(t,1,1,r);
823  p_Setm(t,r);
824  poly pt = p_Add_q(p,p_Neg(t,r),r);
825 
826  for (int i=0; i<IDELEMS(I); i++)
827  {
828  if (p_EqualPolys(I->m[i],pt,r))
829  {
830  p_Delete(&pt,r);
831  return true;
832  }
833  }
834  p_Delete(&pt,r);
835  return false;
836 }
837 
838 int tropicalStrategy::findPositionOfUniformizingBinomial(const ideal I, const ring r) const
839 {
841 
842  // if the valuation is non-trivial then checks if the first generator is p-t
843  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
844  poly p = p_One(r);
845  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
846  poly t = p_One(r);
847  p_SetExp(t,1,1,r);
848  p_Setm(t,r);
849  poly pt = p_Add_q(p,p_Neg(t,r),r);
850 
851  for (int i=0; i<IDELEMS(I); i++)
852  {
853  if (p_EqualPolys(I->m[i],pt,r))
854  {
855  p_Delete(&pt,r);
856  return i;
857  }
858  }
859  p_Delete(&pt,r);
860  return -1;
861 }
862 
863 bool tropicalStrategy::checkForUniformizingParameter(const ideal inI, const ring r) const
864 {
865  // if the valuation is trivial,
866  // then there is no special condition the first generator has to fullfill
867  if (isValuationTrivial())
868  return true;
869 
870  // if the valuation is non-trivial then checks if the first generator is p
871  if (inI->m[0]==NULL)
872  return false;
873  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
874  poly p = p_One(r);
875  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
876 
877  for (int i=0; i<IDELEMS(inI); i++)
878  {
879  if (p_EqualPolys(inI->m[i],p,r))
880  {
881  p_Delete(&p,r);
882  return true;
883  }
884  }
885  p_Delete(&p,r);
886  return false;
887 }
#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...
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:2971
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:519
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:95
ring getOriginalRing() const
returns the polynomial ring over the field with valuation
bool isOrderingLocalInT(const ring r)
return P p
Definition: myNF.cc:203
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:472
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:80
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:542
#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:580
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:1474
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:93
#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:407
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:804
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
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:10
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:32
poly * m
Definition: matpol.h:19
static int rBlocks(ring r)
Definition: ring.h:556
const ring r
Definition: syzextra.cc:208
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:3938
BOOLEAN linealitySpace(leftv res, leftv args)
Definition: bbcone.cc:881
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:551
poly p_One(const ring r)
Definition: p_polys.cc:1312
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:3351
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:464
int j
Definition: myNF.cc:70
void deleteOrdering(ring r)
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:403
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:1321
ideal getStartingIdeal() const
returns the input ideal
int scDimInt(ideal S, ideal Q)
Definition: hdegree.cc:72
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
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:742
#define rTest(r)
Definition: ring.h:775
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:43
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:123
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:24
void mp_Delete(matrix *a, const ring r)
Definition: matpol.cc:792
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:725
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:4320
gfan::ZVector valued_adjustWeightUnderHomogeneity(const gfan::ZVector &e, const gfan::ZVector &w)
void rChangeCurrRing(ring r)
Definition: polys.cc:12
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:47
static FORCE_INLINE coeffs nCopyCoeff(const coeffs r)
"copy" coeffs, i.e. increment ref
Definition: coeffs.h:433
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:843
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
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:1612
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:483
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:477
#define NULL
Definition: omList.c:10
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:455
void rDelete(ring r)
unconditionally deletes fields in r
Definition: ring.cc:448
static BOOLEAN rField_is_Ring_Z(const ring r)
Definition: ring.h:474
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:169
#define pNext(p)
Definition: monomials.h:43
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:228
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:57
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:1013
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete &#39;p&#39;
Definition: coeffs.h:459
ring originalRing
polynomial ring over a field with valuation
polyrec * poly
Definition: hilb.h:10
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:877
static Poly * h
Definition: janet.cc:978
const poly b
Definition: syzextra.cc:213
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:496
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1020
static int sign(int x)
Definition: ring.cc:3328
#define omAlloc0(size)
Definition: omAllocDecl.h:211
int l
Definition: cfEzgcd.cc:94
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)
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:334
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