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

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