source: git/libpolys/polys/monomials/p_polys.h @ 8e45403

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