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