source: git/libpolys/polys/monomials/p_polys.cc @ 975db18

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