source: git/Singular/polys1.cc @ ce9c2d

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