source: git/kernel/polys.cc @ 47b2b2d

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