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

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