source: git/Singular/dyn_modules/gfanlib/tropical.cc @ 7a4e77

jengelh-datetimespielwiese
Last change on this file since 7a4e77 was 7a4e77, checked in by Yue Ren <ren@…>, 9 years ago
new: reduceInitially(poly &h, const poly g)
  • Property mode set to 100644
File size: 27.8 KB
Line 
1#include <kernel/polys.h>
2#include <kernel/kstd1.h>
3#include <libpolys/coeffs/longrat.h>
4#include <libpolys/polys/clapsing.h>
5#include <bbcone.h>
6
7#include <initialReduction.h>
8
9poly initial(poly p)
10{
11  poly g = p;
12  poly h = p_Head(g,currRing);
13  poly f = h;
14  long d = p_Deg(g,currRing);
15  pIter(g);
16  while ((g != NULL) && (p_Deg(g,currRing) == d))
17  {
18    pNext(h) = p_Head(g,currRing);
19    pIter(h);
20    pIter(g);
21  }
22  return(f);
23}
24
25
26BOOLEAN initial(leftv res, leftv args)
27{
28  leftv u = args;
29  if ((u != NULL) && (u->Typ() == POLY_CMD))
30  {
31    leftv v = u->next;
32    if (v == NULL)
33    {
34      poly p = (poly) u->Data();
35      res->rtyp = POLY_CMD;
36      res->data = (void*) initial(p);
37      return FALSE;
38    }
39  }
40  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
41  {
42    leftv v = u->next;
43    if (v == NULL)
44    {
45      ideal I = (ideal) u->Data();
46      ideal inI = idInit(IDELEMS(I));
47      for (int i=0; i<IDELEMS(I); i++)
48        inI->m[i]=initial(I->m[i]);
49      res->rtyp = IDEAL_CMD;
50      res->data = (void*) inI;
51      return FALSE;
52    }
53  }
54  WerrorS("initial: unexpected parameters");
55  return TRUE;
56}
57
58
59BOOLEAN homogeneitySpace(leftv res, leftv args)
60{
61  leftv u = args;
62  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
63  {
64    leftv v = u->next;
65    if (v == NULL)
66    {
67      int n = currRing->N;
68      ideal I = (ideal) u->Data();
69      poly g;
70      int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
71      int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
72      gfan::ZVector leadexpw = gfan::ZVector(n);
73      gfan::ZVector tailexpw = gfan::ZVector(n);
74      gfan::ZMatrix equations = gfan::ZMatrix(0,n);
75      for (int i=0; i<IDELEMS(I); i++)
76      {
77        g = (poly) I->m[i]; pGetExpV(g,leadexpv);
78        leadexpw = intStar2ZVector(n, leadexpv);
79        pIter(g);
80        while (g != NULL)
81        {
82          pGetExpV(g,tailexpv);
83          tailexpw = intStar2ZVector(n, tailexpv);
84          equations.appendRow(leadexpw-tailexpw);
85          pIter(g);
86        }
87      }
88      gfan::ZCone* gCone = new gfan::ZCone(gfan::ZMatrix(0, equations.getWidth()),equations);
89      omFreeSize(leadexpv,(n+1)*sizeof(int));
90      omFreeSize(tailexpv,(n+1)*sizeof(int));
91
92      res->rtyp = coneID;
93      res->data = (void*) gCone;
94      return FALSE;
95    }
96  }
97  WerrorS("homogeneitySpace: unexpected parameters");
98  return TRUE;
99}
100
101
102BOOLEAN groebnerCone(leftv res, leftv args)
103{
104  leftv u = args;
105  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
106  {
107    leftv v = u->next;
108    if (v == NULL)
109    {
110      int n = currRing->N;
111      ideal I = (ideal) u->Data();
112      poly g = NULL;
113      int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
114      int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
115      gfan::ZVector leadexpw = gfan::ZVector(n);
116      gfan::ZVector tailexpw = gfan::ZVector(n);
117      gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
118      gfan::ZMatrix equations = gfan::ZMatrix(0,n);
119      long d;
120      for (int i=0; i<IDELEMS(I); i++)
121      {
122        g = (poly) I->m[i]; pGetExpV(g,leadexpv);
123        leadexpw = intStar2ZVector(n, leadexpv);
124        pIter(g);
125        d = p_Deg(g,currRing);
126        while ((g != NULL) && (p_Deg(g,currRing) == d))
127        {
128          pGetExpV(g,tailexpv);
129          tailexpw = intStar2ZVector(n, tailexpv);
130          equations.appendRow(leadexpw-tailexpw);
131          pIter(g);
132        }
133
134        if (g != NULL)
135        {
136          while (g != NULL)
137          {
138            pGetExpV(g,tailexpv);
139            tailexpw = intStar2ZVector(n, tailexpv);
140            inequalities.appendRow(leadexpw-tailexpw);
141            pIter(g);
142          }
143        }
144      }
145      gfan::ZCone* gCone = new gfan::ZCone(inequalities,equations);
146      omFreeSize(leadexpv,(n+1)*sizeof(int));
147      omFreeSize(tailexpv,(n+1)*sizeof(int));
148
149      res->rtyp = coneID;
150      res->data = (void*) gCone;
151      return FALSE;
152    }
153  }
154  WerrorS("groebnerCone: unexpected parameters");
155  return TRUE;
156}
157
158
159gfan::ZCone* maximalGroebnerCone(const ring &r, const ideal &I)
160{
161  int n = rVar(r);
162  poly g = NULL;
163  int* leadexpv = (int*) omAlloc((n+1)*sizeof(int));
164  int* tailexpv = (int*) omAlloc((n+1)*sizeof(int));
165  gfan::ZVector leadexpw = gfan::ZVector(n);
166  gfan::ZVector tailexpw = gfan::ZVector(n);
167  gfan::ZMatrix inequalities = gfan::ZMatrix(0,n);
168  for (int i=0; i<IDELEMS(I); i++)
169  {
170    g = (poly) I->m[i]; pGetExpV(g,leadexpv);
171    leadexpw = intStar2ZVector(n, leadexpv);
172    pIter(g);
173    while (g != NULL)
174    {
175      pGetExpV(g,tailexpv);
176      tailexpw = intStar2ZVector(n, tailexpv);
177      inequalities.appendRow(leadexpw-tailexpw);
178      pIter(g);
179    }
180  }
181  omFreeSize(leadexpv,(n+1)*sizeof(int));
182  omFreeSize(tailexpv,(n+1)*sizeof(int));
183  return new gfan::ZCone(inequalities,gfan::ZMatrix(0, inequalities.getWidth()));
184}
185
186
187BOOLEAN maximalGroebnerCone(leftv res, leftv args)
188{
189  leftv u = args;
190  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
191  {
192    leftv v = u->next;
193    if (v == NULL)
194    {
195      ideal I = (ideal) u->Data();
196      res->rtyp = coneID;
197      res->data = (void*) maximalGroebnerCone(currRing, I);
198      return FALSE;
199    }
200  }
201  WerrorS("maximalGroebnerCone: unexpected parameters");
202  return TRUE;
203}
204
205/***
206 * If 1, replaces all occuring t with prime p,
207 * where theoretically feasable.
208 * Also computes GCD over integers than
209 * over univariate polynomials
210 **/
211#define T_TO_P 0
212
213/***
214 * Suppose I=g_0,...,g_{n-1}.
215 * This function uses bubble sort to sorts g_1,...,g_{n-1}
216 * such that lm(g_1)>...>lm(g_{n-1})
217 **/
218static inline void sortingLaterGenerators(ideal I)
219{
220  poly cache; int i; int n=IDELEMS(I); int newn;
221  do
222  {
223    newn=0;
224    for (i=2; i<n; i++)
225    {
226      if (pLmCmp(I->m[i-1],I->m[i])<0)
227      {
228        cache=I->m[i-1];
229        I->m[i-1]=I->m[i];
230        I->m[i]=cache;
231        newn = i;
232      }
233    }
234    n=newn;
235  } while(n);
236}
237
238
239// /***
240//  * replaces coefficients in g of lowest degree in t
241//  * divisible by p with a suitable power of t
242//  **/
243// static bool preduce(poly g, const number p)
244// {
245//   // go along g and look for relevant terms.
246//   // monomials in x which have been checked are tracked in done.
247//   // because we assume the leading coefficient of g is 1,
248//   // the leading term does not need to be considered.
249//   poly done = p_LmInit(g,currRing);
250//   p_SetExp(done,1,0,currRing); p_SetCoeff(done,n_Init(1,currRing->cf),currRing);
251//   p_Setm(done,currRing);
252//   poly doneEnd = done;
253//   poly doneCache;
254//   number dnumber; long d;
255//   poly subst; number ppower;
256//   while(pNext(g))
257//   {
258//     // check whether next term needs to be reduced:
259//     // first, check whether monomial in x has come up yet
260//     for (doneCache=done; doneCache; pIter(doneCache))
261//     {
262//       if (p_LmDivisibleBy(doneCache,pNext(g),currRing))
263//         break;
264//     }
265//     if (!doneCache) // if for loop did not terminate prematurely,
266//                     // then the monomial in x is new
267//     {
268//       // second, check whether coefficient is divisible by p
269//       if (n_DivBy(p_GetCoeff(pNext(g),currRing->cf),p,currRing->cf))
270//       {
271//         // reduce the term with respect to p-t:
272//         // divide coefficient by p, remove old term,
273//         // add t multiple of old term
274//         dnumber = n_Div(p_GetCoeff(pNext(g),currRing->cf),p,currRing->cf);
275//         d = n_Int(dnumber,currRing->cf);
276//         n_Delete(&dnumber,currRing->cf);
277//         if (!d)
278//         {
279//           p_Delete(&done,currRing);
280//           WerrorS("initialReduction: overflow in a t-exponent");
281//           return true;
282//         }
283//         subst=p_LmInit(pNext(g),currRing);
284//         if (d>0)
285//         {
286//           p_AddExp(subst,1,d,currRing);
287//           p_SetCoeff(subst,n_Init(1,currRing->cf),currRing);
288//         }
289//         else
290//         {
291//           p_AddExp(subst,1,-d,currRing);
292//           p_SetCoeff(subst,n_Init(-1,currRing->cf),currRing);
293//         }
294//         p_Setm(subst,currRing); pTest(subst);
295//         pNext(g)=p_LmDeleteAndNext(pNext(g),currRing);
296//         pNext(g)=p_Add_q(pNext(g),subst,currRing);
297//         pTest(pNext(g));
298//       }
299//       // either way, add monomial in x to done
300//       pNext(doneEnd)=p_LmInit(pNext(g),currRing);
301//       pIter(doneEnd);
302//       p_SetExp(doneEnd,1,0,currRing);
303//       p_SetCoeff(doneEnd,n_Init(1,currRing->cf),currRing);
304//       p_Setm(doneEnd,currRing);
305//     }
306//     pIter(g);
307//   }
308//   p_Delete(&done,currRing);
309//   return false;
310// }
311
312
313// #ifndef NDEBUG
314// BOOLEAN preduce(leftv res, leftv args)
315// {
316//   leftv u = args;
317//   if ((u != NULL) && (u->Typ() == POLY_CMD))
318//   {
319//     poly g; bool b;
320//     number p = n_Init(3,currRing->cf);
321//     omUpdateInfo();
322//     Print("usedBytes=%ld\n",om_Info.UsedBytes);
323//     g = (poly) u->CopyD();
324//     b = preduce(g,p);
325//     p_Delete(&g,currRing);
326//     if (b) return TRUE;
327//     omUpdateInfo();
328//     Print("usedBytes=%ld\n",om_Info.UsedBytes);
329//     n_Delete(&p,currRing->cf);
330//     res->rtyp = NONE;
331//     res->data = NULL;
332//     return FALSE;
333//   }
334//   return TRUE;
335// }
336// #endif //NDEBUG
337
338
339// /***
340//  * Returns the highest term in g with the matching x-monomial to m
341//  * or, if it does not exist, the NULL pointer
342//  **/
343// static poly highestMatchingX(poly g, const poly m)
344// {
345//   pTest(g); pTest(m);
346//   poly xalpha=p_LmInit(m,currRing);
347
348//   // go along g and find the first term
349//   // with the same monomial in x as xalpha
350//   while (g)
351//   {
352//     p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing);
353//     p_Setm(xalpha,currRing);
354//     if (p_LmEqual(g,xalpha,currRing))
355//       break;
356//     pIter(g);
357//   }
358
359//   // gCache now either points at the wanted term
360//   // or is NULL
361//   p_Delete(&xalpha,currRing);
362//   pTest(g);
363//   return g;
364// }
365
366
367// /***
368//  * Given g with lm(g)=t^\beta x^\alpha, returns g_\alpha
369//  **/
370// static poly powerSeriesCoeff(poly g)
371// {
372//   // the first term obviously is part of our output
373//   // so we copy it...
374//   poly xalpha=p_LmInit(g,currRing);
375//   poly coeff=p_One(currRing);
376//   p_SetCoeff(coeff,n_Copy(p_GetCoeff(g,currRing),currRing->cf),currRing);
377//   p_SetExp(coeff,1,p_GetExp(g,1,currRing),currRing);
378//   p_Setm(coeff,currRing);
379
380//   // ..and proceed with the remaining terms,
381//   // appending the relevant terms to coeff via coeffCache
382//   poly coeffCache=coeff;
383//   pIter(g);
384//   while (g)
385//   {
386//     p_SetExp(xalpha,1,p_GetExp(g,1,currRing),currRing);
387//     p_Setm(xalpha,currRing);
388//     if (p_LmEqual(g,xalpha,currRing))
389//     {
390//       pNext(coeffCache)=p_Init(currRing);
391//       pIter(coeffCache);
392//       p_SetCoeff(coeffCache,n_Copy(p_GetCoeff(g,currRing),currRing->cf),currRing);
393//       p_SetExp(coeffCache,1,p_GetExp(g,1,currRing),currRing);
394//       p_Setm(coeffCache,currRing);
395//     }
396//     pIter(g);
397//   }
398
399//   p_Delete(&xalpha,currRing);
400//   return coeff;
401// }
402
403
404// /***
405//  * divides g by t^b knowing that each term of g
406//  * is divisible by t^b, i.e. no divisibility checks
407//  * needed
408//  **/
409// static void divideByT(poly g, const long b)
410// {
411//   while (g)
412//   {
413//     p_SetExp(g,1,p_GetExp(g,1,currRing)-b,currRing);
414//     p_Setm(g,currRing);
415//     pIter(g);
416//   }
417// }
418
419
420// static void divideByGcd(poly &g)
421// {
422//   // first determine all g_\alpha
423//   // as alpha runs over all exponent vectors
424//   // of monomials in x occuring in g
425
426//   // gAlpha will store all g_\alpha,
427//   // the first term will, for comparison purposes for now,
428//   // also keep their monomial in x.
429//   // we assume that the weight on the x are positive
430//   // so that the x won't make the monomial smaller
431//   ideal gAlphaFront = idInit(pLength(g));
432//   gAlphaFront->m[0] = p_Head(g,currRing);
433//   p_SetExp(gAlphaFront->m[0],1,0,currRing);
434//   p_Setm(gAlphaFront->m[0],currRing);
435//   ideal gAlphaEnd = idInit(pLength(g));
436//   gAlphaEnd->m[0] = gAlphaFront->m[0];
437//   unsigned long gAlpha_sev[pLength(g)];
438//   gAlpha_sev[0] = p_GetShortExpVector(g,currRing);
439//   long beta[pLength(g)];
440//   beta[0] = p_GetExp(g,1,currRing);
441//   int count=0;
442
443//   poly current = pNext(g); unsigned long current_not_sev;
444//   int i;
445//   while (current)
446//   {
447//     // see whether the monomial in x of current already came up
448//     // since everything is homogeneous in x and the ordering is local in t,
449//     // this comes down to a divisibility test in two stages
450//     current_not_sev = ~p_GetShortExpVector(current,currRing);
451//     for(i=0; i<=count; i++)
452//     {
453//       // first stage using short exponent vectors
454//       // second stage a proper test
455//       if (p_LmShortDivisibleBy(gAlphaFront->m[i],gAlpha_sev[i],current,current_not_sev, currRing)
456//             && p_LmDivisibleBy(gAlphaFront->m[i],current,currRing))
457//       {
458//         // if it already occured, add the term to the respective entry
459//         // without the x part
460//         pNext(gAlphaEnd->m[i])=p_Init(currRing);
461//         pIter(gAlphaEnd->m[i]);
462//         p_SetExp(gAlphaEnd->m[i],1,p_GetExp(current,1,currRing)-beta[i],currRing);
463//         p_SetCoeff(gAlphaEnd->m[i],n_Copy(p_GetCoeff(current,currRing),currRing->cf),currRing);
464//         p_Setm(gAlphaEnd->m[i],currRing);
465//         break;
466//       }
467//     }
468//     if (i>count)
469//     {
470//       // if it is new, create a new entry for the term
471//       // including the monomial in x
472//       count++;
473//       gAlphaFront->m[count] = p_Head(current,currRing);
474//       beta[count] = p_GetExp(current,1,currRing);
475//       gAlphaEnd->m[count] = gAlphaFront->m[count];
476//       gAlpha_sev[count] = p_GetShortExpVector(current,currRing);
477//     }
478
479//     pIter(current);
480//   }
481
482//   if (count)
483//   {
484//     // now remove all the x in the leading terms
485//     // so that gAlpha only contais terms in t
486//     int j; long d;
487//     for (i=0; i<=count; i++)
488//     {
489//       for (j=2; j<=currRing->N; j++)
490//         p_SetExp(gAlphaFront->m[i],j,0,currRing);
491//       p_Setm(gAlphaFront->m[i],currRing);
492//       gAlphaEnd->m[i]=NULL;
493//     }
494
495//     // now compute the gcd over all g_\alpha
496//     // and set the input to null as they are deleted in the process
497//     poly gcd = singclap_gcd(gAlphaFront->m[0],gAlphaFront->m[1],currRing);
498//     gAlphaFront->m[0] = NULL;
499//     gAlphaFront->m[1] = NULL;
500//     for(i=2; i<=count; i++)
501//     {
502//       gcd = singclap_gcd(gcd,gAlphaFront->m[i],currRing);
503//       gAlphaFront->m[i] = NULL;
504//     }
505//     // divide g by the gcd
506//     poly h = singclap_pdivide(g,gcd,currRing);
507//     p_Delete(&gcd,currRing);
508//     p_Delete(&g,currRing);
509//     g = h;
510
511//     id_Delete(&gAlphaFront,currRing);
512//     id_Delete(&gAlphaEnd,currRing);
513//   }
514//   else
515//   {
516//     // if g contains only one monomial in x,
517//     // then we can get rid of all the higher t
518//     p_Delete(&g,currRing);
519//     g = gAlphaFront->m[0];
520//     pIter(gAlphaFront->m[0]);
521//     pNext(g)=NULL;
522//     gAlphaEnd->m[0] = NULL;
523//     id_Delete(&gAlphaFront,currRing);
524//     id_Delete(&gAlphaEnd,currRing);
525//   }
526// }
527
528
529// /***
530//  * 1. For each \alpha\in\NN^n, changes (c_\beta t^\beta + c_{\beta+1} t^{\beta+1} + ...)
531//  *    to (c_\beta + c_{\beta+1}*p + ...) t^\beta
532//  * 2. Computes the gcd over all (c_\beta + c_{\beta+1}*p + ...) and divides g by it
533//  **/
534// static void simplify(poly g, const number p)
535// {
536//   // go along g and look for relevant terms.
537//   // monomials in x which have been checked are tracked in done.
538//   poly done = p_LmInit(g,currRing);
539//   p_SetCoeff(done,n_Init(1,currRing->cf),currRing);
540//   p_Setm(done,currRing);
541//   poly doneEnd = done;
542//   poly doneCurrent;
543
544//   poly subst; number substCoeff, substCoeffCache;
545//   unsigned long d;
546
547//   poly gCurrent = g;
548//   while(pNext(gCurrent))
549//   {
550//     // check whether next term needs to be reduced:
551//     // first, check whether monomial in x has come up yet
552//     for (doneCurrent=done; doneCurrent; pIter(doneCurrent))
553//     {
554//       if (p_LmDivisibleBy(doneCurrent,pNext(gCurrent),currRing))
555//         break;
556//     }
557//     // if the monomial in x already occured, then we want to replace
558//     // as many t with p as theoretically feasable
559//     if (doneCurrent)
560//     {
561//       // suppose pNext(gCurrent)=3*t5x and doneCurrent=t3x
562//       // then we want to replace pNext(gCurrent) with 3p2*t3x
563//       // subst = ?*t3x
564//       subst = p_LmInit(doneCurrent,currRing);
565//       // substcoeff = p2
566//       n_Power(p,p_GetExp(subst,1,currRing)-p_GetExp(doneCurrent,1,currRing),&substCoeff,currRing->cf);
567//       // substcoeff = 3p2
568//       n_InpMult(substCoeff,p_GetCoeff(pNext(gCurrent),currRing),currRing->cf);
569//       // subst = 3p2*t3x
570//       p_SetCoeff(subst,substCoeff,currRing);
571//       p_Setm(subst,currRing); pTest(subst);
572
573//       // g = g - pNext(gCurrent) + subst
574//       pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing);
575//       g=p_Add_q(g,subst,currRing);
576//       pTest(pNext(gbeginning));
577//     }
578//     else
579//     {
580//       // if the monomial in x is brand new,
581//       // then we check whether the coefficient is divisible by p
582//       if (n_DivBy(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf))
583//       {
584//         // reduce the term with respect to p-t:
585//         // suppose pNext(gCurrent)=4p3*tx
586//         // then we want to replace it with 4*t4x
587//         // divide 4p3 repeatedly by p until it is not divisible anymore,
588//         // keeping track on the final value 4
589//         // and the number of times it has been divided 3
590//         substCoeff = n_Div(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf);
591//         d = 1;
592//         while (n_DivBy(substCoeff,p,currRing->cf))
593//         {
594//           substCoeffCache = n_Div(p_GetCoeff(pNext(gCurrent),currRing->cf),p,currRing->cf);
595//           n_Delete(&substCoeff,currRing->cf);
596//           substCoeff = substCoeffCache;
597//           d++;
598//           assume(d>0); // check for overflow
599//         }
600
601//         // subst = ?*tx
602//         subst=p_LmInit(pNext(gCurrent),currRing);
603//         // subst = ?*t4x
604//         p_AddExp(subst,1,d,currRing);
605//         // subst = 4*t4x
606//         p_SetCoeff(subst,substCoeffCache,currRing);
607//         p_Setm(subst,currRing); pTest(subst);
608
609//         // g = g - pNext(gCurrent) + subst
610//         pNext(gCurrent)=p_LmDeleteAndNext(pNext(gCurrent),currRing);
611//         pNext(gCurrent)=p_Add_q(pNext(gCurrent),subst,currRing);
612//         pTest(pNext(gCurrent));
613//       }
614
615//       // either way, add monomial in x with minimal t to done
616//       pNext(doneEnd)=p_LmInit(pNext(gCurrent),currRing);
617//       pIter(doneEnd);
618//       p_SetCoeff(doneEnd,n_Init(1,currRing->cf),currRing);
619//       p_Setm(doneEnd,currRing);
620//     }
621//     pIter(gCurrent);
622//   }
623//   p_Delete(&done,currRing);
624// }
625
626
627// #ifndef NDEBUG
628// BOOLEAN divideByGcd(leftv res, leftv args)
629// {
630//   leftv u = args;
631//   if ((u != NULL) && (u->Typ() == POLY_CMD))
632//   {
633//     poly g;
634//     omUpdateInfo();
635//     Print("usedBytes1=%ld\n",om_Info.UsedBytes);
636//     g = (poly) u->CopyD();
637//     divideByGcd(g);
638//     p_Delete(&g,currRing);
639//     omUpdateInfo();
640//     Print("usedBytes1=%ld\n",om_Info.UsedBytes);
641//     res->rtyp = NONE;
642//     res->data = NULL;
643//     return FALSE;
644//   }
645//   return TRUE;
646// }
647// #endif //NDEBUG
648
649
650// BOOLEAN initialReduction(leftv res, leftv args)
651// {
652//   leftv u = args;
653//   if ((u != NULL) && (u->Typ() == IDEAL_CMD))
654//   {
655//     leftv v = u->next;
656//     if (v == NULL)
657//     {
658//       ideal I = (ideal) u->Data();
659
660//       /***
661//        * Suppose I=g_0,...,g_n and g_0=p-t.
662//        * Step 1: sort elements g_1,...,g_{n-1} such that lm(g_1)>...>lm(g_{n-1})
663//        * (the following algorithm is a bubble sort)
664//        **/
665//       sortingLaterGenerators(I);
666
667//       /***
668//        * Step 2: replace coefficient of terms of lowest t-degree divisible by p with t
669//        **/
670//       number p = p_GetCoeff(I->m[0],currRing);
671//       for (int i=1; i<IDELEMS(I); i++)
672//       {
673//         if (preduce(I->m[i],p))
674//           return TRUE;
675//       }
676
677//       /***
678//        * Step 3: the first pass. removing terms with the same monomials in x as lt(g_i)
679//        * out of g_j for i<j
680//        **/
681//       int i,j;
682//       poly uniti, unitj;
683//       for (i=1; i<IDELEMS(I)-1; i++)
684//       {
685//         for (j=i+1; j<IDELEMS(I); j++)
686//         {
687//           unitj = highestMatchingX(I->m[j],I->m[i]);
688//           if (unitj)
689//           {
690//             unitj = powerSeriesCoeff(unitj);
691//             divideByT(unitj,p_GetExp(I->m[i],1,currRing));
692//             uniti = powerSeriesCoeff(I->m[i]);
693//             divideByT(uniti,p_GetExp(I->m[i],1,currRing));
694//             pTest(unitj); pTest(uniti); pTest(I->m[j]); pTest(I->m[i]);
695//             I->m[j]=p_Add_q(p_Mult_q(uniti,I->m[j],currRing),
696//                             p_Neg(p_Mult_q(unitj,p_Copy(I->m[i],currRing),currRing),currRing),
697//                             currRing);
698//             divideByGcd(I->m[j]);
699//           }
700//         }
701//       }
702//       for (int i=1; i<IDELEMS(I); i++)
703//       {
704//         if (preduce(I->m[i],p))
705//           return TRUE;
706//       }
707
708//       /***
709//        * Step 4: the second pass. removing terms divisible by lt(g_j) out of g_i for i<j
710//        **/
711//       for (i=1; i<IDELEMS(I)-1; i++)
712//       {
713//         for (j=i+1; j<IDELEMS(I); j++)
714//         {
715//           uniti = highestMatchingX(I->m[i],I->m[j]);
716//           if (uniti && p_GetExp(uniti,1,currRing)>=p_GetExp(I->m[j],1,currRing))
717//           {
718//             uniti = powerSeriesCoeff(uniti);
719//             divideByT(uniti,p_GetExp(I->m[j],1,currRing));
720//             unitj = powerSeriesCoeff(I->m[j]);
721//             divideByT(unitj,p_GetExp(I->m[j],1,currRing));
722//             I->m[i] = p_Add_q(p_Mult_q(unitj,I->m[i],currRing),
723//                               p_Neg(p_Mult_q(uniti,p_Copy(I->m[j],currRing),currRing),currRing),
724//                               currRing);
725//             divideByGcd(I->m[j]);
726//           }
727//         }
728//       }
729//       for (int i=1; i<IDELEMS(I); i++)
730//       {
731//         if (preduce(I->m[i],p))
732//           return TRUE;
733//       }
734
735//       res->rtyp = NONE;
736//       res->data = NULL;
737//       IDDATA((idhdl)u->data) = (char*) I;
738//       return FALSE;
739//     }
740//   }
741//   WerrorS("initialReduction: unexpected parameters");
742//   return TRUE;
743// }
744
745
746#if 0
747/***
748 * Given a general ring r with any ordering,
749 * changes the ordering to a(v),ws(-w)
750 **/
751bool changetoAWSRing(ring r, gfan::ZVector v, gfan::ZVector w)
752{
753  omFree(r->order);
754  r->order  = (int*) omAlloc0(4*sizeof(int));
755  omFree(r->block0);
756  r->block0 = (int*) omAlloc0(4*sizeof(int));
757  omFree(r->block1);
758  r->block1 = (int*) omAlloc0(4*sizeof(int));
759  for (int i=0; r->wvhdl[i]; i++)
760  { omFree(r->wvhdl[i]); }
761  omFree(r->wvhdl);
762  r->wvhdl  = (int**) omAlloc0(4*sizeof(int*));
763
764  bool ok = false;
765  r->order[0]  = ringorder_a;
766  r->block0[0] = 1;
767  r->block1[0] = r->N;
768  r->wvhdl[0]  = ZVectorToIntStar(v,ok);
769  r->order[1]  = ringorder_ws;
770  r->block0[1] = 1;
771  r->block1[1] = r->N;
772  r->wvhdl[1]  = ZVectorToIntStar(w,ok);
773  r->order[2]=ringorder_C;
774  return ok;
775}
776
777
778/***
779 * Given a ring with ordering a(v'),ws(w'),
780 * changes the weights to v,w
781 **/
782bool changeAWSWeights(ring r, gfan::ZVector v, gfan::ZVector w)
783{
784  omFree(r->wvhdl[0]);
785  omFree(r->wvhdl[1]);
786  bool ok = false;
787  r->wvhdl[0]  = ZVectorToIntStar(v,ok);
788  r->wvhdl[1]  = ZVectorToIntStar(w,ok);
789  return ok;
790}
791
792
793// /***
794//  * Creates an int* representing the transposition of the last two variables
795//  **/
796// static inline int* createPermutationVectorForSaturation(static const ring &r)
797// {
798//   int* w = (int*) omAlloc0((rVar(r)+1)*sizeof(int));
799//   for (int i=1; i<=rVar(r)-2; i++)
800//     w[i] = i;
801//   w[rVar(r)-1] = rVar(r);
802//   w[rVar(r)] = rVar(r)-1;
803// }
804
805
806/***
807 * Creates an int* representing the permutation
808 * 1 -> 1, ..., i-1 -> i-1, i -> n, i+1 -> n-1, ... , n -> i
809 **/
810static inline int* createPermutationVectorForSaturation(const ring &r, const int i)
811{
812  int* sigma = (int*) omAlloc0((rVar(r)+1)*sizeof(int));
813  int j;
814  for (j=1; j<i; j++)
815    sigma[j] = j;
816  for (; j<=rVar(r); j++)
817    sigma[j] = rVar(r)-j+i;
818  return(sigma);
819}
820
821
822/***
823 * Changes the int* representing the permutation
824 * 1 -> 1, ..., i -> i, i+1 -> n, i+2 -> n-1, ... , n -> i+1
825 * to an int* representing the permutation
826 * 1 -> 1, ..., i-1 -> i-1, i -> n, i+1 -> n-1, ... , n -> i
827 **/
828static void changePermutationVectorForSaturation(int* sigma, const ring &r, const int i)
829{
830  for (int j=i; j<rVar(r); j++)
831    sigma[j] = rVar(r)-j+i;
832  sigma[rVar(r)] = i;
833}
834
835
836/***
837 * returns a ring in which the weights of the ring variables are permuted
838 * if handed over a poly in which the variables are permuted, this is basically
839 * as good as permuting the variables of the ring itself.
840 **/
841static ring permuteWeighstOfRingVariables(const ring &r, const int* const sigma)
842{
843  ring s = rCopy0(r);
844  for (int j=0; j<rVar(r); j++)
845  {
846    s->wvhdl[0][j] = r->wvhdl[0][sigma[j+1]];
847    s->wvhdl[1][j] = r->wvhdl[1][sigma[j+1]];
848  }
849  rComplete(s,1);
850  return s;
851}
852
853
854/***
855 * creates a ring s that is a copy of r except with ordering wp(w)
856 **/
857static inline ring createInitialRingForSaturation(const ring &r, const gfan::ZVector &w, bool &ok)
858{
859  assume(rVar(r) == (int) w.size());
860
861  ring s = rCopy0(r); int i;
862  for (i=0; s->order[i]; i++)
863    omFreeSize(s->wvhdl[i],rVar(r)*sizeof(int));
864  i++;
865  omFreeSize(s->order,i*sizeof(int));
866  s->order  = (int*) omAlloc0(3*sizeof(int));
867  omFreeSize(s->block0,i*sizeof(int));
868  s->block0 = (int*) omAlloc0(3*sizeof(int));
869  omFreeSize(s->block1,i*sizeof(int));
870  s->block1 = (int*) omAlloc0(3*sizeof(int));
871  omFreeSize(s->wvhdl,i*sizeof(int*));
872  s->wvhdl  = (int**) omAlloc0(3*sizeof(int*));
873
874  s->order[0]  = ringorder_wp;
875  s->block0[0] = 1;
876  s->block1[0] = rVar(r);
877  s->wvhdl[0]  = ZVectorToIntStar(w,ok);
878  s->order[1]=ringorder_C;
879
880  rComplete(s,1);
881  return s;
882}
883
884
885/***
886 * Given an weighted homogeneous ideal I with respect to weight w
887 * that in standard basis form with respect to the ordering ws(-w),
888 * derives the standard basis of I:<x_n>^\infty
889 * and returns a long k such that I:<x_n>^\infty=I:<x_n>^k
890 **/
891static long deriveStandardBasisOfSaturation(ideal &I, ring &r)
892{
893  long k=0, l; poly current;
894  for (int i=0; i<IDELEMS(I); i++)
895  {
896    current = I->m[i];
897    l = p_GetExp(current,rVar(r),r);
898    if (k<l) k=l;
899    while (current)
900    {
901      p_SubExp(current,rVar(r),l,r); p_Setm(current,r);
902      pIter(current);
903    }
904  }
905  return k;
906}
907
908
909/***
910 * Given a weighted homogeneous ideal I with respect to weight w
911 * with constant first element,
912 * returns NULL if I does not contain a monomial
913 * otherwise returns the monomial contained in I
914 **/
915poly containsMonomial(const ideal &I, const gfan::ZVector &w)
916{
917  assume(rField_is_Ring_Z(currRing));
918
919  // first we switch to the ground field currRing->cf / I->m[0]
920  ring r = rCopy0(currRing);
921  nKillChar(r->cf);
922  r->cf = nInitChar(n_Zp,(void*)(long)n_Int(p_GetCoeff(I->m[0],currRing),currRing->cf));
923  rComplete(r);
924
925  ideal J = id_Copy(I, currRing); poly cache; number temp;
926  for (int i=0; i<IDELEMS(I); i++)
927  {
928    cache = J->m[i];
929    while (cache)
930    {
931      // TODO: temp = npMapGMP(p_GetCoeff(cache,currRing),currRing->cf,r->cf);
932      p_SetCoeff(cache,temp,r); pIter(cache);
933    }
934  }
935
936
937  J = kStd(J,NULL,isHomog,NULL);
938
939  bool b = false;
940  ring s = createInitialRingForSaturation(currRing, w, b);
941  if (b)
942  {
943    WerrorS("containsMonomial: overflow in weight vector");
944    return NULL;
945  }
946
947  return NULL;
948}
949
950
951gfan::ZCone* startingCone(ideal I)
952{
953  I = kStd(I,NULL,isNotHomog,NULL);
954  gfan::ZCone* zc = maximalGroebnerCone(currRing,I);
955  gfan::ZMatrix rays = zc->extremeRays();
956  gfan::ZVector v;
957  for (int i=0; i<rays.getHeight(); i++)
958  {
959    v = rays[i];
960  }
961  return zc;
962}
963#endif
964
965
966void tropical_setup(SModulFunctions* p)
967{
968  p->iiAddCproc("","groebnerCone",FALSE,groebnerCone);
969  p->iiAddCproc("","maximalGroebnerCone",FALSE,maximalGroebnerCone);
970  p->iiAddCproc("","initial",FALSE,initial);
971#ifndef NDEBUG
972  // p->iiAddCproc("","divideByGcd",FALSE,divideByGcd);
973  p->iiAddCproc("","pReduce",FALSE,pReduce);
974  p->iiAddCproc("","reduceInitially0",FALSE,reduceInitially0);
975  p->iiAddCproc("","reduceInitially1",FALSE,reduceInitially1);
976#endif //NDEBUG
977  // p->iiAddCproc("","initialReduction",FALSE,initialReduction);
978  p->iiAddCproc("","homogeneitySpace",FALSE,homogeneitySpace);
979}
Note: See TracBrowser for help on using the repository browser.