source: git/kernel/polys.cc @ 18ff4c

spielwiese
Last change on this file since 18ff4c was 7ba059, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes: NF, pSize, normalize etc. git-svn-id: file:///usr/local/Singular/svn/trunk@10241 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.26 2007-07-24 12:29:31 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  pDiv_nn(result,inv);
177#else
178  pMult_nn(result,inv);
179#endif
180  nDelete(&inv);
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*/
202char * p_Read(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  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      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(char *st, BOOLEAN &ok)
275{
276  poly p;
277  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        c = nDiv(pGetCoeff(p1), k);
733        nDelete(&pGetCoeff(p1));
734        pSetCoeff0(p1, c);
735        h = pNext(p1);
736        while (h != NULL)
737        {
738          c = nDiv(pGetCoeff(h), k);
739          nDelete(&pGetCoeff(h));
740          pSetCoeff0(h, c);
741          pIter(h);
742        }
743      }
744      nDelete(&k);
745    }
746    return;
747  }
748  else
749#endif
750  if (p1!=NULL)
751  {
752    if (pNext(p1)==NULL)
753    {
754      pSetCoeff(p1,nInit(1));
755      return;
756    }
757    if (!nIsOne(pGetCoeff(p1)))
758    {
759      nNormalize(pGetCoeff(p1));
760      k = pGetCoeff(p1);
761      c = nInit(1);
762      pSetCoeff0(p1,c);
763      h = pNext(p1);
764      while (h!=NULL)
765      {
766        c=nDiv(pGetCoeff(h),k);
767        // no need to normalize: Z/p, R
768        // normalize already in nDiv: Q_a, Z/p_a
769        // remains: Q
770        if (rField_is_Q() && (!nIsOne(c))) nNormalize(c);
771        pSetCoeff(h,c);
772        pIter(h);
773      }
774      nDelete(&k);
775    }
776    else
777    {
778      if (nNormalize != nDummy2)
779      {
780        h = pNext(p1);
781        while (h!=NULL)
782        {
783          nNormalize(pGetCoeff(h));
784          pIter(h);
785        }
786      }
787    }
788  }
789}
790
791/*2
792*normalize all coefficients
793*/
794void p_Normalize(poly p, ring r)
795{
796  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
797  while (p!=NULL)
798  {
799    if (currRing==r) {nTest(pGetCoeff(p));}
800    n_Normalize(pGetCoeff(p),r);
801    pIter(p);
802  }
803}
804
805// splits p into polys with Exp(n) == 0 and Exp(n) != 0
806// Poly with Exp(n) != 0 is reversed
807static void pSplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero)
808{
809  if (p == NULL)
810  {
811    *non_zero = NULL;
812    *zero = NULL;
813    return;
814  }
815  spolyrec sz;
816  poly z, n_z, next;
817  z = &sz;
818  n_z = NULL;
819
820  while(p != NULL)
821  {
822    next = pNext(p);
823    if (pGetExp(p, n) == 0)
824    {
825      pNext(z) = p;
826      pIter(z);
827    }
828    else
829    {
830      pNext(p) = n_z;
831      n_z = p;
832    }
833    p = next;
834  }
835  pNext(z) = NULL;
836  *zero = pNext(&sz);
837  *non_zero = n_z;
838  return;
839}
840
841/*3
842* substitute the n-th variable by 1 in p
843* destroy p
844*/
845static poly pSubst1 (poly p,int n)
846{
847  poly qq=NULL, result = NULL;
848  poly zero=NULL, non_zero=NULL;
849
850  // reverse, so that add is likely to be linear
851  pSplitAndReversePoly(p, n, &non_zero, &zero);
852
853  while (non_zero != NULL)
854  {
855    assume(pGetExp(non_zero, n) != 0);
856    qq = non_zero;
857    pIter(non_zero);
858    qq->next = NULL;
859    pSetExp(qq,n,0);
860    pSetm(qq);
861    result = pAdd(result,qq);
862  }
863  p = pAdd(result, zero);
864  pTest(p);
865  return p;
866}
867
868/*3
869* substitute the n-th variable by number e in p
870* destroy p
871*/
872static poly pSubst2 (poly p,int n, number e)
873{
874  assume( ! nIsZero(e) );
875  poly qq,result = NULL;
876  number nn, nm;
877  poly zero, non_zero;
878
879  // reverse, so that add is likely to be linear
880  pSplitAndReversePoly(p, n, &non_zero, &zero);
881
882  while (non_zero != NULL)
883  {
884    assume(pGetExp(non_zero, n) != 0);
885    qq = non_zero;
886    pIter(non_zero);
887    qq->next = NULL;
888    nPower(e, pGetExp(qq, n), &nn);
889    nm = nMult(nn, pGetCoeff(qq));
890    pSetCoeff(qq, nm);
891    nDelete(&nn);
892    pSetExp(qq, n, 0);
893    pSetm(qq);
894    result = pAdd(result,qq);
895  }
896  p = pAdd(result, zero);
897  pTest(p);
898  return p;
899}
900
901
902/* delete monoms whose n-th exponent is different from zero */
903poly pSubst0(poly p, int n)
904{
905  spolyrec res;
906  poly h = &res;
907  pNext(h) = p;
908
909  while (pNext(h)!=NULL)
910  {
911    if (pGetExp(pNext(h),n)!=0)
912    {
913      pDeleteLm(&pNext(h));
914    }
915    else
916    {
917      pIter(h);
918    }
919  }
920  pTest(pNext(&res));
921  return pNext(&res);
922}
923
924/*2
925* substitute the n-th variable by e in p
926* destroy p
927*/
928poly pSubst(poly p, int n, poly e)
929{
930  if (e == NULL) return pSubst0(p, n);
931
932  if (pIsConstant(e))
933  {
934    if (nIsOne(pGetCoeff(e))) return pSubst1(p,n);
935    else return pSubst2(p, n, pGetCoeff(e));
936  }
937
938#ifdef HAVE_PLURAL
939  if (rIsPluralRing(currRing))
940  {
941    return nc_pSubst(p,n,e);
942  }
943#endif
944
945  int exponent,i;
946  poly h, res, m;
947  int *me,*ee;
948  number nu,nu1;
949
950  me=(int *)omAlloc((pVariables+1)*sizeof(int));
951  ee=(int *)omAlloc((pVariables+1)*sizeof(int));
952  if (e!=NULL) pGetExpV(e,ee);
953  res=NULL;
954  h=p;
955  while (h!=NULL)
956  {
957    if ((e!=NULL) || (pGetExp(h,n)==0))
958    {
959      m=pHead(h);
960      pGetExpV(m,me);
961      exponent=me[n];
962      me[n]=0;
963      for(i=pVariables;i>0;i--)
964        me[i]+=exponent*ee[i];
965      pSetExpV(m,me);
966      if (e!=NULL)
967      {
968        nPower(pGetCoeff(e),exponent,&nu);
969        nu1=nMult(pGetCoeff(m),nu);
970        nDelete(&nu);
971        pSetCoeff(m,nu1);
972      }
973      res=pAdd(res,m);
974    }
975    pDeleteLm(&h);
976  }
977  omFreeSize((ADDRESS)me,(pVariables+1)*sizeof(int));
978  omFreeSize((ADDRESS)ee,(pVariables+1)*sizeof(int));
979  return res;
980}
981
982/* Returns TRUE if
983     * LM(p) | LM(lcm)
984     * LC(p) | LC(lcm) only if ring
985     * Exists i, j:
986         * LE(p, i)  != LE(lcm, i)
987         * LE(p1, i) != LE(lcm, i)   ==> LCM(p1, p) != lcm
988         * LE(p, j)  != LE(lcm, j)
989         * LE(p2, j) != LE(lcm, j)   ==> LCM(p2, p) != lcm
990*/
991BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
992{
993  int k, j;
994
995  if (lcm==NULL) return FALSE;
996
997  for (j=pVariables; j; j--)
998    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
999  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1000  for (j=pVariables; j; j--)
1001  {
1002    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1003    {
1004      if (pGetExp(p,j)!=pGetExp(lcm,j))
1005      {
1006        for (k=pVariables; k>j; k--)
1007        {
1008          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1009          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1010            return TRUE;
1011        }
1012        for (k=j-1; k; k--)
1013        {
1014          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1015          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1016            return TRUE;
1017        }
1018        return FALSE;
1019      }
1020    }
1021    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1022    {
1023      if (pGetExp(p,j)!=pGetExp(lcm,j))
1024      {
1025        for (k=pVariables; k>j; k--)
1026        {
1027          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1028          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1029            return TRUE;
1030        }
1031        for (k=j-1; k!=0 ; k--)
1032        {
1033          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1034          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1035            return TRUE;
1036        }
1037        return FALSE;
1038      }
1039    }
1040  }
1041  return FALSE;
1042}
1043
1044int pSize(poly p)
1045{
1046  int count = 0;
1047  while ( p != NULL )
1048  {
1049    count+= nSize( pGetCoeff( p ) );
1050    pIter( p );
1051  }
1052  return count;
1053}
1054
Note: See TracBrowser for help on using the repository browser.