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 && (n_IsUnit(p_GetCoeff(I->m[i],currRing->cf),currRing->cf))
34 )
35 return -1;
36 ideal vv = id_Head(I,currRing);
37 if (i != -1) pDelete(&vv->m[i]);
38 d = scDimInt(vv, currRing->qideal);
39 if (rField_is_Z(currRing) && (i==-1)) d++;
40 idDelete(&vv);
41 return d;
42 }
43 else
44 d = scDimInt(I,currRing->qideal);
45 if (origin != r)
46 rChangeCurrRing(origin);
47 return d;
48}
49
50static void swapElements(ideal I, ideal J)
51{
52 assume(IDELEMS(I)==IDELEMS(J));
53
54 for (int i=IDELEMS(I)-1; i>=0; i--)
55 {
56 poly cache = I->m[i];
57 I->m[i] = J->m[i];
58 J->m[i] = cache;
59 }
60}
61
62static bool noExtraReduction(ideal I, ring r, number /*p*/)
63{
64 int n = rVar(r);
65 gfan::ZVector allOnes(n);
66 for (int i=0; i<n; i++)
67 allOnes[i] = 1;
68 ring rShortcut = rCopy0(r);
69
70 rRingOrder_t* order = rShortcut->order;
71 int* block0 = rShortcut->block0;
72 int* block1 = rShortcut->block1;
73 int** wvhdl = rShortcut->wvhdl;
74
75 int h = rBlocks(r);
76 rShortcut->order = (rRingOrder_t*) omAlloc0((h+2)*sizeof(rRingOrder_t));
77 rShortcut->block0 = (int*) omAlloc0((h+2)*sizeof(int));
78 rShortcut->block1 = (int*) omAlloc0((h+2)*sizeof(int));
79 rShortcut->wvhdl = (int**) omAlloc0((h+2)*sizeof(int*));
80 rShortcut->order[0] = ringorder_a;
81 rShortcut->block0[0] = 1;
82 rShortcut->block1[0] = n;
83 bool overflow;
84 rShortcut->wvhdl[0] = ZVectorToIntStar(allOnes,overflow);
85 for (int i=1; i<=h; i++)
86 {
87 rShortcut->order[i] = order[i-1];
88 rShortcut->block0[i] = block0[i-1];
89 rShortcut->block1[i] = block1[i-1];
90 rShortcut->wvhdl[i] = wvhdl[i-1];
91 }
92 //rShortcut->order[h+1] = (rRingOrder_t)0; -- done by omAlloc0
93 //rShortcut->block0[h+1] = 0;
94 //rShortcut->block1[h+1] = 0;
95 //rShortcut->wvhdl[h+1] = NULL;
96
97 rComplete(rShortcut);
98 rTest(rShortcut);
99
100 omFree(order);
101 omFree(block0);
102 omFree(block1);
103 omFree(wvhdl);
104
105 int k = IDELEMS(I);
106 ideal IShortcut = idInit(k);
107 nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
108 for (int i=0; i<k; i++)
109 {
110 if(I->m[i]!=NULL)
111 {
112 IShortcut->m[i] = p_PermPoly(I->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
113 }
114 }
115
116 ideal JShortcut = gfanlib_kStd_wrapper(IShortcut,rShortcut);
117
118 ideal J = idInit(k);
119 nMapFunc outofShortcut = n_SetMap(rShortcut->cf,r->cf);
120 for (int i=0; i<k; i++)
121 J->m[i] = p_PermPoly(JShortcut->m[i],NULL,rShortcut,r,outofShortcut,NULL,0);
122
123 assume(areIdealsEqual(J,r,I,r));
124 swapElements(I,J);
125 id_Delete(&IShortcut,rShortcut);
126 id_Delete(&JShortcut,rShortcut);
127 rDelete(rShortcut);
128 id_Delete(&J,r);
129 return false;
130}
131
132/**
133 * Initializes all relevant structures and information for the trivial valuation case,
134 * i.e. computing a tropical variety without any valuation.
135 */
161
162/**
163 * Given a polynomial ring r over the rational numbers and a weighted ordering,
164 * returns a polynomial ring s over the integers with one extra variable, which is weighted -1.
165 */
166static ring constructStartingRing(ring r)
167{
169
170 ring s = rCopy0(r,FALSE,FALSE);
171 nKillChar(s->cf);
172 s->cf = nInitChar(n_Z,NULL);
173
174 int n = rVar(s)+1;
175 s->N = n;
176 char** oldNames = s->names;
177 s->names = (char**) omAlloc((n+1)*sizeof(char**));
178 s->names[0] = omStrDup("t");
179 for (int i=1; i<n; i++)
180 s->names[i] = oldNames[i-1];
181 omFree(oldNames);
182
183 s->order = (rRingOrder_t*) omAlloc0(3*sizeof(rRingOrder_t));
184 s->block0 = (int*) omAlloc0(3*sizeof(int));
185 s->block1 = (int*) omAlloc0(3*sizeof(int));
186 s->wvhdl = (int**) omAlloc0(3*sizeof(int**));
187 s->order[0] = ringorder_ws;
188 s->block0[0] = 1;
189 s->block1[0] = n;
190 s->wvhdl[0] = (int*) omAlloc(n*sizeof(int));
191 s->wvhdl[0][0] = 1;
192 if (r->order[0] == ringorder_lp)
193 {
194 s->wvhdl[0][1] = 1;
195 }
196 else if (r->order[0] == ringorder_ls)
197 {
198 s->wvhdl[0][1] = -1;
199 }
200 else if (r->order[0] == ringorder_dp)
201 {
202 for (int i=1; i<n; i++)
203 s->wvhdl[0][i] = -1;
204 }
205 else if (r->order[0] == ringorder_ds)
206 {
207 for (int i=1; i<n; i++)
208 s->wvhdl[0][i] = 1;
209 }
210 else if (r->order[0] == ringorder_ws)
211 {
212 for (int i=1; i<n; i++)
213 s->wvhdl[0][i] = r->wvhdl[0][i-1];
214 }
215 else
216 {
217 for (int i=1; i<n; i++)
218 s->wvhdl[0][i] = -r->wvhdl[0][i-1];
219 }
220 s->order[1] = ringorder_C;
221
222 rComplete(s);
223 rTest(s);
224 return s;
225}
226
227static ideal constructStartingIdeal(ideal originalIdeal, ring originalRing, number uniformizingParameter, ring startingRing)
228{
229 // construct p-t
230 poly g = p_One(startingRing);
231 p_SetCoeff(g,uniformizingParameter,startingRing);
232 pNext(g) = p_One(startingRing);
233 p_SetExp(pNext(g),1,1,startingRing);
234 p_SetCoeff(pNext(g),n_Init(-1,startingRing->cf),startingRing);
235 p_Setm(pNext(g),startingRing);
236 ideal pt = idInit(1);
237 pt->m[0] = g;
238
239 // map originalIdeal from originalRing into startingRing
240 int k = IDELEMS(originalIdeal);
241 ideal J = idInit(k+1);
242 nMapFunc nMap = n_SetMap(originalRing->cf,startingRing->cf);
243 int n = rVar(originalRing);
244 int* shiftByOne = (int*) omAlloc((n+1)*sizeof(int));
245 for (int i=1; i<=n; i++)
246 shiftByOne[i]=i+1;
247 for (int i=0; i<k; i++)
248 {
249 if(originalIdeal->m[i]!=NULL)
250 {
251 J->m[i] = p_PermPoly(originalIdeal->m[i],shiftByOne,originalRing,startingRing,nMap,NULL,0);
252 }
253 }
254 omFreeSize(shiftByOne,(n+1)*sizeof(int));
255
256 ring origin = currRing;
257 rChangeCurrRing(startingRing);
258 ideal startingIdeal = kNF(pt,startingRing->qideal,J); // mathematically redundant,
259 rChangeCurrRing(origin); // but helps with upcoming std computation
260 // ideal startingIdeal = J; J = NULL;
261 assume(startingIdeal->m[k]==NULL);
262 startingIdeal->m[k] = pt->m[0];
263 startingIdeal = gfanlib_kStd_wrapper(startingIdeal,startingRing);
264
265 id_Delete(&J,startingRing);
266 pt->m[0] = NULL;
267 id_Delete(&pt,startingRing);
268 return startingIdeal;
269}
270
271/***
272 * Initializes all relevant structures and information for the valued case,
273 * i.e. computing a tropical variety over the rational numbers with p-adic valuation
274 **/
275tropicalStrategy::tropicalStrategy(ideal J, number q, ring s):
279 linealitySpace(gfan::ZCone()), // to come, see below
280 startingRing(NULL), // to come, see below
281 startingIdeal(NULL), // to come, see below
282 uniformizingParameter(NULL), // to come, see below
283 shortcutRing(NULL), // to come, see below
284 onlyLowerHalfSpace(true),
288{
289 /* assume that the ground field of the originalRing is Q */
291
292 /* replace Q with Z for the startingRing
293 * and add an extra variable for tracking the uniformizing parameter */
295
296 /* map the uniformizing parameter into the new coefficient domain */
299
300 /* map the input ideal into the new polynomial ring */
303
305
306 /* construct the shorcut ring */
307 shortcutRing = rCopy0(startingRing,FALSE); // do not copy q-ideal
312}
313
315 originalRing(rCopy(currentStrategy.getOriginalRing())),
316 originalIdeal(id_Copy(currentStrategy.getOriginalIdeal(),currentStrategy.getOriginalRing())),
317 expectedDimension(currentStrategy.getExpectedDimension()),
318 linealitySpace(currentStrategy.getHomogeneitySpace()),
319 startingRing(rCopy(currentStrategy.getStartingRing())),
320 startingIdeal(id_Copy(currentStrategy.getStartingIdeal(),currentStrategy.getStartingRing())),
327{
332 if (currentStrategy.getUniformizingParameter())
333 {
334 uniformizingParameter = n_Copy(currentStrategy.getUniformizingParameter(),startingRing->cf);
335 n_Test(uniformizingParameter,startingRing->cf);
336 }
337 if (currentStrategy.getShortcutRing())
338 {
339 shortcutRing = rCopy(currentStrategy.getShortcutRing());
340 rTest(shortcutRing);
341 }
342}
343
360
384
385void tropicalStrategy::putUniformizingBinomialInFront(ideal I, const ring r, const number q) const
386{
387 poly p = p_One(r);
388 p_SetCoeff(p,q,r);
389 poly t = p_One(r);
390 p_SetExp(t,1,1,r);
391 p_Setm(t,r);
392 poly pt = p_Add_q(p,p_Neg(t,r),r);
393
394 int k = IDELEMS(I);
395 int l;
396 for (l=0; l<k; l++)
397 {
398 if (p_EqualPolys(I->m[l],pt,r))
399 break;
400 }
401 p_Delete(&pt,r);
402
403 if (l>1)
404 {
405 pt = I->m[l];
406 for (int i=l; i>0; i--)
407 I->m[l] = I->m[l-1];
408 I->m[0] = pt;
409 pt = NULL;
410 }
411 return;
412}
413
414bool tropicalStrategy::reduce(ideal I, const ring r) const
415{
416 rTest(r);
417 id_Test(I,r);
418
419 nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
420 number p = NULL;
422 p = identity(uniformizingParameter,startingRing->cf,r->cf);
423 bool b = extraReductionAlgorithm(I,r,p);
424 // putUniformizingBinomialInFront(I,r,p);
425 if (p!=NULL) 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 FALSE
Definition auxiliary.h:97
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:4086
return false
Definition cfModGcd.cc:85
g
Definition cfModGcd.cc:4098
CanonicalForm b
Definition cfModGcd.cc:4111
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 ...
gfan::ZCone getHomogeneitySpace() const
returns the homogeneity space of the preimage ideal
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:455
#define n_Test(a, r)
BOOLEAN n_Test(number a, const coeffs r)
Definition coeffs.h:713
@ 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:519
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:701
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition numbers.cc:406
static FORCE_INLINE coeffs nCopyCoeff(const coeffs r)
"copy" coeffs, i.e. increment ref
Definition coeffs.h:437
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition coeffs.h:459
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:80
void nKillChar(coeffs r)
undo all initialisations
Definition numbers.cc:556
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition coeffs.h:472
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
if(!FE_OPT_NO_SHELL_FLAG)
Definition fehelp.cc:1000
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:3224
void mp_Delete(matrix *a, const ring r)
Definition matpol.cc:874
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
ip_smatrix * matrix
Definition matpol.h:43
#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)
#define omFreeSize(addr, size)
#define omAlloc(size)
#define omFree(addr)
#define omAlloc0(size)
#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:4211
poly p_Div_nn(poly p, const number n, const ring r)
Definition p_polys.cc:1506
poly p_One(const ring r)
Definition p_polys.cc:1314
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition p_polys.cc:4621
static poly p_Neg(poly p, const ring r)
Definition p_polys.h:1109
static poly p_Add_q(poly p, poly q, const ring r)
Definition p_polys.h:938
static poly p_Mult_q(poly p, poly q, const ring r)
Definition p_polys.h:1120
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:490
static void p_Setm(poly p, const ring r)
Definition p_polys.h:235
static number p_SetCoeff(poly p, number n, ring r)
Definition p_polys.h:414
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:471
static void p_Delete(poly *p, const ring r)
Definition p_polys.h:903
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition p_polys.h:848
void rChangeCurrRing(ring r)
Definition polys.cc:16
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:187
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:3526
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition ring.cc:1426
void rDelete(ring r)
unconditionally deletes fields in r
Definition ring.cc:454
static int sign(int x)
Definition ring.cc:3503
ring rCopy(ring r)
Definition ring.cc:1736
static BOOLEAN rField_is_Z(const ring r)
Definition ring.h:515
static BOOLEAN rField_is_Zp(const ring r)
Definition ring.h:506
static int rBlocks(const ring r)
Definition ring.h:574
rRingOrder_t
order stuff
Definition ring.h:69
@ ringorder_lp
Definition ring.h:78
@ ringorder_a
Definition ring.h:71
@ ringorder_C
Definition ring.h:74
@ ringorder_ds
Definition ring.h:86
@ ringorder_dp
Definition ring.h:79
@ ringorder_ws
Definition ring.h:88
@ ringorder_ls
degree, ip
Definition ring.h:85
static BOOLEAN rField_is_Q(const ring r)
Definition ring.h:512
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition ring.h:598
#define rTest(r)
Definition ring.h:794
#define rField_is_Ring(R)
Definition ring.h:491
ideal idInit(int idsize, int rank)
initialise an ideal / module
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)
#define id_Test(A, lR)
#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