source: git/kernel/polys1.cc @ a41f3aa

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