source: git/kernel/polys1.cc @ b130fb

spielwiese
Last change on this file since b130fb was b130fb, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: gcc-3.4 git-svn-id: file:///usr/local/Singular/svn/trunk@7733 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 32.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: polys1.cc,v 1.8 2005-02-17 09:42:22 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
923void pCleardenom_n(poly ph,number &c)
924{
925  number d, h;
926  poly p;
927
928  p = ph;
929  if(pNext(p)==NULL)
930  {
931    c=nInvers(pGetCoeff(p));
932    pSetCoeff(p,nInit(1));
933  }
934  else
935  {
936    h = nInit(1);
937    while (p!=NULL)
938    {
939      nNormalize(pGetCoeff(p));
940      d=nLcm(h,pGetCoeff(p),currRing);
941      nDelete(&h);
942      h=d;
943      pIter(p);
944    }
945    c=h;
946    /* contains the 1/lcm of all denominators */
947    if(!nIsOne(h))
948    {
949      p = ph;
950      while (p!=NULL)
951      {
952        /* should be:
953        * number hh;
954        * nGetDenom(p->coef,&hh);
955        * nMult(&h,&hh,&d);
956        * nNormalize(d);
957        * nDelete(&hh);
958        * nMult(d,p->coef,&hh);
959        * nDelete(&d);
960        * nDelete(&(p->coef));
961        * p->coef =hh;
962        */
963        d=nMult(h,pGetCoeff(p));
964        nNormalize(d);
965        pSetCoeff(p,d);
966        pIter(p);
967      }
968      if (nGetChar()==1)
969      {
970        loop
971        {
972          h = nInit(1);
973          p=ph;
974          while (p!=NULL)
975          {
976            d=nLcm(h,pGetCoeff(p),currRing);
977            nDelete(&h);
978            h=d;
979            pIter(p);
980          }
981          /* contains the 1/lcm of all denominators */
982          if(!nIsOne(h))
983          {
984            p = ph;
985            while (p!=NULL)
986            {
987              /* should be:
988              * number hh;
989              * nGetDenom(p->coef,&hh);
990              * nMult(&h,&hh,&d);
991              * nNormalize(d);
992              * nDelete(&hh);
993              * nMult(d,p->coef,&hh);
994              * nDelete(&d);
995              * nDelete(&(p->coef));
996              * p->coef =hh;
997              */
998              d=nMult(h,pGetCoeff(p));
999              nNormalize(d);
1000              pSetCoeff(p,d);
1001              pIter(p);
1002            }
1003            number t=nMult(c,h);
1004            nDelete(&c);
1005            c=t;
1006          }
1007          else
1008          {
1009            break;
1010          }
1011          nDelete(&h);
1012        }
1013      }
1014    }
1015  }
1016}
1017
1018/*2
1019*tests if p is homogeneous with respect to the actual weigths
1020*/
1021BOOLEAN pIsHomogeneous (poly p)
1022{
1023  poly qp=p;
1024  int o;
1025
1026  if ((p == NULL) || (pNext(p) == NULL)) return TRUE;
1027  pFDegProc d=(pLexOrder ? pTotaldegree : pFDeg );
1028  o = d(p,currRing);
1029  do
1030  {
1031    if (d(qp,currRing) != o) return FALSE;
1032    pIter(qp);
1033  }
1034  while (qp != NULL);
1035  return TRUE;
1036}
1037
1038// orders monoms of poly using merge sort (ususally faster than
1039// insertion sort). ASSUMES that pSetm was performed on monoms
1040poly pOrdPolyMerge(poly p)
1041{
1042  poly qq,pp,result=NULL;
1043
1044  if (p == NULL) return NULL;
1045
1046  loop
1047  {
1048    qq = p;
1049    loop
1050    {
1051      if (pNext(p) == NULL)
1052      {
1053        result=pAdd(result, qq);
1054        pTest(result);
1055        return result;
1056      }
1057      if (pLmCmp(p,pNext(p)) != 1)
1058      {
1059        pp = p;
1060        pIter(p);
1061        pNext(pp) = NULL;
1062        result = pAdd(result, qq);
1063        break;
1064      }
1065      pIter(p);
1066    }
1067  }
1068}
1069
1070// orders monoms of poly using insertion sort, performs pSetm on each monom
1071poly pOrdPolyInsertSetm(poly p)
1072{
1073  poly qq,result = NULL;
1074
1075#if 0
1076  while (p != NULL)
1077  {
1078    qq = p;
1079    pIter(p);
1080    qq->next = NULL;
1081    pSetm(qq);
1082    result = pAdd(result,qq);
1083    pTest(result);
1084  }
1085#else
1086  while (p != NULL)
1087  {
1088    qq = p;
1089    pIter(p);
1090    qq->next = result;
1091    result = qq;
1092    pSetm(qq);
1093  }
1094  p = result;
1095  result = NULL;
1096  while (p != NULL)
1097  {
1098    qq = p;
1099    pIter(p);
1100    qq->next = NULL;
1101    result = pAdd(result, qq);
1102  }
1103  pTest(result);
1104#endif
1105  return result;
1106}
1107
1108/*2
1109*returns a re-ordered copy of a polynomial, with permutation of the variables
1110*/
1111poly pPermPoly (poly p, int * perm, ring oldRing, nMapFunc nMap,
1112   int *par_perm, int OldPar)
1113{
1114  int OldpVariables = oldRing->N;
1115  poly result = NULL;
1116  poly result_last = NULL;
1117  poly aq=NULL; /* the map coefficient */
1118  poly qq; /* the mapped monomial */
1119
1120  while (p != NULL)
1121  {
1122    if ((OldPar==0)||(rField_is_GF(oldRing)))
1123    {
1124      qq = pInit();
1125      number n=nMap(pGetCoeff(p));
1126      if ((currRing->minpoly!=NULL)
1127      && ((rField_is_Zp_a()) || (rField_is_Q_a())))
1128      {
1129        nNormalize(n);
1130      }
1131      pGetCoeff(qq)=n;
1132    // coef may be zero:  pTest(qq);
1133    }
1134    else
1135    {
1136      qq=pOne();
1137      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar, oldRing);
1138      if ((currRing->minpoly!=NULL)
1139      && ((rField_is_Zp_a()) || (rField_is_Q_a())))
1140      {
1141        poly tmp=aq;
1142        while (tmp!=NULL)
1143        {
1144          number n=pGetCoeff(tmp);
1145          nNormalize(n);
1146          pGetCoeff(tmp)=n;
1147          pIter(tmp);
1148        }
1149      }
1150      pTest(aq);
1151    }
1152    if (rRing_has_Comp(currRing)) pSetComp(qq, p_GetComp(p,oldRing));
1153    if (nIsZero(pGetCoeff(qq)))
1154    {
1155      pDeleteLm(&qq);
1156    }
1157    else
1158    {
1159      int i;
1160      int mapped_to_par=0;
1161      for(i=1; i<=OldpVariables; i++)
1162      {
1163        int e=p_GetExp(p,i,oldRing);
1164        if (e!=0)
1165        {
1166          if (perm==NULL)
1167          {
1168            pSetExp(qq,i, e);
1169          }
1170          else if (perm[i]>0)
1171            pAddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/);
1172          else if (perm[i]<0)
1173          {
1174            if (rField_is_GF())
1175            {
1176              number c=pGetCoeff(qq);
1177              number ee=nfPar(1);
1178              number eee;nfPower(ee,e,&eee); //nfDelete(ee,currRing);
1179              ee=nfMult(c,eee);
1180              //nfDelete(c,currRing);nfDelete(eee,currRing);
1181              pSetCoeff0(qq,ee);
1182            }
1183            else
1184            {
1185              lnumber c=(lnumber)pGetCoeff(qq);
1186              if (c->z->next==NULL)
1187                napAddExp(c->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
1188              else /* more difficult: we have really to multiply: */
1189              {
1190                lnumber mmc=(lnumber)naInit(1);
1191                napSetExp(mmc->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
1192                napSetm(mmc->z);
1193                pGetCoeff(qq)=naMult((number)c,(number)mmc);
1194                nDelete((number *)&c);
1195                nDelete((number *)&mmc); 
1196              }
1197              mapped_to_par=1;
1198            }
1199          }
1200          else
1201          {
1202            /* this variable maps to 0 !*/
1203            pDeleteLm(&qq);
1204            break;
1205          }
1206        }
1207      }
1208      if (mapped_to_par
1209      && (currRing->minpoly!=NULL))
1210      {
1211        number n=pGetCoeff(qq);
1212        nNormalize(n);
1213        pGetCoeff(qq)=n;
1214      }
1215    }
1216    pIter(p);
1217#if 1
1218    if (qq!=NULL)
1219    {
1220      pSetm(qq);
1221      pTest(aq);
1222      pTest(qq);
1223      if (aq!=NULL) qq=pMult(aq,qq);
1224      aq = qq;
1225      while (pNext(aq) != NULL) pIter(aq);
1226      if (result_last==NULL)
1227      {
1228        result=qq;
1229      }
1230      else
1231      {
1232        pNext(result_last)=qq;
1233      }
1234      result_last=aq;
1235      aq = NULL;
1236    }
1237    else if (aq!=NULL)
1238    {
1239      pDelete(&aq);
1240    }
1241  }
1242  result=pSortAdd(result);
1243#else
1244  //  if (qq!=NULL)
1245  //  {
1246  //    pSetm(qq);
1247  //    pTest(qq);
1248  //    pTest(aq);
1249  //    if (aq!=NULL) qq=pMult(aq,qq);
1250  //    aq = qq;
1251  //    while (pNext(aq) != NULL) pIter(aq);
1252  //    pNext(aq) = result;
1253  //    aq = NULL;
1254  //    result = qq;
1255  //  }
1256  //  else if (aq!=NULL)
1257  //  {
1258  //    pDelete(&aq);
1259  //  }
1260  //}
1261  //p = result;
1262  //result = NULL;
1263  //while (p != NULL)
1264  //{
1265  //  qq = p;
1266  //  pIter(p);
1267  //  qq->next = NULL;
1268  //  result = pAdd(result, qq);
1269  //}
1270#endif
1271  pTest(result);
1272  return result;
1273}
1274
1275#if 0
1276/*2
1277*returns a re-ordered copy of a polynomial, with permutation of the variables
1278*/
1279poly p_PermPoly (poly p, int * perm, ring oldRing,
1280   int *par_perm, int OldPar, ring newRing)
1281{
1282  int OldpVariables = oldRing->N;
1283  poly result = NULL;
1284  poly result_last = NULL;
1285  poly aq=NULL; /* the map coefficient */
1286  poly qq; /* the mapped monomial */
1287
1288  while (p != NULL)
1289  {
1290    if (OldPar==0)
1291    {
1292      qq = pInit();
1293      number n=newRing->cf->nMap(pGetCoeff(p));
1294      if ((newRing->minpoly!=NULL)
1295      && ((rField_is_Zp_a(newRing)) || (rField_is_Q_a(newRing))))
1296      {
1297        newRing->cf->nNormalize(n);
1298      }
1299      pGetCoeff(qq)=n;
1300    // coef may be zero:  pTest(qq);
1301    }
1302    else
1303    {
1304      qq=p_ISet(1, newRing);
1305      aq=naPermNumber(pGetCoeff(p),par_perm,OldPar, oldRing);
1306      if ((newRing->minpoly!=NULL)
1307      && ((rField_is_Zp_a(newRing)) || (rField_is_Q_a(newRing))))
1308      {
1309        poly tmp=aq;
1310        while (tmp!=NULL)
1311        {
1312          number n=pGetCoeff(tmp);
1313          newRing->cf->nNormalize(n);
1314          pGetCoeff(tmp)=n;
1315          pIter(tmp);
1316        }
1317      }
1318      //pTest(aq);
1319    }
1320    p_SetComp(qq, p_GetComp(p,oldRing), newRing);
1321    if (newRing->cf->nIsZero(pGetCoeff(qq)))
1322    {
1323      p_DeleteLm(&qq, newRing);
1324    }
1325    else
1326    {
1327      int i;
1328      int mapped_to_par=0;
1329      for(i=1; i<=OldpVariables; i++)
1330      {
1331        int e=p_GetExp(p,i,oldRing);
1332        if (e!=0)
1333        {
1334          if (perm==NULL)
1335          {
1336            p_SetExp(qq,i, e, newRing);
1337          }
1338          else if (perm[i]>0)
1339            p_AddExp(qq,perm[i], e/*p_GetExp( p,i,oldRing)*/, newRing);
1340          else if (perm[i]<0)
1341          {
1342            lnumber c=(lnumber)pGetCoeff(qq);
1343            napAddExp(c->z,-perm[i],e/*p_GetExp( p,i,oldRing)*/);
1344            mapped_to_par=1;
1345          }
1346          else
1347          {
1348            /* this variable maps to 0 !*/
1349            p_DeleteLm(&qq, newRing);
1350            break;
1351          }
1352        }
1353      }
1354      if (mapped_to_par
1355      && (newRing->minpoly!=NULL))
1356      {
1357        number n=pGetCoeff(qq);
1358        newRing->cf->nNormalize(n);
1359        pGetCoeff(qq)=n;
1360      }
1361    }
1362    pIter(p);
1363    if (qq!=NULL)
1364    {
1365      p_Setm(qq, newRing);
1366      //pTest(aq);
1367      //pTest(qq);
1368      if (aq!=NULL) qq=pMult(aq,qq);
1369      aq = qq;
1370      while (pNext(aq) != NULL) pIter(aq);
1371      if (result_last==NULL)
1372      {
1373        result=qq;
1374      }
1375      else
1376      {
1377        pNext(result_last)=qq;
1378      }
1379      result_last=aq;
1380      aq = NULL;
1381    }
1382    else if (aq!=NULL)
1383    {
1384      p_Delete(&aq, newRing);
1385    }
1386  }
1387  result=pOrdPolyMerge(result);
1388  //pTest(result);
1389  return result;
1390}
1391#endif
1392
1393poly ppJet(poly p, int m)
1394{
1395  poly r=NULL;
1396  poly t=NULL;
1397
1398  while (p!=NULL)
1399  {
1400    if (pTotaldegree(p)<=m)
1401    {
1402      if (r==NULL)
1403        r=pHead(p);
1404      else
1405      if (t==NULL)
1406      {
1407        pNext(r)=pHead(p);
1408        t=pNext(r);
1409      }
1410      else
1411      {
1412        pNext(t)=pHead(p);
1413        pIter(t);
1414      }
1415    }
1416    pIter(p);
1417  }
1418  return r;
1419}
1420
1421poly pJet(poly p, int m)
1422{
1423  poly t=NULL;
1424
1425  while((p!=NULL) && (pTotaldegree(p)>m)) pLmDelete(&p);
1426  if (p==NULL) return NULL;
1427  poly r=p;
1428  while (pNext(p)!=NULL)
1429  {
1430    if (pTotaldegree(pNext(p))>m)
1431    {
1432      pLmDelete(&pNext(p));
1433    }
1434    else
1435      pIter(p);
1436  }
1437  return r;
1438}
1439
1440poly ppJetW(poly p, int m, short *w)
1441{
1442  poly r=NULL;
1443  poly t=NULL;
1444  while (p!=NULL)
1445  {
1446    if (totaldegreeWecart_IV(p,currRing,w)<=m)
1447    {
1448      if (r==NULL)
1449        r=pHead(p);
1450      else
1451      if (t==NULL)
1452      {
1453        pNext(r)=pHead(p);
1454        t=pNext(r);
1455      }
1456      else
1457      {
1458        pNext(t)=pHead(p);
1459        pIter(t);
1460      }
1461    }
1462    pIter(p);
1463  }
1464  return r;
1465}
1466
1467poly pJetW(poly p, int m, short *w)
1468{
1469  poly t=NULL;
1470  while((p!=NULL) && (totaldegreeWecart_IV(p,currRing,w)>m)) pLmDelete(&p);
1471  if (p==NULL) return NULL;
1472  poly r=p;
1473  while (pNext(p)!=NULL)
1474  {
1475    if (totaldegreeWecart_IV(pNext(p),currRing,w)>m)
1476    {
1477      pLmDelete(&pNext(p));
1478    }
1479    else
1480      pIter(p);
1481  }
1482  return r;
1483}
1484
1485int pMinDeg(poly p,intvec *w)
1486{
1487  if(p==NULL)
1488    return -1;
1489  int d=-1;
1490  while(p!=NULL)
1491  {
1492    int d0=0;
1493    for(int j=0;j<pVariables;j++)
1494      if(w==NULL||j>=w->length())
1495        d0+=pGetExp(p,j+1);
1496      else
1497        d0+=(*w)[j]*pGetExp(p,j+1);
1498    if(d0<d||d==-1)
1499      d=d0;
1500    pIter(p);
1501  }
1502  return d;
1503}
1504
1505poly pSeries(int n,poly p,poly u, intvec *w)
1506{
1507  short *ww=iv2array(w);
1508  if(p!=NULL)
1509  {
1510    if(u==NULL)
1511      p=pJetW(p,n,ww);
1512    else
1513      p=pJetW(pMult(p,pInvers(n-pMinDeg(p,w),u,w)),n,ww);
1514  }
1515  omFreeSize((ADDRESS)ww,(pVariables+1)*sizeof(short));
1516  return p;
1517}
1518
1519poly pInvers(int n,poly u,intvec *w)
1520{
1521  short *ww=iv2array(w);
1522  if(n<0)
1523    return NULL;
1524  number u0=nInvers(pGetCoeff(u));
1525  poly v=pNSet(u0);
1526  if(n==0)
1527    return v;
1528  poly u1=pJetW(pSub(pOne(),pMult_nn(u,u0)),n,ww);
1529  if(u1==NULL)
1530    return v;
1531  poly v1=pMult_nn(pCopy(u1),u0);
1532  v=pAdd(v,pCopy(v1));
1533  for(int i=n/pMinDeg(u1,w);i>1;i--)
1534  {
1535    v1=pJetW(pMult(v1,pCopy(u1)),n,ww);
1536    v=pAdd(v,pCopy(v1));
1537  }
1538  pDelete(&u1);
1539  pDelete(&v1);
1540  omFreeSize((ADDRESS)ww,(pVariables+1)*sizeof(short));
1541  return v;
1542}
1543
1544long pDegW(poly p, short *w)
1545{
1546  long r=-LONG_MAX;
1547
1548  while (p!=NULL)
1549  {
1550    r=si_max(r, totaldegreeWecart_IV(p,currRing,w));
1551    pIter(p);
1552  }
1553  return r;
1554}
1555
1556/*-----------type conversions ----------------------------*/
1557/*2
1558* input: a set of polys (len elements: p[0]..p[len-1])
1559* output: a vector
1560* p will not be changed
1561*/
1562poly  pPolys2Vec(polyset p, int len)
1563{
1564  poly v=NULL;
1565  poly h;
1566  int i;
1567
1568  for (i=len-1; i>=0; i--)
1569  {
1570    if (p[i])
1571    {
1572      h=pCopy(p[i]);
1573      pSetCompP(h,i+1);
1574      v=pAdd(v,h);
1575    }
1576  }
1577 return v;
1578}
1579
1580/*2
1581* convert a vector to a set of polys,
1582* allocates the polyset, (entries 0..(*len)-1)
1583* the vector will not be changed
1584*/
1585void  pVec2Polys(poly v, polyset *p, int *len)
1586{
1587  poly h;
1588  int k;
1589
1590  *len=pMaxComp(v);
1591  if (*len==0) *len=1;
1592  *p=(polyset)omAlloc0((*len)*sizeof(poly));
1593  while (v!=NULL)
1594  {
1595    h=pHead(v);
1596    k=pGetComp(h);
1597    pSetComp(h,0);
1598    (*p)[k-1]=pAdd((*p)[k-1],h);
1599    pIter(v);
1600  }
1601}
1602
1603int pVar(poly m)
1604{
1605  if (m==NULL) return 0;
1606  if (pNext(m)!=NULL) return 0;
1607  int i,e=0;
1608  for (i=pVariables; i>0; i--)
1609  {
1610    if (pGetExp(m,i)==1)
1611    {
1612      if (e==0) e=i;
1613      else return 0;
1614    }
1615  }
1616  return e;
1617}
1618
1619/*----------utilities for syzygies--------------*/
1620//BOOLEAN   pVectorHasUnitM(poly p, int * k)
1621//{
1622//  while (p!=NULL)
1623//  {
1624//    if (pLmIsConstantComp(p))
1625//    {
1626//      *k = pGetComp(p);
1627//      return TRUE;
1628//    }
1629//    else pIter(p);
1630//  }
1631//  return FALSE;
1632//}
1633
1634BOOLEAN   pVectorHasUnitB(poly p, int * k)
1635{
1636  poly q=p,qq;
1637  int i;
1638
1639  while (q!=NULL)
1640  {
1641    if (pLmIsConstantComp(q))
1642    {
1643      i = pGetComp(q);
1644      qq = p;
1645      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1646      if (qq == q)
1647      {
1648        *k = i;
1649        return TRUE;
1650      }
1651      else
1652        pIter(q);
1653    }
1654    else pIter(q);
1655  }
1656  return FALSE;
1657}
1658
1659void   pVectorHasUnit(poly p, int * k, int * len)
1660{
1661  poly q=p,qq;
1662  int i,j=0;
1663
1664  *len = 0;
1665  while (q!=NULL)
1666  {
1667    if (pLmIsConstantComp(q))
1668    {
1669      i = pGetComp(q);
1670      qq = p;
1671      while ((qq != q) && (pGetComp(qq) != i)) pIter(qq);
1672      if (qq == q)
1673      {
1674       j = 0;
1675       while (qq!=NULL)
1676       {
1677         if (pGetComp(qq)==i) j++;
1678        pIter(qq);
1679       }
1680       if ((*len == 0) || (j<*len))
1681       {
1682         *len = j;
1683         *k = i;
1684       }
1685      }
1686    }
1687    pIter(q);
1688  }
1689}
1690
1691/*2
1692* returns TRUE if p1 = p2
1693*/
1694BOOLEAN pEqualPolys(poly p1,poly p2)
1695{
1696  while ((p1 != NULL) && (p2 != NULL))
1697  {
1698    if (! pLmEqual(p1, p2))
1699      return FALSE;
1700    if (! nEqual(pGetCoeff(p1), pGetCoeff(p2)))
1701      return FALSE;
1702    pIter(p1);
1703    pIter(p2);
1704  }
1705  return (p1==p2);
1706}
1707
1708/*2
1709*returns TRUE if p1 is a skalar multiple of p2
1710*assume p1 != NULL and p2 != NULL
1711*/
1712BOOLEAN pComparePolys(poly p1,poly p2)
1713{
1714  number n,nn;
1715  int i;
1716  pAssume(p1 != NULL && p2 != NULL);
1717
1718  if (!pLmEqual(p1,p2)) //compare leading mons
1719      return FALSE;
1720  if ((pNext(p1)==NULL) && (pNext(p2)!=NULL))
1721     return FALSE;
1722  if ((pNext(p2)==NULL) && (pNext(p1)!=NULL))
1723     return FALSE;
1724  if (pLength(p1) != pLength(p2))
1725    return FALSE;
1726  n=nDiv(pGetCoeff(p1),pGetCoeff(p2));
1727  while ((p1 != NULL) /*&& (p2 != NULL)*/)
1728  {
1729    if ( ! pLmEqual(p1, p2))
1730    {
1731        nDelete(&n);
1732        return FALSE;
1733    }
1734    if (!nEqual(pGetCoeff(p1),nn=nMult(pGetCoeff(p2),n)))
1735    {
1736      nDelete(&n);
1737      nDelete(&nn);
1738      return FALSE;
1739    }
1740    nDelete(&nn);
1741    pIter(p1);
1742    pIter(p2);
1743  }
1744  nDelete(&n);
1745  return TRUE;
1746}
Note: See TracBrowser for help on using the repository browser.