source: git/kernel/polys.cc @ 5ef9d3

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