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