source: git/kernel/polys.cc @ d681e8

fieker-DuValspielwiese
Last change on this file since d681e8 was d681e8, checked in by Oliver Wienand <wienand@…>, 17 years ago
kInline.cc, kutil.h: new: k_GetStringLeadTerms kutil.cc: lots of KDEBUG messages inside of HAVE_RINGS switches added ntl xgcd import new procs: enterOneStrongPoly initenterstrongPairs rmodulon.cc: nrnLcm small error Question: Is x == modul Then return 0 faster than return x % modul? nrnComp missing case corrected intDiv(0,0) = 0, intDiv(0,1) = 1, intDiv(0, a) = modul DIV a rmodulo2m.cc: intDiv(0,0) = 0, intDiv(0,1) = 1, intDiv(0, a) = modul DIV a git-svn-id: file:///usr/local/Singular/svn/trunk@10051 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 19.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: polys.cc,v 1.18 2007-05-19 13:22:23 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#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*/
270poly pHomogen (poly p, int varnum)
271{
272  poly q=NULL, qn;
273  int  o,ii;
274  sBucket_pt bp;
275
276  if (p!=NULL)
277  {
278    if ((varnum < 1) || (varnum > pVariables))
279    {
280      return NULL;
281    }
282    o=pWTotaldegree(p);
283    q=pNext(p);
284    while (q != NULL)
285    {
286      ii=pWTotaldegree(q);
287      if (ii>o) o=ii;
288      pIter(q);
289    }
290    q = pCopy(p);
291    bp = sBucketCreate(currRing);
292    while (q != NULL)
293    {
294      ii = o-pWTotaldegree(q);
295      if (ii!=0)
296      {
297        pAddExp(q,varnum, (Exponent_t)ii);
298        pSetm(q);
299      }
300      qn = pNext(q);
301      pNext(q) = NULL;
302      sBucket_Add_p(bp, q, 1);
303      q = qn;
304    }
305    sBucketDestroyAdd(bp, &q, &ii);
306  }
307  return q;
308}
309
310/*2
311*replaces the maximal powers of the leading monomial of p2 in p1 by
312*the same powers of n, utility for dehomogenization
313*/
314poly pDehomogen (poly p1,poly p2,number n)
315{
316  polyset P;
317  int     SizeOfSet=5;
318  int     i;
319  poly    p;
320  number  nn;
321
322  P = (polyset)omAlloc0(5*sizeof(poly));
323  //for (i=0; i<5; i++)
324  //{
325  //  P[i] = NULL;
326  //}
327  pCancelPolyByMonom(p1,p2,&P,&SizeOfSet);
328  p = P[0];
329  //P[0] = NULL ;// for safety, may be removed later
330  for (i=1; i<SizeOfSet; i++)
331  {
332    if (P[i] != NULL)
333    {
334      nPower(n,i,&nn);
335      pMult_nn(P[i],nn);
336      p = pAdd(p,P[i]);
337      //P[i] =NULL; // for safety, may be removed later
338      nDelete(&nn);
339    }
340  }
341  omFreeSize((ADDRESS)P,SizeOfSet*sizeof(poly));
342  return p;
343}
344
345/*4
346*Returns the exponent of the maximal power of the leading monomial of
347*p2 in that of p1
348*/
349static int pGetMaxPower (poly p1,poly p2)
350{
351  int     i,k,res = 32000; /*a very large integer*/
352
353  if (p1 == NULL) return 0;
354  for (i=1; i<=pVariables; i++)
355  {
356    if ( pGetExp(p2,i) != 0)
357    {
358      k =  pGetExp(p1,i) /  pGetExp(p2,i);
359      if (k < res) res = k;
360    }
361  }
362  return res;
363}
364
365/*2
366*Returns as i-th entry of P the coefficient of the (i-1) power of
367*the leading monomial of p2 in p1
368*/
369void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet)
370{
371  int   maxPow;
372  poly  p,qp,Coeff;
373
374  if (*P == NULL)
375  {
376    *P = (polyset) omAlloc(5*sizeof(poly));
377    *SizeOfSet = 5;
378  }
379  p = pCopy(p1);
380  while (p != NULL)
381  {
382    qp = p->next;
383    p->next = NULL;
384    maxPow = pGetMaxPower(p,p2);
385    Coeff = pDivByMonom(p,p2);
386    if (maxPow > *SizeOfSet)
387    {
388      pEnlargeSet(P,*SizeOfSet,maxPow+1-*SizeOfSet);
389      *SizeOfSet = maxPow+1;
390    }
391    (*P)[maxPow] = pAdd((*P)[maxPow],Coeff);
392    pDelete(&p);
393    p = qp;
394  }
395}
396
397/*2
398*returns the leading monomial of p1 divided by the maximal power of that
399*of p2
400*/
401poly pDivByMonom (poly p1,poly p2)
402{
403  int     k, i;
404
405  if (p1 == NULL) return NULL;
406  k = pGetMaxPower(p1,p2);
407  if (k == 0)
408    return pHead(p1);
409  else
410  {
411    number n;
412    poly p = pInit();
413
414    p->next = NULL;
415    for (i=1; i<=pVariables; i++)
416    {
417       pSetExp(p,i, pGetExp(p1,i)-k* pGetExp(p2,i));
418    }
419    nPower(p2->coef,k,&n);
420    pSetCoeff0(p,nDiv(p1->coef,n));
421    nDelete(&n);
422    pSetm(p);
423    return p;
424  }
425}
426/*----------utilities for syzygies--------------*/
427poly pTakeOutComp(poly * p, int k)
428{
429  poly q = *p,qq=NULL,result = NULL;
430
431  if (q==NULL) return NULL;
432  if (pGetComp(q)==k)
433  {
434    result = q;
435    do
436    {
437      pSetComp(q,0);
438      pSetmComp(q);
439      qq = q;
440      pIter(q);
441    }
442    while ((q!=NULL) && (pGetComp(q)==k));
443    *p = q;
444    pNext(qq) = NULL;
445  }
446  if (q==NULL) return result;
447  if (pGetComp(q) > k)
448  {
449    pDecrComp(q);
450    pSetmComp(q);
451  }
452  poly pNext_q;
453  while ((pNext_q=pNext(q))!=NULL)
454  {
455    if (pGetComp(pNext_q)==k)
456    {
457      if (result==NULL)
458      {
459        result = pNext_q;
460        qq = result;
461      }
462      else
463      {
464        pNext(qq) = pNext_q;
465        pIter(qq);
466      }
467      pNext(q) = pNext(pNext_q);
468      pNext(qq) =NULL;
469      pSetComp(qq,0);
470      pSetmComp(qq);
471    }
472    else
473    {
474      /*pIter(q);*/ q=pNext_q;
475      if (pGetComp(q) > k)
476      {
477        pDecrComp(q);
478        pSetmComp(q);
479      }
480    }
481  }
482  return result;
483}
484
485// Splits *p into two polys: *q which consists of all monoms with
486// component == comp and *p of all other monoms *lq == pLength(*q)
487void pTakeOutComp(poly *r_p, Exponent_t comp, poly *r_q, int *lq)
488{
489  spolyrec pp, qq;
490  poly p, q, p_prev;
491  int l = 0;
492
493#ifdef HAVE_ASSUME
494  int lp = pLength(*r_p);
495#endif
496
497  pNext(&pp) = *r_p;
498  p = *r_p;
499  p_prev = &pp;
500  q = &qq;
501
502  while(p != NULL)
503  {
504    while (pGetComp(p) == comp)
505    {
506      pNext(q) = p;
507      pIter(q);
508      pSetComp(p, 0);
509      pSetmComp(p);
510      pIter(p);
511      l++;
512      if (p == NULL)
513      {
514        pNext(p_prev) = NULL;
515        goto Finish;
516      }
517    }
518    pNext(p_prev) = p;
519    p_prev = p;
520    pIter(p);
521  }
522
523  Finish:
524  pNext(q) = NULL;
525  *r_p = pNext(&pp);
526  *r_q = pNext(&qq);
527  *lq = l;
528#ifdef HAVE_ASSUME
529  assume(pLength(*r_p) + pLength(*r_q) == lp);
530#endif
531  pTest(*r_p);
532  pTest(*r_q);
533}
534
535void pDecrOrdTakeOutComp(poly *r_p, Exponent_t comp, Order_t order,
536                         poly *r_q, int *lq)
537{
538  spolyrec pp, qq;
539  poly p, q, p_prev;
540  int l = 0;
541
542  pNext(&pp) = *r_p;
543  p = *r_p;
544  p_prev = &pp;
545  q = &qq;
546
547#ifdef HAVE_ASSUME
548  if (p != NULL)
549  {
550    while (pNext(p) != NULL)
551    {
552      assume(pGetOrder(p) >= pGetOrder(pNext(p)));
553      pIter(p);
554    }
555  }
556  p = *r_p;
557#endif
558
559  while (p != NULL && pGetOrder(p) > order) pIter(p);
560
561  while(p != NULL && pGetOrder(p) == order)
562  {
563    while (pGetComp(p) == comp)
564    {
565      pNext(q) = p;
566      pIter(q);
567      pIter(p);
568      pSetComp(p, 0);
569      pSetmComp(p);
570      l++;
571      if (p == NULL || pGetOrder(p) != order)
572      {
573        pNext(p_prev) = p;
574        goto Finish;
575      }
576    }
577    pNext(p_prev) = p;
578    p_prev = p;
579    pIter(p);
580  }
581
582  Finish:
583  pNext(q) = NULL;
584  *r_p = pNext(&pp);
585  *r_q = pNext(&qq);
586  *lq = l;
587}
588
589#if 1
590poly pTakeOutComp1(poly * p, int k)
591{
592  poly q = *p;
593
594  if (q==NULL) return NULL;
595
596  poly qq=NULL,result = NULL;
597
598  if (pGetComp(q)==k)
599  {
600    result = q; /* *p */
601    while ((q!=NULL) && (pGetComp(q)==k))
602    {
603      pSetComp(q,0);
604      pSetmComp(q);
605      qq = q;
606      pIter(q);
607    }
608    *p = q;
609    pNext(qq) = NULL;
610  }
611  if (q==NULL) return result;
612//  if (pGetComp(q) > k) pGetComp(q)--;
613  while (pNext(q)!=NULL)
614  {
615    if (pGetComp(pNext(q))==k)
616    {
617      if (result==NULL)
618      {
619        result = pNext(q);
620        qq = result;
621      }
622      else
623      {
624        pNext(qq) = pNext(q);
625        pIter(qq);
626      }
627      pNext(q) = pNext(pNext(q));
628      pNext(qq) =NULL;
629      pSetComp(qq,0);
630      pSetmComp(qq);
631    }
632    else
633    {
634      pIter(q);
635//      if (pGetComp(q) > k) pGetComp(q)--;
636    }
637  }
638  return result;
639}
640#endif
641
642void pDeleteComp(poly * p,int k)
643{
644  poly q;
645
646  while ((*p!=NULL) && (pGetComp(*p)==k)) pDeleteLm(p);
647  if (*p==NULL) return;
648  q = *p;
649  if (pGetComp(q)>k)
650  {
651    pDecrComp(q);
652    pSetmComp(q);
653  }
654  while (pNext(q)!=NULL)
655  {
656    if (pGetComp(pNext(q))==k)
657      pDeleteLm(&(pNext(q)));
658    else
659    {
660      pIter(q);
661      if (pGetComp(q)>k)
662      {
663        pDecrComp(q);
664        pSetmComp(q);
665      }
666    }
667  }
668}
669/*----------end of utilities for syzygies--------------*/
670
671/*2
672* pair has no common factor ? or is no polynomial
673*/
674BOOLEAN pHasNotCF(poly p1, poly p2)
675{
676
677  if (!TEST_OPT_IDLIFT)
678  {
679    if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
680      return FALSE;
681  }
682  int i = 1;
683  loop
684  {
685    if ((pGetExp(p1, i) > 0) && (pGetExp(p2, i) > 0))   return FALSE;
686    if (i == pVariables)                                return TRUE;
687    i++;
688  }
689}
690
691#ifdef HAVE_RINGS  //HACK TODO Oliver
692number nGetUnit(number k) {
693  number unit = nIntDiv(k, nGcd(k, 0, currRing));
694  number gcd = nGcd(unit, 0, currRing);
695  if (!nIsOne(gcd))
696  {
697    number tmp = nMult(unit, unit);
698    number gcd_new = nGcd(tmp, 0, currRing);
699    while (gcd_new != gcd)
700    {
701      gcd = gcd_new;
702      tmp = nMult(tmp, unit);
703      gcd_new = nGcd(tmp, 0, currRing);
704    }
705    unit = nAdd(unit, nIntDiv(0, gcd_new));
706  }
707  return unit;
708}
709#endif
710
711/*2
712*divides p1 by its leading coefficient
713*/
714void pNorm(poly p1)
715{
716  poly h;
717  number k, c;
718#ifdef HAVE_RINGS
719  if (rField_is_Ring(currRing))
720  {
721    if (p1!=NULL)
722    {
723      k = nGetUnit(pGetCoeff(p1));
724      if (!nIsOne(k))
725      {
726        c = nDiv(pGetCoeff(p1), k);
727        pSetCoeff0(p1, c);
728        h = pNext(p1);
729        while (h != NULL)
730        {
731          c = nDiv(pGetCoeff(h), k);
732          pSetCoeff(h, c);
733          pIter(h);
734        }
735        nDelete(&k);
736      }
737     return;
738    }
739  }
740  else
741#endif
742  if (p1!=NULL)
743  {
744    if (pNext(p1)==NULL)
745    {
746      pSetCoeff(p1,nInit(1));
747      return;
748    }
749    if (!nIsOne(pGetCoeff(p1)))
750    {
751      nNormalize(pGetCoeff(p1));
752      k = pGetCoeff(p1);
753      c = nInit(1);
754      pSetCoeff0(p1,c);
755      h = pNext(p1);
756      while (h!=NULL)
757      {
758        c=nDiv(pGetCoeff(h),k);
759        if (!nIsOne(c)) nNormalize(c);
760        pSetCoeff(h,c);
761        pIter(h);
762      }
763      nDelete(&k);
764    }
765    else
766    {
767      if (nNormalize != nDummy2)
768      {
769        h = pNext(p1);
770        while (h!=NULL)
771        {
772          nNormalize(pGetCoeff(h));
773          pIter(h);
774        }
775      }
776    }
777  }
778}
779
780/*2
781*normalize all coefficients
782*/
783void p_Normalize(poly p, ring r)
784{
785  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
786  while (p!=NULL)
787  {
788    if (currRing==r) {nTest(pGetCoeff(p));}
789    n_Normalize(pGetCoeff(p),r);
790    pIter(p);
791  }
792}
793
794// splits p into polys with Exp(n) == 0 and Exp(n) != 0
795// Poly with Exp(n) != 0 is reversed
796static void pSplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero)
797{
798  if (p == NULL)
799  {
800    *non_zero = NULL;
801    *zero = NULL;
802    return;
803  }
804  spolyrec sz;
805  poly z, n_z, next;
806  z = &sz;
807  n_z = NULL;
808
809  while(p != NULL)
810  {
811    next = pNext(p);
812    if (pGetExp(p, n) == 0)
813    {
814      pNext(z) = p;
815      pIter(z);
816    }
817    else
818    {
819      pNext(p) = n_z;
820      n_z = p;
821    }
822    p = next;
823  }
824  pNext(z) = NULL;
825  *zero = pNext(&sz);
826  *non_zero = n_z;
827  return;
828}
829
830/*3
831* substitute the n-th variable by 1 in p
832* destroy p
833*/
834static poly pSubst1 (poly p,int n)
835{
836  poly qq=NULL, result = NULL;
837  poly zero=NULL, non_zero=NULL;
838
839  // reverse, so that add is likely to be linear
840  pSplitAndReversePoly(p, n, &non_zero, &zero);
841
842  while (non_zero != NULL)
843  {
844    assume(pGetExp(non_zero, n) != 0);
845    qq = non_zero;
846    pIter(non_zero);
847    qq->next = NULL;
848    pSetExp(qq,n,0);
849    pSetm(qq);
850    result = pAdd(result,qq);
851  }
852  p = pAdd(result, zero);
853  pTest(p);
854  return p;
855}
856
857/*3
858* substitute the n-th variable by number e in p
859* destroy p
860*/
861static poly pSubst2 (poly p,int n, number e)
862{
863  assume( ! nIsZero(e) );
864  poly qq,result = NULL;
865  number nn, nm;
866  poly zero, non_zero;
867
868  // reverse, so that add is likely to be linear
869  pSplitAndReversePoly(p, n, &non_zero, &zero);
870
871  while (non_zero != NULL)
872  {
873    assume(pGetExp(non_zero, n) != 0);
874    qq = non_zero;
875    pIter(non_zero);
876    qq->next = NULL;
877    nPower(e, pGetExp(qq, n), &nn);
878    nm = nMult(nn, pGetCoeff(qq));
879    pSetCoeff(qq, nm);
880    nDelete(&nn);
881    pSetExp(qq, n, 0);
882    pSetm(qq);
883    result = pAdd(result,qq);
884  }
885  p = pAdd(result, zero);
886  pTest(p);
887  return p;
888}
889
890
891/* delete monoms whose n-th exponent is different from zero */
892poly pSubst0(poly p, int n)
893{
894  spolyrec res;
895  poly h = &res;
896  pNext(h) = p;
897
898  while (pNext(h)!=NULL)
899  {
900    if (pGetExp(pNext(h),n)!=0)
901    {
902      pDeleteLm(&pNext(h));
903    }
904    else
905    {
906      pIter(h);
907    }
908  }
909  pTest(pNext(&res));
910  return pNext(&res);
911}
912
913/*2
914* substitute the n-th variable by e in p
915* destroy p
916*/
917poly pSubst(poly p, int n, poly e)
918{
919  if (e == NULL) return pSubst0(p, n);
920
921  if (pIsConstant(e))
922  {
923    if (nIsOne(pGetCoeff(e))) return pSubst1(p,n);
924    else return pSubst2(p, n, pGetCoeff(e));
925  }
926
927#ifdef HAVE_PLURAL
928  if (rIsPluralRing(currRing))
929  {
930    return nc_pSubst(p,n,e);
931  }
932#endif
933
934  int exponent,i;
935  poly h, res, m;
936  int *me,*ee;
937  number nu,nu1;
938
939  me=(int *)omAlloc((pVariables+1)*sizeof(int));
940  ee=(int *)omAlloc((pVariables+1)*sizeof(int));
941  if (e!=NULL) pGetExpV(e,ee);
942  res=NULL;
943  h=p;
944  while (h!=NULL)
945  {
946    if ((e!=NULL) || (pGetExp(h,n)==0))
947    {
948      m=pHead(h);
949      pGetExpV(m,me);
950      exponent=me[n];
951      me[n]=0;
952      for(i=pVariables;i>0;i--)
953        me[i]+=exponent*ee[i];
954      pSetExpV(m,me);
955      if (e!=NULL)
956      {
957        nPower(pGetCoeff(e),exponent,&nu);
958        nu1=nMult(pGetCoeff(m),nu);
959        nDelete(&nu);
960        pSetCoeff(m,nu1);
961      }
962      res=pAdd(res,m);
963    }
964    pDeleteLm(&h);
965  }
966  omFreeSize((ADDRESS)me,(pVariables+1)*sizeof(int));
967  omFreeSize((ADDRESS)ee,(pVariables+1)*sizeof(int));
968  return res;
969}
970
971/* Returns TRUE if
972     * LM(p) | LM(lcm)
973     * LC(p) | LC(lcm) only if ring
974     * Exists i, j:
975         * LE(p, i)  != LE(lcm, i)
976         * LE(p1, i) != LE(lcm, i)   ==> LCM(p1, p) != lcm
977         * LE(p, j)  != LE(lcm, j)
978         * LE(p2, j) != LE(lcm, j)   ==> LCM(p2, p) != lcm
979*/
980BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
981{
982  int k, j;
983
984  if (lcm==NULL) return FALSE;
985
986  for (j=pVariables; j; j--)
987    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
988  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
989  for (j=pVariables; j; j--)
990  {
991    if (pGetExp(p1,j)!=pGetExp(lcm,j))
992    {
993      if (pGetExp(p,j)!=pGetExp(lcm,j))
994      {
995        for (k=pVariables; k>j; k--)
996        {
997          if ((pGetExp(p,k)!=pGetExp(lcm,k))
998          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
999            return TRUE;
1000        }
1001        for (k=j-1; k; k--)
1002        {
1003          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1004          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1005            return TRUE;
1006        }
1007        return FALSE;
1008      }
1009    }
1010    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1011    {
1012      if (pGetExp(p,j)!=pGetExp(lcm,j))
1013      {
1014        for (k=pVariables; k>j; k--)
1015        {
1016          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1017          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1018            return TRUE;
1019        }
1020        for (k=j-1; k!=0 ; k--)
1021        {
1022          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1023          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1024            return TRUE;
1025        }
1026        return FALSE;
1027      }
1028    }
1029  }
1030  return FALSE;
1031}
Note: See TracBrowser for help on using the repository browser.