source: git/Singular/dyn_modules/gfanlib/tropicalStrategy.cc @ 08b7c2a

spielwiese
Last change on this file since 08b7c2a was 08b7c2a, checked in by Yue Ren <ren@…>, 9 years ago
chg: removed redundant witness computations
  • Property mode set to 100644
File size: 30.1 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
458std::pair<poly,int> tropicalStrategy::checkInitialIdealForMonomial(const ideal I, const ring r, const gfan::ZVector &w) const
459{
460  // quick check whether I already contains an ideal
461  int k = idSize(I);
462  for (int i=0; i<k; i++)
463  {
464    poly g = I->m[i];
465    if (pNext(g)==NULL && (isValuationTrivial() || n_IsOne(p_GetCoeff(g,r),r->cf)))
466      return std::pair<poly,int>(g,i);
467  }
468
469  ring rShortcut;
470  ideal inIShortcut;
471  if (w.size()>0)
472  {
473    // if needed, prepend extra weight for homogeneity
474    // switch to residue field if valuation is non trivial
475    rShortcut = getShortcutRingPrependingWeight(r,w);
476
477    // compute the initial ideal and map it into the constructed ring
478    // if switched to residue field, remove possibly 0 elements
479    ideal inI = initial(I,r,w);
480    inIShortcut = idInit(k);
481    nMapFunc intoShortcut = n_SetMap(r->cf,rShortcut->cf);
482    for (int i=0; i<k; i++)
483      inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,intoShortcut,NULL,0);
484    if (isValuationNonTrivial())
485      idSkipZeroes(inIShortcut);
486    id_Delete(&inI,r);
487  }
488  else
489  {
490    rShortcut = r;
491    inIShortcut = I;
492  }
493
494  // check initial ideal for monomial and
495  // if it exsists, return a copy of the monomial in the input ring
496  poly p = checkForMonomialViaSuddenSaturation(inIShortcut,rShortcut);
497  poly monomial = NULL;
498  if (p!=NULL)
499  {
500    monomial=p_One(r);
501    for (int i=1; i<=rVar(r); i++)
502      p_SetExp(monomial,i,p_GetExp(p,i,rShortcut),r);
503    p_Delete(&p,rShortcut);
504  }
505
506
507  if (w.size()>0)
508  {
509    // if needed, cleanup
510    id_Delete(&inIShortcut,rShortcut);
511    rDelete(rShortcut);
512  }
513  return std::pair<poly,int>(monomial,-1);
514}
515
516ring tropicalStrategy::copyAndChangeCoefficientRing(const ring r) const
517{
518  ring rShortcut = rCopy0(r);
519  nKillChar(rShortcut->cf);
520  rShortcut->cf = nCopyCoeff(shortcutRing->cf);
521  rComplete(rShortcut);
522  rTest(rShortcut);
523  return rShortcut;
524}
525
526ideal tropicalStrategy::computeWitness(const ideal inJ, const ideal inI, const ideal I, const ring r) const
527{
528  // if the valuation is trivial and the ring and ideal have not been extended,
529  // then it is sufficient to return the difference between the elements of inJ
530  // and their normal forms with respect to I and r
531  if (isValuationTrivial())
532    return witness(inJ,I,r);
533  // if the valuation is non-trivial and the ring and ideal have been extended,
534  // then we can make a shortcut through the residue field
535  else
536  {
537    assume(idSize(inI)==idSize(I));
538    int uni = findPositionOfUniformizingBinomial(I,r);
539    assume(uni>=0);
540    /**
541     * change ground ring into finite field
542     * and map the data into it
543     */
544    ring rShortcut = copyAndChangeCoefficientRing(r);
545
546    int k = idSize(inJ);
547    int l = idSize(I);
548    ideal inJShortcut = idInit(k);
549    ideal inIShortcut = idInit(l);
550    nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
551    for (int i=0; i<k; i++)
552      inJShortcut->m[i] = p_PermPoly(inJ->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
553    for (int j=0; j<l; j++)
554      inIShortcut->m[j] = p_PermPoly(inI->m[j],NULL,r,rShortcut,takingResidues,NULL,0);
555    id_Test(inJShortcut,rShortcut);
556    id_Test(inIShortcut,rShortcut);
557
558    /**
559     * Compute a division with remainder over the finite field
560     * and map the result back to r
561     */
562    matrix QShortcut = divisionDiscardingRemainder(inJShortcut,inIShortcut,rShortcut);
563    matrix Q = mpNew(l,k);
564    nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
565    for (int ij=k*l-1; ij>=0; ij--)
566      Q->m[ij] = p_PermPoly(QShortcut->m[ij],NULL,rShortcut,r,takingRepresentatives,NULL,0);
567
568    nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
569    number p = identity(uniformizingParameter,startingRing->cf,r->cf);
570
571    /**
572     * Compute the normal forms
573     */
574    ideal J = idInit(k);
575    for (int j=0; j<k; j++)
576    {
577      poly q0 = p_Copy(inJ->m[j],r);
578      for (int i=0; i<l; i++)
579      {
580        poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
581        poly inIi = p_Copy(inI->m[i],r);
582        q0 = p_Add_q(q0,p_Neg(p_Mult_q(qij,inIi,r),r),r);
583      }
584      q0 = p_Div_nn(q0,p,r);
585      poly q0g0 = p_Mult_q(q0,p_Copy(I->m[uni],r),r);
586      // q0 = NULL;
587      poly qigi = NULL;
588      for (int i=0; i<l; i++)
589      {
590        poly qij = p_Copy(MATELEM(Q,i+1,j+1),r);
591        // poly inIi = p_Copy(I->m[i],r);
592        poly Ii = p_Copy(I->m[i],r);
593        qigi = p_Add_q(qigi,p_Mult_q(qij,Ii,r),r);
594      }
595      J->m[j] = p_Add_q(q0g0,qigi,r);
596    }
597
598    id_Delete(&inIShortcut,rShortcut);
599    id_Delete(&inJShortcut,rShortcut);
600    mp_Delete(&QShortcut,rShortcut);
601    rDelete(rShortcut);
602    mp_Delete(&Q,r);
603    n_Delete(&p,r->cf);
604    return J;
605  }
606}
607
608ideal tropicalStrategy::computeStdOfInitialIdeal(const ideal inI, const ring r) const
609{
610  // if valuation trivial, then compute std as usual
611  if (isValuationTrivial())
612    return gfanlib_kStd_wrapper(inI,r);
613
614  // if valuation non-trivial, then uniformizing parameter is in ideal
615  // so switch to residue field first and compute standard basis over the residue field
616  ring rShortcut = copyAndChangeCoefficientRing(r);
617  nMapFunc takingResidues = n_SetMap(r->cf,rShortcut->cf);
618  int k = idSize(inI);
619  ideal inIShortcut = idInit(k);
620  for (int i=0; i<k; i++)
621    inIShortcut->m[i] = p_PermPoly(inI->m[i],NULL,r,rShortcut,takingResidues,NULL,0);
622  ideal inJShortcut = gfanlib_kStd_wrapper(inIShortcut,rShortcut);
623
624  // and lift the result back to the ring with valuation
625  nMapFunc takingRepresentatives = n_SetMap(rShortcut->cf,r->cf);
626  k = idSize(inJShortcut);
627  ideal inJ = idInit(k+1);
628  inJ->m[0] = p_One(r);
629  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
630  p_SetCoeff(inJ->m[0],identity(uniformizingParameter,startingRing->cf,r->cf),r);
631  for (int i=0; i<k; i++)
632    inJ->m[i+1] = p_PermPoly(inJShortcut->m[i],NULL,rShortcut,r,takingRepresentatives,NULL,0);
633
634  id_Delete(&inJShortcut,rShortcut);
635  id_Delete(&inIShortcut,rShortcut);
636  rDelete(rShortcut);
637  return inJ;
638}
639
640ideal tropicalStrategy::computeLift(const ideal inJs, const ring s, const ideal inIr, const ideal Ir, const ring r) const
641{
642  int k = idSize(inJs);
643  ideal inJr = idInit(k);
644  nMapFunc identitysr = n_SetMap(s->cf,r->cf);
645  for (int i=0; i<k; i++)
646    inJr->m[i] = p_PermPoly(inJs->m[i],NULL,s,r,identitysr,NULL,0);
647
648  ideal Jr = computeWitness(inJr,inIr,Ir,r);
649  nMapFunc identityrs = n_SetMap(r->cf,s->cf);
650  ideal Js = idInit(k);
651  for (int i=0; i<k; i++)
652    Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identityrs,NULL,0);
653  return Js;
654}
655
656static void deleteOrdering(ring r)
657{
658  if (r->order != NULL)
659  {
660    int i=rBlocks(r);
661    assume(r->block0 != NULL && r->block1 != NULL && r->wvhdl != NULL);
662    /* delete order */
663    omFreeSize((ADDRESS)r->order,i*sizeof(int));
664    omFreeSize((ADDRESS)r->block0,i*sizeof(int));
665    omFreeSize((ADDRESS)r->block1,i*sizeof(int));
666    /* delete weights */
667    for (int j=0; j<i; j++)
668      if (r->wvhdl[j]!=NULL)
669        omFree(r->wvhdl[j]);
670    omFreeSize((ADDRESS)r->wvhdl,i*sizeof(int *));
671  }
672  else
673    assume(r->block0 == NULL && r->block1 == NULL && r->wvhdl == NULL);
674  return;
675}
676
677ring tropicalStrategy::copyAndChangeOrderingWP(const ring r, const gfan::ZVector w, const gfan::ZVector v) const
678{
679  // copy shortcutRing and change to desired ordering
680  bool ok;
681  ring s = rCopy0(r);
682  int n = rVar(s);
683  deleteOrdering(s);
684  gfan::ZVector wAdjusted = adjustWeightForHomogeneity(w);
685  gfan::ZVector vAdjusted = adjustWeightUnderHomogeneity(v,wAdjusted);
686  s->order = (int*) omAlloc0(5*sizeof(int));
687  s->block0 = (int*) omAlloc0(5*sizeof(int));
688  s->block1 = (int*) omAlloc0(5*sizeof(int));
689  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
690  s->order[0] = ringorder_a;
691  s->block0[0] = 1;
692  s->block1[0] = n;
693  s->wvhdl[0] = ZVectorToIntStar(wAdjusted,ok);
694  s->order[1] = ringorder_a;
695  s->block0[1] = 1;
696  s->block1[1] = n;
697  s->wvhdl[1] = ZVectorToIntStar(vAdjusted,ok);
698  s->order[2] = ringorder_dp;
699  s->block0[2] = 1;
700  s->block1[2] = n;
701  s->order[3] = ringorder_C;
702  rComplete(s);
703  rTest(s);
704
705  return s;
706}
707
708ring tropicalStrategy::copyAndChangeOrderingLS(const ring r, const gfan::ZVector w, const gfan::ZVector v) const
709{
710  // copy shortcutRing and change to desired ordering
711  bool ok;
712  ring s = rCopy0(r);
713  int n = rVar(s);
714  deleteOrdering(s);
715  s->order = (int*) omAlloc0(5*sizeof(int));
716  s->block0 = (int*) omAlloc0(5*sizeof(int));
717  s->block1 = (int*) omAlloc0(5*sizeof(int));
718  s->wvhdl = (int**) omAlloc0(5*sizeof(int**));
719  s->order[0] = ringorder_a;
720  s->block0[0] = 1;
721  s->block1[0] = n;
722  s->wvhdl[0] = ZVectorToIntStar(w,ok);
723  s->order[1] = ringorder_a;
724  s->block0[1] = 1;
725  s->block1[1] = n;
726  s->wvhdl[1] = ZVectorToIntStar(v,ok);
727  s->order[2] = ringorder_dp;
728  s->block0[2] = 1;
729  s->block1[2] = n;
730  s->order[3] = ringorder_C;
731  rComplete(s);
732  rTest(s);
733
734  return s;
735}
736
737std::pair<ideal,ring> tropicalStrategy::computeFlip(const ideal Ir, const ring r,
738                                                    const gfan::ZVector interiorPoint,
739                                                    const gfan::ZVector facetNormal) const
740{
741  assume(isValuationTrivial() || interiorPoint[0].sign()<0);
742  assume(checkForUniformizingBinomial(Ir,r));
743
744  // get a generating system of the initial ideal
745  // and compute a standard basis with respect to adjacent ordering
746  ideal inIr = initial(Ir,r,interiorPoint);
747  ring sAdjusted = copyAndChangeOrderingWP(r,interiorPoint,facetNormal);
748  nMapFunc identity = n_SetMap(r->cf,sAdjusted->cf);
749  int k = idSize(Ir);
750  ideal inIsAdjusted = idInit(k);
751  for (int i=0; i<k; i++)
752    inIsAdjusted->m[i] = p_PermPoly(inIr->m[i],NULL,r,sAdjusted,identity,NULL,0);
753  ideal inJsAdjusted = computeStdOfInitialIdeal(inIsAdjusted,sAdjusted);
754
755  // find witnesses of the new standard basis elements of the initial ideal
756  // with the help of the old standard basis of the ideal
757  k = idSize(inJsAdjusted);
758  ideal inJr = idInit(k);
759  identity = n_SetMap(sAdjusted->cf,r->cf);
760  for (int i=0; i<k; i++)
761    inJr->m[i] = p_PermPoly(inJsAdjusted->m[i],NULL,sAdjusted,r,identity,NULL,0);
762
763  ideal Jr = computeWitness(inJr,inIr,Ir,r);
764  ring s = copyAndChangeOrderingLS(r,interiorPoint,facetNormal);
765  identity = n_SetMap(r->cf,s->cf);
766  ideal Js = idInit(k);
767  for (int i=0; i<k; i++)
768    Js->m[i] = p_PermPoly(Jr->m[i],NULL,r,s,identity,NULL,0);
769
770  // this->reduce(Jr,r);
771  // cleanup
772  id_Delete(&inIsAdjusted,sAdjusted);
773  id_Delete(&inJsAdjusted,sAdjusted);
774  rDelete(sAdjusted);
775  id_Delete(&inIr,r);
776  id_Delete(&Jr,r);
777  id_Delete(&inJr,r);
778
779  assume(isValuationTrivial() || isOrderingLocalInT(s));
780  return std::make_pair(Js,s);
781}
782
783
784bool tropicalStrategy::checkForUniformizingBinomial(const ideal I, const ring r) const
785{
786  // if the valuation is trivial,
787  // then there is no special condition the first generator has to fullfill
788  if (isValuationTrivial())
789    return true;
790
791  // if the valuation is non-trivial then checks if the first generator is p-t
792  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
793  poly p = p_One(r);
794  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
795  poly t = p_One(r);
796  p_SetExp(t,1,1,r);
797  p_Setm(t,r);
798  poly pt = p_Add_q(p,p_Neg(t,r),r);
799
800  for (int i=0; i<idSize(I); i++)
801  {
802    if (p_EqualPolys(I->m[i],pt,r))
803    {
804      p_Delete(&pt,r);
805      return true;
806    }
807  }
808  p_Delete(&pt,r);
809  return false;
810}
811
812int tropicalStrategy::findPositionOfUniformizingBinomial(const ideal I, const ring r) const
813{
814  assume(isValuationNonTrivial());
815
816  // if the valuation is non-trivial then checks if the first generator is p-t
817  nMapFunc identity = n_SetMap(startingRing->cf,r->cf);
818  poly p = p_One(r);
819  p_SetCoeff(p,identity(uniformizingParameter,startingRing->cf,r->cf),r);
820  poly t = p_One(r);
821  p_SetExp(t,1,1,r);
822  p_Setm(t,r);
823  poly pt = p_Add_q(p,p_Neg(t,r),r);
824
825  for (int i=0; i<idSize(I); i++)
826  {
827    if (p_EqualPolys(I->m[i],pt,r))
828    {
829      p_Delete(&pt,r);
830      return i;
831    }
832  }
833  p_Delete(&pt,r);
834  return -1;
835}
836
837bool tropicalStrategy::checkForUniformizingParameter(const ideal inI, const ring r) const
838{
839  // if the valuation is trivial,
840  // then there is no special condition the first generator has to fullfill
841  if (isValuationTrivial())
842    return true;
843
844  // if the valuation is non-trivial then checks if the first generator is p
845  if (inI->m[0]==NULL)
846    return false;
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
851  for (int i=0; i<idSize(inI); i++)
852  {
853    if (p_EqualPolys(inI->m[i],p,r))
854    {
855      p_Delete(&p,r);
856      return true;
857    }
858  }
859  p_Delete(&p,r);
860  return false;
861}
862
863#ifndef NDEBUG
864tropicalStrategy::tropicalStrategy():
865  originalRing(NULL),
866  originalIdeal(NULL),
867  expectedDimension(NULL),
868  linealitySpace(gfan::ZCone()), // to come, see below
869  startingRing(NULL),            // to come, see below
870  startingIdeal(NULL),           // to come, see below
871  uniformizingParameter(NULL),   // to come, see below
872  shortcutRing(NULL),            // to come, see below
873  onlyLowerHalfSpace(false),
874  weightAdjustingAlgorithm1(NULL),
875  weightAdjustingAlgorithm2(NULL),
876  extraReductionAlgorithm(NULL)
877{
878}
879
880tropicalStrategy tropicalStrategy::debugStrategy(const ideal startIdeal, number unifParameter, ring startRing)
881{
882  tropicalStrategy debug;
883  debug.originalRing = rCopy(startRing);
884  debug.originalIdeal = id_Copy(startIdeal,startRing);
885  debug.startingRing = rCopy(startRing);
886  debug.startingIdeal = id_Copy(startIdeal,startRing);
887  debug.uniformizingParameter = n_Copy(unifParameter,startRing->cf);
888
889  debug.shortcutRing = rCopy0(startRing);
890  nKillChar(debug.shortcutRing->cf);
891  debug.shortcutRing->cf = nInitChar(n_Zp,(void*)(long)IsPrime(n_Int(unifParameter,startRing->cf)));
892  rComplete(debug.shortcutRing);
893  rTest(debug.shortcutRing);
894
895  debug.onlyLowerHalfSpace = true;
896  debug.weightAdjustingAlgorithm1 = valued_adjustWeightForHomogeneity;
897  debug.weightAdjustingAlgorithm2 = valued_adjustWeightUnderHomogeneity;
898  debug.extraReductionAlgorithm = ppreduceInitially;
899
900  return debug;
901}
902
903BOOLEAN computeWitnessDebug(leftv res, leftv args)
904{
905  leftv u = args;
906  if ((u!=NULL) && (u->Typ()==IDEAL_CMD))
907  {
908    leftv v = u->next;
909    if ((v!=NULL) && (v->Typ()==IDEAL_CMD))
910    {
911      leftv w = v->next;
912      if ((w!=NULL) && (w->Typ()==IDEAL_CMD))
913      {
914        leftv x = w->next;
915        if ((x!=NULL) && (x->Typ()==NUMBER_CMD))
916        {
917          omUpdateInfo();
918          Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
919
920          ideal inJ = (ideal) u->CopyD();
921          ideal inI = (ideal) v->CopyD();
922          ideal I = (ideal) w->CopyD();
923          number p = (number) x->CopyD();
924          tropicalStrategy debug = tropicalStrategy::debugStrategy(I,p,currRing);
925          ideal J = debug.computeWitness(inJ,inI,I,currRing);
926          id_Delete(&inJ,currRing);
927          id_Delete(&inI,currRing);
928          id_Delete(&I,currRing);
929          n_Delete(&p,currRing->cf);
930          res->rtyp = IDEAL_CMD;
931          res->data = (char*) J;
932          return FALSE;
933        }
934      }
935    }
936  }
937  return TRUE;
938}
939
940BOOLEAN computeFlipDebug(leftv res, leftv args)
941{
942  leftv u = args;
943  if ((u!=NULL) && (u->Typ()==IDEAL_CMD))
944  {
945    leftv v = u->next;
946    if ((v!=NULL) && (v->Typ()==NUMBER_CMD))
947    {
948      leftv w = v->next;
949      if ((w!=NULL) && (w->Typ()==BIGINTMAT_CMD))
950      {
951        leftv x = w->next;
952        if ((x!=NULL) && (x->Typ()==BIGINTMAT_CMD))
953        {
954          omUpdateInfo();
955          Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
956
957          ideal I = (ideal) u->CopyD();
958          number p = (number) v->CopyD();
959          bigintmat* interiorPoint0 = (bigintmat*) w->CopyD();
960          bigintmat* facetNormal0 = (bigintmat*) x->CopyD();
961          tropicalStrategy debug = tropicalStrategy::debugStrategy(I,p,currRing);
962
963          gfan::ZVector* interiorPoint = bigintmatToZVector(interiorPoint0);
964          gfan::ZVector* facetNormal = bigintmatToZVector(facetNormal0);
965          std::pair<ideal,ring> Js = debug.computeFlip(I,currRing,*interiorPoint,*facetNormal);
966          ideal J = Js.first;
967          ring s = Js.second;
968
969          id_Delete(&J,s);
970          rDelete(s);
971
972          id_Delete(&I,currRing);
973          n_Delete(&p,currRing->cf);
974          delete interiorPoint0;
975          delete facetNormal0;
976          delete interiorPoint;
977          delete facetNormal;
978
979          res->rtyp = NONE;
980          res->data = NULL;
981          return FALSE;
982        }
983      }
984    }
985  }
986  WerrorS("computeFlipDebug: unexpected parameters");
987  return TRUE;
988}
989#endif
Note: See TracBrowser for help on using the repository browser.