source: git/libpolys/polys/monomials/p_polys.h @ 805d0b1

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