source: git/libpolys/polys/monomials/p_polys.cc @ 69f579

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