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