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

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