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

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