source: git/libpolys/polys/monomials/p_polys.cc @ 5a72fe

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