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