source: git/Singular/polys1.cc @ 2f436b

spielwiese
Last change on this file since 2f436b was 2f436b, checked in by Olaf Bachmann <obachman@…>, 23 years ago
* version 1-3-13: sparsemat improivements git-svn-id: file:///usr/local/Singular/svn/trunk@5003 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 22.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: polys1.cc,v 1.61 2000-12-31 15:14:42 obachman 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 "tok.h"
16#include "numbers.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
31/*-----------------------------------------------------------*/
32/*
33* the module weights for std
34*/
35static pFDegProc pOldFDeg;
36static pLDegProc pOldLDeg;
37static intvec * pModW;
38static BOOLEAN pOldLexOrder;
39
40static long pModDeg(poly p, ring r = currRing)
41{
42  return pOldFDeg(p, r)+(*pModW)[p_GetComp(p, r)-1];
43}
44
45void pSetModDeg(intvec *w)
46{
47  if (w!=NULL)
48  {
49    pModW = w;
50    pOldFDeg = pFDeg;
51    pOldLDeg = pLDeg;
52    pOldLexOrder = pLexOrder;
53    pSetDegProcs(pModDeg);
54    pLexOrder = TRUE;
55  }
56  else
57  {
58    pModW = NULL;
59    pRestoreDegProcs(pOldFDeg, pOldLDeg);
60    pLexOrder = pOldLexOrder;
61  }
62}
63
64
65/*2
66* subtract p2 from p1, p1 and p2 are destroyed
67* do not put attention on speed: the procedure is only used in the interpreter
68*/
69poly pSub(poly p1, poly p2)
70{
71  return pAdd(p1, pNeg(p2));
72}
73
74/*3
75*  create binomial coef.
76*/
77static number* pnBin(int exp)
78{
79  int e, i, h;
80  number x, y, *bin=NULL;
81
82  x = nInit(exp);
83  if (nIsZero(x))
84  {
85    nDelete(&x);
86    return bin;
87  }
88  h = (exp >> 1) + 1;
89  bin = (number *)omAlloc0(h*sizeof(number));
90  bin[1] = x;
91  if (exp < 4)
92    return bin;
93  i = exp - 1;
94  for (e=2; e<h; e++)
95  {
96//    if(!nIsZero(bin[e-1]))
97//    {
98      x = nInit(i);
99      i--;
100      y = nMult(x,bin[e-1]);
101      nDelete(&x);
102      x = nInit(e);
103      bin[e] = nIntDiv(y,x);
104      nDelete(&x);
105      nDelete(&y);
106//    }
107//    else
108//    {
109//      bin[e] = nInit(binom(exp,e));
110//    }
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 pMonPower(poly p, int exp)
133{
134  int i;
135
136  if(!nIsOne(pGetCoeff(p)))
137  {
138    number x, y;
139    y = pGetCoeff(p);
140    nPower(y,exp,&x);
141    nDelete(&y);
142    pSetCoeff0(p,x);
143  }
144  for (i=pVariables; i!=0; i--)
145  {
146    pMultExp(p,i, exp);
147  }
148  pSetm(p);
149  return p;
150}
151
152/*3
153* compute for monomials p*q
154* destroys p, keeps q
155*/
156static void pMonMult(poly p, poly q)
157{
158  number x, y;
159  int i;
160
161  y = pGetCoeff(p);
162  x = nMult(y,pGetCoeff(q));
163  nDelete(&y);
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  pExpVectorAdd(p,q);
171}
172
173/*3
174* compute for monomials p*q
175* keeps p, q
176*/
177static poly pMonMultC(poly p, poly q)
178{
179  number x;
180  int i;
181  poly r = pInit();
182
183  x = nMult(pGetCoeff(p),pGetCoeff(q));
184  pSetCoeff0(r,x);
185  pExpVectorSum(r,p, q);
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 pTwoMonPower(poly p, int exp)
195{
196  int eh, e, al;
197  poly *a;
198  poly tail, b, res, h;
199  number x;
200  number *bin = pnBin(exp);
201
202  tail = pNext(p);
203  if (bin == NULL)
204  {
205    pMonPower(p,exp);
206    pMonPower(tail,exp);
207#ifdef PDEBUG
208    pTest(p);
209#endif
210    return p;
211  }
212  eh = exp >> 1;
213  al = (exp + 1) * sizeof(poly);
214  a = (poly *)omAlloc(al);
215  a[1] = p;
216  for (e=1; e<exp; e++)
217  {
218    a[e+1] = pMonMultC(a[e],p);
219  }
220  res = a[exp];
221  b = pHead(tail);
222  for (e=exp-1; e>eh; e--)
223  {
224    h = a[e];
225    x = nMult(bin[exp-e],pGetCoeff(h));
226    pSetCoeff(h,x);
227    pMonMult(h,b);
228    res = pNext(res) = h;
229    pMonMult(b,tail);
230  }
231  for (e=eh; e!=0; e--)
232  {
233    h = a[e];
234    x = nMult(bin[e],pGetCoeff(h));
235    pSetCoeff(h,x);
236    pMonMult(h,b);
237    res = pNext(res) = h;
238    pMonMult(b,tail);
239  }
240  pDeleteLm(&tail);
241  pNext(res) = b;
242  pNext(b) = NULL;
243  res = a[exp];
244  omFreeSize((ADDRESS)a, al);
245  pnFreeBin(bin, exp);
246//  tail=res;
247// while((tail!=NULL)&&(pNext(tail)!=NULL))
248// {
249//   if(nIsZero(pGetCoeff(pNext(tail))))
250//   {
251//     pDeleteLm(&pNext(tail));
252//   }
253//   else
254//     pIter(tail);
255// }
256#ifdef PDEBUG
257  pTest(res);
258#endif
259  return res;
260}
261
262static poly pPow(poly p, int i)
263{
264  poly rc = pCopy(p);
265  i -= 2;
266  do
267  {
268    rc = pMult(rc,pCopy(p));
269    pNormalize(rc);
270    i--;
271  }
272  while (i != 0);
273  return pMult(rc,p);
274}
275
276/*2
277* returns the i-th power of p
278* p will be destroyed
279*/
280poly pPower(poly p, int i)
281{
282  poly rc=NULL;
283
284  if (i==0)
285  {
286    pDelete(&p);
287    return pOne();
288  }
289
290  if(p!=NULL)
291  {
292    if ( (i > 0) && ((unsigned long ) i > (currRing->bitmask)))
293    {
294      Werror("exponent %d is too large, max. is %d",i,currRing->bitmask);
295      return NULL;
296    }
297    switch (i)
298    {
299// cannot happen, see above
300//      case 0:
301//      {
302//        rc=pOne();
303//#ifdef DRING
304//        if ((pDRING) && (pdDFlag(p)==1))
305//        {
306//          pdSetDFlag(rc,1);
307//        }
308//#endif
309//        pDelete(&p);
310//        break;
311//      }
312      case 1:
313        rc=p;
314        break;
315      case 2:
316        rc=pMult(pCopy(p),p);
317        break;
318      default:
319        if (i < 0)
320        {
321          pDelete(&p);
322          return NULL;
323        }
324        else
325        {
326          rc = pNext(p);
327          if (rc == NULL)
328            return pMonPower(p,i);
329          /* else: binom ?*/
330          int char_p=rChar(currRing);
331          if (pNext(rc) != NULL)
332            return pPow(p,i);
333          if ((char_p==0) || (i<=char_p))
334            return pTwoMonPower(p,i);
335          poly p_p=pTwoMonPower(pCopy(p),char_p);
336          return pMult(pPower(p,i-char_p),p_p);
337        }
338      /*end default:*/
339    }
340  }
341  return rc;
342}
343
344/*2
345* returns the partial differentiate of a by the k-th variable
346* does not destroy the input
347*/
348poly pDiff(poly a, int k)
349{
350  poly res, f, last;
351  number t;
352
353  last = res = NULL;
354  while (a!=NULL)
355  {
356    if (pGetExp(a,k)!=0)
357    {
358      f = pLmInit(a);
359      t = nInit(pGetExp(a,k));
360      pSetCoeff0(f,nMult(t,pGetCoeff(a)));
361      nDelete(&t);
362      if (nIsZero(pGetCoeff(f)))
363        pDeleteLm(&f);
364      else
365      {
366        pDecrExp(f,k);
367        pSetm(f);
368        if (res==NULL)
369        {
370          res=last=f;
371        }
372        else
373        {
374          pNext(last)=f;
375          last=f;
376        }
377      }
378    }
379    pIter(a);
380  }
381  return res;
382}
383
384static poly pDiffOpM(poly a, poly b,BOOLEAN multiply)
385{
386  int i,j,s;
387  number n,h,hh;
388  poly p=pOne();
389  n=nMult(pGetCoeff(a),pGetCoeff(b));
390  for(i=pVariables;i>0;i--)
391  {
392    s=pGetExp(b,i);
393    if (s<pGetExp(a,i))
394    {
395      nDelete(&n);
396      pDeleteLm(&p);
397      return NULL;
398    }
399    if (multiply)
400    {
401      for(j=pGetExp(a,i); j>0;j--)
402      {
403        h = nInit(s);
404        hh=nMult(n,h);
405        nDelete(&h);
406        nDelete(&n);
407        n=hh;
408        s--;
409      }
410      pSetExp(p,i,s);
411    }
412    else
413    {
414      pSetExp(p,i,s-pGetExp(a,i));
415    }
416  }
417  pSetm(p);
418  /*if (multiply)*/ pSetCoeff(p,n);
419  return p;
420}
421
422poly pDiffOp(poly a, poly b,BOOLEAN multiply)
423{
424  poly result=NULL;
425  poly h;
426  for(;a!=NULL;pIter(a))
427  {
428    for(h=b;h!=NULL;pIter(h))
429    {
430      result=pAdd(result,pDiffOpM(a,h,multiply));
431    }
432  }
433  return result;
434}
435
436
437void pSplit(poly p, poly *h)
438{
439  *h=pNext(p);
440  pNext(p)=NULL;
441}
442
443
444
445int pMaxCompProc(poly p)
446{
447  return pMaxComp(p);
448}
449
450/*2
451* handle memory request for sets of polynomials (ideals)
452* l is the length of *p, increment is the difference (may be negative)
453*/
454void pEnlargeSet(polyset *p, int l, int increment)
455{
456  int i;
457  polyset h;
458
459  h=(polyset)omReallocSize((poly*)*p,l*sizeof(poly),(l+increment)*sizeof(poly));
460  if (increment>0)
461  {
462    //for (i=l; i<l+increment; i++)
463    //  h[i]=NULL;
464    memset(&(h[l]),0,increment*sizeof(poly));
465  }
466  *p=h;
467}
468
469void pContent(poly ph)
470{
471  number h,d;
472  poly p;
473
474  if(pNext(ph)==NULL)
475  {
476    pSetCoeff(ph,nInit(1));
477  }
478  else
479  {
480    nNormalize(pGetCoeff(ph));
481    if(!nGreaterZero(pGetCoeff(ph))) ph = pNeg(ph);
482    h=nCopy(pGetCoeff(ph));
483    p = pNext(ph);
484    while (p!=NULL)
485    {
486      nNormalize(pGetCoeff(p));
487      d=nGcd(h,pGetCoeff(p));
488      nDelete(&h);
489      h = d;
490      if(nIsOne(h))
491      {
492        break;
493      }
494      pIter(p);
495    }
496    p = ph;
497    //number tmp;
498    if(!nIsOne(h))
499    {
500      while (p!=NULL)
501      {
502        //d = nDiv(pGetCoeff(p),h);
503        //tmp = nIntDiv(pGetCoeff(p),h);
504        //if (!nEqual(d,tmp))
505        //{
506        //  StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
507        //  nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
508        //  nWrite(tmp);Print(StringAppendS("\n"));
509        //}
510        //nDelete(&tmp);
511        d = nIntDiv(pGetCoeff(p),h);
512        pSetCoeff(p,d);
513        pIter(p);
514      }
515    }
516    nDelete(&h);
517#ifdef HAVE_FACTORY
518    if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
519    {
520      singclap_divide_content(ph);
521      if(!nGreaterZero(pGetCoeff(ph))) ph = pNeg(ph);
522    }
523#endif
524  }
525}
526
527//void pContent(poly ph)
528//{
529//  number h,d;
530//  poly p;
531//
532//  p = ph;
533//  if(pNext(p)==NULL)
534//  {
535//    pSetCoeff(p,nInit(1));
536//  }
537//  else
538//  {
539//#ifdef PDEBUG
540//    if (!pTest(p)) return;
541//#endif
542//    nNormalize(pGetCoeff(p));
543//    if(!nGreaterZero(pGetCoeff(ph)))
544//    {
545//      ph = pNeg(ph);
546//      nNormalize(pGetCoeff(p));
547//    }
548//    h=pGetCoeff(p);
549//    pIter(p);
550//    while (p!=NULL)
551//    {
552//      nNormalize(pGetCoeff(p));
553//      if (nGreater(h,pGetCoeff(p))) h=pGetCoeff(p);
554//      pIter(p);
555//    }
556//    h=nCopy(h);
557//    p=ph;
558//    while (p!=NULL)
559//    {
560//      d=nGcd(h,pGetCoeff(p));
561//      nDelete(&h);
562//      h = d;
563//      if(nIsOne(h))
564//      {
565//        break;
566//      }
567//      pIter(p);
568//    }
569//    p = ph;
570//    //number tmp;
571//    if(!nIsOne(h))
572//    {
573//      while (p!=NULL)
574//      {
575//        d = nIntDiv(pGetCoeff(p),h);
576//        pSetCoeff(p,d);
577//        pIter(p);
578//      }
579//    }
580//    nDelete(&h);
581//#ifdef HAVE_FACTORY
582//    if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
583//    {
584//      pTest(ph);
585//      singclap_divide_content(ph);
586//      pTest(ph);
587//    }
588//#endif
589//  }
590//}
591
592void pCleardenom(poly ph)
593{
594  number d, h;
595  poly p;
596
597  p = ph;
598  if(pNext(p)==NULL)
599  {
600    pSetCoeff(p,nInit(1));
601  }
602  else
603  {
604    h = nInit(1);
605    while (p!=NULL)
606    {
607      nNormalize(pGetCoeff(p));
608      d=nLcm(h,pGetCoeff(p));
609      nDelete(&h);
610      h=d;
611      pIter(p);
612    }
613    /* contains the 1/lcm of all denominators */
614    if(!nIsOne(h))
615    {
616      p = ph;
617      while (p!=NULL)
618      {
619        /* should be:
620        * number hh;
621        * nGetDenom(p->coef,&hh);
622        * nMult(&h,&hh,&d);
623        * nNormalize(d);
624        * nDelete(&hh);
625        * nMult(d,p->coef,&hh);
626        * nDelete(&d);
627        * nDelete(&(p->coef));
628        * p->coef =hh;
629        */
630        d=nMult(h,pGetCoeff(p));
631        nNormalize(d);
632        pSetCoeff(p,d);
633        pIter(p);
634      }
635      nDelete(&h);
636      if (nGetChar()==1)
637      {
638        loop
639        {
640          h = nInit(1);
641          p=ph;
642          while (p!=NULL)
643          {
644            d=nLcm(h,pGetCoeff(p));
645            nDelete(&h);
646            h=d;
647            pIter(p);
648          }
649          /* contains the 1/lcm of all denominators */
650          if(!nIsOne(h))
651          {
652            p = ph;
653            while (p!=NULL)
654            {
655              /* should be:
656              * number hh;
657              * nGetDenom(p->coef,&hh);
658              * nMult(&h,&hh,&d);
659              * nNormalize(d);
660              * nDelete(&hh);
661              * nMult(d,p->coef,&hh);
662              * nDelete(&d);
663              * nDelete(&(p->coef));
664              * p->coef =hh;
665              */
666              d=nMult(h,pGetCoeff(p));
667              nNormalize(d);
668              pSetCoeff(p,d);
669              pIter(p);
670            }
671            nDelete(&h);
672          }
673          else
674          {
675            nDelete(&h);
676            break;
677          }
678        }
679      }
680    }
681    if (h!=NULL) nDelete(&h);
682    pContent(ph);
683  }
684}
685
686/*2
687*tests if p is homogeneous with respect to the actual weigths
688*/
689BOOLEAN pIsHomogeneous (poly p)
690{
691  poly qp=p;
692  int o;
693
694  if ((p == NULL) || (pNext(p) == NULL)) return TRUE;
695  pFDegProc d=(pLexOrder ? pTotaldegree : pFDeg );
696  o = d(p);
697  do
698  {
699    if (d(qp) != o) return FALSE;
700    pIter(qp);
701  }
702  while (qp != NULL);
703  return TRUE;
704}
705
706
707/*2
708*returns a re-ordered copy of a polynomial, with permutation of the variables
709*/
710poly pPermPoly (poly p, int * perm, ring oldRing, nMapFunc nMap,
711   int *par_perm, int OldPar)
712{
713  int OldpVariables = oldRing->N;
714  poly result = NULL;
715  poly result_last = NULL;
716  poly aq=NULL; /* the map coefficient */
717  poly qq; /* the mapped monomial */
718
719  while (p != NULL)
720  {
721    if (OldPar==0)
722    {
723      qq = pInit();
724      number n=nMap(pGetCoeff(p));
725      if ((currRing->minpoly!=NULL)
726      && ((rField_is_Zp_a()) || (rField_is_Q_a())))
727      {
728        nNormalize(n);
729      }
730      pGetCoeff(qq)=n;
731    // coef may be zero:  pTest(qq);
732    }
733    else
734    {
735      qq=pOne();
736      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar, oldRing);
737      if ((currRing->minpoly!=NULL)
738      && ((rField_is_Zp_a()) || (rField_is_Q_a())))
739      {
740        poly tmp=aq;
741        while (tmp!=NULL)
742        {
743          number n=pGetCoeff(tmp);
744          nNormalize(n);
745          pGetCoeff(tmp)=n;
746          pIter(tmp);
747        }
748      }
749      pTest(aq);
750    }
751    pSetComp(qq, p_GetComp(p,oldRing));
752    if (nIsZero(pGetCoeff(qq)))
753    {
754      pDeleteLm(&qq);
755    }
756    else
757    {
758      int i;
759      int mapped_to_par=0;
760      for(i=1; i<=OldpVariables; i++)
761      {
762        int e=p_GetExp(p,i,oldRing);
763        if (e!=0)
764        {
765          if (perm==NULL)
766          {
767            pSetExp(qq,i, e);
768          }
769          else if (perm[i]>0)
770            pAddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/);
771          else if (perm[i]<0)
772          {
773            lnumber c=(lnumber)pGetCoeff(qq);
774            napAddExp(c->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
775            mapped_to_par=1;
776          }
777          else
778          {
779            /* this variable maps to 0 !*/
780            pDeleteLm(&qq);
781            break;
782          }
783        }
784      }
785      if (mapped_to_par
786      && (currRing->minpoly!=NULL))
787      {
788        number n=pGetCoeff(qq);
789        nNormalize(n);
790        pGetCoeff(qq)=n;
791      }
792    }
793    pIter(p);
794#if 1
795    if (qq!=NULL)
796    {
797      pSetm(qq);
798      pTest(aq);
799      pTest(qq);
800      if (aq!=NULL) qq=pMult(aq,qq);
801      aq = qq;
802      while (pNext(aq) != NULL) pIter(aq);
803      if (result_last==NULL)
804      {
805        result=qq;
806      }
807      else
808      {
809        pNext(result_last)=qq;
810      }
811      result_last=aq;
812      aq = NULL;
813    }
814    else if (aq!=NULL)
815    {
816      pDelete(&aq);
817    }
818  }
819  result=pSortAdd(result);
820#else
821  //  if (qq!=NULL)
822  //  {
823  //    pSetm(qq);
824  //    pTest(qq);
825  //    pTest(aq);
826  //    if (aq!=NULL) qq=pMult(aq,qq);
827  //    aq = qq;
828  //    while (pNext(aq) != NULL) pIter(aq);
829  //    pNext(aq) = result;
830  //    aq = NULL;
831  //    result = qq;
832  //  }
833  //  else if (aq!=NULL)
834  //  {
835  //    pDelete(&aq);
836  //  }
837  //}
838  //p = result;
839  //result = NULL;
840  //while (p != NULL)
841  //{
842  //  qq = p;
843  //  pIter(p);
844  //  qq->next = NULL;
845  //  result = pAdd(result, qq);
846  //}
847#endif
848  pTest(result);
849  return result;
850}
851
852#if 0
853/*2
854*returns a re-ordered copy of a polynomial, with permutation of the variables
855*/
856poly p_PermPoly (poly p, int * perm, ring oldRing,
857   int *par_perm, int OldPar, ring newRing)
858{
859  int OldpVariables = oldRing->N;
860  poly result = NULL;
861  poly result_last = NULL;
862  poly aq=NULL; /* the map coefficient */
863  poly qq; /* the mapped monomial */
864
865  while (p != NULL)
866  {
867    if (OldPar==0)
868    {
869      qq = pInit();
870      number n=newRing->cf->nMap(pGetCoeff(p));
871      if ((newRing->minpoly!=NULL)
872      && ((rField_is_Zp_a(newRing)) || (rField_is_Q_a(newRing))))
873      {
874        newRing->cf->nNormalize(n);
875      }
876      pGetCoeff(qq)=n;
877    // coef may be zero:  pTest(qq);
878    }
879    else
880    {
881      qq=p_ISet(1, newRing);
882      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar, oldRing);
883      if ((newRing->minpoly!=NULL)
884      && ((rField_is_Zp_a(newRing)) || (rField_is_Q_a(newRing))))
885      {
886        poly tmp=aq;
887        while (tmp!=NULL)
888        {
889          number n=pGetCoeff(tmp);
890          newRing->cf->nNormalize(n);
891          pGetCoeff(tmp)=n;
892          pIter(tmp);
893        }
894      }
895      //pTest(aq);
896    }
897    p_SetComp(qq, p_GetComp(p,oldRing), newRing);
898    if (newRing->cf->nIsZero(pGetCoeff(qq)))
899    {
900      p_DeleteLm(&qq, newRing);
901    }
902    else
903    {
904      int i;
905      int mapped_to_par=0;
906      for(i=1; i<=OldpVariables; i++)
907      {
908        int e=p_GetExp(p,i,oldRing);
909        if (e!=0)
910        {
911          if (perm==NULL)
912          {
913            p_SetExp(qq,i, e, newRing);
914          }
915          else if (perm[i]>0)
916            p_AddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/, newRing);
917          else if (perm[i]<0)
918          {
919            lnumber c=(lnumber)pGetCoeff(qq);
920            napAddExp(c->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
921            mapped_to_par=1;
922          }
923          else
924          {
925            /* this variable maps to 0 !*/
926            p_DeleteLm(&qq, newRing);
927            break;
928          }
929        }
930      }
931      if (mapped_to_par
932      && (newRing->minpoly!=NULL))
933      {
934        number n=pGetCoeff(qq);
935        newRing->cf->nNormalize(n);
936        pGetCoeff(qq)=n;
937      }
938    }
939    pIter(p);
940    if (qq!=NULL)
941    {
942      p_Setm(qq, newRing);
943      //pTest(aq);
944      //pTest(qq);
945      if (aq!=NULL) qq=pMult(aq,qq);
946      aq = qq;
947      while (pNext(aq) != NULL) pIter(aq);
948      if (result_last==NULL)
949      {
950        result=qq;
951      }
952      else
953      {
954        pNext(result_last)=qq;
955      }
956      result_last=aq;
957      aq = NULL;
958    }
959    else if (aq!=NULL)
960    {
961      p_Delete(&aq, newRing);
962    }
963  }
964  result=pOrdPolyMerge(result);
965  //pTest(result);
966  return result;
967}
968#endif
969
970poly pJet(poly p, int m)
971{
972  poly r=NULL;
973  poly t=NULL;
974
975  while (p!=NULL)
976  {
977    if (pTotaldegree(p)<=m)
978    {
979      if (r==NULL)
980        r=pHead(p);
981      else
982      if (t==NULL)
983      {
984        pNext(r)=pHead(p);
985        t=pNext(r);
986      }
987      else
988      {
989        pNext(t)=pHead(p);
990        pIter(t);
991      }
992    }
993    pIter(p);
994  }
995  return r;
996}
997
998poly pJetW(poly p, int m, short *w)
999{
1000  poly r=NULL;
1001  poly t=NULL;
1002  short *wsave=ecartWeights;
1003
1004  ecartWeights=w;
1005
1006  while (p!=NULL)
1007  {
1008    if (totaldegreeWecart(p)<=m)
1009    {
1010      if (r==NULL)
1011        r=pHead(p);
1012      else
1013      if (t==NULL)
1014      {
1015        pNext(r)=pHead(p);
1016        t=pNext(r);
1017      }
1018      else
1019      {
1020        pNext(t)=pHead(p);
1021        pIter(t);
1022      }
1023    }
1024    pIter(p);
1025  }
1026  ecartWeights=wsave;
1027  return r;
1028}
1029
1030int pDegW(poly p, short *w)
1031{
1032  int r=0;
1033  short *wsave=ecartWeights;
1034
1035  ecartWeights=w;
1036
1037  while (p!=NULL)
1038  {
1039    r=max(r, totaldegreeWecart(p));
1040    pIter(p);
1041  }
1042  ecartWeights=wsave;
1043  return r;
1044}
1045
1046/*-----------type conversions ----------------------------*/
1047/*2
1048* input: a set of polys (len elements: p[0]..p[len-1])
1049* output: a vector
1050* p will not be changed
1051*/
1052poly  pPolys2Vec(polyset p, int len)
1053{
1054  poly v=NULL;
1055  poly h;
1056  int i;
1057
1058  for (i=len-1; i>=0; i--)
1059  {
1060    if (p[i])
1061    {
1062      h=pCopy(p[i]);
1063      pSetCompP(h,i+1);
1064      v=pAdd(v,h);
1065    }
1066  }
1067 return v;
1068}
1069
1070/*2
1071* convert a vector to a set of polys,
1072* allocates the polyset, (entries 0..(*len)-1)
1073* the vector will not be changed
1074*/
1075void  pVec2Polys(poly v, polyset *p, int *len)
1076{
1077  poly h;
1078  int k;
1079
1080  *len=pMaxComp(v);
1081  if (*len==0) *len=1;
1082  *p=(polyset)omAlloc0((*len)*sizeof(poly));
1083  while (v!=NULL)
1084  {
1085    h=pHead(v);
1086    k=pGetComp(h);
1087    pSetComp(h,0);
1088    (*p)[k-1]=pAdd((*p)[k-1],h);
1089    pIter(v);
1090  }
1091}
1092
1093int pVar(poly m)
1094{
1095  if (m==NULL) return 0;
1096  if (pNext(m)!=NULL) return 0;
1097  int i,e=0;
1098  for (i=pVariables; i>0; i--)
1099  {
1100    if (pGetExp(m,i)==1)
1101    {
1102      if (e==0) e=i;
1103      else return 0;
1104    }
1105  }
1106  return e;
1107}
1108
1109/*----------utilities for syzygies--------------*/
1110//BOOLEAN   pVectorHasUnitM(poly p, int * k)
1111//{
1112//  while (p!=NULL)
1113//  {
1114//    if (pLmIsConstantComp(p))
1115//    {
1116//      *k = pGetComp(p);
1117//      return TRUE;
1118//    }
1119//    else pIter(p);
1120//  }
1121//  return FALSE;
1122//}
1123
1124BOOLEAN   pVectorHasUnitB(poly p, int * k)
1125{
1126  poly q=p,qq;
1127  int i;
1128
1129  while (q!=NULL)
1130  {
1131    if (pLmIsConstantComp(q))
1132    {
1133      i = pGetComp(q);
1134      qq = p;
1135      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1136      if (qq == q)
1137      {
1138        *k = i;
1139        return TRUE;
1140      }
1141      else
1142        pIter(q);
1143    }
1144    else pIter(q);
1145  }
1146  return FALSE;
1147}
1148
1149void   pVectorHasUnit(poly p, int * k, int * len)
1150{
1151  poly q=p,qq;
1152  int i,j=0;
1153
1154  *len = 0;
1155  while (q!=NULL)
1156  {
1157    if (pLmIsConstantComp(q))
1158    {
1159      i = pGetComp(q);
1160      qq = p;
1161      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1162      if (qq == q)
1163      {
1164       j = 0;
1165       while (qq!=NULL)
1166       {
1167         if (pGetComp(qq)==i) j++;
1168        pIter(qq);
1169       }
1170       if ((*len == 0) || (j<*len))
1171       {
1172         *len = j;
1173         *k = i;
1174       }
1175      }
1176    }
1177    pIter(q);
1178  }
1179}
1180
1181/*2
1182* returns TRUE if p1 = p2
1183*/
1184BOOLEAN pEqualPolys(poly p1,poly p2)
1185{
1186  while ((p1 != NULL) && (p2 != NULL))
1187  {
1188    if (! pLmEqual(p1, p2))
1189      return FALSE;
1190    if (! nEqual(pGetCoeff(p1), pGetCoeff(p2)))
1191      return FALSE;
1192    pIter(p1);
1193    pIter(p2);
1194  }
1195  return (p1==p2);
1196}
1197
1198/*2
1199*returns TRUE if p1 is a skalar multiple of p2
1200*assume p1 != NULL and p2 != NULL
1201*/
1202BOOLEAN pComparePolys(poly p1,poly p2)
1203{
1204  number n,nn;
1205  int i;
1206  pAssume(p1 != NULL && p2 != NULL);
1207
1208  if (!pLmEqual(p1,p2)) //compare leading mons
1209      return FALSE;
1210  if ((pNext(p1)==NULL) && (pNext(p2)!=NULL))
1211     return FALSE;
1212  if ((pNext(p2)==NULL) && (pNext(p1)!=NULL))
1213     return FALSE;
1214  if (pLength(p1) != pLength(p2))
1215    return FALSE;
1216  n=nDiv(pGetCoeff(p1),pGetCoeff(p2));
1217  while ((p1 != NULL) /*&& (p2 != NULL)*/)
1218  {
1219    if ( ! pLmEqual(p1, p2))
1220    {
1221        nDelete(&n);
1222        return FALSE;
1223    }
1224    if (!nEqual(pGetCoeff(p1),nn=nMult(pGetCoeff(p2),n)))
1225    {
1226      nDelete(&n);
1227      nDelete(&nn);
1228      return FALSE;
1229    }
1230    nDelete(&nn);
1231    pIter(p1);
1232    pIter(p2);
1233  }
1234  nDelete(&n);
1235  return TRUE;
1236}
Note: See TracBrowser for help on using the repository browser.