source: git/libpolys/polys/monomials/p_polys.cc @ 920c78

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