source: git/Singular/dyn_modules/customstd/customstd.cc @ 58399e

fieker-DuValspielwiese
Last change on this file since 58399e was 58399e, checked in by Yue <ren@…>, 8 years ago
chg: slimmed down customstd
  • Property mode set to 100755
File size: 16.2 KB
Line 
1#include <Singular/libsingular.h>
2#include <vector>
3#include <iostream>
4
5// global variable potentially storing output
6ideal idealCache=NULL;
7
8std::vector<int> satstdSaturatingVariables;
9
10// //------------------------------------------------------------------------
11// // example of a routine which changes nothing
12// static BOOLEAN display_sp(kStrategy strat)
13// {
14//   // will be call each time a new s-poly is computed (strat->P)
15//   // the algorithm assures that strat->P.p!=NULL, in currRing
16//   // if strat->P.t_p==NULL: strat->P.p->next is in currRing
17//   // otherwise: strat->P.t_p->next==strat->P.p->next, in strat->tailRing
18//   // must return TRUE, if strat->P is changed, FALSE otherwise
19//   PrintS("a new s-poly found: ");
20//   p_Write(strat->P.p,currRing,strat->tailRing);
21//   return FALSE;
22// }
23// static BOOLEAN std_with_display(leftv res, leftv args)
24// {
25//   if (args!=NULL)
26//   {
27//     if (args->Typ()==IDEAL_CMD)
28//     {
29//       ideal I=(ideal)args->Data();
30//       I=kStd(I,currRing->qideal,testHomog,NULL,NULL,0,0,NULL,display_sp);
31//       idSkipZeroes(I);
32//       res->data=(char*)I;
33//       res->rtyp=IDEAL_CMD;
34//       return FALSE;
35//     }
36//   }
37//   WerrorS("expected: std_with_display(`idea;`)");
38//   return TRUE;
39// }
40
41//------------------------------------------------------------------------
42// routine that simplifies the new element by dividing it with the maximal possible
43// partially saturating the ideal with respect to all variables doing so
44static BOOLEAN sat_vars_sp(kStrategy strat)
45{
46  BOOLEAN b = FALSE; // set b to TRUE, if spoly was changed,
47                     // let it remain FALSE otherwise
48  if (strat->P.t_p==NULL)
49  {
50    poly p=strat->P.p;
51    if (pNext(p)==NULL)
52    {
53      // if a term is contained in the ideal, abort std computation
54      // and store the whole ring in idealCache to be returned
55      while ((strat->Ll >= 0))
56        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
57      idealCache = idInit(1);
58      idealCache->m[0] = p_One(currRing);
59      return FALSE;
60    }
61
62    // iterate over all terms of p and
63    // compute the minimum mm of all exponent vectors
64    int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
65    int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
66    p_GetExpV(p,mm,currRing);
67    bool nonTrivialSaturationToBeDone=false;
68    for (p=pNext(p); p!=NULL; pIter(p))
69    {
70      nonTrivialSaturationToBeDone=false;
71      p_GetExpV(p,m0,currRing);
72      for(int i=0; i<satstdSaturatingVariables.size(); i++)
73      {
74        int li = satstdSaturatingVariables[i];
75        mm[li]=si_min(mm[li],m0[li]);
76        if (mm[li]>0) nonTrivialSaturationToBeDone=true;
77      }
78      // abort if the minimum is zero in each component
79      if (nonTrivialSaturationToBeDone==false) break;
80    }
81    if (nonTrivialSaturationToBeDone==true)
82    {
83      // std::cout << "simplifying!" << std::endl;
84      p=p_Copy(strat->P.p,currRing);
85      strat->P.p=p;
86      while(p!=NULL)
87      {
88        for (int i=1; i<satstdSaturatingVariables.size(); i++)
89        {
90          int li = satstdSaturatingVariables[i];
91          p_SubExp(p,li,mm[li],currRing);
92        }
93        p_Setm(p,currRing);
94        pIter(p);
95      }
96      b = TRUE;
97    }
98    omFree(mm);
99    omFree(m0);
100  }
101  else
102  {
103    poly p=strat->P.t_p;
104    if (pNext(p)==NULL)
105    {
106      // if a term is contained in the ideal, abort std computation
107      // and store the output in idealCache to be returned
108      while ((strat->Ll >= 0))
109        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
110      idealCache = idInit(1);
111      idealCache->m[0] = p_One(currRing);
112      return FALSE;
113    }
114
115    // iterate over all terms of p and
116    // compute the minimum mm of all exponent vectors
117    int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
118    int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
119    p_GetExpV(p,mm,strat->tailRing);
120    bool nonTrivialSaturationToBeDone=false;
121    for (p = pNext(p); p!=NULL; pIter(p))
122    {
123      nonTrivialSaturationToBeDone=false;
124      p_GetExpV(p,m0,strat->tailRing);
125      for(int i=0; i<satstdSaturatingVariables.size(); i++)
126      {
127        int li = satstdSaturatingVariables[i];
128        mm[li]=si_min(mm[li],m0[li]);
129        if (mm[li]>0) nonTrivialSaturationToBeDone = true;
130      }
131      // abort if the minimum is zero in each component
132      if (nonTrivialSaturationToBeDone==false) break;
133    }
134    if (nonTrivialSaturationToBeDone==true)
135    {
136      p=p_Copy(strat->P.t_p,strat->tailRing);
137      strat->P.t_p=p;
138      strat->P.p=NULL;
139      while(p!=NULL)
140      {
141        for(int i=1;i<=rVar(currRing);i++)
142          p_SubExp(p,i,mm[i],strat->tailRing);
143        p_Setm(p,strat->tailRing);
144        pIter(p);
145      }
146      strat->P.GetP();
147      b = TRUE;
148    }
149    omFree(mm);
150    omFree(m0);
151  }
152  return b; // return TRUE if sp was changed, FALSE if not
153}
154
155// returns 1<=i<=rVar(r) if x is the i-th variable,
156// return 0 otherwise
157static int getVariableIndex(const poly x, const ring r)
158{
159  // return 0 if x is not a monomial
160  if ((x->next!=NULL) || (!n_IsOne(p_GetCoeff(x,r),r->cf)))
161    return 0;
162
163  // find the first variable index with non-zero exponent
164  int i=1;
165  int l=0;
166  for (; i<=rVar(r); i++)
167  {
168    l = p_GetExp(x,i,r);
169    if (l>0) break;
170  }
171  // return 0 if no such l exist or l is bigger than one
172  if (l!=1)
173    return 0;
174
175  // check that remaining variables have zero exponent
176  for (i++; i<=rVar(r); i++)
177  {
178    if (p_GetExp(x,i,r)>0)
179      return 0;
180  }
181
182  return l;
183}
184
185//------------------------------------------------------------------------
186// routine that simplifies the ideal dividing each generator by the maximal monomial dividing it
187// in particular returns 1 if a generator is a term
188// to be used before starting saturation with respect to all variables
189static void satSimplify(ideal I, const ring r)
190{
191  idSkipZeroes(I);
192  int k = IDELEMS(I);
193  int *mm=(int*)omAlloc((1+rVar(r))*sizeof(int));
194  int *m0=(int*)omAlloc((1+rVar(r))*sizeof(int));
195  for (int i=0; i<k; i++)
196  {
197    poly p = I->m[i];
198    if (p != NULL)
199    {
200      // check whether p is a term, return 1 if true
201      if (p->next == NULL)
202      {
203        p_Delete(&I->m[0],r);
204        I->m[0] = p_One(r);
205        for (int j=1; j<k; j++)
206        {
207          p_Delete(&I->m[j],r);
208          I->m[j] = NULL;
209        }
210        idSkipZeroes(I);
211        omFree(mm);
212        omFree(m0);
213        return;
214      }
215
216      // check whether p is divisible by a monomial
217      // divide if true
218      p_GetExpV(p,mm,r);
219      bool satNecessary=false;
220      for (; p!=NULL; pIter(p))
221      {
222        satNecessary=false;
223        p_GetExpV(p,m0,r);
224        for(int i=1;i<=rVar(r);i++)
225        {
226          mm[i]=si_min(mm[i],m0[i]);
227          if (mm[i]>0)
228            satNecessary=true;
229        }
230        if (satNecessary==false)
231          break;
232      }
233      if (satNecessary==true)
234      {
235        for (p=I->m[i]; p!=NULL; pIter(p))
236        {
237          for (int i=1; i<=rVar(r); i++)
238            p_SubExp(p,i,mm[i],r);
239          p_Setm(p,r);
240        }
241      }
242    }
243  }
244  omFree(mm);
245  omFree(m0);
246}
247
248static BOOLEAN satstd(leftv res, leftv args)
249{
250  leftv u = args;
251  if ((u!=NULL) && (u->Typ()==IDEAL_CMD))
252  {
253    leftv v = u->next;
254
255    if (v==NULL)
256    {
257      int n = rVar(currRing);
258      satstdSaturatingVariables = std::vector<int>(n);
259      for (int i=0; i<n; i++)
260        satstdSaturatingVariables[i] = i+1;
261    }
262    else
263    {
264      if (v->Typ()==IDEAL_CMD)
265      {
266        ideal J = (ideal) v->Data();
267
268        int k = idSize(J);
269        satstdSaturatingVariables = std::vector<int>(k);
270        for (int i=0; i<k; i++)
271        {
272          poly x = J->m[i];
273          int li = getVariableIndex(x,currRing);
274          if (li>0)
275            satstdSaturatingVariables[i]=li;
276          else
277          {
278            WerrorS("satstd: second argument only ideals generated by variables supported for now");
279            return FALSE;
280          }
281        }
282      }
283      else
284      {
285        WerrorS("satstd: unexpected parameters");
286        return TRUE;
287      }
288    }
289
290    ideal I = (ideal) u->Data();
291    satSimplify(I,currRing);
292
293    idealCache = NULL;
294    I=kStd(I,currRing->qideal,testHomog,NULL,NULL,0,0,NULL,sat_vars_sp);
295    satstdSaturatingVariables = std::vector<int>();
296
297    res->rtyp=IDEAL_CMD;
298    if (idealCache)
299    {
300      id_Delete(&I,currRing);
301      res->data = (char*) idealCache;
302      idealCache = NULL;
303    }
304    else
305    {
306      idSkipZeroes(I);
307      res->data=(char*)I;
308    }
309    return FALSE;
310  }
311  WerrorS("satstd: unexpected parameters");
312  return TRUE;
313}
314
315// static BOOLEAN abortIfMonomial_sp(kStrategy strat)
316// {
317//   BOOLEAN b = FALSE; // set b to TRUE, if spoly was changed,
318//                      // let it remain FALSE otherwise
319//   if (strat->P.t_p==NULL)
320//   {
321//     poly p=strat->P.p;
322//     if (pNext(p)==NULL)
323//     {
324//       // if a term is contained in the ideal, abort std computation
325//       // and store the output in idealCache to be returned
326//       while ((strat->Ll >= 0))
327//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
328//       std::cout << "aborting!" << std::endl;
329//       return FALSE;
330//     }
331//   }
332//   else
333//   {
334//     poly p=strat->P.t_p;
335//     if (pNext(p)==NULL)
336//     {
337//       // if a term is contained in the ideal, abort std computation
338//       // and store the output in idealCache to be returned
339//       while ((strat->Ll >= 0))
340//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
341//       std::cout << "aborting!" << std::endl;
342//       return FALSE;
343//     }
344//   }
345//   return b; // return TRUE if sp was changed, FALSE if not
346// }
347// static BOOLEAN abortifmonomialstd(leftv res, leftv args)
348// {
349//   if (args!=NULL)
350//   {
351//     if ((args->Typ()==IDEAL_CMD) && (args->next==NULL))
352//     {
353//       ideal I=(ideal)args->Data();
354//       idealCache = NULL;
355//       I=kStd(I,currRing->qideal,testHomog,NULL,NULL,0,0,NULL,abortIfMonomial_sp);
356//       res->rtyp=IDEAL_CMD;
357//       if (idealCache)
358//         res->data=(char*)idealCache;
359//       else
360//       {
361//         idSkipZeroes(I);
362//         res->data=(char*)I;
363//       }
364//       return FALSE;
365//     }
366//   }
367//   WerrorS("abortifmonomialstd: unexpected parameters");
368//   return TRUE;
369// }
370
371
372// static long wDeg(const poly p, const ring r)
373// {
374//   if (r->order[0] == ringorder_lp)
375//     return p_GetExp(p,1,currRing);
376//   if (r->order[0] == ringorder_ls)
377//     return -p_GetExp(p,1,currRing);
378
379//   if (r->order[0] == ringorder_dp)
380//   {
381//     long d = 0;
382//     for (int i=1; i<=rVar(r); i++)
383//       d = d + p_GetExp(p,i,r);
384//     return d;
385//   }
386//   if (r->order[0] == ringorder_wp || r->order[0] == ringorder_a)
387//   {
388//     long d = 0;
389//     for (int i=r->block0[0]; i<=r->block1[0]; i++)
390//       d = d + p_GetExp(p,i,r)*r->wvhdl[0][i-1];
391//     return d;
392//   }
393//   if (r->order[0] == ringorder_ws)
394//   {
395//     long d = 0;
396//     for (int i=r->block0[0]; i<=r->block1[0]; i++)
397//       d = d - p_GetExp(p,i,r)*r->wvhdl[0][i-1];
398//     return d;
399//   }
400// }
401
402// static bool isInitialFormMonomial(const poly g, const ring r)
403// {
404//   if (g->next==NULL)
405//     return true;
406//   return wDeg(g,r)!=wDeg(g->next,r);
407// }
408
409// //------------------------------------------------------------------------
410// // routine that checks whether the initial form is a monomial,
411// // breaks computation if it finds one, writing the element into idealCache
412// static BOOLEAN sat_sp_initial(kStrategy strat)
413// {
414//   BOOLEAN b = FALSE; // set b to TRUE, if spoly was changed,
415//                      // let it remain FALSE otherwise
416//   if (strat->P.t_p==NULL)
417//   {
418//     poly p=strat->P.p;
419//     if (pNext(p)==NULL)
420//     {
421//       // if a term is contained in the ideal, abort std computation
422//       // and store the output in idealCache to be returned
423//       while ((strat->Ll >= 0))
424//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
425//       idealCache = idInit(1);
426//       idealCache->m[0] = p_One(currRing);
427//       return FALSE;
428//     }
429//     if (isInitialFormMonomial(p,currRing))
430//     {
431//       while ((strat->Ll >= 0))
432//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
433//       idealCache = idInit(1);
434//       idealCache->m[0] = p_Copy(p,currRing);
435//     }
436//     int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
437//     int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
438//     p_GetExpV(p,mm,currRing);
439//     int m_null=0;
440//     while(p!=NULL)
441//     {
442//       m_null=0;
443//       p_GetExpV(p,m0,currRing);
444//       for(int i=1;i<=rVar(currRing);i++)
445//       {
446//         mm[i]=si_min(mm[i],m0[i]);
447//         if (mm[i]>0) m_null++;
448//       }
449//       if (m_null==0) break;
450//       pIter(p);
451//     }
452//     if (m_null>0)
453//     {
454//       std::cout << "simplifying!" << std::endl;
455//       p=p_Copy(strat->P.p,currRing);
456//       strat->P.p=p;
457//       while(p!=NULL)
458//       {
459//         for(int i=1;i<=rVar(currRing);i++)
460//           p_SubExp(p,i,mm[i],currRing);
461//         p_Setm(p,currRing);
462//         pIter(p);
463//       }
464//       b = TRUE;
465//     }
466//     omFree(mm);
467//     omFree(m0);
468//   }
469//   else
470//   {
471//     poly p=strat->P.t_p;
472//     if (pNext(p)==NULL)
473//     {
474//       // if a term is contained in the ideal, abort std computation
475//       // and store the output in idealCache to be returned
476//       while ((strat->Ll >= 0))
477//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
478//       idealCache = idInit(1);
479//       idealCache->m[0] = p_One(currRing);
480//       return FALSE;
481//     }
482//     if (isInitialFormMonomial(p,strat->tailRing))
483//     {
484//       while ((strat->Ll >= 0))
485//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
486//       nMapFunc identity = n_SetMap(strat->tailRing,currRing);
487//       idealCache = idInit(1);
488//       idealCache->m[0] = p_PermPoly(p,NULL,strat->tailRing,currRing,identity,NULL,0);
489//     }
490//     int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
491//     int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
492//     p_GetExpV(p,mm,strat->tailRing);
493//     int m_null=0;
494//     while(p!=NULL)
495//     {
496//       m_null=0;
497//       p_GetExpV(p,m0,strat->tailRing);
498//       for(int i=1;i<=rVar(currRing);i++)
499//       {
500//         mm[i]=si_min(mm[i],m0[i]);
501//         if (mm[i]>0) m_null++;
502//       }
503//       if (m_null==0) break;
504//       pIter(p);
505//     }
506//     if (m_null>0)
507//     {
508//       std::cout << "simplifying!" << std::endl;
509//       p=p_Copy(strat->P.t_p,strat->tailRing);
510//       strat->P.t_p=p;
511//       strat->P.p=NULL;
512//       while(p!=NULL)
513//       {
514//         for(int i=1;i<=rVar(currRing);i++)
515//           p_SubExp(p,i,mm[i],strat->tailRing);
516//         p_Setm(p,strat->tailRing);
517//         pIter(p);
518//       }
519//       strat->P.GetP();
520//       b = TRUE;
521//     }
522//     omFree(mm);
523//     omFree(m0);
524//   }
525//   return b; // return TRUE if sp was changed, FALSE if not
526// }
527// static BOOLEAN satstdWithInitialCheck(leftv res, leftv args)
528// {
529//   if (args!=NULL)
530//   {
531//     if ((args->Typ()==IDEAL_CMD) && (args->next==NULL))
532//     {
533//       ideal I=(ideal)args->Data();
534//       idealCache = NULL;
535//       I=kStd(I,currRing->qideal,testHomog,NULL,NULL,0,0,NULL,sat_sp_initial);
536//       res->rtyp=IDEAL_CMD;
537//       if (idealCache)
538//         res->data=(char*)idealCache;
539//       else
540//         res->data=(char*)I;
541//       return FALSE;
542//     }
543//   }
544//   WerrorS("satstdWithInitialCheck: unexpected parameters");
545//   return TRUE;
546// }
547
548
549
550//------------------------------------------------------------------------
551// initialisation of the module
552extern "C" int SI_MOD_INIT(std_demo)(SModulFunctions* p)
553{
554  // p->iiAddCproc("std_demo","std_with_display",FALSE,std_with_display);
555  p->iiAddCproc("customstd","satstd",FALSE,satstd);
556  // p->iiAddCproc("std_demo","satstdWithInitialCheck",FALSE,satstdWithInitialCheck);
557  // p->iiAddCproc("std_demo","abortifmonomialstd",FALSE,abortifmonomialstd);
558  // PrintS("init of std_demo - type `listvar(Std_demo);` to its contents\n");
559  return (MAX_TOK);
560}
Note: See TracBrowser for help on using the repository browser.