source: git/libpolys/polys/monomials/p_polys.cc @ 644b31

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