source: git/libpolys/polys/monomials/p_polys.cc @ 805d0b1

jengelh-datetimespielwiese
Last change on this file since 805d0b1 was 805d0b1, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
improved (standalone!) header for polynomials (polys/monomials/monomials.h) chg: moved the basic primitives for accessing and setting fields of a poly from p_polys.h into monomials.h NOTE: these basic functions either do NOT need a ring or are macros (or both) chg: no need to include coeffs/coeffs.h and polys/monomials/ring.h in polys/monomials/monomials.h NOTE: just forward-declare structures and typedef pointers at them. chg: minor cleanup of includes in p_polys.cc
  • Property mode set to 100644
File size: 84.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    p_polys.cc
6 *  Purpose: implementation of ring independent poly procedures?
7 *  Author:  obachman (Olaf Bachmann)
8 *  Created: 8/00
9 *******************************************************************/
10
11
12
13#include "config.h"
14#include <misc/auxiliary.h>
15
16#include <ctype.h>
17#include <omalloc/omalloc.h>
18#include <misc/options.h>
19#include <misc/intvec.h>
20
21#include <coeffs/longrat.h> // ???
22#include <coeffs/ffields.h>
23
24#define TRANSEXT_PRIVATES
25
26#include <polys/ext_fields/transext.h>
27#include <polys/ext_fields/algext.h>
28
29#include <polys/weight.h>
30#include <polys/simpleideals.h>
31
32#include "ring.h"
33#include "p_polys.h"
34
35#include <polys/templates/p_MemCmp.h>
36#include <polys/templates/p_MemAdd.h>
37#include <polys/templates/p_MemCopy.h>
38
39
40// #include <???/ideals.h>
41// #include <???/int64vec.h>
42
43#ifndef NDEBUG
44// #include <???/febase.h>
45#endif
46
47#ifdef HAVE_PLURAL
48#include "nc/nc.h"
49#include "nc/sca.h"
50#endif
51
52#include "coeffrings.h"
53#ifdef HAVE_FACTORY
54#include "clapsing.h"
55#endif
56
57/***************************************************************
58 *
59 * Completing what needs to be set for the monomial
60 *
61 ***************************************************************/
62// this is special for the syz stuff
63static int* _components = NULL;
64static long* _componentsShifted = NULL;
65static int _componentsExternal = 0;
66
67BOOLEAN pSetm_error=0;
68
69#ifndef NDEBUG
70# define MYTEST 0
71#else /* ifndef NDEBUG */
72# define MYTEST 0
73#endif /* ifndef NDEBUG */
74
75void p_Setm_General(poly p, const ring r)
76{
77  p_LmCheckPolyRing(p, r);
78  int pos=0;
79  if (r->typ!=NULL)
80  {
81    loop
82    {
83      long ord=0;
84      sro_ord* o=&(r->typ[pos]);
85      switch(o->ord_typ)
86      {
87        case ro_dp:
88        {
89          int a,e;
90          a=o->data.dp.start;
91          e=o->data.dp.end;
92          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r);
93          p->exp[o->data.dp.place]=ord;
94          break;
95        }
96        case ro_wp_neg:
97          ord=POLY_NEGWEIGHT_OFFSET;
98          // no break;
99        case ro_wp:
100        {
101          int a,e;
102          a=o->data.wp.start;
103          e=o->data.wp.end;
104          int *w=o->data.wp.weights;
105#if 1
106          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r)*w[i-a];
107#else
108          long ai;
109          int ei,wi;
110          for(int i=a;i<=e;i++)
111          {
112             ei=p_GetExp(p,i,r);
113             wi=w[i-a];
114             ai=ei*wi;
115             if (ai/ei!=wi) pSetm_error=TRUE;
116             ord+=ai;
117             if (ord<ai) pSetm_error=TRUE;
118          }
119#endif
120          p->exp[o->data.wp.place]=ord;
121          break;
122        }
123        case ro_am:
124        {
125          ord = POLY_NEGWEIGHT_OFFSET;
126          const short a=o->data.am.start;
127          const short e=o->data.am.end;
128          const int * w=o->data.am.weights;
129#if 1
130          for(short i=a; i<=e; i++, w++)
131            ord += ((*w) * p_GetExp(p,i,r));
132#else
133          long ai;
134          int ei,wi;
135          for(short i=a;i<=e;i++)
136          {
137             ei=p_GetExp(p,i,r);
138             wi=w[i-a];
139             ai=ei*wi;
140             if (ai/ei!=wi) pSetm_error=TRUE;
141             ord += ai;
142             if (ord<ai) pSetm_error=TRUE;
143          }
144#endif
145          const int c = p_GetComp(p,r);
146
147          const short len_gen= o->data.am.len_gen;
148         
149          if ((c > 0) && (c <= len_gen))
150          {
151            assume( w == o->data.am.weights_m );
152            assume( w[0] == len_gen );
153            ord += w[c];
154          }
155         
156          p->exp[o->data.am.place] = ord;
157          break;
158        }
159      case ro_wp64:
160        {
161          int64 ord=0;
162          int a,e;
163          a=o->data.wp64.start;
164          e=o->data.wp64.end;
165          int64 *w=o->data.wp64.weights64;
166          int64 ei,wi,ai;
167          for(int i=a;i<=e;i++)
168          {
169            //Print("exp %d w %d \n",p_GetExp(p,i,r),(int)w[i-a]);
170            //ord+=((int64)p_GetExp(p,i,r))*w[i-a];
171            ei=(int64)p_GetExp(p,i,r);
172            wi=w[i-a];
173            ai=ei*wi;
174            if(ei!=0 && ai/ei!=wi)
175            {
176              pSetm_error=TRUE;
177              #if SIZEOF_LONG == 4
178              Print("ai %lld, wi %lld\n",ai,wi);
179              #else
180              Print("ai %ld, wi %ld\n",ai,wi);
181              #endif
182            }
183            ord+=ai;
184            if (ord<ai)
185            {
186              pSetm_error=TRUE;
187              #if SIZEOF_LONG == 4
188              Print("ai %lld, ord %lld\n",ai,ord);
189              #else
190              Print("ai %ld, ord %ld\n",ai,ord);
191              #endif
192            }
193          }
194          int64 mask=(int64)0x7fffffff;
195          long a_0=(long)(ord&mask); //2^31
196          long a_1=(long)(ord >>31 ); /*(ord/(mask+1));*/
197
198          //Print("mask: %x,  ord: %d,  a_0: %d,  a_1: %d\n"
199          //,(int)mask,(int)ord,(int)a_0,(int)a_1);
200                    //Print("mask: %d",mask);
201
202          p->exp[o->data.wp64.place]=a_1;
203          p->exp[o->data.wp64.place+1]=a_0;
204//            if(p_Setm_error) Print("***************************\n
205//                                    ***************************\n
206//                                    **WARNING: overflow error**\n
207//                                    ***************************\n
208//                                    ***************************\n");
209          break;
210        }
211        case ro_cp:
212        {
213          int a,e;
214          a=o->data.cp.start;
215          e=o->data.cp.end;
216          int pl=o->data.cp.place;
217          for(int i=a;i<=e;i++) { p->exp[pl]=p_GetExp(p,i,r); pl++; }
218          break;
219        }
220        case ro_syzcomp:
221        {
222          int c=p_GetComp(p,r);
223          long sc = c;
224          int* Components = (_componentsExternal ? _components :
225                             o->data.syzcomp.Components);
226          long* ShiftedComponents = (_componentsExternal ? _componentsShifted:
227                                     o->data.syzcomp.ShiftedComponents);
228          if (ShiftedComponents != NULL)
229          {
230            assume(Components != NULL);
231            assume(c == 0 || Components[c] != 0);
232            sc = ShiftedComponents[Components[c]];
233            assume(c == 0 || sc != 0);
234          }
235          p->exp[o->data.syzcomp.place]=sc;
236          break;
237        }
238        case ro_syz:
239        {
240          const unsigned long c = p_GetComp(p, r);
241          const short place = o->data.syz.place;
242          const int limit = o->data.syz.limit;
243
244          if (c > limit)
245            p->exp[place] = o->data.syz.curr_index;
246          else if (c > 0)
247          {
248            assume( (1 <= c) && (c <= limit) );
249            p->exp[place]= o->data.syz.syz_index[c];
250          }
251          else
252          {
253            assume(c == 0);
254            p->exp[place]= 0;
255          }
256          break;
257        }
258        // Prefix for Induced Schreyer ordering
259        case ro_isTemp: // Do nothing?? (to be removed into suffix later on...?)
260        {
261          assume(p != NULL);
262
263#ifndef NDEBUG
264#if MYTEST
265          Print("p_Setm_General: ro_isTemp ord: pos: %d, p: ", pos);  p_DebugPrint(p, r, r, 1);
266#endif
267#endif
268          int c = p_GetComp(p, r);
269
270          assume( c >= 0 );
271
272          // Let's simulate case ro_syz above....
273          // Should accumulate (by Suffix) and be a level indicator
274          const int* const pVarOffset = o->data.isTemp.pVarOffset;
275
276          assume( pVarOffset != NULL );
277
278          // TODO: Can this be done in the suffix???
279          for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
280          {
281            const int vo = pVarOffset[i];
282            if( vo != -1) // TODO: optimize: can be done once!
283            {
284              // Hans! Please don't break it again! p_SetExp(p, ..., r, vo) is correct:
285              p_SetExp(p, p_GetExp(p, i, r), r, vo); // copy put them verbatim
286              // Hans! Please don't break it again! p_GetExp(p, r, vo) is correct:
287              assume( p_GetExp(p, r, vo) == p_GetExp(p, i, r) ); // copy put them verbatim
288            }
289          }
290#ifndef NDEBUG
291          for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
292          {
293            const int vo = pVarOffset[i];
294            if( vo != -1) // TODO: optimize: can be done once!
295            {
296              // Hans! Please don't break it again! p_GetExp(p, r, vo) is correct:
297              assume( p_GetExp(p, r, vo) == p_GetExp(p, i, r) ); // copy put them verbatim
298            }
299          }
300#if MYTEST
301//          if( p->exp[o->data.isTemp.start] > 0 )
302            PrintS("after Values: "); p_DebugPrint(p, r, r, 1);
303#endif
304#endif
305          break;
306        }
307
308        // Suffix for Induced Schreyer ordering
309        case ro_is:
310        {
311#ifndef NDEBUG
312#if MYTEST
313          Print("p_Setm_General: ro_is ord: pos: %d, p: ", pos);  p_DebugPrint(p, r, r, 1);
314#endif
315#endif
316
317          assume(p != NULL);
318
319          int c = p_GetComp(p, r);
320
321          assume( c >= 0 );
322          const ideal F = o->data.is.F;
323          const int limit = o->data.is.limit;
324          assume( limit >= 0 );
325          const int start = o->data.is.start;
326
327          if( F != NULL && c > limit )
328          {
329#ifndef NDEBUG
330#if MYTEST
331            Print("p_Setm_General: ro_is : in rSetm: pos: %d, c: %d >  limit: %d\n", c, pos, limit); // p_DebugPrint(p, r, r, 1);
332            PrintS("preComputed Values: ");
333            p_DebugPrint(p, r, r, 1);
334#endif
335#endif
336//          if( c > limit ) // BUG???
337            p->exp[start] = 1;
338//          else
339//            p->exp[start] = 0;
340
341
342            c -= limit;
343            assume( c > 0 );
344            c--;
345
346            if( c >= IDELEMS(F) )
347              break;
348
349            assume( c < IDELEMS(F) ); // What about others???
350           
351            const poly pp = F->m[c]; // get reference monomial!!!
352
353            if(pp == NULL)
354              break;
355           
356            assume(pp != NULL);
357
358#ifndef NDEBUG
359#if MYTEST
360            Print("Respective F[c - %d: %d] pp: ", limit, c);
361            p_DebugPrint(pp, r, r, 1);
362#endif
363#endif
364
365            const int end = o->data.is.end;
366            assume(start <= end);
367
368
369//          const int st = o->data.isTemp.start;
370
371#ifndef NDEBUG
372            Print("p_Setm_General: is(-Temp-) :: c: %d, limit: %d, [st:%d] ===>>> %ld\n", c, limit, start, p->exp[start]);
373#endif
374
375            // p_ExpVectorAdd(p, pp, r);
376
377            for( int i = start; i <= end; i++) // v[0] may be here...
378              p->exp[i] += pp->exp[i]; // !!!!!!!! ADD corresponding LT(F)
379
380            // p_MemAddAdjust(p, ri);
381            if (r->NegWeightL_Offset != NULL)
382            {
383              for (int i=r->NegWeightL_Size-1; i>=0; i--)
384              {
385                const int _i = r->NegWeightL_Offset[i];
386                if( start <= _i && _i <= end )
387                  p->exp[_i] -= POLY_NEGWEIGHT_OFFSET;
388              }
389            }
390
391
392#ifndef NDEBUG
393            const int* const pVarOffset = o->data.is.pVarOffset;
394
395            assume( pVarOffset != NULL );
396
397            for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
398            {
399              const int vo = pVarOffset[i];
400              if( vo != -1) // TODO: optimize: can be done once!
401                // Hans! Please don't break it again! p_GetExp(p/pp, r, vo) is correct:
402                assume( p_GetExp(p, r, vo) == (p_GetExp(p, i, r) + p_GetExp(pp, r, vo)) );
403            }
404            // TODO: how to check this for computed values???
405#if MYTEST
406            PrintS("Computed Values: "); p_DebugPrint(p, r, r, 1);
407#endif
408#endif
409          } else
410          {
411            p->exp[start] = 0; //!!!!????? where?????
412
413            const int* const pVarOffset = o->data.is.pVarOffset;
414
415            // What about v[0] - component: it will be added later by
416            // suffix!!!
417            // TODO: Test it!
418            const int vo = pVarOffset[0];
419            if( vo != -1 )
420              p->exp[vo] = c; // initial component v[0]!
421
422#ifndef NDEBUG
423#if MYTEST
424            Print("ELSE p_Setm_General: ro_is :: c: %d <= limit: %d, vo: %d, exp: %d\n", c, limit, vo, p->exp[vo]);
425            p_DebugPrint(p, r, r, 1);
426#endif
427#endif
428          }
429
430          break;
431        }
432        default:
433          dReportError("wrong ord in rSetm:%d\n",o->ord_typ);
434          return;
435      }
436      pos++;
437      if (pos == r->OrdSize) return;
438    }
439  }
440}
441
442void p_Setm_Syz(poly p, ring r, int* Components, long* ShiftedComponents)
443{
444  _components = Components;
445  _componentsShifted = ShiftedComponents;
446  _componentsExternal = 1;
447  p_Setm_General(p, r);
448  _componentsExternal = 0;
449}
450
451// dummy for lp, ls, etc
452void p_Setm_Dummy(poly p, const ring r)
453{
454  p_LmCheckPolyRing(p, r);
455}
456
457// for dp, Dp, ds, etc
458void p_Setm_TotalDegree(poly p, const ring r)
459{
460  p_LmCheckPolyRing(p, r);
461  p->exp[r->pOrdIndex] = p_Totaldegree(p, r);
462}
463
464// for wp, Wp, ws, etc
465void p_Setm_WFirstTotalDegree(poly p, const ring r)
466{
467  p_LmCheckPolyRing(p, r);
468  p->exp[r->pOrdIndex] = p_WFirstTotalDegree(p, r);
469}
470
471p_SetmProc p_GetSetmProc(ring r)
472{
473  // covers lp, rp, ls,
474  if (r->typ == NULL) return p_Setm_Dummy;
475
476  if (r->OrdSize == 1)
477  {
478    if (r->typ[0].ord_typ == ro_dp &&
479        r->typ[0].data.dp.start == 1 &&
480        r->typ[0].data.dp.end == r->N &&
481        r->typ[0].data.dp.place == r->pOrdIndex)
482      return p_Setm_TotalDegree;
483    if (r->typ[0].ord_typ == ro_wp &&
484        r->typ[0].data.wp.start == 1 &&
485        r->typ[0].data.wp.end == r->N &&
486        r->typ[0].data.wp.place == r->pOrdIndex &&
487        r->typ[0].data.wp.weights == r->firstwv)
488      return p_Setm_WFirstTotalDegree;
489  }
490  return p_Setm_General;
491}
492
493
494/* -------------------------------------------------------------------*/
495/* several possibilities for pFDeg: the degree of the head term       */
496
497/* comptible with ordering */
498long p_Deg(poly a, const ring r)
499{
500  p_LmCheckPolyRing(a, r);
501//  assume(p_GetOrder(a, r) == p_WTotaldegree(a, r)); // WRONG assume!
502  return p_GetOrder(a, r);
503}
504
505// p_WTotalDegree for weighted orderings
506// whose first block covers all variables
507long p_WFirstTotalDegree(poly p, const ring r)
508{
509  int i;
510  long sum = 0;
511
512  for (i=1; i<= r->firstBlockEnds; i++)
513  {
514    sum += p_GetExp(p, i, r)*r->firstwv[i-1];
515  }
516  return sum;
517}
518
519/*2
520* compute the degree of the leading monomial of p
521* with respect to weigths from the ordering
522* the ordering is not compatible with degree so do not use p->Order
523*/
524long p_WTotaldegree(poly p, const ring r)
525{
526  p_LmCheckPolyRing(p, r);
527  int i, k;
528  int pIS = 0;
529  long j =0;
530
531  // iterate through each block:
532  for (i=0;r->order[i]!=0;i++)
533  {
534    int b0=r->block0[i];
535    int b1=r->block1[i];
536    switch(r->order[i])
537    {
538      case ringorder_M:
539        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
540        { // in jedem block:
541          j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/]*r->OrdSgn;
542        }
543        break;
544      case ringorder_wp:
545      case ringorder_ws:
546      case ringorder_Wp:
547      case ringorder_Ws:
548        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
549        { // in jedem block:
550          j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/];
551        }
552        break;
553      case ringorder_lp:
554      case ringorder_ls:
555      case ringorder_rs:
556      case ringorder_dp:
557      case ringorder_ds:
558      case ringorder_Dp:
559      case ringorder_Ds:
560      case ringorder_rp:
561        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
562        {
563          j+= p_GetExp(p,k,r);
564        }
565        break;
566      case ringorder_a64:
567        {
568          int64* w=(int64*)r->wvhdl[i];
569          for (k=0;k<=(b1 /*r->block1[i]*/ - b0 /*r->block0[i]*/);k++)
570          {
571            //there should be added a line which checks if w[k]>2^31
572            j+= p_GetExp(p,k+1, r)*(long)w[k];
573          }
574          //break;
575          return j;
576        }
577      case ringorder_c:
578      case ringorder_C:
579      case ringorder_S:
580      case ringorder_s:
581      case ringorder_aa:
582      case ringorder_IS:
583        break;
584      case ringorder_a:
585        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
586        { // only one line
587          j+= p_GetExp(p,k, r)*r->wvhdl[i][ k- b0 /*r->block0[i]*/];
588        }
589        //break;
590        return j;
591
592#ifndef NDEBUG
593      default:
594        Print("missing order %d in p_WTotaldegree\n",r->order[i]);
595        break;
596#endif
597    }
598  }
599  return  j;
600}
601
602long p_DegW(poly p, const short *w, const ring R)
603{
604  long r=~0L;
605
606  while (p!=NULL)
607  {
608    long t=totaldegreeWecart_IV(p,R,w);
609    if (t>r) r=t;
610    pIter(p);
611  }
612  return r;
613}
614
615int p_Weight(int i, const ring r)
616{
617  if ((r->firstwv==NULL) || (i>r->firstBlockEnds))
618  {
619    return 1;
620  }
621  return r->firstwv[i-1];
622}
623
624long p_WDegree(poly p, const ring r)
625{
626  if (r->firstwv==NULL) return p_Totaldegree(p, r);
627  p_LmCheckPolyRing(p, r);
628  int i;
629  long j =0;
630
631  for(i=1;i<=r->firstBlockEnds;i++)
632    j+=p_GetExp(p, i, r)*r->firstwv[i-1];
633
634  for (;i<=r->N;i++)
635    j+=p_GetExp(p,i, r)*p_Weight(i, r);
636
637  return j;
638}
639
640
641/* ---------------------------------------------------------------------*/
642/* several possibilities for pLDeg: the maximal degree of a monomial in p*/
643/*  compute in l also the pLength of p                                   */
644
645/*2
646* compute the length of a polynomial (in l)
647* and the degree of the monomial with maximal degree: the last one
648*/
649long pLDeg0(poly p,int *l, const ring r)
650{
651  p_CheckPolyRing(p, r);
652  long k= p_GetComp(p, r);
653  int ll=1;
654
655  if (k > 0)
656  {
657    while ((pNext(p)!=NULL) && (p_GetComp(pNext(p), r)==k))
658    {
659      pIter(p);
660      ll++;
661    }
662  }
663  else
664  {
665     while (pNext(p)!=NULL)
666     {
667       pIter(p);
668       ll++;
669     }
670  }
671  *l=ll;
672  return r->pFDeg(p, r);
673}
674
675/*2
676* compute the length of a polynomial (in l)
677* and the degree of the monomial with maximal degree: the last one
678* but search in all components before syzcomp
679*/
680long pLDeg0c(poly p,int *l, const ring r)
681{
682  assume(p!=NULL);
683#ifdef PDEBUG
684  _p_Test(p,r,PDEBUG);
685#endif
686  p_CheckPolyRing(p, r);
687  long o;
688  int ll=1;
689
690  if (! rIsSyzIndexRing(r))
691  {
692    while (pNext(p) != NULL)
693    {
694      pIter(p);
695      ll++;
696    }
697    o = r->pFDeg(p, r);
698  }
699  else
700  {
701    int curr_limit = rGetCurrSyzLimit(r);
702    poly pp = p;
703    while ((p=pNext(p))!=NULL)
704    {
705      if (p_GetComp(p, r)<=curr_limit/*syzComp*/)
706        ll++;
707      else break;
708      pp = p;
709    }
710#ifdef PDEBUG
711    _p_Test(pp,r,PDEBUG);
712#endif
713    o = r->pFDeg(pp, r);
714  }
715  *l=ll;
716  return o;
717}
718
719/*2
720* compute the length of a polynomial (in l)
721* and the degree of the monomial with maximal degree: the first one
722* this works for the polynomial case with degree orderings
723* (both c,dp and dp,c)
724*/
725long pLDegb(poly p,int *l, const ring r)
726{
727  p_CheckPolyRing(p, r);
728  long k= p_GetComp(p, r);
729  long o = r->pFDeg(p, r);
730  int ll=1;
731
732  if (k != 0)
733  {
734    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
735    {
736      ll++;
737    }
738  }
739  else
740  {
741    while ((p=pNext(p)) !=NULL)
742    {
743      ll++;
744    }
745  }
746  *l=ll;
747  return o;
748}
749
750/*2
751* compute the length of a polynomial (in l)
752* and the degree of the monomial with maximal degree:
753* this is NOT the last one, we have to look for it
754*/
755long pLDeg1(poly p,int *l, const ring r)
756{
757  p_CheckPolyRing(p, r);
758  long k= p_GetComp(p, r);
759  int ll=1;
760  long  t,max;
761
762  max=r->pFDeg(p, r);
763  if (k > 0)
764  {
765    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
766    {
767      t=r->pFDeg(p, r);
768      if (t>max) max=t;
769      ll++;
770    }
771  }
772  else
773  {
774    while ((p=pNext(p))!=NULL)
775    {
776      t=r->pFDeg(p, r);
777      if (t>max) max=t;
778      ll++;
779    }
780  }
781  *l=ll;
782  return max;
783}
784
785/*2
786* compute the length of a polynomial (in l)
787* and the degree of the monomial with maximal degree:
788* this is NOT the last one, we have to look for it
789* in all components
790*/
791long pLDeg1c(poly p,int *l, const ring r)
792{
793  p_CheckPolyRing(p, r);
794  int ll=1;
795  long  t,max;
796
797  max=r->pFDeg(p, r);
798  if (rIsSyzIndexRing(r))
799  {
800    long limit = rGetCurrSyzLimit(r);
801    while ((p=pNext(p))!=NULL)
802    {
803      if (p_GetComp(p, r)<=limit)
804      {
805        if ((t=r->pFDeg(p, r))>max) max=t;
806        ll++;
807      }
808      else break;
809    }
810  }
811  else
812  {
813    while ((p=pNext(p))!=NULL)
814    {
815      if ((t=r->pFDeg(p, r))>max) max=t;
816      ll++;
817    }
818  }
819  *l=ll;
820  return max;
821}
822
823// like pLDeg1, only pFDeg == pDeg
824long pLDeg1_Deg(poly p,int *l, const ring r)
825{
826  assume(r->pFDeg == p_Deg);
827  p_CheckPolyRing(p, r);
828  long k= p_GetComp(p, r);
829  int ll=1;
830  long  t,max;
831
832  max=p_GetOrder(p, r);
833  if (k > 0)
834  {
835    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
836    {
837      t=p_GetOrder(p, r);
838      if (t>max) max=t;
839      ll++;
840    }
841  }
842  else
843  {
844    while ((p=pNext(p))!=NULL)
845    {
846      t=p_GetOrder(p, r);
847      if (t>max) max=t;
848      ll++;
849    }
850  }
851  *l=ll;
852  return max;
853}
854
855long pLDeg1c_Deg(poly p,int *l, const ring r)
856{
857  assume(r->pFDeg == p_Deg);
858  p_CheckPolyRing(p, r);
859  int ll=1;
860  long  t,max;
861
862  max=p_GetOrder(p, r);
863  if (rIsSyzIndexRing(r))
864  {
865    long limit = rGetCurrSyzLimit(r);
866    while ((p=pNext(p))!=NULL)
867    {
868      if (p_GetComp(p, r)<=limit)
869      {
870        if ((t=p_GetOrder(p, r))>max) max=t;
871        ll++;
872      }
873      else break;
874    }
875  }
876  else
877  {
878    while ((p=pNext(p))!=NULL)
879    {
880      if ((t=p_GetOrder(p, r))>max) max=t;
881      ll++;
882    }
883  }
884  *l=ll;
885  return max;
886}
887
888// like pLDeg1, only pFDeg == pTotoalDegree
889long pLDeg1_Totaldegree(poly p,int *l, const ring r)
890{
891  p_CheckPolyRing(p, r);
892  long k= p_GetComp(p, r);
893  int ll=1;
894  long  t,max;
895
896  max=p_Totaldegree(p, r);
897  if (k > 0)
898  {
899    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
900    {
901      t=p_Totaldegree(p, r);
902      if (t>max) max=t;
903      ll++;
904    }
905  }
906  else
907  {
908    while ((p=pNext(p))!=NULL)
909    {
910      t=p_Totaldegree(p, r);
911      if (t>max) max=t;
912      ll++;
913    }
914  }
915  *l=ll;
916  return max;
917}
918
919long pLDeg1c_Totaldegree(poly p,int *l, const ring r)
920{
921  p_CheckPolyRing(p, r);
922  int ll=1;
923  long  t,max;
924
925  max=p_Totaldegree(p, r);
926  if (rIsSyzIndexRing(r))
927  {
928    long limit = rGetCurrSyzLimit(r);
929    while ((p=pNext(p))!=NULL)
930    {
931      if (p_GetComp(p, r)<=limit)
932      {
933        if ((t=p_Totaldegree(p, r))>max) max=t;
934        ll++;
935      }
936      else break;
937    }
938  }
939  else
940  {
941    while ((p=pNext(p))!=NULL)
942    {
943      if ((t=p_Totaldegree(p, r))>max) max=t;
944      ll++;
945    }
946  }
947  *l=ll;
948  return max;
949}
950
951// like pLDeg1, only pFDeg == p_WFirstTotalDegree
952long pLDeg1_WFirstTotalDegree(poly p,int *l, const ring r)
953{
954  p_CheckPolyRing(p, r);
955  long k= p_GetComp(p, r);
956  int ll=1;
957  long  t,max;
958
959  max=p_WFirstTotalDegree(p, r);
960  if (k > 0)
961  {
962    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
963    {
964      t=p_WFirstTotalDegree(p, r);
965      if (t>max) max=t;
966      ll++;
967    }
968  }
969  else
970  {
971    while ((p=pNext(p))!=NULL)
972    {
973      t=p_WFirstTotalDegree(p, r);
974      if (t>max) max=t;
975      ll++;
976    }
977  }
978  *l=ll;
979  return max;
980}
981
982long pLDeg1c_WFirstTotalDegree(poly p,int *l, const ring r)
983{
984  p_CheckPolyRing(p, r);
985  int ll=1;
986  long  t,max;
987
988  max=p_WFirstTotalDegree(p, r);
989  if (rIsSyzIndexRing(r))
990  {
991    long limit = rGetCurrSyzLimit(r);
992    while ((p=pNext(p))!=NULL)
993    {
994      if (p_GetComp(p, r)<=limit)
995      {
996        if ((t=p_Totaldegree(p, r))>max) max=t;
997        ll++;
998      }
999      else break;
1000    }
1001  }
1002  else
1003  {
1004    while ((p=pNext(p))!=NULL)
1005    {
1006      if ((t=p_Totaldegree(p, r))>max) max=t;
1007      ll++;
1008    }
1009  }
1010  *l=ll;
1011  return max;
1012}
1013
1014/***************************************************************
1015 *
1016 * Maximal Exponent business
1017 *
1018 ***************************************************************/
1019
1020static inline unsigned long
1021p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r,
1022              unsigned long number_of_exp)
1023{
1024  const unsigned long bitmask = r->bitmask;
1025  unsigned long ml1 = l1 & bitmask;
1026  unsigned long ml2 = l2 & bitmask;
1027  unsigned long max = (ml1 > ml2 ? ml1 : ml2);
1028  unsigned long j = number_of_exp - 1;
1029
1030  if (j > 0)
1031  {
1032    unsigned long mask = bitmask << r->BitsPerExp;
1033    while (1)
1034    {
1035      ml1 = l1 & mask;
1036      ml2 = l2 & mask;
1037      max |= ((ml1 > ml2 ? ml1 : ml2) & mask);
1038      j--;
1039      if (j == 0) break;
1040      mask = mask << r->BitsPerExp;
1041    }
1042  }
1043  return max;
1044}
1045
1046static inline unsigned long
1047p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r)
1048{
1049  return p_GetMaxExpL2(l1, l2, r, r->ExpPerLong);
1050}
1051
1052poly p_GetMaxExpP(poly p, const ring r)
1053{
1054  p_CheckPolyRing(p, r);
1055  if (p == NULL) return p_Init(r);
1056  poly max = p_LmInit(p, r);
1057  pIter(p);
1058  if (p == NULL) return max;
1059  int i, offset;
1060  unsigned long l_p, l_max;
1061  unsigned long divmask = r->divmask;
1062
1063  do
1064  {
1065    offset = r->VarL_Offset[0];
1066    l_p = p->exp[offset];
1067    l_max = max->exp[offset];
1068    // do the divisibility trick to find out whether l has an exponent
1069    if (l_p > l_max ||
1070        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1071      max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
1072
1073    for (i=1; i<r->VarL_Size; i++)
1074    {
1075      offset = r->VarL_Offset[i];
1076      l_p = p->exp[offset];
1077      l_max = max->exp[offset];
1078      // do the divisibility trick to find out whether l has an exponent
1079      if (l_p > l_max ||
1080          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1081        max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
1082    }
1083    pIter(p);
1084  }
1085  while (p != NULL);
1086  return max;
1087}
1088
1089unsigned long p_GetMaxExpL(poly p, const ring r, unsigned long l_max)
1090{
1091  unsigned long l_p, divmask = r->divmask;
1092  int i;
1093
1094  while (p != NULL)
1095  {
1096    l_p = p->exp[r->VarL_Offset[0]];
1097    if (l_p > l_max ||
1098        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1099      l_max = p_GetMaxExpL2(l_max, l_p, r);
1100    for (i=1; i<r->VarL_Size; i++)
1101    {
1102      l_p = p->exp[r->VarL_Offset[i]];
1103      // do the divisibility trick to find out whether l has an exponent
1104      if (l_p > l_max ||
1105          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
1106        l_max = p_GetMaxExpL2(l_max, l_p, r);
1107    }
1108    pIter(p);
1109  }
1110  return l_max;
1111}
1112
1113
1114
1115
1116/***************************************************************
1117 *
1118 * Misc things
1119 *
1120 ***************************************************************/
1121// returns TRUE, if all monoms have the same component
1122BOOLEAN p_OneComp(poly p, const ring r)
1123{
1124  if(p!=NULL)
1125  {
1126    long i = p_GetComp(p, r);
1127    while (pNext(p)!=NULL)
1128    {
1129      pIter(p);
1130      if(i != p_GetComp(p, r)) return FALSE;
1131    }
1132  }
1133  return TRUE;
1134}
1135
1136/*2
1137*test if a monomial /head term is a pure power
1138*/
1139int p_IsPurePower(const poly p, const ring r)
1140{
1141  int i,k=0;
1142
1143  for (i=r->N;i;i--)
1144  {
1145    if (p_GetExp(p,i, r)!=0)
1146    {
1147      if(k!=0) return 0;
1148      k=i;
1149    }
1150  }
1151  return k;
1152}
1153
1154/*2
1155*test if a polynomial is univariate
1156* return -1 for constant,
1157* 0 for not univariate,s
1158* i if dep. on var(i)
1159*/
1160int p_IsUnivariate(poly p, const ring r)
1161{
1162  int i,k=-1;
1163
1164  while (p!=NULL)
1165  {
1166    for (i=r->N;i;i--)
1167    {
1168      if (p_GetExp(p,i, r)!=0)
1169      {
1170        if((k!=-1)&&(k!=i)) return 0;
1171        k=i;
1172      }
1173    }
1174    pIter(p);
1175  }
1176  return k;
1177}
1178
1179// set entry e[i] to 1 if var(i) occurs in p, ignore var(j) if e[j]>0
1180int  p_GetVariables(poly p, int * e, const ring r)
1181{
1182  int i;
1183  int n=0;
1184  while(p!=NULL)
1185  {
1186    n=0;
1187    for(i=r->N; i>0; i--)
1188    {
1189      if(e[i]==0)
1190      {
1191        if (p_GetExp(p,i,r)>0)
1192        {
1193          e[i]=1;
1194          n++;
1195        }
1196      }
1197      else
1198        n++;
1199    }
1200    if (n==r->N) break;
1201    pIter(p);
1202  }
1203  return n;
1204}
1205
1206
1207/*2
1208* returns a polynomial representing the integer i
1209*/
1210poly p_ISet(long i, const ring r)
1211{
1212  poly rc = NULL;
1213  if (i!=0)
1214  {
1215    rc = p_Init(r);
1216    pSetCoeff0(rc,n_Init(i,r->cf));
1217    if (n_IsZero(pGetCoeff(rc),r->cf))
1218      p_LmDelete(&rc,r);
1219  }
1220  return rc;
1221}
1222
1223/*2
1224* an optimized version of p_ISet for the special case 1
1225*/
1226poly p_One(const ring r)
1227{
1228  poly rc = p_Init(r);
1229  pSetCoeff0(rc,n_Init(1,r->cf));
1230  return rc;
1231}
1232
1233void p_Split(poly p, poly *h)
1234{
1235  *h=pNext(p);
1236  pNext(p)=NULL;
1237}
1238
1239/*2
1240* pair has no common factor ? or is no polynomial
1241*/
1242BOOLEAN p_HasNotCF(poly p1, poly p2, const ring r)
1243{
1244
1245  if (p_GetComp(p1,r) > 0 || p_GetComp(p2,r) > 0)
1246    return FALSE;
1247  int i = rVar(r);
1248  loop
1249  {
1250    if ((p_GetExp(p1, i, r) > 0) && (p_GetExp(p2, i, r) > 0))
1251      return FALSE;
1252    i--;
1253    if (i == 0)
1254      return TRUE;
1255  }
1256}
1257
1258/*2
1259* convert monomial given as string to poly, e.g. 1x3y5z
1260*/
1261const char * p_Read(const char *st, poly &rc, const ring r)
1262{
1263  if (r==NULL) { rc=NULL;return st;}
1264  int i,j;
1265  rc = p_Init(r);
1266  const char *s = r->cf->cfRead(st,&(rc->coef),r->cf);
1267  if (s==st)
1268  /* i.e. it does not start with a coeff: test if it is a ringvar*/
1269  {
1270    j = r_IsRingVar(s,r);
1271    if (j >= 0)
1272    {
1273      p_IncrExp(rc,1+j,r);
1274      while (*s!='\0') s++;
1275      goto done;
1276    }
1277  }
1278  while (*s!='\0')
1279  {
1280    char ss[2];
1281    ss[0] = *s++;
1282    ss[1] = '\0';
1283    j = r_IsRingVar(ss,r);
1284    if (j >= 0)
1285    {
1286      const char *s_save=s;
1287      s = eati(s,&i);
1288      if (((unsigned long)i) >  r->bitmask)
1289      {
1290        // exponent to large: it is not a monomial
1291        p_LmDelete(&rc,r);
1292        return s_save;
1293      }
1294      p_AddExp(rc,1+j, (long)i, r);
1295    }
1296    else
1297    {
1298      // 1st char of is not a varname
1299      // We return the parsed polynomial nevertheless. This is needed when
1300      // we are parsing coefficients in a rational function field.
1301      s--;
1302      return s;
1303    }
1304  }
1305done:
1306  if (n_IsZero(pGetCoeff(rc),r->cf)) p_LmDelete(&rc,r);
1307  else
1308  {
1309#ifdef HAVE_PLURAL
1310    // in super-commutative ring
1311    // squares of anti-commutative variables are zeroes!
1312    if(rIsSCA(r))
1313    {
1314      const unsigned int iFirstAltVar = scaFirstAltVar(r);
1315      const unsigned int iLastAltVar  = scaLastAltVar(r);
1316
1317      assume(rc != NULL);
1318
1319      for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
1320        if( p_GetExp(rc, k, r) > 1 )
1321        {
1322          p_LmDelete(&rc, r);
1323          goto finish;
1324        }
1325    }
1326#endif
1327
1328    p_Setm(rc,r);
1329  }
1330finish:
1331  return s;
1332}
1333poly p_mInit(const char *st, BOOLEAN &ok, const ring r)
1334{
1335  poly p;
1336  const char *s=p_Read(st,p,r);
1337  if (*s!='\0')
1338  {
1339    if ((s!=st)&&isdigit(st[0]))
1340    {
1341      errorreported=TRUE;
1342    }
1343    ok=FALSE;
1344    p_Delete(&p,r);
1345    return NULL;
1346  }
1347  #ifdef PDEBUG
1348  _p_Test(p,r,PDEBUG);
1349  #endif
1350  ok=!errorreported;
1351  return p;
1352}
1353
1354/*2
1355* returns a polynomial representing the number n
1356* destroys n
1357*/
1358poly p_NSet(number n, const ring r)
1359{
1360  if (n_IsZero(n,r->cf))
1361  {
1362    n_Delete(&n, r->cf);
1363    return NULL;
1364  }
1365  else
1366  {
1367    poly rc = p_Init(r);
1368    pSetCoeff0(rc,n);
1369    return rc;
1370  }
1371}
1372/*2
1373* assumes that LM(a) = LM(b)*m, for some monomial m,
1374* returns the multiplicant m,
1375* leaves a and b unmodified
1376*/
1377poly p_Divide(poly a, poly b, const ring r)
1378{
1379  assume((p_GetComp(a,r)==p_GetComp(b,r)) || (p_GetComp(b,r)==0));
1380  int i;
1381  poly result = p_Init(r);
1382
1383  for(i=(int)r->N; i; i--)
1384    p_SetExp(result,i, p_GetExp(a,i,r)- p_GetExp(b,i,r),r);
1385  p_SetComp(result, p_GetComp(a,r) - p_GetComp(b,r),r);
1386  p_Setm(result,r);
1387  return result;
1388}
1389
1390poly p_Div_nn(poly p, const number n, const ring r)
1391{
1392  pAssume(!n_IsZero(n,r->cf));
1393  p_Test(p, r);
1394
1395  poly q = p;
1396  while (p != NULL)
1397  {
1398    number nc = pGetCoeff(p);
1399    pSetCoeff0(p, n_Div(nc, n, r->cf));
1400    n_Delete(&nc, r->cf);
1401    pIter(p);
1402  }
1403  p_Test(q, r);
1404  return q;
1405}
1406
1407/*2
1408* divides a by the monomial b, ignores monomials which are not divisible
1409* assumes that b is not NULL, destroyes b
1410*/
1411poly p_DivideM(poly a, poly b, const ring r)
1412{
1413  if (a==NULL) { p_Delete(&b,r); return NULL; }
1414  poly result=a;
1415  poly prev=NULL;
1416  int i;
1417#ifdef HAVE_RINGS
1418  number inv=pGetCoeff(b);
1419#else
1420  number inv=n_Invers(pGetCoeff(b),r->cf);
1421#endif
1422
1423  while (a!=NULL)
1424  {
1425    if (p_DivisibleBy(b,a,r))
1426    {
1427      for(i=(int)r->N; i; i--)
1428         p_SubExp(a,i, p_GetExp(b,i,r),r);
1429      p_SubComp(a, p_GetComp(b,r),r);
1430      p_Setm(a,r);
1431      prev=a;
1432      pIter(a);
1433    }
1434    else
1435    {
1436      if (prev==NULL)
1437      {
1438        p_LmDelete(&result,r);
1439        a=result;
1440      }
1441      else
1442      {
1443        p_LmDelete(&pNext(prev),r);
1444        a=pNext(prev);
1445      }
1446    }
1447  }
1448#ifdef HAVE_RINGS
1449  if (n_IsUnit(inv,r->cf))
1450  {
1451    inv = n_Invers(inv,r->cf);
1452    p_Mult_nn(result,inv,r);
1453    n_Delete(&inv, r->cf);
1454  }
1455  else
1456  {
1457    p_Div_nn(result,inv,r);
1458  }
1459#else
1460  p_Mult_nn(result,inv,r);
1461  n_Delete(&inv, r->cf);
1462#endif
1463  p_Delete(&b, r);
1464  return result;
1465}
1466
1467#ifdef HAVE_RINGS
1468/* TRUE iff LT(f) | LT(g) */
1469BOOLEAN p_DivisibleByRingCase(poly f, poly g, const ring r)
1470{
1471  int exponent;
1472  for(int i = (int)rVar(r); i>0; i--)
1473  {
1474    exponent = p_GetExp(g, i, r) - p_GetExp(f, i, r);
1475    if (exponent < 0) return FALSE;
1476  }
1477  return n_DivBy(pGetCoeff(g), pGetCoeff(f), r->cf);
1478}
1479#endif
1480
1481/*2
1482* returns the LCM of the head terms of a and b in *m
1483*/
1484void p_Lcm(poly a, poly b, poly m, const ring r)
1485{
1486  int i;
1487  for (i=rVar(r); i; i--)
1488  {
1489    p_SetExp(m,i, si_max( p_GetExp(a,i,r), p_GetExp(b,i,r)),r);
1490  }
1491  p_SetComp(m, si_max(p_GetComp(a,r), p_GetComp(b,r)),r);
1492  /* Don't do a pSetm here, otherwise hres/lres chockes */
1493}
1494
1495/* assumes that p and divisor are univariate polynomials in r,
1496   mentioning the same variable;
1497   assumes divisor != NULL;
1498   p may be NULL;
1499   assumes a global monomial ordering in r;
1500   performs polynomial division of p by divisor:
1501     - afterwards p contains the remainder of the division, i.e.,
1502       p_before = result * divisor + p_afterwards;
1503     - if needResult == TRUE, then the method computes and returns 'result',
1504       otherwise NULL is returned (This parametrization can be used when
1505       one is only interested in the remainder of the division. In this
1506       case, the method will be slightly faster.)
1507   leaves divisor unmodified */
1508poly p_PolyDiv(poly &p, poly divisor, BOOLEAN needResult, ring r)
1509{
1510  assume(divisor != NULL);
1511  if (p == NULL) return NULL;
1512
1513  poly result = NULL;
1514  number divisorLC = p_GetCoeff(divisor, r);
1515  int divisorLE = p_GetExp(divisor, 1, r);
1516  while ((p != NULL) && (p_Deg(p, r) >= p_Deg(divisor, r)))
1517  {
1518    /* determine t = LT(p) / LT(divisor) */
1519    poly t = p_ISet(1, r);
1520    number c = n_Div(p_GetCoeff(p, r), divisorLC, r->cf);
1521    p_SetCoeff(t, c, r);
1522    int e = p_GetExp(p, 1, r) - divisorLE;
1523    p_SetExp(t, 1, e, r);
1524    p_Setm(t, r);
1525    if (needResult) result = p_Add_q(result, p_Copy(t, r), r);
1526    p = p_Add_q(p, p_Neg(p_Mult_q(t, p_Copy(divisor, r), r), r), r);
1527  }
1528  return result;
1529}
1530
1531/* returns NULL if p == NULL, otherwise makes p monic by dividing
1532   by its leading coefficient (only done if this is not already 1);
1533   this assumes that we are over a ground field so that division
1534   is well-defined;
1535   modifies p */
1536void p_Monic(poly p, const ring r)
1537{
1538  if (p == NULL) return;
1539  number n = n_Init(1, r->cf);
1540  if (p->next==NULL) { p_SetCoeff(p,n,r); return; }
1541  poly pp = p;
1542  number lc = p_GetCoeff(p, r);
1543  if (n_IsOne(lc, r->cf)) return;
1544  number lcInverse = n_Invers(lc, r->cf);
1545  p_SetCoeff(p, n, r);   // destroys old leading coefficient!
1546  pIter(p);
1547  while (p != NULL)
1548  {
1549    number n = n_Mult(p_GetCoeff(p, r), lcInverse, r->cf);
1550    n_Normalize(n,r->cf);
1551    p_SetCoeff(p, n, r);   // destroys old leading coefficient!
1552    pIter(p);
1553  }
1554  n_Delete(&lcInverse, r->cf);
1555  p = pp;
1556}
1557
1558/* see p_Gcd;
1559   additional assumption: deg(p) >= deg(q);
1560   must destroy p and q (unless one of them is returned) */
1561poly p_GcdHelper(poly &p, poly &q, ring r)
1562{
1563  if (q == NULL)
1564  {
1565    /* We have to make p monic before we return it, so that if the
1566       gcd is a unit in the ground field, we will actually return 1. */
1567    p_Monic(p, r);
1568    return p;
1569  }
1570  else
1571  {
1572    p_PolyDiv(p, q, FALSE, r);
1573    return p_GcdHelper(q, p, r);
1574  }
1575}
1576
1577/* assumes that p and q are univariate polynomials in r,
1578   mentioning the same variable;
1579   assumes a global monomial ordering in r;
1580   assumes that not both p and q are NULL;
1581   returns the gcd of p and q;
1582   leaves p and q unmodified */
1583poly p_Gcd(poly p, poly q, ring r)
1584{
1585  assume((p != NULL) || (q != NULL));
1586
1587  poly a = p; poly b = q;
1588  if (p_Deg(a, r) < p_Deg(b, r)) { a = q; b = p; }
1589  a = p_Copy(a, r); b = p_Copy(b, r);
1590  return p_GcdHelper(a, b, r);
1591}
1592
1593/* see p_ExtGcd;
1594   additional assumption: deg(p) >= deg(q);
1595   must destroy p and q (unless one of them is returned) */
1596poly p_ExtGcdHelper(poly &p, poly &pFactor, poly &q, poly &qFactor,
1597                    ring r)
1598{
1599  if (q == NULL)
1600  {
1601    qFactor = NULL;
1602    pFactor = p_ISet(1, r);
1603    p_SetCoeff(pFactor, n_Invers(p_GetCoeff(p, r), r->cf), r);
1604    p_Monic(p, r);
1605    return p;
1606  }
1607  else
1608  {
1609    poly pDivQ = p_PolyDiv(p, q, TRUE, r);
1610    poly ppFactor = NULL; poly qqFactor = NULL;
1611    poly theGcd = p_ExtGcdHelper(q, qqFactor, p, ppFactor, r);
1612    pFactor = ppFactor;
1613    qFactor = p_Add_q(qqFactor,
1614                      p_Neg(p_Mult_q(pDivQ, p_Copy(ppFactor, r), r), r),
1615                      r);
1616    return theGcd;
1617  }
1618}
1619
1620/* assumes that p and q are univariate polynomials in r,
1621   mentioning the same variable;
1622   assumes a global monomial ordering in r;
1623   assumes that not both p and q are NULL;
1624   returns the gcd of p and q;
1625   moreover, afterwards pFactor and qFactor contain appropriate
1626   factors such that gcd(p, q) = p * pFactor + q * qFactor;
1627   leaves p and q unmodified */
1628poly p_ExtGcd(poly p, poly &pFactor, poly q, poly &qFactor, ring r)
1629{
1630  assume((p != NULL) || (q != NULL));
1631  poly a = p; poly b = q; BOOLEAN aCorrespondsToP = TRUE;
1632  if (p_Deg(a, r) < p_Deg(b, r))
1633  { a = q; b = p; aCorrespondsToP = FALSE; }
1634  a = p_Copy(a, r); b = p_Copy(b, r);
1635  poly aFactor = NULL; poly bFactor = NULL;
1636  poly theGcd = p_ExtGcdHelper(a, aFactor, b, bFactor, r);
1637  if (aCorrespondsToP) { pFactor = aFactor; qFactor = bFactor; }
1638  else                 { pFactor = bFactor; qFactor = aFactor; }
1639  return theGcd;
1640}
1641
1642/*2
1643* returns the partial differentiate of a by the k-th variable
1644* does not destroy the input
1645*/
1646poly p_Diff(poly a, int k, const ring r)
1647{
1648  poly res, f, last;
1649  number t;
1650
1651  last = res = NULL;
1652  while (a!=NULL)
1653  {
1654    if (p_GetExp(a,k,r)!=0)
1655    {
1656      f = p_LmInit(a,r);
1657      t = n_Init(p_GetExp(a,k,r),r->cf);
1658      pSetCoeff0(f,n_Mult(t,pGetCoeff(a),r->cf));
1659      n_Delete(&t,r->cf);
1660      if (n_IsZero(pGetCoeff(f),r->cf))
1661        p_LmDelete(&f,r);
1662      else
1663      {
1664        p_DecrExp(f,k,r);
1665        p_Setm(f,r);
1666        if (res==NULL)
1667        {
1668          res=last=f;
1669        }
1670        else
1671        {
1672          pNext(last)=f;
1673          last=f;
1674        }
1675      }
1676    }
1677    pIter(a);
1678  }
1679  return res;
1680}
1681
1682static poly p_DiffOpM(poly a, poly b,BOOLEAN multiply, const ring r)
1683{
1684  int i,j,s;
1685  number n,h,hh;
1686  poly p=p_One(r);
1687  n=n_Mult(pGetCoeff(a),pGetCoeff(b),r->cf);
1688  for(i=rVar(r);i>0;i--)
1689  {
1690    s=p_GetExp(b,i,r);
1691    if (s<p_GetExp(a,i,r))
1692    {
1693      n_Delete(&n,r->cf);
1694      p_LmDelete(&p,r);
1695      return NULL;
1696    }
1697    if (multiply)
1698    {
1699      for(j=p_GetExp(a,i,r); j>0;j--)
1700      {
1701        h = n_Init(s,r->cf);
1702        hh=n_Mult(n,h,r->cf);
1703        n_Delete(&h,r->cf);
1704        n_Delete(&n,r->cf);
1705        n=hh;
1706        s--;
1707      }
1708      p_SetExp(p,i,s,r);
1709    }
1710    else
1711    {
1712      p_SetExp(p,i,s-p_GetExp(a,i,r),r);
1713    }
1714  }
1715  p_Setm(p,r);
1716  /*if (multiply)*/ p_SetCoeff(p,n,r);
1717  if (n_IsZero(n,r->cf))  p=p_LmDeleteAndNext(p,r); // return NULL as p is a monomial
1718  return p;
1719}
1720
1721poly p_DiffOp(poly a, poly b,BOOLEAN multiply, const ring r)
1722{
1723  poly result=NULL;
1724  poly h;
1725  for(;a!=NULL;pIter(a))
1726  {
1727    for(h=b;h!=NULL;pIter(h))
1728    {
1729      result=p_Add_q(result,p_DiffOpM(a,h,multiply,r),r);
1730    }
1731  }
1732  return result;
1733}
1734/*2
1735* subtract p2 from p1, p1 and p2 are destroyed
1736* do not put attention on speed: the procedure is only used in the interpreter
1737*/
1738poly p_Sub(poly p1, poly p2, const ring r)
1739{
1740  return p_Add_q(p1, p_Neg(p2,r),r);
1741}
1742
1743/*3
1744* compute for a monomial m
1745* the power m^exp, exp > 1
1746* destroys p
1747*/
1748static poly p_MonPower(poly p, int exp, const ring r)
1749{
1750  int i;
1751
1752  if(!n_IsOne(pGetCoeff(p),r->cf))
1753  {
1754    number x, y;
1755    y = pGetCoeff(p);
1756    n_Power(y,exp,&x,r->cf);
1757    n_Delete(&y,r->cf);
1758    pSetCoeff0(p,x);
1759  }
1760  for (i=rVar(r); i!=0; i--)
1761  {
1762    p_MultExp(p,i, exp,r);
1763  }
1764  p_Setm(p,r);
1765  return p;
1766}
1767
1768/*3
1769* compute for monomials p*q
1770* destroys p, keeps q
1771*/
1772static void p_MonMult(poly p, poly q, const ring r)
1773{
1774  number x, y;
1775
1776  y = pGetCoeff(p);
1777  x = n_Mult(y,pGetCoeff(q),r->cf);
1778  n_Delete(&y,r->cf);
1779  pSetCoeff0(p,x);
1780  //for (int i=pVariables; i!=0; i--)
1781  //{
1782  //  pAddExp(p,i, pGetExp(q,i));
1783  //}
1784  //p->Order += q->Order;
1785  p_ExpVectorAdd(p,q,r);
1786}
1787
1788/*3
1789* compute for monomials p*q
1790* keeps p, q
1791*/
1792static poly p_MonMultC(poly p, poly q, const ring rr)
1793{
1794  number x;
1795  poly r = p_Init(rr);
1796
1797  x = n_Mult(pGetCoeff(p),pGetCoeff(q),rr->cf);
1798  pSetCoeff0(r,x);
1799  p_ExpVectorSum(r,p, q, rr);
1800  return r;
1801}
1802
1803/*3
1804*  create binomial coef.
1805*/
1806static number* pnBin(int exp, const ring r)
1807{
1808  int e, i, h;
1809  number x, y, *bin=NULL;
1810
1811  x = n_Init(exp,r->cf);
1812  if (n_IsZero(x,r->cf))
1813  {
1814    n_Delete(&x,r->cf);
1815    return bin;
1816  }
1817  h = (exp >> 1) + 1;
1818  bin = (number *)omAlloc0(h*sizeof(number));
1819  bin[1] = x;
1820  if (exp < 4)
1821    return bin;
1822  i = exp - 1;
1823  for (e=2; e<h; e++)
1824  {
1825      x = n_Init(i,r->cf);
1826      i--;
1827      y = n_Mult(x,bin[e-1],r->cf);
1828      n_Delete(&x,r->cf);
1829      x = n_Init(e,r->cf);
1830      bin[e] = n_IntDiv(y,x,r->cf);
1831      n_Delete(&x,r->cf);
1832      n_Delete(&y,r->cf);
1833  }
1834  return bin;
1835}
1836
1837static void pnFreeBin(number *bin, int exp,const coeffs r)
1838{
1839  int e, h = (exp >> 1) + 1;
1840
1841  if (bin[1] != NULL)
1842  {
1843    for (e=1; e<h; e++)
1844      n_Delete(&(bin[e]),r);
1845  }
1846  omFreeSize((ADDRESS)bin, h*sizeof(number));
1847}
1848
1849/*
1850*  compute for a poly p = head+tail, tail is monomial
1851*          (head + tail)^exp, exp > 1
1852*          with binomial coef.
1853*/
1854static poly p_TwoMonPower(poly p, int exp, const ring r)
1855{
1856  int eh, e;
1857  long al;
1858  poly *a;
1859  poly tail, b, res, h;
1860  number x;
1861  number *bin = pnBin(exp,r);
1862
1863  tail = pNext(p);
1864  if (bin == NULL)
1865  {
1866    p_MonPower(p,exp,r);
1867    p_MonPower(tail,exp,r);
1868#ifdef PDEBUG
1869    p_Test(p,r);
1870#endif
1871    return p;
1872  }
1873  eh = exp >> 1;
1874  al = (exp + 1) * sizeof(poly);
1875  a = (poly *)omAlloc(al);
1876  a[1] = p;
1877  for (e=1; e<exp; e++)
1878  {
1879    a[e+1] = p_MonMultC(a[e],p,r);
1880  }
1881  res = a[exp];
1882  b = p_Head(tail,r);
1883  for (e=exp-1; e>eh; e--)
1884  {
1885    h = a[e];
1886    x = n_Mult(bin[exp-e],pGetCoeff(h),r->cf);
1887    p_SetCoeff(h,x,r);
1888    p_MonMult(h,b,r);
1889    res = pNext(res) = h;
1890    p_MonMult(b,tail,r);
1891  }
1892  for (e=eh; e!=0; e--)
1893  {
1894    h = a[e];
1895    x = n_Mult(bin[e],pGetCoeff(h),r->cf);
1896    p_SetCoeff(h,x,r);
1897    p_MonMult(h,b,r);
1898    res = pNext(res) = h;
1899    p_MonMult(b,tail,r);
1900  }
1901  p_LmDelete(&tail,r);
1902  pNext(res) = b;
1903  pNext(b) = NULL;
1904  res = a[exp];
1905  omFreeSize((ADDRESS)a, al);
1906  pnFreeBin(bin, exp, r->cf);
1907//  tail=res;
1908// while((tail!=NULL)&&(pNext(tail)!=NULL))
1909// {
1910//   if(nIsZero(pGetCoeff(pNext(tail))))
1911//   {
1912//     pLmDelete(&pNext(tail));
1913//   }
1914//   else
1915//     pIter(tail);
1916// }
1917#ifdef PDEBUG
1918  p_Test(res,r);
1919#endif
1920  return res;
1921}
1922
1923static poly p_Pow(poly p, int i, const ring r)
1924{
1925  poly rc = p_Copy(p,r);
1926  i -= 2;
1927  do
1928  {
1929    rc = p_Mult_q(rc,p_Copy(p,r),r);
1930    p_Normalize(rc,r);
1931    i--;
1932  }
1933  while (i != 0);
1934  return p_Mult_q(rc,p,r);
1935}
1936
1937/*2
1938* returns the i-th power of p
1939* p will be destroyed
1940*/
1941poly p_Power(poly p, int i, const ring r)
1942{
1943  poly rc=NULL;
1944
1945  if (i==0)
1946  {
1947    p_Delete(&p,r);
1948    return p_One(r);
1949  }
1950
1951  if(p!=NULL)
1952  {
1953    if ( (i > 0) && ((unsigned long ) i > (r->bitmask)))
1954    {
1955      Werror("exponent %d is too large, max. is %ld",i,r->bitmask);
1956      return NULL;
1957    }
1958    switch (i)
1959    {
1960// cannot happen, see above
1961//      case 0:
1962//      {
1963//        rc=pOne();
1964//        pDelete(&p);
1965//        break;
1966//      }
1967      case 1:
1968        rc=p;
1969        break;
1970      case 2:
1971        rc=p_Mult_q(p_Copy(p,r),p,r);
1972        break;
1973      default:
1974        if (i < 0)
1975        {
1976          p_Delete(&p,r);
1977          return NULL;
1978        }
1979        else
1980        {
1981#ifdef HAVE_PLURAL
1982          if (rIsPluralRing(r)) /* in the NC case nothing helps :-( */
1983          {
1984            int j=i;
1985            rc = p_Copy(p,r);
1986            while (j>1)
1987            {
1988              rc = p_Mult_q(p_Copy(p,r),rc,r);
1989              j--;
1990            }
1991            p_Delete(&p,r);
1992            return rc;
1993          }
1994#endif
1995          rc = pNext(p);
1996          if (rc == NULL)
1997            return p_MonPower(p,i,r);
1998          /* else: binom ?*/
1999          int char_p=rChar(r);
2000          if ((pNext(rc) != NULL)
2001#ifdef HAVE_RINGS
2002             || rField_is_Ring(r)
2003#endif
2004             )
2005            return p_Pow(p,i,r);
2006          if ((char_p==0) || (i<=char_p))
2007            return p_TwoMonPower(p,i,r);
2008          return p_Pow(p,i,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.