source: git/libpolys/polys/monomials/p_polys.cc @ 0b4a87

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