source: git/Singular/polys1.cc @ 6e56de

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