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 = NULL;
424  p = identity(uniformizingParameter,startingRing->cf,r->cf);
425  bool b = extraReductionAlgorithm(I,r,p);
426  // putUniformizingBinomialInFront(I,r,p);
427  if (p!=NULL) n_Delete(&p,r->cf);
428 
429  return b;
430 }
431 
432 void tropicalStrategy::pReduce(ideal I, const ring r) const
433 {
434  rTest(r);
435  id_Test(I,r);
436 
437  if (isValuationTrivial())
438  return;
439 
440  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
441  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
442  ::pReduce(I,p,r);
443  n_Delete(&p,r->cf);
444 
445  return;
446 }
447 
448 ring tropicalStrategy::getShortcutRingPrependingWeight(const ring r, const gfan::ZVector &v) const
449 {
450  ring rShortcut = rCopy0(r,FALSE); // do not copy q-ideal
451 
452  // save old ordering
453  rRingOrder_t* order = rShortcut->order;
454  int* block0 = rShortcut->block0;
455  int* block1 = rShortcut->block1;
456  int** wvhdl = rShortcut->wvhdl;
457 
458  // adjust weight and create new ordering
459  gfan::ZVector w = adjustWeightForHomogeneity(v);
460  int h = rBlocks(r); int n = rVar(r);
461  rShortcut->order = (rRingOrder_t*) omAlloc0((h+2)*sizeof(rRingOrder_t));
462  rShortcut->block0 = (int*) omAlloc0((h+2)*sizeof(int));
463  rShortcut->block1 = (int*) omAlloc0((h+2)*sizeof(int));
464  rShortcut->wvhdl = (int**) omAlloc0((h+2)*sizeof(int*));
465  rShortcut->order[0] = ringorder_a;
466  rShortcut->block0[0] = 1;
467  rShortcut->block1[0] = n;
468  bool overflow;
469  rShortcut->wvhdl[0] = ZVectorToIntStar(w,overflow);
470  for (int i=1; i<=h; i++)
471  {
472  rShortcut->order[i] = order[i-1];
473  rShortcut->block0[i] = block0[i-1];
474  rShortcut->block1[i] = block1[i-1];
475  rShortcut->wvhdl[i] = wvhdl[i-1];
476  }
477 
478  // if valuation non-trivial, change coefficient ring to residue field
479  if (isValuationNonTrivial())
480  {
481  nKillChar(rShortcut->cf);
482  rShortcut->cf = nCopyCoeff(shortcutRing->cf);
483  }
484  rComplete(rShortcut);
485  rTest(rShortcut);
486 
487  // delete old ordering
488  omFree(order);
489  omFree(block0);
490  omFree(block1);
491  omFree(wvhdl);
492 
493  return rShortcut;
494 }
495 
496 std::pair<poly,int> tropicalStrategy::checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector &w) const
497 {
498  // quick check whether I already contains an ideal
499  int k = IDELEMS(I);
500  for (int i=0; i<k; i++)
501  {
502  poly g = I->m[i];
503  if (g!=NULL
504  && pNext(g)==NULL
505  && (isValuationTrivial() || n_IsOne(p_GetCoeff(g,r),r->cf)))
506  return std::pair<poly,int>(g,i);
507  }
508 
509  ring rShortcut;
510  ideal inIShortcut;
511  if (w.size()>0)
512  {
513  // if needed, prepend extra weight for homogeneity
514  // switch to residue field if valuation is non trivial
515  rShortcut = getShortcutRingPrependingWeight(r,w);
516 
517  // compute the initial ideal and map it into the constructed ring
518  // if switched to residue field, remove possibly 0 elements
519  ideal inI = initial(I,r,w);
520  inIShortcut = idInit(k);
521  nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
522  for (int i=0; i<k; i++)
523  inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
524  if (isValuationNonTrivial())
525  idSkipZeroes(inIShortcut);
526  id_Delete(&inI,r);
527  }
528  else
529  {
530  rShortcut = r;
531  inIShortcut = I;
532  }
533 
534  gfan::ZCone C0 = homogeneitySpace(inIShortcut,rShortcut);
535  gfan::ZCone pos = gfan::ZCone::positiveOrthant(C0.ambientDimension());
536  gfan::ZCone C0pos = intersection(C0,pos);
537  C0pos.canonicalize();
538  gfan::ZVector wpos = C0pos.getRelativeInteriorPoint();
540 
541  // check initial ideal for monomial and
542  // if it exsists, return a copy of the monomial in the input ring
543  poly p = searchForMonomialViaStepwiseSaturation(inIShortcut,rShortcut,wpos);
544  poly monomial = NULL;
545  if (p!=NULL)
546  {
547  monomial=p_One(r);
548  for (int i=1; i<=rVar(r); i++)
549  p_SetExp(monomial,i,p_GetExp(p,i,rShortcut),r);
550  p_Setm(monomial,r);
551  p_Delete(&p,rShortcut);
552  }
553 
554 
555  if (w.size()>0)
556  {
557  // if needed, cleanup
558  id_Delete(&inIShortcut,rShortcut);
559  rDelete(rShortcut);
560  }
561  return std::pair<poly,int>(monomial,-1);
562 }
563 
565 {
566  ring rShortcut = rCopy0(r,FALSE); // do not copy q-ideal
567  nKillChar(rShortcut->cf);
568  rShortcut->cf = nCopyCoeff(shortcutRing->cf);
569  rComplete(rShortcut);
570  rTest(rShortcut);
571  return rShortcut;
572 }
573 
574 ideal tropicalStrategy::computeWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
575 {
576  // if the valuation is trivial and the ring and ideal have not been extended,
577  // then it is sufficient to return the difference between the elements of inJ
578  // and their normal forms with respect to I and r
579  if (isValuationTrivial())
580  return witness(inJ,I,r);
581  // if the valuation is non-trivial and the ring and ideal have been extended,
582  // then we can make a shortcut through the residue field
583  else
584  {
585  assume(IDELEMS(inI)==IDELEMS(I));
586  int uni = findPositionOfUniformizingBinomial(I,r);
587  assume(uni>=0);
588  /**
589  * change ground ring into finite field
590  * and map the data into it
591  */
592  ring rShortcut = copyAndChangeCoefficientRing(r);
593 
594  int k = IDELEMS(inJ);
595  int l = IDELEMS(I);
596  ideal inJShortcut = idInit(k);
597  ideal inIShortcut = idInit(l);
598  nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
599  for (int i=0; i<k; i++)
600  inJShortcut->m[i] = p_PermPoly(inJ->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
601  for (int j=0; j<l; j++)
602  inIShortcut->m[j] = p_PermPoly(inI->m[j],NULL,r,rShortcut,takingResidues,NULL,0);
603  id_Test(inJShortcut,rShortcut);
604  id_Test(inIShortcut,rShortcut);
605 
606  /**
607  * Compute a division with remainder over the finite field
608  * and map the result back to r
609  */
610  matrix QShortcut = divisionDiscardingRemainder(inJShortcut,inIShortcut,rShortcut);
611  matrix Q = mpNew(l,k);
612  nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
613  for (int ij=k*l-1; ij>=0; ij--)
614  Q->m[ij] = p_PermPoly(QShortcut->m[ij],NULL,rShortcut,r,takingRepresentatives,NULL,0);
615 
616  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
617  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
618 
619  /**
620  * Compute the normal forms
621  */
622  ideal J = idInit(k);
623  for (int j=0; j<k; j++)
624  {
625  poly q0 = p_Copy(inJ->m[j],r);
626  for (int i=0; i<l; i++)
627  {
628  poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
629  poly inIi = p_Copy(inI->m[i],r);
630  q0 = p_Add_q(q0,p_Neg(p_Mult_q(qij,inIi,r),r),r);
631  }
632  q0 = p_Div_nn(q0,p,r);
633  poly q0g0 = p_Mult_q(q0,p_Copy(I->m[uni],r),r);
634  // q0 = NULL;
635  poly qigi = NULL;
636  for (int i=0; i<l; i++)
637  {
638  poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
639  // poly inIi = p_Copy(I->m[i],r);
640  poly Ii = p_Copy(I->m[i],r);
641  qigi = p_Add_q(qigi,p_Mult_q(qij,Ii,r),r);
642  }
643  J->m[j] = p_Add_q(q0g0,qigi,r);
644  }
645 
646  id_Delete(&inIShortcut,rShortcut);
647  id_Delete(&inJShortcut,rShortcut);
648  mp_Delete(&QShortcut,rShortcut);
649  rDelete(rShortcut);
650  mp_Delete(&Q,r);
651  n_Delete(&p,r->cf);
652  return J;
653  }
654 }
655 
656 ideal tropicalStrategy::computeStdOfInitialIdeal(const ideal inI, const ring r) const
657 {
658  // if valuation trivial, then compute std as usual
659  if (isValuationTrivial())
660  return gfanlib_kStd_wrapper(inI,r);
661 
662  // if valuation non-trivial, then uniformizing parameter is in ideal
663  // so switch to residue field first and compute standard basis over the residue field
664  ring rShortcut = copyAndChangeCoefficientRing(r);
665  nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
666  int k = IDELEMS(inI);
667  ideal inIShortcut = idInit(k);
668  for (int i=0; i<k; i++)
669  inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
670  ideal inJShortcut = gfanlib_kStd_wrapper(inIShortcut,rShortcut);
671 
672  // and lift the result back to the ring with valuation
673  nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
674  k = IDELEMS(inJShortcut);
675  ideal inJ = idInit(k+1);
676  inJ->m[0] = p_One(r);
677  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
678  p_SetCoeff(inJ->m[0],identity(uniformizingParameter,startingRing->cf,r->cf),r);
679  for (int i=0; i<k; i++)
680  inJ->m[i+1] = p_PermPoly(inJShortcut->m[i],NULL,rShortcut,r,takingRepresentatives,NULL,0);
681 
682  id_Delete(&inJShortcut,rShortcut);
683  id_Delete(&inIShortcut,rShortcut);
684  rDelete(rShortcut);
685  return inJ;
686 }
687 
688 ideal tropicalStrategy::computeLift(const ideal inJs, const ring s, const ideal inIr, const ideal Ir, const ring r) const
689 {
690  int k = IDELEMS(inJs);
691  ideal inJr = idInit(k);
692  nMapFunc identitysr = n_SetMap(s->cf,r->cf);
693  for (int i=0; i<k; i++)
694  inJr->m[i] = p_PermPoly(inJs->m[i],NULL,s,r,identitysr,NULL,0);
695 
696  ideal Jr = computeWitness(inJr,inIr,Ir,r);
697  nMapFunc identityrs = n_SetMap(r->cf,s->cf);
698  ideal Js = idInit(k);
699  for (int i=0; i<k; i++)
700  Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identityrs,NULL,0);
701  return Js;
702 }
703 
704 ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
705 {
706  // copy shortcutRing and change to desired ordering
707  bool ok;
708  ring s = rCopy0(r,FALSE,FALSE);
709  int n = rVar(s);
710  gfan::ZVector wAdjusted = adjustWeightForHomogeneity(w);
711  gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(v,wAdjusted);
712  s->order = (rRingOrder_t*) omAlloc0(5*sizeof(rRingOrder_t));
713  s->block0 = (int*) omAlloc0(5*sizeof(int));
714  s->block1 = (int*) omAlloc0(5*sizeof(int));
715  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
716  s->order[0] = ringorder_a;
717  s->block0[0] = 1;
718  s->block1[0] = n;
719  s->wvhdl[0] = ZVectorToIntStar(wAdjusted,ok);
720  s->order[1] = ringorder_a;
721  s->block0[1] = 1;
722  s->block1[1] = n;
723  s->wvhdl[1] = ZVectorToIntStar(vAdjusted,ok);
724  s->order[2] = ringorder_lp;
725  s->block0[2] = 1;
726  s->block1[2] = n;
727  s->order[3] = ringorder_C;
728  rComplete(s);
729  rTest(s);
730 
731  return s;
732 }
733 
734 ring tropicalStrategy::copyAndChangeOrderingLS(const ring r, const gfan::ZVector &w, const gfan::ZVector &v) const
735 {
736  // copy shortcutRing and change to desired ordering
737  bool ok;
738  ring s = rCopy0(r,FALSE,FALSE);
739  int n = rVar(s);
740  s->order = (rRingOrder_t*) omAlloc0(5*sizeof(rRingOrder_t));
741  s->block0 = (int*) omAlloc0(5*sizeof(int));
742  s->block1 = (int*) omAlloc0(5*sizeof(int));
743  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
744  s->order[0] = ringorder_a;
745  s->block0[0] = 1;
746  s->block1[0] = n;
747  s->wvhdl[0] = ZVectorToIntStar(w,ok);
748  s->order[1] = ringorder_a;
749  s->block0[1] = 1;
750  s->block1[1] = n;
751  s->wvhdl[1] = ZVectorToIntStar(v,ok);
752  s->order[2] = ringorder_lp;
753  s->block0[2] = 1;
754  s->block1[2] = n;
755  s->order[3] = ringorder_C;
756  rComplete(s);
757  rTest(s);
758 
759  return s;
760 }
761 
762 std::pair<ideal,ring> tropicalStrategy::computeFlip(const ideal Ir, const ring r,
763  const gfan::ZVector &interiorPoint,
764  const gfan::ZVector &facetNormal) const
765 {
766  assume(isValuationTrivial() || interiorPoint[0].sign()<0);
768  assume(checkWeightVector(Ir,r,interiorPoint));
769 
770  // get a generating system of the initial ideal
771  // and compute a standard basis with respect to adjacent ordering
772  ideal inIr = initial(Ir,r,interiorPoint);
773  ring sAdjusted = copyAndChangeOrderingWP(r,interiorPoint,facetNormal);
774  nMapFunc identity = n_SetMap(r->cf,sAdjusted->cf);
775  int k = IDELEMS(Ir);
776  ideal inIsAdjusted = idInit(k);
777  for (int i=0; i<k; i++)
778  inIsAdjusted->m[i] = p_PermPoly(inIr->m[i],NULL,r,sAdjusted,identity,NULL,0);
779  ideal inJsAdjusted = computeStdOfInitialIdeal(inIsAdjusted,sAdjusted);
780 
781  // find witnesses of the new standard basis elements of the initial ideal
782  // with the help of the old standard basis of the ideal
783  k = IDELEMS(inJsAdjusted);
784  ideal inJr = idInit(k);
785  identity = n_SetMap(sAdjusted->cf,r->cf);
786  for (int i=0; i<k; i++)
787  inJr->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,sAdjusted,r,identity,NULL,0);
788 
789  ideal Jr = computeWitness(inJr,inIr,Ir,r);
790  ring s = copyAndChangeOrderingLS(r,interiorPoint,facetNormal);
791  identity = n_SetMap(r->cf,s->cf);
792  ideal Js = idInit(k);
793  for (int i=0; i<k; i++)
794  Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identity,NULL,0);
795 
796  reduce(Js,s);
797  assume(areIdealsEqual(Js,s,Ir,r));
799  assume(checkWeightVector(Js,s,interiorPoint));
800 
801  // cleanup
802  id_Delete(&inIsAdjusted,sAdjusted);
803  id_Delete(&inJsAdjusted,sAdjusted);
804  rDelete(sAdjusted);
805  id_Delete(&inIr,r);
806  id_Delete(&Jr,r);
807  id_Delete(&inJr,r);
808 
810  return std::make_pair(Js,s);
811 }
812 
813 
814 bool tropicalStrategy::checkForUniformizingBinomial(const ideal I, const ring r) const
815 {
816  // if the valuation is trivial,
817  // then there is no special condition the first generator has to fullfill
818  if (isValuationTrivial())
819  return true;
820 
821  // if the valuation is non-trivial then checks if the first generator is p-t
822  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
823  poly p = p_One(r);
824  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
825  poly t = p_One(r);
826  p_SetExp(t,1,1,r);
827  p_Setm(t,r);
828  poly pt = p_Add_q(p,p_Neg(t,r),r);
829 
830  for (int i=0; i<IDELEMS(I); i++)
831  {
832  if (p_EqualPolys(I->m[i],pt,r))
833  {
834  p_Delete(&pt,r);
835  return true;
836  }
837  }
838  p_Delete(&pt,r);
839  return false;
840 }
841 
842 int tropicalStrategy::findPositionOfUniformizingBinomial(const ideal I, const ring r) const
843 {
845 
846  // if the valuation is non-trivial then checks if the first generator is p-t
847  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
848  poly p = p_One(r);
849  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
850  poly t = p_One(r);
851  p_SetExp(t,1,1,r);
852  p_Setm(t,r);
853  poly pt = p_Add_q(p,p_Neg(t,r),r);
854 
855  for (int i=0; i<IDELEMS(I); i++)
856  {
857  if (p_EqualPolys(I->m[i],pt,r))
858  {
859  p_Delete(&pt,r);
860  return i;
861  }
862  }
863  p_Delete(&pt,r);
864  return -1;
865 }
866 
867 bool tropicalStrategy::checkForUniformizingParameter(const ideal inI, const ring r) const
868 {
869  // if the valuation is trivial,
870  // then there is no special condition the first generator has to fullfill
871  if (isValuationTrivial())
872  return true;
873 
874  // if the valuation is non-trivial then checks if the first generator is p
875  if (inI->m[0]==NULL)
876  return false;
877  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
878  poly p = p_One(r);
879  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
880 
881  for (int i=0; i<IDELEMS(inI); i++)
882  {
883  if (p_EqualPolys(inI->m[i],p,r))
884  {
885  p_Delete(&p,r);
886  return true;
887  }
888  }
889  p_Delete(&p,r);
890  return false;
891 }
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