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

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