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

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