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

fieker-DuValspielwiese
Last change on this file since cf5fba9 was cf5fba9, checked in by Hans Schoenemann <hannes@…>, 8 years ago
fix: name of lib (customstd.so)
  • Property mode set to 100644
File size: 13.1 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
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
207// static BOOLEAN abortIfMonomial_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//       // if a term is contained in the ideal, abort std computation
217//       // and store the output in idealCache to be returned
218//       while ((strat->Ll >= 0))
219//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
220//       std::cout << "aborting!" << std::endl;
221//       return FALSE;
222//     }
223//   }
224//   else
225//   {
226//     poly p=strat->P.t_p;
227//     if (pNext(p)==NULL)
228//     {
229//       // if a term is contained in the ideal, abort std computation
230//       // and store the output in idealCache to be returned
231//       while ((strat->Ll >= 0))
232//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
233//       std::cout << "aborting!" << std::endl;
234//       return FALSE;
235//     }
236//   }
237//   return b; // return TRUE if sp was changed, FALSE if not
238// }
239// static BOOLEAN abortifmonomialstd(leftv res, leftv args)
240// {
241//   if (args!=NULL)
242//   {
243//     if ((args->Typ()==IDEAL_CMD) && (args->next==NULL))
244//     {
245//       ideal I=(ideal)args->Data();
246//       idealCache = NULL;
247//       I=kStd(I,currRing->qideal,testHomog,NULL,NULL,0,0,NULL,abortIfMonomial_sp);
248//       res->rtyp=IDEAL_CMD;
249//       if (idealCache)
250//         res->data=(char*)idealCache;
251//       else
252//       {
253//         idSkipZeroes(I);
254//         res->data=(char*)I;
255//       }
256//       return FALSE;
257//     }
258//   }
259//   WerrorS("abortifmonomialstd: unexpected parameters");
260//   return TRUE;
261// }
262
263
264// static long wDeg(const poly p, const ring r)
265// {
266//   if (r->order[0] == ringorder_lp)
267//     return p_GetExp(p,1,currRing);
268//   if (r->order[0] == ringorder_ls)
269//     return -p_GetExp(p,1,currRing);
270
271//   if (r->order[0] == ringorder_dp)
272//   {
273//     long d = 0;
274//     for (int i=1; i<=rVar(r); i++)
275//       d = d + p_GetExp(p,i,r);
276//     return d;
277//   }
278//   if (r->order[0] == ringorder_wp || r->order[0] == ringorder_a)
279//   {
280//     long d = 0;
281//     for (int i=r->block0[0]; i<=r->block1[0]; i++)
282//       d = d + p_GetExp(p,i,r)*r->wvhdl[0][i-1];
283//     return d;
284//   }
285//   if (r->order[0] == ringorder_ws)
286//   {
287//     long d = 0;
288//     for (int i=r->block0[0]; i<=r->block1[0]; i++)
289//       d = d - p_GetExp(p,i,r)*r->wvhdl[0][i-1];
290//     return d;
291//   }
292// }
293
294// static bool isInitialFormMonomial(const poly g, const ring r)
295// {
296//   if (g->next==NULL)
297//     return true;
298//   return wDeg(g,r)!=wDeg(g->next,r);
299// }
300
301// //------------------------------------------------------------------------
302// // routine that checks whether the initial form is a monomial,
303// // breaks computation if it finds one, writing the element into idealCache
304// static BOOLEAN sat_sp_initial(kStrategy strat)
305// {
306//   BOOLEAN b = FALSE; // set b to TRUE, if spoly was changed,
307//                      // let it remain FALSE otherwise
308//   if (strat->P.t_p==NULL)
309//   {
310//     poly p=strat->P.p;
311//     if (pNext(p)==NULL)
312//     {
313//       // if a term is contained in the ideal, abort std computation
314//       // and store the output in idealCache to be returned
315//       while ((strat->Ll >= 0))
316//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
317//       idealCache = idInit(1);
318//       idealCache->m[0] = p_One(currRing);
319//       return FALSE;
320//     }
321//     if (isInitialFormMonomial(p,currRing))
322//     {
323//       while ((strat->Ll >= 0))
324//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
325//       idealCache = idInit(1);
326//       idealCache->m[0] = p_Copy(p,currRing);
327//     }
328//     int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
329//     int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
330//     p_GetExpV(p,mm,currRing);
331//     int m_null=0;
332//     while(p!=NULL)
333//     {
334//       m_null=0;
335//       p_GetExpV(p,m0,currRing);
336//       for(int i=1;i<=rVar(currRing);i++)
337//       {
338//         mm[i]=si_min(mm[i],m0[i]);
339//         if (mm[i]>0) m_null++;
340//       }
341//       if (m_null==0) break;
342//       pIter(p);
343//     }
344//     if (m_null>0)
345//     {
346//       std::cout << "simplifying!" << std::endl;
347//       p=p_Copy(strat->P.p,currRing);
348//       strat->P.p=p;
349//       while(p!=NULL)
350//       {
351//         for(int i=1;i<=rVar(currRing);i++)
352//           p_SubExp(p,i,mm[i],currRing);
353//         p_Setm(p,currRing);
354//         pIter(p);
355//       }
356//       b = TRUE;
357//     }
358//     omFree(mm);
359//     omFree(m0);
360//   }
361//   else
362//   {
363//     poly p=strat->P.t_p;
364//     if (pNext(p)==NULL)
365//     {
366//       // if a term is contained in the ideal, abort std computation
367//       // and store the output in idealCache to be returned
368//       while ((strat->Ll >= 0))
369//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
370//       idealCache = idInit(1);
371//       idealCache->m[0] = p_One(currRing);
372//       return FALSE;
373//     }
374//     if (isInitialFormMonomial(p,strat->tailRing))
375//     {
376//       while ((strat->Ll >= 0))
377//         deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
378//       nMapFunc identity = n_SetMap(strat->tailRing,currRing);
379//       idealCache = idInit(1);
380//       idealCache->m[0] = p_PermPoly(p,NULL,strat->tailRing,currRing,identity,NULL,0);
381//     }
382//     int *mm=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
383//     int *m0=(int*)omAlloc((1+rVar(currRing))*sizeof(int));
384//     p_GetExpV(p,mm,strat->tailRing);
385//     int m_null=0;
386//     while(p!=NULL)
387//     {
388//       m_null=0;
389//       p_GetExpV(p,m0,strat->tailRing);
390//       for(int i=1;i<=rVar(currRing);i++)
391//       {
392//         mm[i]=si_min(mm[i],m0[i]);
393//         if (mm[i]>0) m_null++;
394//       }
395//       if (m_null==0) break;
396//       pIter(p);
397//     }
398//     if (m_null>0)
399//     {
400//       std::cout << "simplifying!" << std::endl;
401//       p=p_Copy(strat->P.t_p,strat->tailRing);
402//       strat->P.t_p=p;
403//       strat->P.p=NULL;
404//       while(p!=NULL)
405//       {
406//         for(int i=1;i<=rVar(currRing);i++)
407//           p_SubExp(p,i,mm[i],strat->tailRing);
408//         p_Setm(p,strat->tailRing);
409//         pIter(p);
410//       }
411//       strat->P.GetP();
412//       b = TRUE;
413//     }
414//     omFree(mm);
415//     omFree(m0);
416//   }
417//   return b; // return TRUE if sp was changed, FALSE if not
418// }
419// static BOOLEAN satstdWithInitialCheck(leftv res, leftv args)
420// {
421//   if (args!=NULL)
422//   {
423//     if ((args->Typ()==IDEAL_CMD) && (args->next==NULL))
424//     {
425//       ideal I=(ideal)args->Data();
426//       idealCache = NULL;
427//       I=kStd(I,currRing->qideal,testHomog,NULL,NULL,0,0,NULL,sat_sp_initial);
428//       res->rtyp=IDEAL_CMD;
429//       if (idealCache)
430//         res->data=(char*)idealCache;
431//       else
432//         res->data=(char*)I;
433//       return FALSE;
434//     }
435//   }
436//   WerrorS("satstdWithInitialCheck: unexpected parameters");
437//   return TRUE;
438// }
439
440
441
442//------------------------------------------------------------------------
443// initialisation of the module
444extern "C" int SI_MOD_INIT(std_demo)(SModulFunctions* p)
445{
446  // p->iiAddCproc("std_demo","std_with_display",FALSE,std_with_display);
447  p->iiAddCproc("customstd.so","satstd",FALSE,satstd);
448  // p->iiAddCproc("std_demo","satstdWithInitialCheck",FALSE,satstdWithInitialCheck);
449  // p->iiAddCproc("std_demo","abortifmonomialstd",FALSE,abortifmonomialstd);
450  // PrintS("init of std_demo - type `listvar(Std_demo);` to its contents\n");
451  return (MAX_TOK);
452}
Note: See TracBrowser for help on using the repository browser.