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