source: git/Singular/polys1.cc @ 9d72fe

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