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

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