My Project
Loading...
Searching...
No Matches
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 = 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
432void 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
448ring 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
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
496std::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);
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
574ideal 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));
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
656ideal 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
688ideal 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
704ring 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
734ring 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
762std::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
814bool 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
842int 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
867bool 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:544
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:448
#define n_Test(a, r)
BOOLEAN n_Test(number a, const coeffs r)
Definition: coeffs.h:709
@ 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:512
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:697
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:413
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:452
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:535
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:568
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition: coeffs.h:465
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:78
#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
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:3182
void mp_Delete(matrix *a, const ring r)
Definition: matpol.cc:873
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:389
#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:4130
poly p_Div_nn(poly p, const number n, const ring r)
Definition: p_polys.cc:1505
poly p_One(const ring r)
Definition: p_polys.cc:1313
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition: p_polys.cc:4512
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1105
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:934
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1112
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:486
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:231
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:410
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:467
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:899
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:844
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:3450
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1421
void rDelete(ring r)
unconditionally deletes fields in r
Definition: ring.cc:450
static int sign(int x)
Definition: ring.cc:3427
ring rCopy(ring r)
Definition: ring.cc:1731
static BOOLEAN rField_is_Z(const ring r)
Definition: ring.h:509
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:500
static int rBlocks(const ring r)
Definition: ring.h:568
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:506
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:592
#define rTest(r)
Definition: ring.h:782
#define rField_is_Ring(R)
Definition: ring.h:485
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:87
#define Q
Definition: sirandom.c:26
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