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