source: git/libpolys/polys/monomials/p_polys.h @ ae030d

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