# source:git/polys/monomials/p_polys.h@a04c5e

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