source: git/Singular/polys1.cc @ 97a7b44

fieker-DuValspielwiese
Last change on this file since 97a7b44 was 42af82e, checked in by Olaf Bachmann <obachman@…>, 25 years ago
* imap, and pOrdPolyInsertSetm do addition starting with the tail git-svn-id: file:///usr/local/Singular/svn/trunk@3058 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 23.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: polys1.cc,v 1.22 1999-05-26 14:16:50 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 "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#if 0
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#else
1012  while (p != NULL)
1013  {
1014    qq = p;
1015    pIter(p);
1016    qq->next = result;
1017    result = qq;
1018    pSetm(qq);
1019  }
1020  p = result;
1021  result = NULL;
1022  while (p != NULL)
1023  {
1024    qq = p;
1025    pIter(p);
1026    qq->next = NULL;
1027    result = pAdd(result, qq);
1028  }
1029  pTest(result);
1030#endif
1031  return result;
1032}
1033
1034/*2
1035*returns a re-ordered copy of a polynomial, with permutation of the variables
1036*/
1037poly pPermPoly (poly p, int * perm, ring oldRing,
1038   int *par_perm, int OldPar)
1039{
1040  int OldpVariables = oldRing->N;
1041  poly result = NULL;
1042  poly aq=NULL;
1043  poly qq;
1044  int i;
1045
1046  while (p != NULL)
1047  {
1048    if (OldPar==0)
1049    {
1050      qq = pInit();
1051      pGetCoeff(qq)=nMap(pGetCoeff(p));
1052    }
1053    else
1054    {
1055      qq=pOne();
1056      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar);
1057      pTest(aq);
1058    }
1059    pSetComp(qq, pRingGetComp(oldRing,p));
1060    if (nIsZero(pGetCoeff(qq)))
1061    {
1062      pDelete1(&qq);
1063    }
1064    else
1065    {
1066      for(i=1; i<=OldpVariables; i++)
1067      {
1068        if (pRingGetExp(oldRing,p,i)!=0)
1069        {
1070          if (perm==NULL)
1071          {
1072            pSetExp(qq,i, pRingGetExp(oldRing,p,i));
1073          }
1074          else if (perm[i]>0)
1075            pAddExp(qq,perm[i], pRingGetExp(oldRing, p,i));
1076          else if (perm[i]<0)
1077          {
1078            lnumber c=(lnumber)pGetCoeff(qq);
1079            c->z->e[-perm[i]-1]+=pRingGetExp(oldRing, p,i);
1080          }
1081          else
1082          {
1083            /* this variable maps to 0 !*/
1084            pDelete1(&qq);
1085            break;
1086          }
1087        }
1088      }
1089    }
1090    pIter(p);
1091    if (qq!=NULL)
1092    {
1093      pSetm(qq);
1094      pTest(qq);
1095      pTest(aq);
1096      if (aq!=NULL) qq=pMult(aq,qq);
1097      aq = qq;
1098      while (pNext(aq) != NULL) pIter(aq);
1099      pNext(aq) = result;
1100      aq = NULL;
1101      result = qq;
1102    }
1103    else if (aq!=NULL)
1104    {
1105      pDelete(&aq);
1106    }
1107  }
1108  p = result;
1109  result = NULL;
1110  while (p != NULL)
1111  {
1112    qq = p;
1113    pIter(p);
1114    qq->next = NULL;
1115    result = pAdd(result, qq);
1116  }
1117  pTest(result);
1118  return result;
1119}
1120
1121poly pJet(poly p, int m)
1122{
1123  poly r=NULL;
1124  poly t=NULL;
1125
1126  while (p!=NULL)
1127  {
1128    if (pTotaldegree(p)<=m)
1129    {
1130      if (r==NULL)
1131        r=pHead(p);
1132      else
1133      if (t==NULL)
1134      {
1135        pNext(r)=pHead(p);
1136        t=pNext(r);
1137      }
1138      else
1139      {
1140        pNext(t)=pHead(p);
1141        pIter(t);
1142      }
1143    }
1144    pIter(p);
1145  }
1146  return r;
1147}
1148
1149poly pJetW(poly p, int m, short *w)
1150{
1151  poly r=NULL;
1152  poly t=NULL;
1153  short *wsave=ecartWeights;
1154
1155  ecartWeights=w;
1156
1157  while (p!=NULL)
1158  {
1159    if (totaldegreeWecart(p)<=m)
1160    {
1161      if (r==NULL)
1162        r=pHead(p);
1163      else
1164      if (t==NULL)
1165      {
1166        pNext(r)=pHead(p);
1167        t=pNext(r);
1168      }
1169      else
1170      {
1171        pNext(t)=pHead(p);
1172        pIter(t);
1173      }
1174    }
1175    pIter(p);
1176  }
1177  ecartWeights=wsave;
1178  return r;
1179}
1180
1181int pDegW(poly p, short *w)
1182{
1183  int r=0;
1184  short *wsave=ecartWeights;
1185
1186  ecartWeights=w;
1187
1188  while (p!=NULL)
1189  {
1190    r=max(r, totaldegreeWecart(p));
1191    pIter(p);
1192  }
1193  ecartWeights=wsave;
1194  return r;
1195}
1196
1197/*-----------type conversions ----------------------------*/
1198/*2
1199* input: a set of polys (len elements: p[0]..p[len-1])
1200* output: a vector
1201* p will not be changed
1202*/
1203poly  pPolys2Vec(polyset p, int len)
1204{
1205  poly v=NULL;
1206  poly h;
1207  int i;
1208
1209  for (i=len-1; i>=0; i--)
1210  {
1211    if (p[i])
1212    {
1213      h=pCopy(p[i]);
1214      pSetCompP(h,i+1);
1215      v=pAdd(v,h);
1216    }
1217  }
1218 return v;
1219}
1220
1221/*2
1222* convert a vector to a set of polys,
1223* allocates the polyset, (entries 0..(*len)-1)
1224* the vector will not be changed
1225*/
1226void  pVec2Polys(poly v, polyset *p, int *len)
1227{
1228  poly h;
1229  int k;
1230
1231  *len=pMaxComp(v);
1232  if (*len==0) *len=1;
1233  *p=(polyset)Alloc0((*len)*sizeof(poly));
1234  while (v!=NULL)
1235  {
1236    h=pHead(v);
1237    k=pGetComp(h);
1238    pSetComp(h,0);
1239    (*p)[k-1]=pAdd((*p)[k-1],h);
1240    pIter(v);
1241  }
1242}
1243
1244int pVar(poly m)
1245{
1246  if (m==NULL) return 0;
1247  if (pNext(m)!=NULL) return 0;
1248  int i,e=0;
1249  for (i=pVariables; i>0; i--)
1250  {
1251    if (pGetExp(m,i)==1)
1252    {
1253      if (e==0) e=i;
1254      else return 0;
1255    }
1256  }
1257  return e;
1258}
1259
1260/*----------utilities for syzygies--------------*/
1261//BOOLEAN   pVectorHasUnitM(poly p, int * k)
1262//{
1263//  while (p!=NULL)
1264//  {
1265//    if (pIsConstantComp(p))
1266//    {
1267//      *k = pGetComp(p);
1268//      return TRUE;
1269//    }
1270//    else pIter(p);
1271//  }
1272//  return FALSE;
1273//}
1274
1275BOOLEAN   pVectorHasUnitB(poly p, int * k)
1276{
1277  poly q=p,qq;
1278  int i;
1279
1280  while (q!=NULL)
1281  {
1282    if (pIsConstantComp(q))
1283    {
1284      i = pGetComp(q);
1285      qq = p;
1286      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1287      if (qq == q)
1288      {
1289        *k = i;
1290        return TRUE;
1291      }
1292      else
1293        pIter(q);
1294    }
1295    else pIter(q);
1296  }
1297  return FALSE;
1298}
1299
1300void   pVectorHasUnit(poly p, int * k, int * len)
1301{
1302  poly q=p,qq;
1303  int i,j=0;
1304
1305  *len = 0;
1306  while (q!=NULL)
1307  {
1308    if (pIsConstantComp(q))
1309    {
1310      i = pGetComp(q);
1311      qq = p;
1312      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1313      if (qq == q)
1314      {
1315       j = 0;
1316       while (qq!=NULL)
1317       {
1318         if (pGetComp(qq)==i) j++;
1319        pIter(qq);
1320       }
1321       if ((*len == 0) || (j<*len))
1322       {
1323         *len = j;
1324         *k = i;
1325       }
1326      }
1327    }
1328    pIter(q);
1329  }
1330}
1331
1332/*2
1333* returns TRUE if p1 = p2
1334*/
1335BOOLEAN pEqualPolys(poly p1,poly p2)
1336{
1337  while ((p1 != NULL) && (p2 != NULL))
1338  {
1339    /* p1 and p2 are non-NULL, so we may use pComp0 instead of pComp */
1340    if (! pEqual(p1, p2))
1341    {
1342       return FALSE;
1343    }
1344    if (nEqual(pGetCoeff(p1),pGetCoeff(p2)) == FALSE)
1345    {
1346      return FALSE;
1347    }
1348    pIter(p1);
1349    pIter(p2);
1350  }
1351  return (p1==p2);
1352}
1353
1354/*2
1355*returns TRUE if p1 is a skalar multiple of p2
1356*assume p1 != NULL and p2 != NULL
1357*/
1358BOOLEAN pComparePolys(poly p1,poly p2)
1359{
1360  number n,nn;
1361  int i;
1362
1363  if (pLength(p1) != pLength(p2))
1364    return FALSE;
1365  n=nDiv(pGetCoeff(p1),pGetCoeff(p2));
1366  while ((p1 != NULL) /*&& (p2 != NULL)*/)
1367  {
1368    if ( ! pLmEqual(p1, p2))
1369    {
1370        nDelete(&n);
1371        return FALSE;
1372    }
1373    if (!nEqual(pGetCoeff(p1),nn=nMult(pGetCoeff(p2),n)))
1374    {
1375      nDelete(&n);
1376      nDelete(&nn);
1377      return FALSE;
1378    }
1379    nDelete(&nn);
1380    pIter(p1);
1381    pIter(p2);
1382  }
1383  nDelete(&n);
1384  return TRUE;
1385}
Note: See TracBrowser for help on using the repository browser.