My Project
kLiftstd.cc
Go to the documentation of this file.
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 "omalloc/omalloc.h"
14 #include "kernel/GBEngine/kutil.h"
15 #include "misc/options.h"
16 #include "kernel/polys.h"
17 #include "kernel/ideals.h"
18 #include "kernel/GBEngine/kstd1.h"
19 #include "kernel/GBEngine/khstd.h"
20 #include "polys/kbuckets.h"
21 #include "polys/prCopy.h"
22 #include "polys/weight.h"
23 #include "misc/intvec.h"
24 #ifdef HAVE_PLURAL
25 #include "polys/nc/nc.h"
26 #endif
27 
28 static poly kSplitAt(int k,TObject* h,kStrategy strat)
29 {
30  poly p;
31  if (h->t_p==NULL)
32  {
33  if (currRing!=strat->tailRing)
34  {
35  h->t_p=k_LmInit_currRing_2_tailRing(h->p, strat->tailRing);
36  p=h->t_p;
37  }
38  else
39  p=h->p;
40  }
41  else
42  p=h->t_p;
43  if (p->next==NULL) return NULL;
44  const ring tailRing=strat->tailRing;
45  while(p_GetComp(p->next,tailRing)<=k)
46  {
47  pIter(p);
48  if ((p==NULL)||(p->next==NULL))
49  {
50  h->pLength=0; // force re-computation
51  return NULL;
52  }
53  }
54  poly t=p->next;
55  p->next=NULL;
56  h->pLength=0; // force re-computation
57  if ((h->p!=NULL) && (h->t_p!=NULL)
58  && (pNext(h->p)!=pNext(h->t_p)))
59  {
60  pNext(h->p)=pNext(h->t_p);
61  }
62  return t;
63 }
64 static poly kSplitAt(int k,LObject* h,kStrategy strat)
65 {
66  poly p,pr,t=NULL;
67  int l;
68  if (h->bucket!=NULL)
69  {
70  kBucketClear(h->bucket,&p,&l);
71  pr=p;
72  }
73  else
74  {
75  if (h->t_p==NULL)
76  {
77  if (currRing!=strat->tailRing)
78  {
79  h->t_p=k_LmInit_currRing_2_tailRing(h->p, strat->tailRing);
80  p=h->t_p;
81  }
82  else
83  p=h->p;
84  }
85  else
86  p=h->t_p;
87  }
88  const ring tailRing=strat->tailRing;
89  if(p==NULL) return NULL;
90  if (p_GetComp(p,tailRing)>k)
91  {
92  return p;
93  }
94  if (p->next==NULL)
95  {
96  goto finish;
97  }
98  while(p_GetComp(p->next,tailRing)<=k)
99  {
100  pIter(p);
101  if (p->next==NULL) break;
102  }
103  t=p->next;
104  p->next=NULL;
105 finish:
106  if (h->bucket!=NULL)
107  {
108  l=pLength(pr);
109  kBucketInit(h->bucket,pr,l);
110  }
111  else
112  {
113  if ((h->p!=NULL) && (h->t_p!=NULL))
114  {
115  pNext(h->p)=pNext(h->t_p);
116  }
117  }
118 
119  return t;
120 }
121 static void kAppend(poly t,TObject* h)
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 static poly lazyComp(number* A, poly* M,poly* T,int index,poly s,int *l,const ring tailR)
133 {
134  if ((TEST_OPT_PROT) && (index>0)) { Print("<%d>",index+1); mflush(); }
135  kBucket_pt b=kBucketCreate(tailR);
136  kBucketInit(b,s,pLength(s));
137  int cnt=RED_CANONICALIZE;
138  for(int i=0;i<index;i++)
139  {
140  kBucket_Mult_n(b,A[i]);
141  n_Delete(&A[i],tailR->cf);
142  poly tt=T[i];
143  if (tt!=NULL)
144  {
145  cnt--;
146  int dummy=pLength(tt);
147  kBucket_Minus_m_Mult_p(b,M[i],tt,&dummy);
148  }
149  p_Delete(&M[i],tailR);
150  if (UNLIKELY(cnt==0))
151  {
152  cnt=RED_CANONICALIZE;
154  }
155  }
156  poly p;
157  kBucketClear(b,&p,l);
158  kBucketDestroy(&b);
159  return p;
160 }
161 
162 /*2
163 * reduction procedure for the sugar-strategy (honey)
164 * reduces h with elements from T choosing first possible
165 * element in T with respect to the given ecart
166 */
168 {
169  if (strat->tl<0) return 1;
170  assume(h->FDeg == h->pFDeg());
172  poly h_p;
173  int i,j,pass,ei, ii, h_d,ci;
174  unsigned long not_sev;
175  long reddeg,d;
176  #define START_REDUCE 512
177  int red_size=START_REDUCE;
178  number *A=(number*)omAlloc0(red_size*sizeof(number));
179  poly *C=(poly*)omAlloc0(red_size*sizeof(poly));
180  poly *T=(poly*)omAlloc0(red_size*sizeof(poly));
181  const ring tailRing=strat->tailRing;
182 
183  pass = j = 0;
184  d = reddeg = h->GetpFDeg() + h->ecart;
185  h->SetShortExpVector();
186  int li;
187  h_p = h->GetLmTailRing();
188  not_sev = ~ h->sev;
189 
190  // split h into mina part (h) and tail (h_tail)
191  poly h_tail=kSplitAt(strat->syzComp,h,strat);
192  // fix h-pLength
193  h->pLength=0;
194  // remove content
195  //number cont;
196  //p_Content_n(h_p,cont,strat->tailRing);
197  //if (!n_IsOne(cont,strat->tailRing))
198  // h_tail=p_Div_nn(h_tail,cont,tailRing);
199 
200  h->PrepareRed(strat->use_buckets);
201  loop
202  {
203  j=kFindDivisibleByInT(strat, h);
204  if (j < 0)
205  {
206  // lazy computation:
207  int l;
208  poly p=lazyComp(A,C,T,pass,h_tail,&l,strat->tailRing);
209  kBucket_Add_q(h->bucket,p,&l);
210  omFreeSize(A,red_size*sizeof(number));
211  omFreeSize(T,red_size*sizeof(poly));
212  omFreeSize(C,red_size*sizeof(poly));
213  return 1;
214  }
215 
216  ei = strat->T[j].ecart;
217  li = strat->T[j].pLength;
218  ci = nSize(pGetCoeff(strat->T[j].p));
219  ii = j;
220  /*
221  * the polynomial to reduce with (up to the moment) is;
222  * pi with ecart ei (T[ii])
223  */
224  i = j;
225  if (TEST_OPT_LENGTH)
226  {
227  if (li<=0) li=strat->T[j].GetpLength();
228  if (li>1)
229  loop
230  {
231  /*- possible with respect to ecart, minimal nSize -*/
232  i++;
233  if (i > strat->tl)
234  break;
235  //if (ei < h->ecart)
236  // break;
237  if ((((strat->T[i].ecart < ei) && (ei> h->ecart))
238  || ((strat->T[i].ecart <= h->ecart)
239  && (strat->T[i].pLength <= li)
240  && (nSize(pGetCoeff(strat->T[i].p)) <ci)))
241  &&
242  p_LmShortDivisibleBy(strat->T[i].GetLmTailRing(), strat->sevT[i],
243  h_p, not_sev, tailRing))
244  {
245  /*
246  * the polynomial to reduce with is now;
247  */
248  ei = strat->T[i].ecart;
249  li = strat->T[i].pLength;
250  if (li<=0) li=strat->T[i].GetpLength();
251  ii = i;
252  if (li==1) break;
253  }
254  }
255  }
256 
257  /*
258  * end of search: have to reduce with pi
259  */
260 #ifdef KDEBUG
261  if (TEST_OPT_DEBUG)
262  {
263  PrintS("red:");
264  h->wrp();
265  Print("\nwith T[%d]:",ii);
266  strat->T[ii].wrp();
267  }
268 #endif
269  assume(strat->fromT == FALSE);
270 
271  //strat->T[ii].pCleardenom();
272  // split T[ii]:
273  // remember pLength of strat->T[ii]
274  int l_orig=strat->T[ii].pLength;
275  // split strat->T[ii]
276  poly T_tail=kSplitAt(strat->syzComp,&strat->T[ii],strat);
277  h->pLength=0; // force re-computation of length
278  ksReducePoly(h,&(strat->T[ii]),NULL,&A[pass],&C[pass], strat);
279  // restore T[ii]:
280  kAppend(T_tail,&strat->T[ii]);
281  strat->T[ii].pLength=l_orig;
282  // store T_tail
283  T[pass]=T_tail;
284  // delayed computation: A[pass]*tail-M[pass]*T[pass]
285 #ifdef KDEBUG
286  if (TEST_OPT_DEBUG)
287  {
288  PrintS("\nto:");
289  h->wrp();
290  PrintLn();
291  }
292 #endif
293  if(h->IsNull())
294  {
295  // clean up A,C,h_tail:
296  for(int i=0;i<=pass;i++)
297  {
298  n_Delete(&A[i],tailRing->cf);
299  p_Delete(&C[i],tailRing);
300  }
301  p_Delete(&h_tail,tailRing);
302  kDeleteLcm(h);
303  h->Clear();
304  omFreeSize(A,red_size*sizeof(number));
305  omFreeSize(T,red_size*sizeof(poly));
306  omFreeSize(C,red_size*sizeof(poly));
307  return 0;
308  }
309  h->SetShortExpVector();
310  not_sev = ~ h->sev;
311  h_d = h->SetpFDeg();
312  /* compute the ecart */
313  if (ei <= h->ecart)
314  h->ecart = d-h_d;
315  else
316  h->ecart = d-h_d+ei-h->ecart;
317 
318  /*
319  * try to reduce the s-polynomial h
320  *test first whether h should go to the lazyset L
321  *-if the degree jumps
322  *-if the number of pre-defined reductions jumps
323  */
324  pass++;
325  d = h_d + h->ecart;
326  if (pass%RED_CANONICALIZE==0) kBucketCanonicalize(h->bucket);
327  // if cache is to small, double its size:
328  if (pass>=red_size-1)
329  {
330  A=(number*)omRealloc0Size(A,red_size*sizeof(number),2*red_size*sizeof(number));
331  C=(poly*)omRealloc0Size(C,red_size*sizeof(poly),2*red_size*sizeof(poly));
332  T=(poly*)omRealloc0Size(T,red_size*sizeof(poly),2*red_size*sizeof(poly));
333  if(TEST_OPT_PROT) {Print("+%d+",red_size);mflush();}
334  red_size*=2;
335  }
336  }
337 }
#define UNLIKELY(X)
Definition: auxiliary.h:404
#define FALSE
Definition: auxiliary.h:96
int l
Definition: cfEzgcd.cc:100
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
int p
Definition: cfModGcd.cc:4078
CanonicalForm b
Definition: cfModGcd.cc:4103
int syzComp
Definition: kutil.h:354
ring tailRing
Definition: kutil.h:343
TSet T
Definition: kutil.h:326
int tl
Definition: kutil.h:350
unsigned long * sevT
Definition: kutil.h:325
char use_buckets
Definition: kutil.h:383
char fromT
Definition: kutil.h:379
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:455
#define Print
Definition: emacs.cc:80
const CanonicalForm int s
Definition: facAbsFact.cc:51
int j
Definition: facHensel.cc:110
STATIC_VAR jList * T
Definition: janet.cc:30
STATIC_VAR Poly * h
Definition: janet.cc:971
KINLINE poly k_LmInit_currRing_2_tailRing(poly p, ring tailRing, omBin tailBin)
Definition: kInline.h:970
int redLiftstd(LObject *h, kStrategy strat)
Definition: kLiftstd.cc:167
static poly lazyComp(number *A, poly *M, poly *T, int index, poly s, int *l, const ring tailR)
Definition: kLiftstd.cc:132
static void kAppend(poly t, TObject *h)
Definition: kLiftstd.cc:121
#define START_REDUCE
static poly kSplitAt(int k, TObject *h, kStrategy strat)
Definition: kLiftstd.cc:28
void kBucketClear(kBucket_pt bucket, poly *p, int *length)
Definition: kbuckets.cc:521
void kBucket_Minus_m_Mult_p(kBucket_pt bucket, poly m, poly p, int *l, poly spNoether)
Bpoly == Bpoly - m*p; where m is a monom Does not destroy p and m assume (*l <= 0 || pLength(p) == *l...
Definition: kbuckets.cc:722
void kBucket_Mult_n(kBucket_pt bucket, number n)
Multiply Bucket by number ,i.e. Bpoly == n*Bpoly.
Definition: kbuckets.cc:598
void kBucketDestroy(kBucket_pt *bucket_pt)
Definition: kbuckets.cc:216
void kBucketInit(kBucket_pt bucket, poly lm, int length)
Definition: kbuckets.cc:493
kBucket_pt kBucketCreate(const ring bucket_ring)
Creation/Destruction of buckets.
Definition: kbuckets.cc:209
void kBucket_Add_q(kBucket_pt bucket, poly q, int *l)
Add to Bucket a poly ,i.e. Bpoly == q+Bpoly.
Definition: kbuckets.cc:660
int kBucketCanonicalize(kBucket_pt bucket)
Canonicalizes Bpoly, i.e. converts polys of buckets into one poly in one bucket: Returns number of bu...
int ksReducePoly(LObject *PR, TObject *PW, poly spNoether, number *coef, poly *mon, kStrategy strat)
Definition: kspoly.cc:187
int kFindDivisibleByInT(const kStrategy strat, const LObject *L, const int start)
return -1 if no divisor is found number of first divisor in T, otherwise
Definition: kstd2.cc:290
static void kDeleteLcm(LObject *P)
Definition: kutil.h:885
#define RED_CANONICALIZE
Definition: kutil.h:36
class sTObject TObject
Definition: kutil.h:57
class sLObject LObject
Definition: kutil.h:58
#define assume(x)
Definition: mod2.h:387
#define p_GetComp(p, r)
Definition: monomials.h:64
#define pIter(p)
Definition: monomials.h:37
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define nSize(n)
Definition: numbers.h:39
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define omRealloc0Size(addr, o_size, size)
Definition: omAllocDecl.h:221
#define NULL
Definition: omList.c:12
#define TEST_OPT_IDLIFT
Definition: options.h:129
#define TEST_OPT_LENGTH
Definition: options.h:131
#define TEST_OPT_PROT
Definition: options.h:103
#define TEST_OPT_DEBUG
Definition: options.h:108
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
static BOOLEAN p_LmShortDivisibleBy(poly a, unsigned long sev_a, poly b, unsigned long not_sev_b, const ring r)
Definition: p_polys.h:1909
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:873
static unsigned pLength(poly a)
Definition: p_polys.h:191
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatiblity layer for legacy polynomial operations (over currRing)
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
#define mflush()
Definition: reporter.h:58
#define A
Definition: sirandom.c:24
#define M
Definition: sirandom.c:25
#define loop
Definition: structs.h:75