source: git/kernel/polys.cc @ d9ff1b

spielwiese
Last change on this file since d9ff1b was 48884f2, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: pLast -> polys.cc git-svn-id: file:///usr/local/Singular/svn/trunk@12036 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.39 2009-07-29 10:07:22 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    if (currRing==r) {nTest(pGetCoeff(p));}
797    n_Normalize(pGetCoeff(p),r);
798    pIter(p);
799  }
800}
801
802// splits p into polys with Exp(n) == 0 and Exp(n) != 0
803// Poly with Exp(n) != 0 is reversed
804static void pSplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero)
805{
806  if (p == NULL)
807  {
808    *non_zero = NULL;
809    *zero = NULL;
810    return;
811  }
812  spolyrec sz;
813  poly z, n_z, next;
814  z = &sz;
815  n_z = NULL;
816
817  while(p != NULL)
818  {
819    next = pNext(p);
820    if (pGetExp(p, n) == 0)
821    {
822      pNext(z) = p;
823      pIter(z);
824    }
825    else
826    {
827      pNext(p) = n_z;
828      n_z = p;
829    }
830    p = next;
831  }
832  pNext(z) = NULL;
833  *zero = pNext(&sz);
834  *non_zero = n_z;
835  return;
836}
837
838/*3
839* substitute the n-th variable by 1 in p
840* destroy p
841*/
842static poly pSubst1 (poly p,int n)
843{
844  poly qq=NULL, result = NULL;
845  poly zero=NULL, non_zero=NULL;
846
847  // reverse, so that add is likely to be linear
848  pSplitAndReversePoly(p, n, &non_zero, &zero);
849
850  while (non_zero != NULL)
851  {
852    assume(pGetExp(non_zero, n) != 0);
853    qq = non_zero;
854    pIter(non_zero);
855    qq->next = NULL;
856    pSetExp(qq,n,0);
857    pSetm(qq);
858    result = pAdd(result,qq);
859  }
860  p = pAdd(result, zero);
861  pTest(p);
862  return p;
863}
864
865/*3
866* substitute the n-th variable by number e in p
867* destroy p
868*/
869static poly pSubst2 (poly p,int n, number e)
870{
871  assume( ! nIsZero(e) );
872  poly qq,result = NULL;
873  number nn, nm;
874  poly zero, non_zero;
875
876  // reverse, so that add is likely to be linear
877  pSplitAndReversePoly(p, n, &non_zero, &zero);
878
879  while (non_zero != NULL)
880  {
881    assume(pGetExp(non_zero, n) != 0);
882    qq = non_zero;
883    pIter(non_zero);
884    qq->next = NULL;
885    nPower(e, pGetExp(qq, n), &nn);
886    nm = nMult(nn, pGetCoeff(qq));
887#ifdef HAVE_RINGS
888    if (nIsZero(nm))
889    {
890      pLmFree(&qq);
891      nDelete(&nm);
892    }
893    else
894#endif
895    {
896      pSetCoeff(qq, nm);
897      pSetExp(qq, n, 0);
898      pSetm(qq);
899      result = pAdd(result,qq);
900    }
901    nDelete(&nn);
902  }
903  p = pAdd(result, zero);
904  pTest(p);
905  return p;
906}
907
908
909/* delete monoms whose n-th exponent is different from zero */
910poly pSubst0(poly p, int n)
911{
912  spolyrec res;
913  poly h = &res;
914  pNext(h) = p;
915
916  while (pNext(h)!=NULL)
917  {
918    if (pGetExp(pNext(h),n)!=0)
919    {
920      pDeleteLm(&pNext(h));
921    }
922    else
923    {
924      pIter(h);
925    }
926  }
927  pTest(pNext(&res));
928  return pNext(&res);
929}
930
931/*2
932* substitute the n-th variable by e in p
933* destroy p
934*/
935poly pSubst(poly p, int n, poly e)
936{
937  if (e == NULL) return pSubst0(p, n);
938
939  if (pIsConstant(e))
940  {
941    if (nIsOne(pGetCoeff(e))) return pSubst1(p,n);
942    else return pSubst2(p, n, pGetCoeff(e));
943  }
944
945#ifdef HAVE_PLURAL
946  if (rIsPluralRing(currRing))
947  {
948    return nc_pSubst(p,n,e);
949  }
950#endif
951
952  int exponent,i;
953  poly h, res, m;
954  int *me,*ee;
955  number nu,nu1;
956
957  me=(int *)omAlloc((pVariables+1)*sizeof(int));
958  ee=(int *)omAlloc((pVariables+1)*sizeof(int));
959  if (e!=NULL) pGetExpV(e,ee);
960  res=NULL;
961  h=p;
962  while (h!=NULL)
963  {
964    if ((e!=NULL) || (pGetExp(h,n)==0))
965    {
966      m=pHead(h);
967      pGetExpV(m,me);
968      exponent=me[n];
969      me[n]=0;
970      for(i=pVariables;i>0;i--)
971        me[i]+=exponent*ee[i];
972      pSetExpV(m,me);
973      if (e!=NULL)
974      {
975        nPower(pGetCoeff(e),exponent,&nu);
976        nu1=nMult(pGetCoeff(m),nu);
977        nDelete(&nu);
978        pSetCoeff(m,nu1);
979      }
980      res=pAdd(res,m);
981    }
982    pDeleteLm(&h);
983  }
984  omFreeSize((ADDRESS)me,(pVariables+1)*sizeof(int));
985  omFreeSize((ADDRESS)ee,(pVariables+1)*sizeof(int));
986  return res;
987}
988
989/* Returns TRUE if
990     * LM(p) | LM(lcm)
991     * LC(p) | LC(lcm) only if ring
992     * Exists i, j:
993         * LE(p, i)  != LE(lcm, i)
994         * LE(p1, i) != LE(lcm, i)   ==> LCM(p1, p) != lcm
995         * LE(p, j)  != LE(lcm, j)
996         * LE(p2, j) != LE(lcm, j)   ==> LCM(p2, p) != lcm
997*/
998BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
999{
1000  int k, j;
1001
1002  if (lcm==NULL) return FALSE;
1003
1004  for (j=pVariables; j; j--)
1005    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
1006  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1007  for (j=pVariables; j; j--)
1008  {
1009    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1010    {
1011      if (pGetExp(p,j)!=pGetExp(lcm,j))
1012      {
1013        for (k=pVariables; k>j; k--)
1014        {
1015          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1016          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1017            return TRUE;
1018        }
1019        for (k=j-1; k; k--)
1020        {
1021          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1022          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1023            return TRUE;
1024        }
1025        return FALSE;
1026      }
1027    }
1028    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1029    {
1030      if (pGetExp(p,j)!=pGetExp(lcm,j))
1031      {
1032        for (k=pVariables; k>j; k--)
1033        {
1034          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1035          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1036            return TRUE;
1037        }
1038        for (k=j-1; k!=0 ; k--)
1039        {
1040          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1041          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1042            return TRUE;
1043        }
1044        return FALSE;
1045      }
1046    }
1047  }
1048  return FALSE;
1049}
1050#ifdef HAVE_RATGRING
1051BOOLEAN pCompareChainPart (poly p,poly p1,poly p2,poly lcm)
1052{
1053  int k, j;
1054
1055  if (lcm==NULL) return FALSE;
1056
1057  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1058    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
1059  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1060  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1061  {
1062    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1063    {
1064      if (pGetExp(p,j)!=pGetExp(lcm,j))
1065      {
1066        for (k=pVariables; k>j; k--)
1067        for (k=currRing->real_var_end; k>j; k--)
1068        {
1069          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1070          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1071            return TRUE;
1072        }
1073        for (k=j-1; k>=currRing->real_var_start; k--)
1074        {
1075          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1076          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1077            return TRUE;
1078        }
1079        return FALSE;
1080      }
1081    }
1082    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1083    {
1084      if (pGetExp(p,j)!=pGetExp(lcm,j))
1085      {
1086        for (k=currRing->real_var_end; k>j; k--)
1087        {
1088          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1089          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1090            return TRUE;
1091        }
1092        for (k=j-1; k>=currRing->real_var_start; k--)
1093        {
1094          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1095          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1096            return TRUE;
1097        }
1098        return FALSE;
1099      }
1100    }
1101  }
1102  return FALSE;
1103}
1104#endif
1105
1106int pSize(poly p)
1107{
1108  int count = 0;
1109  while ( p != NULL )
1110  {
1111    count+= nSize( pGetCoeff( p ) );
1112    pIter( p );
1113  }
1114  return count;
1115}
1116
1117/*2
1118* returns the length of a (numbers of monomials)
1119* respect syzComp
1120*/
1121poly pLast(poly a, int &l)
1122{
1123  if (a == NULL)
1124  {
1125    l = 0;
1126    return NULL;
1127  }
1128  l = 1;
1129  if (! rIsSyzIndexRing(currRing))
1130  {
1131    while (pNext(a)!=NULL)
1132    {
1133      pIter(a);
1134      l++;
1135    }
1136  }
1137  else
1138  {
1139    int curr_limit = rGetCurrSyzLimit(currRing);
1140    poly pp = a;
1141    while ((a=pNext(a))!=NULL)
1142    {
1143      if (pGetComp(a)<=curr_limit/*syzComp*/)
1144        l++;
1145      else break;
1146      pp = a;
1147    }
1148    a=pp;
1149  }
1150  return a;
1151}
1152
Note: See TracBrowser for help on using the repository browser.