source: git/kernel/polys.cc @ b1f3b55

spielwiese
Last change on this file since b1f3b55 was b1f3b55, checked in by Oliver Wienand <wienand@…>, 16 years ago
polys.cc: --> Memory Fehler git-svn-id: file:///usr/local/Singular/svn/trunk@10787 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 19.7 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: polys.cc,v 1.29 2008-06-24 16:25:02 wienand 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  pDiv_nn(result,inv);
177#else
178  pMult_nn(result,inv);
179  nDelete(&inv);
180#endif
181  pDelete(&b);
182  return result;
183}
184
185/*2
186* returns the LCM of the head terms of a and b in *m
187*/
188void pLcm(poly a, poly b, poly m)
189{
190  int i;
191  for (i=pVariables; i; i--)
192  {
193    pSetExp(m,i, si_max( pGetExp(a,i), pGetExp(b,i)));
194  }
195  pSetComp(m, si_max(pGetComp(a), pGetComp(b)));
196  /* Don't do a pSetm here, otherwise hres/lres chockes */
197}
198
199/*2
200* convert monomial given as string to poly, e.g. 1x3y5z
201*/
202const char * p_Read(const char *st, poly &rc, ring r)
203{
204  if (r==NULL) { rc=NULL;return st;}
205  int i,j;
206  rc = p_Init(r);
207  const char *s = r->cf->nRead(st,&(rc->coef));
208  if (s==st)
209  /* i.e. it does not start with a coeff: test if it is a ringvar*/
210  {
211    j = r_IsRingVar(s,r);
212    if (j >= 0)
213    {
214      p_IncrExp(rc,1+j,r);
215      while (*s!='\0') s++;
216      goto done;
217    }
218  }
219  while (*s!='\0')
220  {
221    char ss[2];
222    ss[0] = *s++;
223    ss[1] = '\0';
224    j = r_IsRingVar(ss,r);
225    if (j >= 0)
226    {
227      const char *s_save=s;
228      s = eati(s,&i);
229      if (((unsigned long)i) >  r->bitmask)
230      {
231        // exponent to large: it is not a monomial
232        p_DeleteLm(&rc,r);
233        return s_save;
234      }
235      p_AddExp(rc,1+j, (Exponent_t)i, r);
236    }
237    else
238    {
239      // 1st char of is not a varname
240      p_DeleteLm(&rc,r);
241      s--;
242      return s;
243    }
244  }
245done:
246  if (r->cf->nIsZero(pGetCoeff(rc))) p_DeleteLm(&rc,r);
247  else
248  {
249#ifdef HAVE_PLURAL
250    // in super-commutative ring
251    // squares of anti-commutative variables are zeroes!
252    if(rIsSCA(r))
253    {
254      const unsigned int iFirstAltVar = scaFirstAltVar(r);
255      const unsigned int iLastAltVar  = scaLastAltVar(r);
256
257      assume(rc != NULL);
258
259      for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
260        if( p_GetExp(rc, k, r) > 1 )
261        {
262          p_DeleteLm(&rc, r);
263          goto finish;
264        }
265    }
266#endif
267   
268    p_Setm(rc,r);
269  }
270finish: 
271  return s;
272}
273
274poly pmInit(const char *st, BOOLEAN &ok)
275{
276  poly p;
277  const char *s=p_Read(st,p,currRing);
278  if (*s!='\0')
279  {
280    if ((s!=st)&&isdigit(st[0]))
281    {
282      errorreported=TRUE;
283    }
284    ok=FALSE;
285    pDelete(&p);
286    return NULL;
287  }
288  ok=!errorreported;
289  return p;
290}
291
292/*2
293*make p homogeneous by multiplying the monomials by powers of x_varnum
294*assume: deg(var(varnum))==1
295*/
296poly pHomogen (poly p, int varnum)
297{
298  poly q=NULL, qn;
299  int  o,ii;
300  sBucket_pt bp;
301
302  if (p!=NULL)
303  {
304    if ((varnum < 1) || (varnum > pVariables))
305    {
306      return NULL;
307    }
308    o=pWTotaldegree(p);
309    q=pNext(p);
310    while (q != NULL)
311    {
312      ii=pWTotaldegree(q);
313      if (ii>o) o=ii;
314      pIter(q);
315    }
316    q = pCopy(p);
317    bp = sBucketCreate(currRing);
318    while (q != NULL)
319    {
320      ii = o-pWTotaldegree(q);
321      if (ii!=0)
322      {
323        pAddExp(q,varnum, (Exponent_t)ii);
324        pSetm(q);
325      }
326      qn = pNext(q);
327      pNext(q) = NULL;
328      sBucket_Add_p(bp, q, 1);
329      q = qn;
330    }
331    sBucketDestroyAdd(bp, &q, &ii);
332  }
333  return q;
334}
335
336/*2
337*replaces the maximal powers of the leading monomial of p2 in p1 by
338*the same powers of n, utility for dehomogenization
339*/
340poly pDehomogen (poly p1,poly p2,number n)
341{
342  polyset P;
343  int     SizeOfSet=5;
344  int     i;
345  poly    p;
346  number  nn;
347
348  P = (polyset)omAlloc0(5*sizeof(poly));
349  //for (i=0; i<5; i++)
350  //{
351  //  P[i] = NULL;
352  //}
353  pCancelPolyByMonom(p1,p2,&P,&SizeOfSet);
354  p = P[0];
355  //P[0] = NULL ;// for safety, may be removed later
356  for (i=1; i<SizeOfSet; i++)
357  {
358    if (P[i] != NULL)
359    {
360      nPower(n,i,&nn);
361      pMult_nn(P[i],nn);
362      p = pAdd(p,P[i]);
363      //P[i] =NULL; // for safety, may be removed later
364      nDelete(&nn);
365    }
366  }
367  omFreeSize((ADDRESS)P,SizeOfSet*sizeof(poly));
368  return p;
369}
370
371/*4
372*Returns the exponent of the maximal power of the leading monomial of
373*p2 in that of p1
374*/
375static int pGetMaxPower (poly p1,poly p2)
376{
377  int     i,k,res = 32000; /*a very large integer*/
378
379  if (p1 == NULL) return 0;
380  for (i=1; i<=pVariables; i++)
381  {
382    if ( pGetExp(p2,i) != 0)
383    {
384      k =  pGetExp(p1,i) /  pGetExp(p2,i);
385      if (k < res) res = k;
386    }
387  }
388  return res;
389}
390
391/*2
392*Returns as i-th entry of P the coefficient of the (i-1) power of
393*the leading monomial of p2 in p1
394*/
395void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet)
396{
397  int   maxPow;
398  poly  p,qp,Coeff;
399
400  if (*P == NULL)
401  {
402    *P = (polyset) omAlloc(5*sizeof(poly));
403    *SizeOfSet = 5;
404  }
405  p = pCopy(p1);
406  while (p != NULL)
407  {
408    qp = p->next;
409    p->next = NULL;
410    maxPow = pGetMaxPower(p,p2);
411    Coeff = pDivByMonom(p,p2);
412    if (maxPow > *SizeOfSet)
413    {
414      pEnlargeSet(P,*SizeOfSet,maxPow+1-*SizeOfSet);
415      *SizeOfSet = maxPow+1;
416    }
417    (*P)[maxPow] = pAdd((*P)[maxPow],Coeff);
418    pDelete(&p);
419    p = qp;
420  }
421}
422
423/*2
424*returns the leading monomial of p1 divided by the maximal power of that
425*of p2
426*/
427poly pDivByMonom (poly p1,poly p2)
428{
429  int     k, i;
430
431  if (p1 == NULL) return NULL;
432  k = pGetMaxPower(p1,p2);
433  if (k == 0)
434    return pHead(p1);
435  else
436  {
437    number n;
438    poly p = pInit();
439
440    p->next = NULL;
441    for (i=1; i<=pVariables; i++)
442    {
443       pSetExp(p,i, pGetExp(p1,i)-k* pGetExp(p2,i));
444    }
445    nPower(p2->coef,k,&n);
446    pSetCoeff0(p,nDiv(p1->coef,n));
447    nDelete(&n);
448    pSetm(p);
449    return p;
450  }
451}
452/*----------utilities for syzygies--------------*/
453poly pTakeOutComp(poly * p, int k)
454{
455  poly q = *p,qq=NULL,result = NULL;
456
457  if (q==NULL) return NULL;
458  if (pGetComp(q)==k)
459  {
460    result = q;
461    do
462    {
463      pSetComp(q,0);
464      pSetmComp(q);
465      qq = q;
466      pIter(q);
467    }
468    while ((q!=NULL) && (pGetComp(q)==k));
469    *p = q;
470    pNext(qq) = NULL;
471  }
472  if (q==NULL) return result;
473  if (pGetComp(q) > k)
474  {
475    pDecrComp(q);
476    pSetmComp(q);
477  }
478  poly pNext_q;
479  while ((pNext_q=pNext(q))!=NULL)
480  {
481    if (pGetComp(pNext_q)==k)
482    {
483      if (result==NULL)
484      {
485        result = pNext_q;
486        qq = result;
487      }
488      else
489      {
490        pNext(qq) = pNext_q;
491        pIter(qq);
492      }
493      pNext(q) = pNext(pNext_q);
494      pNext(qq) =NULL;
495      pSetComp(qq,0);
496      pSetmComp(qq);
497    }
498    else
499    {
500      /*pIter(q);*/ q=pNext_q;
501      if (pGetComp(q) > k)
502      {
503        pDecrComp(q);
504        pSetmComp(q);
505      }
506    }
507  }
508  return result;
509}
510
511// Splits *p into two polys: *q which consists of all monoms with
512// component == comp and *p of all other monoms *lq == pLength(*q)
513void pTakeOutComp(poly *r_p, Exponent_t comp, poly *r_q, int *lq)
514{
515  spolyrec pp, qq;
516  poly p, q, p_prev;
517  int l = 0;
518
519#ifdef HAVE_ASSUME
520  int lp = pLength(*r_p);
521#endif
522
523  pNext(&pp) = *r_p;
524  p = *r_p;
525  p_prev = &pp;
526  q = &qq;
527
528  while(p != NULL)
529  {
530    while (pGetComp(p) == comp)
531    {
532      pNext(q) = p;
533      pIter(q);
534      pSetComp(p, 0);
535      pSetmComp(p);
536      pIter(p);
537      l++;
538      if (p == NULL)
539      {
540        pNext(p_prev) = NULL;
541        goto Finish;
542      }
543    }
544    pNext(p_prev) = p;
545    p_prev = p;
546    pIter(p);
547  }
548
549  Finish:
550  pNext(q) = NULL;
551  *r_p = pNext(&pp);
552  *r_q = pNext(&qq);
553  *lq = l;
554#ifdef HAVE_ASSUME
555  assume(pLength(*r_p) + pLength(*r_q) == lp);
556#endif
557  pTest(*r_p);
558  pTest(*r_q);
559}
560
561void pDecrOrdTakeOutComp(poly *r_p, Exponent_t comp, Order_t order,
562                         poly *r_q, int *lq)
563{
564  spolyrec pp, qq;
565  poly p, q, p_prev;
566  int l = 0;
567
568  pNext(&pp) = *r_p;
569  p = *r_p;
570  p_prev = &pp;
571  q = &qq;
572
573#ifdef HAVE_ASSUME
574  if (p != NULL)
575  {
576    while (pNext(p) != NULL)
577    {
578      assume(pGetOrder(p) >= pGetOrder(pNext(p)));
579      pIter(p);
580    }
581  }
582  p = *r_p;
583#endif
584
585  while (p != NULL && pGetOrder(p) > order) pIter(p);
586
587  while(p != NULL && pGetOrder(p) == order)
588  {
589    while (pGetComp(p) == comp)
590    {
591      pNext(q) = p;
592      pIter(q);
593      pIter(p);
594      pSetComp(p, 0);
595      pSetmComp(p);
596      l++;
597      if (p == NULL || pGetOrder(p) != order)
598      {
599        pNext(p_prev) = p;
600        goto Finish;
601      }
602    }
603    pNext(p_prev) = p;
604    p_prev = p;
605    pIter(p);
606  }
607
608  Finish:
609  pNext(q) = NULL;
610  *r_p = pNext(&pp);
611  *r_q = pNext(&qq);
612  *lq = l;
613}
614
615#if 1
616poly pTakeOutComp1(poly * p, int k)
617{
618  poly q = *p;
619
620  if (q==NULL) return NULL;
621
622  poly qq=NULL,result = NULL;
623
624  if (pGetComp(q)==k)
625  {
626    result = q; /* *p */
627    while ((q!=NULL) && (pGetComp(q)==k))
628    {
629      pSetComp(q,0);
630      pSetmComp(q);
631      qq = q;
632      pIter(q);
633    }
634    *p = q;
635    pNext(qq) = NULL;
636  }
637  if (q==NULL) return result;
638//  if (pGetComp(q) > k) pGetComp(q)--;
639  while (pNext(q)!=NULL)
640  {
641    if (pGetComp(pNext(q))==k)
642    {
643      if (result==NULL)
644      {
645        result = pNext(q);
646        qq = result;
647      }
648      else
649      {
650        pNext(qq) = pNext(q);
651        pIter(qq);
652      }
653      pNext(q) = pNext(pNext(q));
654      pNext(qq) =NULL;
655      pSetComp(qq,0);
656      pSetmComp(qq);
657    }
658    else
659    {
660      pIter(q);
661//      if (pGetComp(q) > k) pGetComp(q)--;
662    }
663  }
664  return result;
665}
666#endif
667
668void pDeleteComp(poly * p,int k)
669{
670  poly q;
671
672  while ((*p!=NULL) && (pGetComp(*p)==k)) pDeleteLm(p);
673  if (*p==NULL) return;
674  q = *p;
675  if (pGetComp(q)>k)
676  {
677    pDecrComp(q);
678    pSetmComp(q);
679  }
680  while (pNext(q)!=NULL)
681  {
682    if (pGetComp(pNext(q))==k)
683      pDeleteLm(&(pNext(q)));
684    else
685    {
686      pIter(q);
687      if (pGetComp(q)>k)
688      {
689        pDecrComp(q);
690        pSetmComp(q);
691      }
692    }
693  }
694}
695/*----------end of utilities for syzygies--------------*/
696
697/*2
698* pair has no common factor ? or is no polynomial
699*/
700BOOLEAN pHasNotCF(poly p1, poly p2)
701{
702
703  if (!TEST_OPT_IDLIFT)
704  {
705    if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
706      return FALSE;
707  }
708  int i = 1;
709  loop
710  {
711    if ((pGetExp(p1, i) > 0) && (pGetExp(p2, i) > 0))   return FALSE;
712    if (i == pVariables)                                return TRUE;
713    i++;
714  }
715}
716
717/*2
718*divides p1 by its leading coefficient
719*/
720void pNorm(poly p1)
721{
722  poly h;
723  number k, c;
724#ifdef HAVE_RINGS
725  if (rField_is_Ring(currRing))
726  {
727    if ((p1!=NULL) && rField_has_Units(currRing))
728    {
729      k = nGetUnit(pGetCoeff(p1));
730      if (!nIsOne(k))
731      {
732        k = nInvers(k);
733        c = nMult(pGetCoeff(p1), k); 
734        nDelete(&pGetCoeff(p1));
735        pSetCoeff0(p1, c);
736        h = pNext(p1);
737        while (h != NULL)
738        {
739          c = nMult(pGetCoeff(h), k);
740          nDelete(&pGetCoeff(h));
741          pSetCoeff0(h, c);
742          pIter(h);
743        }
744      }
745      nDelete(&k);
746    }
747    return;
748  }
749  else
750#endif
751  if (p1!=NULL)
752  {
753    if (pNext(p1)==NULL)
754    {
755      pSetCoeff(p1,nInit(1));
756      return;
757    }
758    if (!nIsOne(pGetCoeff(p1)))
759    {
760      nNormalize(pGetCoeff(p1));
761      k = pGetCoeff(p1);
762      c = nInit(1);
763      pSetCoeff0(p1,c);
764      h = pNext(p1);
765      while (h!=NULL)
766      {
767        c=nDiv(pGetCoeff(h),k);
768        // no need to normalize: Z/p, R
769        // normalize already in nDiv: Q_a, Z/p_a
770        // remains: Q
771        if (rField_is_Q() && (!nIsOne(c))) nNormalize(c);
772        pSetCoeff(h,c);
773        pIter(h);
774      }
775      nDelete(&k);
776    }
777    else
778    {
779      if (nNormalize != nDummy2)
780      {
781        h = pNext(p1);
782        while (h!=NULL)
783        {
784          nNormalize(pGetCoeff(h));
785          pIter(h);
786        }
787      }
788    }
789  }
790}
791
792/*2
793*normalize all coefficients
794*/
795void p_Normalize(poly p, ring r)
796{
797  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
798  while (p!=NULL)
799  {
800    if (currRing==r) {nTest(pGetCoeff(p));}
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    pSetCoeff(qq, nm);
892    nDelete(&nn);
893    pSetExp(qq, n, 0);
894    pSetm(qq);
895    result = pAdd(result,qq);
896  }
897  p = pAdd(result, zero);
898  pTest(p);
899  return p;
900}
901
902
903/* delete monoms whose n-th exponent is different from zero */
904poly pSubst0(poly p, int n)
905{
906  spolyrec res;
907  poly h = &res;
908  pNext(h) = p;
909
910  while (pNext(h)!=NULL)
911  {
912    if (pGetExp(pNext(h),n)!=0)
913    {
914      pDeleteLm(&pNext(h));
915    }
916    else
917    {
918      pIter(h);
919    }
920  }
921  pTest(pNext(&res));
922  return pNext(&res);
923}
924
925/*2
926* substitute the n-th variable by e in p
927* destroy p
928*/
929poly pSubst(poly p, int n, poly e)
930{
931  if (e == NULL) return pSubst0(p, n);
932
933  if (pIsConstant(e))
934  {
935    if (nIsOne(pGetCoeff(e))) return pSubst1(p,n);
936    else return pSubst2(p, n, pGetCoeff(e));
937  }
938
939#ifdef HAVE_PLURAL
940  if (rIsPluralRing(currRing))
941  {
942    return nc_pSubst(p,n,e);
943  }
944#endif
945
946  int exponent,i;
947  poly h, res, m;
948  int *me,*ee;
949  number nu,nu1;
950
951  me=(int *)omAlloc((pVariables+1)*sizeof(int));
952  ee=(int *)omAlloc((pVariables+1)*sizeof(int));
953  if (e!=NULL) pGetExpV(e,ee);
954  res=NULL;
955  h=p;
956  while (h!=NULL)
957  {
958    if ((e!=NULL) || (pGetExp(h,n)==0))
959    {
960      m=pHead(h);
961      pGetExpV(m,me);
962      exponent=me[n];
963      me[n]=0;
964      for(i=pVariables;i>0;i--)
965        me[i]+=exponent*ee[i];
966      pSetExpV(m,me);
967      if (e!=NULL)
968      {
969        nPower(pGetCoeff(e),exponent,&nu);
970        nu1=nMult(pGetCoeff(m),nu);
971        nDelete(&nu);
972        pSetCoeff(m,nu1);
973      }
974      res=pAdd(res,m);
975    }
976    pDeleteLm(&h);
977  }
978  omFreeSize((ADDRESS)me,(pVariables+1)*sizeof(int));
979  omFreeSize((ADDRESS)ee,(pVariables+1)*sizeof(int));
980  return res;
981}
982
983/* Returns TRUE if
984     * LM(p) | LM(lcm)
985     * LC(p) | LC(lcm) only if ring
986     * Exists i, j:
987         * LE(p, i)  != LE(lcm, i)
988         * LE(p1, i) != LE(lcm, i)   ==> LCM(p1, p) != lcm
989         * LE(p, j)  != LE(lcm, j)
990         * LE(p2, j) != LE(lcm, j)   ==> LCM(p2, p) != lcm
991*/
992BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
993{
994  int k, j;
995
996  if (lcm==NULL) return FALSE;
997
998  for (j=pVariables; j; j--)
999    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
1000  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1001  for (j=pVariables; j; j--)
1002  {
1003    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1004    {
1005      if (pGetExp(p,j)!=pGetExp(lcm,j))
1006      {
1007        for (k=pVariables; k>j; k--)
1008        {
1009          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1010          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1011            return TRUE;
1012        }
1013        for (k=j-1; k; k--)
1014        {
1015          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1016          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1017            return TRUE;
1018        }
1019        return FALSE;
1020      }
1021    }
1022    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1023    {
1024      if (pGetExp(p,j)!=pGetExp(lcm,j))
1025      {
1026        for (k=pVariables; k>j; k--)
1027        {
1028          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1029          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1030            return TRUE;
1031        }
1032        for (k=j-1; k!=0 ; k--)
1033        {
1034          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1035          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1036            return TRUE;
1037        }
1038        return FALSE;
1039      }
1040    }
1041  }
1042  return FALSE;
1043}
1044
1045int pSize(poly p)
1046{
1047  int count = 0;
1048  while ( p != NULL )
1049  {
1050    count+= nSize( pGetCoeff( p ) );
1051    pIter( p );
1052  }
1053  return count;
1054}
1055
Note: See TracBrowser for help on using the repository browser.