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

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