source: git/Singular/dyn_modules/gfanlib/tropicalStrategy.cc @ e5618c

spielwiese
Last change on this file since e5618c was e5618c, checked in by Yue Ren <ren@…>, 9 years ago
chg: final version for Singular release
  • Property mode set to 100644
File size: 29.7 KB
Line 
1#include <tropicalStrategy.h>
2#include <adjustWeights.h>
3#include <ppinitialReduction.h>
4#include <ttinitialReduction.h>
5#include <tropical.h>
6#include <std_wrapper.h>
7#include <tropicalCurves.h>
8
9// for various commands in dim(ideal I, ring r):
10#include <kernel/ideals.h>
11#include <kernel/combinatorics/stairc.h>
12#include <kernel/GBEngine/kstd1.h>
13#include <Singular/ipshell.h> // for isPrime(int i)
14
15/***
16 * Computes the dimension of an ideal I in ring r
17 * Copied from jjDim in iparith.cc
18 **/
19int dim(ideal I, ring r)
20{
21  ring origin = currRing;
22  if (origin != r)
23    rChangeCurrRing(r);
24  int d;
25  if (rField_is_Ring(currRing))
26  {
27    int i = idPosConstant(I);
28    if ((i != -1) && (n_IsUnit(p_GetCoeff(I->m[i],currRing->cf),currRing->cf)))
29      return -1;
30    ideal vv = id_Head(I,currRing);
31    if (i != -1) pDelete(&vv->m[i]);
32    d = scDimInt(vv, currRing->qideal);
33    if (rField_is_Ring_Z(currRing) && (i==-1)) d++;
34    idDelete(&vv);
35    return d;
36  }
37  else
38    d = scDimInt(I,currRing->qideal);
39  if (origin != r)
40    rChangeCurrRing(origin);
41  return d;
42}
43
44static void swapElements(ideal I, ideal J)
45{
46  assume(idSize(I)==idSize(J));
47
48  for (int i=0; i<idSize(I); i++)
49  {
50    poly cache = I->m[i];
51    I->m[i] = J->m[i];
52    J->m[i] = cache;
53  }
54
55  return;
56}
57
58static bool noExtraReduction(ideal I, ring r, number /*p*/)
59{
60  int n = rVar(r);
61  gfan::ZVector allOnes(n);
62  for (int i=0; i<n; i++)
63    allOnes[i] = 1;
64  ring rShortcut = rCopy0(r);
65
66  int* order = rShortcut->order;
67  int* block0 = rShortcut->block0;
68  int* block1 = rShortcut->block1;
69  int** wvhdl = rShortcut->wvhdl;
70
71  int h = rBlocks(r);
72  rShortcut->order = (int*) omAlloc0((h+1)*sizeof(int));
73  rShortcut->block0 = (int*) omAlloc0((h+1)*sizeof(int));
74  rShortcut->block1 = (int*) omAlloc0((h+1)*sizeof(int));
75  rShortcut->wvhdl = (int**) omAlloc0((h+1)*sizeof(int*));
76  rShortcut->order[0] = ringorder_a;
77  rShortcut->block0[0] = 1;
78  rShortcut->block1[0] = n;
79  bool overflow;
80  rShortcut->wvhdl[0] = ZVectorToIntStar(allOnes,overflow);
81  for (int i=1; i<=h; i++)
82  {
83    rShortcut->order[i] = order[i-1];
84    rShortcut->block0[i] = block0[i-1];
85    rShortcut->block1[i] = block1[i-1];
86    rShortcut->wvhdl[i] = wvhdl[i-1];
87  }
88
89  rComplete(rShortcut);
90  rTest(rShortcut);
91
92  omFree(order);
93  omFree(block0);
94  omFree(block1);
95  omFree(wvhdl);
96
97  int k = idSize(I);
98  ideal IShortcut = idInit(k);
99  nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
100  for (int i=0; i<k; i++)
101    IShortcut->m[i] = p_PermPoly(I->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
102
103  ideal JShortcut = gfanlib_kStd_wrapper(IShortcut,rShortcut);
104
105  ideal J = idInit(k);
106  nMapFunc outofShortcut = n_SetMap(rShortcut->cf,r->cf);
107  for (int i=0; i<k; i++)
108    J->m[i] = p_PermPoly(JShortcut->m[i],NULL,rShortcut,r,outofShortcut,NULL,0);
109
110  swapElements(I,J);
111  id_Delete(&IShortcut,rShortcut);
112  id_Delete(&JShortcut,rShortcut);
113  rDelete(rShortcut);
114  id_Delete(&J,r);
115  return false;
116}
117
118/**
119 * Initializes all relevant structures and information for the trivial valuation case,
120 * i.e. computing a tropical variety without any valuation.
121 */
122tropicalStrategy::tropicalStrategy(const ideal I, const ring r,
123                                   const bool completelyHomogeneous,
124                                   const bool completeSpace):
125  originalRing(rCopy(r)),
126  originalIdeal(id_Copy(I,r)),
127  expectedDimension(dim(originalIdeal,originalRing)),
128  linealitySpace(homogeneitySpace(originalIdeal,originalRing)),
129  startingRing(rCopy(originalRing)),
130  startingIdeal(id_Copy(originalIdeal,originalRing)),
131  uniformizingParameter(NULL),
132  shortcutRing(NULL),
133  onlyLowerHalfSpace(false),
134  weightAdjustingAlgorithm1(nonvalued_adjustWeightForHomogeneity),
135  weightAdjustingAlgorithm2(nonvalued_adjustWeightUnderHomogeneity),
136  extraReductionAlgorithm(noExtraReduction)
137{
138  assume(rField_is_Q(r) || rField_is_Zp(r));
139  if (!completelyHomogeneous)
140  {
141    weightAdjustingAlgorithm1 = valued_adjustWeightForHomogeneity;
142    weightAdjustingAlgorithm2 = valued_adjustWeightUnderHomogeneity;
143  }
144  if (!completeSpace)
145    onlyLowerHalfSpace = true;
146}
147
148/**
149 * Given a polynomial ring r over the rational numbers and a weighted ordering,
150 * returns a polynomial ring s over the integers with one extra variable, which is weighted -1.
151 */
152static ring constructStartingRing(ring r)
153{
154  assume(rField_is_Q(r));
155
156  ring s = rCopy0(r,FALSE,FALSE);
157  nKillChar(s->cf);
158  s->cf = nInitChar(n_Z,NULL);
159
160  int n = rVar(s)+1;
161  s->N = n;
162  char** oldNames = s->names;
163  s->names = (char**) omAlloc((n+1)*sizeof(char**));
164  s->names[0] = omStrDup("t");
165  for (int i=1; i<n; i++)
166    s->names[i] = oldNames[i-1];
167  omFree(oldNames);
168
169  s->order = (int*) omAlloc0(3*sizeof(int));
170  s->block0 = (int*) omAlloc0(3*sizeof(int));
171  s->block1 = (int*) omAlloc0(3*sizeof(int));
172  s->wvhdl = (int**) omAlloc0(3*sizeof(int**));
173  s->order[0] = ringorder_ws;
174  s->block0[0] = 1;
175  s->block1[0] = n;
176  s->wvhdl[0] = (int*) omAlloc(n*sizeof(int));
177  s->wvhdl[0][0] = 1;
178  for (int i=1; i<n; i++)
179    s->wvhdl[0][i] = -(r->wvhdl[0][i-1]);
180  s->order[1] = ringorder_C;
181
182  rComplete(s);
183  rTest(s);
184  return s;
185}
186
187static ring writeOrderingAsWP(ring r)
188{
189  assume(r->order[0]==ringorder_wp || r->order[0]==ringorder_dp);
190  if (r->order[0]==ringorder_dp)
191  {
192    ring s = rCopy0(r,FALSE,TRUE);
193    rComplete(s);
194    rTest(s);
195    return s;
196  }
197  return rCopy(r);
198}
199
200static ideal constructStartingIdeal(ideal originalIdeal, ring originalRing, number uniformizingParameter, ring startingRing)
201{
202  // construct p-t
203  poly g = p_One(startingRing);
204  p_SetCoeff(g,uniformizingParameter,startingRing);
205  pNext(g) = p_One(startingRing);
206  p_SetExp(pNext(g),1,1,startingRing);
207  p_SetCoeff(pNext(g),n_Init(-1,startingRing->cf),startingRing);
208  p_Setm(pNext(g),startingRing);
209  ideal pt = idInit(1);
210  pt->m[0] = g;
211
212  // map originalIdeal from originalRing into startingRing
213  int k = idSize(originalIdeal);
214  ideal J = idInit(k+1);
215  nMapFunc nMap = n_SetMap(originalRing->cf,startingRing->cf);
216  int n = rVar(originalRing);
217  int* shiftByOne = (int*) omAlloc((n+1)*sizeof(int));
218  for (int i=1; i<=n; i++)
219    shiftByOne[i]=i+1;
220  for (int i=0; i<k; i++)
221    J->m[i] = p_PermPoly(originalIdeal->m[i],shiftByOne,originalRing,startingRing,nMap,NULL,0);
222  omFreeSize(shiftByOne,(n+1)*sizeof(int));
223
224  ring origin = currRing;
225  rChangeCurrRing(startingRing);
226  ideal startingIdeal = kNF(pt,startingRing->qideal,J);
227  rChangeCurrRing(origin);
228  assume(startingIdeal->m[k]==NULL);
229  startingIdeal->m[k] = pt->m[0];
230  startingIdeal = gfanlib_kStd_wrapper(startingIdeal,startingRing);
231
232  id_Delete(&J,startingRing);
233  pt->m[0] = NULL;
234  id_Delete(&pt,startingRing);
235  return startingIdeal;
236}
237
238/***
239 * Initializes all relevant structures and information for the valued case,
240 * i.e. computing a tropical variety over the rational numbers with p-adic valuation
241 **/
242tropicalStrategy::tropicalStrategy(ideal J, number q, ring s):
243  originalRing(rCopy(s)),
244  originalIdeal(id_Copy(J,s)),
245  expectedDimension(dim(originalIdeal,originalRing)+1),
246  linealitySpace(gfan::ZCone()), // to come, see below
247  startingRing(NULL),            // to come, see below
248  startingIdeal(NULL),           // to come, see below
249  uniformizingParameter(NULL),   // to come, see below
250  shortcutRing(NULL),            // to come, see below
251  onlyLowerHalfSpace(true),
252  weightAdjustingAlgorithm1(valued_adjustWeightForHomogeneity),
253  weightAdjustingAlgorithm2(valued_adjustWeightUnderHomogeneity),
254  extraReductionAlgorithm(ppreduceInitially)
255{
256  /* assume that the ground field of the originalRing is Q */
257  assume(rField_is_Q(s));
258
259  /* replace Q with Z for the startingRing
260   * and add an extra variable for tracking the uniformizing parameter */
261  startingRing = constructStartingRing(originalRing);
262
263  /* map the uniformizing parameter into the new coefficient domain */
264  nMapFunc nMap = n_SetMap(originalRing->cf,startingRing->cf);
265  uniformizingParameter = nMap(q,originalRing->cf,startingRing->cf);
266
267  /* map the input ideal into the new polynomial ring */
268  startingIdeal = constructStartingIdeal(J,s,uniformizingParameter,startingRing);
269
270  linealitySpace = homogeneitySpace(startingIdeal,startingRing);
271
272  /* construct the shorcut ring */
273  shortcutRing = rCopy0(startingRing);
274  nKillChar(shortcutRing->cf);
275  shortcutRing->cf = nInitChar(n_Zp,(void*)(long)IsPrime(n_Int(uniformizingParameter,startingRing->cf)));
276  rComplete(shortcutRing);
277  rTest(shortcutRing);
278}
279
280tropicalStrategy::tropicalStrategy(const tropicalStrategy& currentStrategy):
281  originalRing(rCopy(currentStrategy.getOriginalRing())),
282  originalIdeal(id_Copy(currentStrategy.getOriginalIdeal(),currentStrategy.getOriginalRing())),
283  expectedDimension(currentStrategy.getExpectedDimension()),
284  linealitySpace(currentStrategy.getHomogeneitySpace()),
285  startingRing(rCopy(currentStrategy.getStartingRing())),
286  startingIdeal(id_Copy(currentStrategy.getStartingIdeal(),currentStrategy.getStartingRing())),
287  uniformizingParameter(NULL),
288  shortcutRing(NULL),
289  onlyLowerHalfSpace(currentStrategy.restrictToLowerHalfSpace()),
290  weightAdjustingAlgorithm1(currentStrategy.weightAdjustingAlgorithm1),
291  weightAdjustingAlgorithm2(currentStrategy.weightAdjustingAlgorithm2),
292  extraReductionAlgorithm(currentStrategy.extraReductionAlgorithm)
293{
294  if (originalRing) rTest(originalRing);
295  if (originalIdeal) id_Test(originalIdeal,originalRing);
296  if (startingRing) rTest(startingRing);
297  if (startingIdeal) id_Test(startingIdeal,startingRing);
298  if (currentStrategy.getUniformizingParameter())
299  {
300    uniformizingParameter = n_Copy(currentStrategy.getUniformizingParameter(),startingRing->cf);
301    n_Test(uniformizingParameter,startingRing->cf);
302  }
303  if (currentStrategy.getShortcutRing())
304  {
305    shortcutRing = rCopy(currentStrategy.getShortcutRing());
306    rTest(shortcutRing);
307  }
308}
309
310tropicalStrategy::~tropicalStrategy()
311{
312  if (originalRing) rTest(originalRing);
313  if (originalIdeal) id_Test(originalIdeal,originalRing);
314  if (startingRing) rTest(startingRing);
315  if (startingIdeal) id_Test(startingIdeal,startingRing);
316  if (uniformizingParameter) n_Test(uniformizingParameter,startingRing->cf);
317  if (shortcutRing) rTest(shortcutRing);
318
319  id_Delete(&originalIdeal,originalRing);
320  rDelete(originalRing);
321  if (startingIdeal) id_Delete(&startingIdeal,startingRing);
322  if (uniformizingParameter) n_Delete(&uniformizingParameter,startingRing->cf);
323  if (startingRing) rDelete(startingRing);
324  if (shortcutRing) rDelete(shortcutRing);
325}
326
327tropicalStrategy& tropicalStrategy::operator=(const tropicalStrategy& currentStrategy)
328{
329  originalRing = rCopy(currentStrategy.getOriginalRing());
330  originalIdeal = id_Copy(currentStrategy.getOriginalIdeal(),currentStrategy.getOriginalRing());
331  expectedDimension = currentStrategy.getExpectedDimension();
332  startingRing = rCopy(currentStrategy.getStartingRing());
333  startingIdeal = id_Copy(currentStrategy.getStartingIdeal(),currentStrategy.getStartingRing());
334  uniformizingParameter = n_Copy(currentStrategy.getUniformizingParameter(),startingRing->cf);
335  shortcutRing = rCopy(currentStrategy.getShortcutRing());
336  onlyLowerHalfSpace = currentStrategy.restrictToLowerHalfSpace();
337  weightAdjustingAlgorithm1 = currentStrategy.weightAdjustingAlgorithm1;
338  weightAdjustingAlgorithm2 = currentStrategy.weightAdjustingAlgorithm2;
339  extraReductionAlgorithm = currentStrategy.extraReductionAlgorithm;
340
341  if (originalRing) rTest(originalRing);
342  if (originalIdeal) id_Test(originalIdeal,originalRing);
343  if (startingRing) rTest(startingRing);
344  if (startingIdeal) id_Test(startingIdeal,startingRing);
345  if (uniformizingParameter) n_Test(uniformizingParameter,startingRing->cf);
346  if (shortcutRing) rTest(shortcutRing);
347
348  return *this;
349}
350
351void tropicalStrategy::putUniformizingBinomialInFront(ideal I, const ring r, const number q) const
352{
353  poly p = p_One(r);
354  p_SetCoeff(p,q,r);
355  poly t = p_One(r);
356  p_SetExp(t,1,1,r);
357  p_Setm(t,r);
358  poly pt = p_Add_q(p,p_Neg(t,r),r);
359
360  int k = idSize(I);
361  int l;
362  for (l=0; l<k; l++)
363  {
364    if (p_EqualPolys(I->m[l],pt,r))
365      break;
366  }
367  p_Delete(&pt,r);
368
369  if (l>1)
370  {
371    pt = I->m[l];
372    for (int i=l; i>0; i--)
373      I->m[l] = I->m[l-1];
374    I->m[0] = pt;
375    pt = NULL;
376  }
377  return;
378}
379
380bool tropicalStrategy::reduce(ideal I, const ring r) const
381{
382  rTest(r);
383  id_Test(I,r);
384
385  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
386  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
387  bool b = extraReductionAlgorithm(I,r,p);
388  // putUniformizingBinomialInFront(I,r,p);
389  n_Delete(&p,r->cf);
390
391  return b;
392}
393
394void tropicalStrategy::pReduce(ideal I, const ring r) const
395{
396  rTest(r);
397  id_Test(I,r);
398
399  if (isValuationTrivial())
400    return;
401
402  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
403  number p = identity(uniformizingParameter,startingRing->cf,r->cf);
404  ::pReduce(I,p,r);
405  n_Delete(&p,r->cf);
406
407  return;
408}
409
410ring tropicalStrategy::getShortcutRingPrependingWeight(const ring r, const gfan::ZVector v) const
411{
412  ring rShortcut = rCopy0(r);
413
414  // save old ordering
415  int* order = rShortcut->order;
416  int* block0 = rShortcut->block0;
417  int* block1 = rShortcut->block1;
418  int** wvhdl = rShortcut->wvhdl;
419
420  // adjust weight and create new ordering
421  gfan::ZVector w = adjustWeightForHomogeneity(v);
422  int h = rBlocks(r); int n = rVar(r);
423  rShortcut->order = (int*) omAlloc0((h+1)*sizeof(int));
424  rShortcut->block0 = (int*) omAlloc0((h+1)*sizeof(int));
425  rShortcut->block1 = (int*) omAlloc0((h+1)*sizeof(int));
426  rShortcut->wvhdl = (int**) omAlloc0((h+1)*sizeof(int*));
427  rShortcut->order[0] = ringorder_a;
428  rShortcut->block0[0] = 1;
429  rShortcut->block1[0] = n;
430  bool overflow;
431  rShortcut->wvhdl[0] = ZVectorToIntStar(w,overflow);
432  for (int i=1; i<=h; i++)
433  {
434    rShortcut->order[i] = order[i-1];
435    rShortcut->block0[i] = block0[i-1];
436    rShortcut->block1[i] = block1[i-1];
437    rShortcut->wvhdl[i] = wvhdl[i-1];
438  }
439
440  // if valuation non-trivial, change coefficient ring to residue field
441  if (isValuationNonTrivial())
442  {
443    nKillChar(rShortcut->cf);
444    rShortcut->cf = nCopyCoeff(shortcutRing->cf);
445  }
446  rComplete(rShortcut);
447  rTest(rShortcut);
448
449  // delete old ordering
450  omFree(order);
451  omFree(block0);
452  omFree(block1);
453  omFree(wvhdl);
454
455  return rShortcut;
456}
457
458poly tropicalStrategy::checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector w) const
459{
460  // prepend extra weight for homogeneity
461  // switch to residue field if valuation is non trivial
462  ring rShortcut = getShortcutRingPrependingWeight(r,w);
463
464  // compute the initial ideal and map it into the constructed ring
465  // if switched to residue field, remove possibly 0 elements
466  ideal inI = initial(I,r,w);
467  int k = idSize(inI);
468  ideal inIShortcut = idInit(k);
469  nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
470  for (int i=0; i<k; i++)
471    inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
472  if (isValuationNonTrivial())
473    idSkipZeroes(inIShortcut);
474
475  // check initial ideal for monomial and
476  // if it exsists, return a copy of the monomial in the input ring
477  poly p = checkForMonomialViaSuddenSaturation(inIShortcut,rShortcut);
478  poly monomial = NULL; int n = rVar(r);
479  if (p!=NULL)
480  {
481    monomial=p_One(r);
482    for (int i=1; i<=n; i++)
483      p_SetExp(monomial,i,p_GetExp(p,i,rShortcut),r);
484    p_Delete(&p,rShortcut);
485  }
486  id_Delete(&inI,r);
487  id_Delete(&inIShortcut,rShortcut);
488  rDelete(rShortcut);
489  return monomial;
490}
491
492ring tropicalStrategy::copyAndChangeCoefficientRing(const ring r) const
493{
494  ring rShortcut = rCopy0(r);
495  nKillChar(rShortcut->cf);
496  rShortcut->cf = nCopyCoeff(shortcutRing->cf);
497  rComplete(rShortcut);
498  rTest(rShortcut);
499  return rShortcut;
500}
501
502ideal tropicalStrategy::computeWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
503{
504  // if the valuation is trivial and the ring and ideal have not been extended,
505  // then it is sufficient to return the difference between the elements of inJ
506  // and their normal forms with respect to I and r
507  if (isValuationTrivial())
508    return witness(inJ,I,r);
509  // if the valuation is non-trivial and the ring and ideal have been extended,
510  // then we can make a shortcut through the residue field
511  else
512  {
513    assume(idSize(inI)==idSize(I));
514    int uni = findPositionOfUniformizingBinomial(I,r);
515    assume(uni>=0);
516    /**
517     * change ground ring into finite field
518     * and map the data into it
519     */
520    ring rShortcut = copyAndChangeCoefficientRing(r);
521
522    int k = idSize(inJ);
523    int l = idSize(I);
524    ideal inJShortcut = idInit(k);
525    ideal inIShortcut = idInit(l);
526    nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
527    for (int i=0; i<k; i++)
528      inJShortcut->m[i] = p_PermPoly(inJ->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
529    for (int j=0; j<l; j++)
530      inIShortcut->m[j] = p_PermPoly(inI->m[j],NULL,r,rShortcut,takingResidues,NULL,0);
531    id_Test(inJShortcut,rShortcut);
532    id_Test(inIShortcut,rShortcut);
533
534    /**
535     * Compute a division with remainder over the finite field
536     * and map the result back to r
537     */
538    matrix QShortcut = divisionDiscardingRemainder(inJShortcut,inIShortcut,rShortcut);
539    matrix Q = mpNew(l,k);
540    nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
541    for (int ij=k*l-1; ij>=0; ij--)
542      Q->m[ij] = p_PermPoly(QShortcut->m[ij],NULL,rShortcut,r,takingRepresentatives,NULL,0);
543
544    nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
545    number p = identity(uniformizingParameter,startingRing->cf,r->cf);
546
547    /**
548     * Compute the normal forms
549     */
550    ideal J = idInit(k);
551    for (int j=0; j<k; j++)
552    {
553      poly q0 = p_Copy(inJ->m[j],r);
554      for (int i=0; i<l; i++)
555      {
556        poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
557        poly inIi = p_Copy(inI->m[i],r);
558        q0 = p_Add_q(q0,p_Neg(p_Mult_q(qij,inIi,r),r),r);
559      }
560      q0 = p_Div_nn(q0,p,r);
561      poly q0g0 = p_Mult_q(q0,p_Copy(I->m[uni],r),r);
562      // q0 = NULL;
563      poly qigi = NULL;
564      for (int i=0; i<l; i++)
565      {
566        poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
567        // poly inIi = p_Copy(I->m[i],r);
568        poly Ii = p_Copy(I->m[i],r);
569        qigi = p_Add_q(qigi,p_Mult_q(qij,Ii,r),r);
570      }
571      J->m[j] = p_Add_q(q0g0,qigi,r);
572    }
573
574    id_Delete(&inIShortcut,rShortcut);
575    id_Delete(&inJShortcut,rShortcut);
576    mp_Delete(&QShortcut,rShortcut);
577    rDelete(rShortcut);
578    mp_Delete(&Q,r);
579    n_Delete(&p,r->cf);
580    return J;
581  }
582}
583
584ideal tropicalStrategy::computeStdOfInitialIdeal(const ideal inI, const ring r) const
585{
586  // if valuation trivial, then compute std as usual
587  if (isValuationTrivial())
588    return gfanlib_kStd_wrapper(inI,r);
589
590  // if valuation non-trivial, then uniformizing parameter is in ideal
591  // so switch to residue field first and compute standard basis over the residue field
592  ring rShortcut = copyAndChangeCoefficientRing(r);
593  nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
594  int k = idSize(inI);
595  ideal inIShortcut = idInit(k);
596  for (int i=0; i<k; i++)
597    inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
598  ideal inJShortcut = gfanlib_kStd_wrapper(inIShortcut,rShortcut);
599
600  // and lift the result back to the ring with valuation
601  nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
602  k = idSize(inJShortcut);
603  ideal inJ = idInit(k+1);
604  inJ->m[0] = p_One(r);
605  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
606  p_SetCoeff(inJ->m[0],identity(uniformizingParameter,startingRing->cf,r->cf),r);
607  for (int i=0; i<k; i++)
608    inJ->m[i+1] = p_PermPoly(inJShortcut->m[i],NULL,rShortcut,r,takingRepresentatives,NULL,0);
609
610  id_Delete(&inJShortcut,rShortcut);
611  id_Delete(&inIShortcut,rShortcut);
612  rDelete(rShortcut);
613  return inJ;
614}
615
616ideal tropicalStrategy::computeLift(const ideal inJs, const ring s, const ideal inIr, const ideal Ir, const ring r) const
617{
618  int k = idSize(inJs);
619  ideal inJr = idInit(k);
620  nMapFunc identitysr = n_SetMap(s->cf,r->cf);
621  for (int i=0; i<k; i++)
622    inJr->m[i] = p_PermPoly(inJs->m[i],NULL,s,r,identitysr,NULL,0);
623
624  ideal Jr = computeWitness(inJr,inIr,Ir,r);
625  nMapFunc identityrs = n_SetMap(r->cf,s->cf);
626  ideal Js = idInit(k);
627  for (int i=0; i<k; i++)
628    Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identityrs,NULL,0);
629  return Js;
630}
631
632static void deleteOrdering(ring r)
633{
634  if (r->order != NULL)
635  {
636    int i=rBlocks(r);
637    assume(r->block0 != NULL && r->block1 != NULL && r->wvhdl != NULL);
638    /* delete order */
639    omFreeSize((ADDRESS)r->order,i*sizeof(int));
640    omFreeSize((ADDRESS)r->block0,i*sizeof(int));
641    omFreeSize((ADDRESS)r->block1,i*sizeof(int));
642    /* delete weights */
643    for (int j=0; j<i; j++)
644      if (r->wvhdl[j]!=NULL)
645        omFree(r->wvhdl[j]);
646    omFreeSize((ADDRESS)r->wvhdl,i*sizeof(int *));
647  }
648  else
649    assume(r->block0 == NULL && r->block1 == NULL && r->wvhdl == NULL);
650  return;
651}
652
653ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector w, const gfan::ZVector v) const
654{
655  // copy shortcutRing and change to desired ordering
656  bool ok;
657  ring s = rCopy0(r);
658  int n = rVar(s);
659  deleteOrdering(s);
660  gfan::ZVector wAdjusted = adjustWeightForHomogeneity(w);
661  gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(v,wAdjusted);
662  s->order = (int*) omAlloc0(5*sizeof(int));
663  s->block0 = (int*) omAlloc0(5*sizeof(int));
664  s->block1 = (int*) omAlloc0(5*sizeof(int));
665  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
666  s->order[0] = ringorder_a;
667  s->block0[0] = 1;
668  s->block1[0] = n;
669  s->wvhdl[0] = ZVectorToIntStar(wAdjusted,ok);
670  s->order[1] = ringorder_a;
671  s->block0[1] = 1;
672  s->block1[1] = n;
673  s->wvhdl[1] = ZVectorToIntStar(vAdjusted,ok);
674  s->order[2] = ringorder_dp;
675  s->block0[2] = 1;
676  s->block1[2] = n;
677  s->order[3] = ringorder_C;
678  rComplete(s);
679  rTest(s);
680
681  return s;
682}
683
684ring tropicalStrategy::copyAndChangeOrderingLS(const ring r, const gfan::ZVector w, const gfan::ZVector v) const
685{
686  // copy shortcutRing and change to desired ordering
687  bool ok;
688  ring s = rCopy0(r);
689  int n = rVar(s);
690  deleteOrdering(s);
691  s->order = (int*) omAlloc0(5*sizeof(int));
692  s->block0 = (int*) omAlloc0(5*sizeof(int));
693  s->block1 = (int*) omAlloc0(5*sizeof(int));
694  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
695  s->order[0] = ringorder_a;
696  s->block0[0] = 1;
697  s->block1[0] = n;
698  s->wvhdl[0] = ZVectorToIntStar(w,ok);
699  s->order[1] = ringorder_a;
700  s->block0[1] = 1;
701  s->block1[1] = n;
702  s->wvhdl[1] = ZVectorToIntStar(v,ok);
703  s->order[2] = ringorder_dp;
704  s->block0[2] = 1;
705  s->block1[2] = n;
706  s->order[3] = ringorder_C;
707  rComplete(s);
708  rTest(s);
709
710  return s;
711}
712
713std::pair<ideal,ring> tropicalStrategy::computeFlip(const ideal Ir, const ring r,
714                                                    const gfan::ZVector interiorPoint,
715                                                    const gfan::ZVector facetNormal) const
716{
717  assume(isValuationTrivial() || interiorPoint[0].sign()<0);
718  assume(checkForUniformizingBinomial(Ir,r));
719
720  // get a generating system of the initial ideal
721  // and compute a standard basis with respect to adjacent ordering
722  ideal inIr = initial(Ir,r,interiorPoint);
723  ring sAdjusted = copyAndChangeOrderingWP(r,interiorPoint,facetNormal);
724  nMapFunc identity = n_SetMap(r->cf,sAdjusted->cf);
725  int k = idSize(Ir);
726  ideal inIsAdjusted = idInit(k);
727  for (int i=0; i<k; i++)
728    inIsAdjusted->m[i] = p_PermPoly(inIr->m[i],NULL,r,sAdjusted,identity,NULL,0);
729  ideal inJsAdjusted = computeStdOfInitialIdeal(inIsAdjusted,sAdjusted);
730
731  // find witnesses of the new standard basis elements of the initial ideal
732  // with the help of the old standard basis of the ideal
733  k = idSize(inJsAdjusted);
734  ideal inJr = idInit(k);
735  identity = n_SetMap(sAdjusted->cf,r->cf);
736  for (int i=0; i<k; i++)
737    inJr->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,sAdjusted,r,identity,NULL,0);
738
739  ideal Jr = computeWitness(inJr,inIr,Ir,r);
740  ring s = copyAndChangeOrderingLS(r,interiorPoint,facetNormal);
741  identity = n_SetMap(r->cf,s->cf);
742  ideal Js = idInit(k);
743  for (int i=0; i<k; i++)
744    Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identity,NULL,0);
745
746  // this->reduce(Jr,r);
747  // cleanup
748  id_Delete(&inIsAdjusted,sAdjusted);
749  id_Delete(&inJsAdjusted,sAdjusted);
750  rDelete(sAdjusted);
751  id_Delete(&inIr,r);
752  id_Delete(&Jr,r);
753  id_Delete(&inJr,r);
754
755  assume(isValuationTrivial() || isOrderingLocalInT(s));
756  return std::make_pair(Js,s);
757}
758
759
760bool tropicalStrategy::checkForUniformizingBinomial(const ideal I, const ring r) const
761{
762  // if the valuation is trivial,
763  // then there is no special condition the first generator has to fullfill
764  if (isValuationTrivial())
765    return true;
766
767  // if the valuation is non-trivial then checks if the first generator is p-t
768  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
769  poly p = p_One(r);
770  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
771  poly t = p_One(r);
772  p_SetExp(t,1,1,r);
773  p_Setm(t,r);
774  poly pt = p_Add_q(p,p_Neg(t,r),r);
775
776  for (int i=0; i<idSize(I); i++)
777  {
778    if (p_EqualPolys(I->m[i],pt,r))
779    {
780      p_Delete(&pt,r);
781      return true;
782    }
783  }
784  p_Delete(&pt,r);
785  return false;
786}
787
788int tropicalStrategy::findPositionOfUniformizingBinomial(const ideal I, const ring r) const
789{
790  assume(isValuationNonTrivial());
791
792  // if the valuation is non-trivial then checks if the first generator is p-t
793  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
794  poly p = p_One(r);
795  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
796  poly t = p_One(r);
797  p_SetExp(t,1,1,r);
798  p_Setm(t,r);
799  poly pt = p_Add_q(p,p_Neg(t,r),r);
800
801  for (int i=0; i<idSize(I); i++)
802  {
803    if (p_EqualPolys(I->m[i],pt,r))
804    {
805      p_Delete(&pt,r);
806      return i;
807    }
808  }
809  p_Delete(&pt,r);
810  return -1;
811}
812
813bool tropicalStrategy::checkForUniformizingParameter(const ideal inI, const ring r) const
814{
815  // if the valuation is trivial,
816  // then there is no special condition the first generator has to fullfill
817  if (isValuationTrivial())
818    return true;
819
820  // if the valuation is non-trivial then checks if the first generator is p
821  if (inI->m[0]==NULL)
822    return false;
823  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
824  poly p = p_One(r);
825  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
826
827  for (int i=0; i<idSize(inI); i++)
828  {
829    if (p_EqualPolys(inI->m[i],p,r))
830    {
831      p_Delete(&p,r);
832      return true;
833    }
834  }
835  p_Delete(&p,r);
836  return false;
837}
838
839#ifndef NDEBUG
840tropicalStrategy::tropicalStrategy():
841  originalRing(NULL),
842  originalIdeal(NULL),
843  expectedDimension(NULL),
844  linealitySpace(gfan::ZCone()), // to come, see below
845  startingRing(NULL),            // to come, see below
846  startingIdeal(NULL),           // to come, see below
847  uniformizingParameter(NULL),   // to come, see below
848  shortcutRing(NULL),            // to come, see below
849  onlyLowerHalfSpace(false),
850  weightAdjustingAlgorithm1(NULL),
851  weightAdjustingAlgorithm2(NULL),
852  extraReductionAlgorithm(NULL)
853{
854}
855
856tropicalStrategy tropicalStrategy::debugStrategy(const ideal startIdeal, number unifParameter, ring startRing)
857{
858  tropicalStrategy debug;
859  debug.originalRing = rCopy(startRing);
860  debug.originalIdeal = id_Copy(startIdeal,startRing);
861  debug.startingRing = rCopy(startRing);
862  debug.startingIdeal = id_Copy(startIdeal,startRing);
863  debug.uniformizingParameter = n_Copy(unifParameter,startRing->cf);
864
865  debug.shortcutRing = rCopy0(startRing);
866  nKillChar(debug.shortcutRing->cf);
867  debug.shortcutRing->cf = nInitChar(n_Zp,(void*)(long)IsPrime(n_Int(unifParameter,startRing->cf)));
868  rComplete(debug.shortcutRing);
869  rTest(debug.shortcutRing);
870
871  debug.onlyLowerHalfSpace = true;
872  debug.weightAdjustingAlgorithm1 = valued_adjustWeightForHomogeneity;
873  debug.weightAdjustingAlgorithm2 = valued_adjustWeightUnderHomogeneity;
874  debug.extraReductionAlgorithm = ppreduceInitially;
875
876  return debug;
877}
878
879BOOLEAN computeWitnessDebug(leftv res, leftv args)
880{
881  leftv u = args;
882  if ((u!=NULL) && (u->Typ()==IDEAL_CMD))
883  {
884    leftv v = u->next;
885    if ((v!=NULL) && (v->Typ()==IDEAL_CMD))
886    {
887      leftv w = v->next;
888      if ((w!=NULL) && (w->Typ()==IDEAL_CMD))
889      {
890        leftv x = w->next;
891        if ((x!=NULL) && (x->Typ()==NUMBER_CMD))
892        {
893          omUpdateInfo();
894          Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
895
896          ideal inJ = (ideal) u->CopyD();
897          ideal inI = (ideal) v->CopyD();
898          ideal I = (ideal) w->CopyD();
899          number p = (number) x->CopyD();
900          tropicalStrategy debug = tropicalStrategy::debugStrategy(I,p,currRing);
901          ideal J = debug.computeWitness(inJ,inI,I,currRing);
902          id_Delete(&inJ,currRing);
903          id_Delete(&inI,currRing);
904          id_Delete(&I,currRing);
905          n_Delete(&p,currRing->cf);
906          res->rtyp = IDEAL_CMD;
907          res->data = (char*) J;
908          return FALSE;
909        }
910      }
911    }
912  }
913  return TRUE;
914}
915
916BOOLEAN computeFlipDebug(leftv res, leftv args)
917{
918  leftv u = args;
919  if ((u!=NULL) && (u->Typ()==IDEAL_CMD))
920  {
921    leftv v = u->next;
922    if ((v!=NULL) && (v->Typ()==NUMBER_CMD))
923    {
924      leftv w = v->next;
925      if ((w!=NULL) && (w->Typ()==BIGINTMAT_CMD))
926      {
927        leftv x = w->next;
928        if ((x!=NULL) && (x->Typ()==BIGINTMAT_CMD))
929        {
930          omUpdateInfo();
931          Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
932
933          ideal I = (ideal) u->CopyD();
934          number p = (number) v->CopyD();
935          bigintmat* interiorPoint0 = (bigintmat*) w->CopyD();
936          bigintmat* facetNormal0 = (bigintmat*) x->CopyD();
937          tropicalStrategy debug = tropicalStrategy::debugStrategy(I,p,currRing);
938
939          gfan::ZVector* interiorPoint = bigintmatToZVector(interiorPoint0);
940          gfan::ZVector* facetNormal = bigintmatToZVector(facetNormal0);
941          std::pair<ideal,ring> Js = debug.computeFlip(I,currRing,*interiorPoint,*facetNormal);
942          ideal J = Js.first;
943          ring s = Js.second;
944
945          id_Delete(&J,s);
946          rDelete(s);
947
948          id_Delete(&I,currRing);
949          n_Delete(&p,currRing->cf);
950          delete interiorPoint0;
951          delete facetNormal0;
952          delete interiorPoint;
953          delete facetNormal;
954
955          res->rtyp = NONE;
956          res->data = NULL;
957          return FALSE;
958        }
959      }
960    }
961  }
962  WerrorS("computeFlipDebug: unexpected parameters");
963  return TRUE;
964}
965#endif
Note: See TracBrowser for help on using the repository browser.