source: git/Singular/dyn_modules/customstd/customstd.cc @ dbb5c7e

spielwiese
Last change on this file since dbb5c7e was dbb5c7e, checked in by Yue <ren@…>, 8 years ago
chg: made monomialabortstd avaialble in the interpreter
  • Property mode set to 100644
File size: 12.4 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
44BOOLEAN 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
52    // iterate over all terms of p and
53    // compute the minimum mm of all exponent vectors
54    int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
55    int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
56    p_GetExpV(p,mm,currRing);
57    bool nonTrivialSaturationToBeDone=true;
58    for (p=pNext(p); p!=NULL; pIter(p))
59    {
60      nonTrivialSaturationToBeDone=false;
61      p_GetExpV(p,m0,currRing);
62      for (int i=satstdSaturatingVariables.size()-1; i>=0; i--)
63      {
64        int li = satstdSaturatingVariables[i];
65        mm[li]=si_min(mm[li],m0[li]);
66        if (mm[li]>0) nonTrivialSaturationToBeDone=true;
67      }
68      // abort if the minimum is zero in each component
69      if (nonTrivialSaturationToBeDone==false) break;
70    }
71    if (nonTrivialSaturationToBeDone==true)
72    {
73      // std::cout << "simplifying!" << std::endl;
74      p=p_Copy(strat->P.p,currRing);
75      memset(&strat->P,0,sizeof(strat->P));
76      strat->P.tailRing = strat->tailRing;
77      strat->P.p=p;
78      while(p!=NULL)
79      {
80        for (int i=satstdSaturatingVariables.size()-1; i>=0; i--)
81        {
82          int li = satstdSaturatingVariables[i];
83          p_SubExp(p,li,mm[li],currRing);
84        }
85        p_Setm(p,currRing);
86        pIter(p);
87      }
88      b = TRUE;
89    }
90    omFree(mm);
91    omFree(m0);
92  }
93  else
94  {
95    poly p=strat->P.t_p;
96
97    // iterate over all terms of p and
98    // compute the minimum mm of all exponent vectors
99    int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
100    int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
101    p_GetExpV(p,mm,strat->tailRing);
102    bool nonTrivialSaturationToBeDone=true;
103    for (p = pNext(p); p!=NULL; pIter(p))
104    {
105      nonTrivialSaturationToBeDone=false;
106      p_GetExpV(p,m0,strat->tailRing);
107      for(int i=satstdSaturatingVariables.size()-1; i>=0; i--)
108      {
109        int li = satstdSaturatingVariables[i];
110        mm[li]=si_min(mm[li],m0[li]);
111        if (mm[li]>0) nonTrivialSaturationToBeDone = true;
112      }
113      // abort if the minimum is zero in each component
114      if (!nonTrivialSaturationToBeDone) break;
115    }
116    if (nonTrivialSaturationToBeDone)
117    {
118      p=p_Copy(strat->P.t_p,strat->tailRing);
119      memset(&strat->P,0,sizeof(strat->P));
120      strat->P.tailRing = strat->tailRing;
121      strat->P.t_p=p;
122      while(p!=NULL)
123      {
124        for(int i=satstdSaturatingVariables.size()-1; i>=0; i--)
125        {
126          int li = satstdSaturatingVariables[i];
127          p_SubExp(p,li,mm[li],strat->tailRing);
128        }
129        p_Setm(p,strat->tailRing);
130        pIter(p);
131      }
132      strat->P.GetP();
133      b = TRUE;
134    }
135    omFree(mm);
136    omFree(m0);
137  }
138  return b; // return TRUE if sp was changed, FALSE if not
139}
140
141static BOOLEAN satstd(leftv res, leftv args)
142{
143  leftv u = args;
144  if ((u!=NULL) && (u->Typ()==IDEAL_CMD))
145  {
146    leftv v = u->next;
147
148    if (v==NULL)
149    {
150      int n = rVar(currRing);
151      satstdSaturatingVariables = std::vector<int>(n);
152      for (int i=n-1; i>=0; i--)
153        satstdSaturatingVariables[i] = i+1;
154    }
155    else
156    {
157      if (v->Typ()==IDEAL_CMD)
158      {
159        ideal J = (ideal) v->Data();
160
161        int k = IDELEMS(J);
162        satstdSaturatingVariables = std::vector<int>(k);
163        for (int i=0; i<k; i++)
164        {
165          poly x = J->m[i];
166          int li = p_Var(x,currRing);
167          if (li>0)
168            satstdSaturatingVariables[i]=li;
169          else
170          {
171            WerrorS("satstd: second argument only ideals generated by variables supported for now");
172            return TRUE;
173          }
174        }
175      }
176      else
177      {
178        WerrorS("satstd: unexpected parameters");
179        return TRUE;
180      }
181    }
182
183    ideal I = (ideal) u->Data();
184
185    idealCache = NULL;
186    I=kStd(I,currRing->qideal,testHomog,NULL,NULL,0,0,NULL,sat_vars_sp);
187    satstdSaturatingVariables = std::vector<int>();
188
189    res->rtyp=IDEAL_CMD;
190    if (idealCache)
191    {
192      id_Delete(&I,currRing);
193      res->data = (char*) idealCache;
194      idealCache = NULL;
195    }
196    else
197    {
198      idSkipZeroes(I);
199      res->data=(char*)I;
200    }
201    return FALSE;
202  }
203  WerrorS("satstd: unexpected parameters");
204  return TRUE;
205}
206
207BOOLEAN abort_if_monomial_sp(kStrategy strat)
208{
209  BOOLEAN b = FALSE; // set b to TRUE, if spoly was changed,
210                     // let it remain FALSE otherwise
211  if (strat->P.t_p==NULL)
212  {
213    poly p=strat->P.p;
214    if (pNext(p)==NULL)
215    {
216      while ((strat->Ll >= 0))
217        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
218      return FALSE;
219    }
220  }
221  else
222  {
223    poly p=strat->P.t_p;
224    if (pNext(p)==NULL)
225    {
226      while ((strat->Ll >= 0))
227        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
228      return FALSE;
229    }
230  }
231  return b; // return TRUE if sp was changed, FALSE if not
232}
233
234BOOLEAN monomialabortstd(leftv res, leftv args)
235{
236  if (args!=NULL)
237  {
238    if ((args->Typ()==IDEAL_CMD) && (args->next==NULL))
239    {
240      ideal I=(ideal)args->Data();
241      I=kStd(I,currRing->qideal,testHomog,NULL,NULL,0,0,NULL,abort_if_monomial_sp);
242      idSkipZeroes(I);
243      res->rtyp=IDEAL_CMD;
244      res->data=(char*)I;
245      return FALSE;
246    }
247  }
248  WerrorS("monomialabortstd: unexpected parameters");
249  return TRUE;
250}
251
252// static long wDeg(const poly p, const ring r)
253// {
254//   if (r->order[0] == ringorder_lp)
255//     return p_GetExp(p,1,currRing);
256//   if (r->order[0] == ringorder_ls)
257//     return -p_GetExp(p,1,currRing);
258
259//   if (r->order[0] == ringorder_dp)
260//   {
261//     long d = 0;
262//     for (int i=1; i<=rVar(r); i++)
263//       d = d + p_GetExp(p,i,r);
264//     return d;
265//   }
266//   if (r->order[0] == ringorder_wp || r->order[0] == ringorder_a)
267//   {
268//     long d = 0;
269//     for (int i=r->block0[0]; i<=r->block1[0]; i++)
270//       d = d + p_GetExp(p,i,r)*r->wvhdl[0][i-1];
271//     return d;
272//   }
273//   if (r->order[0] == ringorder_ws)
274//   {
275//     long d = 0;
276//     for (int i=r->block0[0]; i<=r->block1[0]; i++)
277//       d = d - p_GetExp(p,i,r)*r->wvhdl[0][i-1];
278//     return d;
279//   }
280// }
281
282// static bool isInitialFormMonomial(const poly g, const ring r)
283// {
284//   if (g->next==NULL)
285//     return true;
286//   return wDeg(g,r)!=wDeg(g->next,r);
287// }
288
289// //------------------------------------------------------------------------
290// // routine that checks whether the initial form is a monomial,
291// // breaks computation if it finds one, writing the element into idealCache
292// static BOOLEAN sat_sp_initial(kStrategy strat)
293// {
294//   BOOLEAN b = FALSE; // set b to TRUE, if spoly was changed,
295//                      // let it remain FALSE otherwise
296//   if (strat->P.t_p==NULL)
297//   {
298//     poly p=strat->P.p;
299//     if (pNext(p)==NULL)
300//     {
301//       // if a term is contained in the ideal, abort std computation
302//       // and store the output in idealCache to be returned
303//       while ((strat->Ll >= 0))
304//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
305//       idealCache = idInit(1);
306//       idealCache->m[0] = p_One(currRing);
307//       return FALSE;
308//     }
309//     if (isInitialFormMonomial(p,currRing))
310//     {
311//       while ((strat->Ll >= 0))
312//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
313//       idealCache = idInit(1);
314//       idealCache->m[0] = p_Copy(p,currRing);
315//     }
316//     int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
317//     int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
318//     p_GetExpV(p,mm,currRing);
319//     int m_null=0;
320//     while(p!=NULL)
321//     {
322//       m_null=0;
323//       p_GetExpV(p,m0,currRing);
324//       for(int i=1;i<=rVar(currRing);i++)
325//       {
326//         mm[i]=si_min(mm[i],m0[i]);
327//         if (mm[i]>0) m_null++;
328//       }
329//       if (m_null==0) break;
330//       pIter(p);
331//     }
332//     if (m_null>0)
333//     {
334//       std::cout << "simplifying!" << std::endl;
335//       p=p_Copy(strat->P.p,currRing);
336//       strat->P.p=p;
337//       while(p!=NULL)
338//       {
339//         for(int i=1;i<=rVar(currRing);i++)
340//           p_SubExp(p,i,mm[i],currRing);
341//         p_Setm(p,currRing);
342//         pIter(p);
343//       }
344//       b = TRUE;
345//     }
346//     omFree(mm);
347//     omFree(m0);
348//   }
349//   else
350//   {
351//     poly p=strat->P.t_p;
352//     if (pNext(p)==NULL)
353//     {
354//       // if a term is contained in the ideal, abort std computation
355//       // and store the output in idealCache to be returned
356//       while ((strat->Ll >= 0))
357//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
358//       idealCache = idInit(1);
359//       idealCache->m[0] = p_One(currRing);
360//       return FALSE;
361//     }
362//     if (isInitialFormMonomial(p,strat->tailRing))
363//     {
364//       while ((strat->Ll >= 0))
365//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
366//       nMapFunc identity = n_SetMap(strat->tailRing,currRing);
367//       idealCache = idInit(1);
368//       idealCache->m[0] = p_PermPoly(p,NULL,strat->tailRing,currRing,identity,NULL,0);
369//     }
370//     int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
371//     int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
372//     p_GetExpV(p,mm,strat->tailRing);
373//     int m_null=0;
374//     while(p!=NULL)
375//     {
376//       m_null=0;
377//       p_GetExpV(p,m0,strat->tailRing);
378//       for(int i=1;i<=rVar(currRing);i++)
379//       {
380//         mm[i]=si_min(mm[i],m0[i]);
381//         if (mm[i]>0) m_null++;
382//       }
383//       if (m_null==0) break;
384//       pIter(p);
385//     }
386//     if (m_null>0)
387//     {
388//       std::cout << "simplifying!" << std::endl;
389//       p=p_Copy(strat->P.t_p,strat->tailRing);
390//       strat->P.t_p=p;
391//       strat->P.p=NULL;
392//       while(p!=NULL)
393//       {
394//         for(int i=1;i<=rVar(currRing);i++)
395//           p_SubExp(p,i,mm[i],strat->tailRing);
396//         p_Setm(p,strat->tailRing);
397//         pIter(p);
398//       }
399//       strat->P.GetP();
400//       b = TRUE;
401//     }
402//     omFree(mm);
403//     omFree(m0);
404//   }
405//   return b; // return TRUE if sp was changed, FALSE if not
406// }
407// static BOOLEAN satstdWithInitialCheck(leftv res, leftv args)
408// {
409//   if (args!=NULL)
410//   {
411//     if ((args->Typ()==IDEAL_CMD) && (args->next==NULL))
412//     {
413//       ideal I=(ideal)args->Data();
414//       idealCache = NULL;
415//       I=kStd(I,currRing->qideal,testHomog,NULL,NULL,0,0,NULL,sat_sp_initial);
416//       res->rtyp=IDEAL_CMD;
417//       if (idealCache)
418//         res->data=(char*)idealCache;
419//       else
420//         res->data=(char*)I;
421//       return FALSE;
422//     }
423//   }
424//   WerrorS("satstdWithInitialCheck: unexpected parameters");
425//   return TRUE;
426// }
427
428
429
430//------------------------------------------------------------------------
431// initialisation of the module
432extern "C" int SI_MOD_INIT(customstd)(SModulFunctions* p)
433{
434  // p->iiAddCproc("std_demo","std_with_display",FALSE,std_with_display);
435  p->iiAddCproc("customstd.so","satstd",FALSE,satstd);
436  // p->iiAddCproc("std_demo","satstdWithInitialCheck",FALSE,satstdWithInitialCheck);
437  p->iiAddCproc("std_demo","monomialabortstd",FALSE,monomialabortstd);
438  // PrintS("init of std_demo - type `listvar(Std_demo);` to its contents\n");
439  return (MAX_TOK);
440}
Note: See TracBrowser for help on using the repository browser.