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

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