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