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 
23 #include <polys/PolyEnumerator.h>
24 
25 #define TRANSEXT_PRIVATES
26 
29 
30 #include <polys/weight.h>
31 #include <polys/simpleideals.h>
32 
33 #include "ring.h"
34 #include "p_polys.h"
35 
39 
40 
41 // #include <???/ideals.h>
42 // #include <???/int64vec.h>
43 
44 #ifndef SING_NDEBUG
45 // #include <???/febase.h>
46 #endif
47 
48 #ifdef HAVE_PLURAL
49 #include "nc/nc.h"
50 #include "nc/sca.h"
51 #endif
52 
53 #include "coeffrings.h"
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);
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)) 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) Print("***************************\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 */
1224 int p_IsPurePower(const poly p, const ring r)
1225 {
1226 #ifdef HAVE_RINGS
1227  if (rField_is_Ring(r))
1228  {
1229  if (p == NULL) return 0;
1230  if (!n_IsUnit(pGetCoeff(p), r->cf)) return 0;
1231  }
1232 #endif
1233  int i,k=0;
1234 
1235  for (i=r->N;i;i--)
1236  {
1237  if (p_GetExp(p,i, r)!=0)
1238  {
1239  if(k!=0) return 0;
1240  k=i;
1241  }
1242  }
1243  return k;
1244 }
1245 
1246 /*2
1247 *test if a polynomial is univariate
1248 * return -1 for constant,
1249 * 0 for not univariate,s
1250 * i if dep. on var(i)
1251 */
1252 int p_IsUnivariate(poly p, const ring r)
1253 {
1254  int i,k=-1;
1255 
1256  while (p!=NULL)
1257  {
1258  for (i=r->N;i;i--)
1259  {
1260  if (p_GetExp(p,i, r)!=0)
1261  {
1262  if((k!=-1)&&(k!=i)) return 0;
1263  k=i;
1264  }
1265  }
1266  pIter(p);
1267  }
1268  return k;
1269 }
1270 
1271 // set entry e[i] to 1 if var(i) occurs in p, ignore var(j) if e[j]>0
1272 int p_GetVariables(poly p, int * e, const ring r)
1273 {
1274  int i;
1275  int n=0;
1276  while(p!=NULL)
1277  {
1278  n=0;
1279  for(i=r->N; i>0; i--)
1280  {
1281  if(e[i]==0)
1282  {
1283  if (p_GetExp(p,i,r)>0)
1284  {
1285  e[i]=1;
1286  n++;
1287  }
1288  }
1289  else
1290  n++;
1291  }
1292  if (n==r->N) break;
1293  pIter(p);
1294  }
1295  return n;
1296 }
1297 
1298 
1299 /*2
1300 * returns a polynomial representing the integer i
1301 */
1302 poly p_ISet(long i, const ring r)
1303 {
1304  poly rc = NULL;
1305  if (i!=0)
1306  {
1307  rc = p_Init(r);
1308  pSetCoeff0(rc,n_Init(i,r->cf));
1309  if (n_IsZero(pGetCoeff(rc),r->cf))
1310  p_LmDelete(&rc,r);
1311  }
1312  return rc;
1313 }
1314 
1315 /*2
1316 * an optimized version of p_ISet for the special case 1
1317 */
1318 poly p_One(const ring r)
1319 {
1320  poly rc = p_Init(r);
1321  pSetCoeff0(rc,n_Init(1,r->cf));
1322  return rc;
1323 }
1324 
1326 {
1327  *h=pNext(p);
1328  pNext(p)=NULL;
1329 }
1330 
1331 /*2
1332 * pair has no common factor ? or is no polynomial
1333 */
1334 BOOLEAN p_HasNotCF(poly p1, poly p2, const ring r)
1335 {
1336 
1337  if (p_GetComp(p1,r) > 0 || p_GetComp(p2,r) > 0)
1338  return FALSE;
1339  int i = rVar(r);
1340  loop
1341  {
1342  if ((p_GetExp(p1, i, r) > 0) && (p_GetExp(p2, i, r) > 0))
1343  return FALSE;
1344  i--;
1345  if (i == 0)
1346  return TRUE;
1347  }
1348 }
1349 
1350 /*2
1351 * convert monomial given as string to poly, e.g. 1x3y5z
1352 */
1353 const char * p_Read(const char *st, poly &rc, const ring r)
1354 {
1355  if (r==NULL) { rc=NULL;return st;}
1356  int i,j;
1357  rc = p_Init(r);
1358  const char *s = n_Read(st,&(p_GetCoeff(rc, r)),r->cf);
1359  if (s==st)
1360  /* i.e. it does not start with a coeff: test if it is a ringvar*/
1361  {
1362  j = r_IsRingVar(s,r->names,r->N);
1363  if (j >= 0)
1364  {
1365  p_IncrExp(rc,1+j,r);
1366  while (*s!='\0') s++;
1367  goto done;
1368  }
1369  }
1370  while (*s!='\0')
1371  {
1372  char ss[2];
1373  ss[0] = *s++;
1374  ss[1] = '\0';
1375  j = r_IsRingVar(ss,r->names,r->N);
1376  if (j >= 0)
1377  {
1378  const char *s_save=s;
1379  s = eati(s,&i);
1380  if (((unsigned long)i) > r->bitmask)
1381  {
1382  // exponent to large: it is not a monomial
1383  p_LmDelete(&rc,r);
1384  return s_save;
1385  }
1386  p_AddExp(rc,1+j, (long)i, r);
1387  }
1388  else
1389  {
1390  // 1st char of is not a varname
1391  // We return the parsed polynomial nevertheless. This is needed when
1392  // we are parsing coefficients in a rational function field.
1393  s--;
1394  break;
1395  }
1396  }
1397 done:
1398  if (n_IsZero(pGetCoeff(rc),r->cf)) p_LmDelete(&rc,r);
1399  else
1400  {
1401 #ifdef HAVE_PLURAL
1402  // in super-commutative ring
1403  // squares of anti-commutative variables are zeroes!
1404  if(rIsSCA(r))
1405  {
1406  const unsigned int iFirstAltVar = scaFirstAltVar(r);
1407  const unsigned int iLastAltVar = scaLastAltVar(r);
1408 
1409  assume(rc != NULL);
1410 
1411  for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
1412  if( p_GetExp(rc, k, r) > 1 )
1413  {
1414  p_LmDelete(&rc, r);
1415  goto finish;
1416  }
1417  }
1418 #endif
1419 
1420  p_Setm(rc,r);
1421  }
1422 finish:
1423  return s;
1424 }
1425 poly p_mInit(const char *st, BOOLEAN &ok, const ring r)
1426 {
1427  poly p;
1428  const char *s=p_Read(st,p,r);
1429  if (*s!='\0')
1430  {
1431  if ((s!=st)&&isdigit(st[0]))
1432  {
1434  }
1435  ok=FALSE;
1436  p_Delete(&p,r);
1437  return NULL;
1438  }
1439  p_Test(p,r);
1440  ok=!errorreported;
1441  return p;
1442 }
1443 
1444 /*2
1445 * returns a polynomial representing the number n
1446 * destroys n
1447 */
1448 poly p_NSet(number n, const ring r)
1449 {
1450  if (n_IsZero(n,r->cf))
1451  {
1452  n_Delete(&n, r->cf);
1453  return NULL;
1454  }
1455  else
1456  {
1457  poly rc = p_Init(r);
1458  pSetCoeff0(rc,n);
1459  return rc;
1460  }
1461 }
1462 /*2
1463 * assumes that LM(a) = LM(b)*m, for some monomial m,
1464 * returns the multiplicant m,
1465 * leaves a and b unmodified
1466 */
1467 poly p_Divide(poly a, poly b, const ring r)
1468 {
1469  assume((p_GetComp(a,r)==p_GetComp(b,r)) || (p_GetComp(b,r)==0));
1470  int i;
1471  poly result = p_Init(r);
1472 
1473  for(i=(int)r->N; i; i--)
1474  p_SetExp(result,i, p_GetExp(a,i,r)- p_GetExp(b,i,r),r);
1475  p_SetComp(result, p_GetComp(a,r) - p_GetComp(b,r),r);
1476  p_Setm(result,r);
1477  return result;
1478 }
1479 
1480 poly p_Div_nn(poly p, const number n, const ring r)
1481 {
1482  pAssume(!n_IsZero(n,r->cf));
1483  p_Test(p, r);
1484 
1485  poly q = p;
1486  while (p != NULL)
1487  {
1488  number nc = pGetCoeff(p);
1489  pSetCoeff0(p, n_Div(nc, n, r->cf));
1490  n_Delete(&nc, r->cf);
1491  pIter(p);
1492  }
1493  p_Test(q, r);
1494  return q;
1495 }
1496 
1497 /*2
1498 * divides a by the monomial b, ignores monomials which are not divisible
1499 * assumes that b is not NULL, destroyes b
1500 */
1501 poly p_DivideM(poly a, poly b, const ring r)
1502 {
1503  if (a==NULL) { p_Delete(&b,r); return NULL; }
1504  poly result=a;
1505  poly prev=NULL;
1506  int i;
1507 #ifdef HAVE_RINGS
1508  number inv=pGetCoeff(b);
1509 #else
1510  number inv=n_Invers(pGetCoeff(b),r->cf);
1511 #endif
1512 
1513  while (a!=NULL)
1514  {
1515  if (p_DivisibleBy(b,a,r))
1516  {
1517  for(i=(int)r->N; i; i--)
1518  p_SubExp(a,i, p_GetExp(b,i,r),r);
1519  p_SubComp(a, p_GetComp(b,r),r);
1520  p_Setm(a,r);
1521  prev=a;
1522  pIter(a);
1523  }
1524  else
1525  {
1526  if (prev==NULL)
1527  {
1528  p_LmDelete(&result,r);
1529  a=result;
1530  }
1531  else
1532  {
1533  p_LmDelete(&pNext(prev),r);
1534  a=pNext(prev);
1535  }
1536  }
1537  }
1538 #ifdef HAVE_RINGS
1539  if (n_IsUnit(inv,r->cf))
1540  {
1541  inv = n_Invers(inv,r->cf);
1542  p_Mult_nn(result,inv,r);
1543  n_Delete(&inv, r->cf);
1544  }
1545  else
1546  {
1547  p_Div_nn(result,inv,r);
1548  }
1549 #else
1550  p_Mult_nn(result,inv,r);
1551  n_Delete(&inv, r->cf);
1552 #endif
1553  p_Delete(&b, r);
1554  return result;
1555 }
1556 
1557 #ifdef HAVE_RINGS
1558 /* TRUE iff LT(f) | LT(g) */
1560 {
1561  int exponent;
1562  for(int i = (int)rVar(r); i>0; i--)
1563  {
1564  exponent = p_GetExp(g, i, r) - p_GetExp(f, i, r);
1565  if (exponent < 0) return FALSE;
1566  }
1567  return n_DivBy(pGetCoeff(g), pGetCoeff(f), r->cf);
1568 }
1569 #endif
1570 
1571 // returns the LCM of the head terms of a and b in *m
1572 void p_Lcm(const poly a, const poly b, poly m, const ring r)
1573 {
1574  for (int i=rVar(r); i; --i)
1575  p_SetExp(m,i, si_max( p_GetExp(a,i,r), p_GetExp(b,i,r)),r);
1576 
1577  p_SetComp(m, si_max(p_GetComp(a,r), p_GetComp(b,r)),r);
1578  /* Don't do a pSetm here, otherwise hres/lres chockes */
1579 }
1580 
1581 
1582 
1583 #ifdef HAVE_RATGRING
1584 /*2
1585 * returns the rational LCM of the head terms of a and b
1586 * without coefficient!!!
1587 */
1588 poly p_LcmRat(const poly a, const poly b, const long lCompM, const ring r)
1589 {
1590  poly m = // p_One( r);
1591  p_Init(r);
1592 
1593 // const int (currRing->N) = r->N;
1594 
1595  // for (int i = (currRing->N); i>=r->real_var_start; i--)
1596  for (int i = r->real_var_end; i>=r->real_var_start; i--)
1597  {
1598  const int lExpA = p_GetExp (a, i, r);
1599  const int lExpB = p_GetExp (b, i, r);
1600 
1601  p_SetExp (m, i, si_max(lExpA, lExpB), r);
1602  }
1603 
1604  p_SetComp (m, lCompM, r);
1605  p_Setm(m,r);
1606  n_New(&(p_GetCoeff(m, r)), r);
1607 
1608  return(m);
1609 };
1610 
1611 void p_LmDeleteAndNextRat(poly *p, int ishift, ring r)
1612 {
1613  /* modifies p*/
1614  // Print("start: "); Print(" "); p_wrp(*p,r);
1615  p_LmCheckPolyRing2(*p, r);
1616  poly q = p_Head(*p,r);
1617  const long cmp = p_GetComp(*p, r);
1618  while ( ( (*p)!=NULL ) && ( p_Comp_k_n(*p, q, ishift+1, r) ) && (p_GetComp(*p, r) == cmp) )
1619  {
1620  p_LmDelete(p,r);
1621  // Print("while: ");p_wrp(*p,r);Print(" ");
1622  }
1623  // p_wrp(*p,r);Print(" ");
1624  // PrintS("end\n");
1625  p_LmDelete(&q,r);
1626 }
1627 
1628 
1629 /* returns x-coeff of p, i.e. a poly in x, s.t. corresponding xd-monomials
1630 have the same D-part and the component 0
1631 does not destroy p
1632 */
1633 poly p_GetCoeffRat(poly p, int ishift, ring r)
1634 {
1635  poly q = pNext(p);
1636  poly res; // = p_Head(p,r);
1637  res = p_GetExp_k_n(p, ishift+1, r->N, r); // does pSetm internally
1638  p_SetCoeff(res,n_Copy(p_GetCoeff(p,r),r),r);
1639  poly s;
1640  long cmp = p_GetComp(p, r);
1641  while ( (q!= NULL) && (p_Comp_k_n(p, q, ishift+1, r)) && (p_GetComp(q, r) == cmp) )
1642  {
1643  s = p_GetExp_k_n(q, ishift+1, r->N, r);
1644  p_SetCoeff(s,n_Copy(p_GetCoeff(q,r),r),r);
1645  res = p_Add_q(res,s,r);
1646  q = pNext(q);
1647  }
1648  cmp = 0;
1649  p_SetCompP(res,cmp,r);
1650  return res;
1651 }
1652 
1653 
1654 
1655 void p_ContentRat(poly &ph, const ring r)
1656 // changes ph
1657 // for rat coefficients in K(x1,..xN)
1658 {
1659  // init array of RatLeadCoeffs
1660  // poly p_GetCoeffRat(poly p, int ishift, ring r);
1661 
1662  int len=pLength(ph);
1663  poly *C = (poly *)omAlloc0((len+1)*sizeof(poly)); //rat coeffs
1664  poly *LM = (poly *)omAlloc0((len+1)*sizeof(poly)); // rat lead terms
1665  int *D = (int *)omAlloc0((len+1)*sizeof(int)); //degrees of coeffs
1666  int *L = (int *)omAlloc0((len+1)*sizeof(int)); //lengths of coeffs
1667  int k = 0;
1668  poly p = p_Copy(ph, r); // ph will be needed below
1669  int mintdeg = p_Totaldegree(p, r);
1670  int minlen = len;
1671  int dd = 0; int i;
1672  int HasConstantCoef = 0;
1673  int is = r->real_var_start - 1;
1674  while (p!=NULL)
1675  {
1676  LM[k] = p_GetExp_k_n(p,1,is, r); // need LmRat istead of p_HeadRat(p, is, currRing); !
1677  C[k] = p_GetCoeffRat(p, is, r);
1678  D[k] = p_Totaldegree(C[k], r);
1679  mintdeg = si_min(mintdeg,D[k]);
1680  L[k] = pLength(C[k]);
1681  minlen = si_min(minlen,L[k]);
1682  if (p_IsConstant(C[k], r))
1683  {
1684  // C[k] = const, so the content will be numerical
1685  HasConstantCoef = 1;
1686  // smth like goto cleanup and return(pContent(p));
1687  }
1688  p_LmDeleteAndNextRat(&p, is, r);
1689  k++;
1690  }
1691 
1692  // look for 1 element of minimal degree and of minimal length
1693  k--;
1694  poly d;
1695  int mindeglen = len;
1696  if (k<=0) // this poly is not a ratgring poly -> pContent
1697  {
1698  p_Delete(&C[0], r);
1699  p_Delete(&LM[0], r);
1700  p_Content(ph, r);
1701  goto cleanup;
1702  }
1703 
1704  int pmindeglen;
1705  for(i=0; i<=k; i++)
1706  {
1707  if (D[i] == mintdeg)
1708  {
1709  if (L[i] < mindeglen)
1710  {
1711  mindeglen=L[i];
1712  pmindeglen = i;
1713  }
1714  }
1715  }
1716  d = p_Copy(C[pmindeglen], r);
1717  // there are dd>=1 mindeg elements
1718  // and pmideglen is the coordinate of one of the smallest among them
1719 
1720  // poly g = singclap_gcd(p_Copy(p,r),p_Copy(q,r));
1721  // return naGcd(d,d2,currRing);
1722 
1723  // adjoin pContentRat here?
1724  for(i=0; i<=k; i++)
1725  {
1726  d=singclap_gcd(d,p_Copy(C[i], r), r);
1727  if (p_Totaldegree(d, r)==0)
1728  {
1729  // cleanup, pContent, return
1730  p_Delete(&d, r);
1731  for(;k>=0;k--)
1732  {
1733  p_Delete(&C[k], r);
1734  p_Delete(&LM[k], r);
1735  }
1736  p_Content(ph, r);
1737  goto cleanup;
1738  }
1739  }
1740  for(i=0; i<=k; i++)
1741  {
1742  poly h=singclap_pdivide(C[i],d, r);
1743  p_Delete(&C[i], r);
1744  C[i]=h;
1745  }
1746 
1747  // zusammensetzen,
1748  p=NULL; // just to be sure
1749  for(i=0; i<=k; i++)
1750  {
1751  p = p_Add_q(p, p_Mult_q(C[i],LM[i], r), r);
1752  C[i]=NULL; LM[i]=NULL;
1753  }
1754  p_Delete(&ph, r); // do not need it anymore
1755  ph = p;
1756  // aufraeumen, return
1757 cleanup:
1758  omFree(C);
1759  omFree(LM);
1760  omFree(D);
1761  omFree(L);
1762 }
1763 
1764 
1765 #endif
1766 
1767 
1768 /* assumes that p and divisor are univariate polynomials in r,
1769  mentioning the same variable;
1770  assumes divisor != NULL;
1771  p may be NULL;
1772  assumes a global monomial ordering in r;
1773  performs polynomial division of p by divisor:
1774  - afterwards p contains the remainder of the division, i.e.,
1775  p_before = result * divisor + p_afterwards;
1776  - if needResult == TRUE, then the method computes and returns 'result',
1777  otherwise NULL is returned (This parametrization can be used when
1778  one is only interested in the remainder of the division. In this
1779  case, the method will be slightly faster.)
1780  leaves divisor unmodified */
1781 poly p_PolyDiv(poly &p, const poly divisor, const BOOLEAN needResult, const ring r)
1782 {
1783  assume(divisor != NULL);
1784  if (p == NULL) return NULL;
1785 
1786  poly result = NULL;
1787  number divisorLC = p_GetCoeff(divisor, r);
1788  int divisorLE = p_GetExp(divisor, 1, r);
1789  while ((p != NULL) && (p_Deg(p, r) >= p_Deg(divisor, r)))
1790  {
1791  /* determine t = LT(p) / LT(divisor) */
1792  poly t = p_ISet(1, r);
1793  number c = n_Div(p_GetCoeff(p, r), divisorLC, r->cf);
1794  n_Normalize(c,r->cf);
1795  p_SetCoeff(t, c, r);
1796  int e = p_GetExp(p, 1, r) - divisorLE;
1797  p_SetExp(t, 1, e, r);
1798  p_Setm(t, r);
1799  if (needResult) result = p_Add_q(result, p_Copy(t, r), r);
1800  p = p_Add_q(p, p_Neg(p_Mult_q(t, p_Copy(divisor, r), r), r), r);
1801  }
1802  return result;
1803 }
1804 
1805 /*2
1806 * returns the partial differentiate of a by the k-th variable
1807 * does not destroy the input
1808 */
1809 poly p_Diff(poly a, int k, const ring r)
1810 {
1811  poly res, f, last;
1812  number t;
1813 
1814  last = res = NULL;
1815  while (a!=NULL)
1816  {
1817  if (p_GetExp(a,k,r)!=0)
1818  {
1819  f = p_LmInit(a,r);
1820  t = n_Init(p_GetExp(a,k,r),r->cf);
1821  pSetCoeff0(f,n_Mult(t,pGetCoeff(a),r->cf));
1822  n_Delete(&t,r->cf);
1823  if (n_IsZero(pGetCoeff(f),r->cf))
1824  p_LmDelete(&f,r);
1825  else
1826  {
1827  p_DecrExp(f,k,r);
1828  p_Setm(f,r);
1829  if (res==NULL)
1830  {
1831  res=last=f;
1832  }
1833  else
1834  {
1835  pNext(last)=f;
1836  last=f;
1837  }
1838  }
1839  }
1840  pIter(a);
1841  }
1842  return res;
1843 }
1844 
1845 static poly p_DiffOpM(poly a, poly b,BOOLEAN multiply, const ring r)
1846 {
1847  int i,j,s;
1848  number n,h,hh;
1849  poly p=p_One(r);
1850  n=n_Mult(pGetCoeff(a),pGetCoeff(b),r->cf);
1851  for(i=rVar(r);i>0;i--)
1852  {
1853  s=p_GetExp(b,i,r);
1854  if (s<p_GetExp(a,i,r))
1855  {
1856  n_Delete(&n,r->cf);
1857  p_LmDelete(&p,r);
1858  return NULL;
1859  }
1860  if (multiply)
1861  {
1862  for(j=p_GetExp(a,i,r); j>0;j--)
1863  {
1864  h = n_Init(s,r->cf);
1865  hh=n_Mult(n,h,r->cf);
1866  n_Delete(&h,r->cf);
1867  n_Delete(&n,r->cf);
1868  n=hh;
1869  s--;
1870  }
1871  p_SetExp(p,i,s,r);
1872  }
1873  else
1874  {
1875  p_SetExp(p,i,s-p_GetExp(a,i,r),r);
1876  }
1877  }
1878  p_Setm(p,r);
1879  /*if (multiply)*/ p_SetCoeff(p,n,r);
1880  if (n_IsZero(n,r->cf)) p=p_LmDeleteAndNext(p,r); // return NULL as p is a monomial
1881  return p;
1882 }
1883 
1884 poly p_DiffOp(poly a, poly b,BOOLEAN multiply, const ring r)
1885 {
1886  poly result=NULL;
1887  poly h;
1888  for(;a!=NULL;pIter(a))
1889  {
1890  for(h=b;h!=NULL;pIter(h))
1891  {
1892  result=p_Add_q(result,p_DiffOpM(a,h,multiply,r),r);
1893  }
1894  }
1895  return result;
1896 }
1897 /*2
1898 * subtract p2 from p1, p1 and p2 are destroyed
1899 * do not put attention on speed: the procedure is only used in the interpreter
1900 */
1901 poly p_Sub(poly p1, poly p2, const ring r)
1902 {
1903  return p_Add_q(p1, p_Neg(p2,r),r);
1904 }
1905 
1906 /*3
1907 * compute for a monomial m
1908 * the power m^exp, exp > 1
1909 * destroys p
1910 */
1911 static poly p_MonPower(poly p, int exp, const ring r)
1912 {
1913  int i;
1914 
1915  if(!n_IsOne(pGetCoeff(p),r->cf))
1916  {
1917  number x, y;
1918  y = pGetCoeff(p);
1919  n_Power(y,exp,&x,r->cf);
1920  n_Delete(&y,r->cf);
1921  pSetCoeff0(p,x);
1922  }
1923  for (i=rVar(r); i!=0; i--)
1924  {
1925  p_MultExp(p,i, exp,r);
1926  }
1927  p_Setm(p,r);
1928  return p;
1929 }
1930 
1931 /*3
1932 * compute for monomials p*q
1933 * destroys p, keeps q
1934 */
1935 static void p_MonMult(poly p, poly q, const ring r)
1936 {
1937  number x, y;
1938 
1939  y = pGetCoeff(p);
1940  x = n_Mult(y,pGetCoeff(q),r->cf);
1941  n_Delete(&y,r->cf);
1942  pSetCoeff0(p,x);
1943  //for (int i=pVariables; i!=0; i--)
1944  //{
1945  // pAddExp(p,i, pGetExp(q,i));
1946  //}
1947  //p->Order += q->Order;
1948  p_ExpVectorAdd(p,q,r);
1949 }
1950 
1951 /*3
1952 * compute for monomials p*q
1953 * keeps p, q
1954 */
1955 static poly p_MonMultC(poly p, poly q, const ring rr)
1956 {
1957  number x;
1958  poly r = p_Init(rr);
1959 
1960  x = n_Mult(pGetCoeff(p),pGetCoeff(q),rr->cf);
1961  pSetCoeff0(r,x);
1962  p_ExpVectorSum(r,p, q, rr);
1963  return r;
1964 }
1965 
1966 /*3
1967 * create binomial coef.
1968 */
1969 static number* pnBin(int exp, const ring r)
1970 {
1971  int e, i, h;
1972  number x, y, *bin=NULL;
1973 
1974  x = n_Init(exp,r->cf);
1975  if (n_IsZero(x,r->cf))
1976  {
1977  n_Delete(&x,r->cf);
1978  return bin;
1979  }
1980  h = (exp >> 1) + 1;
1981  bin = (number *)omAlloc0(h*sizeof(number));
1982  bin[1] = x;
1983  if (exp < 4)
1984  return bin;
1985  i = exp - 1;
1986  for (e=2; e<h; e++)
1987  {
1988  x = n_Init(i,r->cf);
1989  i--;
1990  y = n_Mult(x,bin[e-1],r->cf);
1991  n_Delete(&x,r->cf);
1992  x = n_Init(e,r->cf);
1993  bin[e] = n_ExactDiv(y,x,r->cf);
1994  n_Delete(&x,r->cf);
1995  n_Delete(&y,r->cf);
1996  }
1997  return bin;
1998 }
1999 
2000 static void pnFreeBin(number *bin, int exp,const coeffs r)
2001 {
2002  int e, h = (exp >> 1) + 1;
2003 
2004  if (bin[1] != NULL)
2005  {
2006  for (e=1; e<h; e++)
2007  n_Delete(&(bin[e]),r);
2008  }
2009  omFreeSize((ADDRESS)bin, h*sizeof(number));
2010 }
2011 
2012 /*
2013 * compute for a poly p = head+tail, tail is monomial
2014 * (head + tail)^exp, exp > 1
2015 * with binomial coef.
2016 */
2017 static poly p_TwoMonPower(poly p, int exp, const ring r)
2018 {
2019  int eh, e;
2020  long al;
2021  poly *a;
2022  poly tail, b, res, h;
2023  number x;
2024  number *bin = pnBin(exp,r);
2025 
2026  tail = pNext(p);
2027  if (bin == NULL)
2028  {
2029  p_MonPower(p,exp,r);
2030  p_MonPower(tail,exp,r);
2031  p_Test(p,r);
2032  return p;
2033  }
2034  eh = exp >> 1;
2035  al = (exp + 1) * sizeof(poly);
2036  a = (poly *)omAlloc(al);
2037  a[1] = p;
2038  for (e=1; e<exp; e++)
2039  {
2040  a[e+1] = p_MonMultC(a[e],p,r);
2041  }
2042  res = a[exp];
2043  b = p_Head(tail,r);
2044  for (e=exp-1; e>eh; e--)
2045  {
2046  h = a[e];
2047  x = n_Mult(bin[exp-e],pGetCoeff(h),r->cf);
2048  p_SetCoeff(h,x,r);
2049  p_MonMult(h,b,r);
2050  res = pNext(res) = h;
2051  p_MonMult(b,tail,r);
2052  }
2053  for (e=eh; e!=0; e--)
2054  {
2055  h = a[e];
2056  x = n_Mult(bin[e],pGetCoeff(h),r->cf);
2057  p_SetCoeff(h,x,r);
2058  p_MonMult(h,b,r);
2059  res = pNext(res) = h;
2060  p_MonMult(b,tail,r);
2061  }
2062  p_LmDelete(&tail,r);
2063  pNext(res) = b;
2064  pNext(b) = NULL;
2065  res = a[exp];
2066  omFreeSize((ADDRESS)a, al);
2067  pnFreeBin(bin, exp, r->cf);
2068 // tail=res;
2069 // while((tail!=NULL)&&(pNext(tail)!=NULL))
2070 // {
2071 // if(nIsZero(pGetCoeff(pNext(tail))))
2072 // {
2073 // pLmDelete(&pNext(tail));
2074 // }
2075 // else
2076 // pIter(tail);
2077 // }
2078  p_Test(res,r);
2079  return res;
2080 }
2081 
2082 static poly p_Pow(poly p, int i, const ring r)
2083 {
2084  poly rc = p_Copy(p,r);
2085  i -= 2;
2086  do
2087  {
2088  rc = p_Mult_q(rc,p_Copy(p,r),r);
2089  p_Normalize(rc,r);
2090  i--;
2091  }
2092  while (i != 0);
2093  return p_Mult_q(rc,p,r);
2094 }
2095 
2096 /*2
2097 * returns the i-th power of p
2098 * p will be destroyed
2099 */
2100 poly p_Power(poly p, int i, const ring r)
2101 {
2102  poly rc=NULL;
2103 
2104  if (i==0)
2105  {
2106  p_Delete(&p,r);
2107  return p_One(r);
2108  }
2109 
2110  if(p!=NULL)
2111  {
2112  if ( (i > 0) && ((unsigned long ) i > (r->bitmask)))
2113  {
2114  Werror("exponent %d is too large, max. is %ld",i,r->bitmask);
2115  return NULL;
2116  }
2117  switch (i)
2118  {
2119 // cannot happen, see above
2120 // case 0:
2121 // {
2122 // rc=pOne();
2123 // pDelete(&p);
2124 // break;
2125 // }
2126  case 1:
2127  rc=p;
2128  break;
2129  case 2:
2130  rc=p_Mult_q(p_Copy(p,r),p,r);
2131  break;
2132  default:
2133  if (i < 0)
2134  {
2135  p_Delete(&p,r);
2136  return NULL;
2137  }
2138  else
2139  {
2140 #ifdef HAVE_PLURAL
2141  if (rIsPluralRing(r)) /* in the NC case nothing helps :-( */
2142  {
2143  int j=i;
2144  rc = p_Copy(p,r);
2145  while (j>1)
2146  {
2147  rc = p_Mult_q(p_Copy(p,r),rc,r);
2148  j--;
2149  }
2150  p_Delete(&p,r);
2151  return rc;
2152  }
2153 #endif
2154  rc = pNext(p);
2155  if (rc == NULL)
2156  return p_MonPower(p,i,r);
2157  /* else: binom ?*/
2158  int char_p=rChar(r);
2159  if ((pNext(rc) != NULL)
2160 #ifdef HAVE_RINGS
2161  || rField_is_Ring(r)
2162 #endif
2163  )
2164  return p_Pow(p,i,r);
2165  if ((char_p==0) || (i<=char_p))
2166  return p_TwoMonPower(p,i,r);
2167  return p_Pow(p,i,r);
2168  }
2169  /*end default:*/
2170  }
2171  }
2172  return rc;
2173 }
2174 
2175 /* --------------------------------------------------------------------------------*/
2176 /* content suff */
2177 
2178 static number p_InitContent(poly ph, const ring r);
2179 
2180 #define CLEARENUMERATORS 1
2181 
2182 void p_Content(poly ph, const ring r)
2183 {
2184  assume( ph != NULL );
2185 
2186  assume( r != NULL ); assume( r->cf != NULL );
2187 
2188 
2189 #if CLEARENUMERATORS
2190  if( 0 )
2191  {
2192  const coeffs C = r->cf;
2193  // experimentall (recursive enumerator treatment) of alg. Ext!
2194  CPolyCoeffsEnumerator itr(ph);
2195  n_ClearContent(itr, r->cf);
2196 
2197  p_Test(ph, r); n_Test(pGetCoeff(ph), C);
2198  assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
2199 
2200  // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2201  return;
2202  }
2203 #endif
2204 
2205 
2206 #ifdef HAVE_RINGS
2207  if (rField_is_Ring(r))
2208  {
2209  if (rField_has_Units(r))
2210  {
2211  number k = n_GetUnit(pGetCoeff(ph),r->cf);
2212  if (!n_IsOne(k,r->cf))
2213  {
2214  number tmpGMP = k;
2215  k = n_Invers(k,r->cf);
2216  n_Delete(&tmpGMP,r->cf);
2217  poly h = pNext(ph);
2218  p_SetCoeff(ph, n_Mult(pGetCoeff(ph), k,r->cf),r);
2219  while (h != NULL)
2220  {
2221  p_SetCoeff(h, n_Mult(pGetCoeff(h), k,r->cf),r);
2222  pIter(h);
2223  }
2224 // assume( n_GreaterZero(pGetCoeff(ph),r->cf) );
2225 // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2226  }
2227  n_Delete(&k,r->cf);
2228  }
2229  return;
2230  }
2231 #endif
2232  number h,d;
2233  poly p;
2234 
2235  if(TEST_OPT_CONTENTSB) return;
2236  if(pNext(ph)==NULL)
2237  {
2238  p_SetCoeff(ph,n_Init(1,r->cf),r);
2239  }
2240  else
2241  {
2242  assume( pNext(ph) != NULL );
2243 #if CLEARENUMERATORS
2244  if( nCoeff_is_Q(r->cf) )
2245  {
2246  // experimentall (recursive enumerator treatment) of alg. Ext!
2247  CPolyCoeffsEnumerator itr(ph);
2248  n_ClearContent(itr, r->cf);
2249 
2250  p_Test(ph, r); n_Test(pGetCoeff(ph), r->cf);
2251  assume(n_GreaterZero(pGetCoeff(ph), r->cf)); // ??
2252 
2253  // if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2254  return;
2255  }
2256 #endif
2257 
2258  n_Normalize(pGetCoeff(ph),r->cf);
2259  if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2260  if (rField_is_Q(r)||(getCoeffType(r->cf)==n_transExt)) // should not be used anymore if CLEARENUMERATORS is 1
2261  {
2262  h=p_InitContent(ph,r);
2263  p=ph;
2264  }
2265  else
2266  {
2267  h=n_Copy(pGetCoeff(ph),r->cf);
2268  p = pNext(ph);
2269  }
2270  while (p!=NULL)
2271  {
2272  n_Normalize(pGetCoeff(p),r->cf);
2273  d=n_SubringGcd(h,pGetCoeff(p),r->cf);
2274  n_Delete(&h,r->cf);
2275  h = d;
2276  if(n_IsOne(h,r->cf))
2277  {
2278  break;
2279  }
2280  pIter(p);
2281  }
2282  p = ph;
2283  //number tmp;
2284  if(!n_IsOne(h,r->cf))
2285  {
2286  while (p!=NULL)
2287  {
2288  //d = nDiv(pGetCoeff(p),h);
2289  //tmp = nExactDiv(pGetCoeff(p),h);
2290  //if (!nEqual(d,tmp))
2291  //{
2292  // StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
2293  // nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
2294  // nWrite(tmp);Print(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
2295  //}
2296  //nDelete(&tmp);
2297  d = n_ExactDiv(pGetCoeff(p),h,r->cf);
2298  p_SetCoeff(p,d,r);
2299  pIter(p);
2300  }
2301  }
2302  n_Delete(&h,r->cf);
2303  if (rField_is_Q_a(r))
2304  {
2305  // special handling for alg. ext.:
2306  if (getCoeffType(r->cf)==n_algExt)
2307  {
2308  h = n_Init(1, r->cf->extRing->cf);
2309  p=ph;
2310  while (p!=NULL)
2311  { // each monom: coeff in Q_a
2312  poly c_n_n=(poly)pGetCoeff(p);
2313  poly c_n=c_n_n;
2314  while (c_n!=NULL)
2315  { // each monom: coeff in Q
2316  d=n_NormalizeHelper(h,pGetCoeff(c_n),r->cf->extRing->cf);
2317  n_Delete(&h,r->cf->extRing->cf);
2318  h=d;
2319  pIter(c_n);
2320  }
2321  pIter(p);
2322  }
2323  /* h contains the 1/lcm of all denominators in c_n_n*/
2324  //n_Normalize(h,r->cf->extRing->cf);
2325  if(!n_IsOne(h,r->cf->extRing->cf))
2326  {
2327  p=ph;
2328  while (p!=NULL)
2329  { // each monom: coeff in Q_a
2330  poly c_n=(poly)pGetCoeff(p);
2331  while (c_n!=NULL)
2332  { // each monom: coeff in Q
2333  d=n_Mult(h,pGetCoeff(c_n),r->cf->extRing->cf);
2334  n_Normalize(d,r->cf->extRing->cf);
2335  n_Delete(&pGetCoeff(c_n),r->cf->extRing->cf);
2336  pGetCoeff(c_n)=d;
2337  pIter(c_n);
2338  }
2339  pIter(p);
2340  }
2341  }
2342  n_Delete(&h,r->cf->extRing->cf);
2343  }
2344  /*else
2345  {
2346  // special handling for rat. functions.:
2347  number hzz =NULL;
2348  p=ph;
2349  while (p!=NULL)
2350  { // each monom: coeff in Q_a (Z_a)
2351  fraction f=(fraction)pGetCoeff(p);
2352  poly c_n=NUM(f);
2353  if (hzz==NULL)
2354  {
2355  hzz=n_Copy(pGetCoeff(c_n),r->cf->extRing->cf);
2356  pIter(c_n);
2357  }
2358  while ((c_n!=NULL)&&(!n_IsOne(hzz,r->cf->extRing->cf)))
2359  { // each monom: coeff in Q (Z)
2360  d=n_Gcd(hzz,pGetCoeff(c_n),r->cf->extRing->cf);
2361  n_Delete(&hzz,r->cf->extRing->cf);
2362  hzz=d;
2363  pIter(c_n);
2364  }
2365  pIter(p);
2366  }
2367  // hzz contains the gcd of all numerators in f
2368  h=n_Invers(hzz,r->cf->extRing->cf);
2369  n_Delete(&hzz,r->cf->extRing->cf);
2370  n_Normalize(h,r->cf->extRing->cf);
2371  if(!n_IsOne(h,r->cf->extRing->cf))
2372  {
2373  p=ph;
2374  while (p!=NULL)
2375  { // each monom: coeff in Q_a (Z_a)
2376  fraction f=(fraction)pGetCoeff(p);
2377  NUM(f)=p_Mult_nn(NUM(f),h,r->cf->extRing);
2378  p_Normalize(NUM(f),r->cf->extRing);
2379  pIter(p);
2380  }
2381  }
2382  n_Delete(&h,r->cf->extRing->cf);
2383  }*/
2384  }
2385  }
2386  if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2387 }
2388 
2389 // Not yet?
2390 #if 1 // currently only used by Singular/janet
2391 void p_SimpleContent(poly ph, int smax, const ring r)
2392 {
2393  if(TEST_OPT_CONTENTSB) return;
2394  if (ph==NULL) return;
2395  if (pNext(ph)==NULL)
2396  {
2397  p_SetCoeff(ph,n_Init(1,r->cf),r);
2398  return;
2399  }
2400  if ((pNext(pNext(ph))==NULL)||(!rField_is_Q(r)))
2401  {
2402  return;
2403  }
2404  number d=p_InitContent(ph,r);
2405  if (n_Size(d,r->cf)<=smax)
2406  {
2407  //if (TEST_OPT_PROT) PrintS("G");
2408  return;
2409  }
2410 
2411 
2412  poly p=ph;
2413  number h=d;
2414  if (smax==1) smax=2;
2415  while (p!=NULL)
2416  {
2417 #if 0
2418  d=n_Gcd(h,pGetCoeff(p),r->cf);
2419  n_Delete(&h,r->cf);
2420  h = d;
2421 #else
2422  STATISTIC(n_Gcd); nlInpGcd(h,pGetCoeff(p),r->cf); // FIXME? TODO? // extern void nlInpGcd(number &a, number b, const coeffs r);
2423 #endif
2424  if(n_Size(h,r->cf)<smax)
2425  {
2426  //if (TEST_OPT_PROT) PrintS("g");
2427  return;
2428  }
2429  pIter(p);
2430  }
2431  p = ph;
2432  if (!n_GreaterZero(pGetCoeff(p),r->cf)) h=n_InpNeg(h,r->cf);
2433  if(n_IsOne(h,r->cf)) return;
2434  //if (TEST_OPT_PROT) PrintS("c");
2435  while (p!=NULL)
2436  {
2437 #if 1
2438  d = n_ExactDiv(pGetCoeff(p),h,r->cf);
2439  p_SetCoeff(p,d,r);
2440 #else
2441  STATISTIC(n_ExactDiv); nlInpExactDiv(pGetCoeff(p),h,r->cf); // no such function... ?
2442 #endif
2443  pIter(p);
2444  }
2445  n_Delete(&h,r->cf);
2446 }
2447 #endif
2448 
2449 static number p_InitContent(poly ph, const ring r)
2450 // only for coefficients in Q and rational functions
2451 #if 0
2452 {
2454  assume(ph!=NULL);
2455  assume(pNext(ph)!=NULL);
2456  assume(rField_is_Q(r));
2457  if (pNext(pNext(ph))==NULL)
2458  {
2459  return n_GetNumerator(pGetCoeff(pNext(ph)),r->cf);
2460  }
2461  poly p=ph;
2462  number n1=n_GetNumerator(pGetCoeff(p),r->cf);
2463  pIter(p);
2464  number n2=n_GetNumerator(pGetCoeff(p),r->cf);
2465  pIter(p);
2466  number d;
2467  number t;
2468  loop
2469  {
2470  nlNormalize(pGetCoeff(p),r->cf);
2471  t=n_GetNumerator(pGetCoeff(p),r->cf);
2472  if (nlGreaterZero(t,r->cf))
2473  d=nlAdd(n1,t,r->cf);
2474  else
2475  d=nlSub(n1,t,r->cf);
2476  nlDelete(&t,r->cf);
2477  nlDelete(&n1,r->cf);
2478  n1=d;
2479  pIter(p);
2480  if (p==NULL) break;
2481  nlNormalize(pGetCoeff(p),r->cf);
2482  t=n_GetNumerator(pGetCoeff(p),r->cf);
2483  if (nlGreaterZero(t,r->cf))
2484  d=nlAdd(n2,t,r->cf);
2485  else
2486  d=nlSub(n2,t,r->cf);
2487  nlDelete(&t,r->cf);
2488  nlDelete(&n2,r->cf);
2489  n2=d;
2490  pIter(p);
2491  if (p==NULL) break;
2492  }
2493  d=nlGcd(n1,n2,r->cf);
2494  nlDelete(&n1,r->cf);
2495  nlDelete(&n2,r->cf);
2496  return d;
2497 }
2498 #else
2499 {
2500  number d=pGetCoeff(ph);
2501  int s;
2502  int s2=-1;
2503  if(rField_is_Q(r))
2504  {
2505  if (SR_HDL(d)&SR_INT) return d;
2506  s=mpz_size1(d->z);
2507  }
2508  else
2509  s=n_Size(d,r);
2510  number d2=d;
2511  loop
2512  {
2513  pIter(ph);
2514  if(ph==NULL)
2515  {
2516  if (s2==-1) return n_Copy(d,r->cf);
2517  break;
2518  }
2519  if (rField_is_Q(r))
2520  {
2521  if (SR_HDL(pGetCoeff(ph))&SR_INT)
2522  {
2523  s2=s;
2524  d2=d;
2525  s=0;
2526  d=pGetCoeff(ph);
2527  if (s2==0) break;
2528  }
2529  else if (mpz_size1((pGetCoeff(ph)->z))<=s)
2530  {
2531  s2=s;
2532  d2=d;
2533  d=pGetCoeff(ph);
2534  s=mpz_size1(d->z);
2535  }
2536  }
2537  else
2538  {
2539  int ns=n_Size(pGetCoeff(ph),r);
2540  if (ns<=3)
2541  {
2542  s2=s;
2543  d2=d;
2544  d=pGetCoeff(ph);
2545  s=ns;
2546  if (s2<=3) break;
2547  }
2548  else if (ns<s)
2549  {
2550  s2=s;
2551  d2=d;
2552  d=pGetCoeff(ph);
2553  s=ns;
2554  }
2555  }
2556  }
2557  return n_SubringGcd(d,d2,r->cf);
2558 }
2559 #endif
2560 
2561 //void pContent(poly ph)
2562 //{
2563 // number h,d;
2564 // poly p;
2565 //
2566 // p = ph;
2567 // if(pNext(p)==NULL)
2568 // {
2569 // pSetCoeff(p,nInit(1));
2570 // }
2571 // else
2572 // {
2573 //#ifdef PDEBUG
2574 // if (!pTest(p)) return;
2575 //#endif
2576 // nNormalize(pGetCoeff(p));
2577 // if(!nGreaterZero(pGetCoeff(ph)))
2578 // {
2579 // ph = pNeg(ph);
2580 // nNormalize(pGetCoeff(p));
2581 // }
2582 // h=pGetCoeff(p);
2583 // pIter(p);
2584 // while (p!=NULL)
2585 // {
2586 // nNormalize(pGetCoeff(p));
2587 // if (nGreater(h,pGetCoeff(p))) h=pGetCoeff(p);
2588 // pIter(p);
2589 // }
2590 // h=nCopy(h);
2591 // p=ph;
2592 // while (p!=NULL)
2593 // {
2594 // d=n_Gcd(h,pGetCoeff(p));
2595 // nDelete(&h);
2596 // h = d;
2597 // if(nIsOne(h))
2598 // {
2599 // break;
2600 // }
2601 // pIter(p);
2602 // }
2603 // p = ph;
2604 // //number tmp;
2605 // if(!nIsOne(h))
2606 // {
2607 // while (p!=NULL)
2608 // {
2609 // d = nExactDiv(pGetCoeff(p),h);
2610 // pSetCoeff(p,d);
2611 // pIter(p);
2612 // }
2613 // }
2614 // nDelete(&h);
2615 // if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
2616 // {
2617 // pTest(ph);
2618 // singclap_divide_content(ph);
2619 // pTest(ph);
2620 // }
2621 // }
2622 //}
2623 #if 0
2624 void p_Content(poly ph, const ring r)
2625 {
2626  number h,d;
2627  poly p;
2628 
2629  if(pNext(ph)==NULL)
2630  {
2631  p_SetCoeff(ph,n_Init(1,r->cf),r);
2632  }
2633  else
2634  {
2635  n_Normalize(pGetCoeff(ph),r->cf);
2636  if(!n_GreaterZero(pGetCoeff(ph),r->cf)) ph = p_Neg(ph,r);
2637  h=n_Copy(pGetCoeff(ph),r->cf);
2638  p = pNext(ph);
2639  while (p!=NULL)
2640  {
2641  n_Normalize(pGetCoeff(p),r->cf);
2642  d=n_Gcd(h,pGetCoeff(p),r->cf);
2643  n_Delete(&h,r->cf);
2644  h = d;
2645  if(n_IsOne(h,r->cf))
2646  {
2647  break;
2648  }
2649  pIter(p);
2650  }
2651  p = ph;
2652  //number tmp;
2653  if(!n_IsOne(h,r->cf))
2654  {
2655  while (p!=NULL)
2656  {
2657  //d = nDiv(pGetCoeff(p),h);
2658  //tmp = nExactDiv(pGetCoeff(p),h);
2659  //if (!nEqual(d,tmp))
2660  //{
2661  // StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
2662  // nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
2663  // nWrite(tmp);Print(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
2664  //}
2665  //nDelete(&tmp);
2666  d = n_ExactDiv(pGetCoeff(p),h,r->cf);
2667  p_SetCoeff(p,d,r->cf);
2668  pIter(p);
2669  }
2670  }
2671  n_Delete(&h,r->cf);
2672  //if ( (n_GetChar(r) == 1) || (n_GetChar(r) < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
2673  //{
2674  // singclap_divide_content(ph);
2675  // if(!n_GreaterZero(pGetCoeff(ph),r)) ph = p_Neg(ph,r);
2676  //}
2677  }
2678 }
2679 #endif
2680 /* ---------------------------------------------------------------------------*/
2681 /* cleardenom suff */
2682 poly p_Cleardenom(poly p, const ring r)
2683 {
2684  if( p == NULL )
2685  return NULL;
2686 
2687  assume( r != NULL ); assume( r->cf != NULL ); const coeffs C = r->cf;
2688 
2689 #if CLEARENUMERATORS
2690  if( 0 )
2691  {
2692  CPolyCoeffsEnumerator itr(p);
2693 
2694  n_ClearDenominators(itr, C);
2695 
2696  n_ClearContent(itr, C); // divide out the content
2697 
2698  p_Test(p, r); n_Test(pGetCoeff(p), C);
2699  assume(n_GreaterZero(pGetCoeff(p), C)); // ??
2700 // if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2701 
2702  return p;
2703  }
2704 #endif
2705 
2706 
2707  number d, h;
2708 
2709 #ifdef HAVE_RINGS
2710  if (rField_is_Ring(r))
2711  {
2712  p_Content(p,r);
2713  if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2714  return p;
2715  }
2716 #endif
2717 
2719  {
2720  if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2721  return p;
2722  }
2723 
2724  assume(p != NULL);
2725 
2726  if(pNext(p)==NULL)
2727  {
2728  /*
2729  if (TEST_OPT_CONTENTSB)
2730  {
2731  number n=n_GetDenom(pGetCoeff(p),r->cf);
2732  if (!n_IsOne(n,r->cf))
2733  {
2734  number nn=n_Mult(pGetCoeff(p),n,r->cf);
2735  n_Normalize(nn,r->cf);
2736  p_SetCoeff(p,nn,r);
2737  }
2738  n_Delete(&n,r->cf);
2739  }
2740  else
2741  */
2742  p_SetCoeff(p,n_Init(1,r->cf),r);
2743 
2744  /*assume( n_GreaterZero(pGetCoeff(p),C) );
2745  if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2746  */
2747  return p;
2748  }
2749 
2750  assume(pNext(p)!=NULL);
2751  poly start=p;
2752 
2753 #if 0 && CLEARENUMERATORS
2754 //CF: does not seem to work that well..
2755 
2756  if( nCoeff_is_Q(C) || nCoeff_is_Q_a(C) )
2757  {
2758  CPolyCoeffsEnumerator itr(p);
2759 
2760  n_ClearDenominators(itr, C);
2761 
2762  n_ClearContent(itr, C); // divide out the content
2763 
2764  p_Test(p, r); n_Test(pGetCoeff(p), C);
2765  assume(n_GreaterZero(pGetCoeff(p), C)); // ??
2766 // if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2767 
2768  return start;
2769  }
2770 #endif
2771 
2772  if(1)
2773  {
2774  h = n_Init(1,r->cf);
2775  while (p!=NULL)
2776  {
2777  n_Normalize(pGetCoeff(p),r->cf);
2778  d=n_NormalizeHelper(h,pGetCoeff(p),r->cf);
2779  n_Delete(&h,r->cf);
2780  h=d;
2781  pIter(p);
2782  }
2783  /* contains the 1/lcm of all denominators */
2784  if(!n_IsOne(h,r->cf))
2785  {
2786  p = start;
2787  while (p!=NULL)
2788  {
2789  /* should be: // NOTE: don't use ->coef!!!!
2790  * number hh;
2791  * nGetDenom(p->coef,&hh);
2792  * nMult(&h,&hh,&d);
2793  * nNormalize(d);
2794  * nDelete(&hh);
2795  * nMult(d,p->coef,&hh);
2796  * nDelete(&d);
2797  * nDelete(&(p->coef));
2798  * p->coef =hh;
2799  */
2800  d=n_Mult(h,pGetCoeff(p),r->cf);
2801  n_Normalize(d,r->cf);
2802  p_SetCoeff(p,d,r);
2803  pIter(p);
2804  }
2805  n_Delete(&h,r->cf);
2806  }
2807  n_Delete(&h,r->cf);
2808  p=start;
2809 
2810  p_Content(p,r);
2811 #ifdef HAVE_RATGRING
2812  if (rIsRatGRing(r))
2813  {
2814  /* quick unit detection in the rational case is done in gr_nc_bba */
2815  p_ContentRat(p, r);
2816  start=p;
2817  }
2818 #endif
2819  }
2820 
2821  if(!n_GreaterZero(pGetCoeff(p),C)) p = p_Neg(p,r);
2822 
2823  return start;
2824 }
2825 
2826 void p_Cleardenom_n(poly ph,const ring r,number &c)
2827 {
2828  const coeffs C = r->cf;
2829  number d, h;
2830 
2831  assume( ph != NULL );
2832 
2833  poly p = ph;
2834 
2835 #if CLEARENUMERATORS
2836  if( 0 )
2837  {
2838  CPolyCoeffsEnumerator itr(ph);
2839 
2840  n_ClearDenominators(itr, d, C); // multiply with common denom. d
2841  n_ClearContent(itr, h, C); // divide by the content h
2842 
2843  c = n_Div(d, h, C); // d/h
2844 
2845  n_Delete(&d, C);
2846  n_Delete(&h, C);
2847 
2848  n_Test(c, C);
2849 
2850  p_Test(ph, r); n_Test(pGetCoeff(ph), C);
2851  assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
2852 /*
2853  if(!n_GreaterZero(pGetCoeff(ph),C))
2854  {
2855  ph = p_Neg(ph,r);
2856  c = n_InpNeg(c, C);
2857  }
2858 */
2859  return;
2860  }
2861 #endif
2862 
2863 
2864  if( pNext(p) == NULL )
2865  {
2866  c=n_Invers(pGetCoeff(p), C);
2867  p_SetCoeff(p, n_Init(1, C), r);
2868 
2869  assume( n_GreaterZero(pGetCoeff(ph),C) );
2870  if(!n_GreaterZero(pGetCoeff(ph),C))
2871  {
2872  ph = p_Neg(ph,r);
2873  c = n_InpNeg(c, C);
2874  }
2875 
2876  return;
2877  }
2878 
2879  assume( pNext(p) != NULL );
2880 
2881 #if CLEARENUMERATORS
2882  if( nCoeff_is_Q(C) || nCoeff_is_Q_a(C) )
2883  {
2884  CPolyCoeffsEnumerator itr(ph);
2885 
2886  n_ClearDenominators(itr, d, C); // multiply with common denom. d
2887  n_ClearContent(itr, h, C); // divide by the content h
2888 
2889  c = n_Div(d, h, C); // d/h
2890 
2891  n_Delete(&d, C);
2892  n_Delete(&h, C);
2893 
2894  n_Test(c, C);
2895 
2896  p_Test(ph, r); n_Test(pGetCoeff(ph), C);
2897  assume(n_GreaterZero(pGetCoeff(ph), C)); // ??
2898 /*
2899  if(!n_GreaterZero(pGetCoeff(ph),C))
2900  {
2901  ph = p_Neg(ph,r);
2902  c = n_InpNeg(c, C);
2903  }
2904 */
2905  return;
2906  }
2907 #endif
2908 
2909 
2910 
2911 
2912  if(1)
2913  {
2914  h = n_Init(1,r->cf);
2915  while (p!=NULL)
2916  {
2917  n_Normalize(pGetCoeff(p),r->cf);
2918  d=n_NormalizeHelper(h,pGetCoeff(p),r->cf);
2919  n_Delete(&h,r->cf);
2920  h=d;
2921  pIter(p);
2922  }
2923  c=h;
2924  /* contains the 1/lcm of all denominators */
2925  if(!n_IsOne(h,r->cf))
2926  {
2927  p = ph;
2928  while (p!=NULL)
2929  {
2930  /* should be: // NOTE: don't use ->coef!!!!
2931  * number hh;
2932  * nGetDenom(p->coef,&hh);
2933  * nMult(&h,&hh,&d);
2934  * nNormalize(d);
2935  * nDelete(&hh);
2936  * nMult(d,p->coef,&hh);
2937  * nDelete(&d);
2938  * nDelete(&(p->coef));
2939  * p->coef =hh;
2940  */
2941  d=n_Mult(h,pGetCoeff(p),r->cf);
2942  n_Normalize(d,r->cf);
2943  p_SetCoeff(p,d,r);
2944  pIter(p);
2945  }
2946  if (rField_is_Q_a(r))
2947  {
2948  loop
2949  {
2950  h = n_Init(1,r->cf);
2951  p=ph;
2952  while (p!=NULL)
2953  {
2954  d=n_NormalizeHelper(h,pGetCoeff(p),r->cf);
2955  n_Delete(&h,r->cf);
2956  h=d;
2957  pIter(p);
2958  }
2959  /* contains the 1/lcm of all denominators */
2960  if(!n_IsOne(h,r->cf))
2961  {
2962  p = ph;
2963  while (p!=NULL)
2964  {
2965  /* should be: // NOTE: don't use ->coef!!!!
2966  * number hh;
2967  * nGetDenom(p->coef,&hh);
2968  * nMult(&h,&hh,&d);
2969  * nNormalize(d);
2970  * nDelete(&hh);
2971  * nMult(d,p->coef,&hh);
2972  * nDelete(&d);
2973  * nDelete(&(p->coef));
2974  * p->coef =hh;
2975  */
2976  d=n_Mult(h,pGetCoeff(p),r->cf);
2977  n_Normalize(d,r->cf);
2978  p_SetCoeff(p,d,r);
2979  pIter(p);
2980  }
2981  number t=n_Mult(c,h,r->cf);
2982  n_Delete(&c,r->cf);
2983  c=t;
2984  }
2985  else
2986  {
2987  break;
2988  }
2989  n_Delete(&h,r->cf);
2990  }
2991  }
2992  }
2993  }
2994 
2995  if(!n_GreaterZero(pGetCoeff(ph),C))
2996  {
2997  ph = p_Neg(ph,r);
2998  c = n_InpNeg(c, C);
2999  }
3000 
3001 }
3002 
3003  // normalization: for poly over Q: make poly primitive, integral
3004  // Qa make poly integral with leading
3005  // coefficient minimal in N
3006  // Q(t) make poly primitive, integral
3007 
3008 void p_ProjectiveUnique(poly ph, const ring r)
3009 {
3010  if( ph == NULL )
3011  return;
3012 
3013  assume( r != NULL ); assume( r->cf != NULL ); const coeffs C = r->cf;
3014 
3015  number h;
3016  poly p;
3017 
3018 #ifdef HAVE_RINGS
3019  if (rField_is_Ring(r))
3020  {
3021  p_Content(ph,r);
3022  if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
3023  assume( n_GreaterZero(pGetCoeff(ph),C) );
3024  return;
3025  }
3026 #endif
3027 
3029  {
3030  assume( n_GreaterZero(pGetCoeff(ph),C) );
3031  if(!n_GreaterZero(pGetCoeff(ph),C)) ph = p_Neg(ph,r);
3032  return;
3033  }
3034  p = ph;
3035 
3036  assume(p != NULL);
3037 
3038  if(pNext(p)==NULL) // a monomial
3039  {
3040  p_SetCoeff(p, n_Init(1, C), r);
3041  return;
3042  }
3043 
3044  assume(pNext(p)!=NULL);
3045 
3046  if(!rField_is_Q(r) && !nCoeff_is_transExt(C))
3047  {
3048  h = p_GetCoeff(p, C);
3049  number hInv = n_Invers(h, C);
3050  pIter(p);
3051  while (p!=NULL)
3052  {
3053  p_SetCoeff(p, n_Mult(p_GetCoeff(p, C), hInv, C), r);
3054  pIter(p);
3055  }
3056  n_Delete(&hInv, C);
3057  p = ph;
3058  p_SetCoeff(p, n_Init(1, C), r);
3059  }
3060 
3061  p_Cleardenom(ph, r); //performs also a p_Content
3062 
3063 
3064  /* normalize ph over a transcendental extension s.t.
3065  lead (ph) is > 0 if extRing->cf == Q
3066  or lead (ph) is monic if extRing->cf == Zp*/
3067  if (nCoeff_is_transExt(C))
3068  {
3069  p= ph;
3070  h= p_GetCoeff (p, C);
3071  fraction f = (fraction) h;
3072  number n=p_GetCoeff (NUM (f),C->extRing->cf);
3073  if (rField_is_Q (C->extRing))
3074  {
3075  if (!n_GreaterZero(n,C->extRing->cf))
3076  {
3077  p=p_Neg (p,r);
3078  }
3079  }
3080  else if (rField_is_Zp(C->extRing))
3081  {
3082  if (!n_IsOne (n, C->extRing->cf))
3083  {
3084  n=n_Invers (n,C->extRing->cf);
3085  nMapFunc nMap;
3086  nMap= n_SetMap (C->extRing->cf, C);
3087  number ninv= nMap (n,C->extRing->cf, C);
3088  p=p_Mult_nn (p, ninv, r);
3089  n_Delete (&ninv, C);
3090  n_Delete (&n, C->extRing->cf);
3091  }
3092  }
3093  p= ph;
3094  }
3095 
3096  return;
3097 }
3098 
3099 #if 0 /*unused*/
3100 number p_GetAllDenom(poly ph, const ring r)
3101 {
3102  number d=n_Init(1,r->cf);
3103  poly p = ph;
3104 
3105  while (p!=NULL)
3106  {
3107  number h=n_GetDenom(pGetCoeff(p),r->cf);
3108  if (!n_IsOne(h,r->cf))
3109  {
3110  number dd=n_Mult(d,h,r->cf);
3111  n_Delete(&d,r->cf);
3112  d=dd;
3113  }
3114  n_Delete(&h,r->cf);
3115  pIter(p);
3116  }
3117  return d;
3118 }
3119 #endif
3120 
3121 int p_Size(poly p, const ring r)
3122 {
3123  int count = 0;
3124  if (r->cf->has_simple_Alloc)
3125  return pLength(p);
3126  while ( p != NULL )
3127  {
3128  count+= n_Size( pGetCoeff( p ), r->cf );
3129  pIter( p );
3130  }
3131  return count;
3132 }
3133 
3134 /*2
3135 *make p homogeneous by multiplying the monomials by powers of x_varnum
3136 *assume: deg(var(varnum))==1
3137 */
3138 poly p_Homogen (poly p, int varnum, const ring r)
3139 {
3140  pFDegProc deg;
3141  if (r->pLexOrder && (r->order[0]==ringorder_lp))
3142  deg=p_Totaldegree;
3143  else
3144  deg=r->pFDeg;
3145 
3146  poly q=NULL, qn;
3147  int o,ii;
3148  sBucket_pt bp;
3149 
3150  if (p!=NULL)
3151  {
3152  if ((varnum < 1) || (varnum > rVar(r)))
3153  {
3154  return NULL;
3155  }
3156  o=deg(p,r);
3157  q=pNext(p);
3158  while (q != NULL)
3159  {
3160  ii=deg(q,r);
3161  if (ii>o) o=ii;
3162  pIter(q);
3163  }
3164  q = p_Copy(p,r);
3165  bp = sBucketCreate(r);
3166  while (q != NULL)
3167  {
3168  ii = o-deg(q,r);
3169  if (ii!=0)
3170  {
3171  p_AddExp(q,varnum, (long)ii,r);
3172  p_Setm(q,r);
3173  }
3174  qn = pNext(q);
3175  pNext(q) = NULL;
3176  sBucket_Add_p(bp, q, 1);
3177  q = qn;
3178  }
3179  sBucketDestroyAdd(bp, &q, &ii);
3180  }
3181  return q;
3182 }
3183 
3184 /*2
3185 *tests if p is homogeneous with respect to the actual weigths
3186 */
3188 {
3189  poly qp=p;
3190  int o;
3191 
3192  if ((p == NULL) || (pNext(p) == NULL)) return TRUE;
3193  pFDegProc d;
3194  if (r->pLexOrder && (r->order[0]==ringorder_lp))
3195  d=p_Totaldegree;
3196  else
3197  d=r->pFDeg;
3198  o = d(p,r);
3199  do
3200  {
3201  if (d(qp,r) != o) return FALSE;
3202  pIter(qp);
3203  }
3204  while (qp != NULL);
3205  return TRUE;
3206 }
3207 
3208 /*----------utilities for syzygies--------------*/
3209 BOOLEAN p_VectorHasUnitB(poly p, int * k, const ring r)
3210 {
3211  poly q=p,qq;
3212  int i;
3213 
3214  while (q!=NULL)
3215  {
3216  if (p_LmIsConstantComp(q,r))
3217  {
3218  i = p_GetComp(q,r);
3219  qq = p;
3220  while ((qq != q) && (p_GetComp(qq,r) != i)) pIter(qq);
3221  if (qq == q)
3222  {
3223  *k = i;
3224  return TRUE;
3225  }
3226  else
3227  pIter(q);
3228  }
3229  else pIter(q);
3230  }
3231  return FALSE;
3232 }
3233 
3234 void p_VectorHasUnit(poly p, int * k, int * len, const ring r)
3235 {
3236  poly q=p,qq;
3237  int i,j=0;
3238 
3239  *len = 0;
3240  while (q!=NULL)
3241  {
3242  if (p_LmIsConstantComp(q,r))
3243  {
3244  i = p_GetComp(q,r);
3245  qq = p;
3246  while ((qq != q) && (p_GetComp(qq,r) != i)) pIter(qq);
3247  if (qq == q)
3248  {
3249  j = 0;
3250  while (qq!=NULL)
3251  {
3252  if (p_GetComp(qq,r)==i) j++;
3253  pIter(qq);
3254  }
3255  if ((*len == 0) || (j<*len))
3256  {
3257  *len = j;
3258  *k = i;
3259  }
3260  }
3261  }
3262  pIter(q);
3263  }
3264 }
3265 
3266 poly p_TakeOutComp1(poly * p, int k, const ring r)
3267 {
3268  poly q = *p;
3269 
3270  if (q==NULL) return NULL;
3271 
3272  poly qq=NULL,result = NULL;
3273 
3274  if (p_GetComp(q,r)==k)
3275  {
3276  result = q; /* *p */
3277  while ((q!=NULL) && (p_GetComp(q,r)==k))
3278  {
3279  p_SetComp(q,0,r);
3280  p_SetmComp(q,r);
3281  qq = q;
3282  pIter(q);
3283  }
3284  *p = q;
3285  pNext(qq) = NULL;
3286  }
3287  if (q==NULL) return result;
3288 // if (pGetComp(q) > k) pGetComp(q)--;
3289  while (pNext(q)!=NULL)
3290  {
3291  if (p_GetComp(pNext(q),r)==k)
3292  {
3293  if (result==NULL)
3294  {
3295  result = pNext(q);
3296  qq = result;
3297  }
3298  else
3299  {
3300  pNext(qq) = pNext(q);
3301  pIter(qq);
3302  }
3303  pNext(q) = pNext(pNext(q));
3304  pNext(qq) =NULL;
3305  p_SetComp(qq,0,r);
3306  p_SetmComp(qq,r);
3307  }
3308  else
3309  {
3310  pIter(q);
3311 // if (pGetComp(q) > k) pGetComp(q)--;
3312  }
3313  }
3314  return result;
3315 }
3316 
3317 poly p_TakeOutComp(poly * p, int k, const ring r)
3318 {
3319  poly q = *p,qq=NULL,result = NULL;
3320 
3321  if (q==NULL) return NULL;
3322  BOOLEAN use_setmcomp=rOrd_SetCompRequiresSetm(r);
3323  if (p_GetComp(q,r)==k)
3324  {
3325  result = q;
3326  do
3327  {
3328  p_SetComp(q,0,r);
3329  if (use_setmcomp) p_SetmComp(q,r);
3330  qq = q;
3331  pIter(q);
3332  }
3333  while ((q!=NULL) && (p_GetComp(q,r)==k));
3334  *p = q;
3335  pNext(qq) = NULL;
3336  }
3337  if (q==NULL) return result;
3338  if (p_GetComp(q,r) > k)
3339  {
3340  p_SubComp(q,1,r);
3341  if (use_setmcomp) p_SetmComp(q,r);
3342  }
3343  poly pNext_q;
3344  while ((pNext_q=pNext(q))!=NULL)
3345  {
3346  if (p_GetComp(pNext_q,r)==k)
3347  {
3348  if (result==NULL)
3349  {
3350  result = pNext_q;
3351  qq = result;
3352  }
3353  else
3354  {
3355  pNext(qq) = pNext_q;
3356  pIter(qq);
3357  }
3358  pNext(q) = pNext(pNext_q);
3359  pNext(qq) =NULL;
3360  p_SetComp(qq,0,r);
3361  if (use_setmcomp) p_SetmComp(qq,r);
3362  }
3363  else
3364  {
3365  /*pIter(q);*/ q=pNext_q;
3366  if (p_GetComp(q,r) > k)
3367  {
3368  p_SubComp(q,1,r);
3369  if (use_setmcomp) p_SetmComp(q,r);
3370  }
3371  }
3372  }
3373  return result;
3374 }
3375 
3376 // Splits *p into two polys: *q which consists of all monoms with
3377 // component == comp and *p of all other monoms *lq == pLength(*q)
3378 void p_TakeOutComp(poly *r_p, long comp, poly *r_q, int *lq, const ring r)
3379 {
3380  spolyrec pp, qq;
3381  poly p, q, p_prev;
3382  int l = 0;
3383 
3384 #ifdef HAVE_ASSUME
3385  int lp = pLength(*r_p);
3386 #endif
3387 
3388  pNext(&pp) = *r_p;
3389  p = *r_p;
3390  p_prev = &pp;
3391  q = &qq;
3392 
3393  while(p != NULL)
3394  {
3395  while (p_GetComp(p,r) == comp)
3396  {
3397  pNext(q) = p;
3398  pIter(q);
3399  p_SetComp(p, 0,r);
3400  p_SetmComp(p,r);
3401  pIter(p);
3402  l++;
3403  if (p == NULL)
3404  {
3405  pNext(p_prev) = NULL;
3406  goto Finish;
3407  }
3408  }
3409  pNext(p_prev) = p;
3410  p_prev = p;
3411  pIter(p);
3412  }
3413 
3414  Finish:
3415  pNext(q) = NULL;
3416  *r_p = pNext(&pp);
3417  *r_q = pNext(&qq);
3418  *lq = l;
3419 #ifdef HAVE_ASSUME
3420  assume(pLength(*r_p) + pLength(*r_q) == lp);
3421 #endif
3422  p_Test(*r_p,r);
3423  p_Test(*r_q,r);
3424 }
3425 
3426 void p_DeleteComp(poly * p,int k, const ring r)
3427 {
3428  poly q;
3429 
3430  while ((*p!=NULL) && (p_GetComp(*p,r)==k)) p_LmDelete(p,r);
3431  if (*p==NULL) return;
3432  q = *p;
3433  if (p_GetComp(q,r)>k)
3434  {
3435  p_SubComp(q,1,r);
3436  p_SetmComp(q,r);
3437  }
3438  while (pNext(q)!=NULL)
3439  {
3440  if (p_GetComp(pNext(q),r)==k)
3441  p_LmDelete(&(pNext(q)),r);
3442  else
3443  {
3444  pIter(q);
3445  if (p_GetComp(q,r)>k)
3446  {
3447  p_SubComp(q,1,r);
3448  p_SetmComp(q,r);
3449  }
3450  }
3451  }
3452 }
3453 
3454 /*2
3455 * convert a vector to a set of polys,
3456 * allocates the polyset, (entries 0..(*len)-1)
3457 * the vector will not be changed
3458 */
3459 void p_Vec2Polys(poly v, poly* *p, int *len, const ring r)
3460 {
3461  poly h;
3462  int k;
3463 
3464  *len=p_MaxComp(v,r);
3465  if (*len==0) *len=1;
3466  *p=(poly*)omAlloc0((*len)*sizeof(poly));
3467  while (v!=NULL)
3468  {
3469  h=p_Head(v,r);
3470  k=p_GetComp(h,r);
3471  p_SetComp(h,0,r);
3472  (*p)[k-1]=p_Add_q((*p)[k-1],h,r);
3473  pIter(v);
3474  }
3475 }
3476 
3477 //
3478 // resets the pFDeg and pLDeg: if pLDeg is not given, it is
3479 // set to currRing->pLDegOrig, i.e. to the respective LDegProc which
3480 // only uses pFDeg (and not p_Deg, or pTotalDegree, etc)
3481 void pSetDegProcs(ring r, pFDegProc new_FDeg, pLDegProc new_lDeg)
3482 {
3483  assume(new_FDeg != NULL);
3484  r->pFDeg = new_FDeg;
3485 
3486  if (new_lDeg == NULL)
3487  new_lDeg = r->pLDegOrig;
3488 
3489  r->pLDeg = new_lDeg;
3490 }
3491 
3492 // restores pFDeg and pLDeg:
3493 void pRestoreDegProcs(ring r, pFDegProc old_FDeg, pLDegProc old_lDeg)
3494 {
3495  assume(old_FDeg != NULL && old_lDeg != NULL);
3496  r->pFDeg = old_FDeg;
3497  r->pLDeg = old_lDeg;
3498 }
3499 
3500 /*-------- several access procedures to monomials -------------------- */
3501 /*
3502 * the module weights for std
3503 */
3507 
3508 static long pModDeg(poly p, ring r)
3509 {
3510  long d=pOldFDeg(p, r);
3511  int c=p_GetComp(p, r);
3512  if ((c>0) && ((r->pModW)->range(c-1))) d+= (*(r->pModW))[c-1];
3513  return d;
3514  //return pOldFDeg(p, r)+(*pModW)[p_GetComp(p, r)-1];
3515 }
3516 
3517 void p_SetModDeg(intvec *w, ring r)
3518 {
3519  if (w!=NULL)
3520  {
3521  r->pModW = w;
3522  pOldFDeg = r->pFDeg;
3523  pOldLDeg = r->pLDeg;
3524  pOldLexOrder = r->pLexOrder;
3525  pSetDegProcs(r,pModDeg);
3526  r->pLexOrder = TRUE;
3527  }
3528  else
3529  {
3530  r->pModW = NULL;
3531  pRestoreDegProcs(r,pOldFDeg, pOldLDeg);
3532  r->pLexOrder = pOldLexOrder;
3533  }
3534 }
3535 
3536 /*2
3537 * handle memory request for sets of polynomials (ideals)
3538 * l is the length of *p, increment is the difference (may be negative)
3539 */
3540 void pEnlargeSet(poly* *p, int l, int increment)
3541 {
3542  poly* h;
3543 
3544  if (*p==NULL)
3545  {
3546  if (increment==0) return;
3547  h=(poly*)omAlloc0(increment*sizeof(poly));
3548  }
3549  else
3550  {
3551  h=(poly*)omReallocSize((poly*)*p,l*sizeof(poly),(l+increment)*sizeof(poly));
3552  if (increment>0)
3553  {
3554  //for (i=l; i<l+increment; i++)
3555  // h[i]=NULL;
3556  memset(&(h[l]),0,increment*sizeof(poly));
3557  }
3558  }
3559  *p=h;
3560 }
3561 
3562 /*2
3563 *divides p1 by its leading coefficient
3564 */
3565 void p_Norm(poly p1, const ring r)
3566 {
3567 #ifdef HAVE_RINGS
3568  if (rField_is_Ring(r))
3569  {
3570  if (!n_IsUnit(pGetCoeff(p1), r->cf)) return;
3571  // Werror("p_Norm not possible in the case of coefficient rings.");
3572  }
3573  else
3574 #endif
3575  if (p1!=NULL)
3576  {
3577  if (pNext(p1)==NULL)
3578  {
3579  p_SetCoeff(p1,n_Init(1,r->cf),r);
3580  return;
3581  }
3582  poly h;
3583  if (!n_IsOne(pGetCoeff(p1),r->cf))
3584  {
3585  number k, c;
3586  n_Normalize(pGetCoeff(p1),r->cf);
3587  k = pGetCoeff(p1);
3588  c = n_Init(1,r->cf);
3589  pSetCoeff0(p1,c);
3590  h = pNext(p1);
3591  while (h!=NULL)
3592  {
3593  c=n_Div(pGetCoeff(h),k,r->cf);
3594  // no need to normalize: Z/p, R
3595  // normalize already in nDiv: Q_a, Z/p_a
3596  // remains: Q
3597  if (rField_is_Q(r) && (!n_IsOne(c,r->cf))) n_Normalize(c,r->cf);
3598  p_SetCoeff(h,c,r);
3599  pIter(h);
3600  }
3601  n_Delete(&k,r->cf);
3602  }
3603  else
3604  {
3605  //if (r->cf->cfNormalize != nDummy2) //TODO: OPTIMIZE
3606  {
3607  h = pNext(p1);
3608  while (h!=NULL)
3609  {
3610  n_Normalize(pGetCoeff(h),r->cf);
3611  pIter(h);
3612  }
3613  }
3614  }
3615  }
3616 }
3617 
3618 /*2
3619 *normalize all coefficients
3620 */
3621 void p_Normalize(poly p,const ring r)
3622 {
3623  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
3624  while (p!=NULL)
3625  {
3626  // no test befor n_Normalize: n_Normalize should fix problems
3627  n_Normalize(pGetCoeff(p),r->cf);
3628  pIter(p);
3629  }
3630 }
3631 
3632 // splits p into polys with Exp(n) == 0 and Exp(n) != 0
3633 // Poly with Exp(n) != 0 is reversed
3634 static void p_SplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero, const ring r)
3635 {
3636  if (p == NULL)
3637  {
3638  *non_zero = NULL;
3639  *zero = NULL;
3640  return;
3641  }
3642  spolyrec sz;
3643  poly z, n_z, next;
3644  z = &sz;
3645  n_z = NULL;
3646 
3647  while(p != NULL)
3648  {
3649  next = pNext(p);
3650  if (p_GetExp(p, n,r) == 0)
3651  {
3652  pNext(z) = p;
3653  pIter(z);
3654  }
3655  else
3656  {
3657  pNext(p) = n_z;
3658  n_z = p;
3659  }
3660  p = next;
3661  }
3662  pNext(z) = NULL;
3663  *zero = pNext(&sz);
3664  *non_zero = n_z;
3665 }
3666 /*3
3667 * substitute the n-th variable by 1 in p
3668 * destroy p
3669 */
3670 static poly p_Subst1 (poly p,int n, const ring r)
3671 {
3672  poly qq=NULL, result = NULL;
3673  poly zero=NULL, non_zero=NULL;
3674 
3675  // reverse, so that add is likely to be linear
3676  p_SplitAndReversePoly(p, n, &non_zero, &zero,r);
3677 
3678  while (non_zero != NULL)
3679  {
3680  assume(p_GetExp(non_zero, n,r) != 0);
3681  qq = non_zero;
3682  pIter(non_zero);
3683  qq->next = NULL;
3684  p_SetExp(qq,n,0,r);
3685  p_Setm(qq,r);
3686  result = p_Add_q(result,qq,r);
3687  }
3688  p = p_Add_q(result, zero,r);
3689  p_Test(p,r);
3690  return p;
3691 }
3692 
3693 /*3
3694 * substitute the n-th variable by number e in p
3695 * destroy p
3696 */
3697 static poly p_Subst2 (poly p,int n, number e, const ring r)
3698 {
3699  assume( ! n_IsZero(e,r->cf) );
3700  poly qq,result = NULL;
3701  number nn, nm;
3702  poly zero, non_zero;
3703 
3704  // reverse, so that add is likely to be linear
3705  p_SplitAndReversePoly(p, n, &non_zero, &zero,r);
3706 
3707  while (non_zero != NULL)
3708  {
3709  assume(p_GetExp(non_zero, n, r) != 0);
3710  qq = non_zero;
3711  pIter(non_zero);
3712  qq->next = NULL;
3713  n_Power(e, p_GetExp(qq, n, r), &nn,r->cf);
3714  nm = n_Mult(nn, pGetCoeff(qq),r->cf);
3715 #ifdef HAVE_RINGS
3716  if (n_IsZero(nm,r->cf))
3717  {
3718  p_LmFree(&qq,r);
3719  n_Delete(&nm,r->cf);
3720  }
3721  else
3722 #endif
3723  {
3724  p_SetCoeff(qq, nm,r);
3725  p_SetExp(qq, n, 0,r);
3726  p_Setm(qq,r);
3727  result = p_Add_q(result,qq,r);
3728  }
3729  n_Delete(&nn,r->cf);
3730  }
3731  p = p_Add_q(result, zero,r);
3732  p_Test(p,r);
3733  return p;
3734 }
3735 
3736 
3737 /* delete monoms whose n-th exponent is different from zero */
3738 static poly p_Subst0(poly p, int n, const ring r)
3739 {
3740  spolyrec res;
3741  poly h = &res;
3742  pNext(h) = p;
3743 
3744  while (pNext(h)!=NULL)
3745  {
3746  if (p_GetExp(pNext(h),n,r)!=0)
3747  {
3748  p_LmDelete(&pNext(h),r);
3749  }
3750  else
3751  {
3752  pIter(h);
3753  }
3754  }
3755  p_Test(pNext(&res),r);
3756  return pNext(&res);
3757 }
3758 
3759 /*2
3760 * substitute the n-th variable by e in p
3761 * destroy p
3762 */
3763 poly p_Subst(poly p, int n, poly e, const ring r)
3764 {
3765  if (e == NULL) return p_Subst0(p, n,r);
3766 
3767  if (p_IsConstant(e,r))
3768  {
3769  if (n_IsOne(pGetCoeff(e),r->cf)) return p_Subst1(p,n,r);
3770  else return p_Subst2(p, n, pGetCoeff(e),r);
3771  }
3772 
3773 #ifdef HAVE_PLURAL
3774  if (rIsPluralRing(r))
3775  {
3776  return nc_pSubst(p,n,e,r);
3777  }
3778 #endif
3779 
3780  int exponent,i;
3781  poly h, res, m;
3782  int *me,*ee;
3783  number nu,nu1;
3784 
3785  me=(int *)omAlloc((rVar(r)+1)*sizeof(int));
3786  ee=(int *)omAlloc((rVar(r)+1)*sizeof(int));
3787  if (e!=NULL) p_GetExpV(e,ee,r);
3788  res=NULL;
3789  h=p;
3790  while (h!=NULL)
3791  {
3792  if ((e!=NULL) || (p_GetExp(h,n,r)==0))
3793  {
3794  m=p_Head(h,r);
3795  p_GetExpV(m,me,r);
3796  exponent=me[n];
3797  me[n]=0;
3798  for(i=rVar(r);i>0;i--)
3799  me[i]+=exponent*ee[i];
3800  p_SetExpV(m,me,r);
3801  if (e!=NULL)
3802  {
3803  n_Power(pGetCoeff(e),exponent,&nu,r->cf);
3804  nu1=n_Mult(pGetCoeff(m),nu,r->cf);
3805  n_Delete(&nu,r->cf);
3806  p_SetCoeff(m,nu1,r);
3807  }
3808  res=p_Add_q(res,m,r);
3809  }
3810  p_LmDelete(&h,r);
3811  }
3812  omFreeSize((ADDRESS)me,(rVar(r)+1)*sizeof(int));
3813  omFreeSize((ADDRESS)ee,(rVar(r)+1)*sizeof(int));
3814  return res;
3815 }
3816 
3817 /*2
3818  * returns a re-ordered convertion of a number as a polynomial,
3819  * with permutation of parameters
3820  * NOTE: this only works for Frank's alg. & trans. fields
3821  */
3822 poly n_PermNumber(const number z, const int *par_perm, const int , const ring src, const ring dst)
3823 {
3824 #if 0
3825  PrintS("\nSource Ring: \n");
3826  rWrite(src);
3827 
3828  if(0)
3829  {
3830  number zz = n_Copy(z, src->cf);
3831  PrintS("z: "); n_Write(zz, src);
3832  n_Delete(&zz, src->cf);
3833  }
3834 
3835  PrintS("\nDestination Ring: \n");
3836  rWrite(dst);
3837 
3838  /*Print("\nOldPar: %d\n", OldPar);
3839  for( int i = 1; i <= OldPar; i++ )
3840  {
3841  Print("par(%d) -> par/var (%d)\n", i, par_perm[i-1]);
3842  }*/
3843 #endif
3844  if( z == NULL )
3845  return NULL;
3846 
3847  const coeffs srcCf = src->cf;
3848  assume( srcCf != NULL );
3849 
3850  assume( !nCoeff_is_GF(srcCf) );
3851  assume( src->cf->extRing!=NULL );
3852 
3853  poly zz = NULL;
3854 
3855  const ring srcExtRing = srcCf->extRing;
3856  assume( srcExtRing != NULL );
3857 
3858  const coeffs dstCf = dst->cf;
3859  assume( dstCf != NULL );
3860 
3861  if( nCoeff_is_algExt(srcCf) ) // nCoeff_is_GF(srcCf)?
3862  {
3863  zz = (poly) z;
3864  if( zz == NULL ) return NULL;
3865  }
3866  else if (nCoeff_is_transExt(srcCf))
3867  {
3868  assume( !IS0(z) );
3869 
3870  zz = NUM((fraction)z);
3871  p_Test (zz, srcExtRing);
3872 
3873  if( zz == NULL ) return NULL;
3874  if( !DENIS1((fraction)z) )
3875  {
3876  if (!p_IsConstant(DEN((fraction)z),srcExtRing))
3877  WarnS("Not defined: Cannot map a rational fraction and make a polynomial out of it! Ignoring the denumerator.");
3878  }
3879  }
3880  else
3881  {
3882  assume (FALSE);
3883  Werror("Number permutation is not implemented for this data yet!");
3884  return NULL;
3885  }
3886 
3887  assume( zz != NULL );
3888  p_Test (zz, srcExtRing);
3889 
3890  nMapFunc nMap = n_SetMap(srcExtRing->cf, dstCf);
3891 
3892  assume( nMap != NULL );
3893 
3894  poly qq;
3895  if ((par_perm == NULL) && (rPar(dst) != 0 && rVar (srcExtRing) > 0))
3896  {
3897  int* perm;
3898  perm=(int *)omAlloc0((rVar(srcExtRing)+1)*sizeof(int));
3899  perm[0]= 0;
3900  for(int i=si_min(rVar(srcExtRing),rPar(dst));i>0;i--)
3901  perm[i]=-i;
3902  qq = p_PermPoly(zz, perm, srcExtRing, dst, nMap, NULL, rVar(srcExtRing)-1);
3903  omFreeSize ((ADDRESS)perm, (rVar(srcExtRing)+1)*sizeof(int));
3904  }
3905  else
3906  qq = p_PermPoly(zz, par_perm-1, srcExtRing, dst, nMap, NULL, rVar (srcExtRing)-1);
3907 
3908  if(nCoeff_is_transExt(srcCf)
3909  && (!DENIS1((fraction)z))
3910  && p_IsConstant(DEN((fraction)z),srcExtRing))
3911  {
3912  number n=nMap(pGetCoeff(DEN((fraction)z)),srcExtRing->cf, dstCf);
3913  qq=p_Div_nn(qq,n,dst);
3914  n_Delete(&n,dstCf);
3915  p_Normalize(qq,dst);
3916  }
3917  p_Test (qq, dst);
3918 
3919  return qq;
3920 }
3921 
3922 
3923 /*2
3924 *returns a re-ordered copy of a polynomial, with permutation of the variables
3925 */
3926 poly p_PermPoly (poly p, const int * perm, const ring oldRing, const ring dst,
3927  nMapFunc nMap, const int *par_perm, int OldPar)
3928 {
3929 #if 0
3930  p_Test(p, oldRing);
3931  PrintS("\np_PermPoly::p: "); p_Write(p, oldRing, oldRing); PrintLn();
3932 #endif
3933  const int OldpVariables = rVar(oldRing);
3934  poly result = NULL;
3935  poly result_last = NULL;
3936  poly aq = NULL; /* the map coefficient */
3937  poly qq; /* the mapped monomial */
3938  assume(dst != NULL);
3939  assume(dst->cf != NULL);
3940  while (p != NULL)
3941  {
3942  // map the coefficient
3943  if ( ((OldPar == 0) || (par_perm == NULL) || rField_is_GF(oldRing)) && (nMap != NULL) )
3944  {
3945  qq = p_Init(dst);
3946  assume( nMap != NULL );
3947  number n = nMap(p_GetCoeff(p, oldRing), oldRing->cf, dst->cf);
3948  n_Test (n,dst->cf);
3949  if ( nCoeff_is_algExt(dst->cf) )
3950  n_Normalize(n, dst->cf);
3951  p_GetCoeff(qq, dst) = n;// Note: n can be a ZERO!!!
3952  }
3953  else
3954  {
3955  qq = p_One(dst);
3956 // aq = naPermNumber(p_GetCoeff(p, oldRing), par_perm, OldPar, oldRing); // no dst???
3957 // poly n_PermNumber(const number z, const int *par_perm, const int P, const ring src, const ring dst)
3958  aq = n_PermNumber(p_GetCoeff(p, oldRing), par_perm, OldPar, oldRing, dst);
3959  p_Test(aq, dst);
3960  if ( nCoeff_is_algExt(dst->cf) )
3961  p_Normalize(aq,dst);
3962  if (aq == NULL)
3963  p_SetCoeff(qq, n_Init(0, dst->cf),dst); // Very dirty trick!!!
3964  p_Test(aq, dst);
3965  }
3966  if (rRing_has_Comp(dst))
3967  p_SetComp(qq, p_GetComp(p, oldRing), dst);
3968  if ( n_IsZero(pGetCoeff(qq), dst->cf) )
3969  {
3970  p_LmDelete(&qq,dst);
3971  qq = NULL;
3972  }
3973  else
3974  {
3975  // map pars:
3976  int mapped_to_par = 0;
3977  for(int i = 1; i <= OldpVariables; i++)
3978  {
3979  int e = p_GetExp(p, i, oldRing);
3980  if (e != 0)
3981  {
3982  if (perm==NULL)
3983  p_SetExp(qq, i, e, dst);
3984  else if (perm[i]>0)
3985  p_AddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/, dst);
3986  else if (perm[i]<0)
3987  {
3988  number c = p_GetCoeff(qq, dst);
3989  if (rField_is_GF(dst))
3990  {
3991  assume( dst->cf->extRing == NULL );
3992  number ee = n_Param(1, dst);
3993  number eee;
3994  n_Power(ee, e, &eee, dst->cf); //nfDelete(ee,dst);
3995  ee = n_Mult(c, eee, dst->cf);
3996  //nfDelete(c,dst);nfDelete(eee,dst);
3997  pSetCoeff0(qq,ee);
3998  }
3999  else if (nCoeff_is_Extension(dst->cf))
4000  {
4001  const int par = -perm[i];
4002  assume( par > 0 );
4003 // WarnS("longalg missing 3");
4004 #if 1
4005  const coeffs C = dst->cf;
4006  assume( C != NULL );
4007  const ring R = C->extRing;
4008  assume( R != NULL );
4009  assume( par <= rVar(R) );
4010  poly pcn; // = (number)c
4011  assume( !n_IsZero(c, C) );
4012  if( nCoeff_is_algExt(C) )
4013  pcn = (poly) c;
4014  else // nCoeff_is_transExt(C)
4015  pcn = NUM((fraction)c);
4016  if (pNext(pcn) == NULL) // c->z
4017  p_AddExp(pcn, -perm[i], e, R);
4018  else /* more difficult: we have really to multiply: */
4019  {
4020  poly mmc = p_ISet(1, R);
4021  p_SetExp(mmc, -perm[i], e, R);
4022  p_Setm(mmc, R);
4023  number nnc;
4024  // convert back to a number: number nnc = mmc;
4025  if( nCoeff_is_algExt(C) )
4026  nnc = (number) mmc;
4027  else // nCoeff_is_transExt(C)
4028  nnc = ntInit(mmc, C);
4029  p_GetCoeff(qq, dst) = n_Mult((number)c, nnc, C);
4030  n_Delete((number *)&c, C);
4031  n_Delete((number *)&nnc, C);
4032  }
4033  mapped_to_par=1;
4034 #endif
4035  }
4036  }
4037  else
4038  {
4039  /* this variable maps to 0 !*/
4040  p_LmDelete(&qq, dst);
4041  break;
4042  }
4043  }
4044  }
4045  if ( mapped_to_par && qq!= NULL && nCoeff_is_algExt(dst->cf) )
4046  {
4047  number n = p_GetCoeff(qq, dst);
4048  n_Normalize(n, dst->cf);
4049  p_GetCoeff(qq, dst) = n;
4050  }
4051  }
4052  pIter(p);
4053 
4054 #if 0
4055  p_Test(aq,dst);
4056  PrintS("\naq: "); p_Write(aq, dst, dst); PrintLn();
4057 #endif
4058 
4059 
4060 #if 1
4061  if (qq!=NULL)
4062  {
4063  p_Setm(qq,dst);
4064 
4065  p_Test(aq,dst);
4066  p_Test(qq,dst);
4067 
4068 #if 0
4069  p_Test(qq,dst);
4070  PrintS("\nqq: "); p_Write(qq, dst, dst); PrintLn();
4071 #endif
4072 
4073  if (aq!=NULL)
4074  qq=p_Mult_q(aq,qq,dst);
4075  aq = qq;
4076  while (pNext(aq) != NULL) pIter(aq);
4077  if (result_last==NULL)
4078  {
4079  result=qq;
4080  }
4081  else
4082  {
4083  pNext(result_last)=qq;
4084  }
4085  result_last=aq;
4086  aq = NULL;
4087  }
4088  else if (aq!=NULL)
4089  {
4090  p_Delete(&aq,dst);
4091  }
4092  }
4093  result=p_SortAdd(result,dst);
4094 #else
4095  // if (qq!=NULL)
4096  // {
4097  // pSetm(qq);
4098  // pTest(qq);
4099  // pTest(aq);
4100  // if (aq!=NULL) qq=pMult(aq,qq);
4101  // aq = qq;
4102  // while (pNext(aq) != NULL) pIter(aq);
4103  // pNext(aq) = result;
4104  // aq = NULL;
4105  // result = qq;
4106  // }
4107  // else if (aq!=NULL)
4108  // {
4109  // pDelete(&aq);
4110  // }
4111  //}
4112  //p = result;
4113  //result = NULL;
4114  //while (p != NULL)
4115  //{
4116  // qq = p;
4117  // pIter(p);
4118  // qq->next = NULL;
4119  // result = pAdd(result, qq);
4120  //}
4121 #endif
4122  p_Test(result,dst);
4123 #if 0
4124  p_Test(result,dst);
4125  PrintS("\nresult: "); p_Write(result,dst,dst); PrintLn();
4126 #endif
4127  return result;
4128 }
4129 /**************************************************************
4130  *
4131  * Jet
4132  *
4133  **************************************************************/
4134 
4135 poly pp_Jet(poly p, int m, const ring R)
4136 {
4137  poly r=NULL;
4138  poly t=NULL;
4139 
4140  while (p!=NULL)
4141  {
4142  if (p_Totaldegree(p,R)<=m)
4143  {
4144  if (r==NULL)
4145  r=p_Head(p,R);
4146  else
4147  if (t==NULL)
4148  {
4149  pNext(r)=p_Head(p,R);
4150  t=pNext(r);
4151  }
4152  else
4153  {
4154  pNext(t)=p_Head(p,R);
4155  pIter(t);
4156  }
4157  }
4158  pIter(p);
4159  }
4160  return r;
4161 }
4162 
4163 poly p_Jet(poly p, int m,const ring R)
4164 {
4165  while((p!=NULL) && (p_Totaldegree(p,R)>m)) p_LmDelete(&p,R);
4166  if (p==NULL) return NULL;
4167  poly r=p;
4168  while (pNext(p)!=NULL)
4169  {
4170  if (p_Totaldegree(pNext(p),R)>m)
4171  {
4172  p_LmDelete(&pNext(p),R);
4173  }
4174  else
4175  pIter(p);
4176  }
4177  return r;
4178 }
4179 
4180 poly pp_JetW(poly p, int m, short *w, const ring R)
4181 {
4182  poly r=NULL;
4183  poly t=NULL;
4184  while (p!=NULL)
4185  {
4186  if (totaldegreeWecart_IV(p,R,w)<=m)
4187  {
4188  if (r==NULL)
4189  r=p_Head(p,R);
4190  else
4191  if (t==NULL)
4192  {
4193  pNext(r)=p_Head(p,R);
4194  t=pNext(r);
4195  }
4196  else
4197  {
4198  pNext(t)=p_Head(p,R);
4199  pIter(t);
4200  }
4201  }
4202  pIter(p);
4203  }
4204  return r;
4205 }
4206 
4207 poly p_JetW(poly p, int m, short *w, const ring R)
4208 {
4209  while((p!=NULL) && (totaldegreeWecart_IV(p,R,w)>m)) p_LmDelete(&p,R);
4210  if (p==NULL) return NULL;
4211  poly r=p;
4212  while (pNext(p)!=NULL)
4213  {
4214  if (totaldegreeWecart_IV(pNext(p),R,w)>m)
4215  {
4216  p_LmDelete(&pNext(p),R);
4217  }
4218  else
4219  pIter(p);
4220  }
4221  return r;
4222 }
4223 
4224 /*************************************************************/
4225 int p_MinDeg(poly p,intvec *w, const ring R)
4226 {
4227  if(p==NULL)
4228  return -1;
4229  int d=-1;
4230  while(p!=NULL)
4231  {
4232  int d0=0;
4233  for(int j=0;j<rVar(R);j++)
4234  if(w==NULL||j>=w->length())
4235  d0+=p_GetExp(p,j+1,R);
4236  else
4237  d0+=(*w)[j]*p_GetExp(p,j+1,R);
4238  if(d0<d||d==-1)
4239  d=d0;
4240  pIter(p);
4241  }
4242  return d;
4243 }
4244 
4245 /***************************************************************/
4246 
4247 poly p_Series(int n,poly p,poly u, intvec *w, const ring R)
4248 {
4249  short *ww=iv2array(w,R);
4250  if(p!=NULL)
4251  {
4252  if(u==NULL)
4253  p=p_JetW(p,n,ww,R);
4254  else
4255  p=p_JetW(p_Mult_q(p,p_Invers(n-p_MinDeg(p,w,R),u,w,R),R),n,ww,R);
4256  }
4257  omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(short));
4258  return p;
4259 }
4260 
4261 poly p_Invers(int n,poly u,intvec *w, const ring R)
4262 {
4263  if(n<0)
4264  return NULL;
4265  number u0=n_Invers(pGetCoeff(u),R->cf);
4266  poly v=p_NSet(u0,R);
4267  if(n==0)
4268  return v;
4269  short *ww=iv2array(w,R);
4270  poly u1=p_JetW(p_Sub(p_One(R),p_Mult_nn(u,u0,R),R),n,ww,R);
4271  if(u1==NULL)
4272  {
4273  omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(short));
4274  return v;
4275  }
4276  poly v1=p_Mult_nn(p_Copy(u1,R),u0,R);
4277  v=p_Add_q(v,p_Copy(v1,R),R);
4278  for(int i=n/p_MinDeg(u1,w,R);i>1;i--)
4279  {
4280  v1=p_JetW(p_Mult_q(v1,p_Copy(u1,R),R),n,ww,R);
4281  v=p_Add_q(v,p_Copy(v1,R),R);
4282  }
4283  p_Delete(&u1,R);
4284  p_Delete(&v1,R);
4285  omFreeSize((ADDRESS)ww,(rVar(R)+1)*sizeof(short));
4286  return v;
4287 }
4288 
4289 BOOLEAN p_EqualPolys(poly p1,poly p2, const ring r)
4290 {
4291  while ((p1 != NULL) && (p2 != NULL))
4292  {
4293  if (! p_LmEqual(p1, p2,r))
4294  return FALSE;
4295  if (! n_Equal(p_GetCoeff(p1,r), p_GetCoeff(p2,r),r->cf ))
4296  return FALSE;
4297  pIter(p1);
4298  pIter(p2);
4299  }
4300  return (p1==p2);
4301 }
4302 
4303 static inline BOOLEAN p_ExpVectorEqual(poly p1, poly p2, const ring r1, const ring r2)
4304 {
4305  assume( r1 == r2 || rSamePolyRep(r1, r2) );
4306 
4307  p_LmCheckPolyRing1(p1, r1);
4308  p_LmCheckPolyRing1(p2, r2);
4309 
4310  int i = r1->ExpL_Size;
4311 
4312  assume( r1->ExpL_Size == r2->ExpL_Size );
4313 
4314  unsigned long *ep = p1->exp;
4315  unsigned long *eq = p2->exp;
4316 
4317  do
4318  {
4319  i--;
4320  if (ep[i] != eq[i]) return FALSE;
4321  }
4322  while (i);
4323 
4324  return TRUE;
4325 }
4326 
4327 BOOLEAN p_EqualPolys(poly p1,poly p2, const ring r1, const ring r2)
4328 {
4329  assume( r1 == r2 || rSamePolyRep(r1, r2) ); // will be used in rEqual!
4330  assume( r1->cf == r2->cf );
4331 
4332  while ((p1 != NULL) && (p2 != NULL))
4333  {
4334  // returns 1 if ExpVector(p)==ExpVector(q): does not compare numbers !!
4335  // #define p_LmEqual(p1, p2, r) p_ExpVectorEqual(p1, p2, r)
4336 
4337  if (! p_ExpVectorEqual(p1, p2, r1, r2))
4338  return FALSE;
4339 
4340  if (! n_Equal(p_GetCoeff(p1,r1), p_GetCoeff(p2,r2), r1->cf ))
4341  return FALSE;
4342 
4343  pIter(p1);
4344  pIter(p2);
4345  }
4346  return (p1==p2);
4347 }
4348 
4349 /*2
4350 *returns TRUE if p1 is a skalar multiple of p2
4351 *assume p1 != NULL and p2 != NULL
4352 */
4353 BOOLEAN p_ComparePolys(poly p1,poly p2, const ring r)
4354 {
4355  number n,nn;
4356  pAssume(p1 != NULL && p2 != NULL);
4357 
4358  if (!p_LmEqual(p1,p2,r)) //compare leading mons
4359  return FALSE;
4360  if ((pNext(p1)==NULL) && (pNext(p2)!=NULL))
4361  return FALSE;
4362  if ((pNext(p2)==NULL) && (pNext(p1)!=NULL))
4363  return FALSE;
4364  if (pLength(p1) != pLength(p2))
4365  return FALSE;
4366 #ifdef HAVE_RINGS
4367  if (rField_is_Ring(r))
4368  {
4369  if (!n_DivBy(p_GetCoeff(p1, r), p_GetCoeff(p2, r), r)) return FALSE;
4370  }
4371 #endif
4372  n=n_Div(p_GetCoeff(p1,r),p_GetCoeff(p2,r),r);
4373  while ((p1 != NULL) /*&& (p2 != NULL)*/)
4374  {
4375  if ( ! p_LmEqual(p1, p2,r))
4376  {
4377  n_Delete(&n, r);
4378  return FALSE;
4379  }
4380  if (!n_Equal(p_GetCoeff(p1, r), nn = n_Mult(p_GetCoeff(p2, r),n, r->cf), r->cf))
4381  {
4382  n_Delete(&n, r);
4383  n_Delete(&nn, r);
4384  return FALSE;
4385  }
4386  n_Delete(&nn, r);
4387  pIter(p1);
4388  pIter(p2);
4389  }
4390  n_Delete(&n, r);
4391  return TRUE;
4392 }
4393 
4394 /*2
4395 * returns the length of a (numbers of monomials)
4396 * respect syzComp
4397 */
4398 poly p_Last(const poly p, int &l, const ring r)
4399 {
4400  if (p == NULL)
4401  {
4402  l = 0;
4403  return NULL;
4404  }
4405  l = 1;
4406  poly a = p;
4407  if (! rIsSyzIndexRing(r))
4408  {
4409  poly next = pNext(a);
4410  while (next!=NULL)
4411  {
4412  a = next;
4413  next = pNext(a);
4414  l++;
4415  }
4416  }
4417  else
4418  {
4419  int curr_limit = rGetCurrSyzLimit(r);
4420  poly pp = a;
4421  while ((a=pNext(a))!=NULL)
4422  {
4423  if (p_GetComp(a,r)<=curr_limit/*syzComp*/)
4424  l++;
4425  else break;
4426  pp = a;
4427  }
4428  a=pp;
4429  }
4430  return a;
4431 }
4432 
4433 int p_Var(poly m,const ring r)
4434 {
4435  if (m==NULL) return 0;
4436  if (pNext(m)!=NULL) return 0;
4437  int i,e=0;
4438  for (i=rVar(r); i>0; i--)
4439  {
4440  int exp=p_GetExp(m,i,r);
4441  if (exp==1)
4442  {
4443  if (e==0) e=i;
4444  else return 0;
4445  }
4446  else if (exp!=0)
4447  {
4448  return 0;
4449  }
4450  }
4451  return e;
4452 }
4453 
4454 /*2
4455 *the minimal index of used variables - 1
4456 */
4457 int p_LowVar (poly p, const ring r)
4458 {
4459  int k,l,lex;
4460 
4461  if (p == NULL) return -1;
4462 
4463  k = 32000;/*a very large dummy value*/
4464  while (p != NULL)
4465  {
4466  l = 1;
4467  lex = p_GetExp(p,l,r);
4468  while ((l < (rVar(r))) && (lex == 0))
4469  {
4470  l++;
4471  lex = p_GetExp(p,l,r);
4472  }
4473  l--;
4474  if (l < k) k = l;
4475  pIter(p);
4476  }
4477  return k;
4478 }
4479 
4480 /*2
4481 * verschiebt die Indizees der Modulerzeugenden um i
4482 */
4483 void p_Shift (poly * p,int i, const ring r)
4484 {
4485  poly qp1 = *p,qp2 = *p;/*working pointers*/
4486  int j = p_MaxComp(*p,r),k = p_MinComp(*p,r);
4487 
4488  if (j+i < 0) return ;
4489  while (qp1 != NULL)
4490  {
4491  if ((p_GetComp(qp1,r)+i > 0) || ((j == -i) && (j == k)))
4492  {
4493  p_AddComp(qp1,i,r);
4494  p_SetmComp(qp1,r);
4495  qp2 = qp1;
4496  pIter(qp1);
4497  }
4498  else
4499  {
4500  if (qp2 == *p)
4501  {
4502  pIter(*p);
4503  p_LmDelete(&qp2,r);
4504  qp2 = *p;
4505  qp1 = *p;
4506  }
4507  else
4508  {
4509  qp2->next = qp1->next;
4510  if (qp1!=NULL) p_LmDelete(&qp1,r);
4511  qp1 = qp2->next;
4512  }
4513  }
4514  }
4515 }
4516 
4517 /***************************************************************
4518  *
4519  * Storage Managament Routines
4520  *
4521  ***************************************************************/
4522 
4523 
4524 static inline unsigned long GetBitFields(const long e,
4525  const unsigned int s, const unsigned int n)
4526 {
4527 #define Sy_bit_L(x) (((unsigned long)1L)<<(x))
4528  unsigned int i = 0;
4529  unsigned long ev = 0L;
4530  assume(n > 0 && s < BIT_SIZEOF_LONG);
4531  do
4532  {
4533  assume(s+i < BIT_SIZEOF_LONG);
4534  if (e > (long) i) ev |= Sy_bit_L(s+i);
4535  else break;
4536  i++;
4537  }
4538  while (i < n);
4539  return ev;
4540 }
4541 
4542 // Short Exponent Vectors are used for fast divisibility tests
4543 // ShortExpVectors "squeeze" an exponent vector into one word as follows:
4544 // Let n = BIT_SIZEOF_LONG / pVariables.
4545 // If n == 0 (i.e. pVariables > BIT_SIZE_OF_LONG), let m == the number
4546 // of non-zero exponents. If (m>BIT_SIZEOF_LONG), then sev = ~0, else
4547 // first m bits of sev are set to 1.
4548 // Otherwise (i.e. pVariables <= BIT_SIZE_OF_LONG)
4549 // represented by a bit-field of length n (resp. n+1 for some
4550 // exponents). If the value of an exponent is greater or equal to n, then
4551 // all of its respective n bits are set to 1. If the value of an exponent
4552 // is smaller than n, say m, then only the first m bits of the respective
4553 // n bits are set to 1, the others are set to 0.
4554 // This way, we have:
4555 // exp1 / exp2 ==> (ev1 & ~ev2) == 0, i.e.,
4556 // if (ev1 & ~ev2) then exp1 does not divide exp2
4557 unsigned long p_GetShortExpVector(const poly p, const ring r)
4558 {
4559  assume(p != NULL);
4560  unsigned long ev = 0; // short exponent vector
4561  unsigned int n = BIT_SIZEOF_LONG / r->N; // number of bits per exp
4562  unsigned int m1; // highest bit which is filled with (n+1)
4563  int i=0,j=1;
4564 
4565  if (n == 0)
4566  {
4567  if (r->N <2*BIT_SIZEOF_LONG)
4568  {
4569  n=1;
4570  m1=0;
4571  }
4572  else
4573  {
4574  for (; j<=r->N; j++)
4575  {
4576  if (p_GetExp(p,j,r) > 0) i++;
4577  if (i == BIT_SIZEOF_LONG) break;
4578  }
4579  if (i>0)
4580  ev = ~(0UL) >> (BIT_SIZEOF_LONG - i);
4581  return ev;
4582  }
4583  }
4584  else
4585  {
4586  m1 = (n+1)*(BIT_SIZEOF_LONG - n*r->N);
4587  }
4588 
4589  n++;
4590  while (i<m1)
4591  {
4592  ev |= GetBitFields(p_GetExp(p, j,r), i, n);
4593  i += n;
4594  j++;
4595  }
4596 
4597  n--;
4598  while (i<BIT_SIZEOF_LONG)
4599  {
4600  ev |= GetBitFields(p_GetExp(p, j,r), i, n);
4601  i += n;
4602  j++;
4603  }
4604  return ev;
4605 }
4606 
4607 
4608 /// p_GetShortExpVector of p * pp
4609 unsigned long p_GetShortExpVector(const poly p, const poly pp, const ring r)
4610 {
4611  assume(p != NULL);
4612  assume(pp != NULL);
4613 
4614  unsigned long ev = 0; // short exponent vector
4615  unsigned int n = BIT_SIZEOF_LONG / r->N; // number of bits per exp
4616  unsigned int m1; // highest bit which is filled with (n+1)
4617  int j=1;
4618  unsigned long i = 0L;
4619 
4620  if (n == 0)
4621  {
4622  if (r->N <2*BIT_SIZEOF_LONG)
4623  {
4624  n=1;
4625  m1=0;
4626  }
4627  else
4628  {
4629  for (; j<=r->N; j++)
4630  {
4631  if (p_GetExp(p,j,r) > 0 || p_GetExp(pp,j,r) > 0) i++;
4632  if (i == BIT_SIZEOF_LONG) break;
4633  }
4634  if (i>0)
4635  ev = ~(0UL) >> (BIT_SIZEOF_LONG - i);
4636  return ev;
4637  }
4638  }
4639  else
4640  {
4641  m1 = (n+1)*(BIT_SIZEOF_LONG - n*r->N);
4642  }
4643 
4644  n++;
4645  while (i<m1)
4646  {
4647  ev |= GetBitFields(p_GetExp(p, j,r) + p_GetExp(pp, j,r), i, n);
4648  i += n;
4649  j++;
4650  }
4651 
4652  n--;
4653  while (i<BIT_SIZEOF_LONG)
4654  {
4655  ev |= GetBitFields(p_GetExp(p, j,r) + p_GetExp(pp, j,r), i, n);
4656  i += n;
4657  j++;
4658  }
4659  return ev;
4660 }
4661 
4662 
4663 
4664 /***************************************************************
4665  *
4666  * p_ShallowDelete
4667  *
4668  ***************************************************************/
4669 #undef LINKAGE
4670 #define LINKAGE
4671 #undef p_Delete__T
4672 #define p_Delete__T p_ShallowDelete
4673 #undef n_Delete__T
4674 #define n_Delete__T(n, r) do {} while (0)
4675 
4677 
#define p_LmCheckPolyRing2(p, r)
Definition: monomials.h:207
#define n_New(n, r)
Definition: coeffs.h:441
int status int void size_t count
Definition: si_signals.h:59
BOOLEAN p_VectorHasUnitB(poly p, int *k, const ring r)
Definition: p_polys.cc:3209
void p_SetModDeg(intvec *w, ring r)
Definition: p_polys.cc:3517
for idElimination, like a, except pFDeg, pWeigths ignore it
Definition: ring.h:693
#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:609
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:533
LINLINE number nlSub(number la, number li, const coeffs r)
Definition: longrat.cc:2518
static unsigned long p_AddComp(poly p, unsigned long v, ring r)
Definition: p_polys.h:443
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:685
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:1809
static poly p_MonPower(poly p, int exp, const ring r)
Definition: p_polys.cc:1911
const char * eati(const char *s, int *i)
Definition: reporter.cc:390
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:1272
#define D(A)
Definition: gentable.cc:119
for int64 weights
Definition: ring.h:673
static BOOLEAN p_LmIsConstantComp(const poly p, const ring r)
Definition: p_polys.h:937
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:516
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:599
void PrintLn()
Definition: reporter.cc:327
#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 poly p_LmDeleteAndNext(poly p, const ring r)
Definition: p_polys.h:721
return
Definition: syzextra.cc:280
Definition: ring.h:61
long pLDeg1c_Totaldegree(poly p, int *l, const ring r)
Definition: p_polys.cc:1004
BOOLEAN p_ComparePolys(poly p1, poly p2, const ring r)
returns TRUE if p1 is a skalar multiple of p2 assume p1 != NULL and p2 != NULL
Definition: p_polys.cc:4353
BOOLEAN nlGreaterZero(number za, const coeffs r)
Definition: longrat.cc:1141
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition: coeffs.h:38
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:167
long pLDeg1c(poly p, int *l, const ring r)
Definition: p_polys.cc:876
#define FALSE
Definition: auxiliary.h:140
static BOOLEAN pOldLexOrder
Definition: p_polys.cc:3506
poly p_Homogen(poly p, int varnum, const ring r)
Definition: p_polys.cc:3138
void p_Split(poly p, poly *h)
Definition: p_polys.cc:1325
return P p
Definition: myNF.cc:203
opposite of ls
Definition: ring.h:694
void p_VectorHasUnit(poly p, int *k, int *len, const ring r)
Definition: p_polys.cc:3234
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:469
BOOLEAN p_IsHomogeneous(poly p, const ring r)
Definition: p_polys.cc:3187
f
Definition: cfModGcd.cc:4022
static int rPar(const ring r)
(r->cf->P)
Definition: ring.h:547
poly p_NSet(number n, const ring r)
returns the poly representing the number n, destroys n
Definition: p_polys.cc:1448
#define pAssume(cond)
Definition: monomials.h:98
static long p_IncrExp(poly p, int v, ring r)
Definition: p_polys.h:587
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:4225
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:236
static BOOLEAN rIsSyzIndexRing(const ring r)
Definition: ring.h:714
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:1845
#define p_GetComp(p, r)
Definition: monomials.h:72
static int rGetCurrSyzLimit(const ring r)
Definition: ring.h:717
static poly last
Definition: hdegree.cc:1075
#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:539
static long pModDeg(poly p, ring r)
Definition: p_polys.cc:3508
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1448
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
LINLINE number nlAdd(number la, number li, const coeffs r)
Definition: longrat.cc:2452
static poly p_Subst1(poly p, int n, const ring r)
Definition: p_polys.cc:3670
int rChar(ring r)
Definition: ring.cc:684
void nlInpGcd(number &a, number b, const coeffs r)
Definition: longrat.cc:2588
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:540
poly singclap_gcd(poly f, poly g, const ring r)
destroys f and g
Definition: clapsing.cc:287
long pLDeg0c(poly p, int *l, const ring r)
Definition: p_polys.cc:769
void sBucketDestroyAdd(sBucket_pt bucket, poly *p, int *length)
Definition: sbuckets.h:72
poly p_Div_nn(poly p, const number n, const ring r)
Definition: p_polys.cc:1480
long int64
Definition: auxiliary.h:112
void p_Lcm(const poly a, const poly b, poly m, const ring r)
Definition: p_polys.cc:1572
static BOOLEAN rField_is_Q_a(const ring r)
Definition: ring.h:488
number nlGcd(number a, number b, const coeffs r)
Definition: longrat.cc:1178
#define TRUE
Definition: auxiliary.h:144
int length() const
Definition: intvec.h:86
static void p_MonMult(poly p, poly q, const ring r)
Definition: p_polys.cc:1935
static long p_Totaldegree(poly p, const ring r)
Definition: p_polys.h:1435
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:579
void sBucket_Add_p(sBucket_pt bucket, poly p, int length)
adds poly p to bucket destroys p!
Definition: sbuckets.cc:206
void * ADDRESS
Definition: auxiliary.h:161
poly p_Subst(poly p, int n, poly e, const ring r)
Definition: p_polys.cc:3763
g
Definition: cfModGcd.cc:4031
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:716
poly p_TakeOutComp1(poly *p, int k, const ring r)
Definition: p_polys.cc:3266
static BOOLEAN rField_is_GF(const ring r)
Definition: ring.h:470
void p_Norm(poly p1, const ring r)
Definition: p_polys.cc:3565
static FORCE_INLINE BOOLEAN nCoeff_is_Q(const coeffs r)
Definition: coeffs.h:824
void p_SimpleContent(poly ph, int smax, const ring r)
Definition: p_polys.cc:2391
static long p_MultExp(poly p, int v, long ee, ring r)
Definition: p_polys.h:617
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:51
long pLDeg1_Deg(poly p, int *l, const ring r)
Definition: p_polys.cc:909
#define WarnS
Definition: emacs.cc:81
Definition: ring.h:66
poly singclap_pdivide(poly f, poly g, const ring r)
Definition: clapsing.cc:547
#define omAlloc(size)
Definition: omAllocDecl.h:210
static bool rIsPluralRing(const ring r)
we must always have this test!
Definition: ring.h:361
static poly p_TwoMonPower(poly p, int exp, const ring r)
Definition: p_polys.cc:2017
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:401
poly p_Sub(poly p1, poly p2, const ring r)
Definition: p_polys.cc:1901
long(* pLDegProc)(poly p, int *length, ring r)
Definition: ring.h:45
static void p_LmFree(poly p, ring)
Definition: p_polys.h:679
Definition: ring.h:64
static int pLength(poly a)
Definition: p_polys.h:189
static void p_SetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1456
union sro_ord::@0 data
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:811
poly pp
Definition: myNF.cc:296
static FORCE_INLINE BOOLEAN nCoeff_is_Q_a(const coeffs r)
Definition: coeffs.h:883
static long p_SubExp(poly p, int v, long ee, ring r)
Definition: p_polys.h:609
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:497
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:801
void p_Cleardenom_n(poly ph, const ring r, number &c)
Definition: p_polys.cc:2826
int p_IsUnivariate(poly p, const ring r)
return i, if poly depends only on var(i)
Definition: p_polys.cc:1252
#define pIter(p)
Definition: monomials.h:44
poly pp_JetW(poly p, int m, short *w, const ring R)
Definition: p_polys.cc:4180
poly res
Definition: myNF.cc:322
static poly p_SortAdd(poly p, const ring r, BOOLEAN revert=FALSE)
Definition: p_polys.h:1147
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:635
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:1353
ro_typ ord_typ
Definition: ring.h:182
void p_ContentRat(poly &ph, const ring r)
Definition: p_polys.cc:1655
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:932
static poly p_Head(poly p, const ring r)
Definition: p_polys.h:819
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:1588
int p_Size(poly p, const ring r)
Definition: p_polys.cc:3121
poly p_Invers(int n, poly u, intvec *w, const ring R)
Definition: p_polys.cc:4261
static int p_Comp_k_n(poly a, poly b, int k, ring r)
Definition: p_polys.h:636
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:911
static void p_SetCompP(poly p, int i, ring r)
Definition: p_polys.h:243
Definition: intvec.h:16
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:49
poly p_One(const ring r)
Definition: p_polys.cc:1318
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:465
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:3481
#define omFree(addr)
Definition: omAllocDecl.h:261
number ntInit(long i, const coeffs cf)
Definition: transext.cc:616
#define assume(x)
Definition: mod2.h:405
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition: p_polys.h:1784
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:1611
static poly p_MonMultC(poly p, poly q, const ring rr)
Definition: p_polys.cc:1955
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:3738
void p_Vec2Polys(poly v, poly **p, int *len, const ring r)
Definition: p_polys.cc:3459
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:72
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:771
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:782
const ring R
Definition: DebugPrint.cc:36
sBucket_pt sBucketCreate(ring r)
Definition: sbuckets.cc:125
static FORCE_INLINE void n_Write(number &n, const coeffs r, const BOOLEAN bShortOut=TRUE)
Definition: coeffs.h:592
#define n_Test(a, r)
BOOLEAN n_Test(number a, const coeffs r)
Definition: coeffs.h:923
poly p_Jet(poly p, int m, const ring R)
Definition: p_polys.cc:4163
static BOOLEAN p_DivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1685
void p_Setm_Dummy(poly p, const ring r)
Definition: p_polys.cc:540
Definition: ring.h:180
All the auxiliary stuff.
static int p_LmCmp(poly p, poly q, const ring r)
Definition: p_polys.h:1472
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:565
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:1681
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:558
static int si_max(const int a, const int b)
Definition: auxiliary.h:166
static FORCE_INLINE BOOLEAN nCoeff_is_transExt(const coeffs r)
TRUE iff r represents a transcendental extension field.
Definition: coeffs.h:919
static number * pnBin(int exp, const ring r)
Definition: p_polys.cc:1969
static number p_InitContent(poly ph, const ring r)
Definition: p_polys.cc:2449
poly p_GetCoeffRat(poly p, int ishift, ring r)
Definition: p_polys.cc:1633
int i
Definition: cfEzgcd.cc:123
Induced (Schreyer) ordering.
Definition: ring.h:695
void PrintS(const char *s)
Definition: reporter.cc:294
static poly p_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:902
void(* p_SetmProc)(poly p, const ring r)
Definition: ring.h:47
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:461
int p_LowVar(poly p, const ring r)
the minimal index of used variables - 1
Definition: p_polys.cc:4457
#define p_LmCheckPolyRing1(p, r)
Definition: monomials.h:185
static long p_MinComp(poly p, ring lmRing, ring tailRing)
Definition: p_polys.h:302
poly p_Divide(poly a, poly b, const ring r)
Definition: p_polys.cc:1467
static void p_ExpVectorAdd(poly p1, poly p2, const ring r)
Definition: p_polys.h:1339
S?
Definition: ring.h:677
BOOLEAN rOrd_SetCompRequiresSetm(const ring r)
return TRUE if p_SetComp requires p_Setm
Definition: ring.cc:1875
static poly p_LmFreeAndNext(poly p, ring)
Definition: p_polys.h:699
BOOLEAN pSetm_error
Definition: p_polys.cc:155
void rWrite(ring r, BOOLEAN details)
Definition: ring.cc:236
void p_Content(poly ph, const ring r)
Definition: p_polys.cc:2182
#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:465
static FORCE_INLINE BOOLEAN nCoeff_is_GF(const coeffs r)
Definition: coeffs.h:837
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:720
static poly p_Pow(poly p, int i, const ring r)
Definition: p_polys.cc:2082
poly p_JetW(poly p, int m, short *w, const ring R)
Definition: p_polys.cc:4207
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:1781
static short scaFirstAltVar(ring r)
Definition: sca.h:18
static poly pReverse(poly p)
Definition: p_polys.h:324
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition: coeffs.h:422
Definition: ring.h:69
poly p_Series(int n, poly p, poly u, intvec *w, const ring R)
Definition: p_polys.cc:4247
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition: p_polys.cc:4289
#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:449
void pRestoreDegProcs(ring r, pFDegProc old_FDeg, pLDegProc old_lDeg)
Definition: p_polys.cc:3493
Definition: ring.h:69
p_SetmProc p_GetSetmProc(ring r)
Definition: p_polys.cc:559
static long p_GetOrder(poly p, ring r)
Definition: p_polys.h:410
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:455
static FORCE_INLINE number n_Farey(number a, number b, const coeffs r)
Definition: coeffs.h:785
BOOLEAN p_HasNotCF(poly p1, poly p2, const ring r)
Definition: p_polys.cc:1334
void p_Shift(poly *p, int i, const ring r)
shifts components of the vector p by i
Definition: p_polys.cc:4483
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3621
#define rRing_has_Comp(r)
Definition: monomials.h:274
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:850
#define p_SetmComp
Definition: p_polys.h:233
void nlNormalize(number &x, const coeffs r)
Definition: longrat.cc:1277
poly p_mInit(const char *st, BOOLEAN &ok, const ring r)
Definition: p_polys.cc:1425
unsigned long p_GetShortExpVector(const poly p, const ring r)
Definition: p_polys.cc:4557
#define p_LmEqual(p1, p2, r)
Definition: p_polys.h:1519
poly p_PermPoly(poly p, const int *perm, const ring oldRing, const ring dst, nMapFunc nMap, const int *par_perm, int OldPar)
Definition: p_polys.cc:3926
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:4524
#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:484
poly p_DivideM(poly a, poly b, const ring r)
Definition: p_polys.cc:1501
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:1559
int p_IsPurePower(const poly p, const ring r)
return i, if head depends only on var(i)
Definition: p_polys.cc:1224
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:631
static pLDegProc pOldLDeg
Definition: p_polys.cc:3505
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:437
#define NULL
Definition: omList.c:10
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3540
long pLDeg1_Totaldegree(poly p, int *l, const ring r)
Definition: p_polys.cc:974
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of &#39;n&#39;
Definition: coeffs.h:452
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:3822
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:35
LINLINE void nlDelete(number *a, const coeffs r)
Definition: longrat.cc:2417
poly p_Last(const poly p, int &l, const ring r)
Definition: p_polys.cc:4398
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of &#39;a&#39; and &#39;b&#39;, i.e., a/b; raises an error if &#39;b&#39; is not invertible in r exceptio...
Definition: coeffs.h:615
long(* pFDegProc)(poly p, ring r)
Definition: ring.h:46
static pFDegProc pOldFDeg
Definition: p_polys.cc:3504
#define SR_INT
Definition: longrat.h:65
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:604
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:461
static poly p_LmInit(poly p, const ring r)
Definition: p_polys.h:1263
#define Sy_bit_L(x)
static FORCE_INLINE number n_ExactDiv(number a, number b, const coeffs r)
assume that there is a canonical subring in cf and we know that division is possible for these a and ...
Definition: coeffs.h:622
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:436
#define pSetCoeff0(p, n)
Definition: monomials.h:67
poly p_DiffOp(poly a, poly b, BOOLEAN multiply, const ring r)
Definition: p_polys.cc:1884
static long p_AddExp(poly p, int v, long ee, ring r)
Definition: p_polys.h:602
#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:1300
static FORCE_INLINE number n_SubringGcd(number a, number b, const coeffs r)
Definition: coeffs.h:687
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:844
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:1018
static void p_LmDelete(poly p, const ring r)
Definition: p_polys.h:707
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:3697
#define BIT_SIZEOF_LONG
Definition: auxiliary.h:124
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:456
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:4135
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:3288
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:495
static void p_ExpVectorSum(poly pr, poly p1, poly p2, const ring r)
Definition: p_polys.h:1353
static void p_SplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero, const ring r)
Definition: p_polys.cc:3634
static BOOLEAN p_ExpVectorEqual(poly p1, poly p2, const ring r1, const ring r2)
Definition: p_polys.cc:4303
static long p_DecrExp(poly p, int v, ring r)
Definition: p_polys.h:594
polyrec * poly
Definition: hilb.h:10
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:884
static BOOLEAN rField_has_Units(const ring r)
Definition: ring.h:443
poly p_TakeOutComp(poly *p, int k, const ring r)
Definition: p_polys.cc:3317
int offset
Definition: libparse.cc:1091
static Poly * h
Definition: janet.cc:978
s?
Definition: ring.h:678
int BOOLEAN
Definition: auxiliary.h:131
static poly p_Init(const ring r, omBin bin)
Definition: p_polys.h:1248
const poly b
Definition: syzextra.cc:213
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2682
int p_Var(poly m, const ring r)
Definition: p_polys.cc:4433
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:571
void p_ProjectiveUnique(poly ph, const ring r)
Definition: p_polys.cc:3008
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1302
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1025
static bool rIsRatGRing(const ring r)
Definition: ring.h:372
poly p_Power(poly p, int i, const ring r)
Definition: p_polys.cc:2100
void Werror(const char *fmt,...)
Definition: reporter.cc:199
void p_DeleteComp(poly *p, int k, const ring r)
Definition: p_polys.cc:3426
#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:939
static long p_MaxComp(poly p, ring lmRing, ring tailRing)
Definition: p_polys.h:281
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:2000