# source:git/libpolys/polys/monomials/p_polys.h@5a3ae8

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