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

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