source: git/kernel/shiftgb.cc @ abe5c8

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