source: git/libpolys/polys/monomials/p_polys.cc @ 628f663

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