source: git/libpolys/polys/monomials/p_polys.h @ 304ad9b

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