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

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