source: git/kernel/polys.cc @ d51339

spielwiese
Last change on this file since d51339 was 2b87ac9, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: optim git-svn-id: file:///usr/local/Singular/svn/trunk@11309 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 21.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: polys.cc,v 1.35 2009-01-15 10:33:24 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  if (nIsUnit(inv))
177  {
178    inv = nInvers(inv);
179    pMult_nn(result,inv);
180    nDelete(&inv);
181  }
182  else
183  {
184    pDiv_nn(result,inv);
185  }
186#else
187  pMult_nn(result,inv);
188  nDelete(&inv);
189#endif
190  pDelete(&b);
191  return result;
192}
193
194/*2
195* returns the LCM of the head terms of a and b in *m
196*/
197void pLcm(poly a, poly b, poly m)
198{
199  int i;
200  for (i=pVariables; i; i--)
201  {
202    pSetExp(m,i, si_max( pGetExp(a,i), pGetExp(b,i)));
203  }
204  pSetComp(m, si_max(pGetComp(a), pGetComp(b)));
205  /* Don't do a pSetm here, otherwise hres/lres chockes */
206}
207
208/*2
209* convert monomial given as string to poly, e.g. 1x3y5z
210*/
211const char * p_Read(const char *st, poly &rc, const ring r)
212{
213  if (r==NULL) { rc=NULL;return st;}
214  int i,j;
215  rc = p_Init(r);
216  const char *s = r->cf->nRead(st,&(rc->coef));
217  if (s==st)
218  /* i.e. it does not start with a coeff: test if it is a ringvar*/
219  {
220    j = r_IsRingVar(s,r);
221    if (j >= 0)
222    {
223      p_IncrExp(rc,1+j,r);
224      while (*s!='\0') s++;
225      goto done;
226    }
227  }
228  while (*s!='\0')
229  {
230    char ss[2];
231    ss[0] = *s++;
232    ss[1] = '\0';
233    j = r_IsRingVar(ss,r);
234    if (j >= 0)
235    {
236      const char *s_save=s;
237      s = eati(s,&i);
238      if (((unsigned long)i) >  r->bitmask)
239      {
240        // exponent to large: it is not a monomial
241        p_DeleteLm(&rc,r);
242        return s_save;
243      }
244      p_AddExp(rc,1+j, (Exponent_t)i, r);
245    }
246    else
247    {
248      // 1st char of is not a varname
249      p_DeleteLm(&rc,r);
250      s--;
251      return s;
252    }
253  }
254done:
255  if (r->cf->nIsZero(pGetCoeff(rc))) p_DeleteLm(&rc,r);
256  else
257  {
258#ifdef HAVE_PLURAL
259    // in super-commutative ring
260    // squares of anti-commutative variables are zeroes!
261    if(rIsSCA(r))
262    {
263      const unsigned int iFirstAltVar = scaFirstAltVar(r);
264      const unsigned int iLastAltVar  = scaLastAltVar(r);
265
266      assume(rc != NULL);
267
268      for(unsigned int k = iFirstAltVar; k <= iLastAltVar; k++)
269        if( p_GetExp(rc, k, r) > 1 )
270        {
271          p_DeleteLm(&rc, r);
272          goto finish;
273        }
274    }
275#endif
276   
277    p_Setm(rc,r);
278  }
279finish: 
280  return s;
281}
282
283poly pmInit(const char *st, BOOLEAN &ok)
284{
285  poly p;
286  const char *s=p_Read(st,p,currRing);
287  if (*s!='\0')
288  {
289    if ((s!=st)&&isdigit(st[0]))
290    {
291      errorreported=TRUE;
292    }
293    ok=FALSE;
294    pDelete(&p);
295    return NULL;
296  }
297  ok=!errorreported;
298  return p;
299}
300
301/*2
302*make p homogeneous by multiplying the monomials by powers of x_varnum
303*assume: deg(var(varnum))==1
304*/
305poly pHomogen (poly p, int varnum)
306{
307  poly q=NULL, qn;
308  int  o,ii;
309  sBucket_pt bp;
310
311  if (p!=NULL)
312  {
313    if ((varnum < 1) || (varnum > pVariables))
314    {
315      return NULL;
316    }
317    o=pWTotaldegree(p);
318    q=pNext(p);
319    while (q != NULL)
320    {
321      ii=pWTotaldegree(q);
322      if (ii>o) o=ii;
323      pIter(q);
324    }
325    q = pCopy(p);
326    bp = sBucketCreate(currRing);
327    while (q != NULL)
328    {
329      ii = o-pWTotaldegree(q);
330      if (ii!=0)
331      {
332        pAddExp(q,varnum, (Exponent_t)ii);
333        pSetm(q);
334      }
335      qn = pNext(q);
336      pNext(q) = NULL;
337      sBucket_Add_p(bp, q, 1);
338      q = qn;
339    }
340    sBucketDestroyAdd(bp, &q, &ii);
341  }
342  return q;
343}
344
345/*2
346*replaces the maximal powers of the leading monomial of p2 in p1 by
347*the same powers of n, utility for dehomogenization
348*/
349poly pDehomogen (poly p1,poly p2,number n)
350{
351  polyset P;
352  int     SizeOfSet=5;
353  int     i;
354  poly    p;
355  number  nn;
356
357  P = (polyset)omAlloc0(5*sizeof(poly));
358  //for (i=0; i<5; i++)
359  //{
360  //  P[i] = NULL;
361  //}
362  pCancelPolyByMonom(p1,p2,&P,&SizeOfSet);
363  p = P[0];
364  //P[0] = NULL ;// for safety, may be removed later
365  for (i=1; i<SizeOfSet; i++)
366  {
367    if (P[i] != NULL)
368    {
369      nPower(n,i,&nn);
370      pMult_nn(P[i],nn);
371      p = pAdd(p,P[i]);
372      //P[i] =NULL; // for safety, may be removed later
373      nDelete(&nn);
374    }
375  }
376  omFreeSize((ADDRESS)P,SizeOfSet*sizeof(poly));
377  return p;
378}
379
380/*4
381*Returns the exponent of the maximal power of the leading monomial of
382*p2 in that of p1
383*/
384static int pGetMaxPower (poly p1,poly p2)
385{
386  int     i,k,res = 32000; /*a very large integer*/
387
388  if (p1 == NULL) return 0;
389  for (i=1; i<=pVariables; i++)
390  {
391    if ( pGetExp(p2,i) != 0)
392    {
393      k =  pGetExp(p1,i) /  pGetExp(p2,i);
394      if (k < res) res = k;
395    }
396  }
397  return res;
398}
399
400/*2
401*Returns as i-th entry of P the coefficient of the (i-1) power of
402*the leading monomial of p2 in p1
403*/
404void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet)
405{
406  int   maxPow;
407  poly  p,qp,Coeff;
408
409  if (*P == NULL)
410  {
411    *P = (polyset) omAlloc(5*sizeof(poly));
412    *SizeOfSet = 5;
413  }
414  p = pCopy(p1);
415  while (p != NULL)
416  {
417    qp = p->next;
418    p->next = NULL;
419    maxPow = pGetMaxPower(p,p2);
420    Coeff = pDivByMonom(p,p2);
421    if (maxPow > *SizeOfSet)
422    {
423      pEnlargeSet(P,*SizeOfSet,maxPow+1-*SizeOfSet);
424      *SizeOfSet = maxPow+1;
425    }
426    (*P)[maxPow] = pAdd((*P)[maxPow],Coeff);
427    pDelete(&p);
428    p = qp;
429  }
430}
431
432/*2
433*returns the leading monomial of p1 divided by the maximal power of that
434*of p2
435*/
436poly pDivByMonom (poly p1,poly p2)
437{
438  int     k, i;
439
440  if (p1 == NULL) return NULL;
441  k = pGetMaxPower(p1,p2);
442  if (k == 0)
443    return pHead(p1);
444  else
445  {
446    number n;
447    poly p = pInit();
448
449    p->next = NULL;
450    for (i=1; i<=pVariables; i++)
451    {
452       pSetExp(p,i, pGetExp(p1,i)-k* pGetExp(p2,i));
453    }
454    nPower(p2->coef,k,&n);
455    pSetCoeff0(p,nDiv(p1->coef,n));
456    nDelete(&n);
457    pSetm(p);
458    return p;
459  }
460}
461/*----------utilities for syzygies--------------*/
462poly pTakeOutComp(poly * p, int k)
463{
464  poly q = *p,qq=NULL,result = NULL;
465
466  if (q==NULL) return NULL;
467  BOOLEAN use_setmcomp=rOrd_SetCompRequiresSetm(currRing);
468  if (pGetComp(q)==k)
469  {
470    result = q;
471    do
472    {
473      pSetComp(q,0);
474      if (use_setmcomp) pSetmComp(q);
475      qq = q;
476      pIter(q);
477    }
478    while ((q!=NULL) && (pGetComp(q)==k));
479    *p = q;
480    pNext(qq) = NULL;
481  }
482  if (q==NULL) return result;
483  if (pGetComp(q) > k)
484  {
485    pDecrComp(q);
486    if (use_setmcomp) pSetmComp(q);
487  }
488  poly pNext_q;
489  while ((pNext_q=pNext(q))!=NULL)
490  {
491    if (pGetComp(pNext_q)==k)
492    {
493      if (result==NULL)
494      {
495        result = pNext_q;
496        qq = result;
497      }
498      else
499      {
500        pNext(qq) = pNext_q;
501        pIter(qq);
502      }
503      pNext(q) = pNext(pNext_q);
504      pNext(qq) =NULL;
505      pSetComp(qq,0);
506      if (use_setmcomp) pSetmComp(qq);
507    }
508    else
509    {
510      /*pIter(q);*/ q=pNext_q;
511      if (pGetComp(q) > k)
512      {
513        pDecrComp(q);
514        if (use_setmcomp) pSetmComp(q);
515      }
516    }
517  }
518  return result;
519}
520
521// Splits *p into two polys: *q which consists of all monoms with
522// component == comp and *p of all other monoms *lq == pLength(*q)
523void pTakeOutComp(poly *r_p, Exponent_t comp, poly *r_q, int *lq)
524{
525  spolyrec pp, qq;
526  poly p, q, p_prev;
527  int l = 0;
528
529#ifdef HAVE_ASSUME
530  int lp = pLength(*r_p);
531#endif
532
533  pNext(&pp) = *r_p;
534  p = *r_p;
535  p_prev = &pp;
536  q = &qq;
537
538  while(p != NULL)
539  {
540    while (pGetComp(p) == comp)
541    {
542      pNext(q) = p;
543      pIter(q);
544      pSetComp(p, 0);
545      pSetmComp(p);
546      pIter(p);
547      l++;
548      if (p == NULL)
549      {
550        pNext(p_prev) = NULL;
551        goto Finish;
552      }
553    }
554    pNext(p_prev) = p;
555    p_prev = p;
556    pIter(p);
557  }
558
559  Finish:
560  pNext(q) = NULL;
561  *r_p = pNext(&pp);
562  *r_q = pNext(&qq);
563  *lq = l;
564#ifdef HAVE_ASSUME
565  assume(pLength(*r_p) + pLength(*r_q) == lp);
566#endif
567  pTest(*r_p);
568  pTest(*r_q);
569}
570
571void pDecrOrdTakeOutComp(poly *r_p, Exponent_t comp, Order_t order,
572                         poly *r_q, int *lq)
573{
574  spolyrec pp, qq;
575  poly p, q, p_prev;
576  int l = 0;
577
578  pNext(&pp) = *r_p;
579  p = *r_p;
580  p_prev = &pp;
581  q = &qq;
582
583#ifdef HAVE_ASSUME
584  if (p != NULL)
585  {
586    while (pNext(p) != NULL)
587    {
588      assume(pGetOrder(p) >= pGetOrder(pNext(p)));
589      pIter(p);
590    }
591  }
592  p = *r_p;
593#endif
594
595  while (p != NULL && pGetOrder(p) > order) pIter(p);
596
597  while(p != NULL && pGetOrder(p) == order)
598  {
599    while (pGetComp(p) == comp)
600    {
601      pNext(q) = p;
602      pIter(q);
603      pIter(p);
604      pSetComp(p, 0);
605      pSetmComp(p);
606      l++;
607      if (p == NULL || pGetOrder(p) != order)
608      {
609        pNext(p_prev) = p;
610        goto Finish;
611      }
612    }
613    pNext(p_prev) = p;
614    p_prev = p;
615    pIter(p);
616  }
617
618  Finish:
619  pNext(q) = NULL;
620  *r_p = pNext(&pp);
621  *r_q = pNext(&qq);
622  *lq = l;
623}
624
625#if 1
626poly pTakeOutComp1(poly * p, int k)
627{
628  poly q = *p;
629
630  if (q==NULL) return NULL;
631
632  poly qq=NULL,result = NULL;
633
634  if (pGetComp(q)==k)
635  {
636    result = q; /* *p */
637    while ((q!=NULL) && (pGetComp(q)==k))
638    {
639      pSetComp(q,0);
640      pSetmComp(q);
641      qq = q;
642      pIter(q);
643    }
644    *p = q;
645    pNext(qq) = NULL;
646  }
647  if (q==NULL) return result;
648//  if (pGetComp(q) > k) pGetComp(q)--;
649  while (pNext(q)!=NULL)
650  {
651    if (pGetComp(pNext(q))==k)
652    {
653      if (result==NULL)
654      {
655        result = pNext(q);
656        qq = result;
657      }
658      else
659      {
660        pNext(qq) = pNext(q);
661        pIter(qq);
662      }
663      pNext(q) = pNext(pNext(q));
664      pNext(qq) =NULL;
665      pSetComp(qq,0);
666      pSetmComp(qq);
667    }
668    else
669    {
670      pIter(q);
671//      if (pGetComp(q) > k) pGetComp(q)--;
672    }
673  }
674  return result;
675}
676#endif
677
678void pDeleteComp(poly * p,int k)
679{
680  poly q;
681
682  while ((*p!=NULL) && (pGetComp(*p)==k)) pDeleteLm(p);
683  if (*p==NULL) return;
684  q = *p;
685  if (pGetComp(q)>k)
686  {
687    pDecrComp(q);
688    pSetmComp(q);
689  }
690  while (pNext(q)!=NULL)
691  {
692    if (pGetComp(pNext(q))==k)
693      pDeleteLm(&(pNext(q)));
694    else
695    {
696      pIter(q);
697      if (pGetComp(q)>k)
698      {
699        pDecrComp(q);
700        pSetmComp(q);
701      }
702    }
703  }
704}
705/*----------end of utilities for syzygies--------------*/
706
707/*2
708* pair has no common factor ? or is no polynomial
709*/
710BOOLEAN pHasNotCF(poly p1, poly p2)
711{
712
713  if (!TEST_OPT_IDLIFT)
714  {
715    if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
716      return FALSE;
717  }
718  int i = 1;
719  loop
720  {
721    if ((pGetExp(p1, i) > 0) && (pGetExp(p2, i) > 0))   return FALSE;
722    if (i == pVariables)                                return TRUE;
723    i++;
724  }
725}
726
727/*2
728*divides p1 by its leading coefficient
729*/
730void pNorm(poly p1)
731{
732#ifdef HAVE_RINGS
733  if (rField_is_Ring(currRing))
734  {
735    Werror("pNorm not possible in the case of coefficient rings.");
736  }
737  else
738#endif
739  if (p1!=NULL)
740  {
741    if (pNext(p1)==NULL)
742    {
743      pSetCoeff(p1,nInit(1));
744      return;
745    }
746    poly h;
747    if (!nIsOne(pGetCoeff(p1)))
748    {
749      number k, c;
750      nNormalize(pGetCoeff(p1));
751      k = pGetCoeff(p1);
752      c = nInit(1);
753      pSetCoeff0(p1,c);
754      h = pNext(p1);
755      while (h!=NULL)
756      {
757        c=nDiv(pGetCoeff(h),k);
758        // no need to normalize: Z/p, R
759        // normalize already in nDiv: Q_a, Z/p_a
760        // remains: Q
761        if (rField_is_Q() && (!nIsOne(c))) nNormalize(c);
762        pSetCoeff(h,c);
763        pIter(h);
764      }
765      nDelete(&k);
766    }
767    else
768    {
769      if (nNormalize != nDummy2)
770      {
771        h = pNext(p1);
772        while (h!=NULL)
773        {
774          nNormalize(pGetCoeff(h));
775          pIter(h);
776        }
777      }
778    }
779  }
780}
781
782/*2
783*normalize all coefficients
784*/
785void p_Normalize(poly p,const ring r)
786{
787  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
788  while (p!=NULL)
789  {
790    if (currRing==r) {nTest(pGetCoeff(p));}
791    n_Normalize(pGetCoeff(p),r);
792    pIter(p);
793  }
794}
795
796// splits p into polys with Exp(n) == 0 and Exp(n) != 0
797// Poly with Exp(n) != 0 is reversed
798static void pSplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero)
799{
800  if (p == NULL)
801  {
802    *non_zero = NULL;
803    *zero = NULL;
804    return;
805  }
806  spolyrec sz;
807  poly z, n_z, next;
808  z = &sz;
809  n_z = NULL;
810
811  while(p != NULL)
812  {
813    next = pNext(p);
814    if (pGetExp(p, n) == 0)
815    {
816      pNext(z) = p;
817      pIter(z);
818    }
819    else
820    {
821      pNext(p) = n_z;
822      n_z = p;
823    }
824    p = next;
825  }
826  pNext(z) = NULL;
827  *zero = pNext(&sz);
828  *non_zero = n_z;
829  return;
830}
831
832/*3
833* substitute the n-th variable by 1 in p
834* destroy p
835*/
836static poly pSubst1 (poly p,int n)
837{
838  poly qq=NULL, result = NULL;
839  poly zero=NULL, non_zero=NULL;
840
841  // reverse, so that add is likely to be linear
842  pSplitAndReversePoly(p, n, &non_zero, &zero);
843
844  while (non_zero != NULL)
845  {
846    assume(pGetExp(non_zero, n) != 0);
847    qq = non_zero;
848    pIter(non_zero);
849    qq->next = NULL;
850    pSetExp(qq,n,0);
851    pSetm(qq);
852    result = pAdd(result,qq);
853  }
854  p = pAdd(result, zero);
855  pTest(p);
856  return p;
857}
858
859/*3
860* substitute the n-th variable by number e in p
861* destroy p
862*/
863static poly pSubst2 (poly p,int n, number e)
864{
865  assume( ! nIsZero(e) );
866  poly qq,result = NULL;
867  number nn, nm;
868  poly zero, non_zero;
869
870  // reverse, so that add is likely to be linear
871  pSplitAndReversePoly(p, n, &non_zero, &zero);
872
873  while (non_zero != NULL)
874  {
875    assume(pGetExp(non_zero, n) != 0);
876    qq = non_zero;
877    pIter(non_zero);
878    qq->next = NULL;
879    nPower(e, pGetExp(qq, n), &nn);
880    nm = nMult(nn, pGetCoeff(qq));
881#ifdef HAVE_RINGS
882    if (nIsZero(nm))
883    {
884      pLmFree(&qq);
885      nDelete(&nm);
886    }
887    else
888#endif
889    {
890      pSetCoeff(qq, nm);
891      pSetExp(qq, n, 0);
892      pSetm(qq);
893      result = pAdd(result,qq);
894    }
895    nDelete(&nn);
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}
1044BOOLEAN pCompareChainPart (poly p,poly p1,poly p2,poly lcm)
1045{
1046  int k, j;
1047
1048  if (lcm==NULL) return FALSE;
1049
1050  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1051    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
1052  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1053  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1054  {
1055    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1056    {
1057      if (pGetExp(p,j)!=pGetExp(lcm,j))
1058      {
1059        for (k=pVariables; k>j; k--)
1060        for (k=currRing->real_var_end; k>j; k--)
1061        {
1062          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1063          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1064            return TRUE;
1065        }
1066        for (k=j-1; k>=currRing->real_var_start; k--)
1067        {
1068          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1069          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1070            return TRUE;
1071        }
1072        return FALSE;
1073      }
1074    }
1075    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1076    {
1077      if (pGetExp(p,j)!=pGetExp(lcm,j))
1078      {
1079        for (k=currRing->real_var_end; k>j; k--)
1080        {
1081          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1082          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1083            return TRUE;
1084        }
1085        for (k=j-1; k>=currRing->real_var_start; k--)
1086        {
1087          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1088          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1089            return TRUE;
1090        }
1091        return FALSE;
1092      }
1093    }
1094  }
1095  return FALSE;
1096}
1097
1098int pSize(poly p)
1099{
1100  int count = 0;
1101  while ( p != NULL )
1102  {
1103    count+= nSize( pGetCoeff( p ) );
1104    pIter( p );
1105  }
1106  return count;
1107}
1108
Note: See TracBrowser for help on using the repository browser.