source: git/libpolys/polys/monomials/p_polys.h @ 48ca29

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