source: git/libpolys/polys/monomials/p_polys.h @ 5948a8

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