# source:git/libpolys/polys/monomials/p_polys.h@71ba5b8

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