source: git/Singular/polys1.cc @ 8a150b

spielwiese
Last change on this file since 8a150b was 8a150b, checked in by Hans Schönemann <hannes@…>, 25 years ago
* hannes: added long reals git-svn-id: file:///usr/local/Singular/svn/trunk@3012 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.21 1999-04-29 11:38: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 "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  = pInit();
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            return pMonPower(p,i);
495          /* else: binom */
496          int char_p=rChar(currRing);
497          if (pNext(rc) != NULL)
498            return pPow(p,i);
499          if ((char_p==0) || (i<=char_p))
500            return pTwoMonPower(p,i);
501          poly p_p=pTwoMonPower(pCopy(p),char_p);
502          return pMult(pPower(p,i-char_p),p_p);
503        }
504      /*end default:*/
505    }
506  }
507  return rc;
508}
509
510/*2
511* returns the partial differentiate of a by the k-th variable
512* does not destroy the input
513*/
514poly pDiff(poly a, int k)
515{
516  poly res, f, last;
517  number t;
518
519  last = res = NULL;
520  while (a!=NULL)
521  {
522    if (pGetExp(a,k)!=0)
523    {
524      f = pNew();
525      pCopy2(f,a);
526      pNext(f)=NULL;
527      t = nInit(pGetExp(a,k));
528      pSetCoeff0(f,nMult(t,pGetCoeff(a)));
529      nDelete(&t);
530      if (nIsZero(pGetCoeff(f)))
531        pDelete1(&f);
532      else
533      {
534        pDecrExp(f,k);
535        pSetm(f);
536        if (res==NULL)
537        {
538          res=last=f;
539        }
540        else
541        {
542          pNext(last)=f;
543          last=f;
544        }
545        //res = pAdd(res, f);
546      }
547    }
548    pIter(a);
549  }
550  return res;
551}
552
553static poly pDiffOpM(poly a, poly b,BOOLEAN multiply)
554{
555  int i,j,s;
556  number n,h,hh;
557  poly p=pOne();
558  n=nMult(pGetCoeff(a),pGetCoeff(b));
559  for(i=pVariables;i>0;i--)
560  {
561    s=pGetExp(b,i);
562    if (s<pGetExp(a,i))
563    {
564      nDelete(&n);
565      pDelete1(&p);
566      return NULL;
567    }
568    if (multiply)
569    {
570      for(j=pGetExp(a,i); j>0;j--)
571      {
572        h = nInit(s);
573        hh=nMult(n,h);
574        nDelete(&h);
575        nDelete(&n);
576        n=hh;
577        s--;
578      }
579      pSetExp(p,i,s);
580    }
581    else
582    {
583      pSetExp(p,i,s-pGetExp(a,i));
584    }
585  }
586  pSetm(p);
587  /*if (multiply)*/ pSetCoeff(p,n);
588  return p;
589}
590
591poly pDiffOp(poly a, poly b,BOOLEAN multiply)
592{
593  poly result=NULL;
594  poly h;
595  for(;a!=NULL;pIter(a))
596  {
597    for(h=b;h!=NULL;pIter(h))
598    {
599      result=pAdd(result,pDiffOpM(a,h,multiply));
600    }
601  }
602  return result;
603}
604
605/*2
606* returns the length of a (numbers of monomials)
607*/
608int pLength(poly a)
609{
610  int l = 0;
611
612  while (a!=NULL)
613  {
614    pIter(a);
615    l++;
616  }
617  return l;
618}
619
620
621void pSplit(poly p, poly *h)
622{
623  *h=pNext(p);
624  pNext(p)=NULL;
625}
626
627/*2
628* returns maximal column number in the modul element a (or 0)
629*/
630int pMaxComp(poly p)
631{
632  int result,i;
633
634  if(p==NULL) return 0;
635  result = pGetComp(p);
636  while (pNext(p)!=NULL)
637  {
638    pIter(p);
639    i = pGetComp(p);
640    if (i>result) result = i;
641  }
642  return result;
643}
644
645/*2
646* returns minimal column number in the modul element a (or 0)
647*/
648int pMinComp(poly p)
649{
650  int result,i;
651
652  if(p==NULL) return 0;
653  result = pGetComp(p);
654  while (pNext(p)!=NULL)
655  {
656    pIter(p);
657    i = pGetComp(p);
658    if (i<result) result = i;
659  }
660  return result;
661}
662
663/*2
664* returns TRUE, if all monoms have the same component
665*/
666BOOLEAN pOneComp(poly p)
667{
668  if(p!=NULL)
669  {
670    int i = pGetComp(p);
671    while (pNext(p)!=NULL)
672    {
673      pIter(p);
674      if(i != pGetComp(p)) return FALSE;
675    }
676  }
677  return TRUE;
678}
679
680/*2
681* multiplies the polynomial a by the column generator with number i
682*/
683void pSetCompP(poly p, int i)
684{
685
686  if ((p!=NULL) && (pGetComp(p)==0))
687  {
688    do
689    {
690      pSetComp(p, (Exponent_t)i);
691      pIter(p);
692    } while (p!=NULL);
693  }
694}
695
696/*2
697* handle memory request for sets of polynomials (ideals)
698* l is the length of *p, increment is the difference (may be negative)
699*/
700void pEnlargeSet(polyset *p, int l, int increment)
701{
702  int i;
703  polyset h;
704
705  h=(polyset)ReAlloc((poly*)*p,l*sizeof(poly),(l+increment)*sizeof(poly));
706  if (increment>0)
707  {
708    //for (i=l; i<l+increment; i++)
709    //  h[i]=NULL;
710    memset(&(h[l]),0,increment*sizeof(poly));
711  }
712  *p=h;
713}
714
715/*2
716* returns a polynomial representing the integer i
717*/
718poly pISet(int i)
719{
720  poly rc = NULL;
721  if (i!=0)
722  {
723    rc = pInit();
724    pSetCoeff0(rc,nInit(i));
725    if (nIsZero(pGetCoeff(rc)))
726      pDelete1(&rc);
727#ifdef TEST_MAC_ORDER
728    else if (bNoAdd)
729      pSetm(rc);
730#endif
731  }
732  return rc;
733}
734
735void pContent(poly ph)
736{
737  number h,d;
738  poly p;
739
740  if(pNext(ph)==NULL)
741  {
742    pSetCoeff(ph,nInit(1));
743  }
744  else
745  {
746#ifdef PDEBUG
747    if (!pTest(ph)) return;
748#endif
749    nNormalize(pGetCoeff(ph));
750    if(!nGreaterZero(pGetCoeff(ph))) ph = pNeg(ph);
751    h=nCopy(pGetCoeff(ph));
752    p = pNext(ph);
753    while (p!=NULL)
754    {
755      nNormalize(pGetCoeff(p));
756      d=nGcd(h,pGetCoeff(p));
757      nDelete(&h);
758      h = d;
759      if(nIsOne(h))
760      {
761        break;
762      }
763      pIter(p);
764    }
765    p = ph;
766    //number tmp;
767    if(!nIsOne(h))
768    {
769      while (p!=NULL)
770      {
771        //d = nDiv(pGetCoeff(p),h);
772        //tmp = nIntDiv(pGetCoeff(p),h);
773        //if (!nEqual(d,tmp))
774        //{
775        //  StringSetS("** div0:");nWrite(pGetCoeff(p));StringAppendS("/");
776        //  nWrite(h);StringAppendS("=");nWrite(d);StringAppendS(" int:");
777        //  nWrite(tmp);Print(StringAppendS("\n"));
778        //}
779        //nDelete(&tmp);
780        d = nIntDiv(pGetCoeff(p),h);
781        pSetCoeff(p,d);
782        pIter(p);
783      }
784    }
785    nDelete(&h);
786#ifdef HAVE_FACTORY
787    if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
788    {
789      singclap_divide_content(ph);
790      if(!nGreaterZero(pGetCoeff(ph))) ph = pNeg(ph);
791    }
792#endif
793  }
794  pTest(ph);
795}
796
797//void pContent(poly ph)
798//{
799//  number h,d;
800//  poly p;
801//
802//  p = ph;
803//  if(pNext(p)==NULL)
804//  {
805//    pSetCoeff(p,nInit(1));
806//  }
807//  else
808//  {
809//#ifdef PDEBUG
810//    if (!pTest(p)) return;
811//#endif
812//    nNormalize(pGetCoeff(p));
813//    if(!nGreaterZero(pGetCoeff(ph)))
814//    {
815//      ph = pNeg(ph);
816//      nNormalize(pGetCoeff(p));
817//    }
818//    h=pGetCoeff(p);
819//    pIter(p);
820//    while (p!=NULL)
821//    {
822//      nNormalize(pGetCoeff(p));
823//      if (nGreater(h,pGetCoeff(p))) h=pGetCoeff(p);
824//      pIter(p);
825//    }
826//    h=nCopy(h);
827//    p=ph;
828//    while (p!=NULL)
829//    {
830//      d=nGcd(h,pGetCoeff(p));
831//      nDelete(&h);
832//      h = d;
833//      if(nIsOne(h))
834//      {
835//        break;
836//      }
837//      pIter(p);
838//    }
839//    p = ph;
840//    //number tmp;
841//    if(!nIsOne(h))
842//    {
843//      while (p!=NULL)
844//      {
845//        d = nIntDiv(pGetCoeff(p),h);
846//        pSetCoeff(p,d);
847//        pIter(p);
848//      }
849//    }
850//    nDelete(&h);
851//#ifdef HAVE_FACTORY
852//    if ( (nGetChar() == 1) || (nGetChar() < 0) ) /* Q[a],Q(a),Zp[a],Z/p(a) */
853//    {
854//      pTest(ph);
855//      singclap_divide_content(ph);
856//      pTest(ph);
857//    }
858//#endif
859//  }
860//}
861
862void pCleardenom(poly ph)
863{
864  number d, h;
865  poly p;
866
867  p = ph;
868  if(pNext(p)==NULL)
869  {
870    pSetCoeff(p,nInit(1));
871  }
872  else
873  {
874    h = nInit(1);
875    while (p!=NULL)
876    {
877      nNormalize(pGetCoeff(p));
878      d=nLcm(h,pGetCoeff(p));
879      nDelete(&h);
880      h=d;
881      pIter(p);
882    }
883    /* contains the 1/lcm of all denominators */
884    if(!nIsOne(h))
885    {
886      p = ph;
887      while (p!=NULL)
888      {
889        /* should be:
890        * number hh;
891        * nGetDenom(p->coef,&hh);
892        * nMult(&h,&hh,&d);
893        * nNormalize(d);
894        * nDelete(&hh);
895        * nMult(d,p->coef,&hh);
896        * nDelete(&d);
897        * nDelete(&(p->coef));
898        * p->coef =hh;
899        */
900        d=nMult(h,pGetCoeff(p));
901        nNormalize(d);
902        pSetCoeff(p,d);
903        pIter(p);
904      }
905      nDelete(&h);
906      if (nGetChar()==1)
907      {
908        h = nInit(1);
909        p=ph;
910        while (p!=NULL)
911        {
912          d=nLcm(h,pGetCoeff(p));
913          nDelete(&h);
914          h=d;
915          pIter(p);
916        }
917        /* contains the 1/lcm of all denominators */
918        if(!nIsOne(h))
919        {
920          p = ph;
921          while (p!=NULL)
922          {
923            /* should be:
924            * number hh;
925            * nGetDenom(p->coef,&hh);
926            * nMult(&h,&hh,&d);
927            * nNormalize(d);
928            * nDelete(&hh);
929            * nMult(d,p->coef,&hh);
930            * nDelete(&d);
931            * nDelete(&(p->coef));
932            * p->coef =hh;
933            */
934            d=nMult(h,pGetCoeff(p));
935            nNormalize(d);
936            pSetCoeff(p,d);
937            pIter(p);
938          }
939          nDelete(&h);
940        }
941      }
942    }
943    pContent(ph);
944  }
945}
946
947/*2
948*tests if p is homogeneous with respect to the actual weigths
949*/
950BOOLEAN pIsHomogeneous (poly p)
951{
952  poly qp=p;
953  int o;
954
955  if ((p == NULL) || (pNext(p) == NULL)) return TRUE;
956  pFDegProc d=(pLexOrder ? pTotaldegree : pFDeg );
957  o = d(p);
958  do
959  {
960    if (d(qp) != o) return FALSE;
961    pIter(qp);
962  }
963  while (qp != NULL);
964  return TRUE;
965}
966
967// orders monoms of poly using merge sort (ususally faster than
968// insertion sort). ASSUMES that pSetm was performed on monoms
969// (i.e. that Order field is set correctly)
970poly pOrdPolyMerge(poly p)
971{
972  poly qq,pp,result=NULL;
973
974  if (p == NULL) return NULL;
975
976  for (;;)
977  {
978    qq = p;
979    for (;;)
980    {
981      if (pNext(p) == NULL) return pAdd(result, qq);
982      if (pComp(p,pNext(p)) != 1)
983      {
984        pp = p;
985        pIter(p);
986        pNext(pp) = NULL;
987        result = pAdd(result, qq);
988        break;
989      }
990      pIter(p);
991    }
992  }
993}
994
995// orders monoms of poly using insertion sort, performs pSetm on each
996// monom (i.e. sets Order field)
997poly pOrdPolyInsertSetm(poly p)
998{
999  poly qq,result = NULL;
1000
1001  while (p != NULL)
1002  {
1003    qq = p;
1004    pIter(p);
1005    qq->next = NULL;
1006    pSetm(qq);
1007    result = pAdd(result,qq);
1008    pTest(result);
1009  }
1010  return result;
1011}
1012
1013/*2
1014*returns a re-ordered copy of a polynomial, with permutation of the variables
1015*/
1016poly pPermPoly (poly p, int * perm, ring oldRing,
1017   int *par_perm, int OldPar)
1018{
1019  int OldpVariables = oldRing->N;
1020  poly result = NULL;
1021  poly aq=NULL;
1022  poly qq;
1023  int i;
1024
1025  while (p != NULL)
1026  {
1027    if (OldPar==0)
1028    {
1029      qq = pInit();
1030      pGetCoeff(qq)=nMap(pGetCoeff(p));
1031    }
1032    else
1033    {
1034      qq=pOne();
1035      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar);
1036      pTest(aq);
1037    }
1038    pSetComp(qq, pRingGetComp(oldRing,p));
1039    if (nIsZero(pGetCoeff(qq)))
1040    {
1041      pDelete1(&qq);
1042    }
1043    else
1044    {
1045      for(i=1; i<=OldpVariables; i++)
1046      {
1047        if (pRingGetExp(oldRing,p,i)!=0)
1048        {
1049          if (perm==NULL)
1050          {
1051            pSetExp(qq,i, pRingGetExp(oldRing,p,i));
1052          }
1053          else if (perm[i]>0)
1054            pAddExp(qq,perm[i], pRingGetExp(oldRing, p,i));
1055          else if (perm[i]<0)
1056          {
1057            lnumber c=(lnumber)pGetCoeff(qq);
1058            c->z->e[-perm[i]-1]+=pRingGetExp(oldRing, p,i);
1059          }
1060          else
1061          {
1062            /* this variable maps to 0 !*/
1063            pDelete1(&qq);
1064            break;
1065          }
1066        }
1067      }
1068    }
1069    pIter(p);
1070    if (qq!=NULL)
1071    {
1072      pSetm(qq);
1073      pTest(qq);
1074      pTest(aq);
1075      if (aq!=NULL)
1076      {
1077        qq=pMult(aq,qq);
1078        aq=NULL;
1079      }
1080      result = pAdd(result,qq);
1081      pTest(result);
1082    }
1083    else if (aq!=NULL)
1084    {
1085      pDelete(&aq);
1086    }
1087  }
1088  return result;
1089}
1090
1091poly pJet(poly p, int m)
1092{
1093  poly r=NULL;
1094  poly t=NULL;
1095
1096  while (p!=NULL)
1097  {
1098    if (pTotaldegree(p)<=m)
1099    {
1100      if (r==NULL)
1101        r=pHead(p);
1102      else
1103      if (t==NULL)
1104      {
1105        pNext(r)=pHead(p);
1106        t=pNext(r);
1107      }
1108      else
1109      {
1110        pNext(t)=pHead(p);
1111        pIter(t);
1112      }
1113    }
1114    pIter(p);
1115  }
1116  return r;
1117}
1118
1119poly pJetW(poly p, int m, short *w)
1120{
1121  poly r=NULL;
1122  poly t=NULL;
1123  short *wsave=ecartWeights;
1124
1125  ecartWeights=w;
1126
1127  while (p!=NULL)
1128  {
1129    if (totaldegreeWecart(p)<=m)
1130    {
1131      if (r==NULL)
1132        r=pHead(p);
1133      else
1134      if (t==NULL)
1135      {
1136        pNext(r)=pHead(p);
1137        t=pNext(r);
1138      }
1139      else
1140      {
1141        pNext(t)=pHead(p);
1142        pIter(t);
1143      }
1144    }
1145    pIter(p);
1146  }
1147  ecartWeights=wsave;
1148  return r;
1149}
1150
1151int pDegW(poly p, short *w)
1152{
1153  int r=0;
1154  short *wsave=ecartWeights;
1155
1156  ecartWeights=w;
1157
1158  while (p!=NULL)
1159  {
1160    r=max(r, totaldegreeWecart(p));
1161    pIter(p);
1162  }
1163  ecartWeights=wsave;
1164  return r;
1165}
1166
1167/*-----------type conversions ----------------------------*/
1168/*2
1169* input: a set of polys (len elements: p[0]..p[len-1])
1170* output: a vector
1171* p will not be changed
1172*/
1173poly  pPolys2Vec(polyset p, int len)
1174{
1175  poly v=NULL;
1176  poly h;
1177  int i;
1178
1179  for (i=len-1; i>=0; i--)
1180  {
1181    if (p[i])
1182    {
1183      h=pCopy(p[i]);
1184      pSetCompP(h,i+1);
1185      v=pAdd(v,h);
1186    }
1187  }
1188 return v;
1189}
1190
1191/*2
1192* convert a vector to a set of polys,
1193* allocates the polyset, (entries 0..(*len)-1)
1194* the vector will not be changed
1195*/
1196void  pVec2Polys(poly v, polyset *p, int *len)
1197{
1198  poly h;
1199  int k;
1200
1201  *len=pMaxComp(v);
1202  if (*len==0) *len=1;
1203  *p=(polyset)Alloc0((*len)*sizeof(poly));
1204  while (v!=NULL)
1205  {
1206    h=pHead(v);
1207    k=pGetComp(h);
1208    pSetComp(h,0);
1209    (*p)[k-1]=pAdd((*p)[k-1],h);
1210    pIter(v);
1211  }
1212}
1213
1214int pVar(poly m)
1215{
1216  if (m==NULL) return 0;
1217  if (pNext(m)!=NULL) return 0;
1218  int i,e=0;
1219  for (i=pVariables; i>0; i--)
1220  {
1221    if (pGetExp(m,i)==1)
1222    {
1223      if (e==0) e=i;
1224      else return 0;
1225    }
1226  }
1227  return e;
1228}
1229
1230/*----------utilities for syzygies--------------*/
1231//BOOLEAN   pVectorHasUnitM(poly p, int * k)
1232//{
1233//  while (p!=NULL)
1234//  {
1235//    if (pIsConstantComp(p))
1236//    {
1237//      *k = pGetComp(p);
1238//      return TRUE;
1239//    }
1240//    else pIter(p);
1241//  }
1242//  return FALSE;
1243//}
1244
1245BOOLEAN   pVectorHasUnitB(poly p, int * k)
1246{
1247  poly q=p,qq;
1248  int i;
1249
1250  while (q!=NULL)
1251  {
1252    if (pIsConstantComp(q))
1253    {
1254      i = pGetComp(q);
1255      qq = p;
1256      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1257      if (qq == q)
1258      {
1259        *k = i;
1260        return TRUE;
1261      }
1262      else
1263        pIter(q);
1264    }
1265    else pIter(q);
1266  }
1267  return FALSE;
1268}
1269
1270void   pVectorHasUnit(poly p, int * k, int * len)
1271{
1272  poly q=p,qq;
1273  int i,j=0;
1274
1275  *len = 0;
1276  while (q!=NULL)
1277  {
1278    if (pIsConstantComp(q))
1279    {
1280      i = pGetComp(q);
1281      qq = p;
1282      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1283      if (qq == q)
1284      {
1285       j = 0;
1286       while (qq!=NULL)
1287       {
1288         if (pGetComp(qq)==i) j++;
1289        pIter(qq);
1290       }
1291       if ((*len == 0) || (j<*len))
1292       {
1293         *len = j;
1294         *k = i;
1295       }
1296      }
1297    }
1298    pIter(q);
1299  }
1300}
1301
1302/*2
1303* returns TRUE if p1 = p2
1304*/
1305BOOLEAN pEqualPolys(poly p1,poly p2)
1306{
1307  while ((p1 != NULL) && (p2 != NULL))
1308  {
1309    /* p1 and p2 are non-NULL, so we may use pComp0 instead of pComp */
1310    if (! pEqual(p1, p2))
1311    {
1312       return FALSE;
1313    }
1314    if (nEqual(pGetCoeff(p1),pGetCoeff(p2)) == FALSE)
1315    {
1316      return FALSE;
1317    }
1318    pIter(p1);
1319    pIter(p2);
1320  }
1321  return (p1==p2);
1322}
1323
1324/*2
1325*returns TRUE if p1 is a skalar multiple of p2
1326*assume p1 != NULL and p2 != NULL
1327*/
1328BOOLEAN pComparePolys(poly p1,poly p2)
1329{
1330  number n,nn;
1331  int i;
1332
1333  if (pLength(p1) != pLength(p2))
1334    return FALSE;
1335  n=nDiv(pGetCoeff(p1),pGetCoeff(p2));
1336  while ((p1 != NULL) /*&& (p2 != NULL)*/)
1337  {
1338    if ( ! pLmEqual(p1, p2))
1339    {
1340        nDelete(&n);
1341        return FALSE;
1342    }
1343    if (!nEqual(pGetCoeff(p1),nn=nMult(pGetCoeff(p2),n)))
1344    {
1345      nDelete(&n);
1346      nDelete(&nn);
1347      return FALSE;
1348    }
1349    nDelete(&nn);
1350    pIter(p1);
1351    pIter(p2);
1352  }
1353  nDelete(&n);
1354  return TRUE;
1355}
Note: See TracBrowser for help on using the repository browser.