source: git/kernel/polys1.cc @ 3c6379

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