source: git/Singular/polys1.cc @ b98976

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