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

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