source: git/kernel/polys1.cc @ 0885fd

spielwiese
Last change on this file since 0885fd was 0885fd, checked in by Hans Schoenemann <hannes@…>, 13 years ago
janet: pSimpleContent git-svn-id: file:///usr/local/Singular/svn/trunk@13706 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 36.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5
6/*
7* ABSTRACT - all basic methods to manipulate polynomials:
8* independent of representation
9*/
10
11/* includes */
12#include <string.h>
13#include <kernel/mod2.h>
14#include <kernel/options.h>
15#include <kernel/numbers.h>
16#include <kernel/ffields.h>
17#include <omalloc/omalloc.h>
18#include <kernel/febase.h>
19#include <kernel/weight.h>
20#include <kernel/intvec.h>
21#include <kernel/longalg.h>
22#include <kernel/ring.h>
23#include <kernel/ideals.h>
24#include <kernel/polys.h>
25//#include "ipid.h"
26#ifdef HAVE_FACTORY
27#include <kernel/clapsing.h>
28#endif
29
30#ifdef HAVE_RATGRING
31#include <kernel/ratgring.h>
32#endif
33
34/*-------- several access procedures to monomials -------------------- */
35/*
36* the module weights for std
37*/
38static pFDegProc pOldFDeg;
39static pLDegProc pOldLDeg;
40static intvec * pModW;
41static BOOLEAN pOldLexOrder;
42
43static long pModDeg(poly p, ring r = currRing)
44{
45  long d=pOldFDeg(p, r);
46  int c=p_GetComp(p, r);
47  if ((c>0) && (pModW->range(c-1))) d+= (*pModW)[c-1];
48  return d;
49  //return pOldFDeg(p, r)+(*pModW)[p_GetComp(p, r)-1];
50}
51
52void pSetModDeg(intvec *w)
53{
54  if (w!=NULL)
55  {
56    pModW = w;
57    pOldFDeg = pFDeg;
58    pOldLDeg = pLDeg;
59    pOldLexOrder = pLexOrder;
60    pSetDegProcs(pModDeg);
61    pLexOrder = TRUE;
62  }
63  else
64  {
65    pModW = NULL;
66    pRestoreDegProcs(pOldFDeg, pOldLDeg);
67    pLexOrder = pOldLexOrder;
68  }
69}
70
71
72/*2
73* subtract p2 from p1, p1 and p2 are destroyed
74* do not put attention on speed: the procedure is only used in the interpreter
75*/
76poly pSub(poly p1, poly p2)
77{
78  return pAdd(p1, pNeg(p2));
79}
80
81/*3
82*  create binomial coef.
83*/
84static number* pnBin(int exp)
85{
86  int e, i, h;
87  number x, y, *bin=NULL;
88
89  x = nInit(exp);
90  if (nIsZero(x))
91  {
92    nDelete(&x);
93    return bin;
94  }
95  h = (exp >> 1) + 1;
96  bin = (number *)omAlloc0(h*sizeof(number));
97  bin[1] = x;
98  if (exp < 4)
99    return bin;
100  i = exp - 1;
101  for (e=2; e<h; e++)
102  {
103      x = nInit(i);
104      i--;
105      y = nMult(x,bin[e-1]);
106      nDelete(&x);
107      x = nInit(e);
108      bin[e] = nIntDiv(y,x);
109      nDelete(&x);
110      nDelete(&y);
111  }
112  return bin;
113}
114
115static void pnFreeBin(number *bin, int exp)
116{
117  int e, h = (exp >> 1) + 1;
118
119  if (bin[1] != NULL)
120  {
121    for (e=1; e<h; e++)
122      nDelete(&(bin[e]));
123  }
124  omFreeSize((ADDRESS)bin, h*sizeof(number));
125}
126
127/*3
128* compute for a monomial m
129* the power m^exp, exp > 1
130* destroys p
131*/
132static poly p_MonPower(poly p, int exp, const ring r)
133{
134  int i;
135
136  if(!n_IsOne(pGetCoeff(p),r))
137  {
138    number x, y;
139    y = pGetCoeff(p);
140    n_Power(y,exp,&x,r);
141    n_Delete(&y,r);
142    pSetCoeff0(p,x);
143  }
144  for (i=rVar(r); i!=0; i--)
145  {
146    p_MultExp(p,i, exp,r);
147  }
148  p_Setm(p,r);
149  return p;
150}
151
152/*3
153* compute for monomials p*q
154* destroys p, keeps q
155*/
156static void p_MonMult(poly p, poly q, const ring r)
157{
158  number x, y;
159  int i;
160
161  y = pGetCoeff(p);
162  x = n_Mult(y,pGetCoeff(q),r);
163  n_Delete(&y,r);
164  pSetCoeff0(p,x);
165  //for (i=pVariables; i!=0; i--)
166  //{
167  //  pAddExp(p,i, pGetExp(q,i));
168  //}
169  //p->Order += q->Order;
170  p_ExpVectorAdd(p,q,r);
171}
172
173/*3
174* compute for monomials p*q
175* keeps p, q
176*/
177static poly p_MonMultC(poly p, poly q, const ring rr)
178{
179  number x;
180  int i;
181  poly r = p_Init(rr);
182
183  x = n_Mult(pGetCoeff(p),pGetCoeff(q),rr);
184  pSetCoeff0(r,x);
185  p_ExpVectorSum(r,p, q, rr);
186  return r;
187}
188
189/*
190*  compute for a poly p = head+tail, tail is monomial
191*          (head + tail)^exp, exp > 1
192*          with binomial coef.
193*/
194static poly p_TwoMonPower(poly p, int exp, const ring r)
195{
196  int eh, e;
197  long al;
198  poly *a;
199  poly tail, b, res, h;
200  number x;
201  number *bin = pnBin(exp);
202
203  tail = pNext(p);
204  if (bin == NULL)
205  {
206    p_MonPower(p,exp,r);
207    p_MonPower(tail,exp,r);
208#ifdef PDEBUG
209    p_Test(p,r);
210#endif
211    return p;
212  }
213  eh = exp >> 1;
214  al = (exp + 1) * sizeof(poly);
215  a = (poly *)omAlloc(al);
216  a[1] = p;
217  for (e=1; e<exp; e++)
218  {
219    a[e+1] = p_MonMultC(a[e],p,r);
220  }
221  res = a[exp];
222  b = p_Head(tail,r);
223  for (e=exp-1; e>eh; e--)
224  {
225    h = a[e];
226    x = n_Mult(bin[exp-e],pGetCoeff(h),r);
227    p_SetCoeff(h,x,r);
228    p_MonMult(h,b,r);
229    res = pNext(res) = h;
230    p_MonMult(b,tail,r);
231  }
232  for (e=eh; e!=0; e--)
233  {
234    h = a[e];
235    x = n_Mult(bin[e],pGetCoeff(h),r);
236    p_SetCoeff(h,x,r);
237    p_MonMult(h,b,r);
238    res = pNext(res) = h;
239    p_MonMult(b,tail,r);
240  }
241  p_LmDelete(&tail,r);
242  pNext(res) = b;
243  pNext(b) = NULL;
244  res = a[exp];
245  omFreeSize((ADDRESS)a, al);
246  pnFreeBin(bin, exp);
247//  tail=res;
248// while((tail!=NULL)&&(pNext(tail)!=NULL))
249// {
250//   if(nIsZero(pGetCoeff(pNext(tail))))
251//   {
252//     pLmDelete(&pNext(tail));
253//   }
254//   else
255//     pIter(tail);
256// }
257#ifdef PDEBUG
258  p_Test(res,r);
259#endif
260  return res;
261}
262
263static poly p_Pow(poly p, int i, const ring r)
264{
265  poly rc = p_Copy(p,r);
266  i -= 2;
267  do
268  {
269    rc = p_Mult_q(rc,p_Copy(p,r),r);
270    p_Normalize(rc,r);
271    i--;
272  }
273  while (i != 0);
274  return p_Mult_q(rc,p,r);
275}
276
277/*2
278* returns the i-th power of p
279* p will be destroyed
280*/
281poly p_Power(poly p, int i, const ring r)
282{
283  poly rc=NULL;
284
285  if (i==0)
286  {
287    p_Delete(&p,r);
288    return p_One(r);
289  }
290
291  if(p!=NULL)
292  {
293    if ( (i > 0) && ((unsigned long ) i > (r->bitmask)))
294    {
295      Werror("exponent %d is too large, max. is %ld",i,r->bitmask);
296      return NULL;
297    }
298    switch (i)
299    {
300// cannot happen, see above
301//      case 0:
302//      {
303//        rc=pOne();
304//        pDelete(&p);
305//        break;
306//      }
307      case 1:
308        rc=p;
309        break;
310      case 2:
311        rc=p_Mult_q(p_Copy(p,r),p,r);
312        break;
313      default:
314        if (i < 0)
315        {
316          p_Delete(&p,r);
317          return NULL;
318        }
319        else
320        {
321#ifdef HAVE_PLURAL
322          if (rIsPluralRing(r)) /* in the NC case nothing helps :-( */
323          {
324            int j=i;
325            rc = p_Copy(p,r);
326            while (j>1)
327            {
328              rc = p_Mult_q(p_Copy(p,r),rc,r);
329              j--;
330            }
331            p_Delete(&p,r);
332            return rc;
333          }
334#endif
335          rc = pNext(p);
336          if (rc == NULL)
337            return p_MonPower(p,i,r);
338          /* else: binom ?*/
339          int char_p=rChar(r);
340          if ((pNext(rc) != NULL)
341#ifdef HAVE_RINGS
342             || rField_is_Ring(r)
343#endif
344             )
345            return p_Pow(p,i,r);
346          if ((char_p==0) || (i<=char_p))
347            return p_TwoMonPower(p,i,r);
348          poly p_p=p_TwoMonPower(p_Copy(p,r),char_p,r);
349          return p_Mult_q(p_Power(p,i-char_p,r),p_p,r);
350        }
351      /*end default:*/
352    }
353  }
354  return rc;
355}
356
357/*2
358* returns the partial differentiate of a by the k-th variable
359* does not destroy the input
360*/
361poly pDiff(poly a, int k)
362{
363  poly res, f, last;
364  number t;
365
366  last = res = NULL;
367  while (a!=NULL)
368  {
369    if (pGetExp(a,k)!=0)
370    {
371      f = pLmInit(a);
372      t = nInit(pGetExp(a,k));
373      pSetCoeff0(f,nMult(t,pGetCoeff(a)));
374      nDelete(&t);
375      if (nIsZero(pGetCoeff(f)))
376        pLmDelete(&f);
377      else
378      {
379        pDecrExp(f,k);
380        pSetm(f);
381        if (res==NULL)
382        {
383          res=last=f;
384        }
385        else
386        {
387          pNext(last)=f;
388          last=f;
389        }
390      }
391    }
392    pIter(a);
393  }
394  return res;
395}
396
397static poly pDiffOpM(poly a, poly b,BOOLEAN multiply)
398{
399  int i,j,s;
400  number n,h,hh;
401  poly p=pOne();
402  n=nMult(pGetCoeff(a),pGetCoeff(b));
403  for(i=pVariables;i>0;i--)
404  {
405    s=pGetExp(b,i);
406    if (s<pGetExp(a,i))
407    {
408      nDelete(&n);
409      pLmDelete(&p);
410      return NULL;
411    }
412    if (multiply)
413    {
414      for(j=pGetExp(a,i); j>0;j--)
415      {
416        h = nInit(s);
417        hh=nMult(n,h);
418        nDelete(&h);
419        nDelete(&n);
420        n=hh;
421        s--;
422      }
423      pSetExp(p,i,s);
424    }
425    else
426    {
427      pSetExp(p,i,s-pGetExp(a,i));
428    }
429  }
430  pSetm(p);
431  /*if (multiply)*/ pSetCoeff(p,n);
432  if (nIsZero(n))  p=pLmDeleteAndNext(p); // return NULL as p is a monomial
433  return p;
434}
435
436poly pDiffOp(poly a, poly b,BOOLEAN multiply)
437{
438  poly result=NULL;
439  poly h;
440  for(;a!=NULL;pIter(a))
441  {
442    for(h=b;h!=NULL;pIter(h))
443    {
444      result=pAdd(result,pDiffOpM(a,h,multiply));
445    }
446  }
447  return result;
448}
449
450
451void pSplit(poly p, poly *h)
452{
453  *h=pNext(p);
454  pNext(p)=NULL;
455}
456
457
458
459int pMaxCompProc(poly p)
460{
461  return pMaxComp(p);
462}
463
464/*2
465* handle memory request for sets of polynomials (ideals)
466* l is the length of *p, increment is the difference (may be negative)
467*/
468void pEnlargeSet(polyset *p, int l, int increment)
469{
470  int i;
471  polyset h;
472
473  h=(polyset)omReallocSize((poly*)*p,l*sizeof(poly),(l+increment)*sizeof(poly));
474  if (increment>0)
475  {
476    //for (i=l; i<l+increment; i++)
477    //  h[i]=NULL;
478    memset(&(h[l]),0,increment*sizeof(poly));
479  }
480  *p=h;
481}
482
483number pInitContent(poly ph);
484number pInitContent_a(poly ph);
485
486void p_Content(poly ph, const ring r)
487{
488#ifdef HAVE_RINGS
489  if (rField_is_Ring(r))
490  {
491    if ((ph!=NULL) && rField_has_Units(r))
492    {
493      number k = nGetUnit(pGetCoeff(ph));
494      if (!nIsOne(k))
495      {
496        number tmpGMP = k;
497        k = nInvers(k);
498        nDelete(&tmpGMP);
499        poly h = pNext(ph);
500        pSetCoeff(ph, nMult(pGetCoeff(ph), k));
501        while (h != NULL)
502        {
503          pSetCoeff(h, nMult(pGetCoeff(h), k));
504          pIter(h);
505        }
506      }
507      nDelete(&k);
508    }
509    return;
510  }
511#endif
512  number h,d;
513  poly p;
514
515  if(TEST_OPT_CONTENTSB) return;
516  if(pNext(ph)==NULL)
517  {
518    pSetCoeff(ph,nInit(1));
519  }
520  else
521  {
522    nNormalize(pGetCoeff(ph));
523    if(!nGreaterZero(pGetCoeff(ph))) ph = pNeg(ph);
524    if (rField_is_Q())
525    {
526      h=pInitContent(ph);
527      p=ph;
528    }
529    else if ((rField_is_Extension(r))
530    && ((rPar(r)>1)||(r->minpoly==NULL)))
531    {
532      h=pInitContent_a(ph);
533      p=ph;
534    }
535    else
536    {
537      h=nCopy(pGetCoeff(ph));
538      p = pNext(ph);
539    }
540    while (p!=NULL)
541    {
542      nNormalize(pGetCoeff(p));
543      d=nGcd(h,pGetCoeff(p),r);
544      nDelete(&h);
545      h = d;
546      if(nIsOne(h))
547      {
548        break;
549      }
550      pIter(p);
551    }
552    p = ph;
553    //number tmp;
554    if(!nIsOne(h))
555    {
556      while (p!=NULL)
557      {
558        //d = nDiv(pGetCoeff(p),h);
559        //tmp = nIntDiv(pGetCoeff(p),h);
560        //if (!nEqual(d,tmp))
561        //{
562        //  StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
563        //  nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
564        //  nWrite(tmp);Print(StringAppendS("\n"));
565        //}
566        //nDelete(&tmp);
567        d = nIntDiv(pGetCoeff(p),h);
568        pSetCoeff(p,d);
569        pIter(p);
570      }
571    }
572    nDelete(&h);
573#ifdef HAVE_FACTORY
574    if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
575    {
576      singclap_divide_content(ph);
577      if(!nGreaterZero(pGetCoeff(ph))) ph = pNeg(ph);
578    }
579#endif
580    if (rField_is_Q_a(r))
581    {
582      number hzz = nlInit(1, r);
583      h = nlInit(1, r);
584      p=ph;
585      while (p!=NULL)
586      { // each monom: coeff in Q_a
587        lnumber c_n_n=(lnumber)pGetCoeff(p);
588        napoly c_n=c_n_n->z;
589        while (c_n!=NULL)
590        { // each monom: coeff in Q
591          d=nlLcm(hzz,pGetCoeff(c_n),r->algring);
592          n_Delete(&hzz,r->algring);
593          hzz=d;
594          pIter(c_n);
595        }
596        c_n=c_n_n->n;
597        while (c_n!=NULL)
598        { // each monom: coeff in Q
599          d=nlLcm(h,pGetCoeff(c_n),r->algring);
600          n_Delete(&h,r->algring);
601          h=d;
602          pIter(c_n);
603        }
604        pIter(p);
605      }
606      /* hzz contains the 1/lcm of all denominators in c_n_n->z*/
607      /* h contains the 1/lcm of all denominators in c_n_n->n*/
608      number htmp=nlInvers(h);
609      number hzztmp=nlInvers(hzz);
610      number hh=nlMult(hzz,h);
611      nlDelete(&hzz,r->algring);
612      nlDelete(&h,r->algring);
613      number hg=nlGcd(hzztmp,htmp,r->algring);
614      nlDelete(&hzztmp,r->algring);
615      nlDelete(&htmp,r->algring);
616      h=nlMult(hh,hg);
617      nlDelete(&hg,r->algring);
618      nlDelete(&hh,r->algring);
619      nlNormalize(h);
620      if(!nlIsOne(h))
621      {
622        p=ph;
623        while (p!=NULL)
624        { // each monom: coeff in Q_a
625          lnumber c_n_n=(lnumber)pGetCoeff(p);
626          napoly c_n=c_n_n->z;
627          while (c_n!=NULL)
628          { // each monom: coeff in Q
629            d=nlMult(h,pGetCoeff(c_n));
630            nlNormalize(d);
631            nlDelete(&pGetCoeff(c_n),r->algring);
632            pGetCoeff(c_n)=d;
633            pIter(c_n);
634          }
635          c_n=c_n_n->n;
636          while (c_n!=NULL)
637          { // each monom: coeff in Q
638            d=nlMult(h,pGetCoeff(c_n));
639            nlNormalize(d);
640            nlDelete(&pGetCoeff(c_n),r->algring);
641            pGetCoeff(c_n)=d;
642            pIter(c_n);
643          }
644          pIter(p);
645        }
646      }
647      nlDelete(&h,r->algring);
648    }
649  }
650}
651
652void pSimpleContent(poly ph,int smax)
653{
654  if(TEST_OPT_CONTENTSB) return;
655  if (ph==NULL) return;
656  if (pNext(ph)==NULL)
657  {
658    pSetCoeff(ph,nInit(1));
659    return;
660  }
661  if ((pNext(pNext(ph))==NULL)||(!rField_is_Q()))
662  {
663    return;
664  }
665  number d=pInitContent(ph);
666  if (nlSize(d)<=smax)
667  {
668    //if (TEST_OPT_PROT) PrintS("G");
669    return;
670  }
671  poly p=ph;
672  number h=d;
673  if (smax==1) smax=2;
674  while (p!=NULL)
675  {
676#if 0
677    d=nlGcd(h,pGetCoeff(p),currRing);
678    nlDelete(&h,currRing);
679    h = d;
680#else
681    nlInpGcd(h,pGetCoeff(p),currRing);
682#endif
683    if(nlSize(h)<smax)
684    {
685      //if (TEST_OPT_PROT) PrintS("g");
686      return;
687    }
688    pIter(p);
689  }
690  p = ph;
691  if (!nlGreaterZero(pGetCoeff(p))) h=nlNeg(h);
692  if(nlIsOne(h)) return;
693  //if (TEST_OPT_PROT) PrintS("c");
694  while (p!=NULL)
695  {
696#if 1
697    d = nlIntDiv(pGetCoeff(p),h);
698    pSetCoeff(p,d);
699#else
700    nlInpIntDiv(pGetCoeff(p),h,currRing);
701#endif
702    pIter(p);
703  }
704  nlDelete(&h,currRing);
705}
706
707number pInitContent(poly ph)
708// only for coefficients in Q
709#if 0
710{
711  assume(!TEST_OPT_CONTENTSB);
712  assume(ph!=NULL);
713  assume(pNext(ph)!=NULL);
714  assume(rField_is_Q());
715  if (pNext(pNext(ph))==NULL)
716  {
717    return nlGetNom(pGetCoeff(pNext(ph)),currRing);
718  }
719  poly p=ph;
720  number n1=nlGetNom(pGetCoeff(p),currRing);
721  pIter(p);
722  number n2=nlGetNom(pGetCoeff(p),currRing);
723  pIter(p);
724  number d;
725  number t;
726  loop
727  {
728    nlNormalize(pGetCoeff(p));
729    t=nlGetNom(pGetCoeff(p),currRing);
730    if (nlGreaterZero(t))
731      d=nlAdd(n1,t);
732    else
733      d=nlSub(n1,t);
734    nlDelete(&t,currRing);
735    nlDelete(&n1,currRing);
736    n1=d;
737    pIter(p);
738    if (p==NULL) break;
739    nlNormalize(pGetCoeff(p));
740    t=nlGetNom(pGetCoeff(p),currRing);
741    if (nlGreaterZero(t))
742      d=nlAdd(n2,t);
743    else
744      d=nlSub(n2,t);
745    nlDelete(&t,currRing);
746    nlDelete(&n2,currRing);
747    n2=d;
748    pIter(p);
749    if (p==NULL) break;
750  }
751  d=nlGcd(n1,n2,currRing);
752  nlDelete(&n1,currRing);
753  nlDelete(&n2,currRing);
754  return d;
755}
756#else
757{
758  number d=pGetCoeff(ph);
759  if(SR_HDL(d)&SR_INT) return d;
760  int s=mpz_size1(d->z);
761  int s2=-1;
762  number d2;
763  loop
764  {
765    pIter(ph);
766    if(ph==NULL)
767    {
768      if (s2==-1) return nlCopy(d);
769      break;
770    }
771    if (SR_HDL(pGetCoeff(ph))&SR_INT)
772    {
773      s2=s;
774      d2=d;
775      s=0;
776      d=pGetCoeff(ph);
777      if (s2==0) break;
778    }
779    else
780    if (mpz_size1((pGetCoeff(ph)->z))<=s)
781    {
782      s2=s;
783      d2=d;
784      d=pGetCoeff(ph);
785      s=mpz_size1(d->z);
786    }
787  }
788  return nlGcd(d,d2,currRing);
789}
790#endif
791
792number pInitContent_a(poly ph)
793// only for coefficients in K(a) anf K(a,...)
794{
795  number d=pGetCoeff(ph);
796  int s=naParDeg(d);
797  if (s /* naParDeg(d)*/ <=1) return naCopy(d);
798  int s2=-1;
799  number d2;
800  int ss;
801  loop
802  {
803    pIter(ph);
804    if(ph==NULL)
805    {
806      if (s2==-1) return naCopy(d);
807      break;
808    }
809    if ((ss=naParDeg(pGetCoeff(ph)))<s)
810    {
811      s2=s;
812      d2=d;
813      s=ss;
814      d=pGetCoeff(ph);
815      if (s2<=1) break;
816    }
817  }
818  return naGcd(d,d2,currRing);
819}
820
821
822//void pContent(poly ph)
823//{
824//  number h,d;
825//  poly p;
826//
827//  p = ph;
828//  if(pNext(p)==NULL)
829//  {
830//    pSetCoeff(p,nInit(1));
831//  }
832//  else
833//  {
834//#ifdef PDEBUG
835//    if (!pTest(p)) return;
836//#endif
837//    nNormalize(pGetCoeff(p));
838//    if(!nGreaterZero(pGetCoeff(ph)))
839//    {
840//      ph = pNeg(ph);
841//      nNormalize(pGetCoeff(p));
842//    }
843//    h=pGetCoeff(p);
844//    pIter(p);
845//    while (p!=NULL)
846//    {
847//      nNormalize(pGetCoeff(p));
848//      if (nGreater(h,pGetCoeff(p))) h=pGetCoeff(p);
849//      pIter(p);
850//    }
851//    h=nCopy(h);
852//    p=ph;
853//    while (p!=NULL)
854//    {
855//      d=nGcd(h,pGetCoeff(p));
856//      nDelete(&h);
857//      h = d;
858//      if(nIsOne(h))
859//      {
860//        break;
861//      }
862//      pIter(p);
863//    }
864//    p = ph;
865//    //number tmp;
866//    if(!nIsOne(h))
867//    {
868//      while (p!=NULL)
869//      {
870//        d = nIntDiv(pGetCoeff(p),h);
871//        pSetCoeff(p,d);
872//        pIter(p);
873//      }
874//    }
875//    nDelete(&h);
876//#ifdef HAVE_FACTORY
877//    if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
878//    {
879//      pTest(ph);
880//      singclap_divide_content(ph);
881//      pTest(ph);
882//    }
883//#endif
884//  }
885//}
886#if 0
887void p_Content(poly ph, ring r)
888{
889  number h,d;
890  poly p;
891
892  if(pNext(ph)==NULL)
893  {
894    pSetCoeff(ph,n_Init(1,r));
895  }
896  else
897  {
898    n_Normalize(pGetCoeff(ph),r);
899    if(!n_GreaterZero(pGetCoeff(ph),r)) ph = p_Neg(ph,r);
900    h=n_Copy(pGetCoeff(ph),r);
901    p = pNext(ph);
902    while (p!=NULL)
903    {
904      n_Normalize(pGetCoeff(p),r);
905      d=n_Gcd(h,pGetCoeff(p),r);
906      n_Delete(&h,r);
907      h = d;
908      if(n_IsOne(h,r))
909      {
910        break;
911      }
912      pIter(p);
913    }
914    p = ph;
915    //number tmp;
916    if(!n_IsOne(h,r))
917    {
918      while (p!=NULL)
919      {
920        //d = nDiv(pGetCoeff(p),h);
921        //tmp = nIntDiv(pGetCoeff(p),h);
922        //if (!nEqual(d,tmp))
923        //{
924        //  StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
925        //  nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
926        //  nWrite(tmp);Print(StringAppendS("\n"));
927        //}
928        //nDelete(&tmp);
929        d = n_IntDiv(pGetCoeff(p),h,r);
930        p_SetCoeff(p,d,r);
931        pIter(p);
932      }
933    }
934    n_Delete(&h,r);
935#ifdef HAVE_FACTORY
936    //if ( (n_GetChar(r) == 1) || (n_GetChar(r) < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
937    //{
938    //  singclap_divide_content(ph);
939    //  if(!n_GreaterZero(pGetCoeff(ph),r)) ph = p_Neg(ph,r);
940    //}
941#endif
942  }
943}
944#endif
945
946poly p_Cleardenom(poly ph, const ring r)
947{
948  poly start=ph;
949  number d, h;
950  poly p;
951
952#ifdef HAVE_RINGS
953  if (rField_is_Ring(r))
954  {
955    p_Content(ph,r);
956    return start;
957  }
958#endif
959  if (rField_is_Zp(r) && TEST_OPT_INTSTRATEGY) return start;
960  p = ph;
961  if(pNext(p)==NULL)
962  {
963    if (TEST_OPT_CONTENTSB)
964    {
965      number n=nGetDenom(pGetCoeff(p));
966      if (!nIsOne(n))
967      {
968        number nn=nMult(pGetCoeff(p),n);
969        nNormalize(nn);
970        pSetCoeff(p,nn);
971      }
972      nDelete(&n);
973    }
974    else
975      pSetCoeff(p,nInit(1));
976  }
977  else
978  {
979    h = nInit(1);
980    while (p!=NULL)
981    {
982      nNormalize(pGetCoeff(p));
983      d=nLcm(h,pGetCoeff(p),currRing);
984      nDelete(&h);
985      h=d;
986      pIter(p);
987    }
988    /* contains the 1/lcm of all denominators */
989    if(!nIsOne(h))
990    {
991      p = ph;
992      while (p!=NULL)
993      {
994        /* should be:
995        * number hh;
996        * nGetDenom(p->coef,&hh);
997        * nMult(&h,&hh,&d);
998        * nNormalize(d);
999        * nDelete(&hh);
1000        * nMult(d,p->coef,&hh);
1001        * nDelete(&d);
1002        * nDelete(&(p->coef));
1003        * p->coef =hh;
1004        */
1005        d=nMult(h,pGetCoeff(p));
1006        nNormalize(d);
1007        pSetCoeff(p,d);
1008        pIter(p);
1009      }
1010      nDelete(&h);
1011      if (nGetChar()==1)
1012      {
1013        loop
1014        {
1015          h = nInit(1);
1016          p=ph;
1017          while (p!=NULL)
1018          {
1019            d=nLcm(h,pGetCoeff(p),currRing);
1020            nDelete(&h);
1021            h=d;
1022            pIter(p);
1023          }
1024          /* contains the 1/lcm of all denominators */
1025          if(!nIsOne(h))
1026          {
1027            p = ph;
1028            while (p!=NULL)
1029            {
1030              /* should be:
1031              * number hh;
1032              * nGetDenom(p->coef,&hh);
1033              * nMult(&h,&hh,&d);
1034              * nNormalize(d);
1035              * nDelete(&hh);
1036              * nMult(d,p->coef,&hh);
1037              * nDelete(&d);
1038              * nDelete(&(p->coef));
1039              * p->coef =hh;
1040              */
1041              d=nMult(h,pGetCoeff(p));
1042              nNormalize(d);
1043              pSetCoeff(p,d);
1044              pIter(p);
1045            }
1046            nDelete(&h);
1047          }
1048          else
1049          {
1050            nDelete(&h);
1051            break;
1052          }
1053        }
1054      }
1055    }
1056    if (h!=NULL) nDelete(&h);
1057 
1058    p_Content(ph,r);
1059#ifdef HAVE_RATGRING
1060    if (rIsRatGRing(r))
1061    {
1062      /* quick unit detection in the rational case is done in gr_nc_bba */
1063      pContentRat(ph);
1064      start=ph;
1065    }
1066#endif
1067  }
1068  return start;
1069}
1070
1071void p_Cleardenom_n(poly ph,const ring r,number &c)
1072{
1073  number d, h;
1074  poly p;
1075
1076  p = ph;
1077  if(pNext(p)==NULL)
1078  {
1079    c=nInvers(pGetCoeff(p));
1080    pSetCoeff(p,nInit(1));
1081  }
1082  else
1083  {
1084    h = nInit(1);
1085    while (p!=NULL)
1086    {
1087      nNormalize(pGetCoeff(p));
1088      d=nLcm(h,pGetCoeff(p),r);
1089      nDelete(&h);
1090      h=d;
1091      pIter(p);
1092    }
1093    c=h;
1094    /* contains the 1/lcm of all denominators */
1095    if(!nIsOne(h))
1096    {
1097      p = ph;
1098      while (p!=NULL)
1099      {
1100        /* should be:
1101        * number hh;
1102        * nGetDenom(p->coef,&hh);
1103        * nMult(&h,&hh,&d);
1104        * nNormalize(d);
1105        * nDelete(&hh);
1106        * nMult(d,p->coef,&hh);
1107        * nDelete(&d);
1108        * nDelete(&(p->coef));
1109        * p->coef =hh;
1110        */
1111        d=nMult(h,pGetCoeff(p));
1112        nNormalize(d);
1113        pSetCoeff(p,d);
1114        pIter(p);
1115      }
1116      if (nGetChar()==1)
1117      {
1118        loop
1119        {
1120          h = nInit(1);
1121          p=ph;
1122          while (p!=NULL)
1123          {
1124            d=nLcm(h,pGetCoeff(p),r);
1125            nDelete(&h);
1126            h=d;
1127            pIter(p);
1128          }
1129          /* contains the 1/lcm of all denominators */
1130          if(!nIsOne(h))
1131          {
1132            p = ph;
1133            while (p!=NULL)
1134            {
1135              /* should be:
1136              * number hh;
1137              * nGetDenom(p->coef,&hh);
1138              * nMult(&h,&hh,&d);
1139              * nNormalize(d);
1140              * nDelete(&hh);
1141              * nMult(d,p->coef,&hh);
1142              * nDelete(&d);
1143              * nDelete(&(p->coef));
1144              * p->coef =hh;
1145              */
1146              d=nMult(h,pGetCoeff(p));
1147              nNormalize(d);
1148              pSetCoeff(p,d);
1149              pIter(p);
1150            }
1151            number t=nMult(c,h);
1152            nDelete(&c);
1153            c=t;
1154          }
1155          else
1156          {
1157            break;
1158          }
1159          nDelete(&h);
1160        }
1161      }
1162    }
1163  }
1164}
1165
1166number p_GetAllDenom(poly ph, const ring r)
1167{
1168  number d=n_Init(1,r);
1169  poly p = ph;
1170
1171  while (p!=NULL)
1172  {
1173    number h=n_GetDenom(pGetCoeff(p),r);
1174    if (!n_IsOne(h,r))
1175    {
1176      number dd=n_Mult(d,h,r);
1177      n_Delete(&d,r);
1178      d=dd;
1179    }
1180    n_Delete(&h,r);
1181    pIter(p);
1182  }
1183  return d;
1184}
1185
1186/*2
1187*tests if p is homogeneous with respect to the actual weigths
1188*/
1189BOOLEAN pIsHomogeneous (poly p)
1190{
1191  poly qp=p;
1192  int o;
1193
1194  if ((p == NULL) || (pNext(p) == NULL)) return TRUE;
1195  pFDegProc d;
1196  if (pLexOrder && (currRing->order[0]==ringorder_lp))
1197    d=p_Totaldegree;
1198  else 
1199    d=pFDeg;
1200  o = d(p,currRing);
1201  do
1202  {
1203    if (d(qp,currRing) != o) return FALSE;
1204    pIter(qp);
1205  }
1206  while (qp != NULL);
1207  return TRUE;
1208}
1209
1210// orders monoms of poly using merge sort (ususally faster than
1211// insertion sort). ASSUMES that pSetm was performed on monoms
1212poly pOrdPolyMerge(poly p)
1213{
1214  poly qq,pp,result=NULL;
1215
1216  if (p == NULL) return NULL;
1217
1218  loop
1219  {
1220    qq = p;
1221    loop
1222    {
1223      if (pNext(p) == NULL)
1224      {
1225        result=pAdd(result, qq);
1226        pTest(result);
1227        return result;
1228      }
1229      if (pLmCmp(p,pNext(p)) != 1)
1230      {
1231        pp = p;
1232        pIter(p);
1233        pNext(pp) = NULL;
1234        result = pAdd(result, qq);
1235        break;
1236      }
1237      pIter(p);
1238    }
1239  }
1240}
1241
1242// orders monoms of poly using insertion sort, performs pSetm on each monom
1243poly pOrdPolyInsertSetm(poly p)
1244{
1245  poly qq,result = NULL;
1246
1247#if 0
1248  while (p != NULL)
1249  {
1250    qq = p;
1251    pIter(p);
1252    qq->next = NULL;
1253    pSetm(qq);
1254    result = pAdd(result,qq);
1255    pTest(result);
1256  }
1257#else
1258  while (p != NULL)
1259  {
1260    qq = p;
1261    pIter(p);
1262    qq->next = result;
1263    result = qq;
1264    pSetm(qq);
1265  }
1266  p = result;
1267  result = NULL;
1268  while (p != NULL)
1269  {
1270    qq = p;
1271    pIter(p);
1272    qq->next = NULL;
1273    result = pAdd(result, qq);
1274  }
1275  pTest(result);
1276#endif
1277  return result;
1278}
1279
1280/*2
1281*returns a re-ordered copy of a polynomial, with permutation of the variables
1282*/
1283poly pPermPoly (poly p, int * perm, const ring oldRing, nMapFunc nMap,
1284   int *par_perm, int OldPar)
1285{
1286  int OldpVariables = oldRing->N;
1287  poly result = NULL;
1288  poly result_last = NULL;
1289  poly aq=NULL; /* the map coefficient */
1290  poly qq; /* the mapped monomial */
1291
1292  while (p != NULL)
1293  {
1294    if ((OldPar==0)||(rField_is_GF(oldRing)))
1295    {
1296      qq = pInit();
1297      number n=nMap(pGetCoeff(p));
1298      if ((currRing->minpoly!=NULL)
1299      && ((rField_is_Zp_a()) || (rField_is_Q_a())))
1300      {
1301        nNormalize(n);
1302      }
1303      pGetCoeff(qq)=n;
1304    // coef may be zero:  pTest(qq);
1305    }
1306    else
1307    {
1308      qq=pOne();
1309      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar, oldRing);
1310      if ((currRing->minpoly!=NULL)
1311      && ((rField_is_Zp_a()) || (rField_is_Q_a())))
1312      {
1313        poly tmp=aq;
1314        while (tmp!=NULL)
1315        {
1316          number n=pGetCoeff(tmp);
1317          nNormalize(n);
1318          pGetCoeff(tmp)=n;
1319          pIter(tmp);
1320        }
1321      }
1322      pTest(aq);
1323    }
1324    if (rRing_has_Comp(currRing)) pSetComp(qq, p_GetComp(p,oldRing));
1325    if (nIsZero(pGetCoeff(qq)))
1326    {
1327      pLmDelete(&qq);
1328    }
1329    else
1330    {
1331      int i;
1332      int mapped_to_par=0;
1333      for(i=1; i<=OldpVariables; i++)
1334      {
1335        int e=p_GetExp(p,i,oldRing);
1336        if (e!=0)
1337        {
1338          if (perm==NULL)
1339          {
1340            pSetExp(qq,i, e);
1341          }
1342          else if (perm[i]>0)
1343            pAddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/);
1344          else if (perm[i]<0)
1345          {
1346            if (rField_is_GF())
1347            {
1348              number c=pGetCoeff(qq);
1349              number ee=nfPar(1);
1350              number eee;nfPower(ee,e,&eee); //nfDelete(ee,currRing);
1351              ee=nfMult(c,eee);
1352              //nfDelete(c,currRing);nfDelete(eee,currRing);
1353              pSetCoeff0(qq,ee);
1354            }
1355            else
1356            {
1357              lnumber c=(lnumber)pGetCoeff(qq);
1358              if (c->z->next==NULL)
1359                napAddExp(c->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
1360              else /* more difficult: we have really to multiply: */
1361              {
1362                lnumber mmc=(lnumber)naInit(1,currRing);
1363                napSetExp(mmc->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
1364                napSetm(mmc->z);
1365                pGetCoeff(qq)=naMult((number)c,(number)mmc);
1366                nDelete((number *)&c);
1367                nDelete((number *)&mmc); 
1368              }
1369              mapped_to_par=1;
1370            }
1371          }
1372          else
1373          {
1374            /* this variable maps to 0 !*/
1375            pLmDelete(&qq);
1376            break;
1377          }
1378        }
1379      }
1380      if (mapped_to_par
1381      && (currRing->minpoly!=NULL))
1382      {
1383        number n=pGetCoeff(qq);
1384        nNormalize(n);
1385        pGetCoeff(qq)=n;
1386      }
1387    }
1388    pIter(p);
1389#if 1
1390    if (qq!=NULL)
1391    {
1392      pSetm(qq);
1393      pTest(aq);
1394      pTest(qq);
1395      if (aq!=NULL) qq=pMult(aq,qq);
1396      aq = qq;
1397      while (pNext(aq) != NULL) pIter(aq);
1398      if (result_last==NULL)
1399      {
1400        result=qq;
1401      }
1402      else
1403      {
1404        pNext(result_last)=qq;
1405      }
1406      result_last=aq;
1407      aq = NULL;
1408    }
1409    else if (aq!=NULL)
1410    {
1411      pDelete(&aq);
1412    }
1413  }
1414  result=pSortAdd(result);
1415#else
1416  //  if (qq!=NULL)
1417  //  {
1418  //    pSetm(qq);
1419  //    pTest(qq);
1420  //    pTest(aq);
1421  //    if (aq!=NULL) qq=pMult(aq,qq);
1422  //    aq = qq;
1423  //    while (pNext(aq) != NULL) pIter(aq);
1424  //    pNext(aq) = result;
1425  //    aq = NULL;
1426  //    result = qq;
1427  //  }
1428  //  else if (aq!=NULL)
1429  //  {
1430  //    pDelete(&aq);
1431  //  }
1432  //}
1433  //p = result;
1434  //result = NULL;
1435  //while (p != NULL)
1436  //{
1437  //  qq = p;
1438  //  pIter(p);
1439  //  qq->next = NULL;
1440  //  result = pAdd(result, qq);
1441  //}
1442#endif
1443  pTest(result);
1444  return result;
1445}
1446
1447#if 0
1448/*2
1449*returns a re-ordered copy of a polynomial, with permutation of the variables
1450*/
1451poly p_PermPoly (poly p, int * perm, ring oldRing,
1452   int *par_perm, int OldPar, ring newRing)
1453{
1454  int OldpVariables = oldRing->N;
1455  poly result = NULL;
1456  poly result_last = NULL;
1457  poly aq=NULL; /* the map coefficient */
1458  poly qq; /* the mapped monomial */
1459
1460  while (p != NULL)
1461  {
1462    if (OldPar==0)
1463    {
1464      qq = pInit();
1465      number n=newRing->cf->nMap(pGetCoeff(p));
1466      if ((newRing->minpoly!=NULL)
1467      && ((rField_is_Zp_a(newRing)) || (rField_is_Q_a(newRing))))
1468      {
1469        newRing->cf->nNormalize(n);
1470      }
1471      pGetCoeff(qq)=n;
1472    // coef may be zero:  pTest(qq);
1473    }
1474    else
1475    {
1476      qq=p_ISet(1, newRing);
1477      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar, oldRing);
1478      if ((newRing->minpoly!=NULL)
1479      && ((rField_is_Zp_a(newRing)) || (rField_is_Q_a(newRing))))
1480      {
1481        poly tmp=aq;
1482        while (tmp!=NULL)
1483        {
1484          number n=pGetCoeff(tmp);
1485          newRing->cf->nNormalize(n);
1486          pGetCoeff(tmp)=n;
1487          pIter(tmp);
1488        }
1489      }
1490      //pTest(aq);
1491    }
1492    p_SetComp(qq, p_GetComp(p,oldRing), newRing);
1493    if (newRing->cf->nIsZero(pGetCoeff(qq)))
1494    {
1495      p_DeleteLm(&qq, newRing);
1496    }
1497    else
1498    {
1499      int i;
1500      int mapped_to_par=0;
1501      for(i=1; i<=OldpVariables; i++)
1502      {
1503        int e=p_GetExp(p,i,oldRing);
1504        if (e!=0)
1505        {
1506          if (perm==NULL)
1507          {
1508            p_SetExp(qq,i, e, newRing);
1509          }
1510          else if (perm[i]>0)
1511            p_AddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/, newRing);
1512          else if (perm[i]<0)
1513          {
1514            lnumber c=(lnumber)pGetCoeff(qq);
1515            napAddExp(c->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
1516            mapped_to_par=1;
1517          }
1518          else
1519          {
1520            /* this variable maps to 0 !*/
1521            p_DeleteLm(&qq, newRing);
1522            break;
1523          }
1524        }
1525      }
1526      if (mapped_to_par
1527      && (newRing->minpoly!=NULL))
1528      {
1529        number n=pGetCoeff(qq);
1530        newRing->cf->nNormalize(n);
1531        pGetCoeff(qq)=n;
1532      }
1533    }
1534    pIter(p);
1535    if (qq!=NULL)
1536    {
1537      p_Setm(qq, newRing);
1538      //pTest(aq);
1539      //pTest(qq);
1540      if (aq!=NULL) qq=pMult(aq,qq);
1541      aq = qq;
1542      while (pNext(aq) != NULL) pIter(aq);
1543      if (result_last==NULL)
1544      {
1545        result=qq;
1546      }
1547      else
1548      {
1549        pNext(result_last)=qq;
1550      }
1551      result_last=aq;
1552      aq = NULL;
1553    }
1554    else if (aq!=NULL)
1555    {
1556      p_Delete(&aq, newRing);
1557    }
1558  }
1559  result=pOrdPolyMerge(result);
1560  //pTest(result);
1561  return result;
1562}
1563#endif
1564
1565poly ppJet(poly p, int m)
1566{
1567  poly r=NULL;
1568  poly t=NULL;
1569
1570  while (p!=NULL)
1571  {
1572    if (p_Totaldegree(p,currRing)<=m)
1573    {
1574      if (r==NULL)
1575        r=pHead(p);
1576      else
1577      if (t==NULL)
1578      {
1579        pNext(r)=pHead(p);
1580        t=pNext(r);
1581      }
1582      else
1583      {
1584        pNext(t)=pHead(p);
1585        pIter(t);
1586      }
1587    }
1588    pIter(p);
1589  }
1590  return r;
1591}
1592
1593poly pJet(poly p, int m)
1594{
1595  poly t=NULL;
1596
1597  while((p!=NULL) && (p_Totaldegree(p,currRing)>m)) pLmDelete(&p);
1598  if (p==NULL) return NULL;
1599  poly r=p;
1600  while (pNext(p)!=NULL)
1601  {
1602    if (p_Totaldegree(pNext(p),currRing)>m)
1603    {
1604      pLmDelete(&pNext(p));
1605    }
1606    else
1607      pIter(p);
1608  }
1609  return r;
1610}
1611
1612poly ppJetW(poly p, int m, short *w)
1613{
1614  poly r=NULL;
1615  poly t=NULL;
1616  while (p!=NULL)
1617  {
1618    if (totaldegreeWecart_IV(p,currRing,w)<=m)
1619    {
1620      if (r==NULL)
1621        r=pHead(p);
1622      else
1623      if (t==NULL)
1624      {
1625        pNext(r)=pHead(p);
1626        t=pNext(r);
1627      }
1628      else
1629      {
1630        pNext(t)=pHead(p);
1631        pIter(t);
1632      }
1633    }
1634    pIter(p);
1635  }
1636  return r;
1637}
1638
1639poly pJetW(poly p, int m, short *w)
1640{
1641  poly t=NULL;
1642  while((p!=NULL) && (totaldegreeWecart_IV(p,currRing,w)>m)) pLmDelete(&p);
1643  if (p==NULL) return NULL;
1644  poly r=p;
1645  while (pNext(p)!=NULL)
1646  {
1647    if (totaldegreeWecart_IV(pNext(p),currRing,w)>m)
1648    {
1649      pLmDelete(&pNext(p));
1650    }
1651    else
1652      pIter(p);
1653  }
1654  return r;
1655}
1656
1657int pMinDeg(poly p,intvec *w)
1658{
1659  if(p==NULL)
1660    return -1;
1661  int d=-1;
1662  while(p!=NULL)
1663  {
1664    int d0=0;
1665    for(int j=0;j<pVariables;j++)
1666      if(w==NULL||j>=w->length())
1667        d0+=pGetExp(p,j+1);
1668      else
1669        d0+=(*w)[j]*pGetExp(p,j+1);
1670    if(d0<d||d==-1)
1671      d=d0;
1672    pIter(p);
1673  }
1674  return d;
1675}
1676
1677poly pSeries(int n,poly p,poly u, intvec *w)
1678{
1679  short *ww=iv2array(w);
1680  if(p!=NULL)
1681  {
1682    if(u==NULL)
1683      p=pJetW(p,n,ww);
1684    else
1685      p=pJetW(pMult(p,pInvers(n-pMinDeg(p,w),u,w)),n,ww);
1686  }
1687  omFreeSize((ADDRESS)ww,(pVariables+1)*sizeof(short));
1688  return p;
1689}
1690
1691poly pInvers(int n,poly u,intvec *w)
1692{
1693  short *ww=iv2array(w);
1694  if(n<0)
1695    return NULL;
1696  number u0=nInvers(pGetCoeff(u));
1697  poly v=pNSet(u0);
1698  if(n==0)
1699    return v;
1700  poly u1=pJetW(pSub(pOne(),pMult_nn(u,u0)),n,ww);
1701  if(u1==NULL)
1702    return v;
1703  poly v1=pMult_nn(pCopy(u1),u0);
1704  v=pAdd(v,pCopy(v1));
1705  for(int i=n/pMinDeg(u1,w);i>1;i--)
1706  {
1707    v1=pJetW(pMult(v1,pCopy(u1)),n,ww);
1708    v=pAdd(v,pCopy(v1));
1709  }
1710  pDelete(&u1);
1711  pDelete(&v1);
1712  omFreeSize((ADDRESS)ww,(pVariables+1)*sizeof(short));
1713  return v;
1714}
1715
1716long pDegW(poly p, const short *w)
1717{
1718  long r=-LONG_MAX;
1719
1720  while (p!=NULL)
1721  {
1722    long t=totaldegreeWecart_IV(p,currRing,w);
1723    if (t>r) r=t;
1724    pIter(p);
1725  }
1726  return r;
1727}
1728
1729/*-----------type conversions ----------------------------*/
1730/*2
1731* input: a set of polys (len elements: p[0]..p[len-1])
1732* output: a vector
1733* p will not be changed
1734*/
1735poly  pPolys2Vec(polyset p, int len)
1736{
1737  poly v=NULL;
1738  poly h;
1739  int i;
1740
1741  for (i=len-1; i>=0; i--)
1742  {
1743    if (p[i])
1744    {
1745      h=pCopy(p[i]);
1746      pSetCompP(h,i+1);
1747      v=pAdd(v,h);
1748    }
1749  }
1750 return v;
1751}
1752
1753/*2
1754* convert a vector to a set of polys,
1755* allocates the polyset, (entries 0..(*len)-1)
1756* the vector will not be changed
1757*/
1758void  pVec2Polys(poly v, polyset *p, int *len)
1759{
1760  poly h;
1761  int k;
1762
1763  *len=pMaxComp(v);
1764  if (*len==0) *len=1;
1765  *p=(polyset)omAlloc0((*len)*sizeof(poly));
1766  while (v!=NULL)
1767  {
1768    h=pHead(v);
1769    k=pGetComp(h);
1770    pSetComp(h,0);
1771    (*p)[k-1]=pAdd((*p)[k-1],h);
1772    pIter(v);
1773  }
1774}
1775
1776int p_Var(poly m,const ring r)
1777{
1778  if (m==NULL) return 0;
1779  if (pNext(m)!=NULL) return 0;
1780  int i,e=0;
1781  for (i=r->N; i>0; i--)
1782  {
1783    int exp=p_GetExp(m,i,r);
1784    if (exp==1)
1785    {
1786      if (e==0) e=i;
1787      else return 0;
1788    }
1789    else if (exp!=0)
1790    {
1791      return 0;
1792    }
1793  }
1794  return e;
1795}
1796
1797/*----------utilities for syzygies--------------*/
1798//BOOLEAN   pVectorHasUnitM(poly p, int * k)
1799//{
1800//  while (p!=NULL)
1801//  {
1802//    if (pLmIsConstantComp(p))
1803//    {
1804//      *k = pGetComp(p);
1805//      return TRUE;
1806//    }
1807//    else pIter(p);
1808//  }
1809//  return FALSE;
1810//}
1811
1812BOOLEAN   pVectorHasUnitB(poly p, int * k)
1813{
1814  poly q=p,qq;
1815  int i;
1816
1817  while (q!=NULL)
1818  {
1819    if (pLmIsConstantComp(q))
1820    {
1821      i = pGetComp(q);
1822      qq = p;
1823      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1824      if (qq == q)
1825      {
1826        *k = i;
1827        return TRUE;
1828      }
1829      else
1830        pIter(q);
1831    }
1832    else pIter(q);
1833  }
1834  return FALSE;
1835}
1836
1837void   pVectorHasUnit(poly p, int * k, int * len)
1838{
1839  poly q=p,qq;
1840  int i,j=0;
1841
1842  *len = 0;
1843  while (q!=NULL)
1844  {
1845    if (pLmIsConstantComp(q))
1846    {
1847      i = pGetComp(q);
1848      qq = p;
1849      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1850      if (qq == q)
1851      {
1852       j = 0;
1853       while (qq!=NULL)
1854       {
1855         if (pGetComp(qq)==i) j++;
1856        pIter(qq);
1857       }
1858       if ((*len == 0) || (j<*len))
1859       {
1860         *len = j;
1861         *k = i;
1862       }
1863      }
1864    }
1865    pIter(q);
1866  }
1867}
1868
1869/*2
1870* returns TRUE if p1 = p2
1871*/
1872BOOLEAN p_EqualPolys(poly p1,poly p2, const ring r)
1873{
1874  while ((p1 != NULL) && (p2 != NULL))
1875  {
1876    if (! p_LmEqual(p1, p2,r))
1877      return FALSE;
1878    if (! n_Equal(p_GetCoeff(p1,r), p_GetCoeff(p2,r),r ))
1879      return FALSE;
1880    pIter(p1);
1881    pIter(p2);
1882  }
1883  return (p1==p2);
1884}
1885
1886/*2
1887*returns TRUE if p1 is a skalar multiple of p2
1888*assume p1 != NULL and p2 != NULL
1889*/
1890BOOLEAN pComparePolys(poly p1,poly p2)
1891{
1892  number n,nn;
1893  int i;
1894  pAssume(p1 != NULL && p2 != NULL);
1895
1896  if (!pLmEqual(p1,p2)) //compare leading mons
1897      return FALSE;
1898  if ((pNext(p1)==NULL) && (pNext(p2)!=NULL))
1899     return FALSE;
1900  if ((pNext(p2)==NULL) && (pNext(p1)!=NULL))
1901     return FALSE;
1902  if (pLength(p1) != pLength(p2))
1903    return FALSE;
1904  n=nDiv(pGetCoeff(p1),pGetCoeff(p2));
1905  while ((p1 != NULL) /*&& (p2 != NULL)*/)
1906  {
1907    if ( ! pLmEqual(p1, p2))
1908    {
1909        nDelete(&n);
1910        return FALSE;
1911    }
1912    if (!nEqual(pGetCoeff(p1),nn=nMult(pGetCoeff(p2),n)))
1913    {
1914      nDelete(&n);
1915      nDelete(&nn);
1916      return FALSE;
1917    }
1918    nDelete(&nn);
1919    pIter(p1);
1920    pIter(p2);
1921  }
1922  nDelete(&n);
1923  return TRUE;
1924}
Note: See TracBrowser for help on using the repository browser.