source: git/Singular/dyn_modules/gfanlib/tropicalStrategy.cc @ 5d0cbd

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