1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /***************************************************************
5 * File: p_polys.cc
6 * Purpose: implementation of ring independent poly procedures?
7 * Author: obachman (Olaf Bachmann)
8 * Created: 8/00
9 *******************************************************************/
10
11 #include <ctype.h>
12
13 #include "misc/auxiliary.h"
14
15 #include "misc/options.h"
16 #include "misc/intvec.h"
17
18
19 #include "coeffs/longrat.h" // snumber is needed...
20 #include "coeffs/numbers.h" // ndCopyMap
21
22 #include "polys/PolyEnumerator.h"
23
24 #define TRANSEXT_PRIVATES
25
26 #include "polys/ext_fields/transext.h"
27 #include "polys/ext_fields/algext.h"
28
29 #include "polys/weight.h"
30 #include "polys/simpleideals.h"
31
32 #include "ring.h"
33 #include "p_polys.h"
34
35 #include "polys/templates/p_MemCmp.h"
36 #include "polys/templates/p_MemAdd.h"
37 #include "polys/templates/p_MemCopy.h"
38
39
40 #ifdef HAVE_PLURAL
41 #include "nc/nc.h"
42 #include "nc/sca.h"
43 #endif
44
45 #ifdef HAVE_SHIFTBBA
46 #include "polys/shiftop.h"
47 #endif
48
49 #include "clapsing.h"
50
51 /*
52 * lift ideal with coeffs over Z (mod N) to Q via Farey
53 */
p_Farey(poly p,number N,const ring r)54 poly p_Farey(poly p, number N, const ring r)
55 {
56 poly h=p_Copy(p,r);
57 poly hh=h;
58 while(h!=NULL)
59 {
60 number c=pGetCoeff(h);
61 pSetCoeff0(h,n_Farey(c,N,r->cf));
62 n_Delete(&c,r->cf);
63 pIter(h);
64 }
65 while((hh!=NULL)&&(n_IsZero(pGetCoeff(hh),r->cf)))
66 {
67 p_LmDelete(&hh,r);
68 }
69 h=hh;
70 while((h!=NULL) && (pNext(h)!=NULL))
71 {
72 if(n_IsZero(pGetCoeff(pNext(h)),r->cf))
73 {
74 p_LmDelete(&pNext(h),r);
75 }
76 else pIter(h);
77 }
78 return hh;
79 }
80 /*2
81 * xx,q: arrays of length 0..rl-1
82 * xx[i]: SB mod q[i]
83 * assume: char=0
84 * assume: q[i]!=0
85 * x: work space
86 * destroys xx
87 */
p_ChineseRemainder(poly * xx,number * x,number * q,int rl,CFArray & inv_cache,const ring R)88 poly p_ChineseRemainder(poly *xx, number *x,number *q, int rl, CFArray &inv_cache,const ring R)
89 {
90 poly r,h,hh;
91 int j;
92 poly res_p=NULL;
93 loop
94 {
95 /* search the lead term */
96 r=NULL;
97 for(j=rl-1;j>=0;j--)
98 {
99 h=xx[j];
100 if ((h!=NULL)
101 &&((r==NULL)||(p_LmCmp(r,h,R)==-1)))
102 r=h;
103 }
104 /* nothing found -> return */
105 if (r==NULL) break;
106 /* create the monomial in h */
107 h=p_Head(r,R);
108 /* collect the coeffs in x[..]*/
109 for(j=rl-1;j>=0;j--)
110 {
111 hh=xx[j];
112 if ((hh!=NULL) && (p_LmCmp(h,hh,R)==0))
113 {
114 x[j]=pGetCoeff(hh);
115 hh=p_LmFreeAndNext(hh,R);
116 xx[j]=hh;
117 }
118 else
119 x[j]=n_Init(0, R->cf);
120 }
121 number n=n_ChineseRemainderSym(x,q,rl,TRUE,inv_cache,R->cf);
122 for(j=rl-1;j>=0;j--)
123 {
124 x[j]=NULL; // n_Init(0...) takes no memory
125 }
126 if (n_IsZero(n,R->cf)) p_Delete(&h,R);
127 else
128 {
129 //Print("new mon:");pWrite(h);
130 p_SetCoeff(h,n,R);
131 pNext(h)=res_p;
132 res_p=h; // building res_p in reverse order!
133 }
134 }
135 res_p=pReverse(res_p);
136 p_Test(res_p, R);
137 return res_p;
138 }
139 /***************************************************************
140 *
141 * Completing what needs to be set for the monomial
142 *
143 ***************************************************************/
144 // this is special for the syz stuff
145 STATIC_VAR int* _components = NULL;
146 STATIC_VAR long* _componentsShifted = NULL;
147 STATIC_VAR int _componentsExternal = 0;
148
149 VAR BOOLEAN pSetm_error=0;
150
151 #ifndef SING_NDEBUG
152 # define MYTEST 0
153 #else /* ifndef SING_NDEBUG */
154 # define MYTEST 0
155 #endif /* ifndef SING_NDEBUG */
156
p_Setm_General(poly p,const ring r)157 void p_Setm_General(poly p, const ring r)
158 {
159 p_LmCheckPolyRing(p, r);
160 int pos=0;
161 if (r->typ!=NULL)
162 {
163 loop
164 {
165 unsigned long ord=0;
166 sro_ord* o=&(r->typ[pos]);
167 switch(o->ord_typ)
168 {
169 case ro_dp:
170 {
171 int a,e;
172 a=o->data.dp.start;
173 e=o->data.dp.end;
174 for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r);
175 p->exp[o->data.dp.place]=ord;
176 break;
177 }
178 case ro_wp_neg:
179 ord=POLY_NEGWEIGHT_OFFSET;
180 // no break;
181 case ro_wp:
182 {
183 int a,e;
184 a=o->data.wp.start;
185 e=o->data.wp.end;
186 int *w=o->data.wp.weights;
187 #if 1
188 for(int i=a;i<=e;i++) ord+=((unsigned long)p_GetExp(p,i,r))*((unsigned long)w[i-a]);
189 #else
190 long ai;
191 int ei,wi;
192 for(int i=a;i<=e;i++)
193 {
194 ei=p_GetExp(p,i,r);
195 wi=w[i-a];
196 ai=ei*wi;
197 if (ai/ei!=wi) pSetm_error=TRUE;
198 ord+=ai;
199 if (ord<ai) pSetm_error=TRUE;
200 }
201 #endif
202 p->exp[o->data.wp.place]=ord;
203 break;
204 }
205 case ro_am:
206 {
207 ord = POLY_NEGWEIGHT_OFFSET;
208 const short a=o->data.am.start;
209 const short e=o->data.am.end;
210 const int * w=o->data.am.weights;
211 #if 1
212 for(short i=a; i<=e; i++, w++)
213 ord += ((*w) * p_GetExp(p,i,r));
214 #else
215 long ai;
216 int ei,wi;
217 for(short i=a;i<=e;i++)
218 {
219 ei=p_GetExp(p,i,r);
220 wi=w[i-a];
221 ai=ei*wi;
222 if (ai/ei!=wi) pSetm_error=TRUE;
223 ord += ai;
224 if (ord<ai) pSetm_error=TRUE;
225 }
226 #endif
227 const int c = p_GetComp(p,r);
228
229 const short len_gen= o->data.am.len_gen;
230
231 if ((c > 0) && (c <= len_gen))
232 {
233 assume( w == o->data.am.weights_m );
234 assume( w[0] == len_gen );
235 ord += w[c];
236 }
237
238 p->exp[o->data.am.place] = ord;
239 break;
240 }
241 case ro_wp64:
242 {
243 int64 ord=0;
244 int a,e;
245 a=o->data.wp64.start;
246 e=o->data.wp64.end;
247 int64 *w=o->data.wp64.weights64;
248 int64 ei,wi,ai;
249 for(int i=a;i<=e;i++)
250 {
251 //Print("exp %d w %d \n",p_GetExp(p,i,r),(int)w[i-a]);
252 //ord+=((int64)p_GetExp(p,i,r))*w[i-a];
253 ei=(int64)p_GetExp(p,i,r);
254 wi=w[i-a];
255 ai=ei*wi;
256 if(ei!=0 && ai/ei!=wi)
257 {
258 pSetm_error=TRUE;
259 #if SIZEOF_LONG == 4
260 Print("ai %lld, wi %lld\n",ai,wi);
261 #else
262 Print("ai %ld, wi %ld\n",ai,wi);
263 #endif
264 }
265 ord+=ai;
266 if (ord<ai)
267 {
268 pSetm_error=TRUE;
269 #if SIZEOF_LONG == 4
270 Print("ai %lld, ord %lld\n",ai,ord);
271 #else
272 Print("ai %ld, ord %ld\n",ai,ord);
273 #endif
274 }
275 }
276 int64 mask=(int64)0x7fffffff;
277 long a_0=(long)(ord&mask); //2^31
278 long a_1=(long)(ord >>31 ); /*(ord/(mask+1));*/
279
280 //Print("mask: %x, ord: %d, a_0: %d, a_1: %d\n"
281 //,(int)mask,(int)ord,(int)a_0,(int)a_1);
282 //Print("mask: %d",mask);
283
284 p->exp[o->data.wp64.place]=a_1;
285 p->exp[o->data.wp64.place+1]=a_0;
286 // if(p_Setm_error) PrintS("***************************\n"
287 // "***************************\n"
288 // "**WARNING: overflow error**\n"
289 // "***************************\n"
290 // "***************************\n");
291 break;
292 }
293 case ro_cp:
294 {
295 int a,e;
296 a=o->data.cp.start;
297 e=o->data.cp.end;
298 int pl=o->data.cp.place;
299 for(int i=a;i<=e;i++) { p->exp[pl]=p_GetExp(p,i,r); pl++; }
300 break;
301 }
302 case ro_syzcomp:
303 {
304 long c=__p_GetComp(p,r);
305 long sc = c;
306 int* Components = (_componentsExternal ? _components :
307 o->data.syzcomp.Components);
308 long* ShiftedComponents = (_componentsExternal ? _componentsShifted:
309 o->data.syzcomp.ShiftedComponents);
310 if (ShiftedComponents != NULL)
311 {
312 assume(Components != NULL);
313 assume(c == 0 || Components[c] != 0);
314 sc = ShiftedComponents[Components[c]];
315 assume(c == 0 || sc != 0);
316 }
317 p->exp[o->data.syzcomp.place]=sc;
318 break;
319 }
320 case ro_syz:
321 {
322 const unsigned long c = __p_GetComp(p, r);
323 const short place = o->data.syz.place;
324 const int limit = o->data.syz.limit;
325
326 if (c > (unsigned long)limit)
327 p->exp[place] = o->data.syz.curr_index;
328 else if (c > 0)
329 {
330 assume( (1 <= c) && (c <= (unsigned long)limit) );
331 p->exp[place]= o->data.syz.syz_index[c];
332 }
333 else
334 {
335 assume(c == 0);
336 p->exp[place]= 0;
337 }
338 break;
339 }
340 // Prefix for Induced Schreyer ordering
341 case ro_isTemp: // Do nothing?? (to be removed into suffix later on...?)
342 {
343 assume(p != NULL);
344
345 #ifndef SING_NDEBUG
346 #if MYTEST
347 Print("p_Setm_General: ro_isTemp ord: pos: %d, p: ", pos); p_wrp(p, r);
348 #endif
349 #endif
350 int c = p_GetComp(p, r);
351
352 assume( c >= 0 );
353
354 // Let's simulate case ro_syz above....
355 // Should accumulate (by Suffix) and be a level indicator
356 const int* const pVarOffset = o->data.isTemp.pVarOffset;
357
358 assume( pVarOffset != NULL );
359
360 // TODO: Can this be done in the suffix???
361 for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
362 {
363 const int vo = pVarOffset[i];
364 if( vo != -1) // TODO: optimize: can be done once!
365 {
366 // Hans! Please don't break it again! p_SetExp(p, ..., r, vo) is correct:
367 p_SetExp(p, p_GetExp(p, i, r), r, vo); // copy put them verbatim
368 // Hans! Please don't break it again! p_GetExp(p, r, vo) is correct:
369 assume( p_GetExp(p, r, vo) == p_GetExp(p, i, r) ); // copy put them verbatim
370 }
371 }
372 #ifndef SING_NDEBUG
373 for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
374 {
375 const int vo = pVarOffset[i];
376 if( vo != -1) // TODO: optimize: can be done once!
377 {
378 // Hans! Please don't break it again! p_GetExp(p, r, vo) is correct:
379 assume( p_GetExp(p, r, vo) == p_GetExp(p, i, r) ); // copy put them verbatim
380 }
381 }
382 #if MYTEST
383 // if( p->exp[o->data.isTemp.start] > 0 )
384 PrintS("after Values: "); p_wrp(p, r);
385 #endif
386 #endif
387 break;
388 }
389
390 // Suffix for Induced Schreyer ordering
391 case ro_is:
392 {
393 #ifndef SING_NDEBUG
394 #if MYTEST
395 Print("p_Setm_General: ro_is ord: pos: %d, p: ", pos); p_wrp(p, r);
396 #endif
397 #endif
398
399 assume(p != NULL);
400
401 int c = p_GetComp(p, r);
402
403 assume( c >= 0 );
404 const ideal F = o->data.is.F;
405 const int limit = o->data.is.limit;
406 assume( limit >= 0 );
407 const int start = o->data.is.start;
408
409 if( F != NULL && c > limit )
410 {
411 #ifndef SING_NDEBUG
412 #if MYTEST
413 Print("p_Setm_General: ro_is : in rSetm: pos: %d, c: %d > limit: %d\n", c, pos, limit);
414 PrintS("preComputed Values: ");
415 p_wrp(p, r);
416 #endif
417 #endif
418 // if( c > limit ) // BUG???
419 p->exp[start] = 1;
420 // else
421 // p->exp[start] = 0;
422
423
424 c -= limit;
425 assume( c > 0 );
426 c--;
427
428 if( c >= IDELEMS(F) )
429 break;
430
431 assume( c < IDELEMS(F) ); // What about others???
432
433 const poly pp = F->m[c]; // get reference monomial!!!
434
435 if(pp == NULL)
436 break;
437
438 assume(pp != NULL);
439
440 #ifndef SING_NDEBUG
441 #if MYTEST
442 Print("Respective F[c - %d: %d] pp: ", limit, c);
443 p_wrp(pp, r);
444 #endif
445 #endif
446
447 const int end = o->data.is.end;
448 assume(start <= end);
449
450
451 // const int st = o->data.isTemp.start;
452
453 #ifndef SING_NDEBUG
454 #if MYTEST
455 Print("p_Setm_General: is(-Temp-) :: c: %d, limit: %d, [st:%d] ===>>> %ld\n", c, limit, start, p->exp[start]);
456 #endif
457 #endif
458
459 // p_ExpVectorAdd(p, pp, r);
460
461 for( int i = start; i <= end; i++) // v[0] may be here...
462 p->exp[i] += pp->exp[i]; // !!!!!!!! ADD corresponding LT(F)
463
464 // p_MemAddAdjust(p, ri);
465 if (r->NegWeightL_Offset != NULL)
466 {
467 for (int i=r->NegWeightL_Size-1; i>=0; i--)
468 {
469 const int _i = r->NegWeightL_Offset[i];
470 if( start <= _i && _i <= end )
471 p->exp[_i] -= POLY_NEGWEIGHT_OFFSET;
472 }
473 }
474
475
476 #ifndef SING_NDEBUG
477 const int* const pVarOffset = o->data.is.pVarOffset;
478
479 assume( pVarOffset != NULL );
480
481 for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
482 {
483 const int vo = pVarOffset[i];
484 if( vo != -1) // TODO: optimize: can be done once!
485 // Hans! Please don't break it again! p_GetExp(p/pp, r, vo) is correct:
486 assume( p_GetExp(p, r, vo) == (p_GetExp(p, i, r) + p_GetExp(pp, r, vo)) );
487 }
488 // TODO: how to check this for computed values???
489 #if MYTEST
490 PrintS("Computed Values: "); p_wrp(p, r);
491 #endif
492 #endif
493 } else
494 {
495 p->exp[start] = 0; //!!!!????? where?????
496
497 const int* const pVarOffset = o->data.is.pVarOffset;
498
499 // What about v[0] - component: it will be added later by
500 // suffix!!!
501 // TODO: Test it!
502 const int vo = pVarOffset[0];
503 if( vo != -1 )
504 p->exp[vo] = c; // initial component v[0]!
505
506 #ifndef SING_NDEBUG
507 #if MYTEST
508 Print("ELSE p_Setm_General: ro_is :: c: %d <= limit: %d, vo: %d, exp: %d\n", c, limit, vo, p->exp[vo]);
509 p_wrp(p, r);
510 #endif
511 #endif
512 }
513
514 break;
515 }
516 default:
517 dReportError("wrong ord in rSetm:%d\n",o->ord_typ);
518 return;
519 }
520 pos++;
521 if (pos == r->OrdSize) return;
522 }
523 }
524 }
525
p_Setm_Syz(poly p,ring r,int * Components,long * ShiftedComponents)526 void p_Setm_Syz(poly p, ring r, int* Components, long* ShiftedComponents)
527 {
528 _components = Components;
529 _componentsShifted = ShiftedComponents;
530 _componentsExternal = 1;
531 p_Setm_General(p, r);
532 _componentsExternal = 0;
533 }
534
535 // dummy for lp, ls, etc
p_Setm_Dummy(poly p,const ring r)536 void p_Setm_Dummy(poly p, const ring r)
537 {
538 p_LmCheckPolyRing(p, r);
539 }
540
541 // for dp, Dp, ds, etc
p_Setm_TotalDegree(poly p,const ring r)542 void p_Setm_TotalDegree(poly p, const ring r)
543 {
544 p_LmCheckPolyRing(p, r);
545 p->exp[r->pOrdIndex] = p_Totaldegree(p, r);
546 }
547
548 // for wp, Wp, ws, etc
p_Setm_WFirstTotalDegree(poly p,const ring r)549 void p_Setm_WFirstTotalDegree(poly p, const ring r)
550 {
551 p_LmCheckPolyRing(p, r);
552 p->exp[r->pOrdIndex] = p_WFirstTotalDegree(p, r);
553 }
554
p_GetSetmProc(const ring r)555 p_SetmProc p_GetSetmProc(const ring r)
556 {
557 // covers lp, rp, ls,
558 if (r->typ == NULL) return p_Setm_Dummy;
559
560 if (r->OrdSize == 1)
561 {
562 if (r->typ[0].ord_typ == ro_dp &&
563 r->typ[0].data.dp.start == 1 &&
564 r->typ[0].data.dp.end == r->N &&
565 r->typ[0].data.dp.place == r->pOrdIndex)
566 return p_Setm_TotalDegree;
567 if (r->typ[0].ord_typ == ro_wp &&
568 r->typ[0].data.wp.start == 1 &&
569 r->typ[0].data.wp.end == r->N &&
570 r->typ[0].data.wp.place == r->pOrdIndex &&
571 r->typ[0].data.wp.weights == r->firstwv)
572 return p_Setm_WFirstTotalDegree;
573 }
574 return p_Setm_General;
575 }
576
577
578 /* -------------------------------------------------------------------*/
579 /* several possibilities for pFDeg: the degree of the head term */
580
581 /* comptible with ordering */
p_Deg(poly a,const ring r)582 long p_Deg(poly a, const ring r)
583 {
584 p_LmCheckPolyRing(a, r);
585 // assume(p_GetOrder(a, r) == p_WTotaldegree(a, r)); // WRONG assume!
586 return p_GetOrder(a, r);
587 }
588
589 // p_WTotalDegree for weighted orderings
590 // whose first block covers all variables
p_WFirstTotalDegree(poly p,const ring r)591 long p_WFirstTotalDegree(poly p, const ring r)
592 {
593 int i;
594 long sum = 0;
595
596 for (i=1; i<= r->firstBlockEnds; i++)
597 {
598 sum += p_GetExp(p, i, r)*r->firstwv[i-1];
599 }
600 return sum;
601 }
602
603 /*2
604 * compute the degree of the leading monomial of p
605 * with respect to weigths from the ordering
606 * the ordering is not compatible with degree so do not use p->Order
607 */
p_WTotaldegree(poly p,const ring r)608 long p_WTotaldegree(poly p, const ring r)
609 {
610 p_LmCheckPolyRing(p, r);
611 int i, k;
612 long j =0;
613
614 // iterate through each block:
615 for (i=0;r->order[i]!=0;i++)
616 {
617 int b0=r->block0[i];
618 int b1=r->block1[i];
619 switch(r->order[i])
620 {
621 case ringorder_M:
622 for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
623 { // in jedem block:
624 j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/]*r->OrdSgn;
625 }
626 break;
627 case ringorder_am:
628 b1=si_min(b1,r->N);
629 /* no break, continue as ringorder_a*/
630 case ringorder_a:
631 for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
632 { // only one line
633 j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/];
634 }
635 return j*r->OrdSgn;
636 case ringorder_wp:
637 case ringorder_ws:
638 case ringorder_Wp:
639 case ringorder_Ws:
640 for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
641 { // in jedem block:
642 j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/];
643 }
644 break;
645 case ringorder_lp:
646 case ringorder_ls:
647 case ringorder_rs:
648 case ringorder_dp:
649 case ringorder_ds:
650 case ringorder_Dp:
651 case ringorder_Ds:
652 case ringorder_rp:
653 for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
654 {
655 j+= p_GetExp(p,k,r);
656 }
657 break;
658 case ringorder_a64:
659 {
660 int64* w=(int64*)r->wvhdl[i];
661 for (k=0;k<=(b1 /*r->block1[i]*/ - b0 /*r->block0[i]*/);k++)
662 {
663 //there should be added a line which checks if w[k]>2^31
664 j+= p_GetExp(p,k+1, r)*(long)w[k];
665 }
666 //break;
667 return j;
668 }
669 case ringorder_c: /* nothing to do*/
670 case ringorder_C: /* nothing to do*/
671 case ringorder_S: /* nothing to do*/
672 case ringorder_s: /* nothing to do*/
673 case ringorder_IS: /* nothing to do */
674 case ringorder_unspec: /* to make clang happy, does not occur*/
675 case ringorder_no: /* to make clang happy, does not occur*/
676 case ringorder_L: /* to make clang happy, does not occur*/
677 case ringorder_aa: /* ignored by p_WTotaldegree*/
678 break;
679 /* no default: all orderings covered */
680 }
681 }
682 return j;
683 }
684
p_DegW(poly p,const int * w,const ring R)685 long p_DegW(poly p, const int *w, const ring R)
686 {
687 p_Test(p, R);
688 assume( w != NULL );
689 long r=-LONG_MAX;
690
691 while (p!=NULL)
692 {
693 long t=totaldegreeWecart_IV(p,R,w);
694 if (t>r) r=t;
695 pIter(p);
696 }
697 return r;
698 }
699
p_Weight(int i,const ring r)700 int p_Weight(int i, const ring r)
701 {
702 if ((r->firstwv==NULL) || (i>r->firstBlockEnds))
703 {
704 return 1;
705 }
706 return r->firstwv[i-1];
707 }
708
p_WDegree(poly p,const ring r)709 long p_WDegree(poly p, const ring r)
710 {
711 if (r->firstwv==NULL) return p_Totaldegree(p, r);
712 p_LmCheckPolyRing(p, r);
713 int i;
714 long j =0;
715
716 for(i=1;i<=r->firstBlockEnds;i++)
717 j+=p_GetExp(p, i, r)*r->firstwv[i-1];
718
719 for (;i<=rVar(r);i++)
720 j+=p_GetExp(p,i, r)*p_Weight(i, r);
721
722 return j;
723 }
724
725
726 /* ---------------------------------------------------------------------*/
727 /* several possibilities for pLDeg: the maximal degree of a monomial in p*/
728 /* compute in l also the pLength of p */
729
730 /*2
731 * compute the length of a polynomial (in l)
732 * and the degree of the monomial with maximal degree: the last one
733 */
pLDeg0(poly p,int * l,const ring r)734 long pLDeg0(poly p,int *l, const ring r)
735 {
736 p_CheckPolyRing(p, r);
737 long k= p_GetComp(p, r);
738 int ll=1;
739
740 if (k > 0)
741 {
742 while ((pNext(p)!=NULL) && (__p_GetComp(pNext(p), r)==k))
743 {
744 pIter(p);
745 ll++;
746 }
747 }
748 else
749 {
750 while (pNext(p)!=NULL)
751 {
752 pIter(p);
753 ll++;
754 }
755 }
756 *l=ll;
757 return r->pFDeg(p, r);
758 }
759
760 /*2
761 * compute the length of a polynomial (in l)
762 * and the degree of the monomial with maximal degree: the last one
763 * but search in all components before syzcomp
764 */
pLDeg0c(poly p,int * l,const ring r)765 long pLDeg0c(poly p,int *l, const ring r)
766 {
767 assume(p!=NULL);
768 p_Test(p,r);
769 p_CheckPolyRing(p, r);
770 long o;
771 int ll=1;
772
773 if (! rIsSyzIndexRing(r))
774 {
775 while (pNext(p) != NULL)
776 {
777 pIter(p);
778 ll++;
779 }
780 o = r->pFDeg(p, r);
781 }
782 else
783 {
784 int curr_limit = rGetCurrSyzLimit(r);
785 poly pp = p;
786 while ((p=pNext(p))!=NULL)
787 {
788 if (__p_GetComp(p, r)<=curr_limit/*syzComp*/)
789 ll++;
790 else break;
791 pp = p;
792 }
793 p_Test(pp,r);
794 o = r->pFDeg(pp, r);
795 }
796 *l=ll;
797 return o;
798 }
799
800 /*2
801 * compute the length of a polynomial (in l)
802 * and the degree of the monomial with maximal degree: the first one
803 * this works for the polynomial case with degree orderings
804 * (both c,dp and dp,c)
805 */
pLDegb(poly p,int * l,const ring r)806 long pLDegb(poly p,int *l, const ring r)
807 {
808 p_CheckPolyRing(p, r);
809 long k= p_GetComp(p, r);
810 long o = r->pFDeg(p, r);
811 int ll=1;
812
813 if (k != 0)
814 {
815 while (((p=pNext(p))!=NULL) && (__p_GetComp(p, r)==k))
816 {
817 ll++;
818 }
819 }
820 else
821 {
822 while ((p=pNext(p)) !=NULL)
823 {
824 ll++;
825 }
826 }
827 *l=ll;
828 return o;
829 }
830
831 /*2
832 * compute the length of a polynomial (in l)
833 * and the degree of the monomial with maximal degree:
834 * this is NOT the last one, we have to look for it
835 */
pLDeg1(poly p,int * l,const ring r)836 long pLDeg1(poly p,int *l, const ring r)
837 {
838 p_CheckPolyRing(p, r);
839 long k= p_GetComp(p, r);
840 int ll=1;
841 long t,max;
842
843 max=r->pFDeg(p, r);
844 if (k > 0)
845 {
846 while (((p=pNext(p))!=NULL) && (__p_GetComp(p, r)==k))
847 {
848 t=r->pFDeg(p, r);
849 if (t>max) max=t;
850 ll++;
851 }
852 }
853 else
854 {
855 while ((p=pNext(p))!=NULL)
856 {
857 t=r->pFDeg(p, r);
858 if (t>max) max=t;
859 ll++;
860 }
861 }
862 *l=ll;
863 return max;
864 }
865
866 /*2
867 * compute the length of a polynomial (in l)
868 * and the degree of the monomial with maximal degree:
869 * this is NOT the last one, we have to look for it
870 * in all components
871 */
pLDeg1c(poly p,int * l,const ring r)872 long pLDeg1c(poly p,int *l, const ring r)
873 {
874 p_CheckPolyRing(p, r);
875 int ll=1;
876 long t,max;
877
878 max=r->pFDeg(p, r);
879 if (rIsSyzIndexRing(r))
880 {
881 long limit = rGetCurrSyzLimit(r);
882 while ((p=pNext(p))!=NULL)
883 {
884 if (__p_GetComp(p, r)<=limit)
885 {
886 if ((t=r->pFDeg(p, r))>max) max=t;
887 ll++;
888 }
889 else break;
890 }
891 }
892 else
893 {
894 while ((p=pNext(p))!=NULL)
895 {
896 if ((t=r->pFDeg(p, r))>max) max=t;
897 ll++;
898 }
899 }
900 *l=ll;
901 return max;
902 }
903
904 // like pLDeg1, only pFDeg == pDeg
pLDeg1_Deg(poly p,int * l,const ring r)905 long pLDeg1_Deg(poly p,int *l, const ring r)
906 {
907 assume(r->pFDeg == p_Deg);
908 p_CheckPolyRing(p, r);
909 long k= p_GetComp(p, r);
910 int ll=1;
911 long t,max;
912
913 max=p_GetOrder(p, r);
914 if (k > 0)
915 {
916 while (((p=pNext(p))!=NULL) && (__p_GetComp(p, r)==k))
917 {
918 t=p_GetOrder(p, r);
919 if (t>max) max=t;
920 ll++;
921 }
922 }
923 else
924 {
925 while ((p=pNext(p))!=NULL)
926 {
927 t=p_GetOrder(p, r);
928 if (t>max) max=t;
929 ll++;
930 }
931 }
932 *l=ll;
933 return max;
934 }
935
pLDeg1c_Deg(poly p,int * l,const ring r)936 long pLDeg1c_Deg(poly p,int *l, const ring r)
937 {
938 assume(r->pFDeg == p_Deg);
939 p_CheckPolyRing(p, r);
940 int ll=1;
941 long t,max;
942
943 max=p_GetOrder(p, r);
944 if (rIsSyzIndexRing(r))
945 {
946 long limit = rGetCurrSyzLimit(r);
947 while ((p=pNext(p))!=NULL)
948 {
949 if (__p_GetComp(p, r)<=limit)
950 {
951 if ((t=p_GetOrder(p, r))>max) max=t;
952 ll++;
953 }
954 else break;
955 }
956 }
957 else
958 {
959 while ((p=pNext(p))!=NULL)
960 {
961 if ((t=p_GetOrder(p, r))>max) max=t;
962 ll++;
963 }
964 }
965 *l=ll;
966 return max;
967 }
968
969 // like pLDeg1, only pFDeg == pTotoalDegree
pLDeg1_Totaldegree(poly p,int * l,const ring r)970 long pLDeg1_Totaldegree(poly p,int *l, const ring r)
971 {
972 p_CheckPolyRing(p, r);
973 long k= p_GetComp(p, r);
974 int ll=1;
975 long t,max;
976
977 max=p_Totaldegree(p, r);
978 if (k > 0)
979 {
980 while (((p=pNext(p))!=NULL) && (__p_GetComp(p, r)==k))
981 {
982 t=p_Totaldegree(p, r);
983 if (t>max) max=t;
984 ll++;
985 }
986 }
987 else
988 {
989 while ((p=pNext(p))!=NULL)
990 {
991 t=p_Totaldegree(p, r);
992 if (t>max) max=t;
993 ll++;
994 }
995 }
996 *l=ll;
997 return max;
998 }
999
pLDeg1c_Totaldegree(poly p,int * l,const ring r)1000 long pLDeg1c_Totaldegree(poly p,int *l, const ring r)
1001 {
1002 p_CheckPolyRing(p, r);
1003 int ll=1;
1004 long t,max;
1005
1006 max=p_Totaldegree(p, r);
1007 if (rIsSyzIndexRing(r))
1008 {
1009 long limit = rGetCurrSyzLimit(r);
1010 while ((p=pNext(p))!=NULL)
1011 {
1012 if (__p_GetComp(p, r)<=limit)
1013 {
1014 if ((t=p_Totaldegree(p, r))>max) max=t;
1015 ll++;
1016 }
1017 else break;
1018 }
1019 }
1020 else
1021 {
1022 while ((p=pNext(p))!=NULL)
1023 {
1024 if ((t=p_Totaldegree(p, r))>max) max=t;
1025 ll++;
1026 }
1027 }
1028 *l=ll;
1029 return max;
1030 }
1031
1032 // like pLDeg1, only pFDeg == p_WFirstTotalDegree
pLDeg1_WFirstTotalDegree(poly p,int * l,const ring r)1033 long pLDeg1_WFirstTotalDegree(poly p,int *l, const ring r)
1034 {
1035 p_CheckPolyRing(p, r);
1036 long k= p_GetComp(p, r);
1037 int ll=1;
1038 long t,max;
1039
1040 max=p_WFirstTotalDegree(p, r);
1041 if (k > 0)
1042 {
1043 while (((p=pNext(p))!=NULL) && (__p_GetComp(p, r)==k))
1044 {
1045 t=p_WFirstTotalDegree(p, r);
1046 if (t>max) max=t;
1047 ll++;
1048 }
1049 }
1050 else
1051 {
1052 while ((p=pNext(p))!=NULL)
1053 {
1054 t=p_WFirstTotalDegree(p, r);
1055 if (t>max) max=t;
1056 ll++;
1057 }
1058 }
1059 *l=ll;
1060 return max;
1061 }
1062
pLDeg1c_WFirstTotalDegree(poly p,int * l,const ring r)1063 long pLDeg1c_WFirstTotalDegree(poly p,int *l, const ring r)
1064 {
1065 p_CheckPolyRing(p, r);
1066 int ll=1;
1067 long t,max;
1068
1069 max=p_WFirstTotalDegree(p, r);
1070 if (rIsSyzIndexRing(r))
1071 {
1072 long limit = rGetCurrSyzLimit(r);
1073 while ((p=pNext(p))!=NULL)
1074 {
1075 if (__p_GetComp(p, r)<=limit)
1076 {
1077 if ((t=p_Totaldegree(p, r))>max) max=t;
1078 ll++;
1079 }
1080 else break;
1081 }
1082 }
1083 else
1084 {
1085 while ((p=pNext(p))!=NULL)
1086 {
1087 if ((t=p_Totaldegree(p, r))>max) max=t;
1088 ll++;
1089 }
1090 }
1091 *l=ll;
1092 return max;
1093 }
1094
1095 /***************************************************************
1096 *
1097 * Maximal Exponent business
1098 *
1099 ***************************************************************/
1100
1101 static inline unsigned long
p_GetMaxExpL2(unsigned long l1,unsigned long l2,const ring r,unsigned long number_of_exp)1102 p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r,
1103 unsigned long number_of_exp)
1104 {
1105 const unsigned long bitmask = r->bitmask;
1106 unsigned long ml1 = l1 & bitmask;
1107 unsigned long ml2 = l2 & bitmask;
1108 unsigned long max = (ml1 > ml2 ? ml1 : ml2);
1109 unsigned long j = number_of_exp - 1;
1110
1111 if (j > 0)
1112 {
1113 unsigned long mask = bitmask << r->BitsPerExp;
1114 while (1)
1115 {
1116 ml1 = l1 & mask;
1117 ml2 = l2 & mask;
1118 max |= ((ml1 > ml2 ? ml1 : ml2) & mask);
1119 j--;
1120 if (j == 0) break;
1121 mask = mask << r->BitsPerExp;
1122 }
1123 }
1124 return max;
1125 }
1126
1127 static inline unsigned long
p_GetMaxExpL2(unsigned long l1,unsigned long l2,const ring r)1128 p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r)
1129 {
1130 return p_GetMaxExpL2(l1, l2, r, r->ExpPerLong);
1131 }
1132
p_GetMaxExpP(poly p,const ring r)1133 poly p_GetMaxExpP(poly p, const ring r)
1134 {
1135 p_CheckPolyRing(p, r);
1136 if (p == NULL) return p_Init(r);
1137 poly max = p_LmInit(p, r);
1138 pIter(p);
1139 if (p == NULL) return max;
1140 int i, offset;
1141 unsigned long l_p, l_max;
1142 unsigned long divmask = r->divmask;
1143
1144 do
1145 {
1146 offset = r->VarL_Offset[0];
1147 l_p = p->exp[offset];
1148 l_max = max->exp[offset];
1149 // do the divisibility trick to find out whether l has an exponent
1150 if (l_p > l_max ||
1151 (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1152 max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
1153
1154 for (i=1; i<r->VarL_Size; i++)
1155 {
1156 offset = r->VarL_Offset[i];
1157 l_p = p->exp[offset];
1158 l_max = max->exp[offset];
1159 // do the divisibility trick to find out whether l has an exponent
1160 if (l_p > l_max ||
1161 (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1162 max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
1163 }
1164 pIter(p);
1165 }
1166 while (p != NULL);
1167 return max;
1168 }
1169
p_GetMaxExpL(poly p,const ring r,unsigned long l_max)1170 unsigned long p_GetMaxExpL(poly p, const ring r, unsigned long l_max)
1171 {
1172 unsigned long l_p, divmask = r->divmask;
1173 int i;
1174
1175 while (p != NULL)
1176 {
1177 l_p = p->exp[r->VarL_Offset[0]];
1178 if (l_p > l_max ||
1179 (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1180 l_max = p_GetMaxExpL2(l_max, l_p, r);
1181 for (i=1; i<r->VarL_Size; i++)
1182 {
1183 l_p = p->exp[r->VarL_Offset[i]];
1184 // do the divisibility trick to find out whether l has an exponent
1185 if (l_p > l_max ||
1186 (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1187 l_max = p_GetMaxExpL2(l_max, l_p, r);
1188 }
1189 pIter(p);
1190 }
1191 return l_max;
1192 }
1193
1194
1195
1196
1197 /***************************************************************
1198 *
1199 * Misc things
1200 *
1201 ***************************************************************/
1202 // returns TRUE, if all monoms have the same component
p_OneComp(poly p,const ring r)1203 BOOLEAN p_OneComp(poly p, const ring r)
1204 {
1205 if(p!=NULL)
1206 {
1207 long i = p_GetComp(p, r);
1208 while (pNext(p)!=NULL)
1209 {
1210 pIter(p);
1211 if(i != p_GetComp(p, r)) return FALSE;
1212 }
1213 }
1214 return TRUE;
1215 }
1216
1217 /*2
1218 *test if a monomial /head term is a pure power,
1219 * i.e. depends on only one variable
1220 */
p_IsPurePower(const poly p,const ring r)1221 int p_IsPurePower(const poly p, const ring r)
1222 {
1223 int i,k=0;
1224
1225 for (i=r->N;i;i--)
1226 {
1227 if (p_GetExp(p,i, r)!=0)
1228 {
1229 if(k!=0) return 0;
1230 k=i;
1231 }
1232 }
1233 return k;
1234 }
1235
1236 /*2
1237 *test if a polynomial is univariate
1238 * return -1 for constant,
1239 * 0 for not univariate,s
1240 * i if dep. on var(i)
1241 */
p_IsUnivariate(poly p,const ring r)1242 int p_IsUnivariate(poly p, const ring r)
1243 {
1244 int i,k=-1;
1245
1246 while (p!=NULL)
1247 {
1248 for (i=r->N;i;i--)
1249 {
1250 if (p_GetExp(p,i, r)!=0)
1251 {
1252 if((k!=-1)&&(k!=i)) return 0;
1253 k=i;
1254 }
1255 }
1256 pIter(p);
1257 }
1258 return k;
1259 }
1260
1261 // set entry e[i] to 1 if var(i) occurs in p, ignore var(j) if e[j]>0
p_GetVariables(poly p,int * e,const ring r)1262 int p_GetVariables(poly p, int * e, const ring r)
1263 {
1264 int i;
1265 int n=0;
1266 while(p!=NULL)
1267 {
1268 n=0;
1269 for(i=r->N; i>0; i--)
1270 {
1271 if(e[i]==0)
1272 {
1273 if (p_GetExp(p,i,r)>0)
1274 {
1275 e[i]=1;
1276 n++;
1277 }
1278 }
1279 else
1280 n++;
1281 }
1282 if (n==r->N) break;
1283 pIter(p);
1284 }
1285 return n;
1286 }
1287
1288
1289 /*2
1290 * returns a polynomial representing the integer i
1291 */
p_ISet(long i,const ring r)1292 poly p_ISet(long i, const ring r)
1293 {
1294 poly rc = NULL;
1295 if (i!=0)
1296 {
1297 rc = p_Init(r);
1298 pSetCoeff0(rc,n_Init(i,r->cf));
1299 if (n_IsZero(pGetCoeff(rc),r->cf))
1300 p_LmDelete(&rc,r);
1301 }
1302 return rc;
1303 }
1304
1305 /*2
1306 * an optimized version of p_ISet for the special case 1
1307 */
p_One(const ring r)1308 poly p_One(const ring r)
1309 {
1310 poly rc = p_Init(r);
1311 pSetCoeff0(rc,n_Init(1,r->cf));
1312 return rc;
1313 }
1314
p_Split(poly p,poly * h)1315 void p_Split(poly p, poly *h)
1316 {
1317 *h=pNext(p);
1318 pNext(p)=NULL;
1319 }
1320
1321 /*2
1322 * pair has no common factor ? or is no polynomial
1323 */
p_HasNotCF(poly p1,poly p2,const ring r)1324 BOOLEAN p_HasNotCF(poly p1, poly p2, const ring r)
1325 {
1326
1327 if (p_GetComp(p1,r) > 0 || p_GetComp(p2,r) > 0)
1328 return FALSE;
1329 int i = rVar(r);
1330 loop
1331 {
1332 if ((p_GetExp(p1, i, r) > 0) && (p_GetExp(p2, i, r) > 0))
1333 return FALSE;
1334 i--;
1335 if (i == 0)
1336 return TRUE;
1337 }
1338 }
1339
p_HasNotCFRing(poly p1,poly p2,const ring r)1340 BOOLEAN p_HasNotCFRing(poly p1, poly p2, const ring r)
1341 {
1342
1343 if (p_GetComp(p1,r) > 0 || p_GetComp(p2,r) > 0)
1344 return FALSE;
1345 int i = rVar(r);
1346 loop
1347 {
1348 if ((p_GetExp(p1, i, r) > 0) && (p_GetExp(p2, i, r) > 0))
1349 return FALSE;
1350 i--;
1351 if (i == 0) {
1352 if (n_DivBy(pGetCoeff(p1), pGetCoeff(p2), r->cf) ||
1353 n_DivBy(pGetCoeff(p2), pGetCoeff(p1), r->cf)) {
1354 return FALSE;
1355 } else {
1356 return TRUE;
1357 }
1358 }
1359 }
1360 }
1361
1362 /*2
1363 * convert monomial given as string to poly, e.g. 1x3y5z
1364 */
p_Read(const char * st,poly & rc,const ring r)1365 const char * p_Read(const char *st, poly &rc, const ring r)
1366 {
1367 if (r==NULL) { rc=NULL;return st;}
1368 int i,j;
1369 rc = p_Init(r);
1370 const char *s = n_Read(st,&(p_GetCoeff(rc, r)),r->cf);
1371 if (s==st)
1372 /* i.e. it does not start with a coeff: test if it is a ringvar*/
1373 {
1374 j = r_IsRingVar(s,r->names,r->N);
1375 if (j >= 0)
1376 {
1377 p_IncrExp(rc,1+j,r);
1378 while (*s!='\0') s++;
1379 goto done;
1380 }
1381 }
1382 while (*s!='\0')
1383 {
1384 char ss[2];
1385 ss[0] = *s++;
1386 ss[1] = '\0';
1387 j = r_IsRingVar(ss,r->names,r->N);
1388 if (j >= 0)
1389 {
1390 const char *s_save=s;
1391 s = eati(s,&i);
1392 if (((unsigned long)i) > r->bitmask/2)
1393 {
1394 // exponent to large: it is not a monomial
1395 p_LmDelete(&rc,r);
1396 return s_save;
1397 }
1398 p_AddExp(rc,1+j, (long)i, r);
1399 }
1400 else
1401 {
1402 // 1st char of is not a varname
1403 // We return the parsed polynomial nevertheless. This is needed when
1404 // we are parsing coefficients in a rational function field.
1405 s--;
1406 break;
1407 }
1408 }
1409 done:
1410 if (n_IsZero(pGetCoeff(rc),r->cf)) p_LmDelete(&rc,r);
1411 else
1412 {
1413 #ifdef HAVE_PLURAL
1414 // in super-commutative ring
1415 // squares of anti-commutative variables are zeroes!
1416 if(rIsSCA(r))
1417 {
1418 const unsigned int iFirstAltVar = scaFirstAltVar(r);
1419 const unsigned int iLastAltVar = scaLastAltVar(r);
1420
1421 assume(rc != NULL);
1422
1423 for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
1424 if( p_GetExp(rc, k, r) > 1 )
1425 {
1426 p_LmDelete(&rc, r);
1427 goto finish;
1428 }
1429 }
1430 #endif
1431
1432 p_Setm(rc,r);
1433 }
1434 finish:
1435 return s;
1436 }
p_mInit(const char * st,BOOLEAN & ok,const ring r)1437 poly p_mInit(const char *st, BOOLEAN &ok, const ring r)
1438 {
1439 poly p;
1440 const char *s=p_Read(st,p,r);
1441 if (*s!='\0')
1442 {
1443 if ((s!=st)&&isdigit(st[0]))
1444 {
1445 errorreported=TRUE;
1446 }
1447 ok=FALSE;
1448 p_Delete(&p,r);
1449 return NULL;
1450 }
1451 p_Test(p,r);
1452 ok=!errorreported;
1453 return p;
1454 }
1455
1456 /*2
1457 * returns a polynomial representing the number n
1458 * destroys n
1459 */
p_NSet(number n,const ring r)1460 poly p_NSet(number n, const ring r)
1461 {
1462 if (n_IsZero(n,r->cf))
1463 {
1464 n_Delete(&n, r->cf);
1465 return NULL;
1466 }
1467 else
1468 {
1469 poly rc = p_Init(r);
1470 pSetCoeff0(rc,n);
1471 return rc;
1472 }
1473 }
1474 /*2
1475 * assumes that LM(a) = LM(b)*m, for some monomial m,
1476 * returns the multiplicant m,
1477 * leaves a and b unmodified
1478 */
p_MDivide(poly a,poly b,const ring r)1479 poly p_MDivide(poly a, poly b, const ring r)
1480 {
1481 assume((p_GetComp(a,r)==p_GetComp(b,r)) || (p_GetComp(b,r)==0));
1482 int i;
1483 poly result = p_Init(r);
1484
1485 for(i=(int)r->N; i; i--)
1486 p_SetExp(result,i, p_GetExp(a,i,r)- p_GetExp(b,i,r),r);
1487 p_SetComp(result, p_GetComp(a,r) - p_GetComp(b,r),r);
1488 p_Setm(result,r);
1489 return result;
1490 }
1491
p_Div_nn(poly p,const number n,const ring r)1492 poly p_Div_nn(poly p, const number n, const ring r)
1493 {
1494 pAssume(!n_IsZero(n,r->cf));
1495 p_Test(p, r);
1496 poly result = p;
1497 poly prev = NULL;
1498 while (p!=NULL)
1499 {
1500 number nc = n_Div(pGetCoeff(p),n,r->cf);
1501 if (!n_IsZero(nc,r->cf))
1502 {
1503 p_SetCoeff(p,nc,r);
1504 prev=p;
1505 pIter(p);
1506 }
1507 else
1508 {
1509 if (prev==NULL)
1510 {
1511 p_LmDelete(&result,r);
1512 p=result;
1513 }
1514 else
1515 {
1516 p_LmDelete(&pNext(prev),r);
1517 p=pNext(prev);
1518 }
1519 }
1520 }
1521 p_Test(result,r);
1522 return(result);
1523 }
1524
p_Div_mm(poly p,const poly m,const ring r)1525 poly p_Div_mm(poly p, const poly m, const ring r)
1526 {
1527 p_Test(p, r);
1528 p_Test(m, r);
1529 poly result = p;
1530 poly prev = NULL;
1531 number n=pGetCoeff(m);
1532 while (p!=NULL)
1533 {
1534 number nc = n_Div(pGetCoeff(p),n,r->cf);
1535 n_Normalize(nc,r->cf);
1536 if (!n_IsZero(nc,r->cf))
1537 {
1538 p_SetCoeff(p,nc,r);
1539 prev=p;
1540 p_ExpVectorSub(p,m,r);
1541 pIter(p);
1542 }
1543 else
1544 {
1545 if (prev==NULL)
1546 {
1547 p_LmDelete(&result,r);
1548 p=result;
1549 }
1550 else
1551 {
1552 p_LmDelete(&pNext(prev),r);
1553 p=pNext(prev);
1554 }
1555 }
1556 }
1557 p_Test(result,r);
1558 return(result);
1559 }
1560
1561 /*2
1562 * divides a by the monomial b, ignores monomials which are not divisible
1563 * assumes that b is not NULL, destroyes b
1564 */
p_DivideM(poly a,poly b,const ring r)1565 poly p_DivideM(poly a, poly b, const ring r)
1566 {
1567 if (a==NULL) { p_Delete(&b,r); return NULL; }
1568 poly result=a;
1569
1570 if(!p_IsConstant(b,r))
1571 {
1572 if (rIsNCRing(r))
1573 {
1574 WerrorS("p_DivideM not implemented for non-commuative rings");
1575 return NULL;
1576 }
1577 poly prev=NULL;
1578 while (a!=NULL)
1579 {
1580 if (p_DivisibleBy(b,a,r))
1581 {
1582 p_ExpVectorSub(a,b,r);
1583 prev=a;
1584 pIter(a);
1585 }
1586 else
1587 {
1588 if (prev==NULL)
1589 {
1590 p_LmDelete(&result,r);
1591 a=result;
1592 }
1593 else
1594 {
1595 p_LmDelete(&pNext(prev),r);
1596 a=pNext(prev);
1597 }
1598 }
1599 }
1600 }
1601 if (result!=NULL)
1602 {
1603 number inv=pGetCoeff(b);
1604 //if ((!rField_is_Ring(r)) || n_IsUnit(inv,r->cf))
1605 if (rField_is_Zp(r))
1606 {
1607 inv = n_Invers(inv,r->cf);
1608 __p_Mult_nn(result,inv,r);
1609 n_Delete(&inv, r->cf);
1610 }
1611 else
1612 {
1613 result = p_Div_nn(result,inv,r);
1614 }
1615 }
1616 p_Delete(&b, r);
1617 return result;
1618 }
1619
pp_DivideM(poly a,poly b,const ring r)1620 poly pp_DivideM(poly a, poly b, const ring r)
1621 {
1622 if (a==NULL) { return NULL; }
1623 // TODO: better implementation without copying a,b
1624 return p_DivideM(p_Copy(a,r),p_Head(b,r),r);
1625 }
1626
1627 #ifdef HAVE_RINGS
1628 /* TRUE iff LT(f) | LT(g) */
p_DivisibleByRingCase(poly f,poly g,const ring r)1629 BOOLEAN p_DivisibleByRingCase(poly f, poly g, const ring r)
1630 {
1631 int exponent;
1632 for(int i = (int)rVar(r); i>0; i--)
1633 {
1634 exponent = p_GetExp(g, i, r) - p_GetExp(f, i, r);
1635 if (exponent < 0) return FALSE;
1636 }
1637 return n_DivBy(pGetCoeff(g), pGetCoeff(f), r->cf);
1638 }
1639 #endif
1640
1641 // returns the LCM of the head terms of a and b in *m
p_Lcm(const poly a,const poly b,poly m,const ring r)1642 void p_Lcm(const poly a, const poly b, poly m, const ring r)
1643 {
1644 for (int i=r->N; i; --i)
1645 p_SetExp(m,i, si_max( p_GetExp(a,i,r), p_GetExp(b,i,r)),r);
1646
1647 p_SetComp(m, si_max(p_GetComp(a,r), p_GetComp(b,r)),r);
1648 /* Don't do a pSetm here, otherwise hres/lres chockes */
1649 }
1650
p_Lcm(const poly a,const poly b,const ring r)1651 poly p_Lcm(const poly a, const poly b, const ring r)
1652 {
1653 poly m=p_Init(r);
1654 p_Lcm(a, b, m, r);
1655 p_Setm(m,r);
1656 return(m);
1657 }
1658
1659 #ifdef HAVE_RATGRING
1660 /*2
1661 * returns the rational LCM of the head terms of a and b
1662 * without coefficient!!!
1663 */
p_LcmRat(const poly a,const poly b,const long lCompM,const ring r)1664 poly p_LcmRat(const poly a, const poly b, const long lCompM, const ring r)
1665 {
1666 poly m = // p_One( r);
1667 p_Init(r);
1668
1669 // const int (currRing->N) = r->N;
1670
1671 // for (int i = (currRing->N); i>=r->real_var_start; i--)
1672 for (int i = r->real_var_end; i>=r->real_var_start; i--)
1673 {
1674 const int lExpA = p_GetExp (a, i, r);
1675 const int lExpB = p_GetExp (b, i, r);
1676
1677 p_SetExp (m, i, si_max(lExpA, lExpB), r);
1678 }
1679
1680 p_SetComp (m, lCompM, r);
1681 p_Setm(m,r);
1682 n_New(&(p_GetCoeff(m, r)), r);
1683
1684 return(m);
1685 };
1686
p_LmDeleteAndNextRat(poly * p,int ishift,ring r)1687 void p_LmDeleteAndNextRat(poly *p, int ishift, ring r)
1688 {
1689 /* modifies p*/
1690 // Print("start: "); Print(" "); p_wrp(*p,r);
1691 p_LmCheckPolyRing2(*p, r);
1692 poly q = p_Head(*p,r);
1693 const long cmp = p_GetComp(*p, r);
1694 while ( ( (*p)!=NULL ) && ( p_Comp_k_n(*p, q, ishift+1, r) ) && (p_GetComp(*p, r) == cmp) )
1695 {
1696 p_LmDelete(p,r);
1697 // Print("while: ");p_wrp(*p,r);Print(" ");
1698 }
1699 // p_wrp(*p,r);Print(" ");
1700 // PrintS("end\n");
1701 p_LmDelete(&q,r);
1702 }
1703
1704
1705 /* returns x-coeff of p, i.e. a poly in x, s.t. corresponding xd-monomials
1706 have the same D-part and the component 0
1707 does not destroy p
1708 */
p_GetCoeffRat(poly p,int ishift,ring r)1709 poly p_GetCoeffRat(poly p, int ishift, ring r)
1710 {
1711 poly q = pNext(p);
1712 poly res; // = p_Head(p,r);
1713 res = p_GetExp_k_n(p, ishift+1, r->N, r); // does pSetm internally
1714 p_SetCoeff(res,n_Copy(p_GetCoeff(p,r),r),r);
1715 poly s;
1716 long cmp = p_GetComp(p, r);
1717 while ( (q!= NULL) && (p_Comp_k_n(p, q, ishift+1, r)) && (p_GetComp(q, r) == cmp) )
1718 {
1719 s = p_GetExp_k_n(q, ishift+1, r->N, r);
1720 p_SetCoeff(s,n_Copy(p_GetCoeff(q,r),r),r);
1721 res = p_Add_q(res,s,r);
1722 q = pNext(q);
1723 }
1724 cmp = 0;
1725 p_SetCompP(res,cmp,r);
1726 return res;
1727 }
1728
1729
1730
p_ContentRat(poly & ph,const ring r)1731 void p_ContentRat(poly &ph, const ring r)
1732 // changes ph
1733 // for rat coefficients in K(x1,..xN)
1734 {
1735 // init array of RatLeadCoeffs
1736 // poly p_GetCoeffRat(poly p, int ishift, ring r);
1737
1738 int len=pLength(ph);
1739 poly *C = (poly *)omAlloc0((len+1)*sizeof(poly)); //rat coeffs
1740 poly *LM = (poly *)omAlloc0((len+1)*sizeof(poly)); // rat lead terms
1741 int *D = (int *)omAlloc0((len+1)*sizeof(int)); //degrees of coeffs
1742 int *L = (int *)omAlloc0((len+1)*sizeof(int)); //lengths of coeffs
1743 int k = 0;
1744 poly p = p_Copy(ph, r); // ph will be needed below
1745 int mintdeg = p_Totaldegree(p, r);
1746 int minlen = len;
1747 int dd = 0; int i;
1748 int HasConstantCoef = 0;
1749 int is = r->real_var_start - 1;
1750 while (p!=NULL)
1751 {
1752 LM[k] = p_GetExp_k_n(p,1,is, r); // need LmRat istead of p_HeadRat(p, is, currRing); !
1753 C[k] = p_GetCoeffRat(p, is, r);
1754 D[k] = p_Totaldegree(C[k], r);
1755 mintdeg = si_min(mintdeg,D[k]);
1756 L[k] = pLength(C[k]);
1757 minlen = si_min(minlen,L[k]);
1758 if (p_IsConstant(C[k], r))
1759 {
1760 // C[k] = const, so the content will be numerical
1761 HasConstantCoef = 1;
1762 // smth like goto cleanup and return(pContent(p));
1763 }
1764 p_LmDeleteAndNextRat(&p, is, r);
1765 k++;
1766 }
1767
1768 // look for 1 element of minimal degree and of minimal length
1769 k--;
1770 poly d;
1771 int mindeglen = len;
1772 if (k<=0) // this poly is not a ratgring poly -> pContent
1773 {
1774 p_Delete(&C[0], r);
1775 p_Delete(&LM[0], r);
1776 p_ContentForGB(ph, r);
1777 goto cleanup;
1778 }
1779
1780 int pmindeglen;
1781 for(i=0; i<=k; i++)
1782 {
1783 if (D[i] == mintdeg)
1784 {
1785 if (L[i] < mindeglen)
1786 {
1787 mindeglen=L[i];
1788 pmindeglen = i;
1789 }
1790 }
1791 }
1792 d = p_Copy(C[pmindeglen], r);
1793 // there are dd>=1 mindeg elements
1794 // and pmideglen is the coordinate of one of the smallest among them
1795
1796 // poly g = singclap_gcd(p_Copy(p,r),p_Copy(q,r));
1797 // return naGcd(d,d2,currRing);
1798
1799 // adjoin pContentRat here?
1800 for(i=0; i<=k; i++)
1801 {
1802 d=singclap_gcd(d,p_Copy(C[i], r), r);
1803 if (p_Totaldegree(d, r)==0)
1804 {
1805 // cleanup, pContent, return
1806 p_Delete(&d, r);
1807 for(;k>=0;k--)
1808 {
1809 p_Delete(&C[k], r);
1810 p_Delete(&LM[k], r);
1811 }
1812 p_ContentForGB(ph, r);
1813 goto cleanup;
1814 }
1815 }
1816 for(i=0; i<=k; i++)
1817 {
1818 poly h=singclap_pdivide(C[i],d, r);
1819 p_Delete(&C[i], r);
1820 C[i]=h;
1821 }
1822
1823 // zusammensetzen,
1824 p=NULL; // just to be sure
1825 for(i=0; i<=k; i++)
1826 {
1827 p = p_Add_q(p, p_Mult_q(C[i],LM[i], r), r);
1828 C[i]=NULL; LM[i]=NULL;
1829 }
1830 p_Delete(&ph, r); // do not need it anymore
1831 ph = p;
1832 // aufraeumen, return
1833 cleanup:
1834 omFree(C);
1835 omFree(LM);
1836 omFree(D);
1837 omFree(L);
1838 }
1839
1840
1841 #endif
1842
1843
1844 /* assumes that p and divisor are univariate polynomials in r,
1845 mentioning the same variable;
1846 assumes divisor != NULL;
1847 p may be NULL;
1848 assumes a global monomial ordering in r;
1849 performs polynomial division of p by divisor:
1850 - afterwards p contains the remainder of the division, i.e.,
1851 p_before = result * divisor + p_afterwards;
1852 - if needResult == TRUE, then the method computes and returns 'result',
1853 otherwise NULL is returned (This parametrization can be used when
1854 one is only interested in the remainder of the division. In this
1855 case, the method will be slightly faster.)
1856 leaves divisor unmodified */
p_PolyDiv(poly & p,const poly divisor,const BOOLEAN needResult,const ring r)1857 poly p_PolyDiv(poly &p, const poly divisor, const BOOLEAN needResult, const ring r)
1858 {
1859 assume(divisor != NULL);
1860 if (p == NULL) return NULL;
1861
1862 poly result = NULL;
1863 number divisorLC = p_GetCoeff(divisor, r);
1864 int divisorLE = p_GetExp(divisor, 1, r);
1865 while ((p != NULL) && (p_Deg(p, r) >= p_Deg(divisor, r)))
1866 {
1867 /* determine t = LT(p) / LT(divisor) */
1868 poly t = p_ISet(1, r);
1869 number c = n_Div(p_GetCoeff(p, r), divisorLC, r->cf);
1870 n_Normalize(c,r->cf);
1871 p_SetCoeff(t, c, r);
1872 int e = p_GetExp(p, 1, r) - divisorLE;
1873 p_SetExp(t, 1, e, r);
1874 p_Setm(t, r);
1875 if (needResult) result = p_Add_q(result, p_Copy(t, r), r);
1876 p = p_Add_q(p, p_Neg(p_Mult_q(t, p_Copy(divisor, r), r), r), r);
1877 }
1878 return result;
1879 }
1880
1881 /*2
1882 * returns the partial differentiate of a by the k-th variable
1883 * does not destroy the input
1884 */
p_Diff(poly a,int k,const ring r)1885 poly p_Diff(poly a, int k, const ring r)
1886 {
1887 poly res, f, last;
1888 number t;
1889
1890 last = res = NULL;
1891 while (a!=NULL)
1892 {
1893 if (p_GetExp(a,k,r)!=0)
1894 {
1895 f = p_LmInit(a,r);
1896 t = n_Init(p_GetExp(a,k,r),r->cf);
1897 pSetCoeff0(f,n_Mult(t,pGetCoeff(a),r->cf));
1898 n_Delete(&t,r->cf);
1899 if (n_IsZero(pGetCoeff(f),r->cf))
1900 p_LmDelete(&f,r);
1901 else
1902 {
1903 p_DecrExp(f,k,r);
1904 p_Setm(f,r);
1905 if (res==NULL)
1906 {
1907 res=last=f;
1908 }
1909 else
1910 {
1911 pNext(last)=f;
1912 last=f;
1913 }
1914 }
1915 }
1916 pIter(a);
1917 }
1918 return res;
1919 }
1920
p_DiffOpM(poly a,poly b,BOOLEAN multiply,const ring r)1921 static poly p_DiffOpM(poly a, poly b,BOOLEAN multiply, const ring r)
1922 {
1923 int i,j,s;
1924 number n,h,hh;
1925 poly p=p_One(r);
1926 n=n_Mult(pGetCoeff(a),pGetCoeff(b),r->cf);
1927 for(i=rVar(r);i>0;i--)
1928 {
1929 s=p_GetExp(b,i,r);
1930 if (s<p_GetExp(a,i,r))
1931 {
1932 n_Delete(&n,r->cf);
1933 p_LmDelete(&p,r);
1934 return NULL;
1935 }
1936 if (multiply)
1937 {
1938 for(j=p_GetExp(a,i,r); j>0;j--)
1939 {
1940 h = n_Init(s,r->cf);
1941 hh=n_Mult(n,h,r->cf);
1942 n_Delete(&h,r->cf);
1943 n_Delete(&n,r->cf);
1944 n=hh;
1945 s--;
1946 }
1947 p_SetExp(p,i,s,r);
1948 }
1949 else
1950 {
1951 p_SetExp(p,i,s-p_GetExp(a,i,r),r);
1952 }
1953 }
1954 p_Setm(p,r);
1955 /*if (multiply)*/ p_SetCoeff(p,n,r);
1956 if (n_IsZero(n,r->cf)) p=p_LmDeleteAndNext(p,r); // return NULL as p is a monomial
1957 return p;
1958 }
1959
p_DiffOp(poly a,poly b,BOOLEAN multiply,const ring r)1960 poly p_DiffOp(poly a, poly b,BOOLEAN multiply, const ring r)
1961 {
1962 poly result=NULL;
1963 poly h;
1964 for(;a!=NULL;pIter(a))
1965 {
1966 for(h=b;h!=NULL;pIter(h))
1967 {
1968 result=p_Add_q(result,p_DiffOpM(a,h,multiply,r),r);
1969 }
1970 }
1971 return result;
1972 }
1973 /*2
1974 * subtract p2 from p1, p1 and p2 are destroyed
1975 * do not put attention on speed: the procedure is only used in the interpreter
1976 */
p_Sub(poly p1,poly p2,const ring r)1977 poly p_Sub(poly p1, poly p2, const ring r)
1978 {
1979 return p_Add_q(p1, p_Neg(p2,r),r);
1980 }
1981
1982 /*3
1983 * compute for a monomial m
1984 * the power m^exp, exp > 1
1985 * destroys p
1986 */
p_MonPower(poly p,int exp,const ring r)1987 static poly p_MonPower(poly p, int exp, const ring r)
1988 {
1989 int i;
1990
1991 if(!n_IsOne(pGetCoeff(p),r->cf))
1992 {
1993 number x, y;
1994 y = pGetCoeff(p);
1995 n_Power(y,exp,&x,r->cf);
1996 n_Delete(&y,r->cf);
1997 pSetCoeff0(p,x);
1998 }
1999 for (i=rVar(r); i!=0; i--)
2000 {
2001 p_MultExp(p,i, exp,r);
2002 }
2003 p_Setm(p,r);
2004 return p;
2005 }
2006
2007 /*3
2008 * compute for monomials p*q
2009 * destroys p, keeps q
2010 */
p_MonMult(poly p,poly q,const ring r)2011 static void p_MonMult(poly p, poly q, const ring r)
2012 {
2013 number x, y;
2014
2015 y = pGetCoeff(p);
2016 x = n_Mult(y,pGetCoeff(q),r->cf);
2017 n_Delete(&y,r->cf);
2018 pSetCoeff0(p,x);
2019 //for (int i=pVariables; i!=0; i--)
2020 //{
2021 // pAddExp(p,i, pGetExp(q,i));
2022 //}
2023 //p->Order += q->Order;
2024 p_ExpVectorAdd(p,q,r);
2025 }
2026
2027 /*3
2028 * compute for monomials p*q
2029 * keeps p, q
2030 */
p_MonMultC(poly p,poly q,const ring rr)2031 static poly p_MonMultC(poly p, poly q, const ring rr)
2032 {
2033 number x;
2034 poly r = p_Init(rr);
2035
2036 x = n_Mult(pGetCoeff(p),pGetCoeff(q),rr->cf);
2037 pSetCoeff0(r,x);
2038 p_ExpVectorSum(r,p, q, rr);
2039 return r;
2040 }
2041
2042 /*3
2043 * create binomial coef.
2044 */
pnBin(int exp,const ring r)2045 static number* pnBin(int exp, const ring r)
2046 {
2047 int e, i, h;
2048 number x, y, *bin=NULL;
2049
2050 x = n_Init(exp,r->cf);
2051 if (n_IsZero(x,r->cf))
2052 {
2053 n_Delete(&x,r->cf);
2054 return bin;
2055 }
2056 h = (exp >> 1) + 1;
2057 bin = (number *)omAlloc0(h*sizeof(number));
2058 bin[1] = x;
2059 if (exp < 4)
2060 return bin;
2061 i = exp - 1;
2062 for (e=2; e<h; e++)
2063 {
2064 x = n_Init(i,r->cf);
2065 i--;
2066 y = n_Mult(x,bin[e-1],r->cf);
2067 n_Delete(&x,r->cf);
2068 x = n_Init(e,r->cf);
2069 bin[e] = n_ExactDiv(y,x,r->cf);
2070 n_Delete(&x,r->cf);
2071 n_Delete(&y,r->cf);
2072 }
2073 return bin;
2074 }
2075
pnFreeBin(number * bin,int exp,const coeffs r)2076 static void pnFreeBin(number *bin, int exp,const coeffs r)
2077 {
2078 int e, h = (exp >> 1) + 1;
2079
2080 if (bin[1] != NULL)
2081 {
2082 for (e=1; e<h; e++)
2083 n_Delete(&(bin[e]),r);
2084 }
2085 omFreeSize((ADDRESS)bin, h*sizeof(number));
2086 }
2087
2088 /*
2089 * compute for a poly p = head+tail, tail is monomial
2090 * (head + tail)^exp, exp > 1
2091 * with binomial coef.
2092 */
p_TwoMonPower(poly p,int exp,const ring r)2093 static poly p_TwoMonPower(poly p, int exp, const ring r)
2094 {
2095 int eh, e;
2096 long al;
2097 poly *a;
2098 poly tail, b, res, h;
2099 number x;
2100 number *bin = pnBin(exp,r);
2101
2102 tail = pNext(p);
2103 if (bin == NULL)
2104 {
2105 p_MonPower(p,exp,r);
2106 p_MonPower(tail,exp,r);
2107 p_Test(p,r);
2108 return p;
2109 }
2110 eh = exp >> 1;
2111 al = (exp + 1) * sizeof(poly);
2112 a = (poly *)omAlloc(al);
2113 a[1] = p;
2114 for (e=1; e<exp; e++)
2115 {
2116 a[e+1] = p_MonMultC(a[e],p,r);
2117 }
2118 res = a[exp];
2119 b = p_Head(tail,r);
2120 for (e=exp-1; e>eh; e--)
2121 {
2122 h = a[e];
2123 x = n_Mult(bin[exp-e],pGetCoeff(h),r->cf);
2124 p_SetCoeff(h,x,r);
2125 p_MonMult(h,b,r);
2126 res = pNext(res) = h;
2127 p_MonMult(b,tail,r);
2128 }
2129 for (e=eh; e!=0; e--)
2130 {
2131 h = a[e];
2132 x = n_Mult(bin[e],pGetCoeff(h),r->cf);
2133 p_SetCoeff(h,x,r);
2134 p_MonMult(h,b,r);
2135 res = pNext(res) = h;
2136 p_MonMult(b,tail,r);
2137 }
2138 p_LmDelete(&tail,r);
2139 pNext(res) = b;
2140 pNext(b) = NULL;
2141 res = a[exp];
2142 omFreeSize((ADDRESS)a, al);
2143 pnFreeBin(bin, exp, r->cf);
2144 // tail=res;
2145 // while((tail!=NULL)&&(pNext(tail)!=NULL))
2146 // {
2147 // if(nIsZero(pGetCoeff(pNext(tail))))
2148 // {
2149 // pLmDelete(&pNext(tail));
2150 // }
2151 // else
2152 // pIter(tail);
2153 // }
2154 p_Test(res,r);
2155 return res;
2156 }
2157
p_Pow(poly p,int i,const ring r)2158 static poly p_Pow(poly p, int i, const ring r)
2159 {
2160 poly rc = p_Copy(p,r);
2161 i -= 2;
2162 do
2163 {
2164 rc = p_Mult_q(rc,p_Copy(p,r),r);
2165 p_Normalize(rc,r);
2166 i--;
2167 }
2168 while (i != 0);
2169 return p_Mult_q(rc,p,r);
2170 }
2171
p_Pow_charp(poly p,int i,const ring r)2172 static poly p_Pow_charp(poly p, int i, const ring r)
2173 {
2174 //assume char_p == i
2175 poly h=p;
2176 while(h!=NULL) { p_MonPower(h,i,r);pIter(h);}
2177 return p;
2178 }
2179
2180 /*2
2181 * returns the i-th power of p
2182 * p will be destroyed
2183 */
p_Power(poly p,int i,const ring r)2184 poly p_Power(poly p, int i, const ring r)
2185 {
2186 poly rc=NULL;
2187
2188 if (i==0)
2189 {
2190 p_Delete(&p,r);
2191 return p_One(r);
2192 }
2193
2194 if(p!=NULL)
2195 {
2196 if ( (i > 0) && ((unsigned long ) i > (r->bitmask))
2197 #ifdef HAVE_SHIFTBBA
2198 && (!rIsLPRing(r))
2199 #endif
2200 )
2201 {
2202 Werror("exponent %d is too large, max. is %ld",i,r->bitmask);
2203 return NULL;
2204 }
2205 switch (i)
2206 {
2207 // cannot happen, see above
2208 // case 0:
2209 // {
2210 // rc=pOne();
2211 // pDelete(&p);
2212 // break;
2213 // }
2214 case 1:
2215 rc=p;
2216 break;
2217 case 2:
2218 rc=p_Mult_q(p_Copy(p,r),p,r);
2219 break;
2220 default:
2221 if (i < 0)
2222 {
2223 p_Delete(&p,r);
2224 return NULL;
2225 }
2226 else
2227 {
2228 #ifdef HAVE_PLURAL
2229 if (rIsNCRing(r)) /* in the NC case nothing helps :-( */
2230 {
2231 int j=i;
2232 rc = p_Copy(p,r);
2233 while (j>1)
2234 {
2235 rc = p_Mult_q(p_Copy(p,r),rc,r);
2236 j--;
2237 }
2238 p_Delete(&p,r);
2239 return rc;
2240 }
2241 #endif
2242 rc = pNext(p);
2243 if (rc == NULL)
2244 return p_MonPower(p,i,r);
2245 /* else: binom ?*/
2246 int char_p=rChar(r);
2247 if ((char_p>0) && (i>char_p)
2248 && ((rField_is_Zp(r,char_p)
2249 || (rField_is_Zp_a(r,char_p)))))
2250 {
2251 poly h=p_Pow_charp(p_Copy(p,r),char_p,r);
2252 int rest=i-char_p;
2253 while (rest>=char_p)
2254 {
2255 rest-=char_p;
2256 h=p_Mult_q(h,p_Pow_charp(p_Copy(p,r),char_p,r),r);
2257 }
2258 poly res=h;
2259 if (rest>0)
2260 res=p_Mult_q(p_Power(p_Copy(p,r),rest,r),h,r);
2261 p_Delete(&p,r);
2262 return res;
2263 }
2264 if ((pNext(rc) != NULL)
2265 || rField_is_Ring(r)
2266 )
2267 return p_Pow(p,i,r);
2268 if ((char_p==0) || (i<=char_p))
2269 return p_TwoMonPower(p,i,r);
2270 return p_Pow(p,i,r);
2271 }
2272 /*end default:*/
2273 }
2274 }
2275 return rc;
2276 }
2277
2278 /* --------------------------------------------------------------------------------*/
2279 /* content suff */
2280 //number p_InitContent(poly ph, const ring r);
2281
p_Content(poly ph,const ring r)2282 void p_Content(poly ph, const ring r)
2283 {
2284 if (ph==NULL) return;
2285 const coeffs cf=r->cf;
2286 if (pNext(ph)==NULL)
2287 {
2288 p_SetCoeff(ph,n_Init(1,cf),r);
2289 }
2290 if ((cf->cfSubringGcd==ndGcd)
2291 || (cf->cfGcd==ndGcd)) /* trivial gcd*/
2292 return;
2293 number h;
2294 if ((rField_is_Q(r))
2295 || (rField_is_Q_a(r))
2296 || (rField_is_Zp_a)(r)
2297 || (rField_is_Z(r))
2298 )
2299 {
2300 h=p_InitContent(ph,r); /* first guess of a gcd of all coeffs */
2301 }
2302 else
2303 {
2304 h=n_Copy(pGetCoeff(ph),cf);
2305 }
2306 poly p;
2307 if(n_IsOne(h,cf))
2308 {
2309 goto content_finish;
2310 }
2311 p=ph;
2312 // take the SubringGcd of all coeffs
2313 while (p!=NULL)
2314 {
2315 n_Normalize(pGetCoeff(p),cf);
2316 number d=n_SubringGcd(h,pGetCoeff(p),cf);
2317 n_Delete(&h,cf);
2318 h = d;
2319 if(n_IsOne(h,cf))
2320 {
2321 goto content_finish;
2322 }
2323 pIter(p);
2324 }
2325 // if found<>1, divide by it
2326 p = ph;
2327 while (p!=NULL)
2328 {
2329 number d = n_ExactDiv(pGetCoeff(p),h,cf);
2330 p_SetCoeff(p,d,r);
2331 pIter(p);
2332 }
2333 content_finish:
2334 n_Delete(&h,r->cf);
2335 // and last: check leading sign:
2336 if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2337 }
2338
p_Content_n(poly ph,number & c,const ring r)2339 void p_Content_n(poly ph, number &c,const ring r)
2340 {
2341 if (ph==NULL)
2342 {
2343 c=n_Init(1,r->cf);
2344 return;
2345 }
2346 const coeffs cf=r->cf;
2347 if (pNext(ph)==NULL)
2348 {
2349 c=pGetCoeff(ph);
2350 p_SetCoeff0(ph,n_Init(1,cf),r);
2351 }
2352 if ((cf->cfSubringGcd==ndGcd)
2353 || (cf->cfGcd==ndGcd)) /* trivial gcd*/
2354 {
2355 c=n_Init(1,r->cf);
2356 return;
2357 }
2358 number h;
2359 if ((rField_is_Q(r))
2360 || (rField_is_Q_a(r))
2361 || (rField_is_Zp_a)(r)
2362 || (rField_is_Z(r))
2363 )
2364 {
2365 h=p_InitContent(ph,r); /* first guess of a gcd of all coeffs */
2366 }
2367 else
2368 {
2369 h=n_Copy(pGetCoeff(ph),cf);
2370 }
2371 poly p;
2372 if(n_IsOne(h,cf))
2373 {
2374 goto content_finish;
2375 }
2376 p=ph;
2377 // take the SubringGcd of all coeffs
2378 while (p!=NULL)
2379 {
2380 n_Normalize(pGetCoeff(p),cf);
2381 number d=n_SubringGcd(h,pGetCoeff(p),cf);
2382 n_Delete(&h,cf);
2383 h = d;
2384 if(n_IsOne(h,cf))
2385 {
2386 goto content_finish;
2387 }
2388 pIter(p);
2389 }
2390 // if found<>1, divide by it
2391 p = ph;
2392 while (p!=NULL)
2393 {
2394 number d = n_ExactDiv(pGetCoeff(p),h,cf);
2395 p_SetCoeff(p,d,r);
2396 pIter(p);
2397 }
2398 content_finish:
2399 c=h;
2400 // and last: check leading sign:
2401 if(!n_GreaterZero(pGetCoeff(ph),r->cf))
2402 {
2403 c = n_InpNeg(c,r->cf);
2404 ph = p_Neg(ph,r);
2405 }
2406 }
2407
2408 #define CLEARENUMERATORS 1
2409
p_ContentForGB(poly ph,const ring r)2410 void p_ContentForGB(poly ph, const ring r)
2411 {
2412 if(TEST_OPT_CONTENTSB) return;
2413 assume( ph != NULL );
2414
2415 assume( r != NULL ); assume( r->cf != NULL );
2416
2417
2418 #if CLEARENUMERATORS
2419 if( 0 )
2420 {
2421 const coeffs C = r->cf;
2422 // experimentall (recursive enumerator treatment) of alg. Ext!
2423 CPolyCoeffsEnumerator itr(ph);
2424 n_ClearContent(itr, r->cf);
2425
2426 p_Test(ph, r); n_Test(pGetCoeff(ph), C);
2427 assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
2428
2429 // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2430 return;
2431 }
2432 #endif
2433
2434
2435 #ifdef HAVE_RINGS
2436 if (rField_is_Ring(r))
2437 {
2438 if (rField_has_Units(r))
2439 {
2440 number k = n_GetUnit(pGetCoeff(ph),r->cf);
2441 if (!n_IsOne(k,r->cf))
2442 {
2443 number tmpGMP = k;
2444 k = n_Invers(k,r->cf);
2445 n_Delete(&tmpGMP,r->cf);
2446 poly h = pNext(ph);
2447 p_SetCoeff(ph, n_Mult(pGetCoeff(ph), k,r->cf),r);
2448 while (h != NULL)
2449 {
2450 p_SetCoeff(h, n_Mult(pGetCoeff(h), k,r->cf),r);
2451 pIter(h);
2452 }
2453 // assume( n_GreaterZero(pGetCoeff(ph),r->cf) );
2454 // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2455 }
2456 n_Delete(&k,r->cf);
2457 }
2458 return;
2459 }
2460 #endif
2461 number h,d;
2462 poly p;
2463
2464 if(pNext(ph)==NULL)
2465 {
2466 p_SetCoeff(ph,n_Init(1,r->cf),r);
2467 }
2468 else
2469 {
2470 assume( pNext(ph) != NULL );
2471 #if CLEARENUMERATORS
2472 if( nCoeff_is_Q(r->cf) )
2473 {
2474 // experimentall (recursive enumerator treatment) of alg. Ext!
2475 CPolyCoeffsEnumerator itr(ph);
2476 n_ClearContent(itr, r->cf);
2477
2478 p_Test(ph, r); n_Test(pGetCoeff(ph), r->cf);
2479 assume(n_GreaterZero(pGetCoeff(ph), r->cf)); // ??
2480
2481 // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2482 return;
2483 }
2484 #endif
2485
2486 n_Normalize(pGetCoeff(ph),r->cf);
2487 if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2488 if (rField_is_Q(r)||(getCoeffType(r->cf)==n_transExt)) // should not be used anymore if CLEARENUMERATORS is 1
2489 {
2490 h=p_InitContent(ph,r);
2491 p=ph;
2492 }
2493 else
2494 {
2495 h=n_Copy(pGetCoeff(ph),r->cf);
2496 p = pNext(ph);
2497 }
2498 while (p!=NULL)
2499 {
2500 n_Normalize(pGetCoeff(p),r->cf);
2501 d=n_SubringGcd(h,pGetCoeff(p),r->cf);
2502 n_Delete(&h,r->cf);
2503 h = d;
2504 if(n_IsOne(h,r->cf))
2505 {
2506 break;
2507 }
2508 pIter(p);
2509 }
2510 //number tmp;
2511 if(!n_IsOne(h,r->cf))
2512 {
2513 p = ph;
2514 while (p!=NULL)
2515 {
2516 //d = nDiv(pGetCoeff(p),h);
2517 //tmp = nExactDiv(pGetCoeff(p),h);
2518 //if (!nEqual(d,tmp))
2519 //{
2520 // StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
2521 // nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
2522 // nWrite(tmp);Print(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
2523 //}
2524 //nDelete(&tmp);
2525 d = n_ExactDiv(pGetCoeff(p),h,r->cf);
2526 p_SetCoeff(p,d,r);
2527 pIter(p);
2528 }
2529 }
2530 n_Delete(&h,r->cf);
2531 if (rField_is_Q_a(r))
2532 {
2533 // special handling for alg. ext.:
2534 if (getCoeffType(r->cf)==n_algExt)
2535 {
2536 h = n_Init(1, r->cf->extRing->cf);
2537 p=ph;
2538 while (p!=NULL)
2539 { // each monom: coeff in Q_a
2540 poly c_n_n=(poly)pGetCoeff(p);
2541 poly c_n=c_n_n;
2542 while (c_n!=NULL)
2543 { // each monom: coeff in Q
2544 d=n_NormalizeHelper(h,pGetCoeff(c_n),r->cf->extRing->cf);
2545 n_Delete(&h,r->cf->extRing->cf);
2546 h=d;
2547 pIter(c_n);
2548 }
2549 pIter(p);
2550 }
2551 /* h contains the 1/lcm of all denominators in c_n_n*/
2552 //n_Normalize(h,r->cf->extRing->cf);
2553 if(!n_IsOne(h,r->cf->extRing->cf))
2554 {
2555 p=ph;
2556 while (p!=NULL)
2557 { // each monom: coeff in Q_a
2558 poly c_n=(poly)pGetCoeff(p);
2559 while (c_n!=NULL)
2560 { // each monom: coeff in Q
2561 d=n_Mult(h,pGetCoeff(c_n),r->cf->extRing->cf);
2562 n_Normalize(d,r->cf->extRing->cf);
2563 n_Delete(&pGetCoeff(c_n),r->cf->extRing->cf);
2564 pGetCoeff(c_n)=d;
2565 pIter(c_n);
2566 }
2567 pIter(p);
2568 }
2569 }
2570 n_Delete(&h,r->cf->extRing->cf);
2571 }
2572 /*else
2573 {
2574 // special handling for rat. functions.:
2575 number hzz =NULL;
2576 p=ph;
2577 while (p!=NULL)
2578 { // each monom: coeff in Q_a (Z_a)
2579 fraction f=(fraction)pGetCoeff(p);
2580 poly c_n=NUM(f);
2581 if (hzz==NULL)
2582 {
2583 hzz=n_Copy(pGetCoeff(c_n),r->cf->extRing->cf);
2584 pIter(c_n);
2585 }
2586 while ((c_n!=NULL)&&(!n_IsOne(hzz,r->cf->extRing->cf)))
2587 { // each monom: coeff in Q (Z)
2588 d=n_Gcd(hzz,pGetCoeff(c_n),r->cf->extRing->cf);
2589 n_Delete(&hzz,r->cf->extRing->cf);
2590 hzz=d;
2591 pIter(c_n);
2592 }
2593 pIter(p);
2594 }
2595 // hzz contains the gcd of all numerators in f
2596 h=n_Invers(hzz,r->cf->extRing->cf);
2597 n_Delete(&hzz,r->cf->extRing->cf);
2598 n_Normalize(h,r->cf->extRing->cf);
2599 if(!n_IsOne(h,r->cf->extRing->cf))
2600 {
2601 p=ph;
2602 while (p!=NULL)
2603 { // each monom: coeff in Q_a (Z_a)
2604 fraction f=(fraction)pGetCoeff(p);
2605 NUM(f)=__p_Mult_nn(NUM(f),h,r->cf->extRing);
2606 p_Normalize(NUM(f),r->cf->extRing);
2607 pIter(p);
2608 }
2609 }
2610 n_Delete(&h,r->cf->extRing->cf);
2611 }*/
2612 }
2613 }
2614 if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2615 }
2616
2617 // Not yet?
2618 #if 1 // currently only used by Singular/janet
p_SimpleContent(poly ph,int smax,const ring r)2619 void p_SimpleContent(poly ph, int smax, const ring r)
2620 {
2621 if(TEST_OPT_CONTENTSB) return;
2622 if (ph==NULL) return;
2623 if (pNext(ph)==NULL)
2624 {
2625 p_SetCoeff(ph,n_Init(1,r->cf),r);
2626 return;
2627 }
2628 if (pNext(pNext(ph))==NULL)
2629 {
2630 return;
2631 }
2632 if (!(rField_is_Q(r))
2633 && (!rField_is_Q_a(r))
2634 && (!rField_is_Zp_a(r))
2635 && (!rField_is_Z(r))
2636 )
2637 {
2638 return;
2639 }
2640 number d=p_InitContent(ph,r);
2641 number h=d;
2642 if (n_Size(d,r->cf)<=smax)
2643 {
2644 n_Delete(&h,r->cf);
2645 //if (TEST_OPT_PROT) PrintS("G");
2646 return;
2647 }
2648
2649 poly p=ph;
2650 if (smax==1) smax=2;
2651 while (p!=NULL)
2652 {
2653 #if 1
2654 d=n_SubringGcd(h,pGetCoeff(p),r->cf);
2655 n_Delete(&h,r->cf);
2656 h = d;
2657 #else
2658 n_InpGcd(h,pGetCoeff(p),r->cf);
2659 #endif
2660 if(n_Size(h,r->cf)<smax)
2661 {
2662 //if (TEST_OPT_PROT) PrintS("g");
2663 n_Delete(&h,r->cf);
2664 return;
2665 }
2666 pIter(p);
2667 }
2668 p = ph;
2669 if (!n_GreaterZero(pGetCoeff(p),r->cf)) h=n_InpNeg(h,r->cf);
2670 if(n_IsOne(h,r->cf))
2671 {
2672 n_Delete(&h,r->cf);
2673 return;
2674 }
2675 if (TEST_OPT_PROT) PrintS("c");
2676 while (p!=NULL)
2677 {
2678 #if 1
2679 d = n_ExactDiv(pGetCoeff(p),h,r->cf);
2680 p_SetCoeff(p,d,r);
2681 #else
2682 STATISTIC(n_ExactDiv); nlInpExactDiv(pGetCoeff(p),h,r->cf); // no such function... ?
2683 #endif
2684 pIter(p);
2685 }
2686 n_Delete(&h,r->cf);
2687 }
2688 #endif
2689
p_InitContent(poly ph,const ring r)2690 number p_InitContent(poly ph, const ring r)
2691 // only for coefficients in Q and rational functions
2692 #if 0
2693 {
2694 assume(!TEST_OPT_CONTENTSB);
2695 assume(ph!=NULL);
2696 assume(pNext(ph)!=NULL);
2697 assume(rField_is_Q(r));
2698 if (pNext(pNext(ph))==NULL)
2699 {
2700 return n_GetNumerator(pGetCoeff(pNext(ph)),r->cf);
2701 }
2702 poly p=ph;
2703 number n1=n_GetNumerator(pGetCoeff(p),r->cf);
2704 pIter(p);
2705 number n2=n_GetNumerator(pGetCoeff(p),r->cf);
2706 pIter(p);
2707 number d;
2708 number t;
2709 loop
2710 {
2711 nlNormalize(pGetCoeff(p),r->cf);
2712 t=n_GetNumerator(pGetCoeff(p),r->cf);
2713 if (nlGreaterZero(t,r->cf))
2714 d=nlAdd(n1,t,r->cf);
2715 else
2716 d=nlSub(n1,t,r->cf);
2717 nlDelete(&t,r->cf);
2718 nlDelete(&n1,r->cf);
2719 n1=d;
2720 pIter(p);
2721 if (p==NULL) break;
2722 nlNormalize(pGetCoeff(p),r->cf);
2723 t=n_GetNumerator(pGetCoeff(p),r->cf);
2724 if (nlGreaterZero(t,r->cf))
2725 d=nlAdd(n2,t,r->cf);
2726 else
2727 d=nlSub(n2,t,r->cf);
2728 nlDelete(&t,r->cf);
2729 nlDelete(&n2,r->cf);
2730 n2=d;
2731 pIter(p);
2732 if (p==NULL) break;
2733 }
2734 d=nlGcd(n1,n2,r->cf);
2735 nlDelete(&n1,r->cf);
2736 nlDelete(&n2,r->cf);
2737 return d;
2738 }
2739 #else
2740 {
2741 /* ph has al least 2 terms */
2742 number d=pGetCoeff(ph);
2743 int s=n_Size(d,r->cf);
2744 pIter(ph);
2745 number d2=pGetCoeff(ph);
2746 int s2=n_Size(d2,r->cf);
2747 pIter(ph);
2748 if (ph==NULL)
2749 {
2750 if (s<s2) return n_Copy(d,r->cf);
2751 else return n_Copy(d2,r->cf);
2752 }
2753 do
2754 {
2755 number nd=pGetCoeff(ph);
2756 int ns=n_Size(nd,r->cf);
2757 if (ns<=2)
2758 {
2759 s2=s;
2760 d2=d;
2761 d=nd;
2762 s=ns;
2763 break;
2764 }
2765 else if (ns<s)
2766 {
2767 s2=s;
2768 d2=d;
2769 d=nd;
2770 s=ns;
2771 }
2772 pIter(ph);
2773 }
2774 while(ph!=NULL);
2775 return n_SubringGcd(d,d2,r->cf);
2776 }
2777 #endif
2778
2779 //void pContent(poly ph)
2780 //{
2781 // number h,d;
2782 // poly p;
2783 //
2784 // p = ph;
2785 // if(pNext(p)==NULL)
2786 // {
2787 // pSetCoeff(p,nInit(1));
2788 // }
2789 // else
2790 // {
2791 //#ifdef PDEBUG
2792 // if (!pTest(p)) return;
2793 //#endif
2794 // nNormalize(pGetCoeff(p));
2795 // if(!nGreaterZero(pGetCoeff(ph)))
2796 // {
2797 // ph = pNeg(ph);
2798 // nNormalize(pGetCoeff(p));
2799 // }
2800 // h=pGetCoeff(p);
2801 // pIter(p);
2802 // while (p!=NULL)
2803 // {
2804 // nNormalize(pGetCoeff(p));
2805 // if (nGreater(h,pGetCoeff(p))) h=pGetCoeff(p);
2806 // pIter(p);
2807 // }
2808 // h=nCopy(h);
2809 // p=ph;
2810 // while (p!=NULL)
2811 // {
2812 // d=n_Gcd(h,pGetCoeff(p));
2813 // nDelete(&h);
2814 // h = d;
2815 // if(nIsOne(h))
2816 // {
2817 // break;
2818 // }
2819 // pIter(p);
2820 // }
2821 // p = ph;
2822 // //number tmp;
2823 // if(!nIsOne(h))
2824 // {
2825 // while (p!=NULL)
2826 // {
2827 // d = nExactDiv(pGetCoeff(p),h);
2828 // pSetCoeff(p,d);
2829 // pIter(p);
2830 // }
2831 // }
2832 // nDelete(&h);
2833 // if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
2834 // {
2835 // pTest(ph);
2836 // singclap_divide_content(ph);
2837 // pTest(ph);
2838 // }
2839 // }
2840 //}
2841 #if 0
2842 void p_Content(poly ph, const ring r)
2843 {
2844 number h,d;
2845 poly p;
2846
2847 if(pNext(ph)==NULL)
2848 {
2849 p_SetCoeff(ph,n_Init(1,r->cf),r);
2850 }
2851 else
2852 {
2853 n_Normalize(pGetCoeff(ph),r->cf);
2854 if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2855 h=n_Copy(pGetCoeff(ph),r->cf);
2856 p = pNext(ph);
2857 while (p!=NULL)
2858 {
2859 n_Normalize(pGetCoeff(p),r->cf);
2860 d=n_Gcd(h,pGetCoeff(p),r->cf);
2861 n_Delete(&h,r->cf);
2862 h = d;
2863 if(n_IsOne(h,r->cf))
2864 {
2865 break;
2866 }
2867 pIter(p);
2868 }
2869 p = ph;
2870 //number tmp;
2871 if(!n_IsOne(h,r->cf))
2872 {
2873 while (p!=NULL)
2874 {
2875 //d = nDiv(pGetCoeff(p),h);
2876 //tmp = nExactDiv(pGetCoeff(p),h);
2877 //if (!nEqual(d,tmp))
2878 //{
2879 // StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
2880 // nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
2881 // nWrite(tmp);Print(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
2882 //}
2883 //nDelete(&tmp);
2884 d = n_ExactDiv(pGetCoeff(p),h,r->cf);
2885 p_SetCoeff(p,d,r->cf);
2886 pIter(p);
2887 }
2888 }
2889 n_Delete(&h,r->cf);
2890 //if ( (n_GetChar(r) == 1) || (n_GetChar(r) < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
2891 //{
2892 // singclap_divide_content(ph);
2893 // if(!n_GreaterZero(pGetCoeff(ph),r)) ph = p_Neg(ph,r);
2894 //}
2895 }
2896 }
2897 #endif
2898 /* ---------------------------------------------------------------------------*/
2899 /* cleardenom suff */
p_Cleardenom(poly p,const ring r)2900 poly p_Cleardenom(poly p, const ring r)
2901 {
2902 if( p == NULL )
2903 return NULL;
2904
2905 assume( r != NULL ); assume( r->cf != NULL ); const coeffs C = r->cf;
2906
2907 #if CLEARENUMERATORS
2908 if( 0 )
2909 {
2910 CPolyCoeffsEnumerator itr(p);
2911 n_ClearDenominators(itr, C);
2912 n_ClearContent(itr, C); // divide out the content
2913 p_Test(p, r); n_Test(pGetCoeff(p), C);
2914 assume(n_GreaterZero(pGetCoeff(p), C)); // ??
2915 // if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2916 return p;
2917 }
2918 #endif
2919
2920 number d, h;
2921
2922 if (rField_is_Ring(r))
2923 {
2924 p_ContentForGB(p,r);
2925 if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2926 return p;
2927 }
2928
2929 if (rField_is_Zp(r) && TEST_OPT_INTSTRATEGY)
2930 {
2931 if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2932 return p;
2933 }
2934
2935 assume(p != NULL);
2936
2937 if(pNext(p)==NULL)
2938 {
2939 if (!TEST_OPT_CONTENTSB
2940 && !rField_is_Ring(r))
2941 p_SetCoeff(p,n_Init(1,r->cf),r);
2942 else if(!n_GreaterZero(pGetCoeff(p),C))
2943 p = p_Neg(p,r);
2944 return p;
2945 }
2946
2947 assume(pNext(p)!=NULL);
2948 poly start=p;
2949
2950 #if 0 && CLEARENUMERATORS
2951 //CF: does not seem to work that well..
2952
2953 if( nCoeff_is_Q(C) || nCoeff_is_Q_a(C) )
2954 {
2955 CPolyCoeffsEnumerator itr(p);
2956 n_ClearDenominators(itr, C);
2957 n_ClearContent(itr, C); // divide out the content
2958 p_Test(p, r); n_Test(pGetCoeff(p), C);
2959 assume(n_GreaterZero(pGetCoeff(p), C)); // ??
2960 // if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2961 return start;
2962 }
2963 #endif
2964
2965 if(1)
2966 {
2967 // get lcm of all denominators ----------------------------------
2968 h = n_Init(1,r->cf);
2969 while (p!=NULL)
2970 {
2971 n_Normalize(pGetCoeff(p),r->cf);
2972 d=n_NormalizeHelper(h,pGetCoeff(p),r->cf);
2973 n_Delete(&h,r->cf);
2974 h=d;
2975 pIter(p);
2976 }
2977 /* h now contains the 1/lcm of all denominators */
2978 if(!n_IsOne(h,r->cf))
2979 {
2980 // multiply by the lcm of all denominators
2981 p = start;
2982 while (p!=NULL)
2983 {
2984 d=n_Mult(h,pGetCoeff(p),r->cf);
2985 n_Normalize(d,r->cf);
2986 p_SetCoeff(p,d,r);
2987 pIter(p);
2988 }
2989 }
2990 n_Delete(&h,r->cf);
2991 p=start;
2992
2993 p_ContentForGB(p,r);
2994 #ifdef HAVE_RATGRING
2995 if (rIsRatGRing(r))
2996 {
2997 /* quick unit detection in the rational case is done in gr_nc_bba */
2998 p_ContentRat(p, r);
2999 start=p;
3000 }
3001 #endif
3002 }
3003
3004 if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
3005
3006 return start;
3007 }
3008
p_Cleardenom_n(poly ph,const ring r,number & c)3009 void p_Cleardenom_n(poly ph,const ring r,number &c)
3010 {
3011 const coeffs C = r->cf;
3012 number d, h;
3013
3014 assume( ph != NULL );
3015
3016 poly p = ph;
3017
3018 #if CLEARENUMERATORS
3019 if( 0 )
3020 {
3021 CPolyCoeffsEnumerator itr(ph);
3022
3023 n_ClearDenominators(itr, d, C); // multiply with common denom. d
3024 n_ClearContent(itr, h, C); // divide by the content h
3025
3026 c = n_Div(d, h, C); // d/h
3027
3028 n_Delete(&d, C);
3029 n_Delete(&h, C);
3030
3031 n_Test(c, C);
3032
3033 p_Test(ph, r); n_Test(pGetCoeff(ph), C);
3034 assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
3035 /*
3036 if(!n_GreaterZero(pGetCoeff(ph),C))
3037 {
3038 ph = p_Neg(ph,r);
3039 c = n_InpNeg(c, C);
3040 }
3041 */
3042 return;
3043 }
3044 #endif
3045
3046
3047 if( pNext(p) == NULL )
3048 {
3049 if(!TEST_OPT_CONTENTSB)
3050 {
3051 c=n_Invers(pGetCoeff(p), C);
3052 p_SetCoeff(p, n_Init(1, C), r);
3053 }
3054 else
3055 {
3056 c=n_Init(1,C);
3057 }
3058
3059 if(!n_GreaterZero(pGetCoeff(ph),C))
3060 {
3061 ph = p_Neg(ph,r);
3062 c = n_InpNeg(c, C);
3063 }
3064
3065 return;
3066 }
3067 if (TEST_OPT_CONTENTSB) { c=n_Init(1,C); return; }
3068
3069 assume( pNext(p) != NULL );
3070
3071 #if CLEARENUMERATORS
3072 if( nCoeff_is_Q(C) || nCoeff_is_Q_a(C) )
3073 {
3074 CPolyCoeffsEnumerator itr(ph);
3075
3076 n_ClearDenominators(itr, d, C); // multiply with common denom. d
3077 n_ClearContent(itr, h, C); // divide by the content h
3078
3079 c = n_Div(d, h, C); // d/h
3080
3081 n_Delete(&d, C);
3082 n_Delete(&h, C);
3083
3084 n_Test(c, C);
3085
3086 p_Test(ph, r); n_Test(pGetCoeff(ph), C);
3087 assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
3088 /*
3089 if(!n_GreaterZero(pGetCoeff(ph),C))
3090 {
3091 ph = p_Neg(ph,r);
3092 c = n_InpNeg(c, C);
3093 }
3094 */
3095 return;
3096 }
3097 #endif
3098
3099
3100
3101
3102 if(1)
3103 {
3104 h = n_Init(1,r->cf);
3105 while (p!=NULL)
3106 {
3107 n_Normalize(pGetCoeff(p),r->cf);
3108 d=n_NormalizeHelper(h,pGetCoeff(p),r->cf);
3109 n_Delete(&h,r->cf);
3110 h=d;
3111 pIter(p);
3112 }
3113 c=h;
3114 /* contains the 1/lcm of all denominators */
3115 if(!n_IsOne(h,r->cf))
3116 {
3117 p = ph;
3118 while (p!=NULL)
3119 {
3120 /* should be: // NOTE: don't use ->coef!!!!
3121 * number hh;
3122 * nGetDenom(p->coef,&hh);
3123 * nMult(&h,&hh,&d);
3124 * nNormalize(d);
3125 * nDelete(&hh);
3126 * nMult(d,p->coef,&hh);
3127 * nDelete(&d);
3128 * nDelete(&(p->coef));
3129 * p->coef =hh;
3130 */
3131 d=n_Mult(h,pGetCoeff(p),r->cf);
3132 n_Normalize(d,r->cf);
3133 p_SetCoeff(p,d,r);
3134 pIter(p);
3135 }
3136 if (rField_is_Q_a(r))
3137 {
3138 loop
3139 {
3140 h = n_Init(1,r->cf);
3141 p=ph;
3142 while (p!=NULL)
3143 {
3144 d=n_NormalizeHelper(h,pGetCoeff(p),r->cf);
3145 n_Delete(&h,r->cf);
3146 h=d;
3147 pIter(p);
3148 }
3149 /* contains the 1/lcm of all denominators */
3150 if(!n_IsOne(h,r->cf))
3151 {
3152 p = ph;
3153 while (p!=NULL)
3154 {
3155 /* should be: // NOTE: don't use ->coef!!!!
3156 * number hh;
3157 * nGetDenom(p->coef,&hh);
3158 * nMult(&h,&hh,&d);
3159 * nNormalize(d);
3160 * nDelete(&hh);
3161 * nMult(d,p->coef,&hh);
3162 * nDelete(&d);
3163 * nDelete(&(p->coef));
3164 * p->coef =hh;
3165 */
3166 d=n_Mult(h,pGetCoeff(p),r->cf);
3167 n_Normalize(d,r->cf);
3168 p_SetCoeff(p,d,r);
3169 pIter(p);
3170 }
3171 number t=n_Mult(c,h,r->cf);
3172 n_Delete(&c,r->cf);
3173 c=t;
3174 }
3175 else
3176 {
3177 break;
3178 }
3179 n_Delete(&h,r->cf);
3180 }
3181 }
3182 }
3183 }
3184
3185 if(!n_GreaterZero(pGetCoeff(ph),C))
3186 {
3187 ph = p_Neg(ph,r);
3188 c = n_InpNeg(c, C);
3189 }
3190
3191 }
3192
3193 // normalization: for poly over Q: make poly primitive, integral
3194 // Qa make poly integral with leading
3195 // coefficient minimal in N
3196 // Q(t) make poly primitive, integral
3197
p_ProjectiveUnique(poly ph,const ring r)3198 void p_ProjectiveUnique(poly ph, const ring r)
3199 {
3200 if( ph == NULL )
3201 return;
3202
3203 assume( r != NULL ); assume( r->cf != NULL );
3204 const coeffs C = r->cf;
3205
3206 number h;
3207 poly p;
3208
3209 if (nCoeff_is_Ring(C))
3210 {
3211 p_ContentForGB(ph,r);
3212 if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
3213 assume( n_GreaterZero(pGetCoeff(ph),C) );
3214 return;
3215 }
3216
3217 if (nCoeff_is_Zp(C) && TEST_OPT_INTSTRATEGY)
3218 {
3219 assume( n_GreaterZero(pGetCoeff(ph),C) );
3220 if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
3221 return;
3222 }
3223 p = ph;
3224
3225 assume(p != NULL);
3226
3227 if(pNext(p)==NULL) // a monomial
3228 {
3229 p_SetCoeff(p, n_Init(1, C), r);
3230 return;
3231 }
3232
3233 assume(pNext(p)!=NULL);
3234
3235 if(!nCoeff_is_Q(C) && !nCoeff_is_transExt(C))
3236 {
3237 h = p_GetCoeff(p, C);
3238 number hInv = n_Invers(h, C);
3239 pIter(p);
3240 while (p!=NULL)
3241 {
3242 p_SetCoeff(p, n_Mult(p_GetCoeff(p, C), hInv, C), r);
3243 pIter(p);
3244 }
3245 n_Delete(&hInv, C);
3246 p = ph;
3247 p_SetCoeff(p, n_Init(1, C), r);
3248 }
3249
3250 p_Cleardenom(ph, r); //removes also Content
3251
3252
3253 /* normalize ph over a transcendental extension s.t.
3254 lead (ph) is > 0 if extRing->cf == Q
3255 or lead (ph) is monic if extRing->cf == Zp*/
3256 if (nCoeff_is_transExt(C))
3257 {
3258 p= ph;
3259 h= p_GetCoeff (p, C);
3260 fraction f = (fraction) h;
3261 number n=p_GetCoeff (NUM (f),C->extRing->cf);
3262 if (rField_is_Q (C->extRing))
3263 {
3264 if (!n_GreaterZero(n,C->extRing->cf))
3265 {
3266 p=p_Neg (p,r);
3267 }
3268 }
3269 else if (rField_is_Zp(C->extRing))
3270 {
3271 if (!n_IsOne (n, C->extRing->cf))
3272 {
3273 n=n_Invers (n,C->extRing->cf);
3274 nMapFunc nMap;
3275 nMap= n_SetMap (C->extRing->cf, C);
3276 number ninv= nMap (n,C->extRing->cf, C);
3277 p=__p_Mult_nn (p, ninv, r);
3278 n_Delete (&ninv, C);
3279 n_Delete (&n, C->extRing->cf);
3280 }
3281 }
3282 p= ph;
3283 }
3284
3285 return;
3286 }
3287
3288 #if 0 /*unused*/
3289 number p_GetAllDenom(poly ph, const ring r)
3290 {
3291 number d=n_Init(1,r->cf);
3292 poly p = ph;
3293
3294 while (p!=NULL)
3295 {
3296 number h=n_GetDenom(pGetCoeff(p),r->cf);
3297 if (!n_IsOne(h,r->cf))
3298 {
3299 number dd=n_Mult(d,h,r->cf);
3300 n_Delete(&d,r->cf);
3301 d=dd;
3302 }
3303 n_Delete(&h,r->cf);
3304 pIter(p);
3305 }
3306 return d;
3307 }
3308 #endif
3309
p_Size(poly p,const ring r)3310 int p_Size(poly p, const ring r)
3311 {
3312 int count = 0;
3313 if (r->cf->has_simple_Alloc)
3314 return pLength(p);
3315 while ( p != NULL )
3316 {
3317 count+= n_Size( pGetCoeff( p ), r->cf );
3318 pIter( p );
3319 }
3320 return count;
3321 }
3322
3323 /*2
3324 *make p homogeneous by multiplying the monomials by powers of x_varnum
3325 *assume: deg(var(varnum))==1
3326 */
p_Homogen(poly p,int varnum,const ring r)3327 poly p_Homogen (poly p, int varnum, const ring r)
3328 {
3329 pFDegProc deg;
3330 if (r->pLexOrder && (r->order[0]==ringorder_lp))
3331 deg=p_Totaldegree;
3332 else
3333 deg=r->pFDeg;
3334
3335 poly q=NULL, qn;
3336 int o,ii;
3337 sBucket_pt bp;
3338
3339 if (p!=NULL)
3340 {
3341 if ((varnum < 1) || (varnum > rVar(r)))
3342 {
3343 return NULL;
3344 }
3345 o=deg(p,r);
3346 q=pNext(p);
3347 while (q != NULL)
3348 {
3349 ii=deg(q,r);
3350 if (ii>o) o=ii;
3351 pIter(q);
3352 }
3353 q = p_Copy(p,r);
3354 bp = sBucketCreate(r);
3355 while (q != NULL)
3356 {
3357 ii = o-deg(q,r);
3358 if (ii!=0)
3359 {
3360 p_AddExp(q,varnum, (long)ii,r);
3361 p_Setm(q,r);
3362 }
3363 qn = pNext(q);
3364 pNext(q) = NULL;
3365 sBucket_Add_m(bp, q);
3366 q = qn;
3367 }
3368 sBucketDestroyAdd(bp, &q, &ii);
3369 }
3370 return q;
3371 }
3372
3373 /*2
3374 *tests if p is homogeneous with respect to the actual weigths
3375 */
p_IsHomogeneous(poly p,const ring r)3376 BOOLEAN p_IsHomogeneous (poly p, const ring r)
3377 {
3378 poly qp=p;
3379 int o;
3380
3381 if ((p == NULL) || (pNext(p) == NULL)) return TRUE;
3382 pFDegProc d;
3383 if (r->pLexOrder && (r->order[0]==ringorder_lp))
3384 d=p_Totaldegree;
3385 else
3386 d=r->pFDeg;
3387 o = d(p,r);
3388 do
3389 {
3390 if (d(qp,r) != o) return FALSE;
3391 pIter(qp);
3392 }
3393 while (qp != NULL);
3394 return TRUE;
3395 }
3396
3397 /*----------utilities for syzygies--------------*/
p_VectorHasUnitB(poly p,int * k,const ring r)3398 BOOLEAN p_VectorHasUnitB(poly p, int * k, const ring r)
3399 {
3400 poly q=p,qq;
3401 int i;
3402
3403 while (q!=NULL)
3404 {
3405 if (p_LmIsConstantComp(q,r))
3406 {
3407 i = __p_GetComp(q,r);
3408 qq = p;
3409 while ((qq != q) && (__p_GetComp(qq,r) != i)) pIter(qq);
3410 if (qq == q)
3411 {
3412 *k = i;
3413 return TRUE;
3414 }
3415 }
3416 pIter(q);
3417 }
3418 return FALSE;
3419 }
3420
p_VectorHasUnit(poly p,int * k,int * len,const ring r)3421 void p_VectorHasUnit(poly p, int * k, int * len, const ring r)
3422 {
3423 poly q=p,qq;
3424 int i,j=0;
3425
3426 *len = 0;
3427 while (q!=NULL)
3428 {
3429 if (p_LmIsConstantComp(q,r))
3430 {
3431 i = __p_GetComp(q,r);
3432 qq = p;
3433 while ((qq != q) && (__p_GetComp(qq,r) != i)) pIter(qq);
3434 if (qq == q)
3435 {
3436 j = 0;
3437 while (qq!=NULL)
3438 {
3439 if (__p_GetComp(qq,r)==i) j++;
3440 pIter(qq);
3441 }
3442 if ((*len == 0) || (j<*len))
3443 {
3444 *len = j;
3445 *k = i;
3446 }
3447 }
3448 }
3449 pIter(q);
3450 }
3451 }
3452
p_TakeOutComp1(poly * p,int k,const ring r)3453 poly p_TakeOutComp1(poly * p, int k, const ring r)
3454 {
3455 poly q = *p;
3456
3457 if (q==NULL) return NULL;
3458
3459 poly qq=NULL,result = NULL;
3460
3461 if (__p_GetComp(q,r)==k)
3462 {
3463 result = q; /* *p */
3464 while ((q!=NULL) && (__p_GetComp(q,r)==k))
3465 {
3466 p_SetComp(q,0,r);
3467 p_SetmComp(q,r);
3468 qq = q;
3469 pIter(q);
3470 }
3471 *p = q;
3472 pNext(qq) = NULL;
3473 }
3474 if (q==NULL) return result;
3475 // if (pGetComp(q) > k) pGetComp(q)--;
3476 while (pNext(q)!=NULL)
3477 {
3478 if (__p_GetComp(pNext(q),r)==k)
3479 {
3480 if (result==NULL)
3481 {
3482 result = pNext(q);
3483 qq = result;
3484 }
3485 else
3486 {
3487 pNext(qq) = pNext(q);
3488 pIter(qq);
3489 }
3490 pNext(q) = pNext(pNext(q));
3491 pNext(qq) =NULL;
3492 p_SetComp(qq,0,r);
3493 p_SetmComp(qq,r);
3494 }
3495 else
3496 {
3497 pIter(q);
3498 // if (pGetComp(q) > k) pGetComp(q)--;
3499 }
3500 }
3501 return result;
3502 }
3503
p_TakeOutComp(poly * p,int k,const ring r)3504 poly p_TakeOutComp(poly * p, int k, const ring r)
3505 {
3506 poly q = *p,qq=NULL,result = NULL;
3507
3508 if (q==NULL) return NULL;
3509 BOOLEAN use_setmcomp=rOrd_SetCompRequiresSetm(r);
3510 if (__p_GetComp(q,r)==k)
3511 {
3512 result = q;
3513 do
3514 {
3515 p_SetComp(q,0,r);
3516 if (use_setmcomp) p_SetmComp(q,r);
3517 qq = q;
3518 pIter(q);
3519 }
3520 while ((q!=NULL) && (__p_GetComp(q,r)==k));
3521 *p = q;
3522 pNext(qq) = NULL;
3523 }
3524 if (q==NULL) return result;
3525 if (__p_GetComp(q,r) > k)
3526 {
3527 p_SubComp(q,1,r);
3528 if (use_setmcomp) p_SetmComp(q,r);
3529 }
3530 poly pNext_q;
3531 while ((pNext_q=pNext(q))!=NULL)
3532 {
3533 if (__p_GetComp(pNext_q,r)==k)
3534 {
3535 if (result==NULL)
3536 {
3537 result = pNext_q;
3538 qq = result;
3539 }
3540 else
3541 {
3542 pNext(qq) = pNext_q;
3543 pIter(qq);
3544 }
3545 pNext(q) = pNext(pNext_q);
3546 pNext(qq) =NULL;
3547 p_SetComp(qq,0,r);
3548 if (use_setmcomp) p_SetmComp(qq,r);
3549 }
3550 else
3551 {
3552 /*pIter(q);*/ q=pNext_q;
3553 if (__p_GetComp(q,r) > k)
3554 {
3555 p_SubComp(q,1,r);
3556 if (use_setmcomp) p_SetmComp(q,r);
3557 }
3558 }
3559 }
3560 return result;
3561 }
3562
3563 // Splits *p into two polys: *q which consists of all monoms with
3564 // component == comp and *p of all other monoms *lq == pLength(*q)
p_TakeOutComp(poly * r_p,long comp,poly * r_q,int * lq,const ring r)3565 void p_TakeOutComp(poly *r_p, long comp, poly *r_q, int *lq, const ring r)
3566 {
3567 spolyrec pp, qq;
3568 poly p, q, p_prev;
3569 int l = 0;
3570
3571 #ifndef SING_NDEBUG
3572 int lp = pLength(*r_p);
3573 #endif
3574
3575 pNext(&pp) = *r_p;
3576 p = *r_p;
3577 p_prev = &pp;
3578 q = &qq;
3579
3580 while(p != NULL)
3581 {
3582 while (__p_GetComp(p,r) == comp)
3583 {
3584 pNext(q) = p;
3585 pIter(q);
3586 p_SetComp(p, 0,r);
3587 p_SetmComp(p,r);
3588 pIter(p);
3589 l++;
3590 if (p == NULL)
3591 {
3592 pNext(p_prev) = NULL;
3593 goto Finish;
3594 }
3595 }
3596 pNext(p_prev) = p;
3597 p_prev = p;
3598 pIter(p);
3599 }
3600
3601 Finish:
3602 pNext(q) = NULL;
3603 *r_p = pNext(&pp);
3604 *r_q = pNext(&qq);
3605 *lq = l;
3606 #ifndef SING_NDEBUG
3607 assume(pLength(*r_p) + pLength(*r_q) == lp);
3608 #endif
3609 p_Test(*r_p,r);
3610 p_Test(*r_q,r);
3611 }
3612
p_DeleteComp(poly * p,int k,const ring r)3613 void p_DeleteComp(poly * p,int k, const ring r)
3614 {
3615 poly q;
3616
3617 while ((*p!=NULL) && (__p_GetComp(*p,r)==k)) p_LmDelete(p,r);
3618 if (*p==NULL) return;
3619 q = *p;
3620 if (__p_GetComp(q,r)>k)
3621 {
3622 p_SubComp(q,1,r);
3623 p_SetmComp(q,r);
3624 }
3625 while (pNext(q)!=NULL)
3626 {
3627 if (__p_GetComp(pNext(q),r)==k)
3628 p_LmDelete(&(pNext(q)),r);
3629 else
3630 {
3631 pIter(q);
3632 if (__p_GetComp(q,r)>k)
3633 {
3634 p_SubComp(q,1,r);
3635 p_SetmComp(q,r);
3636 }
3637 }
3638 }
3639 }
3640
p_Vec2Poly(poly v,int k,const ring r)3641 poly p_Vec2Poly(poly v, int k, const ring r)
3642 {
3643 poly h;
3644 poly res=NULL;
3645
3646 while (v!=NULL)
3647 {
3648 if (__p_GetComp(v,r)==k)
3649 {
3650 h=p_Head(v,r);
3651 p_SetComp(h,0,r);
3652 pNext(h)=res;res=h;
3653 }
3654 pIter(v);
3655 }
3656 if (res!=NULL) res=pReverse(res);
3657 return res;
3658 }
3659
3660 /// vector to already allocated array (len>=p_MaxComp(v,r))
3661 // also used for p_Vec2Polys
p_Vec2Array(poly v,poly * p,int len,const ring r)3662 void p_Vec2Array(poly v, poly *p, int len, const ring r)
3663 {
3664 poly h;
3665 int k;
3666
3667 for(int i=len-1;i>=0;i--) p[i]=NULL;
3668 while (v!=NULL)
3669 {
3670 h=p_Head(v,r);
3671 k=__p_GetComp(h,r);
3672 if (k>len) { Werror("wrong rank:%d, should be %d",len,k); }
3673 else
3674 {
3675 p_SetComp(h,0,r);
3676 p_Setm(h,r);
3677 pNext(h)=p[k-1];p[k-1]=h;
3678 }
3679 pIter(v);
3680 }
3681 for(int i=len-1;i>=0;i--)
3682 {
3683 if (p[i]!=NULL) p[i]=pReverse(p[i]);
3684 }
3685 }
3686
3687 /*2
3688 * convert a vector to a set of polys,
3689 * allocates the polyset, (entries 0..(*len)-1)
3690 * the vector will not be changed
3691 */
p_Vec2Polys(poly v,poly ** p,int * len,const ring r)3692 void p_Vec2Polys(poly v, poly* *p, int *len, const ring r)
3693 {
3694 poly h;
3695 int k;
3696
3697 *len=p_MaxComp(v,r);
3698 if (*len==0) *len=1;
3699 *p=(poly*)omAlloc((*len)*sizeof(poly));
3700 p_Vec2Array(v,*p,*len,r);
3701 }
3702
3703 //
3704 // resets the pFDeg and pLDeg: if pLDeg is not given, it is
3705 // set to currRing->pLDegOrig, i.e. to the respective LDegProc which
3706 // only uses pFDeg (and not p_Deg, or pTotalDegree, etc)
pSetDegProcs(ring r,pFDegProc new_FDeg,pLDegProc new_lDeg)3707 void pSetDegProcs(ring r, pFDegProc new_FDeg, pLDegProc new_lDeg)
3708 {
3709 assume(new_FDeg != NULL);
3710 r->pFDeg = new_FDeg;
3711
3712 if (new_lDeg == NULL)
3713 new_lDeg = r->pLDegOrig;
3714
3715 r->pLDeg = new_lDeg;
3716 }
3717
3718 // restores pFDeg and pLDeg:
pRestoreDegProcs(ring r,pFDegProc old_FDeg,pLDegProc old_lDeg)3719 void pRestoreDegProcs(ring r, pFDegProc old_FDeg, pLDegProc old_lDeg)
3720 {
3721 assume(old_FDeg != NULL && old_lDeg != NULL);
3722 r->pFDeg = old_FDeg;
3723 r->pLDeg = old_lDeg;
3724 }
3725
3726 /*-------- several access procedures to monomials -------------------- */
3727 /*
3728 * the module weights for std
3729 */
3730 STATIC_VAR pFDegProc pOldFDeg;
3731 STATIC_VAR pLDegProc pOldLDeg;
3732 STATIC_VAR BOOLEAN pOldLexOrder;
3733
pModDeg(poly p,ring r)3734 static long pModDeg(poly p, ring r)
3735 {
3736 long d=pOldFDeg(p, r);
3737 int c=__p_GetComp(p, r);
3738 if ((c>0) && ((r->pModW)->range(c-1))) d+= (*(r->pModW))[c-1];
3739 return d;
3740 //return pOldFDeg(p, r)+(*pModW)[p_GetComp(p, r)-1];
3741 }
3742
p_SetModDeg(intvec * w,ring r)3743 void p_SetModDeg(intvec *w, ring r)
3744 {
3745 if (w!=NULL)
3746 {
3747 r->pModW = w;
3748 pOldFDeg = r->pFDeg;
3749 pOldLDeg = r->pLDeg;
3750 pOldLexOrder = r->pLexOrder;
3751 pSetDegProcs(r,pModDeg);
3752 r->pLexOrder = TRUE;
3753 }
3754 else
3755 {
3756 r->pModW = NULL;
3757 pRestoreDegProcs(r,pOldFDeg, pOldLDeg);
3758 r->pLexOrder = pOldLexOrder;
3759 }
3760 }
3761
3762 /*2
3763 * handle memory request for sets of polynomials (ideals)
3764 * l is the length of *p, increment is the difference (may be negative)
3765 */
pEnlargeSet(poly ** p,int l,int increment)3766 void pEnlargeSet(poly* *p, int l, int increment)
3767 {
3768 poly* h;
3769
3770 if (*p==NULL)
3771 {
3772 if (increment==0) return;
3773 h=(poly*)omAlloc0(increment*sizeof(poly));
3774 }
3775 else
3776 {
3777 h=(poly*)omReallocSize((poly*)*p,l*sizeof(poly),(l+increment)*sizeof(poly));
3778 if (increment>0)
3779 {
3780 memset(&(h[l]),0,increment*sizeof(poly));
3781 }
3782 }
3783 *p=h;
3784 }
3785
3786 /*2
3787 *divides p1 by its leading coefficient
3788 */
p_Norm(poly p1,const ring r)3789 void p_Norm(poly p1, const ring r)
3790 {
3791 if (rField_is_Ring(r))
3792 {
3793 if (!n_IsUnit(pGetCoeff(p1), r->cf)) return;
3794 // Werror("p_Norm not possible in the case of coefficient rings.");
3795 }
3796 else if (p1!=NULL)
3797 {
3798 if (pNext(p1)==NULL)
3799 {
3800 p_SetCoeff(p1,n_Init(1,r->cf),r);
3801 return;
3802 }
3803 poly h;
3804 if (!n_IsOne(pGetCoeff(p1),r->cf))
3805 {
3806 number k, c;
3807 n_Normalize(pGetCoeff(p1),r->cf);
3808 k = pGetCoeff(p1);
3809 c = n_Init(1,r->cf);
3810 pSetCoeff0(p1,c);
3811 h = pNext(p1);
3812 while (h!=NULL)
3813 {
3814 c=n_Div(pGetCoeff(h),k,r->cf);
3815 // no need to normalize: Z/p, R
3816 // normalize already in nDiv: Q_a, Z/p_a
3817 // remains: Q
3818 if (rField_is_Q(r) && (!n_IsOne(c,r->cf))) n_Normalize(c,r->cf);
3819 p_SetCoeff(h,c,r);
3820 pIter(h);
3821 }
3822 n_Delete(&k,r->cf);
3823 }
3824 else
3825 {
3826 //if (r->cf->cfNormalize != nDummy2) //TODO: OPTIMIZE
3827 {
3828 h = pNext(p1);
3829 while (h!=NULL)
3830 {
3831 n_Normalize(pGetCoeff(h),r->cf);
3832 pIter(h);
3833 }
3834 }
3835 }
3836 }
3837 }
3838
3839 /*2
3840 *normalize all coefficients
3841 */
p_Normalize(poly p,const ring r)3842 void p_Normalize(poly p,const ring r)
3843 {
3844 if ((rField_has_simple_inverse(r)) /* Z/p, GF(p,n), R, long R/C */
3845 || (r->cf->cfNormalize==ndNormalize)) /* Nemo rings, ...*/
3846 return;
3847 while (p!=NULL)
3848 {
3849 // no test befor n_Normalize: n_Normalize should fix problems
3850 n_Normalize(pGetCoeff(p),r->cf);
3851 pIter(p);
3852 }
3853 }
3854
3855 // splits p into polys with Exp(n) == 0 and Exp(n) != 0
3856 // Poly with Exp(n) != 0 is reversed
p_SplitAndReversePoly(poly p,int n,poly * non_zero,poly * zero,const ring r)3857 static void p_SplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero, const ring r)
3858 {
3859 if (p == NULL)
3860 {
3861 *non_zero = NULL;
3862 *zero = NULL;
3863 return;
3864 }
3865 spolyrec sz;
3866 poly z, n_z, next;
3867 z = &sz;
3868 n_z = NULL;
3869
3870 while(p != NULL)
3871 {
3872 next = pNext(p);
3873 if (p_GetExp(p, n,r) == 0)
3874 {
3875 pNext(z) = p;
3876 pIter(z);
3877 }
3878 else
3879 {
3880 pNext(p) = n_z;
3881 n_z = p;
3882 }
3883 p = next;
3884 }
3885 pNext(z) = NULL;
3886 *zero = pNext(&sz);
3887 *non_zero = n_z;
3888 }
3889 /*3
3890 * substitute the n-th variable by 1 in p
3891 * destroy p
3892 */
p_Subst1(poly p,int n,const ring r)3893 static poly p_Subst1 (poly p,int n, const ring r)
3894 {
3895 poly qq=NULL, result = NULL;
3896 poly zero=NULL, non_zero=NULL;
3897
3898 // reverse, so that add is likely to be linear
3899 p_SplitAndReversePoly(p, n, &non_zero, &zero,r);
3900
3901 while (non_zero != NULL)
3902 {
3903 assume(p_GetExp(non_zero, n,r) != 0);
3904 qq = non_zero;
3905 pIter(non_zero);
3906 qq->next = NULL;
3907 p_SetExp(qq,n,0,r);
3908 p_Setm(qq,r);
3909 result = p_Add_q(result,qq,r);
3910 }
3911 p = p_Add_q(result, zero,r);
3912 p_Test(p,r);
3913 return p;
3914 }
3915
3916 /*3
3917 * substitute the n-th variable by number e in p
3918 * destroy p
3919 */
p_Subst2(poly p,int n,number e,const ring r)3920 static poly p_Subst2 (poly p,int n, number e, const ring r)
3921 {
3922 assume( ! n_IsZero(e,r->cf) );
3923 poly qq,result = NULL;
3924 number nn, nm;
3925 poly zero, non_zero;
3926
3927 // reverse, so that add is likely to be linear
3928 p_SplitAndReversePoly(p, n, &non_zero, &zero,r);
3929
3930 while (non_zero != NULL)
3931 {
3932 assume(p_GetExp(non_zero, n, r) != 0);
3933 qq = non_zero;
3934 pIter(non_zero);
3935 qq->next = NULL;
3936 n_Power(e, p_GetExp(qq, n, r), &nn,r->cf);
3937 nm = n_Mult(nn, pGetCoeff(qq),r->cf);
3938 #ifdef HAVE_RINGS
3939 if (n_IsZero(nm,r->cf))
3940 {
3941 p_LmFree(&qq,r);
3942 n_Delete(&nm,r->cf);
3943 }
3944 else
3945 #endif
3946 {
3947 p_SetCoeff(qq, nm,r);
3948 p_SetExp(qq, n, 0,r);
3949 p_Setm(qq,r);
3950 result = p_Add_q(result,qq,r);
3951 }
3952 n_Delete(&nn,r->cf);
3953 }
3954 p = p_Add_q(result, zero,r);
3955 p_Test(p,r);
3956 return p;
3957 }
3958
3959
3960 /* delete monoms whose n-th exponent is different from zero */
p_Subst0(poly p,int n,const ring r)3961 static poly p_Subst0(poly p, int n, const ring r)
3962 {
3963 spolyrec res;
3964 poly h = &res;
3965 pNext(h) = p;
3966
3967 while (pNext(h)!=NULL)
3968 {
3969 if (p_GetExp(pNext(h),n,r)!=0)
3970 {
3971 p_LmDelete(&pNext(h),r);
3972 }
3973 else
3974 {
3975 pIter(h);
3976 }
3977 }
3978 p_Test(pNext(&res),r);
3979 return pNext(&res);
3980 }
3981
3982 /*2
3983 * substitute the n-th variable by e in p
3984 * destroy p
3985 */
p_Subst(poly p,int n,poly e,const ring r)3986 poly p_Subst(poly p, int n, poly e, const ring r)
3987 {
3988 #ifdef HAVE_SHIFTBBA
3989 // also don't even use p_Subst0 for Letterplace
3990 if (rIsLPRing(r))
3991 {
3992 poly subst = p_LPSubst(p, n, e, r);
3993 p_Delete(&p, r);
3994 return subst;
3995 }
3996 #endif
3997
3998 if (e == NULL) return p_Subst0(p, n,r);
3999
4000 if (p_IsConstant(e,r))
4001 {
4002 if (n_IsOne(pGetCoeff(e),r->cf)) return p_Subst1(p,n,r);
4003 else return p_Subst2(p, n, pGetCoeff(e),r);
4004 }
4005
4006 #ifdef HAVE_PLURAL
4007 if (rIsPluralRing(r))
4008 {
4009 return nc_pSubst(p,n,e,r);
4010 }
4011 #endif
4012
4013 int exponent,i;
4014 poly h, res, m;
4015 int *me,*ee;
4016 number nu,nu1;
4017
4018 me=(int *)omAlloc((rVar(r)+1)*sizeof(int));
4019 ee=(int *)omAlloc((rVar(r)+1)*sizeof(int));
4020 if (e!=NULL) p_GetExpV(e,ee,r);
4021 res=NULL;
4022 h=p;
4023 while (h!=NULL)
4024 {
4025 if ((e!=NULL) || (p_GetExp(h,n,r)==0))
4026 {
4027 m=p_Head(h,r);
4028 p_GetExpV(m,me,r);
4029 exponent=me[n];
4030 me[n]=0;
4031 for(i=rVar(r);i>0;i--)
4032 me[i]+=exponent*ee[i];
4033 p_SetExpV(m,me,r);
4034 if (e!=NULL)
4035 {
4036 n_Power(pGetCoeff(e),exponent,&nu,r->cf);
4037 nu1=n_Mult(pGetCoeff(m),nu,r->cf);
4038 n_Delete(&nu,r->cf);
4039 p_SetCoeff(m,nu1,r);
4040 }
4041 res=p_Add_q(res,m,r);
4042 }
4043 p_LmDelete(&h,r);
4044 }
4045 omFreeSize((ADDRESS)me,(rVar(r)+1)*sizeof(int));
4046 omFreeSize((ADDRESS)ee,(rVar(r)+1)*sizeof(int));
4047 return res;
4048 }
4049
4050 /*2
4051 * returns a re-ordered convertion of a number as a polynomial,
4052 * with permutation of parameters
4053 * NOTE: this only works for Frank's alg. & trans. fields
4054 */
n_PermNumber(const number z,const int * par_perm,const int,const ring src,const ring dst)4055 poly n_PermNumber(const number z, const int *par_perm, const int , const ring src, const ring dst)
4056 {
4057 #if 0
4058 PrintS("\nSource Ring: \n");
4059 rWrite(src);
4060
4061 if(0)
4062 {
4063 number zz = n_Copy(z, src->cf);
4064 PrintS("z: "); n_Write(zz, src);
4065 n_Delete(&zz, src->cf);
4066 }
4067
4068 PrintS("\nDestination Ring: \n");
4069 rWrite(dst);
4070
4071 /*Print("\nOldPar: %d\n", OldPar);
4072 for( int i = 1; i <= OldPar; i++ )
4073 {
4074 Print("par(%d) -> par/var (%d)\n", i, par_perm[i-1]);
4075 }*/
4076 #endif
4077 if( z == NULL )
4078 return NULL;
4079
4080 const coeffs srcCf = src->cf;
4081 assume( srcCf != NULL );
4082
4083 assume( !nCoeff_is_GF(srcCf) );
4084 assume( src->cf->extRing!=NULL );
4085
4086 poly zz = NULL;
4087
4088 const ring srcExtRing = srcCf->extRing;
4089 assume( srcExtRing != NULL );
4090
4091 const coeffs dstCf = dst->cf;
4092 assume( dstCf != NULL );
4093
4094 if( nCoeff_is_algExt(srcCf) ) // nCoeff_is_GF(srcCf)?
4095 {
4096 zz = (poly) z;
4097 if( zz == NULL ) return NULL;
4098 }
4099 else if (nCoeff_is_transExt(srcCf))
4100 {
4101 assume( !IS0(z) );
4102
4103 zz = NUM((fraction)z);
4104 p_Test (zz, srcExtRing);
4105
4106 if( zz == NULL ) return NULL;
4107 if( !DENIS1((fraction)z) )
4108 {
4109 if (!p_IsConstant(DEN((fraction)z),srcExtRing))
4110 WarnS("Not defined: Cannot map a rational fraction and make a polynomial out of it! Ignoring the denominator.");
4111 }
4112 }
4113 else
4114 {
4115 assume (FALSE);
4116 WerrorS("Number permutation is not implemented for this data yet!");
4117 return NULL;
4118 }
4119
4120 assume( zz != NULL );
4121 p_Test (zz, srcExtRing);
4122
4123 nMapFunc nMap = n_SetMap(srcExtRing->cf, dstCf);
4124
4125 assume( nMap != NULL );
4126
4127 poly qq;
4128 if ((par_perm == NULL) && (rPar(dst) != 0 && rVar (srcExtRing) > 0))
4129 {
4130 int* perm;
4131 perm=(int *)omAlloc0((rVar(srcExtRing)+1)*sizeof(int));
4132 for(int i=si_min(rVar(srcExtRing),rPar(dst));i>0;i--)
4133 perm[i]=-i;
4134 qq = p_PermPoly(zz, perm, srcExtRing, dst, nMap, NULL, rVar(srcExtRing)-1);
4135 omFreeSize ((ADDRESS)perm, (rVar(srcExtRing)+1)*sizeof(int));
4136 }
4137 else
4138 qq = p_PermPoly(zz, par_perm-1, srcExtRing, dst, nMap, NULL, rVar (srcExtRing)-1);
4139
4140 if(nCoeff_is_transExt(srcCf)
4141 && (!DENIS1((fraction)z))
4142 && p_IsConstant(DEN((fraction)z),srcExtRing))
4143 {
4144 number n=nMap(pGetCoeff(DEN((fraction)z)),srcExtRing->cf, dstCf);
4145 qq=p_Div_nn(qq,n,dst);
4146 n_Delete(&n,dstCf);
4147 p_Normalize(qq,dst);
4148 }
4149 p_Test (qq, dst);
4150
4151 return qq;
4152 }
4153
4154
4155 /*2
4156 *returns a re-ordered copy of a polynomial, with permutation of the variables
4157 */
p_PermPoly(poly p,const int * perm,const ring oldRing,const ring dst,nMapFunc nMap,const int * par_perm,int OldPar,BOOLEAN use_mult)4158 poly p_PermPoly (poly p, const int * perm, const ring oldRing, const ring dst,
4159 nMapFunc nMap, const int *par_perm, int OldPar, BOOLEAN use_mult)
4160 {
4161 #if 0
4162 p_Test(p, oldRing);
4163 PrintS("p_PermPoly::p: "); p_Write(p, oldRing, oldRing);
4164 #endif
4165 const int OldpVariables = rVar(oldRing);
4166 poly result = NULL;
4167 poly result_last = NULL;
4168 poly aq = NULL; /* the map coefficient */
4169 poly qq; /* the mapped monomial */
4170 assume(dst != NULL);
4171 assume(dst->cf != NULL);
4172 #ifdef HAVE_PLURAL
4173 poly tmp_mm=p_One(dst);
4174 #endif
4175 while (p != NULL)
4176 {
4177 // map the coefficient
4178 if ( ((OldPar == 0) || (par_perm == NULL) || rField_is_GF(oldRing) || (nMap==ndCopyMap))
4179 && (nMap != NULL) )
4180 {
4181 qq = p_Init(dst);
4182 assume( nMap != NULL );
4183 number n = nMap(p_GetCoeff(p, oldRing), oldRing->cf, dst->cf);
4184 n_Test (n,dst->cf);
4185 if ( nCoeff_is_algExt(dst->cf) )
4186 n_Normalize(n, dst->cf);
4187 p_GetCoeff(qq, dst) = n;// Note: n can be a ZERO!!!
4188 }
4189 else
4190 {
4191 qq = p_One(dst);
4192 // aq = naPermNumber(p_GetCoeff(p, oldRing), par_perm, OldPar, oldRing); // no dst???
4193 // poly n_PermNumber(const number z, const int *par_perm, const int P, const ring src, const ring dst)
4194 aq = n_PermNumber(p_GetCoeff(p, oldRing), par_perm, OldPar, oldRing, dst);
4195 p_Test(aq, dst);
4196 if ( nCoeff_is_algExt(dst->cf) )
4197 p_Normalize(aq,dst);
4198 if (aq == NULL)
4199 p_SetCoeff(qq, n_Init(0, dst->cf),dst); // Very dirty trick!!!
4200 p_Test(aq, dst);
4201 }
4202 if (rRing_has_Comp(dst))
4203 p_SetComp(qq, p_GetComp(p, oldRing), dst);
4204 if ( n_IsZero(pGetCoeff(qq), dst->cf) )
4205 {
4206 p_LmDelete(&qq,dst);
4207 qq = NULL;
4208 }
4209 else
4210 {
4211 // map pars:
4212 int mapped_to_par = 0;
4213 for(int i = 1; i <= OldpVariables; i++)
4214 {
4215 int e = p_GetExp(p, i, oldRing);
4216 if (e != 0)
4217 {
4218 if (perm==NULL)
4219 p_SetExp(qq, i, e, dst);
4220 else if (perm[i]>0)
4221 {
4222 #ifdef HAVE_PLURAL
4223 if(use_mult)
4224 {
4225 p_SetExp(tmp_mm,perm[i],e,dst);
4226 p_Setm(tmp_mm,dst);
4227 qq=p_Mult_mm(qq,tmp_mm,dst);
4228 p_SetExp(tmp_mm,perm[i],0,dst);
4229
4230 }
4231 else
4232 #endif
4233 p_AddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/, dst);
4234 }
4235 else if (perm[i]<0)
4236 {
4237 number c = p_GetCoeff(qq, dst);
4238 if (rField_is_GF(dst))
4239 {
4240 assume( dst->cf->extRing == NULL );
4241 number ee = n_Param(1, dst);
4242 number eee;
4243 n_Power(ee, e, &eee, dst->cf); //nfDelete(ee,dst);
4244 ee = n_Mult(c, eee, dst->cf);
4245 //nfDelete(c,dst);nfDelete(eee,dst);
4246 pSetCoeff0(qq,ee);
4247 }
4248 else if (nCoeff_is_Extension(dst->cf))
4249 {
4250 const int par = -perm[i];
4251 assume( par > 0 );
4252 // WarnS("longalg missing 3");
4253 #if 1
4254 const coeffs C = dst->cf;
4255 assume( C != NULL );
4256 const ring R = C->extRing;
4257 assume( R != NULL );
4258 assume( par <= rVar(R) );
4259 poly pcn; // = (number)c
4260 assume( !n_IsZero(c, C) );
4261 if( nCoeff_is_algExt(C) )
4262 pcn = (poly) c;
4263 else // nCoeff_is_transExt(C)
4264 pcn = NUM((fraction)c);
4265 if (pNext(pcn) == NULL) // c->z
4266 p_AddExp(pcn, -perm[i], e, R);
4267 else /* more difficult: we have really to multiply: */
4268 {
4269 poly mmc = p_ISet(1, R);
4270 p_SetExp(mmc, -perm[i], e, R);
4271 p_Setm(mmc, R);
4272 number nnc;
4273 // convert back to a number: number nnc = mmc;
4274 if( nCoeff_is_algExt(C) )
4275 nnc = (number) mmc;
4276 else // nCoeff_is_transExt(C)
4277 nnc = ntInit(mmc, C);
4278 p_GetCoeff(qq, dst) = n_Mult((number)c, nnc, C);
4279 n_Delete((number *)&c, C);
4280 n_Delete((number *)&nnc, C);
4281 }
4282 mapped_to_par=1;
4283 #endif
4284 }
4285 }
4286 else
4287 {
4288 /* this variable maps to 0 !*/
4289 p_LmDelete(&qq, dst);
4290 break;
4291 }
4292 }
4293 }
4294 if ( mapped_to_par && (qq!= NULL) && nCoeff_is_algExt(dst->cf) )
4295 {
4296 number n = p_GetCoeff(qq, dst);
4297 n_Normalize(n, dst->cf);
4298 p_GetCoeff(qq, dst) = n;
4299 }
4300 }
4301 pIter(p);
4302
4303 #if 0
4304 p_Test(aq,dst);
4305 PrintS("aq: "); p_Write(aq, dst, dst);
4306 #endif
4307
4308
4309 #if 1
4310 if (qq!=NULL)
4311 {
4312 p_Setm(qq,dst);
4313
4314 p_Test(aq,dst);
4315 p_Test(qq,dst);
4316
4317 #if 0
4318 PrintS("qq: "); p_Write(qq, dst, dst);
4319 #endif
4320
4321 if (aq!=NULL)
4322 qq=p_Mult_q(aq,qq,dst);
4323 aq = qq;
4324 while (pNext(aq) != NULL) pIter(aq);
4325 if (result_last==NULL)
4326 {
4327 result=qq;
4328 }
4329 else
4330 {
4331 pNext(result_last)=qq;
4332 }
4333 result_last=aq;
4334 aq = NULL;
4335 }
4336 else if (aq!=NULL)
4337 {
4338 p_Delete(&aq,dst);
4339 }
4340 }
4341 result=p_SortAdd(result,dst);
4342 #else
4343 // if (qq!=NULL)
4344 // {
4345 // pSetm(qq);
4346 // pTest(qq);
4347 // pTest(aq);
4348 // if (aq!=NULL) qq=pMult(aq,qq);
4349 // aq = qq;
4350 // while (pNext(aq) != NULL) pIter(aq);
4351 // pNext(aq) = result;
4352 // aq = NULL;
4353 // result = qq;
4354 // }
4355 // else if (aq!=NULL)
4356 // {
4357 // pDelete(&aq);
4358 // }
4359 //}
4360 //p = result;
4361 //result = NULL;
4362 //while (p != NULL)
4363 //{
4364 // qq = p;
4365 // pIter(p);
4366 // qq->next = NULL;
4367 // result = pAdd(result, qq);
4368 //}
4369 #endif
4370 p_Test(result,dst);
4371 #if 0
4372 p_Test(result,dst);
4373 PrintS("result: "); p_Write(result,dst,dst);
4374 #endif
4375 #ifdef HAVE_PLURAL
4376 p_LmDelete(&tmp_mm,dst);
4377 #endif
4378 return result;
4379 }
4380 /**************************************************************
4381 *
4382 * Jet
4383 *
4384 **************************************************************/
4385
pp_Jet(poly p,int m,const ring R)4386 poly pp_Jet(poly p, int m, const ring R)
4387 {
4388 poly r=NULL;
4389 poly t=NULL;
4390
4391 while (p!=NULL)
4392 {
4393 if (p_Totaldegree(p,R)<=m)
4394 {
4395 if (r==NULL)
4396 r=p_Head(p,R);
4397 else
4398 if (t==NULL)
4399 {
4400 pNext(r)=p_Head(p,R);
4401 t=pNext(r);
4402 }
4403 else
4404 {
4405 pNext(t)=p_Head(p,R);
4406 pIter(t);
4407 }
4408 }
4409 pIter(p);
4410 }
4411 return r;
4412 }
4413
p_Jet(poly p,int m,const ring R)4414 poly p_Jet(poly p, int m,const ring R)
4415 {
4416 while((p!=NULL) && (p_Totaldegree(p,R)>m)) p_LmDelete(&p,R);
4417 if (p==NULL) return NULL;
4418 poly r=p;
4419 while (pNext(p)!=NULL)
4420 {
4421 if (p_Totaldegree(pNext(p),R)>m)
4422 {
4423 p_LmDelete(&pNext(p),R);
4424 }
4425 else
4426 pIter(p);
4427 }
4428 return r;
4429 }
4430
pp_JetW(poly p,int m,int * w,const ring R)4431 poly pp_JetW(poly p, int m, int *w, const ring R)
4432 {
4433 poly r=NULL;
4434 poly t=NULL;
4435 while (p!=NULL)
4436 {
4437 if (totaldegreeWecart_IV(p,R,w)<=m)
4438 {
4439 if (r==NULL)
4440 r=p_Head(p,R);
4441 else
4442 if (t==NULL)
4443 {
4444 pNext(r)=p_Head(p,R);
4445 t=pNext(r);
4446 }
4447 else
4448 {
4449 pNext(t)=p_Head(p,R);
4450 pIter(t);
4451 }
4452 }
4453 pIter(p);
4454 }
4455 return r;
4456 }
4457
p_JetW(poly p,int m,int * w,const ring R)4458 poly p_JetW(poly p, int m, int *w, const ring R)
4459 {
4460 while((p!=NULL) && (totaldegreeWecart_IV(p,R,w)>m)) p_LmDelete(&p,R);
4461 if (p==NULL) return NULL;
4462 poly r=p;
4463 while (pNext(p)!=NULL)
4464 {
4465 if (totaldegreeWecart_IV(pNext(p),R,w)>m)
4466 {
4467 p_LmDelete(&pNext(p),R);
4468 }
4469 else
4470 pIter(p);
4471 }
4472 return r;
4473 }
4474
4475 /*************************************************************/
p_MinDeg(poly p,intvec * w,const ring R)4476 int p_MinDeg(poly p,intvec *w, const ring R)
4477 {
4478 if(p==NULL)
4479 return -1;
4480 int d=-1;
4481 while(p!=NULL)
4482 {
4483 int d0=0;
4484 for(int j=0;j<rVar(R);j++)
4485 if(w==NULL||j>=w->length())
4486 d0+=p_GetExp(p,j+1,R);
4487 else
4488 d0+=(*w)[j]*p_GetExp(p,j+1,R);
4489 if(d0<d||d==-1)
4490 d=d0;
4491 pIter(p);
4492 }
4493 return d;
4494 }
4495
4496 /***************************************************************/
p_Invers(int n,poly u,intvec * w,const ring R)4497 static poly p_Invers(int n,poly u,intvec *w, const ring R)
4498 {
4499 if(n<0)
4500 return NULL;
4501 number u0=n_Invers(pGetCoeff(u),R->cf);
4502 poly v=p_NSet(u0,R);
4503 if(n==0)
4504 return v;
4505 int *ww=iv2array(w,R);
4506 poly u1=p_JetW(p_Sub(p_One(R),__p_Mult_nn(u,u0,R),R),n,ww,R);
4507 if(u1==NULL)
4508 {
4509 omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(int));
4510 return v;
4511 }
4512 poly v1=__p_Mult_nn(p_Copy(u1,R),u0,R);
4513 v=p_Add_q(v,p_Copy(v1,R),R);
4514 for(int i=n/p_MinDeg(u1,w,R);i>1;i--)
4515 {
4516 v1=p_JetW(p_Mult_q(v1,p_Copy(u1,R),R),n,ww,R);
4517 v=p_Add_q(v,p_Copy(v1,R),R);
4518 }
4519 p_Delete(&u1,R);
4520 p_Delete(&v1,R);
4521 omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(int));
4522 return v;
4523 }
4524
4525
p_Series(int n,poly p,poly u,intvec * w,const ring R)4526 poly p_Series(int n,poly p,poly u, intvec *w, const ring R)
4527 {
4528 int *ww=iv2array(w,R);
4529 if(p!=NULL)
4530 {
4531 if(u==NULL)
4532 p=p_JetW(p,n,ww,R);
4533 else
4534 p=p_JetW(p_Mult_q(p,p_Invers(n-p_MinDeg(p,w,R),u,w,R),R),n,ww,R);
4535 }
4536 omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(int));
4537 return p;
4538 }
4539
p_EqualPolys(poly p1,poly p2,const ring r)4540 BOOLEAN p_EqualPolys(poly p1,poly p2, const ring r)
4541 {
4542 while ((p1 != NULL) && (p2 != NULL))
4543 {
4544 if (! p_LmEqual(p1, p2,r))
4545 return FALSE;
4546 if (! n_Equal(p_GetCoeff(p1,r), p_GetCoeff(p2,r),r->cf ))
4547 return FALSE;
4548 pIter(p1);
4549 pIter(p2);
4550 }
4551 return (p1==p2);
4552 }
4553
p_ExpVectorEqual(poly p1,poly p2,const ring r1,const ring r2)4554 static inline BOOLEAN p_ExpVectorEqual(poly p1, poly p2, const ring r1, const ring r2)
4555 {
4556 assume( r1 == r2 || rSamePolyRep(r1, r2) );
4557
4558 p_LmCheckPolyRing1(p1, r1);
4559 p_LmCheckPolyRing1(p2, r2);
4560
4561 int i = r1->ExpL_Size;
4562
4563 assume( r1->ExpL_Size == r2->ExpL_Size );
4564
4565 unsigned long *ep = p1->exp;
4566 unsigned long *eq = p2->exp;
4567
4568 do
4569 {
4570 i--;
4571 if (ep[i] != eq[i]) return FALSE;
4572 }
4573 while (i);
4574
4575 return TRUE;
4576 }
4577
p_EqualPolys(poly p1,poly p2,const ring r1,const ring r2)4578 BOOLEAN p_EqualPolys(poly p1,poly p2, const ring r1, const ring r2)
4579 {
4580 assume( r1 == r2 || rSamePolyRep(r1, r2) ); // will be used in rEqual!
4581 assume( r1->cf == r2->cf );
4582
4583 while ((p1 != NULL) && (p2 != NULL))
4584 {
4585 // returns 1 if ExpVector(p)==ExpVector(q): does not compare numbers !!
4586 // #define p_LmEqual(p1, p2, r) p_ExpVectorEqual(p1, p2, r)
4587
4588 if (! p_ExpVectorEqual(p1, p2, r1, r2))
4589 return FALSE;
4590
4591 if (! n_Equal(p_GetCoeff(p1,r1), p_GetCoeff(p2,r2), r1->cf ))
4592 return FALSE;
4593
4594 pIter(p1);
4595 pIter(p2);
4596 }
4597 return (p1==p2);
4598 }
4599
4600 /*2
4601 *returns TRUE if p1 is a skalar multiple of p2
4602 *assume p1 != NULL and p2 != NULL
4603 */
p_ComparePolys(poly p1,poly p2,const ring r)4604 BOOLEAN p_ComparePolys(poly p1,poly p2, const ring r)
4605 {
4606 number n,nn;
4607 pAssume(p1 != NULL && p2 != NULL);
4608
4609 if (!p_LmEqual(p1,p2,r)) //compare leading mons
4610 return FALSE;
4611 if ((pNext(p1)==NULL) && (pNext(p2)!=NULL))
4612 return FALSE;
4613 if ((pNext(p2)==NULL) && (pNext(p1)!=NULL))
4614 return FALSE;
4615 if (pLength(p1) != pLength(p2))
4616 return FALSE;
4617 #ifdef HAVE_RINGS
4618 if (rField_is_Ring(r))
4619 {
4620 if (!n_DivBy(pGetCoeff(p1), pGetCoeff(p2), r->cf)) return FALSE;
4621 }
4622 #endif
4623 n=n_Div(pGetCoeff(p1),pGetCoeff(p2),r->cf);
4624 while ((p1 != NULL) /*&& (p2 != NULL)*/)
4625 {
4626 if ( ! p_LmEqual(p1, p2,r))
4627 {
4628 n_Delete(&n, r->cf);
4629 return FALSE;
4630 }
4631 if (!n_Equal(pGetCoeff(p1), nn = n_Mult(pGetCoeff(p2),n, r->cf), r->cf))
4632 {
4633 n_Delete(&n, r->cf);
4634 n_Delete(&nn, r->cf);
4635 return FALSE;
4636 }
4637 n_Delete(&nn, r->cf);
4638 pIter(p1);
4639 pIter(p2);
4640 }
4641 n_Delete(&n, r->cf);
4642 return TRUE;
4643 }
4644
4645 /*2
4646 * returns the length of a (numbers of monomials)
4647 * respect syzComp
4648 */
p_Last(const poly p,int & l,const ring r)4649 poly p_Last(const poly p, int &l, const ring r)
4650 {
4651 if (p == NULL)
4652 {
4653 l = 0;
4654 return NULL;
4655 }
4656 l = 1;
4657 poly a = p;
4658 if (! rIsSyzIndexRing(r))
4659 {
4660 poly next = pNext(a);
4661 while (next!=NULL)
4662 {
4663 a = next;
4664 next = pNext(a);
4665 l++;
4666 }
4667 }
4668 else
4669 {
4670 int curr_limit = rGetCurrSyzLimit(r);
4671 poly pp = a;
4672 while ((a=pNext(a))!=NULL)
4673 {
4674 if (__p_GetComp(a,r)<=curr_limit/*syzComp*/)
4675 l++;
4676 else break;
4677 pp = a;
4678 }
4679 a=pp;
4680 }
4681 return a;
4682 }
4683
p_Var(poly m,const ring r)4684 int p_Var(poly m,const ring r)
4685 {
4686 if (m==NULL) return 0;
4687 if (pNext(m)!=NULL) return 0;
4688 int i,e=0;
4689 for (i=rVar(r); i>0; i--)
4690 {
4691 int exp=p_GetExp(m,i,r);
4692 if (exp==1)
4693 {
4694 if (e==0) e=i;
4695 else return 0;
4696 }
4697 else if (exp!=0)
4698 {
4699 return 0;
4700 }
4701 }
4702 return e;
4703 }
4704
4705 /*2
4706 *the minimal index of used variables - 1
4707 */
p_LowVar(poly p,const ring r)4708 int p_LowVar (poly p, const ring r)
4709 {
4710 int k,l,lex;
4711
4712 if (p == NULL) return -1;
4713
4714 k = 32000;/*a very large dummy value*/
4715 while (p != NULL)
4716 {
4717 l = 1;
4718 lex = p_GetExp(p,l,r);
4719 while ((l < (rVar(r))) && (lex == 0))
4720 {
4721 l++;
4722 lex = p_GetExp(p,l,r);
4723 }
4724 l--;
4725 if (l < k) k = l;
4726 pIter(p);
4727 }
4728 return k;
4729 }
4730
4731 /*2
4732 * verschiebt die Indizees der Modulerzeugenden um i
4733 */
p_Shift(poly * p,int i,const ring r)4734 void p_Shift (poly * p,int i, const ring r)
4735 {
4736 poly qp1 = *p,qp2 = *p;/*working pointers*/
4737 int j = p_MaxComp(*p,r),k = p_MinComp(*p,r);
4738
4739 if (j+i < 0) return ;
4740 BOOLEAN toPoly= ((j == -i) && (j == k));
4741 while (qp1 != NULL)
4742 {
4743 if (toPoly || (__p_GetComp(qp1,r)+i > 0))
4744 {
4745 p_AddComp(qp1,i,r);
4746 p_SetmComp(qp1,r);
4747 qp2 = qp1;
4748 pIter(qp1);
4749 }
4750 else
4751 {
4752 if (qp2 == *p)
4753 {
4754 pIter(*p);
4755 p_LmDelete(&qp2,r);
4756 qp2 = *p;
4757 qp1 = *p;
4758 }
4759 else
4760 {
4761 qp2->next = qp1->next;
4762 if (qp1!=NULL) p_LmDelete(&qp1,r);
4763 qp1 = qp2->next;
4764 }
4765 }
4766 }
4767 }
4768
4769 /***************************************************************
4770 *
4771 * Storage Managament Routines
4772 *
4773 ***************************************************************/
4774
4775
GetBitFields(const long e,const unsigned int s,const unsigned int n)4776 static inline unsigned long GetBitFields(const long e,
4777 const unsigned int s, const unsigned int n)
4778 {
4779 #define Sy_bit_L(x) (((unsigned long)1L)<<(x))
4780 unsigned int i = 0;
4781 unsigned long ev = 0L;
4782 assume(n > 0 && s < BIT_SIZEOF_LONG);
4783 do
4784 {
4785 assume(s+i < BIT_SIZEOF_LONG);
4786 if (e > (long) i) ev |= Sy_bit_L(s+i);
4787 else break;
4788 i++;
4789 }
4790 while (i < n);
4791 return ev;
4792 }
4793
4794 // Short Exponent Vectors are used for fast divisibility tests
4795 // ShortExpVectors "squeeze" an exponent vector into one word as follows:
4796 // Let n = BIT_SIZEOF_LONG / pVariables.
4797 // If n == 0 (i.e. pVariables > BIT_SIZE_OF_LONG), let m == the number
4798 // of non-zero exponents. If (m>BIT_SIZEOF_LONG), then sev = ~0, else
4799 // first m bits of sev are set to 1.
4800 // Otherwise (i.e. pVariables <= BIT_SIZE_OF_LONG)
4801 // represented by a bit-field of length n (resp. n+1 for some
4802 // exponents). If the value of an exponent is greater or equal to n, then
4803 // all of its respective n bits are set to 1. If the value of an exponent
4804 // is smaller than n, say m, then only the first m bits of the respective
4805 // n bits are set to 1, the others are set to 0.
4806 // This way, we have:
4807 // exp1 / exp2 ==> (ev1 & ~ev2) == 0, i.e.,
4808 // if (ev1 & ~ev2) then exp1 does not divide exp2
p_GetShortExpVector(const poly p,const ring r)4809 unsigned long p_GetShortExpVector(const poly p, const ring r)
4810 {
4811 assume(p != NULL);
4812 unsigned long ev = 0; // short exponent vector
4813 unsigned int n = BIT_SIZEOF_LONG / r->N; // number of bits per exp
4814 unsigned int m1; // highest bit which is filled with (n+1)
4815 int i=0,j=1;
4816
4817 if (n == 0)
4818 {
4819 if (r->N <2*BIT_SIZEOF_LONG)
4820 {
4821 n=1;
4822 m1=0;
4823 }
4824 else
4825 {
4826 for (; j<=r->N; j++)
4827 {
4828 if (p_GetExp(p,j,r) > 0) i++;
4829 if (i == BIT_SIZEOF_LONG) break;
4830 }
4831 if (i>0)
4832 ev = ~(0UL) >> (BIT_SIZEOF_LONG - i);
4833 return ev;
4834 }
4835 }
4836 else
4837 {
4838 m1 = (n+1)*(BIT_SIZEOF_LONG - n*r->N);
4839 }
4840
4841 n++;
4842 while (i<m1)
4843 {
4844 ev |= GetBitFields(p_GetExp(p, j,r), i, n);
4845 i += n;
4846 j++;
4847 }
4848
4849 n--;
4850 while (i<BIT_SIZEOF_LONG)
4851 {
4852 ev |= GetBitFields(p_GetExp(p, j,r), i, n);
4853 i += n;
4854 j++;
4855 }
4856 return ev;
4857 }
4858
4859
4860 /// p_GetShortExpVector of p * pp
p_GetShortExpVector(const poly p,const poly pp,const ring r)4861 unsigned long p_GetShortExpVector(const poly p, const poly pp, const ring r)
4862 {
4863 assume(p != NULL);
4864 assume(pp != NULL);
4865
4866 unsigned long ev = 0; // short exponent vector
4867 unsigned int n = BIT_SIZEOF_LONG / r->N; // number of bits per exp
4868 unsigned int m1; // highest bit which is filled with (n+1)
4869 int j=1;
4870 unsigned long i = 0L;
4871
4872 if (n == 0)
4873 {
4874 if (r->N <2*BIT_SIZEOF_LONG)
4875 {
4876 n=1;
4877 m1=0;
4878 }
4879 else
4880 {
4881 for (; j<=r->N; j++)
4882 {
4883 if (p_GetExp(p,j,r) > 0 || p_GetExp(pp,j,r) > 0) i++;
4884 if (i == BIT_SIZEOF_LONG) break;
4885 }
4886 if (i>0)
4887 ev = ~(0UL) >> (BIT_SIZEOF_LONG - i);
4888 return ev;
4889 }
4890 }
4891 else
4892 {
4893 m1 = (n+1)*(BIT_SIZEOF_LONG - n*r->N);
4894 }
4895
4896 n++;
4897 while (i<m1)
4898 {
4899 ev |= GetBitFields(p_GetExp(p, j,r) + p_GetExp(pp, j,r), i, n);
4900 i += n;
4901 j++;
4902 }
4903
4904 n--;
4905 while (i<BIT_SIZEOF_LONG)
4906 {
4907 ev |= GetBitFields(p_GetExp(p, j,r) + p_GetExp(pp, j,r), i, n);
4908 i += n;
4909 j++;
4910 }
4911 return ev;
4912 }
4913
4914
4915
4916 /***************************************************************
4917 *
4918 * p_ShallowDelete
4919 *
4920 ***************************************************************/
4921 #undef LINKAGE
4922 #define LINKAGE
4923 #undef p_Delete__T
4924 #define p_Delete__T p_ShallowDelete
4925 #undef n_Delete__T
4926 #define n_Delete__T(n, r) do {} while (0)
4927
4928 #include "polys/templates/p_Delete__T.cc"
4929
4930 /***************************************************************/
4931 /*
4932 * compare a and b
4933 */
p_Compare(const poly a,const poly b,const ring R)4934 int p_Compare(const poly a, const poly b, const ring R)
4935 {
4936 int r=p_Cmp(a,b,R);
4937 if ((r==0)&&(a!=NULL))
4938 {
4939 number h=n_Sub(pGetCoeff(a),pGetCoeff(b),R->cf);
4940 /* compare lead coeffs */
4941 r = -1+n_IsZero(h,R->cf)+2*n_GreaterZero(h,R->cf); /* -1: <, 0:==, 1: > */
4942 n_Delete(&h,R->cf);
4943 }
4944 else if (a==NULL)
4945 {
4946 if (b==NULL)
4947 {
4948 /* compare 0, 0 */
4949 r=0;
4950 }
4951 else if(p_IsConstant(b,R))
4952 {
4953 /* compare 0, const */
4954 r = 1-2*n_GreaterZero(pGetCoeff(b),R->cf); /* -1: <, 1: > */
4955 }
4956 }
4957 else if (b==NULL)
4958 {
4959 if (p_IsConstant(a,R))
4960 {
4961 /* compare const, 0 */
4962 r = -1+2*n_GreaterZero(pGetCoeff(a),R->cf); /* -1: <, 1: > */
4963 }
4964 }
4965 return(r);
4966 }
4967
p_GcdMon(poly f,poly g,const ring r)4968 poly p_GcdMon(poly f, poly g, const ring r)
4969 {
4970 assume(f!=NULL);
4971 assume(g!=NULL);
4972 assume(pNext(f)==NULL);
4973 poly G=p_Head(f,r);
4974 poly h=g;
4975 int *mf=(int*)omAlloc((r->N+1)*sizeof(int));
4976 p_GetExpV(f,mf,r);
4977 int *mh=(int*)omAlloc((r->N+1)*sizeof(int));
4978 BOOLEAN const_mon;
4979 BOOLEAN one_coeff=n_IsOne(pGetCoeff(G),r->cf);
4980 loop
4981 {
4982 if (h==NULL) break;
4983 if(!one_coeff)
4984 {
4985 number n=n_SubringGcd(pGetCoeff(G),pGetCoeff(h),r->cf);
4986 one_coeff=n_IsOne(n,r->cf);
4987 p_SetCoeff(G,n,r);
4988 }
4989 p_GetExpV(h,mh,r);
4990 const_mon=TRUE;
4991 for(unsigned j=r->N;j!=0;j--)
4992 {
4993 if (mh[j]<mf[j]) mf[j]=mh[j];
4994 if (mf[j]>0) const_mon=FALSE;
4995 }
4996 if (one_coeff && const_mon) break;
4997 pIter(h);
4998 }
4999 mf[0]=0;
5000 p_SetExpV(G,mf,r); // included is p_SetComp, p_Setm
5001 omFreeSize(mf,(r->N+1)*sizeof(int));
5002 omFreeSize(mh,(r->N+1)*sizeof(int));
5003 return G;
5004 }
5005
p_CopyPowerProduct(poly p,const ring r)5006 poly p_CopyPowerProduct(poly p, const ring r)
5007 {
5008 if (p == NULL) return NULL;
5009 p_LmCheckPolyRing1(p, r);
5010 poly np;
5011 omTypeAllocBin(poly, np, r->PolyBin);
5012 p_SetRingOfLm(np, r);
5013 memcpy(np->exp, p->exp, r->ExpL_Size*sizeof(long));
5014 pNext(np) = NULL;
5015 pSetCoeff0(np, n_Init(1, r->cf));
5016 return np;
5017 }
5018
p_MaxExpPerVar(poly p,int i,const ring r)5019 int p_MaxExpPerVar(poly p, int i, const ring r)
5020 {
5021 int m=0;
5022 while(p!=NULL)
5023 {
5024 int mm=p_GetExp(p,i,r);
5025 if (mm>m) m=mm;
5026 pIter(p);
5027 }
5028 return m;
5029 }
5030
5031