source: git/kernel/polys.cc @ c3005ea

spielwiese
Last change on this file since c3005ea was c3005ea, checked in by Hans Schoenemann <hannes@…>, 13 years ago
code cleanup: pDehomgen, idDehomgen git-svn-id: file:///usr/local/Singular/svn/trunk@13700 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 21.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5
6/*
7* ABSTRACT - all basic methods to manipulate polynomials
8*/
9
10/* includes */
11#include <stdio.h>
12#include <string.h>
13#include <ctype.h>
14#include <kernel/mod2.h>
15#include <kernel/options.h>
16#include <omalloc/omalloc.h>
17#include <kernel/febase.h>
18#include <kernel/numbers.h>
19#include <kernel/polys.h>
20#include <kernel/ring.h>
21#include <kernel/sbuckets.h>
22
23#ifdef HAVE_PLURAL
24#include <kernel/gring.h>
25#include <kernel/sca.h>
26#endif
27
28/* ----------- global variables, set by pSetGlobals --------------------- */
29/* computes length and maximal degree of a POLYnomial */
30pLDegProc pLDeg;
31/* computes the degree of the initial term, used for std */
32pFDegProc pFDeg;
33/* the monomial ordering of the head monomials a and b */
34/* returns -1 if a comes before b, 0 if a=b, 1 otherwise */
35
36int pVariables;     // number of variables
37
38/* 1 for polynomial ring, -1 otherwise */
39int     pOrdSgn;
40// it is of type int, not BOOLEAN because it is also in ip
41/* TRUE if the monomial ordering is not compatible with pFDeg */
42BOOLEAN pLexOrder;
43
44/* ----------- global variables, set by procedures from hecke/kstd1 ----- */
45/* the highest monomial below pHEdge */
46poly      ppNoether = NULL;
47
48/* -------------------------------------------------------- */
49/*2
50* change all global variables to fit the description of the new ring
51*/
52
53
54void pSetGlobals(const ring r, BOOLEAN complete)
55{
56  int i;
57  if (ppNoether!=NULL) pDelete(&ppNoether);
58  pVariables = r->N;
59  pOrdSgn = r->OrdSgn;
60  pFDeg=r->pFDeg;
61  pLDeg=r->pLDeg;
62  pLexOrder=r->LexOrder;
63
64  if (complete)
65  {
66    test &= ~ TEST_RINGDEP_OPTS;
67    test |= r->options;
68  }
69}
70
71// resets the pFDeg and pLDeg: if pLDeg is not given, it is
72// set to currRing->pLDegOrig, i.e. to the respective LDegProc which
73// only uses pFDeg (and not pDeg, or pTotalDegree, etc)
74void pSetDegProcs(pFDegProc new_FDeg, pLDegProc new_lDeg)
75{
76  assume(new_FDeg != NULL);
77  pFDeg = new_FDeg;
78  currRing->pFDeg = new_FDeg;
79
80  if (new_lDeg == NULL)
81    new_lDeg = currRing->pLDegOrig;
82
83  pLDeg = new_lDeg;
84  currRing->pLDeg = new_lDeg;
85}
86
87
88// restores pFDeg and pLDeg:
89extern void pRestoreDegProcs(pFDegProc old_FDeg, pLDegProc old_lDeg)
90{
91  assume(old_FDeg != NULL && old_lDeg != NULL);
92  pFDeg = old_FDeg;
93  currRing->pFDeg = old_FDeg;
94  pLDeg = old_lDeg;
95  currRing->pLDeg = old_lDeg;
96}
97
98/*2
99* assumes that the head term of b is a multiple of the head term of a
100* and return the multiplicant *m
101*/
102poly pDivide(poly a, poly b)
103{
104  int i;
105  poly result = pInit();
106
107  for(i=(int)pVariables; i; i--)
108    pSetExp(result,i, pGetExp(a,i)- pGetExp(b,i));
109  pSetComp(result, pGetComp(a) - pGetComp(b));
110  pSetm(result);
111  return result;
112}
113
114#ifdef HAVE_RINGS   //TODO Oliver
115#define pDiv_nn(p, n)              p_Div_nn(p, n, currRing)
116
117poly p_Div_nn(poly p, const number n, const ring r)
118{
119  pAssume(!n_IsZero(n,r));
120  p_Test(p, r);
121
122  poly q = p;
123  while (p != NULL)
124  {
125    number nc = pGetCoeff(p);
126    pSetCoeff0(p, n_Div(nc, n, r));
127    n_Delete(&nc, r);
128    pIter(p);
129  }
130  p_Test(q, r);
131  return q;
132}
133#endif
134
135/*2
136* divides a by the monomial b, ignores monomials which are not divisible
137* assumes that b is not NULL
138*/
139poly pDivideM(poly a, poly b)
140{
141  if (a==NULL) return NULL;
142  poly result=a;
143  poly prev=NULL;
144  int i;
145#ifdef HAVE_RINGS
146  number inv=pGetCoeff(b);
147#else
148  number inv=nInvers(pGetCoeff(b));
149#endif
150
151  while (a!=NULL)
152  {
153    if (pDivisibleBy(b,a))
154    {
155      for(i=(int)pVariables; i; i--)
156         pSubExp(a,i, pGetExp(b,i));
157      pSubComp(a, pGetComp(b));
158      pSetm(a);
159      prev=a;
160      pIter(a);
161    }
162    else
163    {
164      if (prev==NULL)
165      {
166        pLmDelete(&result);
167        a=result;
168      }
169      else
170      {
171        pLmDelete(&pNext(prev));
172        a=pNext(prev);
173      }
174    }
175  }
176#ifdef HAVE_RINGS
177  if (nIsUnit(inv))
178  {
179    inv = nInvers(inv);
180    pMult_nn(result,inv);
181    nDelete(&inv);
182  }
183  else
184  {
185    pDiv_nn(result,inv);
186  }
187#else
188  pMult_nn(result,inv);
189  nDelete(&inv);
190#endif
191  pDelete(&b);
192  return result;
193}
194
195/*2
196* returns the LCM of the head terms of a and b in *m
197*/
198void pLcm(poly a, poly b, poly m)
199{
200  int i;
201  for (i=pVariables; i; i--)
202  {
203    pSetExp(m,i, si_max( pGetExp(a,i), pGetExp(b,i)));
204  }
205  pSetComp(m, si_max(pGetComp(a), pGetComp(b)));
206  /* Don't do a pSetm here, otherwise hres/lres chockes */
207}
208
209/*2
210* convert monomial given as string to poly, e.g. 1x3y5z
211*/
212const char * p_Read(const char *st, poly &rc, const ring r)
213{
214  if (r==NULL) { rc=NULL;return st;}
215  int i,j;
216  rc = p_Init(r);
217  const char *s = r->cf->nRead(st,&(rc->coef));
218  if (s==st)
219  /* i.e. it does not start with a coeff: test if it is a ringvar*/
220  {
221    j = r_IsRingVar(s,r);
222    if (j >= 0)
223    {
224      p_IncrExp(rc,1+j,r);
225      while (*s!='\0') s++;
226      goto done;
227    }
228  }
229  while (*s!='\0')
230  {
231    char ss[2];
232    ss[0] = *s++;
233    ss[1] = '\0';
234    j = r_IsRingVar(ss,r);
235    if (j >= 0)
236    {
237      const char *s_save=s;
238      s = eati(s,&i);
239      if (((unsigned long)i) >  r->bitmask)
240      {
241        // exponent to large: it is not a monomial
242        p_LmDelete(&rc,r);
243        return s_save;
244      }
245      p_AddExp(rc,1+j, (long)i, r);
246    }
247    else
248    {
249      // 1st char of is not a varname
250      p_LmDelete(&rc,r);
251      s--;
252      return s;
253    }
254  }
255done:
256  if (r->cf->nIsZero(pGetCoeff(rc))) p_LmDelete(&rc,r);
257  else
258  {
259#ifdef HAVE_PLURAL
260    // in super-commutative ring
261    // squares of anti-commutative variables are zeroes!
262    if(rIsSCA(r))
263    {
264      const unsigned int iFirstAltVar = scaFirstAltVar(r);
265      const unsigned int iLastAltVar  = scaLastAltVar(r);
266
267      assume(rc != NULL);
268
269      for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
270        if( p_GetExp(rc, k, r) > 1 )
271        {
272          p_LmDelete(&rc, r);
273          goto finish;
274        }
275    }
276#endif
277    p_Setm(rc,r);
278  }
279finish:
280  return s;
281}
282
283BOOLEAN _p_Test(poly p, ring r, int level);
284poly pmInit(const char *st, BOOLEAN &ok)
285{
286  poly p;
287  const char *s=p_Read(st,p,currRing);
288  if (*s!='\0')
289  {
290    if ((s!=st)&&isdigit(st[0]))
291    {
292      errorreported=TRUE;
293    }
294    ok=FALSE;
295    pDelete(&p);
296    return NULL;
297  }
298  #ifdef PDEBUG
299  _p_Test(p,currRing,PDEBUG);
300  #endif
301  ok=!errorreported;
302  return p;
303}
304
305/*2
306*make p homogeneous by multiplying the monomials by powers of x_varnum
307*assume: deg(var(varnum))==1
308*/
309poly pHomogen (poly p, int varnum)
310{
311  pFDegProc deg;
312  if (pLexOrder && (currRing->order[0]==ringorder_lp))
313    deg=p_Totaldegree;
314  else
315    deg=pFDeg;
316
317  poly q=NULL, qn;
318  int  o,ii;
319  sBucket_pt bp;
320
321  if (p!=NULL)
322  {
323    if ((varnum < 1) || (varnum > pVariables))
324    {
325      return NULL;
326    }
327    o=deg(p,currRing);
328    q=pNext(p);
329    while (q != NULL)
330    {
331      ii=deg(q,currRing);
332      if (ii>o) o=ii;
333      pIter(q);
334    }
335    q = pCopy(p);
336    bp = sBucketCreate(currRing);
337    while (q != NULL)
338    {
339      ii = o-deg(q,currRing);
340      if (ii!=0)
341      {
342        pAddExp(q,varnum, (long)ii);
343        pSetm(q);
344      }
345      qn = pNext(q);
346      pNext(q) = NULL;
347      sBucket_Add_p(bp, q, 1);
348      q = qn;
349    }
350    sBucketDestroyAdd(bp, &q, &ii);
351  }
352  return q;
353}
354
355/*4
356*Returns the exponent of the maximal power of the leading monomial of
357*p2 in that of p1
358*/
359static int pGetMaxPower (poly p1,poly p2)
360{
361  int     i,k,res = 32000; /*a very large integer*/
362
363  if (p1 == NULL) return 0;
364  for (i=1; i<=pVariables; i++)
365  {
366    if ( pGetExp(p2,i) != 0)
367    {
368      k =  pGetExp(p1,i) /  pGetExp(p2,i);
369      if (k < res) res = k;
370    }
371  }
372  return res;
373}
374
375/*2
376*returns the leading monomial of p1 divided by the maximal power of that
377*of p2
378*/
379static poly pDivByMonom (poly p1,poly p2)
380{
381  int     k, i;
382
383  if (p1 == NULL) return NULL;
384  k = pGetMaxPower(p1,p2);
385  if (k == 0)
386    return pHead(p1);
387  else
388  {
389    number n;
390    poly p = pInit();
391
392    p->next = NULL;
393    for (i=1; i<=pVariables; i++)
394    {
395       pSetExp(p,i, pGetExp(p1,i)-k* pGetExp(p2,i));
396    }
397    nPower(p2->coef,k,&n);
398    pSetCoeff0(p,nDiv(p1->coef,n));
399    nDelete(&n);
400    pSetm(p);
401    return p;
402  }
403}
404// Returns as i-th entry of P the coefficient of the (i-1) power of
405// the leading monomial of p2 in p1
406//
407static void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet)
408{
409  int   maxPow;
410  poly  p,qp,Coeff;
411
412  if (*P == NULL)
413  {
414    *P = (polyset) omAlloc(5*sizeof(poly));
415    *SizeOfSet = 5;
416  }
417  p = pCopy(p1);
418  while (p != NULL)
419  {
420    qp = p->next;
421    p->next = NULL;
422    maxPow = pGetMaxPower(p,p2);
423    Coeff = pDivByMonom(p,p2);
424    if (maxPow > *SizeOfSet)
425    {
426      pEnlargeSet(P,*SizeOfSet,maxPow+1-*SizeOfSet);
427      *SizeOfSet = maxPow+1;
428    }
429    (*P)[maxPow] = pAdd((*P)[maxPow],Coeff);
430    pDelete(&p);
431    p = qp;
432  }
433}
434
435/*----------utilities for syzygies--------------*/
436poly pTakeOutComp(poly * p, int k)
437{
438  poly q = *p,qq=NULL,result = NULL;
439
440  if (q==NULL) return NULL;
441  BOOLEAN use_setmcomp=rOrd_SetCompRequiresSetm(currRing);
442  if (pGetComp(q)==k)
443  {
444    result = q;
445    do
446    {
447      pSetComp(q,0);
448      if (use_setmcomp) pSetmComp(q);
449      qq = q;
450      pIter(q);
451    }
452    while ((q!=NULL) && (pGetComp(q)==k));
453    *p = q;
454    pNext(qq) = NULL;
455  }
456  if (q==NULL) return result;
457  if (pGetComp(q) > k)
458  {
459    pSubComp(q,1);
460    if (use_setmcomp) pSetmComp(q);
461  }
462  poly pNext_q;
463  while ((pNext_q=pNext(q))!=NULL)
464  {
465    if (pGetComp(pNext_q)==k)
466    {
467      if (result==NULL)
468      {
469        result = pNext_q;
470        qq = result;
471      }
472      else
473      {
474        pNext(qq) = pNext_q;
475        pIter(qq);
476      }
477      pNext(q) = pNext(pNext_q);
478      pNext(qq) =NULL;
479      pSetComp(qq,0);
480      if (use_setmcomp) pSetmComp(qq);
481    }
482    else
483    {
484      /*pIter(q);*/ q=pNext_q;
485      if (pGetComp(q) > k)
486      {
487        pSubComp(q,1);
488        if (use_setmcomp) pSetmComp(q);
489      }
490    }
491  }
492  return result;
493}
494
495// Splits *p into two polys: *q which consists of all monoms with
496// component == comp and *p of all other monoms *lq == pLength(*q)
497void pTakeOutComp(poly *r_p, long comp, poly *r_q, int *lq)
498{
499  spolyrec pp, qq;
500  poly p, q, p_prev;
501  int l = 0;
502
503#ifdef HAVE_ASSUME
504  int lp = pLength(*r_p);
505#endif
506
507  pNext(&pp) = *r_p;
508  p = *r_p;
509  p_prev = &pp;
510  q = &qq;
511
512  while(p != NULL)
513  {
514    while (pGetComp(p) == comp)
515    {
516      pNext(q) = p;
517      pIter(q);
518      pSetComp(p, 0);
519      pSetmComp(p);
520      pIter(p);
521      l++;
522      if (p == NULL)
523      {
524        pNext(p_prev) = NULL;
525        goto Finish;
526      }
527    }
528    pNext(p_prev) = p;
529    p_prev = p;
530    pIter(p);
531  }
532
533  Finish:
534  pNext(q) = NULL;
535  *r_p = pNext(&pp);
536  *r_q = pNext(&qq);
537  *lq = l;
538#ifdef HAVE_ASSUME
539  assume(pLength(*r_p) + pLength(*r_q) == lp);
540#endif
541  pTest(*r_p);
542  pTest(*r_q);
543}
544
545void pDecrOrdTakeOutComp(poly *r_p, long comp, long order,
546                         poly *r_q, int *lq)
547{
548  spolyrec pp, qq;
549  poly p, q, p_prev;
550  int l = 0;
551
552  pNext(&pp) = *r_p;
553  p = *r_p;
554  p_prev = &pp;
555  q = &qq;
556
557#ifdef HAVE_ASSUME
558  if (p != NULL)
559  {
560    while (pNext(p) != NULL)
561    {
562      assume(pGetOrder(p) >= pGetOrder(pNext(p)));
563      pIter(p);
564    }
565  }
566  p = *r_p;
567#endif
568
569  while (p != NULL && pGetOrder(p) > order) pIter(p);
570
571  while(p != NULL && pGetOrder(p) == order)
572  {
573    while (pGetComp(p) == comp)
574    {
575      pNext(q) = p;
576      pIter(q);
577      pIter(p);
578      pSetComp(p, 0);
579      pSetmComp(p);
580      l++;
581      if (p == NULL || pGetOrder(p) != order)
582      {
583        pNext(p_prev) = p;
584        goto Finish;
585      }
586    }
587    pNext(p_prev) = p;
588    p_prev = p;
589    pIter(p);
590  }
591
592  Finish:
593  pNext(q) = NULL;
594  *r_p = pNext(&pp);
595  *r_q = pNext(&qq);
596  *lq = l;
597}
598
599#if 1
600poly pTakeOutComp1(poly * p, int k)
601{
602  poly q = *p;
603
604  if (q==NULL) return NULL;
605
606  poly qq=NULL,result = NULL;
607
608  if (pGetComp(q)==k)
609  {
610    result = q; /* *p */
611    while ((q!=NULL) && (pGetComp(q)==k))
612    {
613      pSetComp(q,0);
614      pSetmComp(q);
615      qq = q;
616      pIter(q);
617    }
618    *p = q;
619    pNext(qq) = NULL;
620  }
621  if (q==NULL) return result;
622//  if (pGetComp(q) > k) pGetComp(q)--;
623  while (pNext(q)!=NULL)
624  {
625    if (pGetComp(pNext(q))==k)
626    {
627      if (result==NULL)
628      {
629        result = pNext(q);
630        qq = result;
631      }
632      else
633      {
634        pNext(qq) = pNext(q);
635        pIter(qq);
636      }
637      pNext(q) = pNext(pNext(q));
638      pNext(qq) =NULL;
639      pSetComp(qq,0);
640      pSetmComp(qq);
641    }
642    else
643    {
644      pIter(q);
645//      if (pGetComp(q) > k) pGetComp(q)--;
646    }
647  }
648  return result;
649}
650#endif
651
652void pDeleteComp(poly * p,int k)
653{
654  poly q;
655
656  while ((*p!=NULL) && (pGetComp(*p)==k)) pLmDelete(p);
657  if (*p==NULL) return;
658  q = *p;
659  if (pGetComp(q)>k)
660  {
661    pSubComp(q,1);
662    pSetmComp(q);
663  }
664  while (pNext(q)!=NULL)
665  {
666    if (pGetComp(pNext(q))==k)
667      pLmDelete(&(pNext(q)));
668    else
669    {
670      pIter(q);
671      if (pGetComp(q)>k)
672      {
673        pSubComp(q,1);
674        pSetmComp(q);
675      }
676    }
677  }
678}
679/*----------end of utilities for syzygies--------------*/
680
681/*2
682* pair has no common factor ? or is no polynomial
683*/
684BOOLEAN pHasNotCF(poly p1, poly p2)
685{
686
687  if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
688    return FALSE;
689  int i = pVariables;
690  loop
691  {
692    if ((pGetExp(p1, i) > 0) && (pGetExp(p2, i) > 0))   return FALSE;
693    i--;
694    if (i == 0)                                         return TRUE;
695  }
696}
697
698/*2
699*divides p1 by its leading coefficient
700*/
701void pNorm(poly p1)
702{
703#ifdef HAVE_RINGS
704  if (rField_is_Ring(currRing))
705  {
706    Werror("pNorm not possible in the case of coefficient rings.");
707  }
708  else
709#endif
710  if (p1!=NULL)
711  {
712    if (pNext(p1)==NULL)
713    {
714      pSetCoeff(p1,nInit(1));
715      return;
716    }
717    poly h;
718    if (!nIsOne(pGetCoeff(p1)))
719    {
720      number k, c;
721      nNormalize(pGetCoeff(p1));
722      k = pGetCoeff(p1);
723      c = nInit(1);
724      pSetCoeff0(p1,c);
725      h = pNext(p1);
726      while (h!=NULL)
727      {
728        c=nDiv(pGetCoeff(h),k);
729        // no need to normalize: Z/p, R
730        // normalize already in nDiv: Q_a, Z/p_a
731        // remains: Q
732        if (rField_is_Q() && (!nIsOne(c))) nNormalize(c);
733        pSetCoeff(h,c);
734        pIter(h);
735      }
736      nDelete(&k);
737    }
738    else
739    {
740      if (nNormalize != nDummy2)
741      {
742        h = pNext(p1);
743        while (h!=NULL)
744        {
745          nNormalize(pGetCoeff(h));
746          pIter(h);
747        }
748      }
749    }
750  }
751}
752
753/*2
754*normalize all coefficients
755*/
756void p_Normalize(poly p,const ring r)
757{
758  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
759  while (p!=NULL)
760  {
761#ifdef LDEBUG
762    if (currRing==r) {nTest(pGetCoeff(p));}
763#endif
764    n_Normalize(pGetCoeff(p),r);
765    pIter(p);
766  }
767}
768
769// splits p into polys with Exp(n) == 0 and Exp(n) != 0
770// Poly with Exp(n) != 0 is reversed
771static void pSplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero)
772{
773  if (p == NULL)
774  {
775    *non_zero = NULL;
776    *zero = NULL;
777    return;
778  }
779  spolyrec sz;
780  poly z, n_z, next;
781  z = &sz;
782  n_z = NULL;
783
784  while(p != NULL)
785  {
786    next = pNext(p);
787    if (pGetExp(p, n) == 0)
788    {
789      pNext(z) = p;
790      pIter(z);
791    }
792    else
793    {
794      pNext(p) = n_z;
795      n_z = p;
796    }
797    p = next;
798  }
799  pNext(z) = NULL;
800  *zero = pNext(&sz);
801  *non_zero = n_z;
802  return;
803}
804
805/*3
806* substitute the n-th variable by 1 in p
807* destroy p
808*/
809static poly pSubst1 (poly p,int n)
810{
811  poly qq=NULL, result = NULL;
812  poly zero=NULL, non_zero=NULL;
813
814  // reverse, so that add is likely to be linear
815  pSplitAndReversePoly(p, n, &non_zero, &zero);
816
817  while (non_zero != NULL)
818  {
819    assume(pGetExp(non_zero, n) != 0);
820    qq = non_zero;
821    pIter(non_zero);
822    qq->next = NULL;
823    pSetExp(qq,n,0);
824    pSetm(qq);
825    result = pAdd(result,qq);
826  }
827  p = pAdd(result, zero);
828  pTest(p);
829  return p;
830}
831
832/*3
833* substitute the n-th variable by number e in p
834* destroy p
835*/
836static poly pSubst2 (poly p,int n, number e)
837{
838  assume( ! nIsZero(e) );
839  poly qq,result = NULL;
840  number nn, nm;
841  poly zero, non_zero;
842
843  // reverse, so that add is likely to be linear
844  pSplitAndReversePoly(p, n, &non_zero, &zero);
845
846  while (non_zero != NULL)
847  {
848    assume(pGetExp(non_zero, n) != 0);
849    qq = non_zero;
850    pIter(non_zero);
851    qq->next = NULL;
852    nPower(e, pGetExp(qq, n), &nn);
853    nm = nMult(nn, pGetCoeff(qq));
854#ifdef HAVE_RINGS
855    if (nIsZero(nm))
856    {
857      pLmFree(&qq);
858      nDelete(&nm);
859    }
860    else
861#endif
862    {
863      pSetCoeff(qq, nm);
864      pSetExp(qq, n, 0);
865      pSetm(qq);
866      result = pAdd(result,qq);
867    }
868    nDelete(&nn);
869  }
870  p = pAdd(result, zero);
871  pTest(p);
872  return p;
873}
874
875
876/* delete monoms whose n-th exponent is different from zero */
877poly pSubst0(poly p, int n)
878{
879  spolyrec res;
880  poly h = &res;
881  pNext(h) = p;
882
883  while (pNext(h)!=NULL)
884  {
885    if (pGetExp(pNext(h),n)!=0)
886    {
887      pLmDelete(&pNext(h));
888    }
889    else
890    {
891      pIter(h);
892    }
893  }
894  pTest(pNext(&res));
895  return pNext(&res);
896}
897
898/*2
899* substitute the n-th variable by e in p
900* destroy p
901*/
902poly pSubst(poly p, int n, poly e)
903{
904  if (e == NULL) return pSubst0(p, n);
905
906  if (pIsConstant(e))
907  {
908    if (nIsOne(pGetCoeff(e))) return pSubst1(p,n);
909    else return pSubst2(p, n, pGetCoeff(e));
910  }
911
912#ifdef HAVE_PLURAL
913  if (rIsPluralRing(currRing))
914  {
915    return nc_pSubst(p,n,e);
916  }
917#endif
918
919  int exponent,i;
920  poly h, res, m;
921  int *me,*ee;
922  number nu,nu1;
923
924  me=(int *)omAlloc((pVariables+1)*sizeof(int));
925  ee=(int *)omAlloc((pVariables+1)*sizeof(int));
926  if (e!=NULL) pGetExpV(e,ee);
927  res=NULL;
928  h=p;
929  while (h!=NULL)
930  {
931    if ((e!=NULL) || (pGetExp(h,n)==0))
932    {
933      m=pHead(h);
934      pGetExpV(m,me);
935      exponent=me[n];
936      me[n]=0;
937      for(i=pVariables;i>0;i--)
938        me[i]+=exponent*ee[i];
939      pSetExpV(m,me);
940      if (e!=NULL)
941      {
942        nPower(pGetCoeff(e),exponent,&nu);
943        nu1=nMult(pGetCoeff(m),nu);
944        nDelete(&nu);
945        pSetCoeff(m,nu1);
946      }
947      res=pAdd(res,m);
948    }
949    pLmDelete(&h);
950  }
951  omFreeSize((ADDRESS)me,(pVariables+1)*sizeof(int));
952  omFreeSize((ADDRESS)ee,(pVariables+1)*sizeof(int));
953  return res;
954}
955
956/* Returns TRUE if
957     * LM(p) | LM(lcm)
958     * LC(p) | LC(lcm) only if ring
959     * Exists i, j:
960         * LE(p, i)  != LE(lcm, i)
961         * LE(p1, i) != LE(lcm, i)   ==> LCM(p1, p) != lcm
962         * LE(p, j)  != LE(lcm, j)
963         * LE(p2, j) != LE(lcm, j)   ==> LCM(p2, p) != lcm
964*/
965BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
966{
967  int k, j;
968
969  if (lcm==NULL) return FALSE;
970
971  for (j=pVariables; j; j--)
972    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
973  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
974  for (j=pVariables; j; j--)
975  {
976    if (pGetExp(p1,j)!=pGetExp(lcm,j))
977    {
978      if (pGetExp(p,j)!=pGetExp(lcm,j))
979      {
980        for (k=pVariables; k>j; k--)
981        {
982          if ((pGetExp(p,k)!=pGetExp(lcm,k))
983          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
984            return TRUE;
985        }
986        for (k=j-1; k; k--)
987        {
988          if ((pGetExp(p,k)!=pGetExp(lcm,k))
989          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
990            return TRUE;
991        }
992        return FALSE;
993      }
994    }
995    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
996    {
997      if (pGetExp(p,j)!=pGetExp(lcm,j))
998      {
999        for (k=pVariables; k>j; k--)
1000        {
1001          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1002          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1003            return TRUE;
1004        }
1005        for (k=j-1; k!=0 ; k--)
1006        {
1007          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1008          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1009            return TRUE;
1010        }
1011        return FALSE;
1012      }
1013    }
1014  }
1015  return FALSE;
1016}
1017#ifdef HAVE_RATGRING
1018BOOLEAN pCompareChainPart (poly p,poly p1,poly p2,poly lcm)
1019{
1020  int k, j;
1021
1022  if (lcm==NULL) return FALSE;
1023
1024  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1025    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
1026  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1027  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1028  {
1029    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1030    {
1031      if (pGetExp(p,j)!=pGetExp(lcm,j))
1032      {
1033        for (k=pVariables; k>j; k--)
1034        for (k=currRing->real_var_end; k>j; k--)
1035        {
1036          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1037          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1038            return TRUE;
1039        }
1040        for (k=j-1; k>=currRing->real_var_start; k--)
1041        {
1042          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1043          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1044            return TRUE;
1045        }
1046        return FALSE;
1047      }
1048    }
1049    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1050    {
1051      if (pGetExp(p,j)!=pGetExp(lcm,j))
1052      {
1053        for (k=currRing->real_var_end; k>j; k--)
1054        {
1055          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1056          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1057            return TRUE;
1058        }
1059        for (k=j-1; k>=currRing->real_var_start; k--)
1060        {
1061          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1062          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1063            return TRUE;
1064        }
1065        return FALSE;
1066      }
1067    }
1068  }
1069  return FALSE;
1070}
1071#endif
1072
1073int pSize(poly p)
1074{
1075  int count = 0;
1076  while ( p != NULL )
1077  {
1078    count+= nSize( pGetCoeff( p ) );
1079    pIter( p );
1080  }
1081  return count;
1082}
1083
1084/*2
1085* returns the length of a (numbers of monomials)
1086* respect syzComp
1087*/
1088poly pLast(poly a, int &l)
1089{
1090  if (a == NULL)
1091  {
1092    l = 0;
1093    return NULL;
1094  }
1095  l = 1;
1096  if (! rIsSyzIndexRing(currRing))
1097  {
1098    while (pNext(a)!=NULL)
1099    {
1100      pIter(a);
1101      l++;
1102    }
1103  }
1104  else
1105  {
1106    int curr_limit = rGetCurrSyzLimit(currRing);
1107    poly pp = a;
1108    while ((a=pNext(a))!=NULL)
1109    {
1110      if (pGetComp(a)<=curr_limit/*syzComp*/)
1111        l++;
1112      else break;
1113      pp = a;
1114    }
1115    a=pp;
1116  }
1117  return a;
1118}
1119
Note: See TracBrowser for help on using the repository browser.