source: git/libpolys/polys/monomials/p_polys.cc @ 16c400

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