source: git/libpolys/polys/monomials/p_polys.h @ 4f0f42

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