My Project
tropicalStrategy.cc
Go to the documentation of this file.
1#include "tropicalStrategy.h"
2#include "singularWishlist.h"
3#include "adjustWeights.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"
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 **/
23int dim(ideal I, ring r)
24{
25 ring origin = currRing;
26 if (origin != 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
52static 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
64static 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 */
138tropicalStrategy::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{
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 */
168static ring constructStartingRing(ring r)
169{
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
229static 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 **/
277tropicalStrategy::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 */
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 */
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
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());
376
383
384 return *this;
385}
386
387void 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
416bool 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
430void 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
446ring 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
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
494std::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);
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
572ideal 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));
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
654ideal 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
686ideal 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
702ring 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
732ring 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
760std::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
812bool 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
840int 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
865bool tropicalStrategy::checkForUniformizingParameter(const ideal inI, const ring r) const
866{
867 // if the valuation is trivial,
868 // then there is no special condition the first generator has to fullfill
869 if (isValuationTrivial())
870 return true;
871
872 // if the valuation is non-trivial then checks if the first generator is p
873 if (inI->m[0]==NULL)
874 return false;
875 nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
876 poly p = p_One(r);
877 p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
878
879 for (int i=0; i<IDELEMS(inI); i++)
880 {
881 if (p_EqualPolys(inI->m[i],p,r))
882 {
883 p_Delete(&p,r);
884 return true;
885 }
886 }
887 p_Delete(&p,r);
888 return false;
889}
gfan::ZVector nonvalued_adjustWeightForHomogeneity(const gfan::ZVector &w)
gfan::ZVector nonvalued_adjustWeightUnderHomogeneity(const gfan::ZVector &e, const gfan::ZVector &)
gfan::ZVector valued_adjustWeightForHomogeneity(const gfan::ZVector &w)
gfan::ZVector valued_adjustWeightUnderHomogeneity(const gfan::ZVector &e, const gfan::ZVector &w)
#define NULL
Definition: auxiliary.h:104
#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 sign(const CanonicalForm &a)
int l
Definition: cfEzgcd.cc:100
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
int p
Definition: cfModGcd.cc:4080
return false
Definition: cfModGcd.cc:84
g
Definition: cfModGcd.cc:4092
CanonicalForm b
Definition: cfModGcd.cc:4105
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:548
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:452
#define n_Test(a, r)
BOOLEAN n_Test(number a, const coeffs r)
Definition: coeffs.h:736
@ n_Zp
\F{p < 2^31}
Definition: coeffs.h:30
@ n_Z
only used if HAVE_RINGS is defined
Definition: coeffs.h:44
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:516
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:723
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:358
static FORCE_INLINE coeffs nCopyCoeff(const coeffs r)
"copy" coeffs, i.e. increment ref
Definition: coeffs.h:430
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:456
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:539
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:74
void nKillChar(coeffs r)
undo all initialisations
Definition: numbers.cc:526
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition: coeffs.h:469
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:3169
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
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:4158
poly p_Div_nn(poly p, const number n, const ring r)
Definition: p_polys.cc:1492
poly p_One(const ring r)
Definition: p_polys.cc:1308
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition: p_polys.cc:4540
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1067
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:896
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1074
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:861
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:812
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:3400
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
ring rCopy(ring r)
Definition: ring.cc:1645
static BOOLEAN rField_is_Ring(const ring r)
Definition: ring.h:486
static int rBlocks(ring r)
Definition: ring.h:570
static BOOLEAN rField_is_Z(const ring r)
Definition: ring.h:511
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:502
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:508
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:594
#define rTest(r)
Definition: ring.h:787
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