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