source: git/libpolys/polys/monomials/p_polys.cc @ a3f0fea

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