source: git/libpolys/polys/monomials/p_polys.h @ bb6c8a

fieker-DuValspielwiese
Last change on this file since bb6c8a was bb6c8a, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: poly-cmp procedures should be in a correct order
  • Property mode set to 100644
File size: 47.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    p_polys.h
6 *  Purpose: declaration of poly stuf which are independent of
7 *           currRing
8 *  Author:  obachman (Olaf Bachmann)
9 *  Created: 9/00
10 *  Version: $Id$
11 *******************************************************************/
12#ifndef P_POLYS_H
13#define P_POLYS_H
14
15#include <polys/monomials/ring.h>
16#include <polys/monomials/monomials.h>
17#include <polys/monomials/polys-impl.h>
18#include <polys/templates/p_Procs.h>
19#include <polys/templates/p_Procs.h>
20#include <polys/sbuckets.h>
21
22/***************************************************************
23 *
24 * Primitives for accessing and setting fields of a poly
25 * poly must be != NULL
26 *
27 ***************************************************************/
28// next
29#define pNext(p)            ((p)->next)
30#define pIter(p)            ((p) = (p)->next)
31
32// coeff
33#define pGetCoeff(p)        ((p)->coef)
34// deletes old coeff before setting the new one
35#define pSetCoeff0(p,n)     (p)->coef=(n)
36#define p_GetCoeff(p,r)     pGetCoeff(p)
37#define p_SetCoeff0(p,n,r)  pSetCoeff0(p,n)
38
39#define __p_GetComp(p, r)   (p)->exp[r->pCompIndex]
40#define p_GetComp(p, r)    ((long) (r->pCompIndex >= 0 ? __p_GetComp(p, r) : 0))
41
42/***************************************************************
43 *
44 * Divisiblity tests, args must be != NULL, except for
45 * pDivisbleBy
46 *
47 ***************************************************************/
48unsigned long p_GetShortExpVector(poly a, ring r);
49
50/***************************************************************
51 *
52 * Misc things on polys
53 *
54 ***************************************************************/
55// return the maximal exponent of p in form of the maximal long var
56unsigned long p_GetMaxExpL(poly p, const ring r, unsigned long l_max = 0);
57// return monomial r such that GetExp(r,i) is maximum of all
58// monomials in p; coeff == 0, next == NULL, ord is not set
59poly p_GetMaxExpP(poly p, ring r);
60
61int p_MinDeg(poly p,intvec *w, const ring R);
62
63long p_DegW(poly p, const short *w, const ring R);
64
65// return TRUE if all monoms have the same component
66BOOLEAN   p_OneComp(poly p, ring r);
67
68// return i, if head depends only on var(i)
69int       p_IsPurePower(const poly p, const ring r);
70
71// return i, if poly depends only on var(i)
72int       p_IsUnivariate(poly p, const ring r);
73
74// set entry e[i] to 1 if var(i) occurs in p, ignore var(j) if e[j]>0
75// return #(e[i]>0)
76int      p_GetVariables(poly p, int * e, const ring r);
77
78// returns the poly representing the integer i
79poly      p_ISet(int i, ring r);
80
81// returns the poly representing the number n, destroys n
82poly      p_NSet(number n, ring r);
83
84/***************************************************************
85 *
86 * Copying/Deletion of polys: args may be NULL
87 *
88 ***************************************************************/
89
90// simply deletes monomials, does not free coeffs
91void p_ShallowDelete(poly *p, const ring r);
92
93
94
95/***************************************************************
96 *
97 * Copying/Deleteion of polys: args may be NULL
98 *  - p/q as arg mean a poly
99 *  - m a monomial
100 *  - n a number
101 *  - pp (resp. qq, mm, nn) means arg is constant
102 *  - p (resp, q, m, n)     means arg is destroyed
103 *
104 ***************************************************************/
105
106poly      p_Sub(poly a, poly b, const ring r);
107
108poly      p_Power(poly p, int i, const ring r);
109
110
111/***************************************************************
112 *
113 * PDEBUG stuff
114 *
115 ***************************************************************/
116#ifdef PDEBUG
117// Returns TRUE if m is monom of p, FALSE otherwise
118BOOLEAN pIsMonomOf(poly p, poly m);
119// Returns TRUE if p and q have common monoms
120BOOLEAN pHaveCommonMonoms(poly p, poly q);
121
122// p_Check* routines return TRUE if everything is ok,
123// else, they report error message and return false
124
125// check if Lm(p) is from ring r
126BOOLEAN p_LmCheckIsFromRing(poly p, ring r);
127// check if Lm(p) != NULL, r != NULL and initialized && Lm(p) is from r
128BOOLEAN p_LmCheckPolyRing(poly p, ring r);
129// check if all monoms of p are from ring r
130BOOLEAN p_CheckIsFromRing(poly p, ring r);
131// check r != NULL and initialized && all monoms of p are from r
132BOOLEAN p_CheckPolyRing(poly p, ring r);
133// check if r != NULL and initialized
134BOOLEAN p_CheckRing(ring r);
135// only do check if cond
136
137
138#define pIfThen(cond, check) do {if (cond) {check;}} while (0)
139
140BOOLEAN _p_Test(poly p, ring r, int level);
141BOOLEAN _p_LmTest(poly p, ring r, int level);
142BOOLEAN _pp_Test(poly p, ring lmRing, ring tailRing, int level);
143
144#define p_Test(p,r)     _p_Test(p, r, PDEBUG)
145#define p_LmTest(p,r)   _p_LmTest(p, r, PDEBUG)
146#define pp_Test(p, lmRing, tailRing)    _pp_Test(p, lmRing, tailRing, PDEBUG)
147
148#else // ! PDEBUG
149
150#define pIsMonomOf(p, q)        (TRUE)
151#define pHaveCommonMonoms(p, q) (TRUE)
152#define p_LmCheckIsFromRing(p,r)  ((void)0)
153#define p_LmCheckPolyRing(p,r)    ((void)0)
154#define p_CheckIsFromRing(p,r)  ((void)0)
155#define p_CheckPolyRing(p,r)    ((void)0)
156#define p_CheckRing(r)          ((void)0)
157#define P_CheckIf(cond, check)  ((void)0)
158
159#define p_Test(p,r)     (1)
160#define p_LmTest(p,r)   (1)
161#define pp_Test(p, lmRing, tailRing) (1)
162
163#endif
164
165/***************************************************************
166 *
167 * Misc stuff
168 *
169 ***************************************************************/
170/*2
171* returns the length of a (numbers of monomials)
172*/
173static inline int pLength(poly a)
174{
175  int l = 0;
176  while (a!=NULL)
177  {
178    pIter(a);
179    l++;
180  }
181  return l;
182}
183
184void      p_Norm(poly p1, const ring r);
185void      p_Normalize(poly p,const ring r);
186
187void      p_Content(poly p, const ring r);
188//void      p_SimpleContent(poly p, int s, const ring r);
189
190poly      p_Cleardenom(poly p, const ring r);
191void      p_Cleardenom_n(poly p, const ring r,number &c);
192number    p_GetAllDenom(poly ph, const ring r);
193
194int       pSize( poly p, const ring r );
195
196// homogenizes p by multiplying certain powers of the varnum-th variable
197poly      p_Homogen (poly p, int varnum, const ring r);
198
199BOOLEAN   p_IsHomogeneous (poly p, const ring r);
200
201static inline void p_Setm(poly p, const ring r);
202p_SetmProc p_GetSetmProc(ring r);
203
204poly      p_Subst(poly p, int n, poly e, const ring r);
205
206// TODO:
207#define p_SetmComp  p_Setm
208
209// component
210static inline  unsigned long p_SetComp(poly p, unsigned long c, ring r)
211{
212  p_LmCheckPolyRing2(p, r);
213  pAssume2(rRing_has_Comp(r));
214  __p_GetComp(p,r) = c;
215  return c;
216}
217// sets component of poly a to i, returns length of p
218static inline   void p_SetCompP(poly p, int i, ring r)
219{
220  if (p != NULL)
221  {
222#ifdef PDEBUG
223    poly q = p;
224    int l = 0;
225#endif
226
227    if (rOrd_SetCompRequiresSetm(r))
228    {
229      do
230      {
231        p_SetComp(p, i, r);
232        p_SetmComp(p, r);
233#ifdef PDEBUG
234        l++;
235#endif
236        pIter(p);
237      }
238      while (p != NULL);
239    }
240    else
241    {
242      do
243      {
244        p_SetComp(p, i, r);
245#ifdef PDEBUG
246        l++;
247#endif
248        pIter(p);
249      }
250      while(p != NULL);
251    }
252#ifdef PDEBUG
253    p_Test(q, r);
254    assume(l == pLength(q));
255#endif
256  }
257}
258
259static inline   void p_SetCompP(poly p, int i, ring lmRing, ring tailRing)
260{
261  if (p != NULL)
262  {
263    p_SetComp(p, i, lmRing);
264    p_SetmComp(p, lmRing);
265    p_SetCompP(pNext(p), i, tailRing);
266  }
267}
268
269// returns maximal column number in the modul element a (or 0)
270static inline long p_MaxComp(poly p, ring lmRing, ring tailRing)
271{
272  long result,i;
273
274  if(p==NULL) return 0;
275  result = p_GetComp(p, lmRing);
276  if (result != 0)
277  {
278    loop
279    {
280      pIter(p);
281      if(p==NULL) break;
282      i = p_GetComp(p, tailRing);
283      if (i>result) result = i;
284    }
285  }
286  return result;
287}
288
289static inline long p_MaxComp(poly p,ring lmRing) {return p_MaxComp(p,lmRing,lmRing);}
290
291static inline   long p_MinComp(poly p, ring lmRing, ring tailRing)
292{
293  long result,i;
294
295  if(p==NULL) return 0;
296  result = p_GetComp(p,lmRing);
297  if (result != 0)
298  {
299    loop
300    {
301      pIter(p);
302      if(p==NULL) break;
303      i = p_GetComp(p,tailRing);
304      if (i<result) result = i;
305    }
306  }
307  return result;
308}
309
310static inline long p_MinComp(poly p,ring lmRing) {return p_MinComp(p,lmRing,lmRing);}
311
312/***************************************************************
313 *
314 * poly things which are independent of ring
315 *
316 ***************************************************************/
317static inline poly      pLast(poly a, int &length);
318inline   poly      pLast(poly a) { int l; return pLast(a, l);}
319static inline poly pReverse(poly p)
320{
321  if (p == NULL || pNext(p) == NULL) return p;
322
323  poly q = pNext(p), // == pNext(p)
324    qn;
325  pNext(p) = NULL;
326  do
327  {
328    qn = pNext(q);
329    pNext(q) = p;
330    p = q;
331    q = qn;
332  }
333  while (qn != NULL);
334  return p;
335}
336void      pEnlargeSet(poly**p, int length, int increment);
337
338
339/***************************************************************
340 *
341 * I/O
342 *
343 ***************************************************************/
344char*     p_String(poly p, ring lmRing, ring tailRing);
345char*     p_String0(poly p, ring lmRing, ring tailRing);
346void      p_Write(poly p, ring lmRing, ring tailRing);
347void      p_Write0(poly p, ring lmRing, ring tailRing);
348void      p_wrp(poly p, ring lmRing, ring tailRing);
349
350/***************************************************************
351 *
352 * Degree stuff -- see p_polys.cc for explainations
353 *
354 ***************************************************************/
355extern pLDegProc pLDeg;
356extern pFDegProc pFDeg;
357long p_WFirstTotalDegree(poly p, ring r);
358long p_WTotaldegree(poly p, const ring r);
359long p_WDegree(poly p,const ring r);
360long pLDeg0(poly p,int *l, ring r);
361long pLDeg0c(poly p,int *l, ring r);
362long pLDegb(poly p,int *l, ring r);
363long pLDeg1(poly p,int *l, ring r);
364long pLDeg1c(poly p,int *l, ring r);
365long pLDeg1_Deg(poly p,int *l, ring r);
366long pLDeg1c_Deg(poly p,int *l, ring r);
367long pLDeg1_Totaldegree(poly p,int *l, ring r);
368long pLDeg1c_Totaldegree(poly p,int *l, ring r);
369long pLDeg1_WFirstTotalDegree(poly p,int *l, ring r);
370long pLDeg1c_WFirstTotalDegree(poly p,int *l, ring r);
371BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r);
372
373long p_Deg(poly a, const ring r);
374
375
376/***************************************************************
377 *
378 * Primitives for accessing and setting fields of a poly
379 *
380 ***************************************************************/
381
382static inline number p_SetCoeff(poly p, number n, ring r)
383{
384  p_LmCheckPolyRing2(p, r);
385  n_Delete(&(p->coef), r->cf);
386  (p)->coef=n;
387  return n;
388}
389
390// order
391static inline long p_GetOrder(poly p, ring r)
392{
393  p_LmCheckPolyRing2(p, r);
394  if (r->typ==NULL) return ((p)->exp[r->pOrdIndex]);
395  int i=0;
396  loop
397  {
398    switch(r->typ[i].ord_typ)
399    {
400      case ro_wp_neg:
401        return (((long)((p)->exp[r->pOrdIndex]))-POLY_NEGWEIGHT_OFFSET);
402      case ro_syzcomp:
403      case ro_syz:
404      case ro_cp:
405        i++;
406        break;
407      //case ro_dp:
408      //case ro_wp:
409      default:
410        return ((p)->exp[r->pOrdIndex]);
411    }
412  }
413}
414
415// Setm
416static inline void p_Setm(poly p, const ring r)
417{
418  p_CheckRing2(r);
419  r->p_Setm(p, r);
420}
421
422
423static inline unsigned long p_AddComp(poly p, unsigned long v, ring r)
424{
425  p_LmCheckPolyRing2(p, r);
426  pAssume2(rRing_has_Comp(r));
427  return __p_GetComp(p,r) += v;
428}
429static inline unsigned long p_SubComp(poly p, unsigned long v, ring r)
430{
431  p_LmCheckPolyRing2(p, r);
432  pAssume2(rRing_has_Comp(r));
433  _pPolyAssume2(__p_GetComp(p,r) >= v,p,r);
434  return __p_GetComp(p,r) -= v;
435}
436
437#ifndef HAVE_EXPSIZES
438
439/// get a single variable exponent
440/// @Note:
441/// the integer VarOffset encodes:
442/// 1. the position of a variable in the exponent vector p->exp (lower 24 bits)
443/// 2. number of bits to shift to the right in the upper 8 bits (which takes at most 6 bits for 64 bit)
444/// Thus VarOffset always has 2 zero higher bits!
445static inline long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
446{
447  pAssume2((VarOffset >> (24 + 6)) == 0);
448#if 0
449  int pos=(VarOffset & 0xffffff);
450  int bitpos=(VarOffset >> 24);
451  unsigned long exp=(p->exp[pos] >> bitmask) & iBitmask;
452  return exp;
453#else
454  return (long)
455         ((p->exp[(VarOffset & 0xffffff)] >> (VarOffset >> 24))
456          & iBitmask);
457#endif
458}
459
460
461/// set a single variable exponent
462/// @Note:
463/// VarOffset encodes the position in p->exp @see p_GetExp
464static inline unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
465{
466  pAssume2(e>=0);
467  pAssume2(e<=iBitmask);
468  pAssume2((VarOffset >> (24 + 6)) == 0);
469
470  // shift e to the left:
471  register int shift = VarOffset >> 24;
472  unsigned long ee = e << shift /*(VarOffset >> 24)*/;
473  // find the bits in the exponent vector
474  register int offset = (VarOffset & 0xffffff);
475  // clear the bits in the exponent vector:
476  p->exp[offset]  &= ~( iBitmask << shift );
477  // insert e with |
478  p->exp[ offset ] |= ee;
479  return e;
480}
481
482
483#else // #ifdef HAVE_EXPSIZES // EXPERIMENTAL!!!
484
485static inline unsigned long BitMask(unsigned long bitmask, int twobits)
486{
487  // bitmask = 00000111111111111
488  // 0 must give bitmask!
489  // 1, 2, 3 - anything like 00011..11
490  pAssume2((twobits >> 2) == 0);
491  static const unsigned long _bitmasks[4] = {-1, 0x7fff, 0x7f, 0x3};
492  return bitmask & _bitmasks[twobits];
493}
494
495
496/// @Note: we may add some more info (6 ) into VarOffset and thus encode
497static inline long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
498{
499  int pos  =(VarOffset & 0xffffff);
500  int hbyte= (VarOffset >> 24); // the highest byte
501  int bitpos = hbyte & 0x3f; // last 6 bits
502  long bitmask = BitMask(iBitmask, hbyte >> 6);
503
504  long exp=(p->exp[pos] >> bitpos) & bitmask;
505  return exp;
506
507}
508
509static inline long p_SetExp(poly p, const long e, const unsigned long iBitmask, const int VarOffset)
510{
511  pAssume2(e>=0);
512  pAssume2(e <= BitMask(iBitmask, VarOffset >> 30));
513
514  // shift e to the left:
515  register int hbyte = VarOffset >> 24;
516  int bitmask = BitMask(iBitmask, hbyte >> 6);
517  register int shift = hbyte & 0x3f;
518  long ee = e << shift;
519  // find the bits in the exponent vector
520  register int offset = (VarOffset & 0xffffff);
521  // clear the bits in the exponent vector:
522  p->exp[offset]  &= ~( bitmask << shift );
523  // insert e with |
524  p->exp[ offset ] |= ee;
525  return e;
526}
527
528#endif // #ifndef HAVE_EXPSIZES
529
530
531static inline long p_GetExp(const poly p, const ring r, const int VarOffset)
532{
533  p_LmCheckPolyRing2(p, r);
534  pAssume2(VarOffset != -1);
535  return p_GetExp(p, r->bitmask, VarOffset);
536}
537
538static inline long p_SetExp(poly p, const long e, const ring r, const int VarOffset)
539{
540  p_LmCheckPolyRing2(p, r);
541  pAssume2(VarOffset != -1);
542  return p_SetExp(p, e, r->bitmask, VarOffset);
543}
544
545
546
547/// get v^th exponent for a monomial
548static inline long p_GetExp(const poly p, const int v, const ring r)
549{
550  p_LmCheckPolyRing2(p, r);
551  pAssume2(v>0 && v <= r->N);
552  pAssume2(r->VarOffset[v] != -1);
553  return p_GetExp(p, r->bitmask, r->VarOffset[v]);
554}
555
556
557/// set v^th exponent for a monomial
558static inline long p_SetExp(poly p, const int v, const long e, const ring r)
559{
560  p_LmCheckPolyRing2(p, r);
561  pAssume2(v>0 && v <= r->N);
562  pAssume2(r->VarOffset[v] != -1);
563  return p_SetExp(p, e, r->bitmask, r->VarOffset[v]);
564}
565
566// the following should be implemented more efficiently
567static inline  long p_IncrExp(poly p, int v, ring r)
568{
569  p_LmCheckPolyRing2(p, r);
570  int e = p_GetExp(p,v,r);
571  e++;
572  return p_SetExp(p,v,e,r);
573}
574static inline  long p_DecrExp(poly p, int v, ring r)
575{
576  p_LmCheckPolyRing2(p, r);
577  int e = p_GetExp(p,v,r);
578  pAssume2(e > 0);
579  e--;
580  return p_SetExp(p,v,e,r);
581}
582static inline  long p_AddExp(poly p, int v, long ee, ring r)
583{
584  p_LmCheckPolyRing2(p, r);
585  int e = p_GetExp(p,v,r);
586  e += ee;
587  return p_SetExp(p,v,e,r);
588}
589static inline  long p_SubExp(poly p, int v, long ee, ring r)
590{
591  p_LmCheckPolyRing2(p, r);
592  long e = p_GetExp(p,v,r);
593  pAssume2(e >= ee);
594  e -= ee;
595  return p_SetExp(p,v,e,r);
596}
597static inline  long p_MultExp(poly p, int v, long ee, ring r)
598{
599  p_LmCheckPolyRing2(p, r);
600  long e = p_GetExp(p,v,r);
601  e *= ee;
602  return p_SetExp(p,v,e,r);
603}
604
605static inline long p_GetExpSum(poly p1, poly p2, int i, ring r)
606{
607  p_LmCheckPolyRing2(p1, r);
608  p_LmCheckPolyRing2(p2, r);
609  return p_GetExp(p1,i,r) + p_GetExp(p2,i,r);
610}
611static inline long p_GetExpDiff(poly p1, poly p2, int i, ring r)
612{
613  return p_GetExp(p1,i,r) - p_GetExp(p2,i,r);
614}
615
616static inline int p_Comp_k_n(poly a, poly b, int k, ring r)
617{
618  if ((a==NULL) || (b==NULL) ) return FALSE;
619  p_LmCheckPolyRing2(a, r);
620  p_LmCheckPolyRing2(b, r);
621  pAssume2(k > 0 && k <= r->N);
622  int i=k;
623  for(;i<=r->N;i++)
624  {
625    if (p_GetExp(a,i,r) != p_GetExp(b,i,r)) return FALSE;
626    //    if (a->exp[(r->VarOffset[i] & 0xffffff)] != b->exp[(r->VarOffset[i] & 0xffffff)]) return FALSE;
627  }
628  return TRUE;
629}
630
631
632/***************************************************************
633 *
634 * Allocation/Initalization/Deletion
635 *
636 ***************************************************************/
637static inline poly p_New(ring r, omBin bin)
638{
639  p_CheckRing2(r);
640  pAssume2(bin != NULL && r->PolyBin->sizeW == bin->sizeW);
641  poly p;
642  omTypeAllocBin(poly, p, bin);
643  p_SetRingOfLm(p, r);
644  return p;
645}
646
647static inline poly p_New(ring r)
648{
649  return p_New(r, r->PolyBin);
650}
651
652static inline void p_LmFree(poly p, ring r)
653{
654  p_LmCheckPolyRing2(p, r);
655  omFreeBinAddr(p);
656}
657static inline void p_LmFree(poly *p, ring r)
658{
659  p_LmCheckPolyRing2(*p, r);
660  poly h = *p;
661  *p = pNext(h);
662  omFreeBinAddr(h);
663}
664static inline poly p_LmFreeAndNext(poly p, ring r)
665{
666  p_LmCheckPolyRing2(p, r);
667  poly pnext = pNext(p);
668  omFreeBinAddr(p);
669  return pnext;
670}
671static inline void p_LmDelete(poly p, const ring r)
672{
673  p_LmCheckPolyRing2(p, r);
674  n_Delete(&pGetCoeff(p), r->cf);
675  omFreeBinAddr(p);
676}
677static inline void p_LmDelete(poly *p, const ring r)
678{
679  p_LmCheckPolyRing2(*p, r);
680  poly h = *p;
681  *p = pNext(h);
682  n_Delete(&pGetCoeff(h), r->cf);
683  omFreeBinAddr(h);
684}
685static inline poly p_LmDeleteAndNext(poly p, const ring r)
686{
687  p_LmCheckPolyRing2(p, r);
688  poly pnext = pNext(p);
689  n_Delete(&pGetCoeff(p), r->cf);
690  omFreeBinAddr(p);
691  return pnext;
692}
693
694/***************************************************************
695 *
696 * Misc routines
697 *
698 ***************************************************************/
699
700
701static inline unsigned long p_GetMaxExp(const poly p, const ring r)
702{
703  return p_GetMaxExp(p_GetMaxExpL(p, r), r);
704}
705
706static inline unsigned long p_GetMaxExp(const unsigned long l, const ring r)
707{
708  unsigned long bitmask = r->bitmask;
709  unsigned long max = (l & bitmask);
710  unsigned long j = r->ExpPerLong - 1;
711
712  if (j > 0)
713  {
714    unsigned long i = r->BitsPerExp;
715    long e;
716    loop
717    {
718      e = ((l >> i) & bitmask);
719      if ((unsigned long) e > max)
720        max = e;
721      j--;
722      if (j==0) break;
723      i += r->BitsPerExp;
724    }
725  }
726  return max;
727}
728
729static inline unsigned long
730p_GetTotalDegree(const unsigned long l, const ring r, const int number_of_exps)
731{
732  const unsigned long bitmask = r->bitmask;
733  unsigned long sum = (l & bitmask);
734  unsigned long j = number_of_exps - 1;
735
736  if (j > 0)
737  {
738    unsigned long i = r->BitsPerExp;
739    loop
740    {
741      sum += ((l >> i) & bitmask);
742      j--;
743      if (j==0) break;
744      i += r->BitsPerExp;
745    }
746  }
747  return sum;
748}
749
750static inline unsigned long
751p_GetTotalDegree(const unsigned long l, const ring r)
752{
753  return p_GetTotalDegree(l, r, r->ExpPerLong);
754}
755
756/***************************************************************
757 *
758 * Dispatcher to r->p_Procs, they do the tests/checks
759 *
760 ***************************************************************/
761// returns a copy of p
762static inline poly p_Copy(poly p, const ring r)
763{
764#ifdef PDEBUG
765  poly pp= r->p_Procs->p_Copy(p, r);
766  p_Test(pp,r);
767  return pp;
768#else
769  return r->p_Procs->p_Copy(p, r);
770#endif
771}
772
773static inline poly p_Head(poly p, const ring r)
774{
775  if (p == NULL) return NULL;
776  p_LmCheckPolyRing1(p, r);
777  poly np;
778  omTypeAllocBin(poly, np, r->PolyBin);
779  p_SetRingOfLm(np, r);
780  p_MemCopy_LengthGeneral(np->exp, p->exp, r->ExpL_Size);
781  pNext(np) = NULL;
782  pSetCoeff0(np, n_Copy(pGetCoeff(p), r->cf));
783  return np;
784}
785
786// returns a copy of p with Lm(p) from lmRing and Tail(p) from tailRing
787static inline poly p_Copy(poly p, const ring lmRing, const ring tailRing)
788{
789#ifndef PDEBUG
790  if (tailRing == lmRing)
791    return tailRing->p_Procs->p_Copy(p, tailRing);
792#endif
793  if (p != NULL)
794  {
795    poly pres = p_Head(p, lmRing);
796    pNext(pres) = tailRing->p_Procs->p_Copy(pNext(p), tailRing);
797    return pres;
798  }
799  else
800    return NULL;
801}
802
803// deletes *p, and sets *p to NULL
804static inline void p_Delete(poly *p, const ring r)
805{
806  r->p_Procs->p_Delete(p, r);
807}
808
809static inline void p_Delete(poly *p,  const ring lmRing, const ring tailRing)
810{
811#ifndef PDEBUG
812  if (tailRing == lmRing)
813  {
814    tailRing->p_Procs->p_Delete(p, tailRing);
815    return;
816  }
817#endif
818  if (*p != NULL)
819  {
820    if (pNext(*p) != NULL)
821      tailRing->p_Procs->p_Delete(&pNext(*p), tailRing);
822    p_LmDelete(p, lmRing);
823  }
824}
825
826// copys monomials of p, allocates new monomials from bin,
827// deletes monomoals of p
828static inline poly p_ShallowCopyDelete(poly p, const ring r, omBin bin)
829{
830  p_LmCheckPolyRing2(p, r);
831  pAssume2(r->PolyBin->sizeW == bin->sizeW);
832  return r->p_Procs->p_ShallowCopyDelete(p, r, bin);
833}
834
835// returns p+q, destroys p and q
836static inline poly p_Add_q(poly p, poly q, const ring r)
837{
838  int shorter;
839  return r->p_Procs->p_Add_q(p, q, shorter, r);
840}
841
842/// like p_Add_q, except that if lp == pLength(lp) lq == pLength(lq) then lp == pLength(p+q)
843static inline poly p_Add_q(poly p, poly q, int &lp, int lq, const ring r)
844{
845  int shorter;
846  poly res = r->p_Procs->p_Add_q(p, q, shorter, r);
847  lp = (lp + lq) - shorter;
848  return res;
849}
850
851// returns p*n, destroys p
852static inline poly p_Mult_nn(poly p, number n, const ring r)
853{
854  if (n_IsOne(n, r->cf))
855    return p;
856  else
857    return r->p_Procs->p_Mult_nn(p, n, r);
858}
859
860static inline poly p_Mult_nn(poly p, number n, const ring lmRing,
861                        const ring tailRing)
862{
863#ifndef PDEBUG
864  if (lmRing == tailRing)
865  {
866    return p_Mult_nn(p, n, tailRing);
867  }
868#endif
869  poly pnext = pNext(p);
870  pNext(p) = NULL;
871  p = lmRing->p_Procs->p_Mult_nn(p, n, lmRing);
872  pNext(p) = tailRing->p_Procs->p_Mult_nn(pnext, n, tailRing);
873  return p;
874}
875
876// returns p*n, does not destroy p
877static inline poly pp_Mult_nn(poly p, number n, const ring r)
878{
879  if (n_IsOne(n, r->cf))
880    return p_Copy(p, r);
881  else
882    return r->p_Procs->pp_Mult_nn(p, n, r);
883}
884
885// test if the monomial is a constant as a vector component
886// i.e., test if all exponents are zero
887static inline BOOLEAN p_LmIsConstantComp(const poly p, const ring r)
888{
889  //p_LmCheckPolyRing(p, r);
890  int i = r->VarL_Size - 1;
891
892  do
893  {
894    if (p->exp[r->VarL_Offset[i]] != 0)
895      return FALSE;
896    i--;
897  }
898  while (i >= 0);
899  return TRUE;
900}
901
902// test if monomial is a constant, i.e. if all exponents and the component
903// is zero
904static inline BOOLEAN p_LmIsConstant(const poly p, const ring r)
905{
906  if (p_LmIsConstantComp(p, r))
907    return (p_GetComp(p, r) == 0);
908  return FALSE;
909}
910
911// returns Copy(p)*m, does neither destroy p nor m
912static inline poly pp_Mult_mm(poly p, poly m, const ring r)
913{
914  if (p_LmIsConstant(m, r))
915    return pp_Mult_nn(p, pGetCoeff(m), r);
916  else
917  {
918    poly last;
919    return r->p_Procs->pp_Mult_mm(p, m, r, last);
920  }
921}
922
923// returns p*m, destroys p, const: m
924static inline poly p_Mult_mm(poly p, poly m, const ring r)
925{
926  if (p_LmIsConstant(m, r))
927    return p_Mult_nn(p, pGetCoeff(m), r);
928  else
929    return r->p_Procs->p_Mult_mm(p, m, r);
930}
931
932// return p - m*Copy(q), destroys p; const: p,m
933static inline poly p_Minus_mm_Mult_qq(poly p, poly m, poly q, const ring r)
934{
935#ifdef HAVE_PLURAL
936  if (rIsPluralRing(r))
937  {
938    int lp, lq;
939    poly spNoether;
940    return nc_p_Minus_mm_Mult_qq(p, m, q, lp, lq, spNoether, r);
941  }
942#endif
943
944  int shorter;
945  poly last;
946
947  return r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, NULL, r, last); // !!!
948}
949
950// like p_Minus_mm_Mult_qq, except that if lp == pLength(lp) lq == pLength(lq)
951// then lp == pLength(p -m*q)
952static inline poly p_Minus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq,
953                                 poly spNoether, const ring r)
954{
955#ifdef HAVE_PLURAL
956  if (rIsPluralRing(r))
957     return nc_p_Minus_mm_Mult_qq(p, m, q, lp, lq, spNoether, r);
958#endif
959
960  int shorter;
961  poly last,res;
962  res = r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, spNoether, r, last);
963  lp = (lp + lq) - shorter;
964  return res;
965}
966
967// returns p*Coeff(m) for such monomials pm of p, for which m is divisble by pm
968static inline poly pp_Mult_Coeff_mm_DivSelect(poly p, const poly m, const ring r)
969{
970  int shorter;
971  return r->p_Procs->pp_Mult_Coeff_mm_DivSelect(p, m, shorter, r);
972}
973
974// returns p*Coeff(m) for such monomials pm of p, for which m is divisble by pm
975// if lp is length of p on input then lp is length of returned poly on output
976static inline poly pp_Mult_Coeff_mm_DivSelect(poly p, int &lp, const poly m, const ring r)
977{
978  int shorter;
979  poly pp = r->p_Procs->pp_Mult_Coeff_mm_DivSelect(p, m, shorter, r);
980  lp -= shorter;
981  return pp;
982}
983
984// returns -p, destroys p
985static inline poly p_Neg(poly p, const ring r)
986{
987  return r->p_Procs->p_Neg(p, r);
988}
989
990extern poly  _p_Mult_q(poly p, poly q, const int copy, const ring r);
991// returns p*q, destroys p and q
992static inline poly p_Mult_q(poly p, poly q, const ring r)
993{
994  if (p == NULL)
995  {
996    r->p_Procs->p_Delete(&q, r);
997    return NULL;
998  }
999  if (q == NULL)
1000  {
1001    r->p_Procs->p_Delete(&p, r);
1002    return NULL;
1003  }
1004
1005  if (pNext(p) == NULL)
1006  {
1007#ifdef HAVE_PLURAL
1008    if (rIsPluralRing(r))
1009      q = nc_mm_Mult_p(p, q, r);
1010    else
1011#endif /* HAVE_PLURAL */
1012      q = r->p_Procs->p_Mult_mm(q, p, r);
1013
1014    r->p_Procs->p_Delete(&p, r);
1015    return q;
1016  }
1017
1018  if (pNext(q) == NULL)
1019  {
1020  // NEEDED
1021#ifdef HAVE_PLURAL
1022/*    if (rIsPluralRing(r))
1023      p = gnc_p_Mult_mm(p, q, r); // ???
1024    else*/
1025#endif /* HAVE_PLURAL */
1026      p = r->p_Procs->p_Mult_mm(p, q, r);
1027
1028    r->p_Procs->p_Delete(&q, r);
1029    return p;
1030  }
1031#ifdef HAVE_PLURAL
1032  if (rIsPluralRing(r))
1033    return _nc_p_Mult_q(p, q, r);
1034  else
1035#endif
1036  return _p_Mult_q(p, q, 0, r);
1037}
1038
1039// returns p*q, does neither destroy p nor q
1040static inline poly pp_Mult_qq(poly p, poly q, const ring r)
1041{
1042  poly last;
1043  if (p == NULL || q == NULL) return NULL;
1044
1045  if (pNext(p) == NULL)
1046  {
1047#ifdef HAVE_PLURAL
1048    if (rIsPluralRing(r))
1049      return nc_mm_Mult_pp(p, q, r);
1050#endif
1051    return r->p_Procs->pp_Mult_mm(q, p, r, last);
1052  }
1053
1054  if (pNext(q) == NULL)
1055  {
1056    return r->p_Procs->pp_Mult_mm(p, q, r, last);
1057  }
1058
1059  poly qq = q;
1060  if (p == q)
1061    qq = p_Copy(q, r);
1062
1063  poly res;
1064#ifdef HAVE_PLURAL
1065  if (rIsPluralRing(r))
1066    res = _nc_pp_Mult_qq(p, qq, r);
1067  else
1068#endif
1069    res = _p_Mult_q(p, qq, 1, r);
1070
1071  if (qq != q)
1072    p_Delete(&qq, r);
1073  return res;
1074}
1075
1076// returns p + m*q destroys p, const: q, m
1077static inline poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, int &lp, int lq,
1078                                const ring r)
1079{
1080#ifdef HAVE_PLURAL
1081  if (rIsPluralRing(r))
1082    return nc_p_Plus_mm_Mult_qq(p, m, q, lp, lq, r);
1083#endif
1084
1085// this should be implemented more efficiently
1086  poly res, last;
1087  int shorter;
1088  number n_old = pGetCoeff(m);
1089  number n_neg = n_Copy(n_old, r->cf);
1090  n_neg = n_Neg(n_neg, r->cf);
1091  pSetCoeff0(m, n_neg);
1092  res = r->p_Procs->p_Minus_mm_Mult_qq(p, m, q, shorter, NULL, r, last);
1093  lp = (lp + lq) - shorter;
1094  pSetCoeff0(m, n_old);
1095  n_Delete(&n_neg, r->cf);
1096  return res;
1097}
1098
1099static inline poly p_Plus_mm_Mult_qq(poly p, poly m, poly q, const ring r)
1100{
1101  int lp = 0, lq = 0;
1102  return p_Plus_mm_Mult_qq(p, m, q, lp, lq, r);
1103}
1104
1105// returns merged p and q, assumes p and q have no monomials which are equal
1106static inline poly p_Merge_q(poly p, poly q, const ring r)
1107{
1108  return r->p_Procs->p_Merge_q(p, q, r);
1109}
1110
1111// like p_SortMerge, except that p may have equal monimals
1112static inline poly p_SortAdd(poly p, const ring r, BOOLEAN revert= FALSE)
1113{
1114  if (revert) p = pReverse(p);
1115  return sBucketSortAdd(p, r);
1116}
1117
1118// sorts p using bucket sort: returns sorted poly
1119// assumes that monomials of p are all different
1120// reverses it first, if revert == TRUE, use this if input p is "almost" sorted
1121// correctly
1122static inline poly p_SortMerge(poly p, const ring r, BOOLEAN revert= FALSE)
1123{
1124  if (revert) p = pReverse(p);
1125  return sBucketSortMerge(p, r);
1126}
1127
1128/***************************************************************
1129 *
1130 * I/O
1131 *
1132 ***************************************************************/
1133static inline char*     p_String(poly p, ring p_ring)
1134{
1135  return p_String(p, p_ring, p_ring);
1136}
1137static inline char*     p_String0(poly p, ring p_ring)
1138{
1139  return p_String0(p, p_ring, p_ring);
1140}
1141static inline void      p_Write(poly p, ring p_ring)
1142{
1143  p_Write(p, p_ring, p_ring);
1144}
1145static inline void      p_Write0(poly p, ring p_ring)
1146{
1147  p_Write0(p, p_ring, p_ring);
1148}
1149static inline void      p_wrp(poly p, ring p_ring)
1150{
1151  p_wrp(p, p_ring, p_ring);
1152}
1153
1154/***************************************************************
1155 *  Purpose: implementation of poly procs which iter over ExpVector
1156 *  Author:  obachman (Olaf Bachmann)
1157 *  Created: 8/00
1158 *  Version: $Id$
1159 *******************************************************************/
1160#include <misc/mylimits.h>
1161#include <polys/templates/p_MemCmp.h>
1162// #include <polys/structs.h>
1163#include <polys/monomials/ring.h>
1164#include <coeffs/coeffs.h>
1165
1166#if PDEBUG > 0
1167
1168#define _p_LmCmpAction(p, q, r, actionE, actionG, actionS)  \
1169do                                                          \
1170{                                                           \
1171  int _cmp = p_LmCmp(p,q,r);                                \
1172  if (_cmp == 0) actionE;                                   \
1173  if (_cmp == 1) actionG;                                   \
1174  actionS;                                                  \
1175}                                                           \
1176while(0)
1177
1178#else
1179
1180#define _p_LmCmpAction(p, q, r, actionE, actionG, actionS)                      \
1181 p_MemCmp_LengthGeneral_OrdGeneral(p->exp, q->exp, r->CmpL_Size, r->ordsgn,    \
1182                                   actionE, actionG, actionS)
1183
1184#endif
1185
1186#define pDivAssume(x)   ((void)0)
1187
1188#include <omalloc/omalloc.h>
1189#include <coeffs/coeffs.h>
1190#include <polys/monomials/p_polys.h>
1191#include <polys/templates/p_MemAdd.h>
1192#include <polys/templates/p_MemCopy.h>
1193
1194/***************************************************************
1195 *
1196 * Allocation/Initalization/Deletion
1197 *
1198 ***************************************************************/
1199// adjustments for negative weights
1200static inline void p_MemAdd_NegWeightAdjust(poly p, const ring r)
1201{
1202  if (r->NegWeightL_Offset != NULL)
1203  {
1204    for (int i=r->NegWeightL_Size-1; i>=0; i--)
1205    {
1206      p->exp[r->NegWeightL_Offset[i]] -= POLY_NEGWEIGHT_OFFSET;
1207    }
1208  }
1209}
1210static inline void p_MemSub_NegWeightAdjust(poly p, const ring r)
1211{
1212  if (r->NegWeightL_Offset != NULL)
1213  {
1214    for (int i=r->NegWeightL_Size-1; i>=0; i--)
1215    {
1216      p->exp[r->NegWeightL_Offset[i]] += POLY_NEGWEIGHT_OFFSET;
1217    }
1218  }
1219}
1220// ExpVextor(d_p) = ExpVector(s_p)
1221static inline void p_ExpVectorCopy(poly d_p, poly s_p, const ring r)
1222{
1223  p_LmCheckPolyRing1(d_p, r);
1224  p_LmCheckPolyRing1(s_p, r);
1225  p_MemCopy_LengthGeneral(d_p->exp, s_p->exp, r->ExpL_Size);
1226}
1227
1228static inline poly p_Init(const ring r, omBin bin)
1229{
1230  p_CheckRing1(r);
1231  pAssume1(bin != NULL && r->PolyBin->sizeW == bin->sizeW);
1232  poly p;
1233  omTypeAlloc0Bin(poly, p, bin);
1234  p_MemAdd_NegWeightAdjust(p, r);
1235  p_SetRingOfLm(p, r);
1236  return p;
1237}
1238static inline poly p_Init(const ring r)
1239{
1240  return p_Init(r, r->PolyBin);
1241}
1242
1243static inline poly p_LmInit(poly p, const ring r)
1244{
1245  p_LmCheckPolyRing1(p, r);
1246  poly np;
1247  omTypeAllocBin(poly, np, r->PolyBin);
1248  p_SetRingOfLm(np, r);
1249  p_MemCopy_LengthGeneral(np->exp, p->exp, r->ExpL_Size);
1250  pNext(np) = NULL;
1251  pSetCoeff0(np, NULL);
1252  return np;
1253}
1254static inline poly p_LmInit(poly s_p, const ring s_r, const ring d_r, omBin d_bin)
1255{
1256  p_LmCheckPolyRing1(s_p, s_r);
1257  p_CheckRing(d_r);
1258  pAssume1(d_r->N <= s_r->N);
1259  poly d_p = p_Init(d_r, d_bin);
1260  for (int i=d_r->N; i>0; i--)
1261  {
1262    p_SetExp(d_p, i, p_GetExp(s_p, i,s_r), d_r);
1263  }
1264  if (rRing_has_Comp(d_r))
1265  {
1266    p_SetComp(d_p, p_GetComp(s_p,s_r), d_r);
1267  }
1268  p_Setm(d_p, d_r);
1269  return d_p;
1270}
1271static inline poly p_LmInit(poly s_p, const ring s_r, const ring d_r)
1272{
1273  pAssume1(d_r != NULL);
1274  return p_LmInit(s_p, s_r, d_r, d_r->PolyBin);
1275}
1276
1277// set all exponents l..k to 0, assume exp. k+1..n and 1..l-1 are in
1278// different blocks
1279// set coeff to 1
1280static inline poly p_GetExp_k_n(poly p, int l, int k, const ring r)
1281{
1282  if (p == NULL) return NULL;
1283  p_LmCheckPolyRing1(p, r);
1284  poly np;
1285  omTypeAllocBin(poly, np, r->PolyBin);
1286  p_SetRingOfLm(np, r);
1287  p_MemCopy_LengthGeneral(np->exp, p->exp, r->ExpL_Size);
1288  pNext(np) = NULL;
1289  pSetCoeff0(np, n_Init(1, r->cf));
1290  int i;
1291  for(i=l;i<=k;i++)
1292  {
1293    //np->exp[(r->VarOffset[i] & 0xffffff)] =0;
1294    p_SetExp(np,i,0,r);
1295  }
1296  p_Setm(np,r);
1297  return np;
1298}
1299
1300// simialar to p_ShallowCopyDelete but does it only for leading monomial
1301static inline poly p_LmShallowCopyDelete(poly p, const ring r, omBin bin)
1302{
1303  p_LmCheckPolyRing1(p, r);
1304  pAssume1(bin->sizeW == r->PolyBin->sizeW);
1305  poly new_p = p_New(r);
1306  p_MemCopy_LengthGeneral(new_p->exp, p->exp, r->ExpL_Size);
1307  pSetCoeff0(new_p, pGetCoeff(p));
1308  pNext(new_p) = pNext(p);
1309  omFreeBinAddr(p);
1310  return new_p;
1311}
1312
1313/***************************************************************
1314 *
1315 * Operation on ExpVectors
1316 *
1317 ***************************************************************/
1318// ExpVector(p1) += ExpVector(p2)
1319static inline void p_ExpVectorAdd(poly p1, poly p2, const ring r)
1320{
1321  p_LmCheckPolyRing1(p1, r);
1322  p_LmCheckPolyRing1(p2, r);
1323#if PDEBUG >= 1
1324  for (int i=1; i<=r->N; i++)
1325    pAssume1((unsigned long) (p_GetExp(p1, i, r) + p_GetExp(p2, i, r)) <= r->bitmask);
1326  pAssume1(p_GetComp(p1, r) == 0 || p_GetComp(p2, r) == 0);
1327#endif
1328
1329  p_MemAdd_LengthGeneral(p1->exp, p2->exp, r->ExpL_Size);
1330  p_MemAdd_NegWeightAdjust(p1, r);
1331}
1332// ExpVector(p1) -= ExpVector(p2)
1333static inline void p_ExpVectorSub(poly p1, poly p2, const ring r)
1334{
1335  p_LmCheckPolyRing1(p1, r);
1336  p_LmCheckPolyRing1(p2, r);
1337#if PDEBUG >= 1
1338  for (int i=1; i<=r->N; i++)
1339    pAssume1(p_GetExp(p1, i, r) >= p_GetExp(p2, i, r));
1340  pAssume1(p_GetComp(p1, r) == 0 || p_GetComp(p2, r) == 0 ||
1341          p_GetComp(p1, r) == p_GetComp(p2, r));
1342#endif
1343
1344  p_MemSub_LengthGeneral(p1->exp, p2->exp, r->ExpL_Size);
1345  p_MemSub_NegWeightAdjust(p1, r);
1346
1347}
1348// ExpVector(p1) += ExpVector(p2) - ExpVector(p3)
1349static inline void p_ExpVectorAddSub(poly p1, poly p2, poly p3, const ring r)
1350{
1351  p_LmCheckPolyRing1(p1, r);
1352  p_LmCheckPolyRing1(p2, r);
1353  p_LmCheckPolyRing1(p3, r);
1354#if PDEBUG >= 1
1355  for (int i=1; i<=r->N; i++)
1356    pAssume1(p_GetExp(p1, i, r) + p_GetExp(p2, i, r) >= p_GetExp(p3, i, r));
1357  pAssume1(p_GetComp(p1, r) == 0 ||
1358           (p_GetComp(p2, r) - p_GetComp(p3, r) == 0) ||
1359           (p_GetComp(p1, r) == p_GetComp(p2, r) - p_GetComp(p3, r)));
1360#endif
1361
1362  p_MemAddSub_LengthGeneral(p1->exp, p2->exp, p3->exp, r->ExpL_Size);
1363  // no need to adjust in case of NegWeights
1364}
1365
1366// ExpVector(pr) = ExpVector(p1) + ExpVector(p2)
1367static inline void p_ExpVectorSum(poly pr, poly p1, poly p2, const ring r)
1368{
1369  p_LmCheckPolyRing1(p1, r);
1370  p_LmCheckPolyRing1(p2, r);
1371  p_LmCheckPolyRing1(pr, r);
1372#if PDEBUG >= 1
1373  for (int i=1; i<=r->N; i++)
1374    pAssume1((unsigned long) (p_GetExp(p1, i, r) + p_GetExp(p2, i, r)) <= r->bitmask);
1375  pAssume1(p_GetComp(p1, r) == 0 || p_GetComp(p2, r) == 0);
1376#endif
1377
1378  p_MemSum_LengthGeneral(pr->exp, p1->exp, p2->exp, r->ExpL_Size);
1379  p_MemAdd_NegWeightAdjust(pr, r);
1380}
1381// ExpVector(pr) = ExpVector(p1) - ExpVector(p2)
1382static inline void p_ExpVectorDiff(poly pr, poly p1, poly p2, const ring r)
1383{
1384  p_LmCheckPolyRing1(p1, r);
1385  p_LmCheckPolyRing1(p2, r);
1386  p_LmCheckPolyRing1(pr, r);
1387#if PDEBUG >= 2
1388  for (int i=1; i<=r->N; i++)
1389    pAssume1(p_GetExp(p1, i, r) >= p_GetExp(p2, i, r));
1390  pAssume1(!rRing_has_Comp(r) || p_GetComp(p1, r) == p_GetComp(p2, r));
1391#endif
1392
1393  p_MemDiff_LengthGeneral(pr->exp, p1->exp, p2->exp, r->ExpL_Size);
1394  p_MemSub_NegWeightAdjust(pr, r);
1395}
1396
1397static inline BOOLEAN p_ExpVectorEqual(poly p1, poly p2, const ring r)
1398{
1399  p_LmCheckPolyRing1(p1, r);
1400  p_LmCheckPolyRing1(p2, r);
1401
1402  int i = r->ExpL_Size;
1403  unsigned long *ep = p1->exp;
1404  unsigned long *eq = p2->exp;
1405
1406  do
1407  {
1408    i--;
1409    if (ep[i] != eq[i]) return FALSE;
1410  }
1411  while (i);
1412  return TRUE;
1413}
1414
1415static inline long p_Totaldegree(poly p, const ring r)
1416{
1417  p_LmCheckPolyRing1(p, r);
1418  unsigned long s = p_GetTotalDegree(p->exp[r->VarL_Offset[0]],
1419                                     r,
1420                                     r->MinExpPerLong);
1421  for (int i=r->VarL_Size-1; i>0; i--)
1422  {
1423    s += p_GetTotalDegree(p->exp[r->VarL_Offset[i]], r);
1424  }
1425  return (long)s;
1426}
1427
1428static inline void p_GetExpV(poly p, int *ev, const ring r)
1429{
1430  p_LmCheckPolyRing1(p, r);
1431  for (int j = r->N; j; j--)
1432      ev[j] = p_GetExp(p, j, r);
1433
1434  ev[0] = p_GetComp(p, r);
1435}
1436static inline void p_SetExpV(poly p, int *ev, const ring r)
1437{
1438  p_LmCheckPolyRing1(p, r);
1439  for (int j = r->N; j; j--)
1440      p_SetExp(p, j, ev[j], r);
1441
1442  p_SetComp(p, ev[0],r);
1443  p_Setm(p, r);
1444}
1445
1446/***************************************************************
1447 *
1448 * Comparison w.r.t. monomial ordering
1449 *
1450 ***************************************************************/
1451static inline int p_LmCmp(poly p, poly q, const ring r)
1452{
1453  p_LmCheckPolyRing1(p, r);
1454  p_LmCheckPolyRing1(q, r);
1455
1456  p_MemCmp_LengthGeneral_OrdGeneral(p->exp, q->exp, r->CmpL_Size, r->ordsgn,
1457                                    return 0, return 1, return -1);
1458}
1459
1460
1461
1462/***************************************************************
1463 *
1464 * Comparisons: they are all done without regarding coeffs
1465 *
1466 ***************************************************************/
1467#define p_LmCmpAction(p, q, r, actionE, actionG, actionS) \
1468  _p_LmCmpAction(p, q, r, actionE, actionG, actionS)
1469
1470// returns 1 if ExpVector(p)==ExpVector(q): does not compare numbers !!
1471#define p_LmEqual(p1, p2, r) p_ExpVectorEqual(p1, p2, r)
1472
1473// pCmp: args may be NULL
1474// returns: (p2==NULL ? 1 : (p1 == NULL ? -1 : p_LmCmp(p1, p2)))
1475static inline int p_Cmp(poly p1, poly p2, ring r)
1476{
1477  if (p2==NULL)
1478    return 1;
1479  if (p1==NULL)
1480    return -1;
1481  return p_LmCmp(p1,p2,r);
1482}
1483
1484
1485/***************************************************************
1486 *
1487 * divisibility
1488 *
1489 ***************************************************************/
1490// return: FALSE, if there exists i, such that a->exp[i] > b->exp[i]
1491//         TRUE, otherwise
1492// (1) Consider long vars, instead of single exponents
1493// (2) Clearly, if la > lb, then FALSE
1494// (3) Suppose la <= lb, and consider first bits of single exponents in l:
1495//     if TRUE, then value of these bits is la ^ lb
1496//     if FALSE, then la-lb causes an "overflow" into one of those bits, i.e.,
1497//               la ^ lb != la - lb
1498static inline BOOLEAN _p_LmDivisibleByNoComp(poly a, poly b, const ring r)
1499{
1500  int i=r->VarL_Size - 1;
1501  unsigned long divmask = r->divmask;
1502  unsigned long la, lb;
1503
1504  if (r->VarL_LowIndex >= 0)
1505  {
1506    i += r->VarL_LowIndex;
1507    do
1508    {
1509      la = a->exp[i];
1510      lb = b->exp[i];
1511      if ((la > lb) ||
1512          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
1513      {
1514        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
1515        return FALSE;
1516      }
1517      i--;
1518    }
1519    while (i>=r->VarL_LowIndex);
1520  }
1521  else
1522  {
1523    do
1524    {
1525      la = a->exp[r->VarL_Offset[i]];
1526      lb = b->exp[r->VarL_Offset[i]];
1527      if ((la > lb) ||
1528          (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask)))
1529      {
1530        pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == FALSE);
1531        return FALSE;
1532      }
1533      i--;
1534    }
1535    while (i>=0);
1536  }
1537#ifdef HAVE_RINGS
1538  pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == nDivBy(p_GetCoeff(b, r), p_GetCoeff(a, r)));
1539  return (!rField_is_Ring(r)) || n_DivBy(p_GetCoeff(b, r), p_GetCoeff(a, r), r->cf);
1540#else
1541  pDivAssume(p_DebugLmDivisibleByNoComp(a, b, r) == TRUE);
1542  return TRUE;
1543#endif
1544}
1545
1546static inline BOOLEAN _p_LmDivisibleByNoComp(poly a, const ring r_a, poly b, const ring r_b)
1547{
1548  int i=r_a->N;
1549  pAssume1(r_a->N == r_b->N);
1550
1551  do
1552  {
1553    if (p_GetExp(a,i,r_a) > p_GetExp(b,i,r_b))
1554      return FALSE;
1555    i--;
1556  }
1557  while (i);
1558#ifdef HAVE_RINGS
1559  return n_DivBy(p_GetCoeff(b, r_b), p_GetCoeff(a, r_a), r_a->cf);
1560#else
1561  return TRUE;
1562#endif
1563}
1564
1565#ifdef HAVE_RATGRING
1566static inline BOOLEAN _p_LmDivisibleByNoCompPart(poly a, const ring r_a, poly b, const ring r_b,const int start, const int end)
1567{
1568  int i=end;
1569  pAssume1(r_a->N == r_b->N);
1570
1571  do
1572  {
1573    if (p_GetExp(a,i,r_a) > p_GetExp(b,i,r_b))
1574      return FALSE;
1575    i--;
1576  }
1577  while (i>=start);
1578#ifdef HAVE_RINGS
1579  return nDivBy(p_GetCoeff(b, r), p_GetCoeff(a, r));
1580#else
1581  return TRUE;
1582#endif
1583}
1584static inline BOOLEAN _p_LmDivisibleByPart(poly a, const ring r_a, poly b, const ring r_b,const int start, const int end)
1585{
1586  if (p_GetComp(a, r_a) == 0 || p_GetComp(a,r_a) == p_GetComp(b,r_b))
1587    return _p_LmDivisibleByNoCompPart(a, r_a, b, r_b,start,end);
1588  return FALSE;
1589}
1590static inline BOOLEAN p_LmDivisibleByPart(poly a, poly b, const ring r,const int start, const int end)
1591{
1592  p_LmCheckPolyRing1(b, r);
1593  pIfThen1(a != NULL, p_LmCheckPolyRing1(b, r));
1594  if (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r))
1595    return _p_LmDivisibleByNoCompPart(a, r, b, r,start, end);
1596  return FALSE;
1597}
1598#endif
1599static inline BOOLEAN _p_LmDivisibleBy(poly a, poly b, const ring r)
1600{
1601  if (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r))
1602    return _p_LmDivisibleByNoComp(a, b, r);
1603  return FALSE;
1604}
1605static inline BOOLEAN _p_LmDivisibleBy(poly a, const ring r_a, poly b, const ring r_b)
1606{
1607  if (p_GetComp(a, r_a) == 0 || p_GetComp(a,r_a) == p_GetComp(b,r_b))
1608    return _p_LmDivisibleByNoComp(a, r_a, b, r_b);
1609  return FALSE;
1610}
1611static inline BOOLEAN p_LmDivisibleByNoComp(poly a, poly b, const ring r)
1612{
1613  p_LmCheckPolyRing1(a, r);
1614  p_LmCheckPolyRing1(b, r);
1615  return _p_LmDivisibleByNoComp(a, b, r);
1616}
1617static inline BOOLEAN p_LmDivisibleBy(poly a, poly b, const ring r)
1618{
1619  p_LmCheckPolyRing1(b, r);
1620  pIfThen1(a != NULL, p_LmCheckPolyRing1(b, r));
1621  if (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r))
1622    return _p_LmDivisibleByNoComp(a, b, r);
1623  return FALSE;
1624}
1625
1626static inline BOOLEAN p_DivisibleBy(poly a, poly b, const ring r)
1627{
1628  pIfThen1(b!=NULL, p_LmCheckPolyRing1(b, r));
1629  pIfThen1(a!=NULL, p_LmCheckPolyRing1(a, r));
1630
1631  if (a != NULL && (p_GetComp(a, r) == 0 || p_GetComp(a,r) == p_GetComp(b,r)))
1632      return _p_LmDivisibleByNoComp(a,b,r);
1633  return FALSE;
1634}
1635static inline BOOLEAN p_DivisibleBy(poly a, const ring r_a, poly b, const ring r_b)
1636{
1637  pIfThen1(b!=NULL, p_LmCheckPolyRing1(b, r_b));
1638  pIfThen1(a!=NULL, p_LmCheckPolyRing1(a, r_a));
1639  if (a != NULL) {
1640      return _p_LmDivisibleBy(a, r_a, b, r_b);
1641  }
1642  return FALSE;
1643}
1644static inline BOOLEAN p_LmDivisibleBy(poly a, const ring r_a, poly b, const ring r_b)
1645{
1646  p_LmCheckPolyRing(a, r_a);
1647  p_LmCheckPolyRing(b, r_b);
1648  return _p_LmDivisibleBy(a, r_a, b, r_b);
1649}
1650static inline BOOLEAN p_LmShortDivisibleBy(poly a, unsigned long sev_a,
1651                                    poly b, unsigned long not_sev_b, const ring r)
1652{
1653  p_LmCheckPolyRing1(a, r);
1654  p_LmCheckPolyRing1(b, r);
1655#ifndef PDIV_DEBUG
1656  _pPolyAssume2(p_GetShortExpVector(a, r) == sev_a, a, r);
1657  _pPolyAssume2(p_GetShortExpVector(b, r) == ~ not_sev_b, b, r);
1658
1659  if (sev_a & not_sev_b)
1660  {
1661    pAssume1(p_LmDivisibleByNoComp(a, b, r) == FALSE);
1662    return FALSE;
1663  }
1664  return p_LmDivisibleBy(a, b, r);
1665#else
1666  return pDebugLmShortDivisibleBy(a, sev_a, r, b, not_sev_b, r);
1667#endif
1668}
1669
1670static inline BOOLEAN p_LmShortDivisibleBy(poly a, unsigned long sev_a, const ring r_a,
1671                                      poly b, unsigned long not_sev_b, const ring r_b)
1672{
1673  p_LmCheckPolyRing1(a, r_a);
1674  p_LmCheckPolyRing1(b, r_b);
1675#ifndef PDIV_DEBUG
1676  _pPolyAssume2(p_GetShortExpVector(a, r_a) == sev_a, a, r_a);
1677  _pPolyAssume2(p_GetShortExpVector(b, r_b) == ~ not_sev_b, b, r_b);
1678
1679  if (sev_a & not_sev_b)
1680  {
1681    pAssume1(_p_LmDivisibleByNoComp(a, r_a, b, r_b) == FALSE);
1682    return FALSE;
1683  }
1684  return _p_LmDivisibleBy(a, r_a, b, r_b);
1685#else
1686  return pDebugLmShortDivisibleBy(a, sev_a, r_a, b, not_sev_b, r_b);
1687#endif
1688}
1689
1690/***************************************************************
1691 *
1692 * Misc things on Lm
1693 *
1694 ***************************************************************/
1695
1696
1697// like the respective p_LmIs* routines, except that p might be empty
1698static inline BOOLEAN p_IsConstantComp(const poly p, const ring r)
1699{
1700  if (p == NULL) return TRUE;
1701  return (pNext(p)==NULL) && p_LmIsConstantComp(p, r);
1702}
1703
1704static inline BOOLEAN p_IsConstant(const poly p, const ring r)
1705{
1706  if (p == NULL) return TRUE;
1707  return (pNext(p)==NULL) && p_LmIsConstant(p, r);
1708}
1709
1710static inline BOOLEAN p_IsConstantPoly(const poly p, const ring r)
1711{
1712  poly pp=p;
1713  while(pp!=NULL)
1714  {
1715    if (! p_LmIsConstantComp(pp, r))
1716      return FALSE;
1717    pIter(pp);
1718  }
1719  return TRUE;
1720}
1721
1722static inline BOOLEAN p_IsUnit(const poly p, const ring r)
1723{
1724  if (p == NULL) return FALSE;
1725#ifdef HAVE_RINGS
1726  if (rField_is_Ring(r))
1727    return (p_LmIsConstant(p, r) && nIsUnit(pGetCoeff(p),r->cf));
1728#endif
1729  return p_LmIsConstant(p, r);
1730}
1731
1732static inline BOOLEAN p_LmExpVectorAddIsOk(const poly p1, const poly p2,
1733                                      const ring r)
1734{
1735  p_LmCheckPolyRing(p1, r);
1736  p_LmCheckPolyRing(p2, r);
1737  unsigned long l1, l2, divmask = r->divmask;
1738  int i;
1739
1740  for (i=0; i<r->VarL_Size; i++)
1741  {
1742    l1 = p1->exp[r->VarL_Offset[i]];
1743    l2 = p2->exp[r->VarL_Offset[i]];
1744    // do the divisiblity trick
1745    if ( (l1 > ULONG_MAX - l2) ||
1746         (((l1 & divmask) ^ (l2 & divmask)) != ((l1 + l2) & divmask)))
1747      return FALSE;
1748  }
1749  return TRUE;
1750}
1751void      p_Split(poly p, poly * r);   /*p => IN(p), r => REST(p) */
1752BOOLEAN p_HasNotCF(poly p1, poly p2, const ring r);
1753poly      p_mInit(const char *s, BOOLEAN &ok, const ring r); /* monom s -> poly, interpreter */
1754const char *    p_Read(const char *s, poly &p,const ring r); /* monom -> poly */
1755poly      p_Divide(poly a, poly b, const ring r);
1756poly      p_DivideM(poly a, poly b, const ring r);
1757void      p_Lcm(poly a, poly b, poly m, const ring r);
1758poly      p_Diff(poly a, int k, const ring r);
1759poly      p_DiffOp(poly a, poly b,BOOLEAN multiply, const ring r);
1760int       p_Weight(int c, const ring r);
1761
1762/* syszygy stuff */
1763BOOLEAN   p_VectorHasUnitB(poly p, int * k, const ring r);
1764void      p_VectorHasUnit(poly p, int * k, int * len, const ring r);
1765poly      p_TakeOutComp1(poly * p, int k, const ring r);
1766// Splits *p into two polys: *q which consists of all monoms with
1767// component == comp and *p of all other monoms *lq == pLength(*q)
1768// On return all components pf *q == 0
1769void p_TakeOutComp(poly *p, long comp, poly *q, int *lq, const ring r);
1770
1771// This is something weird -- Don't use it, unless you know what you are doing
1772poly      p_TakeOutComp(poly * p, int k);
1773
1774void      p_DeleteComp(poly * p,int k, const ring r);
1775
1776/*-------------ring management:----------------------*/
1777void p_SetGlobals(const ring r, BOOLEAN complete = TRUE);
1778
1779// resets the pFDeg and pLDeg: if pLDeg is not given, it is
1780// set to currRing->pLDegOrig, i.e. to the respective LDegProc which
1781// only uses pFDeg (and not pDeg, or pTotalDegree, etc).
1782// If you use this, make sure your procs does not make any assumptions
1783// on ordering and/or OrdIndex -- otherwise they might return wrong results
1784// on strat->tailRing
1785void pSetDegProcs(ring r, pFDegProc new_FDeg, pLDegProc new_lDeg = NULL);
1786// restores pFDeg and pLDeg:
1787void pRestoreDegProcs(ring r, pFDegProc old_FDeg, pLDegProc old_lDeg);
1788
1789/*-------------pComp for syzygies:-------------------*/
1790void p_SetModDeg(intvec *w, ring r);
1791
1792/*------------ Jet ----------------------------------*/
1793poly pp_Jet(poly p, int m, const ring R);
1794poly p_Jet(poly p, int m,const ring R);
1795poly pp_JetW(poly p, int m, short *w, const ring R);
1796poly p_JetW(poly p, int m, short *w, const ring R);
1797
1798
1799poly      p_PermPoly (poly p, int * perm,const ring OldRing, const ring dst,
1800                     nMapFunc nMap, int *par_perm=NULL, int OldPar=0);
1801
1802/*----------------------------------------------------*/
1803poly p_Series(int n,poly p,poly u, intvec *w, const ring R);
1804poly p_Invers(int n,poly u,intvec *w, const ring R);
1805
1806
1807#endif // P_POLYS_H
1808
Note: See TracBrowser for help on using the repository browser.