source: git/kernel/polys1.cc @ 53f716

fieker-DuValspielwiese
Last change on this file since 53f716 was 53f716, checked in by Hans Schoenemann <hannes@…>, 14 years ago
pDeleteLm -> pLmDelete git-svn-id: file:///usr/local/Singular/svn/trunk@12895 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 36.1 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 "mod2.h"
14#include "options.h"
15#include "numbers.h"
16#include "ffields.h"
17#include "omalloc.h"
18#include "febase.h"
19#include "weight.h"
20#include "intvec.h"
21#include "longalg.h"
22#include "ring.h"
23#include "ideals.h"
24#include "polys.h"
25//#include "ipid.h"
26#ifdef HAVE_FACTORY
27#include "clapsing.h"
28#endif
29
30#ifdef HAVE_RATGRING
31#include "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_DeleteLm(&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}
651void pSimpleContent(poly ph,int smax)
652{
653  if(TEST_OPT_CONTENTSB) return;
654  if (ph==NULL) return;
655  if (pNext(ph)==NULL)
656  {
657    pSetCoeff(ph,nInit(1));
658    return;
659  }
660  if ((pNext(pNext(ph))==NULL)||(!rField_is_Q()))
661  {
662    return;
663  }
664  number d=pInitContent(ph);
665  if (nlSize(d)<=smax)
666  {
667    //if (TEST_OPT_PROT) PrintS("G");
668    return;
669  }
670  poly p=ph;
671  number h=d;
672  if (smax==1) smax=2;
673  while (p!=NULL)
674  {
675#if 0
676    d=nlGcd(h,pGetCoeff(p),currRing);
677    nlDelete(&h,currRing);
678    h = d;
679#else
680    nlInpGcd(h,pGetCoeff(p),currRing);
681#endif
682    if(nlSize(h)<smax)
683    {
684      //if (TEST_OPT_PROT) PrintS("g");
685      return;
686    }
687    pIter(p);
688  }
689  p = ph;
690  if (!nlGreaterZero(pGetCoeff(p))) h=nlNeg(h);
691  if(nlIsOne(h)) return;
692  //if (TEST_OPT_PROT) PrintS("c");
693  while (p!=NULL)
694  {
695#if 1
696    d = nlIntDiv(pGetCoeff(p),h);
697    pSetCoeff(p,d);
698#else
699    nlInpIntDiv(pGetCoeff(p),h,currRing);
700#endif
701    pIter(p);
702  }
703  nlDelete(&h,currRing);
704}
705
706number pInitContent(poly ph)
707// only for coefficients in Q
708#if 0
709{
710  assume(!TEST_OPT_CONTENTSB);
711  assume(ph!=NULL);
712  assume(pNext(ph)!=NULL);
713  assume(rField_is_Q());
714  if (pNext(pNext(ph))==NULL)
715  {
716    return nlGetNom(pGetCoeff(pNext(ph)),currRing);
717  }
718  poly p=ph;
719  number n1=nlGetNom(pGetCoeff(p),currRing);
720  pIter(p);
721  number n2=nlGetNom(pGetCoeff(p),currRing);
722  pIter(p);
723  number d;
724  number t;
725  loop
726  {
727    nlNormalize(pGetCoeff(p));
728    t=nlGetNom(pGetCoeff(p),currRing);
729    if (nlGreaterZero(t))
730      d=nlAdd(n1,t);
731    else
732      d=nlSub(n1,t);
733    nlDelete(&t,currRing);
734    nlDelete(&n1,currRing);
735    n1=d;
736    pIter(p);
737    if (p==NULL) break;
738    nlNormalize(pGetCoeff(p));
739    t=nlGetNom(pGetCoeff(p),currRing);
740    if (nlGreaterZero(t))
741      d=nlAdd(n2,t);
742    else
743      d=nlSub(n2,t);
744    nlDelete(&t,currRing);
745    nlDelete(&n2,currRing);
746    n2=d;
747    pIter(p);
748    if (p==NULL) break;
749  }
750  d=nlGcd(n1,n2,currRing);
751  nlDelete(&n1,currRing);
752  nlDelete(&n2,currRing);
753  return d;
754}
755#else
756{
757  number d=pGetCoeff(ph);
758  if(SR_HDL(d)&SR_INT) return d;
759  int s=mpz_size1(d->z);
760  int s2=-1;
761  number d2;
762  loop
763  {
764    pIter(ph);
765    if(ph==NULL)
766    {
767      if (s2==-1) return nlCopy(d);
768      break;
769    }
770    if (SR_HDL(pGetCoeff(ph))&SR_INT)
771    {
772      s2=s;
773      d2=d;
774      s=0;
775      d=pGetCoeff(ph);
776      if (s2==0) break;
777    }
778    else
779    if (mpz_size1((pGetCoeff(ph)->z))<=s)
780    {
781      s2=s;
782      d2=d;
783      d=pGetCoeff(ph);
784      s=mpz_size1(d->z);
785    }
786  }
787  return nlGcd(d,d2,currRing);
788}
789#endif
790
791number pInitContent_a(poly ph)
792// only for coefficients in K(a) anf K(a,...)
793{
794  number d=pGetCoeff(ph);
795  int s=naParDeg(d);
796  if (s /* naParDeg(d)*/ <=1) return naCopy(d);
797  int s2=-1;
798  number d2;
799  int ss;
800  loop
801  {
802    pIter(ph);
803    if(ph==NULL)
804    {
805      if (s2==-1) return naCopy(d);
806      break;
807    }
808    if ((ss=naParDeg(pGetCoeff(ph)))<s)
809    {
810      s2=s;
811      d2=d;
812      s=ss;
813      d=pGetCoeff(ph);
814      if (s2<=1) break;
815    }
816  }
817  return naGcd(d,d2,currRing);
818}
819
820
821//void pContent(poly ph)
822//{
823//  number h,d;
824//  poly p;
825//
826//  p = ph;
827//  if(pNext(p)==NULL)
828//  {
829//    pSetCoeff(p,nInit(1));
830//  }
831//  else
832//  {
833//#ifdef PDEBUG
834//    if (!pTest(p)) return;
835//#endif
836//    nNormalize(pGetCoeff(p));
837//    if(!nGreaterZero(pGetCoeff(ph)))
838//    {
839//      ph = pNeg(ph);
840//      nNormalize(pGetCoeff(p));
841//    }
842//    h=pGetCoeff(p);
843//    pIter(p);
844//    while (p!=NULL)
845//    {
846//      nNormalize(pGetCoeff(p));
847//      if (nGreater(h,pGetCoeff(p))) h=pGetCoeff(p);
848//      pIter(p);
849//    }
850//    h=nCopy(h);
851//    p=ph;
852//    while (p!=NULL)
853//    {
854//      d=nGcd(h,pGetCoeff(p));
855//      nDelete(&h);
856//      h = d;
857//      if(nIsOne(h))
858//      {
859//        break;
860//      }
861//      pIter(p);
862//    }
863//    p = ph;
864//    //number tmp;
865//    if(!nIsOne(h))
866//    {
867//      while (p!=NULL)
868//      {
869//        d = nIntDiv(pGetCoeff(p),h);
870//        pSetCoeff(p,d);
871//        pIter(p);
872//      }
873//    }
874//    nDelete(&h);
875//#ifdef HAVE_FACTORY
876//    if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
877//    {
878//      pTest(ph);
879//      singclap_divide_content(ph);
880//      pTest(ph);
881//    }
882//#endif
883//  }
884//}
885#if 0
886void p_Content(poly ph, ring r)
887{
888  number h,d;
889  poly p;
890
891  if(pNext(ph)==NULL)
892  {
893    pSetCoeff(ph,n_Init(1,r));
894  }
895  else
896  {
897    n_Normalize(pGetCoeff(ph),r);
898    if(!n_GreaterZero(pGetCoeff(ph),r)) ph = p_Neg(ph,r);
899    h=n_Copy(pGetCoeff(ph),r);
900    p = pNext(ph);
901    while (p!=NULL)
902    {
903      n_Normalize(pGetCoeff(p),r);
904      d=n_Gcd(h,pGetCoeff(p),r);
905      n_Delete(&h,r);
906      h = d;
907      if(n_IsOne(h,r))
908      {
909        break;
910      }
911      pIter(p);
912    }
913    p = ph;
914    //number tmp;
915    if(!n_IsOne(h,r))
916    {
917      while (p!=NULL)
918      {
919        //d = nDiv(pGetCoeff(p),h);
920        //tmp = nIntDiv(pGetCoeff(p),h);
921        //if (!nEqual(d,tmp))
922        //{
923        //  StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
924        //  nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
925        //  nWrite(tmp);Print(StringAppendS("\n"));
926        //}
927        //nDelete(&tmp);
928        d = n_IntDiv(pGetCoeff(p),h,r);
929        p_SetCoeff(p,d,r);
930        pIter(p);
931      }
932    }
933    n_Delete(&h,r);
934#ifdef HAVE_FACTORY
935    //if ( (n_GetChar(r) == 1) || (n_GetChar(r) < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
936    //{
937    //  singclap_divide_content(ph);
938    //  if(!n_GreaterZero(pGetCoeff(ph),r)) ph = p_Neg(ph,r);
939    //}
940#endif
941  }
942}
943#endif
944
945poly p_Cleardenom(poly ph, const ring r)
946{
947  poly start=ph;
948  number d, h;
949  poly p;
950
951#ifdef HAVE_RINGS
952  if (rField_is_Ring(r))
953  {
954    p_Content(ph,r);
955    return start;
956  }
957#endif
958  if (rField_is_Zp(r) && TEST_OPT_INTSTRATEGY) return start;
959  p = ph;
960  if(pNext(p)==NULL)
961  {
962    if (TEST_OPT_CONTENTSB)
963    {
964      number n=nGetDenom(pGetCoeff(p));
965      if (!nIsOne(n))
966      {
967        number nn=nMult(pGetCoeff(p),n);
968        nNormalize(nn);
969        pSetCoeff(p,nn);
970      }
971      nDelete(&n);
972    }
973    else
974      pSetCoeff(p,nInit(1));
975  }
976  else
977  {
978    h = nInit(1);
979    while (p!=NULL)
980    {
981      nNormalize(pGetCoeff(p));
982      d=nLcm(h,pGetCoeff(p),currRing);
983      nDelete(&h);
984      h=d;
985      pIter(p);
986    }
987    /* contains the 1/lcm of all denominators */
988    if(!nIsOne(h))
989    {
990      p = ph;
991      while (p!=NULL)
992      {
993        /* should be:
994        * number hh;
995        * nGetDenom(p->coef,&hh);
996        * nMult(&h,&hh,&d);
997        * nNormalize(d);
998        * nDelete(&hh);
999        * nMult(d,p->coef,&hh);
1000        * nDelete(&d);
1001        * nDelete(&(p->coef));
1002        * p->coef =hh;
1003        */
1004        d=nMult(h,pGetCoeff(p));
1005        nNormalize(d);
1006        pSetCoeff(p,d);
1007        pIter(p);
1008      }
1009      nDelete(&h);
1010      if (nGetChar()==1)
1011      {
1012        loop
1013        {
1014          h = nInit(1);
1015          p=ph;
1016          while (p!=NULL)
1017          {
1018            d=nLcm(h,pGetCoeff(p),currRing);
1019            nDelete(&h);
1020            h=d;
1021            pIter(p);
1022          }
1023          /* contains the 1/lcm of all denominators */
1024          if(!nIsOne(h))
1025          {
1026            p = ph;
1027            while (p!=NULL)
1028            {
1029              /* should be:
1030              * number hh;
1031              * nGetDenom(p->coef,&hh);
1032              * nMult(&h,&hh,&d);
1033              * nNormalize(d);
1034              * nDelete(&hh);
1035              * nMult(d,p->coef,&hh);
1036              * nDelete(&d);
1037              * nDelete(&(p->coef));
1038              * p->coef =hh;
1039              */
1040              d=nMult(h,pGetCoeff(p));
1041              nNormalize(d);
1042              pSetCoeff(p,d);
1043              pIter(p);
1044            }
1045            nDelete(&h);
1046          }
1047          else
1048          {
1049            nDelete(&h);
1050            break;
1051          }
1052        }
1053      }
1054    }
1055    if (h!=NULL) nDelete(&h);
1056 
1057    p_Content(ph,r);
1058#ifdef HAVE_RATGRING
1059    if (rIsRatGRing(r))
1060    {
1061      /* quick unit detection in the rational case is done in gr_nc_bba */
1062      pContentRat(ph);
1063      start=ph;
1064    }
1065#endif
1066  }
1067  return start;
1068}
1069
1070void p_Cleardenom_n(poly ph,const ring r,number &c)
1071{
1072  number d, h;
1073  poly p;
1074
1075  p = ph;
1076  if(pNext(p)==NULL)
1077  {
1078    c=nInvers(pGetCoeff(p));
1079    pSetCoeff(p,nInit(1));
1080  }
1081  else
1082  {
1083    h = nInit(1);
1084    while (p!=NULL)
1085    {
1086      nNormalize(pGetCoeff(p));
1087      d=nLcm(h,pGetCoeff(p),r);
1088      nDelete(&h);
1089      h=d;
1090      pIter(p);
1091    }
1092    c=h;
1093    /* contains the 1/lcm of all denominators */
1094    if(!nIsOne(h))
1095    {
1096      p = ph;
1097      while (p!=NULL)
1098      {
1099        /* should be:
1100        * number hh;
1101        * nGetDenom(p->coef,&hh);
1102        * nMult(&h,&hh,&d);
1103        * nNormalize(d);
1104        * nDelete(&hh);
1105        * nMult(d,p->coef,&hh);
1106        * nDelete(&d);
1107        * nDelete(&(p->coef));
1108        * p->coef =hh;
1109        */
1110        d=nMult(h,pGetCoeff(p));
1111        nNormalize(d);
1112        pSetCoeff(p,d);
1113        pIter(p);
1114      }
1115      if (nGetChar()==1)
1116      {
1117        loop
1118        {
1119          h = nInit(1);
1120          p=ph;
1121          while (p!=NULL)
1122          {
1123            d=nLcm(h,pGetCoeff(p),r);
1124            nDelete(&h);
1125            h=d;
1126            pIter(p);
1127          }
1128          /* contains the 1/lcm of all denominators */
1129          if(!nIsOne(h))
1130          {
1131            p = ph;
1132            while (p!=NULL)
1133            {
1134              /* should be:
1135              * number hh;
1136              * nGetDenom(p->coef,&hh);
1137              * nMult(&h,&hh,&d);
1138              * nNormalize(d);
1139              * nDelete(&hh);
1140              * nMult(d,p->coef,&hh);
1141              * nDelete(&d);
1142              * nDelete(&(p->coef));
1143              * p->coef =hh;
1144              */
1145              d=nMult(h,pGetCoeff(p));
1146              nNormalize(d);
1147              pSetCoeff(p,d);
1148              pIter(p);
1149            }
1150            number t=nMult(c,h);
1151            nDelete(&c);
1152            c=t;
1153          }
1154          else
1155          {
1156            break;
1157          }
1158          nDelete(&h);
1159        }
1160      }
1161    }
1162  }
1163}
1164
1165number p_GetAllDenom(poly ph, const ring r)
1166{
1167  number d=n_Init(1,r);
1168  poly p = ph;
1169
1170  while (p!=NULL)
1171  {
1172    number h=n_GetDenom(pGetCoeff(p),r);
1173    if (!n_IsOne(h,r))
1174    {
1175      number dd=n_Mult(d,h,r);
1176      n_Delete(&d,r);
1177      d=dd;
1178    }
1179    n_Delete(&h,r);
1180    pIter(p);
1181  }
1182  return d;
1183}
1184
1185/*2
1186*tests if p is homogeneous with respect to the actual weigths
1187*/
1188BOOLEAN pIsHomogeneous (poly p)
1189{
1190  poly qp=p;
1191  int o;
1192
1193  if ((p == NULL) || (pNext(p) == NULL)) return TRUE;
1194  pFDegProc d;
1195  if (pLexOrder && (currRing->order[0]==ringorder_lp))
1196    d=pTotaldegree;
1197  else 
1198    d=pFDeg;
1199  o = d(p,currRing);
1200  do
1201  {
1202    if (d(qp,currRing) != o) return FALSE;
1203    pIter(qp);
1204  }
1205  while (qp != NULL);
1206  return TRUE;
1207}
1208
1209// orders monoms of poly using merge sort (ususally faster than
1210// insertion sort). ASSUMES that pSetm was performed on monoms
1211poly pOrdPolyMerge(poly p)
1212{
1213  poly qq,pp,result=NULL;
1214
1215  if (p == NULL) return NULL;
1216
1217  loop
1218  {
1219    qq = p;
1220    loop
1221    {
1222      if (pNext(p) == NULL)
1223      {
1224        result=pAdd(result, qq);
1225        pTest(result);
1226        return result;
1227      }
1228      if (pLmCmp(p,pNext(p)) != 1)
1229      {
1230        pp = p;
1231        pIter(p);
1232        pNext(pp) = NULL;
1233        result = pAdd(result, qq);
1234        break;
1235      }
1236      pIter(p);
1237    }
1238  }
1239}
1240
1241// orders monoms of poly using insertion sort, performs pSetm on each monom
1242poly pOrdPolyInsertSetm(poly p)
1243{
1244  poly qq,result = NULL;
1245
1246#if 0
1247  while (p != NULL)
1248  {
1249    qq = p;
1250    pIter(p);
1251    qq->next = NULL;
1252    pSetm(qq);
1253    result = pAdd(result,qq);
1254    pTest(result);
1255  }
1256#else
1257  while (p != NULL)
1258  {
1259    qq = p;
1260    pIter(p);
1261    qq->next = result;
1262    result = qq;
1263    pSetm(qq);
1264  }
1265  p = result;
1266  result = NULL;
1267  while (p != NULL)
1268  {
1269    qq = p;
1270    pIter(p);
1271    qq->next = NULL;
1272    result = pAdd(result, qq);
1273  }
1274  pTest(result);
1275#endif
1276  return result;
1277}
1278
1279/*2
1280*returns a re-ordered copy of a polynomial, with permutation of the variables
1281*/
1282poly pPermPoly (poly p, int * perm, const ring oldRing, nMapFunc nMap,
1283   int *par_perm, int OldPar)
1284{
1285  int OldpVariables = oldRing->N;
1286  poly result = NULL;
1287  poly result_last = NULL;
1288  poly aq=NULL; /* the map coefficient */
1289  poly qq; /* the mapped monomial */
1290
1291  while (p != NULL)
1292  {
1293    if ((OldPar==0)||(rField_is_GF(oldRing)))
1294    {
1295      qq = pInit();
1296      number n=nMap(pGetCoeff(p));
1297      if ((currRing->minpoly!=NULL)
1298      && ((rField_is_Zp_a()) || (rField_is_Q_a())))
1299      {
1300        nNormalize(n);
1301      }
1302      pGetCoeff(qq)=n;
1303    // coef may be zero:  pTest(qq);
1304    }
1305    else
1306    {
1307      qq=pOne();
1308      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar, oldRing);
1309      if ((currRing->minpoly!=NULL)
1310      && ((rField_is_Zp_a()) || (rField_is_Q_a())))
1311      {
1312        poly tmp=aq;
1313        while (tmp!=NULL)
1314        {
1315          number n=pGetCoeff(tmp);
1316          nNormalize(n);
1317          pGetCoeff(tmp)=n;
1318          pIter(tmp);
1319        }
1320      }
1321      pTest(aq);
1322    }
1323    if (rRing_has_Comp(currRing)) pSetComp(qq, p_GetComp(p,oldRing));
1324    if (nIsZero(pGetCoeff(qq)))
1325    {
1326      pLmDelete(&qq);
1327    }
1328    else
1329    {
1330      int i;
1331      int mapped_to_par=0;
1332      for(i=1; i<=OldpVariables; i++)
1333      {
1334        int e=p_GetExp(p,i,oldRing);
1335        if (e!=0)
1336        {
1337          if (perm==NULL)
1338          {
1339            pSetExp(qq,i, e);
1340          }
1341          else if (perm[i]>0)
1342            pAddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/);
1343          else if (perm[i]<0)
1344          {
1345            if (rField_is_GF())
1346            {
1347              number c=pGetCoeff(qq);
1348              number ee=nfPar(1);
1349              number eee;nfPower(ee,e,&eee); //nfDelete(ee,currRing);
1350              ee=nfMult(c,eee);
1351              //nfDelete(c,currRing);nfDelete(eee,currRing);
1352              pSetCoeff0(qq,ee);
1353            }
1354            else
1355            {
1356              lnumber c=(lnumber)pGetCoeff(qq);
1357              if (c->z->next==NULL)
1358                napAddExp(c->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
1359              else /* more difficult: we have really to multiply: */
1360              {
1361                lnumber mmc=(lnumber)naInit(1,currRing);
1362                napSetExp(mmc->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
1363                napSetm(mmc->z);
1364                pGetCoeff(qq)=naMult((number)c,(number)mmc);
1365                nDelete((number *)&c);
1366                nDelete((number *)&mmc); 
1367              }
1368              mapped_to_par=1;
1369            }
1370          }
1371          else
1372          {
1373            /* this variable maps to 0 !*/
1374            pLmDelete(&qq);
1375            break;
1376          }
1377        }
1378      }
1379      if (mapped_to_par
1380      && (currRing->minpoly!=NULL))
1381      {
1382        number n=pGetCoeff(qq);
1383        nNormalize(n);
1384        pGetCoeff(qq)=n;
1385      }
1386    }
1387    pIter(p);
1388#if 1
1389    if (qq!=NULL)
1390    {
1391      pSetm(qq);
1392      pTest(aq);
1393      pTest(qq);
1394      if (aq!=NULL) qq=pMult(aq,qq);
1395      aq = qq;
1396      while (pNext(aq) != NULL) pIter(aq);
1397      if (result_last==NULL)
1398      {
1399        result=qq;
1400      }
1401      else
1402      {
1403        pNext(result_last)=qq;
1404      }
1405      result_last=aq;
1406      aq = NULL;
1407    }
1408    else if (aq!=NULL)
1409    {
1410      pDelete(&aq);
1411    }
1412  }
1413  result=pSortAdd(result);
1414#else
1415  //  if (qq!=NULL)
1416  //  {
1417  //    pSetm(qq);
1418  //    pTest(qq);
1419  //    pTest(aq);
1420  //    if (aq!=NULL) qq=pMult(aq,qq);
1421  //    aq = qq;
1422  //    while (pNext(aq) != NULL) pIter(aq);
1423  //    pNext(aq) = result;
1424  //    aq = NULL;
1425  //    result = qq;
1426  //  }
1427  //  else if (aq!=NULL)
1428  //  {
1429  //    pDelete(&aq);
1430  //  }
1431  //}
1432  //p = result;
1433  //result = NULL;
1434  //while (p != NULL)
1435  //{
1436  //  qq = p;
1437  //  pIter(p);
1438  //  qq->next = NULL;
1439  //  result = pAdd(result, qq);
1440  //}
1441#endif
1442  pTest(result);
1443  return result;
1444}
1445
1446#if 0
1447/*2
1448*returns a re-ordered copy of a polynomial, with permutation of the variables
1449*/
1450poly p_PermPoly (poly p, int * perm, ring oldRing,
1451   int *par_perm, int OldPar, ring newRing)
1452{
1453  int OldpVariables = oldRing->N;
1454  poly result = NULL;
1455  poly result_last = NULL;
1456  poly aq=NULL; /* the map coefficient */
1457  poly qq; /* the mapped monomial */
1458
1459  while (p != NULL)
1460  {
1461    if (OldPar==0)
1462    {
1463      qq = pInit();
1464      number n=newRing->cf->nMap(pGetCoeff(p));
1465      if ((newRing->minpoly!=NULL)
1466      && ((rField_is_Zp_a(newRing)) || (rField_is_Q_a(newRing))))
1467      {
1468        newRing->cf->nNormalize(n);
1469      }
1470      pGetCoeff(qq)=n;
1471    // coef may be zero:  pTest(qq);
1472    }
1473    else
1474    {
1475      qq=p_ISet(1, newRing);
1476      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar, oldRing);
1477      if ((newRing->minpoly!=NULL)
1478      && ((rField_is_Zp_a(newRing)) || (rField_is_Q_a(newRing))))
1479      {
1480        poly tmp=aq;
1481        while (tmp!=NULL)
1482        {
1483          number n=pGetCoeff(tmp);
1484          newRing->cf->nNormalize(n);
1485          pGetCoeff(tmp)=n;
1486          pIter(tmp);
1487        }
1488      }
1489      //pTest(aq);
1490    }
1491    p_SetComp(qq, p_GetComp(p,oldRing), newRing);
1492    if (newRing->cf->nIsZero(pGetCoeff(qq)))
1493    {
1494      p_DeleteLm(&qq, newRing);
1495    }
1496    else
1497    {
1498      int i;
1499      int mapped_to_par=0;
1500      for(i=1; i<=OldpVariables; i++)
1501      {
1502        int e=p_GetExp(p,i,oldRing);
1503        if (e!=0)
1504        {
1505          if (perm==NULL)
1506          {
1507            p_SetExp(qq,i, e, newRing);
1508          }
1509          else if (perm[i]>0)
1510            p_AddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/, newRing);
1511          else if (perm[i]<0)
1512          {
1513            lnumber c=(lnumber)pGetCoeff(qq);
1514            napAddExp(c->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
1515            mapped_to_par=1;
1516          }
1517          else
1518          {
1519            /* this variable maps to 0 !*/
1520            p_DeleteLm(&qq, newRing);
1521            break;
1522          }
1523        }
1524      }
1525      if (mapped_to_par
1526      && (newRing->minpoly!=NULL))
1527      {
1528        number n=pGetCoeff(qq);
1529        newRing->cf->nNormalize(n);
1530        pGetCoeff(qq)=n;
1531      }
1532    }
1533    pIter(p);
1534    if (qq!=NULL)
1535    {
1536      p_Setm(qq, newRing);
1537      //pTest(aq);
1538      //pTest(qq);
1539      if (aq!=NULL) qq=pMult(aq,qq);
1540      aq = qq;
1541      while (pNext(aq) != NULL) pIter(aq);
1542      if (result_last==NULL)
1543      {
1544        result=qq;
1545      }
1546      else
1547      {
1548        pNext(result_last)=qq;
1549      }
1550      result_last=aq;
1551      aq = NULL;
1552    }
1553    else if (aq!=NULL)
1554    {
1555      p_Delete(&aq, newRing);
1556    }
1557  }
1558  result=pOrdPolyMerge(result);
1559  //pTest(result);
1560  return result;
1561}
1562#endif
1563
1564poly ppJet(poly p, int m)
1565{
1566  poly r=NULL;
1567  poly t=NULL;
1568
1569  while (p!=NULL)
1570  {
1571    if (pTotaldegree(p)<=m)
1572    {
1573      if (r==NULL)
1574        r=pHead(p);
1575      else
1576      if (t==NULL)
1577      {
1578        pNext(r)=pHead(p);
1579        t=pNext(r);
1580      }
1581      else
1582      {
1583        pNext(t)=pHead(p);
1584        pIter(t);
1585      }
1586    }
1587    pIter(p);
1588  }
1589  return r;
1590}
1591
1592poly pJet(poly p, int m)
1593{
1594  poly t=NULL;
1595
1596  while((p!=NULL) && (pTotaldegree(p)>m)) pLmDelete(&p);
1597  if (p==NULL) return NULL;
1598  poly r=p;
1599  while (pNext(p)!=NULL)
1600  {
1601    if (pTotaldegree(pNext(p))>m)
1602    {
1603      pLmDelete(&pNext(p));
1604    }
1605    else
1606      pIter(p);
1607  }
1608  return r;
1609}
1610
1611poly ppJetW(poly p, int m, short *w)
1612{
1613  poly r=NULL;
1614  poly t=NULL;
1615  while (p!=NULL)
1616  {
1617    if (totaldegreeWecart_IV(p,currRing,w)<=m)
1618    {
1619      if (r==NULL)
1620        r=pHead(p);
1621      else
1622      if (t==NULL)
1623      {
1624        pNext(r)=pHead(p);
1625        t=pNext(r);
1626      }
1627      else
1628      {
1629        pNext(t)=pHead(p);
1630        pIter(t);
1631      }
1632    }
1633    pIter(p);
1634  }
1635  return r;
1636}
1637
1638poly pJetW(poly p, int m, short *w)
1639{
1640  poly t=NULL;
1641  while((p!=NULL) && (totaldegreeWecart_IV(p,currRing,w)>m)) pLmDelete(&p);
1642  if (p==NULL) return NULL;
1643  poly r=p;
1644  while (pNext(p)!=NULL)
1645  {
1646    if (totaldegreeWecart_IV(pNext(p),currRing,w)>m)
1647    {
1648      pLmDelete(&pNext(p));
1649    }
1650    else
1651      pIter(p);
1652  }
1653  return r;
1654}
1655
1656int pMinDeg(poly p,intvec *w)
1657{
1658  if(p==NULL)
1659    return -1;
1660  int d=-1;
1661  while(p!=NULL)
1662  {
1663    int d0=0;
1664    for(int j=0;j<pVariables;j++)
1665      if(w==NULL||j>=w->length())
1666        d0+=pGetExp(p,j+1);
1667      else
1668        d0+=(*w)[j]*pGetExp(p,j+1);
1669    if(d0<d||d==-1)
1670      d=d0;
1671    pIter(p);
1672  }
1673  return d;
1674}
1675
1676poly pSeries(int n,poly p,poly u, intvec *w)
1677{
1678  short *ww=iv2array(w);
1679  if(p!=NULL)
1680  {
1681    if(u==NULL)
1682      p=pJetW(p,n,ww);
1683    else
1684      p=pJetW(pMult(p,pInvers(n-pMinDeg(p,w),u,w)),n,ww);
1685  }
1686  omFreeSize((ADDRESS)ww,(pVariables+1)*sizeof(short));
1687  return p;
1688}
1689
1690poly pInvers(int n,poly u,intvec *w)
1691{
1692  short *ww=iv2array(w);
1693  if(n<0)
1694    return NULL;
1695  number u0=nInvers(pGetCoeff(u));
1696  poly v=pNSet(u0);
1697  if(n==0)
1698    return v;
1699  poly u1=pJetW(pSub(pOne(),pMult_nn(u,u0)),n,ww);
1700  if(u1==NULL)
1701    return v;
1702  poly v1=pMult_nn(pCopy(u1),u0);
1703  v=pAdd(v,pCopy(v1));
1704  for(int i=n/pMinDeg(u1,w);i>1;i--)
1705  {
1706    v1=pJetW(pMult(v1,pCopy(u1)),n,ww);
1707    v=pAdd(v,pCopy(v1));
1708  }
1709  pDelete(&u1);
1710  pDelete(&v1);
1711  omFreeSize((ADDRESS)ww,(pVariables+1)*sizeof(short));
1712  return v;
1713}
1714
1715long pDegW(poly p, const short *w)
1716{
1717  long r=-LONG_MAX;
1718
1719  while (p!=NULL)
1720  {
1721    long t=totaldegreeWecart_IV(p,currRing,w);
1722    if (t>r) r=t;
1723    pIter(p);
1724  }
1725  return r;
1726}
1727
1728/*-----------type conversions ----------------------------*/
1729/*2
1730* input: a set of polys (len elements: p[0]..p[len-1])
1731* output: a vector
1732* p will not be changed
1733*/
1734poly  pPolys2Vec(polyset p, int len)
1735{
1736  poly v=NULL;
1737  poly h;
1738  int i;
1739
1740  for (i=len-1; i>=0; i--)
1741  {
1742    if (p[i])
1743    {
1744      h=pCopy(p[i]);
1745      pSetCompP(h,i+1);
1746      v=pAdd(v,h);
1747    }
1748  }
1749 return v;
1750}
1751
1752/*2
1753* convert a vector to a set of polys,
1754* allocates the polyset, (entries 0..(*len)-1)
1755* the vector will not be changed
1756*/
1757void  pVec2Polys(poly v, polyset *p, int *len)
1758{
1759  poly h;
1760  int k;
1761
1762  *len=pMaxComp(v);
1763  if (*len==0) *len=1;
1764  *p=(polyset)omAlloc0((*len)*sizeof(poly));
1765  while (v!=NULL)
1766  {
1767    h=pHead(v);
1768    k=pGetComp(h);
1769    pSetComp(h,0);
1770    (*p)[k-1]=pAdd((*p)[k-1],h);
1771    pIter(v);
1772  }
1773}
1774
1775int p_Var(poly m,const ring r)
1776{
1777  if (m==NULL) return 0;
1778  if (pNext(m)!=NULL) return 0;
1779  int i,e=0;
1780  for (i=r->N; i>0; i--)
1781  {
1782    int exp=p_GetExp(m,i,r);
1783    if (exp==1)
1784    {
1785      if (e==0) e=i;
1786      else return 0;
1787    }
1788    else if (exp!=0)
1789    {
1790      return 0;
1791    }
1792  }
1793  return e;
1794}
1795
1796/*----------utilities for syzygies--------------*/
1797//BOOLEAN   pVectorHasUnitM(poly p, int * k)
1798//{
1799//  while (p!=NULL)
1800//  {
1801//    if (pLmIsConstantComp(p))
1802//    {
1803//      *k = pGetComp(p);
1804//      return TRUE;
1805//    }
1806//    else pIter(p);
1807//  }
1808//  return FALSE;
1809//}
1810
1811BOOLEAN   pVectorHasUnitB(poly p, int * k)
1812{
1813  poly q=p,qq;
1814  int i;
1815
1816  while (q!=NULL)
1817  {
1818    if (pLmIsConstantComp(q))
1819    {
1820      i = pGetComp(q);
1821      qq = p;
1822      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1823      if (qq == q)
1824      {
1825        *k = i;
1826        return TRUE;
1827      }
1828      else
1829        pIter(q);
1830    }
1831    else pIter(q);
1832  }
1833  return FALSE;
1834}
1835
1836void   pVectorHasUnit(poly p, int * k, int * len)
1837{
1838  poly q=p,qq;
1839  int i,j=0;
1840
1841  *len = 0;
1842  while (q!=NULL)
1843  {
1844    if (pLmIsConstantComp(q))
1845    {
1846      i = pGetComp(q);
1847      qq = p;
1848      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1849      if (qq == q)
1850      {
1851       j = 0;
1852       while (qq!=NULL)
1853       {
1854         if (pGetComp(qq)==i) j++;
1855        pIter(qq);
1856       }
1857       if ((*len == 0) || (j<*len))
1858       {
1859         *len = j;
1860         *k = i;
1861       }
1862      }
1863    }
1864    pIter(q);
1865  }
1866}
1867
1868/*2
1869* returns TRUE if p1 = p2
1870*/
1871BOOLEAN p_EqualPolys(poly p1,poly p2, const ring r)
1872{
1873  while ((p1 != NULL) && (p2 != NULL))
1874  {
1875    if (! p_LmEqual(p1, p2,r))
1876      return FALSE;
1877    if (! n_Equal(p_GetCoeff(p1,r), p_GetCoeff(p2,r),r ))
1878      return FALSE;
1879    pIter(p1);
1880    pIter(p2);
1881  }
1882  return (p1==p2);
1883}
1884
1885/*2
1886*returns TRUE if p1 is a skalar multiple of p2
1887*assume p1 != NULL and p2 != NULL
1888*/
1889BOOLEAN pComparePolys(poly p1,poly p2)
1890{
1891  number n,nn;
1892  int i;
1893  pAssume(p1 != NULL && p2 != NULL);
1894
1895  if (!pLmEqual(p1,p2)) //compare leading mons
1896      return FALSE;
1897  if ((pNext(p1)==NULL) && (pNext(p2)!=NULL))
1898     return FALSE;
1899  if ((pNext(p2)==NULL) && (pNext(p1)!=NULL))
1900     return FALSE;
1901  if (pLength(p1) != pLength(p2))
1902    return FALSE;
1903  n=nDiv(pGetCoeff(p1),pGetCoeff(p2));
1904  while ((p1 != NULL) /*&& (p2 != NULL)*/)
1905  {
1906    if ( ! pLmEqual(p1, p2))
1907    {
1908        nDelete(&n);
1909        return FALSE;
1910    }
1911    if (!nEqual(pGetCoeff(p1),nn=nMult(pGetCoeff(p2),n)))
1912    {
1913      nDelete(&n);
1914      nDelete(&nn);
1915      return FALSE;
1916    }
1917    nDelete(&nn);
1918    pIter(p1);
1919    pIter(p2);
1920  }
1921  nDelete(&n);
1922  return TRUE;
1923}
Note: See TracBrowser for help on using the repository browser.