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