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

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