source: git/kernel/shiftgb.cc @ fbc7cb

jengelh-datetimespielwiese
Last change on this file since fbc7cb was ba5e9e, checked in by Oleksandr Motsak <motsak@…>, 9 years ago
Changed configure-scripts to generate individual public config files for each package: resources, libpolys, singular (main) fix: sources should include correct corresponding config headers.
  • Property mode set to 100644
File size: 14.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: kernel: utils for shift GB and free GB
6*/
7
8#ifdef HAVE_CONFIG_H
9#include "singularconfig.h"
10#endif /* HAVE_CONFIG_H */
11#include <kernel/mod2.h>
12
13#ifdef HAVE_SHIFTBBA
14#include <kernel/febase.h>
15#include <polys/monomials/ring.h>
16#include <kernel/polys.h>
17#include <coeffs/numbers.h>
18#include <kernel/ideals.h>
19#include <polys/matpol.h>
20#include <polys/kbuckets.h>
21#include <kernel/kstd1.h>
22#include <polys/sbuckets.h>
23#include <polys/operations/p_Mult_q.h>
24#include <kernel/kutil.h>
25#include <kernel/structs.h>
26#include <omalloc/omalloc.h>
27#include <kernel/khstd.h>
28#include <polys/kbuckets.h>
29#include <polys/weight.h>
30#include <misc/intvec.h>
31#include <kernel/structs.h>
32#include <kernel/kInline.h>
33#include <kernel/stairc.h>
34#include <polys/weight.h>
35#include <misc/intvec.h>
36#include <kernel/timer.h>
37#include <kernel/shiftgb.h>
38#include <polys/nc/sca.h>
39
40
41#define freeT(A,v) omFreeSize((ADDRESS)A,(v+1)*sizeof(int))
42
43
44/* TODO: write p* stuff as instances of p_* for all the functions */
45/* p_* functions are new, p* are old */
46
47poly p_LPshiftT(poly p, int sh, int uptodeg, int lV, kStrategy strat, const ring r)
48{
49  /* assume shift takes place, shifts the poly p by sh */
50  /* p is like TObject: lm in currRing = r, tail in tailRing  */
51
52  if (p==NULL) return(p);
53
54  assume(p_LmCheckIsFromRing(p,r));
55  assume(p_CheckIsFromRing(pNext(p),strat->tailRing));
56
57  /* assume sh and uptodeg agree  TODO check */
58
59  if (sh == 0) return(p); /* the zero shift */
60
61  poly q   = NULL;
62  poly s   = p_mLPshift(p, sh, uptodeg, lV, r); // lm in currRing
63  poly pp = pNext(p);
64
65  while (pp != NULL)
66  {
67    q = p_Add_q(q, p_mLPshift(pp,sh,uptodeg,lV,strat->tailRing),strat->tailRing);
68    pIter(pp);
69  }
70  pNext(s) = q;
71  /* int version: returns TRUE if it was successful */
72  return(s);
73}
74
75
76poly p_LPshift(poly p, int sh, int uptodeg, int lV, const ring r)
77{
78  /* assume shift takes place */
79  /* shifts the poly p from the ring r by sh */
80
81  /* assume sh and uptodeg agree TODO check */
82
83  if (p==NULL) return(p);
84  if (sh == 0) return(p); /* the zero shift */
85
86  poly q  = NULL;
87  poly pp = p; // do not take copies
88  while (pp!=NULL)
89  {
90    q = p_Add_q(q, p_mLPshift(pp,sh,uptodeg,lV,r),r);
91    pIter(pp);
92  }
93  return(q);
94}
95
96poly p_mLPshift(poly p, int sh, int uptodeg, int lV, const ring r)
97{
98  /* p is a monomial from the ring r */
99
100  if (sh == 0) return(p); /* the zero shift */
101
102  if (sh < 0 )
103  {
104#ifdef PDEBUG
105    PrintS("pmLPshift: negative shift requested\n");
106#endif
107    return(NULL); /* violation, 2check */
108  }
109
110  int L = p_mLastVblock(p,lV,r);
111  if (L+sh-1 > uptodeg)
112  {
113#ifdef PDEBUG
114    PrintS("p_mLPshift: too big shift requested\n");
115#endif
116    return(NULL); /* violation, 2check */
117  }
118  int *e=(int *)omAlloc0((r->N+1)*sizeof(int));
119  int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
120  p_GetExpV(p,e,r);
121
122  int j;
123  //  for (j=1; j<=r->N; j++)
124  // L*lV gives the last position of the last block
125  for (j=1; j<= L*lV ; j++)
126  {
127    if (e[j]==1)
128    {
129      s[j + (sh*lV)] = e[j]; /* actually 1 */
130#ifdef PDEBUG
131      omCheckAddr(s);
132#endif
133    }
134#ifdef PDEBUG
135    else
136    {
137      if (e[j]!=0)
138      {
139         //         Print("p_mLPshift: ex[%d]=%d\n",j,e[j]);
140      }
141    }
142#endif
143  }
144  poly m = p_One(r);
145  p_SetExpV(m,s,r);
146  freeT(e, r->N);
147  freeT(s, r->N);
148  /*  pSetm(m); */ /* done in the pSetExpV */
149  /* think on the component and coefficient */
150  //  number c = pGetCoeff(p);
151  //  p_SetCoeff0(m,p_GetCoeff(p,r),r);
152  p_SetComp(m,p_GetComp(p,r),r); // component is preserved
153  p_SetCoeff0(m,n_Copy(p_GetCoeff(p,r),r->cf),r);  // coeff is preserved
154  return(m);
155}
156
157poly pLPshift(poly p, int sh, int uptodeg, int lV)
158{
159  /* assume shift takes place */
160  /* shifts the poly p by sh */
161  /* deletes p */
162
163  /* assume sh and uptodeg agree */
164
165  if (sh == 0) return(p); /* the zero shift */
166
167  poly q  = NULL;
168  poly pp = p; // pCopy(p);
169  while (pp!=NULL)
170  {
171    q = p_Add_q(q, pmLPshift(pp,sh,uptodeg,lV),currRing);
172    pIter(pp);
173  }
174  /* delete pp? */
175  /* int version: returns TRUE if it was successful */
176  return(q);
177}
178
179poly pmLPshift(poly p, int sh, int uptodeg, int lV)
180{
181  /* TODO: use a shortcut with p_ version */
182  /* pm is a monomial */
183
184  if (sh == 0) return(p); /* the zero shift */
185
186  if (sh < 0 )
187  {
188#ifdef PDEBUG
189    PrintS("pmLPshift: negative shift requested\n");
190#endif
191    return(NULL); /* violation, 2check */
192  }
193
194  int L = pmLastVblock(p,lV);
195  if (L+sh-1 > uptodeg)
196  {
197#ifdef PDEBUG
198    PrintS("pmLPshift: too big shift requested\n");
199#endif
200    return(NULL); /* violation, 2check */
201  }
202  int *e=(int *)omAlloc0((currRing->N+1)*sizeof(int));
203  int *s=(int *)omAlloc0((currRing->N+1)*sizeof(int));
204  pGetExpV(p,e);
205  number c = pGetCoeff(p);
206  int j;
207  for (j=1; j<=currRing->N; j++)
208  {
209    if (e[j]==1)
210    {
211      s[j + (sh*lV)] = e[j]; /* actually 1 */
212    }
213  }
214  poly m = pOne();
215  pSetExpV(m,s);
216  /*  pSetm(m); */ /* done in the pSetExpV */
217  /* think on the component */
218  pSetCoeff0(m,c);
219  freeT(e, currRing->N);
220  freeT(s, currRing->N);
221  return(m);
222}
223
224int pLastVblock(poly p, int lV)
225{
226  /* returns the number of maximal block */
227  /* appearing among the monomials of p */
228  /* the 0th block is the 1st one */
229  poly q = p; //p_Copy(p,currRing); /* need it ? */
230  int ans = 0;
231  int ansnew = 0;
232  while (q!=NULL)
233  {
234    ansnew = pmLastVblock(q,lV);
235    ans    = si_max(ans,ansnew);
236    pIter(q);
237  }
238  /* do not need to delete q */
239  return(ans);
240}
241
242int pmLastVblock(poly p, int lV)
243{
244  /* for a monomial p, returns the number of the last block */
245  /* where a nonzero exponent is sitting */
246  if (pIsConstantPoly(p))
247  {
248    return(int(0));
249  }
250  int *e=(int *)omAlloc0((currRing->N+1)*sizeof(int));
251  pGetExpV(p,e);
252  int j,b;
253  j = currRing->N;
254  while ( (!e[j]) && (j>=1) ) j--;
255  freeT(e, currRing->N);
256  if (j==0)
257  {
258#ifdef PDEBUG
259    PrintS("pmLastVblock: unexpected zero exponent vector\n");
260#endif
261    return(j);
262  }
263  b = (int)(j/lV) + 1; /* the number of the block, >=1 */
264  return (b);
265}
266
267int p_LastVblockT(poly p, int lV, kStrategy strat, const ring r)
268{
269  /* returns the number of maximal block */
270  /* appearing among the monomials of p */
271  /* the 0th block is the 1st one */
272
273  /* p is like TObject: lm in currRing = r, tail in tailRing  */
274  assume(p_LmCheckIsFromRing(p,r));
275  assume(p_CheckIsFromRing(pNext(p),strat->tailRing));
276
277  int ans = p_mLastVblock(p, lV, r); // Block of LM
278  poly q = pNext(p);
279  int ansnew = 0;
280  while (q != NULL)
281  {
282    ansnew = p_mLastVblock(q, lV, strat->tailRing);
283    ans       = si_max(ans,ansnew);
284    pIter(q);
285  }
286  /* do not need to delete q */
287  return(ans);
288}
289
290int p_LastVblock(poly p, int lV, const ring r)
291{
292  /* returns the number of maximal block */
293  /* appearing among the monomials of p */
294  /* the 0th block is the 1st one */
295  poly q = p; //p_Copy(p,currRing); /* need it ? */
296  int ans = 0;
297  int ansnew = 0;
298  while (q!=NULL)
299  {
300    ansnew = p_mLastVblock(q, lV, r);
301    ans    = si_max(ans,ansnew);
302    pIter(q);
303  }
304  /* do not need to delete q */
305  return(ans);
306}
307
308int p_mLastVblock(poly p, int lV, const ring r)
309{
310  /* for a monomial p, returns the number of the last block */
311  /* where a nonzero exponent is sitting */
312  if (p_LmIsConstant(p,r))
313  {
314    return(0);
315  }
316  int *e=(int *)omAlloc0((r->N+1)*sizeof(int));
317  p_GetExpV(p,e,r);
318  int j,b;
319  j = r->N;
320  while ( (!e[j]) && (j>=1) ) j--;
321  if (j==0)
322  {
323#ifdef PDEBUG
324    PrintS("pmLastVblock: unexpected zero exponent vector\n");
325#endif
326    return(j);
327  }
328  b = (int)((j+lV-1)/lV); /* the number of the block, >=1 */
329  freeT(e,r->N);
330  return (b);
331}
332
333int pFirstVblock(poly p, int lV)
334{
335  /* returns the number of maximal block */
336  /* appearing among the monomials of p */
337  /* the 0th block is the 1st one */
338  poly q = p; //p_Copy(p,currRing); /* need it ? */
339  int ans = 0;
340  int ansnew = 0;
341  while (q!=NULL)
342  {
343    ansnew = pmFirstVblock(q,lV);
344    ans    = si_min(ans,ansnew);
345    pIter(q);
346  }
347  /* do not need to delete q */
348  return(ans);
349}
350
351int pmFirstVblock(poly p, int lV)
352{
353  if (pIsConstantPoly(p))
354  {
355    return(int(0));
356  }
357  /* for a monomial p, returns the number of the first block */
358  /* where a nonzero exponent is sitting */
359  int *e=(int *)omAlloc0((currRing->N+1)*sizeof(int));
360  pGetExpV(p,e);
361  int j,b;
362  j = 1;
363  while ( (!e[j]) && (j<=currRing->N-1) ) j++;
364  if (j==currRing->N + 1)
365  {
366#ifdef PDEBUG
367    PrintS("pmFirstVblock: unexpected zero exponent vector\n");
368#endif
369    return(j);
370  }
371  b = (int)(j/lV)+1; /* the number of the block, 1<= N <= currRing->N  */
372  return (b);
373}
374
375  /* there should be two routines: */
376  /* 1. test place-squarefreeness: in homog this suffices: isInV */
377  /* 2. test the presence of a hole -> in the tail??? */
378
379int isInV(poly p, int lV)
380{
381  /* investigate only the leading monomial of p in currRing */
382  if ( pIsConstant(p) ) return(1);
383  if (lV <= 0) return(0);
384  /* returns 1 iff p is in V */
385  /* that is in each block up to a certain one there is only one nonzero exponent */
386  /* lV = the length of V = the number of orig vars */
387  int *e = (int *)omAlloc0((currRing->N+1)*sizeof(int));
388  int  b = (int)((currRing->N +lV-1)/lV); /* the number of blocks */
389  //int b  = (int)(currRing->N)/lV;
390  int *B = (int *)omAlloc0((b+1)*sizeof(int)); /* the num of elements in a block */
391  pGetExpV(p,e);
392  int i,j;
393  for (j=1; j<=b; j++)
394  {
395    /* we go through all the vars */
396    /* by blocks in lV vars */
397    for (i=(j-1)*lV + 1; i<= j*lV; i++)
398    {
399      if (e[i]) B[j] = B[j]+1;
400    }
401  }
402  //  j = b;
403  //  while ( (!B[j]) && (j>=1)) j--;
404  for (j=b; j>=1; j--)
405  {
406    if (B[j]!=0) break;
407  }
408  /* do not need e anymore */
409  freeT(e, currRing->N);
410
411  if (j==0) goto ret_true;
412//   {
413//     /* it is a zero exp vector, which is in V */
414//     freeT(B, b);
415//     return(1);
416//   }
417  /* now B[j] != 0 and we test place-squarefreeness */
418  for (; j>=1; j--)
419  {
420    if (B[j]!=1)
421    {
422      freeT(B, b);
423      return(0);
424    }
425  }
426 ret_true:
427  freeT(B, b);
428  return(1);
429}
430
431int poly_isInV(poly p, int lV)
432{
433  /* tests whether the whole polynomial p in in V */
434  poly q = p;
435  while (q!=NULL)
436  {
437    if ( !isInV(q,lV) )
438    {
439      return(0);
440    }
441    q = pNext(q);
442  }
443  return(1);
444}
445
446int ideal_isInV(ideal I, int lV)
447{
448  /* tests whether each polynomial of an ideal I lies in in V */
449  int i;
450  int s    = IDELEMS(I)-1;
451  for(i = 0; i <= s; i++)
452  {
453    if ( !poly_isInV(I->m[i],lV) )
454    {
455      return(0);
456    }
457  }
458  return(1);
459}
460
461
462int itoInsert(poly p, int uptodeg, int lV, const ring r)
463{
464  /* for poly in lmCR/tailTR presentation */
465  /* the below situation (commented out) might happen! */
466//   if (r == currRing)
467//   {
468//     "Current ring is not expected in toInsert";
469//     return(0);
470//   }
471  /* compute the number of insertions */
472  int i = p_mLastVblock(p, lV, currRing);
473  if (pNext(p) != NULL)
474  {
475    i = si_max(i, p_LastVblock(pNext(p), lV, r) );
476  }
477  //  i = uptodeg  - i +1;
478  i = uptodeg  - i;
479  //  p_wrp(p,currRing,r); Print("----i:%d",i); PrintLn();
480  return(i);
481}
482
483poly p_ShrinkT(poly p, int lV, kStrategy strat, const ring r)
484//poly p_Shrink(poly p, int uptodeg, int lV, kStrategy strat, const ring r)
485{
486  /* p is like TObject: lm in currRing = r, tail in tailRing  */
487  /* proc shrinks the poly p in ring r */
488  /* lV = the length of V = the number of orig vars */
489  /* check assumes/exceptions */
490  /* r->N is a multiple of lV */
491
492  if (p==NULL) return(p);
493
494  assume(p_LmCheckIsFromRing(p,r));
495  assume(p_CheckIsFromRing(pNext(p),strat->tailRing));
496
497  poly q   = NULL;
498  poly s   = p_mShrink(p, lV, r); // lm in currRing
499  poly pp = pNext(p);
500
501  while (pp != NULL)
502  {
503    //    q = p_Add_q(q, p_mShrink(pp,uptodeg,lV,strat->tailRing),strat->tailRing);
504    q = p_Add_q(q, p_mShrink(pp,lV,strat->tailRing),strat->tailRing);
505    pIter(pp);
506  }
507  pNext(s) = q;
508  return(s);
509}
510
511poly p_Shrink(poly p, int lV, const ring r)
512{
513  /* proc shrinks the poly p in ring r */
514  /* lV = the length of V = the number of orig vars */
515  /* check assumes/exceptions */
516  /* r->N is a multiple of lV */
517
518  if (p==NULL) return(p);
519  assume(p_CheckIsFromRing(p,r));
520  poly q = NULL;
521  poly pp = p;
522
523  while (pp != NULL)
524  {
525    q = p_Add_q(q, p_mShrink(pp,lV,r),r);
526    pIter(pp);
527  }
528  return(q);
529}
530
531poly p_mShrink(poly p, int lV, const ring r)
532{
533  /* shrinks the monomial p in ring r */
534  /* lV = the length of V = the number of orig vars */
535
536  /* check assumes/exceptions */
537  /* r->N is a multiple of lV */
538
539  int *e = (int *)omAlloc0((r->N+1)*sizeof(int));
540  int  b = (int)((r->N +lV-1)/lV); /* the number of blocks */
541  //  int *B = (int *)omAlloc0((b+1)*sizeof(int)); /* the num of elements in a block */
542  int *S = (int *)omAlloc0((r->N+1)*sizeof(int)); /* the shrinked exponent */
543  p_GetExpV(p,e,r);
544  int i,j; int cnt = 1; //counter for blocks in S
545  for (j=1; j<=b; j++)
546  {
547    /* we go through all the vars */
548    /* by blocks in lV vars */
549    for (i=(j-1)*lV + 1; i<= j*lV; i++)
550    {
551      if (e[i]==1)
552      {
553         //      B[j] = B[j]+1; // for control in V?
554         S[(cnt-1)*lV + (i - (j-1)*lV)] = e[i];
555         /* assuming we are in V, can interrupt here */
556         cnt++;
557         //  break; //results in incomplete shrink!
558         i = j*lV; // manual break under assumption p is in V
559      }
560    }
561  }
562#ifdef PDEBUG
563  //  Print("p_mShrink: cnt = [%d], b = %d\n",cnt,b);
564#endif
565  // cnt -1 <= b  must hold!
566  //  freeT(B, b);
567  poly s = p_One(r);
568  p_SetExpV(s,S,r);
569  freeT(e, r->N);
570  freeT(S, r->N);
571  /*  p_Setm(s,r); // done by p_SetExpV */
572  p_SetComp(s,p_GetComp(p,r),r); // component is preserved
573  p_SetCoeff(s,p_GetCoeff(p,r),r);  // coeff is preserved
574#ifdef PDEBUG
575  //  Print("p_mShrink: from "); p_wrp(p,r); Print(" to "); p_wrp(s,r); PrintLn();
576#endif
577  return(s);
578}
579
580/* shiftgb stuff */
581
582
583/*2
584 *if the leading term of p
585 *divides the leading term of some T[i] it will be canceled
586 */
587// static inline void clearSShift (poly p, unsigned long p_sev,int l, int* at, int* k,
588//                            kStrategy strat)
589// {
590//   assume(p_sev == pGetShortExpVector(p));
591//   if (!pLmShortDivisibleBy(p,p_sev, strat->T[*at].p, ~ strat->sevT[*at])) return;
592//   //  if (l>=strat->lenS[*at]) return;
593//   if (TEST_OPT_PROT)
594//     PrintS("!");
595//   mflush();
596//   //pDelete(&strat->S[*at]);
597//   deleteInS((*at),strat);
598//   (*at)--;
599//   (*k)--;
600// //  assume(lenS_correct(strat));
601// }
602
603/* remarks: cleanT : just deletion
604enlargeT: just reallocation */
605
606#endif
Note: See TracBrowser for help on using the repository browser.