source: git/kernel/polys.cc @ 075c01

spielwiese
Last change on this file since 075c01 was 997e23, checked in by Hans Schoenemann <hannes@…>, 14 years ago
debug stuff git-svn-id: file:///usr/local/Singular/svn/trunk@12834 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$ */
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 "options.h"
16#include "omalloc.h"
17#include "febase.h"
18#include "numbers.h"
19#include "polys.h"
20#include "ring.h"
21#include "sbuckets.h"
22
23#ifdef HAVE_PLURAL
24#include "gring.h"
25#include "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        pDeleteLm(&result);
167        a=result;
168      }
169      else
170      {
171        pDeleteLm(&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_DeleteLm(&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_DeleteLm(&rc,r);
251      s--;
252      return s;
253    }
254  }
255done:
256  if (r->cf->nIsZero(pGetCoeff(rc))) p_DeleteLm(&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_DeleteLm(&rc, r);
273          goto finish;
274        }
275    }
276#endif
277   
278    p_Setm(rc,r);
279  }
280finish: 
281  return s;
282}
283
284BOOLEAN _p_Test(poly p, ring r, int level);
285poly pmInit(const char *st, BOOLEAN &ok)
286{
287  poly p;
288  const char *s=p_Read(st,p,currRing);
289  if (*s!='\0')
290  {
291    if ((s!=st)&&isdigit(st[0]))
292    {
293      errorreported=TRUE;
294    }
295    ok=FALSE;
296    pDelete(&p);
297    return NULL;
298  }
299  #ifdef PDEBUG
300  _p_Test(p,currRing,PDEBUG);
301  #endif
302  ok=!errorreported;
303  return p;
304}
305
306/*2
307*make p homogeneous by multiplying the monomials by powers of x_varnum
308*assume: deg(var(varnum))==1
309*/
310poly pHomogen (poly p, int varnum)
311{
312  pFDegProc deg;
313  if (pLexOrder && (currRing->order[0]==ringorder_lp))
314    deg=pTotaldegree;
315  else
316    deg=pFDeg;
317
318  poly q=NULL, qn;
319  int  o,ii;
320  sBucket_pt bp;
321
322  if (p!=NULL)
323  {
324    if ((varnum < 1) || (varnum > pVariables))
325    {
326      return NULL;
327    }
328    o=deg(p,currRing);
329    q=pNext(p);
330    while (q != NULL)
331    {
332      ii=deg(q,currRing);
333      if (ii>o) o=ii;
334      pIter(q);
335    }
336    q = pCopy(p);
337    bp = sBucketCreate(currRing);
338    while (q != NULL)
339    {
340      ii = o-deg(q,currRing);
341      if (ii!=0)
342      {
343        pAddExp(q,varnum, (long)ii);
344        pSetm(q);
345      }
346      qn = pNext(q);
347      pNext(q) = NULL;
348      sBucket_Add_p(bp, q, 1);
349      q = qn;
350    }
351    sBucketDestroyAdd(bp, &q, &ii);
352  }
353  return q;
354}
355
356/*2
357*replaces the maximal powers of the leading monomial of p2 in p1 by
358*the same powers of n, utility for dehomogenization
359*/
360poly pDehomogen (poly p1,poly p2,number n)
361{
362  polyset P;
363  int     SizeOfSet=5;
364  int     i;
365  poly    p;
366  number  nn;
367
368  P = (polyset)omAlloc0(5*sizeof(poly));
369  //for (i=0; i<5; i++)
370  //{
371  //  P[i] = NULL;
372  //}
373  pCancelPolyByMonom(p1,p2,&P,&SizeOfSet);
374  p = P[0];
375  //P[0] = NULL ;// for safety, may be removed later
376  for (i=1; i<SizeOfSet; i++)
377  {
378    if (P[i] != NULL)
379    {
380      nPower(n,i,&nn);
381      pMult_nn(P[i],nn);
382      p = pAdd(p,P[i]);
383      //P[i] =NULL; // for safety, may be removed later
384      nDelete(&nn);
385    }
386  }
387  omFreeSize((ADDRESS)P,SizeOfSet*sizeof(poly));
388  return p;
389}
390
391/*4
392*Returns the exponent of the maximal power of the leading monomial of
393*p2 in that of p1
394*/
395static int pGetMaxPower (poly p1,poly p2)
396{
397  int     i,k,res = 32000; /*a very large integer*/
398
399  if (p1 == NULL) return 0;
400  for (i=1; i<=pVariables; i++)
401  {
402    if ( pGetExp(p2,i) != 0)
403    {
404      k =  pGetExp(p1,i) /  pGetExp(p2,i);
405      if (k < res) res = k;
406    }
407  }
408  return res;
409}
410
411/*2
412*Returns as i-th entry of P the coefficient of the (i-1) power of
413*the leading monomial of p2 in p1
414*/
415void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet)
416{
417  int   maxPow;
418  poly  p,qp,Coeff;
419
420  if (*P == NULL)
421  {
422    *P = (polyset) omAlloc(5*sizeof(poly));
423    *SizeOfSet = 5;
424  }
425  p = pCopy(p1);
426  while (p != NULL)
427  {
428    qp = p->next;
429    p->next = NULL;
430    maxPow = pGetMaxPower(p,p2);
431    Coeff = pDivByMonom(p,p2);
432    if (maxPow > *SizeOfSet)
433    {
434      pEnlargeSet(P,*SizeOfSet,maxPow+1-*SizeOfSet);
435      *SizeOfSet = maxPow+1;
436    }
437    (*P)[maxPow] = pAdd((*P)[maxPow],Coeff);
438    pDelete(&p);
439    p = qp;
440  }
441}
442
443/*2
444*returns the leading monomial of p1 divided by the maximal power of that
445*of p2
446*/
447poly pDivByMonom (poly p1,poly p2)
448{
449  int     k, i;
450
451  if (p1 == NULL) return NULL;
452  k = pGetMaxPower(p1,p2);
453  if (k == 0)
454    return pHead(p1);
455  else
456  {
457    number n;
458    poly p = pInit();
459
460    p->next = NULL;
461    for (i=1; i<=pVariables; i++)
462    {
463       pSetExp(p,i, pGetExp(p1,i)-k* pGetExp(p2,i));
464    }
465    nPower(p2->coef,k,&n);
466    pSetCoeff0(p,nDiv(p1->coef,n));
467    nDelete(&n);
468    pSetm(p);
469    return p;
470  }
471}
472/*----------utilities for syzygies--------------*/
473poly pTakeOutComp(poly * p, int k)
474{
475  poly q = *p,qq=NULL,result = NULL;
476
477  if (q==NULL) return NULL;
478  BOOLEAN use_setmcomp=rOrd_SetCompRequiresSetm(currRing);
479  if (pGetComp(q)==k)
480  {
481    result = q;
482    do
483    {
484      pSetComp(q,0);
485      if (use_setmcomp) pSetmComp(q);
486      qq = q;
487      pIter(q);
488    }
489    while ((q!=NULL) && (pGetComp(q)==k));
490    *p = q;
491    pNext(qq) = NULL;
492  }
493  if (q==NULL) return result;
494  if (pGetComp(q) > k)
495  {
496    pDecrComp(q);
497    if (use_setmcomp) pSetmComp(q);
498  }
499  poly pNext_q;
500  while ((pNext_q=pNext(q))!=NULL)
501  {
502    if (pGetComp(pNext_q)==k)
503    {
504      if (result==NULL)
505      {
506        result = pNext_q;
507        qq = result;
508      }
509      else
510      {
511        pNext(qq) = pNext_q;
512        pIter(qq);
513      }
514      pNext(q) = pNext(pNext_q);
515      pNext(qq) =NULL;
516      pSetComp(qq,0);
517      if (use_setmcomp) pSetmComp(qq);
518    }
519    else
520    {
521      /*pIter(q);*/ q=pNext_q;
522      if (pGetComp(q) > k)
523      {
524        pDecrComp(q);
525        if (use_setmcomp) pSetmComp(q);
526      }
527    }
528  }
529  return result;
530}
531
532// Splits *p into two polys: *q which consists of all monoms with
533// component == comp and *p of all other monoms *lq == pLength(*q)
534void pTakeOutComp(poly *r_p, long comp, poly *r_q, int *lq)
535{
536  spolyrec pp, qq;
537  poly p, q, p_prev;
538  int l = 0;
539
540#ifdef HAVE_ASSUME
541  int lp = pLength(*r_p);
542#endif
543
544  pNext(&pp) = *r_p;
545  p = *r_p;
546  p_prev = &pp;
547  q = &qq;
548
549  while(p != NULL)
550  {
551    while (pGetComp(p) == comp)
552    {
553      pNext(q) = p;
554      pIter(q);
555      pSetComp(p, 0);
556      pSetmComp(p);
557      pIter(p);
558      l++;
559      if (p == NULL)
560      {
561        pNext(p_prev) = NULL;
562        goto Finish;
563      }
564    }
565    pNext(p_prev) = p;
566    p_prev = p;
567    pIter(p);
568  }
569
570  Finish:
571  pNext(q) = NULL;
572  *r_p = pNext(&pp);
573  *r_q = pNext(&qq);
574  *lq = l;
575#ifdef HAVE_ASSUME
576  assume(pLength(*r_p) + pLength(*r_q) == lp);
577#endif
578  pTest(*r_p);
579  pTest(*r_q);
580}
581
582void pDecrOrdTakeOutComp(poly *r_p, long comp, long order,
583                         poly *r_q, int *lq)
584{
585  spolyrec pp, qq;
586  poly p, q, p_prev;
587  int l = 0;
588
589  pNext(&pp) = *r_p;
590  p = *r_p;
591  p_prev = &pp;
592  q = &qq;
593
594#ifdef HAVE_ASSUME
595  if (p != NULL)
596  {
597    while (pNext(p) != NULL)
598    {
599      assume(pGetOrder(p) >= pGetOrder(pNext(p)));
600      pIter(p);
601    }
602  }
603  p = *r_p;
604#endif
605
606  while (p != NULL && pGetOrder(p) > order) pIter(p);
607
608  while(p != NULL && pGetOrder(p) == order)
609  {
610    while (pGetComp(p) == comp)
611    {
612      pNext(q) = p;
613      pIter(q);
614      pIter(p);
615      pSetComp(p, 0);
616      pSetmComp(p);
617      l++;
618      if (p == NULL || pGetOrder(p) != order)
619      {
620        pNext(p_prev) = p;
621        goto Finish;
622      }
623    }
624    pNext(p_prev) = p;
625    p_prev = p;
626    pIter(p);
627  }
628
629  Finish:
630  pNext(q) = NULL;
631  *r_p = pNext(&pp);
632  *r_q = pNext(&qq);
633  *lq = l;
634}
635
636#if 1
637poly pTakeOutComp1(poly * p, int k)
638{
639  poly q = *p;
640
641  if (q==NULL) return NULL;
642
643  poly qq=NULL,result = NULL;
644
645  if (pGetComp(q)==k)
646  {
647    result = q; /* *p */
648    while ((q!=NULL) && (pGetComp(q)==k))
649    {
650      pSetComp(q,0);
651      pSetmComp(q);
652      qq = q;
653      pIter(q);
654    }
655    *p = q;
656    pNext(qq) = NULL;
657  }
658  if (q==NULL) return result;
659//  if (pGetComp(q) > k) pGetComp(q)--;
660  while (pNext(q)!=NULL)
661  {
662    if (pGetComp(pNext(q))==k)
663    {
664      if (result==NULL)
665      {
666        result = pNext(q);
667        qq = result;
668      }
669      else
670      {
671        pNext(qq) = pNext(q);
672        pIter(qq);
673      }
674      pNext(q) = pNext(pNext(q));
675      pNext(qq) =NULL;
676      pSetComp(qq,0);
677      pSetmComp(qq);
678    }
679    else
680    {
681      pIter(q);
682//      if (pGetComp(q) > k) pGetComp(q)--;
683    }
684  }
685  return result;
686}
687#endif
688
689void pDeleteComp(poly * p,int k)
690{
691  poly q;
692
693  while ((*p!=NULL) && (pGetComp(*p)==k)) pDeleteLm(p);
694  if (*p==NULL) return;
695  q = *p;
696  if (pGetComp(q)>k)
697  {
698    pDecrComp(q);
699    pSetmComp(q);
700  }
701  while (pNext(q)!=NULL)
702  {
703    if (pGetComp(pNext(q))==k)
704      pDeleteLm(&(pNext(q)));
705    else
706    {
707      pIter(q);
708      if (pGetComp(q)>k)
709      {
710        pDecrComp(q);
711        pSetmComp(q);
712      }
713    }
714  }
715}
716/*----------end of utilities for syzygies--------------*/
717
718/*2
719* pair has no common factor ? or is no polynomial
720*/
721BOOLEAN pHasNotCF(poly p1, poly p2)
722{
723
724  if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
725    return FALSE;
726  int i = pVariables;
727  loop
728  {
729    if ((pGetExp(p1, i) > 0) && (pGetExp(p2, i) > 0))   return FALSE;
730    i--;
731    if (i == 0)                                         return TRUE;
732  }
733}
734
735/*2
736*divides p1 by its leading coefficient
737*/
738void pNorm(poly p1)
739{
740#ifdef HAVE_RINGS
741  if (rField_is_Ring(currRing))
742  {
743    Werror("pNorm not possible in the case of coefficient rings.");
744  }
745  else
746#endif
747  if (p1!=NULL)
748  {
749    if (pNext(p1)==NULL)
750    {
751      pSetCoeff(p1,nInit(1));
752      return;
753    }
754    poly h;
755    if (!nIsOne(pGetCoeff(p1)))
756    {
757      number k, c;
758      nNormalize(pGetCoeff(p1));
759      k = pGetCoeff(p1);
760      c = nInit(1);
761      pSetCoeff0(p1,c);
762      h = pNext(p1);
763      while (h!=NULL)
764      {
765        c=nDiv(pGetCoeff(h),k);
766        // no need to normalize: Z/p, R
767        // normalize already in nDiv: Q_a, Z/p_a
768        // remains: Q
769        if (rField_is_Q() && (!nIsOne(c))) nNormalize(c);
770        pSetCoeff(h,c);
771        pIter(h);
772      }
773      nDelete(&k);
774    }
775    else
776    {
777      if (nNormalize != nDummy2)
778      {
779        h = pNext(p1);
780        while (h!=NULL)
781        {
782          nNormalize(pGetCoeff(h));
783          pIter(h);
784        }
785      }
786    }
787  }
788}
789
790/*2
791*normalize all coefficients
792*/
793void p_Normalize(poly p,const ring r)
794{
795  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
796  while (p!=NULL)
797  {
798#ifdef LDEBUG
799    if (currRing==r) {nTest(pGetCoeff(p));}
800#endif
801    n_Normalize(pGetCoeff(p),r);
802    pIter(p);
803  }
804}
805
806// splits p into polys with Exp(n) == 0 and Exp(n) != 0
807// Poly with Exp(n) != 0 is reversed
808static void pSplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero)
809{
810  if (p == NULL)
811  {
812    *non_zero = NULL;
813    *zero = NULL;
814    return;
815  }
816  spolyrec sz;
817  poly z, n_z, next;
818  z = &sz;
819  n_z = NULL;
820
821  while(p != NULL)
822  {
823    next = pNext(p);
824    if (pGetExp(p, n) == 0)
825    {
826      pNext(z) = p;
827      pIter(z);
828    }
829    else
830    {
831      pNext(p) = n_z;
832      n_z = p;
833    }
834    p = next;
835  }
836  pNext(z) = NULL;
837  *zero = pNext(&sz);
838  *non_zero = n_z;
839  return;
840}
841
842/*3
843* substitute the n-th variable by 1 in p
844* destroy p
845*/
846static poly pSubst1 (poly p,int n)
847{
848  poly qq=NULL, result = NULL;
849  poly zero=NULL, non_zero=NULL;
850
851  // reverse, so that add is likely to be linear
852  pSplitAndReversePoly(p, n, &non_zero, &zero);
853
854  while (non_zero != NULL)
855  {
856    assume(pGetExp(non_zero, n) != 0);
857    qq = non_zero;
858    pIter(non_zero);
859    qq->next = NULL;
860    pSetExp(qq,n,0);
861    pSetm(qq);
862    result = pAdd(result,qq);
863  }
864  p = pAdd(result, zero);
865  pTest(p);
866  return p;
867}
868
869/*3
870* substitute the n-th variable by number e in p
871* destroy p
872*/
873static poly pSubst2 (poly p,int n, number e)
874{
875  assume( ! nIsZero(e) );
876  poly qq,result = NULL;
877  number nn, nm;
878  poly zero, non_zero;
879
880  // reverse, so that add is likely to be linear
881  pSplitAndReversePoly(p, n, &non_zero, &zero);
882
883  while (non_zero != NULL)
884  {
885    assume(pGetExp(non_zero, n) != 0);
886    qq = non_zero;
887    pIter(non_zero);
888    qq->next = NULL;
889    nPower(e, pGetExp(qq, n), &nn);
890    nm = nMult(nn, pGetCoeff(qq));
891#ifdef HAVE_RINGS
892    if (nIsZero(nm))
893    {
894      pLmFree(&qq);
895      nDelete(&nm);
896    }
897    else
898#endif
899    {
900      pSetCoeff(qq, nm);
901      pSetExp(qq, n, 0);
902      pSetm(qq);
903      result = pAdd(result,qq);
904    }
905    nDelete(&nn);
906  }
907  p = pAdd(result, zero);
908  pTest(p);
909  return p;
910}
911
912
913/* delete monoms whose n-th exponent is different from zero */
914poly pSubst0(poly p, int n)
915{
916  spolyrec res;
917  poly h = &res;
918  pNext(h) = p;
919
920  while (pNext(h)!=NULL)
921  {
922    if (pGetExp(pNext(h),n)!=0)
923    {
924      pDeleteLm(&pNext(h));
925    }
926    else
927    {
928      pIter(h);
929    }
930  }
931  pTest(pNext(&res));
932  return pNext(&res);
933}
934
935/*2
936* substitute the n-th variable by e in p
937* destroy p
938*/
939poly pSubst(poly p, int n, poly e)
940{
941  if (e == NULL) return pSubst0(p, n);
942
943  if (pIsConstant(e))
944  {
945    if (nIsOne(pGetCoeff(e))) return pSubst1(p,n);
946    else return pSubst2(p, n, pGetCoeff(e));
947  }
948
949#ifdef HAVE_PLURAL
950  if (rIsPluralRing(currRing))
951  {
952    return nc_pSubst(p,n,e);
953  }
954#endif
955
956  int exponent,i;
957  poly h, res, m;
958  int *me,*ee;
959  number nu,nu1;
960
961  me=(int *)omAlloc((pVariables+1)*sizeof(int));
962  ee=(int *)omAlloc((pVariables+1)*sizeof(int));
963  if (e!=NULL) pGetExpV(e,ee);
964  res=NULL;
965  h=p;
966  while (h!=NULL)
967  {
968    if ((e!=NULL) || (pGetExp(h,n)==0))
969    {
970      m=pHead(h);
971      pGetExpV(m,me);
972      exponent=me[n];
973      me[n]=0;
974      for(i=pVariables;i>0;i--)
975        me[i]+=exponent*ee[i];
976      pSetExpV(m,me);
977      if (e!=NULL)
978      {
979        nPower(pGetCoeff(e),exponent,&nu);
980        nu1=nMult(pGetCoeff(m),nu);
981        nDelete(&nu);
982        pSetCoeff(m,nu1);
983      }
984      res=pAdd(res,m);
985    }
986    pDeleteLm(&h);
987  }
988  omFreeSize((ADDRESS)me,(pVariables+1)*sizeof(int));
989  omFreeSize((ADDRESS)ee,(pVariables+1)*sizeof(int));
990  return res;
991}
992
993/* Returns TRUE if
994     * LM(p) | LM(lcm)
995     * LC(p) | LC(lcm) only if ring
996     * Exists i, j:
997         * LE(p, i)  != LE(lcm, i)
998         * LE(p1, i) != LE(lcm, i)   ==> LCM(p1, p) != lcm
999         * LE(p, j)  != LE(lcm, j)
1000         * LE(p2, j) != LE(lcm, j)   ==> LCM(p2, p) != lcm
1001*/
1002BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
1003{
1004  int k, j;
1005
1006  if (lcm==NULL) return FALSE;
1007
1008  for (j=pVariables; j; j--)
1009    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
1010  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1011  for (j=pVariables; j; j--)
1012  {
1013    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1014    {
1015      if (pGetExp(p,j)!=pGetExp(lcm,j))
1016      {
1017        for (k=pVariables; k>j; k--)
1018        {
1019          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1020          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1021            return TRUE;
1022        }
1023        for (k=j-1; k; k--)
1024        {
1025          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1026          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1027            return TRUE;
1028        }
1029        return FALSE;
1030      }
1031    }
1032    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1033    {
1034      if (pGetExp(p,j)!=pGetExp(lcm,j))
1035      {
1036        for (k=pVariables; k>j; k--)
1037        {
1038          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1039          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1040            return TRUE;
1041        }
1042        for (k=j-1; k!=0 ; k--)
1043        {
1044          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1045          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1046            return TRUE;
1047        }
1048        return FALSE;
1049      }
1050    }
1051  }
1052  return FALSE;
1053}
1054#ifdef HAVE_RATGRING
1055BOOLEAN pCompareChainPart (poly p,poly p1,poly p2,poly lcm)
1056{
1057  int k, j;
1058
1059  if (lcm==NULL) return FALSE;
1060
1061  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1062    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
1063  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1064  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1065  {
1066    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1067    {
1068      if (pGetExp(p,j)!=pGetExp(lcm,j))
1069      {
1070        for (k=pVariables; k>j; k--)
1071        for (k=currRing->real_var_end; k>j; k--)
1072        {
1073          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1074          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1075            return TRUE;
1076        }
1077        for (k=j-1; k>=currRing->real_var_start; k--)
1078        {
1079          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1080          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1081            return TRUE;
1082        }
1083        return FALSE;
1084      }
1085    }
1086    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1087    {
1088      if (pGetExp(p,j)!=pGetExp(lcm,j))
1089      {
1090        for (k=currRing->real_var_end; k>j; k--)
1091        {
1092          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1093          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1094            return TRUE;
1095        }
1096        for (k=j-1; k>=currRing->real_var_start; k--)
1097        {
1098          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1099          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1100            return TRUE;
1101        }
1102        return FALSE;
1103      }
1104    }
1105  }
1106  return FALSE;
1107}
1108#endif
1109
1110int pSize(poly p)
1111{
1112  int count = 0;
1113  while ( p != NULL )
1114  {
1115    count+= nSize( pGetCoeff( p ) );
1116    pIter( p );
1117  }
1118  return count;
1119}
1120
1121/*2
1122* returns the length of a (numbers of monomials)
1123* respect syzComp
1124*/
1125poly pLast(poly a, int &l)
1126{
1127  if (a == NULL)
1128  {
1129    l = 0;
1130    return NULL;
1131  }
1132  l = 1;
1133  if (! rIsSyzIndexRing(currRing))
1134  {
1135    while (pNext(a)!=NULL)
1136    {
1137      pIter(a);
1138      l++;
1139    }
1140  }
1141  else
1142  {
1143    int curr_limit = rGetCurrSyzLimit(currRing);
1144    poly pp = a;
1145    while ((a=pNext(a))!=NULL)
1146    {
1147      if (pGetComp(a)<=curr_limit/*syzComp*/)
1148        l++;
1149      else break;
1150      pp = a;
1151    }
1152    a=pp;
1153  }
1154  return a;
1155}
1156
Note: See TracBrowser for help on using the repository browser.