source: git/kernel/GBEngine/kLiftstd.cc @ 817917a

fieker-DuValspielwiese
Last change on this file since 817917a was 817917a, checked in by Hans Schoenemann <hannes@…>, 3 years ago
fix: example warkedPreimageStd
  • Property mode set to 100644
File size: 6.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5*  ABSTRACT -  Kernel: alg. of Buchberger
6*/
7
8#include "kernel/mod2.h"
9
10// define if no buckets should be used
11// #define NO_BUCKETS
12
13#include "kernel/GBEngine/kutil.h"
14#include "misc/options.h"
15#include "kernel/polys.h"
16#include "kernel/ideals.h"
17#include "kernel/GBEngine/kstd1.h"
18#include "kernel/GBEngine/khstd.h"
19#include "polys/kbuckets.h"
20#include "polys/prCopy.h"
21#include "polys/weight.h"
22#include "misc/intvec.h"
23#ifdef HAVE_PLURAL
24#include "polys/nc/nc.h"
25#endif
26
27static poly kSplitAt(int k,poly p, const ring r)
28{
29  if((p==NULL) || (p->next==NULL)) return NULL;
30  while(p_GetComp(p->next,r)<=k)
31  {
32    pIter(p);
33    if (p->next==NULL) return NULL;
34  }
35  poly t=p->next;
36  p->next=NULL;
37  return t;
38}
39static poly kSplitAt(int k,TObject* h,int *l,kStrategy strat)
40{
41  poly p;
42  if (h->t_p==NULL) h->GetLmTailRing();
43  p=h->t_p;
44  if ((p==NULL) ||(p->next==NULL)) return NULL;
45  int ll=1;
46  const ring tailRing=strat->tailRing;
47  while(p_GetComp(p->next,tailRing)<=k)
48  {
49    pIter(p);
50    *l=ll;
51    if ((p==NULL)||(p->next==NULL)) return NULL;
52    ll++;
53  }
54  poly t=p->next;
55  p->next=NULL;
56  *l=ll;
57  return t;
58}
59static poly kSplitAt(int k,LObject* h,kStrategy strat)
60{
61  poly p,pr;
62  int l;
63  if (h->bucket!=NULL)
64  {
65    kBucketClear(h->bucket,&p,&l);
66    pr=p;
67  }
68  else
69  { 
70    if (h->t_p!=NULL) h->GetLmTailRing();
71    p=h->t_p;
72  }
73  const ring tailRing=strat->tailRing;
74  if (p_GetComp(p,tailRing)>k)
75  {
76    return p;
77  }
78  if (p->next==NULL) return NULL;
79  while(p_GetComp(p->next,tailRing)<=k) pIter(p);
80  poly t=p->next;
81  p->next=NULL;
82  if (h->bucket!=NULL)
83  {
84    l=pLength(pr);
85    kBucketInit(h->bucket,pr,l);
86  }
87  return t;
88}
89static poly kTailAt(int k,TObject* h,kStrategy strat)
90{
91  poly p;
92  if (h->t_p!=NULL)
93    p=h->t_p;
94  else
95    p=h->p;
96  assume(p_GetComp(p,strat->tailRing)<=k);
97  if (p->next==NULL) return NULL;
98  const ring tailRing=strat->tailRing;
99  while(p_GetComp(p->next,tailRing)<=k) pIter(p);
100  poly t=p->next;
101  return t;
102}
103static void kAppend(poly t,TObject* h)
104{
105  poly p;
106  if (h->t_p!=NULL)
107    p=h->t_p;
108  else
109    p=h->p;
110  while(p->next!=NULL) pIter(p);
111  p->next=t;
112  if ((h->p!=NULL)&&(h->t_p!=NULL)) pNext(h->p)=pNext(h->t_p);
113}
114static void kAppend(poly t,LObject* h)
115{
116  if (h->bucket!=NULL)
117  {
118    int l=-1;
119    kBucket_Add_q(h->bucket,t,&l);
120  }
121  else
122  {
123    poly p;
124    if (h->t_p!=NULL)
125      p=h->t_p;
126    else
127      p=h->p;
128    while(p->next!=NULL) pIter(p);
129    p->next=t;
130    if ((h->p!=NULL)&&(h->t_p!=NULL)) pNext(h->p)=pNext(h->t_p);
131  }
132}
133
134static poly lazyComp(number* A, poly* M,poly* T,int index,poly s,int *l,const ring tailR)
135{
136  if ((TEST_OPT_PROT) && (index>0)) { Print("<%d>",index+1); mflush(); }
137  kBucket_pt b=kBucketCreate(tailR);
138  kBucketInit(b,s,pLength(s));
139  for(int i=0;i<index;i++)
140  {
141    kBucket_Mult_n(b,A[i]);
142    n_Delete(&A[i],tailR->cf);
143    poly tt=T[i];
144    if (tt!=NULL)
145    {
146      int dummy=pLength(tt);
147      kBucket_Minus_m_Mult_p(b,M[i],tt,&dummy);
148    }
149    p_Delete(&M[i],tailR);
150  }
151  poly p;
152  kBucketClear(b,&p,l);
153  kBucketDestroy(&b);
154  return p;
155}
156/*2
157*  reduction procedure for the sugar-strategy (honey)
158* reduces h with elements from T choosing first possible
159* element in T with respect to the given ecart
160*/
161int redLiftstd (LObject* h, kStrategy strat)
162{
163  if (strat->tl<0) return 1;
164  assume(h->FDeg == h->pFDeg());
165  poly h_p;
166  int i,j,at,pass,ei, ii, h_d,ci;
167  unsigned long not_sev;
168  long reddeg,d;
169  number A[500];
170  poly C[500];
171  poly T[500];
172  memset(T,0,sizeof(T));
173  memset(C,0,sizeof(T));
174  const ring tailRing=strat->tailRing;
175
176  pass = j = 0;
177  d = reddeg = h->GetpFDeg() + h->ecart;
178  h->SetShortExpVector();
179  int li;
180  h_p = h->GetLmTailRing();
181  not_sev = ~ h->sev;
182
183  // split h into mina part (h) and tail (h_tail)
184  poly h_tail=kSplitAt(strat->syzComp,h,strat);
185  // fix h-pLength
186  h->pLength=0;
187  // remove content
188  //number cont;
189  //p_Content_n(h_p,cont,strat->tailRing);
190  //if (!n_IsOne(cont,strat->tailRing))
191  //  h_tail=p_Div_nn(h_tail,cont,tailRing);
192
193  h->PrepareRed(strat->use_buckets);
194  loop
195  {
196    j=kFindDivisibleByInT(strat, h);
197    if (j < 0)
198    {
199      // lazy computation:
200      int l;
201      poly p=lazyComp(A,C,T,pass,h_tail,&l,strat->tailRing);
202      kBucket_Add_q(h->bucket,p,&l);
203      return 1;
204    }
205
206    ei = strat->T[j].ecart;
207    li = strat->T[j].pLength;
208    ci = nSize(pGetCoeff(strat->T[j].p));
209    if (li<=0) li=strat->T[j].GetpLength();
210    ii = j;
211    /*
212     * the polynomial to reduce with (up to the moment) is;
213     * pi with ecart ei (T[ii])
214     */
215    i = j;
216    if (TEST_OPT_LENGTH)
217    loop
218    {
219      /*- possible with respect to ecart, minimal nSize -*/
220      i++;
221      if (i > strat->tl)
222        break;
223      //if (ei < h->ecart)
224      //  break;
225      if (li==1)
226        break;
227      if ((((strat->T[i].ecart < ei) && (ei> h->ecart))
228         || ((strat->T[i].ecart <= h->ecart)
229            && (strat->T[i].pLength <= li)
230            && (nSize(pGetCoeff(strat->T[i].p)) <ci)))
231         &&
232          p_LmShortDivisibleBy(strat->T[i].GetLmTailRing(), strat->sevT[i],
233                               h_p, not_sev, tailRing))
234      {
235        /*
236         * the polynomial to reduce with is now;
237         */
238        ei = strat->T[i].ecart;
239        li = strat->T[i].pLength;
240        if (li<=0) li=strat->T[i].GetpLength();
241        ii = i;
242      }
243    }
244
245    /*
246     * end of search: have to reduce with pi
247     */
248#ifdef KDEBUG
249    if (TEST_OPT_DEBUG)
250    {
251      PrintS("red:");
252      h->wrp();
253      Print("\nwith T[%d]:",ii);
254      strat->T[ii].wrp();
255    }
256#endif
257    assume(strat->fromT == FALSE);
258
259    //strat->T[ii].pCleardenom();
260    // split T[ii]:
261    // remember pLength of strat->T[ii]
262    int l_orig=strat->T[ii].pLength;
263    // split strat->T[ii]
264    poly T_tail=kSplitAt(strat->syzComp,&strat->T[ii],&strat->T[ii].pLength,strat);
265    ksReducePoly(h,&(strat->T[ii]),NULL,&A[pass],&C[pass], strat);
266    // restore T[ii]:
267    kAppend(T_tail,&strat->T[ii]);
268    strat->T[ii].pLength=l_orig;
269    // store T_tail
270    T[pass]=T_tail;
271    // delayed computation: A[pass]*tail-M[pass]*T[pass]
272#ifdef KDEBUG
273    if (TEST_OPT_DEBUG)
274    {
275      PrintS("\nto:");
276      h->wrp();
277      PrintLn();
278    }
279#endif
280    if(h->IsNull())
281    {
282      // clean up A,C,h_tail:
283      for(int i=0;i<=pass;i++)
284      {
285        n_Delete(&A[i],tailRing->cf);
286        p_Delete(&C[i],tailRing);
287      }
288      p_Delete(&h_tail,tailRing);
289      kDeleteLcm(h);
290      h->Clear();
291      return 0;
292    }
293    h->SetShortExpVector();
294    not_sev = ~ h->sev;
295    h_d = h->SetpFDeg();
296    /* compute the ecart */
297    if (ei <= h->ecart)
298      h->ecart = d-h_d;
299    else
300      h->ecart = d-h_d+ei-h->ecart;
301
302    /*
303     * try to reduce the s-polynomial h
304     *test first whether h should go to the lazyset L
305     *-if the degree jumps
306     *-if the number of pre-defined reductions jumps
307     */
308    pass++;
309    d = h_d + h->ecart;
310    if (pass%64==0) kBucketCanonicalize(h->bucket);
311  }
312}
Note: See TracBrowser for help on using the repository browser.