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

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