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

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