source: git/kernel/polys.cc @ 0b5e3d

spielwiese
Last change on this file since 0b5e3d was 0b5e3d, checked in by Hans Schönemann <hannes@…>, 14 years ago
Exponent_t removed git-svn-id: file:///usr/local/Singular/svn/trunk@12664 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 21.8 KB
RevLine 
[35aab3]1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
[341696]4/* $Id$ */
[35aab3]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"
[690e21e]15#include "options.h"
[35aab3]16#include "omalloc.h"
17#include "febase.h"
18#include "numbers.h"
19#include "polys.h"
20#include "ring.h"
21
[68349d]22#ifdef HAVE_PLURAL
23#include "gring.h"
[5d7a3b]24#include "sca.h"
[68349d]25#endif
26
[35aab3]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
[8c5988]53void pSetGlobals(const ring r, BOOLEAN complete)
[35aab3]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
[206e158]113#ifdef HAVE_RINGS   //TODO Oliver
[8251cc]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
[35aab3]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;
[009d80]144#ifdef HAVE_RINGS
145  number inv=pGetCoeff(b);
[8251cc]146#else
[35aab3]147  number inv=nInvers(pGetCoeff(b));
[8251cc]148#endif
[35aab3]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  }
[009d80]175#ifdef HAVE_RINGS
[b547dc]176  if (nIsUnit(inv))
177  {
[ec1744b]178    inv = nInvers(inv);
[b547dc]179    pMult_nn(result,inv);
180    nDelete(&inv);
181  }
182  else
183  {
184    pDiv_nn(result,inv);
185  }
[009d80]186#else
[35aab3]187  pMult_nn(result,inv);
188  nDelete(&inv);
[b1f3b55]189#endif
[35aab3]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*/
[107986]211const char * p_Read(const char *st, poly &rc, const ring r)
[35aab3]212{
213  if (r==NULL) { rc=NULL;return st;}
214  int i,j;
215  rc = p_Init(r);
[85e68dd]216  const char *s = r->cf->nRead(st,&(rc->coef));
[35aab3]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    {
[85e68dd]236      const char *s_save=s;
[35aab3]237      s = eati(s,&i);
[ded4ee]238      if (((unsigned long)i) >  r->bitmask)
239      {
[9b10a76]240        // exponent to large: it is not a monomial
[b529d9]241        p_DeleteLm(&rc,r);
[ded4ee]242        return s_save;
243      }
[0b5e3d]244      p_AddExp(rc,1+j, (long)i, r);
[35aab3]245    }
246    else
247    {
[b529d9]248      // 1st char of is not a varname
249      p_DeleteLm(&rc,r);
[35aab3]250      s--;
251      return s;
252    }
253  }
254done:
255  if (r->cf->nIsZero(pGetCoeff(rc))) p_DeleteLm(&rc,r);
256  else
257  {
[5d7a3b]258#ifdef HAVE_PLURAL
[9b10a76]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        }
[5d7a3b]274    }
275#endif
[9b10a76]276   
[35aab3]277    p_Setm(rc,r);
278  }
[9b10a76]279finish: 
[35aab3]280  return s;
281}
282
[690e21e]283BOOLEAN _p_Test(poly p, ring r, int level);
[85e68dd]284poly pmInit(const char *st, BOOLEAN &ok)
[35aab3]285{
286  poly p;
[85e68dd]287  const char *s=p_Read(st,p,currRing);
[35aab3]288  if (*s!='\0')
289  {
290    if ((s!=st)&&isdigit(st[0]))
291    {
292      errorreported=TRUE;
293    }
294    ok=FALSE;
295    pDelete(&p);
296    return NULL;
297  }
[690e21e]298  #ifdef PDEBUG
299  _p_Test(p,currRing,PDEBUG);
300  #endif
[35aab3]301  ok=!errorreported;
302  return p;
303}
304
305/*2
306*make p homogeneous by multiplying the monomials by powers of x_varnum
[5ef9d3]307*assume: deg(var(varnum))==1
[35aab3]308*/
309poly pHomogen (poly p, int varnum)
310{
[37318d]311  pFDegProc deg;
312  if (pLexOrder && (currRing->order[0]==ringorder_lp))
313    deg=pTotaldegree;
314  else
315    deg=pFDeg;
316
[35aab3]317  poly q=NULL, qn;
318  int  o,ii;
319  sBucket_pt bp;
320
321  if (p!=NULL)
322  {
323    if ((varnum < 1) || (varnum > pVariables))
324    {
325      return NULL;
326    }
[37318d]327    o=deg(p,currRing);
[35aab3]328    q=pNext(p);
329    while (q != NULL)
330    {
[37318d]331      ii=deg(q,currRing);
[35aab3]332      if (ii>o) o=ii;
333      pIter(q);
334    }
335    q = pCopy(p);
336    bp = sBucketCreate(currRing);
337    while (q != NULL)
338    {
[37318d]339      ii = o-deg(q,currRing);
[35aab3]340      if (ii!=0)
341      {
[0b5e3d]342        pAddExp(q,varnum, (long)ii);
[35aab3]343        pSetm(q);
344      }
345      qn = pNext(q);
346      pNext(q) = NULL;
347      sBucket_Add_p(bp, q, 1);
348      q = qn;
349    }
350    sBucketDestroyAdd(bp, &q, &ii);
351  }
352  return q;
353}
354
355/*2
356*replaces the maximal powers of the leading monomial of p2 in p1 by
357*the same powers of n, utility for dehomogenization
358*/
359poly pDehomogen (poly p1,poly p2,number n)
360{
361  polyset P;
362  int     SizeOfSet=5;
363  int     i;
364  poly    p;
365  number  nn;
366
367  P = (polyset)omAlloc0(5*sizeof(poly));
368  //for (i=0; i<5; i++)
369  //{
370  //  P[i] = NULL;
371  //}
372  pCancelPolyByMonom(p1,p2,&P,&SizeOfSet);
373  p = P[0];
374  //P[0] = NULL ;// for safety, may be removed later
375  for (i=1; i<SizeOfSet; i++)
376  {
377    if (P[i] != NULL)
378    {
379      nPower(n,i,&nn);
380      pMult_nn(P[i],nn);
381      p = pAdd(p,P[i]);
382      //P[i] =NULL; // for safety, may be removed later
383      nDelete(&nn);
384    }
385  }
386  omFreeSize((ADDRESS)P,SizeOfSet*sizeof(poly));
387  return p;
388}
389
390/*4
391*Returns the exponent of the maximal power of the leading monomial of
392*p2 in that of p1
393*/
394static int pGetMaxPower (poly p1,poly p2)
395{
396  int     i,k,res = 32000; /*a very large integer*/
397
398  if (p1 == NULL) return 0;
399  for (i=1; i<=pVariables; i++)
400  {
401    if ( pGetExp(p2,i) != 0)
402    {
403      k =  pGetExp(p1,i) /  pGetExp(p2,i);
404      if (k < res) res = k;
405    }
406  }
407  return res;
408}
409
410/*2
411*Returns as i-th entry of P the coefficient of the (i-1) power of
412*the leading monomial of p2 in p1
413*/
414void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet)
415{
416  int   maxPow;
417  poly  p,qp,Coeff;
418
419  if (*P == NULL)
420  {
421    *P = (polyset) omAlloc(5*sizeof(poly));
422    *SizeOfSet = 5;
423  }
424  p = pCopy(p1);
425  while (p != NULL)
426  {
427    qp = p->next;
428    p->next = NULL;
429    maxPow = pGetMaxPower(p,p2);
430    Coeff = pDivByMonom(p,p2);
431    if (maxPow > *SizeOfSet)
432    {
433      pEnlargeSet(P,*SizeOfSet,maxPow+1-*SizeOfSet);
434      *SizeOfSet = maxPow+1;
435    }
436    (*P)[maxPow] = pAdd((*P)[maxPow],Coeff);
437    pDelete(&p);
438    p = qp;
439  }
440}
441
442/*2
443*returns the leading monomial of p1 divided by the maximal power of that
444*of p2
445*/
446poly pDivByMonom (poly p1,poly p2)
447{
448  int     k, i;
449
450  if (p1 == NULL) return NULL;
451  k = pGetMaxPower(p1,p2);
452  if (k == 0)
453    return pHead(p1);
454  else
455  {
456    number n;
457    poly p = pInit();
458
459    p->next = NULL;
460    for (i=1; i<=pVariables; i++)
461    {
462       pSetExp(p,i, pGetExp(p1,i)-k* pGetExp(p2,i));
463    }
464    nPower(p2->coef,k,&n);
465    pSetCoeff0(p,nDiv(p1->coef,n));
466    nDelete(&n);
467    pSetm(p);
468    return p;
469  }
470}
471/*----------utilities for syzygies--------------*/
472poly pTakeOutComp(poly * p, int k)
473{
474  poly q = *p,qq=NULL,result = NULL;
475
476  if (q==NULL) return NULL;
[2b87ac9]477  BOOLEAN use_setmcomp=rOrd_SetCompRequiresSetm(currRing);
[35aab3]478  if (pGetComp(q)==k)
479  {
480    result = q;
[0ef842]481    do
[35aab3]482    {
483      pSetComp(q,0);
[2b87ac9]484      if (use_setmcomp) pSetmComp(q);
[35aab3]485      qq = q;
486      pIter(q);
487    }
[0ef842]488    while ((q!=NULL) && (pGetComp(q)==k));
[35aab3]489    *p = q;
490    pNext(qq) = NULL;
491  }
492  if (q==NULL) return result;
493  if (pGetComp(q) > k)
494  {
495    pDecrComp(q);
[2b87ac9]496    if (use_setmcomp) pSetmComp(q);
[35aab3]497  }
498  poly pNext_q;
499  while ((pNext_q=pNext(q))!=NULL)
500  {
501    if (pGetComp(pNext_q)==k)
502    {
503      if (result==NULL)
504      {
505        result = pNext_q;
506        qq = result;
507      }
508      else
509      {
510        pNext(qq) = pNext_q;
511        pIter(qq);
512      }
513      pNext(q) = pNext(pNext_q);
514      pNext(qq) =NULL;
515      pSetComp(qq,0);
[2b87ac9]516      if (use_setmcomp) pSetmComp(qq);
[35aab3]517    }
518    else
519    {
520      /*pIter(q);*/ q=pNext_q;
521      if (pGetComp(q) > k)
522      {
523        pDecrComp(q);
[2b87ac9]524        if (use_setmcomp) pSetmComp(q);
[35aab3]525      }
526    }
527  }
528  return result;
529}
530
531// Splits *p into two polys: *q which consists of all monoms with
532// component == comp and *p of all other monoms *lq == pLength(*q)
[0b5e3d]533void pTakeOutComp(poly *r_p, long comp, poly *r_q, int *lq)
[35aab3]534{
535  spolyrec pp, qq;
536  poly p, q, p_prev;
537  int l = 0;
538
539#ifdef HAVE_ASSUME
540  int lp = pLength(*r_p);
541#endif
542
543  pNext(&pp) = *r_p;
544  p = *r_p;
545  p_prev = &pp;
546  q = &qq;
547
548  while(p != NULL)
549  {
550    while (pGetComp(p) == comp)
551    {
552      pNext(q) = p;
553      pIter(q);
554      pSetComp(p, 0);
555      pSetmComp(p);
556      pIter(p);
557      l++;
558      if (p == NULL)
559      {
560        pNext(p_prev) = NULL;
561        goto Finish;
562      }
563    }
564    pNext(p_prev) = p;
565    p_prev = p;
566    pIter(p);
567  }
568
569  Finish:
570  pNext(q) = NULL;
571  *r_p = pNext(&pp);
572  *r_q = pNext(&qq);
573  *lq = l;
574#ifdef HAVE_ASSUME
575  assume(pLength(*r_p) + pLength(*r_q) == lp);
576#endif
577  pTest(*r_p);
578  pTest(*r_q);
579}
580
[0b5e3d]581void pDecrOrdTakeOutComp(poly *r_p, long comp, long order,
[35aab3]582                         poly *r_q, int *lq)
583{
584  spolyrec pp, qq;
585  poly p, q, p_prev;
586  int l = 0;
587
588  pNext(&pp) = *r_p;
589  p = *r_p;
590  p_prev = &pp;
591  q = &qq;
592
593#ifdef HAVE_ASSUME
594  if (p != NULL)
595  {
596    while (pNext(p) != NULL)
597    {
598      assume(pGetOrder(p) >= pGetOrder(pNext(p)));
599      pIter(p);
600    }
601  }
602  p = *r_p;
603#endif
604
605  while (p != NULL && pGetOrder(p) > order) pIter(p);
606
607  while(p != NULL && pGetOrder(p) == order)
608  {
609    while (pGetComp(p) == comp)
610    {
611      pNext(q) = p;
612      pIter(q);
613      pIter(p);
614      pSetComp(p, 0);
615      pSetmComp(p);
616      l++;
617      if (p == NULL || pGetOrder(p) != order)
618      {
619        pNext(p_prev) = p;
620        goto Finish;
621      }
622    }
623    pNext(p_prev) = p;
624    p_prev = p;
625    pIter(p);
626  }
627
628  Finish:
629  pNext(q) = NULL;
630  *r_p = pNext(&pp);
631  *r_q = pNext(&qq);
632  *lq = l;
633}
634
635#if 1
636poly pTakeOutComp1(poly * p, int k)
637{
638  poly q = *p;
639
640  if (q==NULL) return NULL;
641
642  poly qq=NULL,result = NULL;
643
644  if (pGetComp(q)==k)
645  {
646    result = q; /* *p */
647    while ((q!=NULL) && (pGetComp(q)==k))
648    {
649      pSetComp(q,0);
650      pSetmComp(q);
651      qq = q;
652      pIter(q);
653    }
654    *p = q;
655    pNext(qq) = NULL;
656  }
657  if (q==NULL) return result;
658//  if (pGetComp(q) > k) pGetComp(q)--;
659  while (pNext(q)!=NULL)
660  {
661    if (pGetComp(pNext(q))==k)
662    {
663      if (result==NULL)
664      {
665        result = pNext(q);
666        qq = result;
667      }
668      else
669      {
670        pNext(qq) = pNext(q);
671        pIter(qq);
672      }
673      pNext(q) = pNext(pNext(q));
674      pNext(qq) =NULL;
675      pSetComp(qq,0);
676      pSetmComp(qq);
677    }
678    else
679    {
680      pIter(q);
681//      if (pGetComp(q) > k) pGetComp(q)--;
682    }
683  }
684  return result;
685}
686#endif
687
688void pDeleteComp(poly * p,int k)
689{
690  poly q;
691
692  while ((*p!=NULL) && (pGetComp(*p)==k)) pDeleteLm(p);
693  if (*p==NULL) return;
694  q = *p;
695  if (pGetComp(q)>k)
696  {
697    pDecrComp(q);
698    pSetmComp(q);
699  }
700  while (pNext(q)!=NULL)
701  {
702    if (pGetComp(pNext(q))==k)
703      pDeleteLm(&(pNext(q)));
704    else
705    {
706      pIter(q);
707      if (pGetComp(q)>k)
708      {
709        pDecrComp(q);
710        pSetmComp(q);
711      }
712    }
713  }
714}
715/*----------end of utilities for syzygies--------------*/
716
717/*2
718* pair has no common factor ? or is no polynomial
719*/
720BOOLEAN pHasNotCF(poly p1, poly p2)
721{
722
[1c35568]723  if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
724    return FALSE;
[d2fb9b7]725  int i = pVariables;
[35aab3]726  loop
727  {
728    if ((pGetExp(p1, i) > 0) && (pGetExp(p2, i) > 0))   return FALSE;
[d2fb9b7]729    i--;
730    if (i == 0)                                         return TRUE;
[35aab3]731  }
732}
733
734/*2
735*divides p1 by its leading coefficient
736*/
737void pNorm(poly p1)
738{
[009d80]739#ifdef HAVE_RINGS
740  if (rField_is_Ring(currRing))
[585bbcb]741  {
[521349]742    Werror("pNorm not possible in the case of coefficient rings.");
[585bbcb]743  }
[994445]744  else
[f92547]745#endif
[35aab3]746  if (p1!=NULL)
747  {
[5e8fe91]748    if (pNext(p1)==NULL)
749    {
[98938c]750      pSetCoeff(p1,nInit(1));
[5e8fe91]751      return;
752    }
[521349]753    poly h;
[35aab3]754    if (!nIsOne(pGetCoeff(p1)))
755    {
[521349]756      number k, c;
[35aab3]757      nNormalize(pGetCoeff(p1));
[585bbcb]758      k = pGetCoeff(p1);
[35aab3]759      c = nInit(1);
760      pSetCoeff0(p1,c);
761      h = pNext(p1);
762      while (h!=NULL)
763      {
764        c=nDiv(pGetCoeff(h),k);
[7ba059]765        // no need to normalize: Z/p, R
766        // normalize already in nDiv: Q_a, Z/p_a
767        // remains: Q
768        if (rField_is_Q() && (!nIsOne(c))) nNormalize(c);
[35aab3]769        pSetCoeff(h,c);
770        pIter(h);
771      }
772      nDelete(&k);
773    }
774    else
775    {
776      if (nNormalize != nDummy2)
777      {
778        h = pNext(p1);
779        while (h!=NULL)
780        {
781          nNormalize(pGetCoeff(h));
782          pIter(h);
783        }
784      }
785    }
786  }
787}
788
789/*2
790*normalize all coefficients
791*/
[107986]792void p_Normalize(poly p,const ring r)
[35aab3]793{
794  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
795  while (p!=NULL)
796  {
[47b2b2d]797#ifdef LDEBUG
[35aab3]798    if (currRing==r) {nTest(pGetCoeff(p));}
[47b2b2d]799#endif
[35aab3]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{
[788529d]847  poly qq=NULL, result = NULL;
848  poly zero=NULL, non_zero=NULL;
[35aab3]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));
[67dbdb]890#ifdef HAVE_RINGS
891    if (nIsZero(nm))
892    {
893      pLmFree(&qq);
894      nDelete(&nm);
895    }
896    else
897#endif
898    {
899      pSetCoeff(qq, nm);
900      pSetExp(qq, n, 0);
901      pSetm(qq);
902      result = pAdd(result,qq);
903    }
[35aab3]904    nDelete(&nn);
905  }
906  p = pAdd(result, zero);
907  pTest(p);
908  return p;
909}
910
911
912/* delete monoms whose n-th exponent is different from zero */
913poly pSubst0(poly p, int n)
914{
915  spolyrec res;
916  poly h = &res;
917  pNext(h) = p;
918
919  while (pNext(h)!=NULL)
920  {
921    if (pGetExp(pNext(h),n)!=0)
922    {
923      pDeleteLm(&pNext(h));
924    }
925    else
926    {
927      pIter(h);
928    }
929  }
930  pTest(pNext(&res));
931  return pNext(&res);
932}
933
934/*2
935* substitute the n-th variable by e in p
936* destroy p
937*/
938poly pSubst(poly p, int n, poly e)
939{
940  if (e == NULL) return pSubst0(p, n);
941
942  if (pIsConstant(e))
943  {
944    if (nIsOne(pGetCoeff(e))) return pSubst1(p,n);
945    else return pSubst2(p, n, pGetCoeff(e));
946  }
947
[68349d]948#ifdef HAVE_PLURAL
949  if (rIsPluralRing(currRing))
950  {
951    return nc_pSubst(p,n,e);
952  }
953#endif
954
[35aab3]955  int exponent,i;
956  poly h, res, m;
957  int *me,*ee;
958  number nu,nu1;
959
960  me=(int *)omAlloc((pVariables+1)*sizeof(int));
961  ee=(int *)omAlloc((pVariables+1)*sizeof(int));
962  if (e!=NULL) pGetExpV(e,ee);
963  res=NULL;
964  h=p;
965  while (h!=NULL)
966  {
967    if ((e!=NULL) || (pGetExp(h,n)==0))
968    {
969      m=pHead(h);
970      pGetExpV(m,me);
971      exponent=me[n];
972      me[n]=0;
973      for(i=pVariables;i>0;i--)
974        me[i]+=exponent*ee[i];
975      pSetExpV(m,me);
976      if (e!=NULL)
977      {
978        nPower(pGetCoeff(e),exponent,&nu);
979        nu1=nMult(pGetCoeff(m),nu);
980        nDelete(&nu);
981        pSetCoeff(m,nu1);
982      }
983      res=pAdd(res,m);
984    }
985    pDeleteLm(&h);
986  }
987  omFreeSize((ADDRESS)me,(pVariables+1)*sizeof(int));
988  omFreeSize((ADDRESS)ee,(pVariables+1)*sizeof(int));
989  return res;
990}
991
[a2466f]992/* Returns TRUE if
993     * LM(p) | LM(lcm)
994     * LC(p) | LC(lcm) only if ring
995     * Exists i, j:
996         * LE(p, i)  != LE(lcm, i)
997         * LE(p1, i) != LE(lcm, i)   ==> LCM(p1, p) != lcm
998         * LE(p, j)  != LE(lcm, j)
999         * LE(p2, j) != LE(lcm, j)   ==> LCM(p2, p) != lcm
1000*/
[35aab3]1001BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
1002{
1003  int k, j;
1004
1005  if (lcm==NULL) return FALSE;
1006
1007  for (j=pVariables; j; j--)
1008    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
1009  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1010  for (j=pVariables; j; j--)
1011  {
1012    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1013    {
1014      if (pGetExp(p,j)!=pGetExp(lcm,j))
1015      {
1016        for (k=pVariables; k>j; k--)
1017        {
1018          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1019          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1020            return TRUE;
1021        }
1022        for (k=j-1; k; k--)
1023        {
1024          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1025          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1026            return TRUE;
1027        }
1028        return FALSE;
1029      }
1030    }
1031    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1032    {
1033      if (pGetExp(p,j)!=pGetExp(lcm,j))
1034      {
1035        for (k=pVariables; k>j; k--)
1036        {
1037          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1038          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1039            return TRUE;
1040        }
1041        for (k=j-1; k!=0 ; k--)
1042        {
1043          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1044          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1045            return TRUE;
1046        }
1047        return FALSE;
1048      }
1049    }
1050  }
1051  return FALSE;
1052}
[27e750]1053#ifdef HAVE_RATGRING
[107986]1054BOOLEAN pCompareChainPart (poly p,poly p1,poly p2,poly lcm)
1055{
1056  int k, j;
1057
1058  if (lcm==NULL) return FALSE;
1059
1060  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1061    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
1062  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
1063  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
1064  {
1065    if (pGetExp(p1,j)!=pGetExp(lcm,j))
1066    {
1067      if (pGetExp(p,j)!=pGetExp(lcm,j))
1068      {
1069        for (k=pVariables; k>j; k--)
1070        for (k=currRing->real_var_end; k>j; k--)
1071        {
1072          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1073          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1074            return TRUE;
1075        }
1076        for (k=j-1; k>=currRing->real_var_start; k--)
1077        {
1078          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1079          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
1080            return TRUE;
1081        }
1082        return FALSE;
1083      }
1084    }
1085    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
1086    {
1087      if (pGetExp(p,j)!=pGetExp(lcm,j))
1088      {
1089        for (k=currRing->real_var_end; k>j; k--)
1090        {
1091          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1092          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1093            return TRUE;
1094        }
1095        for (k=j-1; k>=currRing->real_var_start; k--)
1096        {
1097          if ((pGetExp(p,k)!=pGetExp(lcm,k))
1098          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
1099            return TRUE;
1100        }
1101        return FALSE;
1102      }
1103    }
1104  }
1105  return FALSE;
1106}
[27e750]1107#endif
[7ba059]1108
1109int pSize(poly p)
1110{
1111  int count = 0;
1112  while ( p != NULL )
1113  {
1114    count+= nSize( pGetCoeff( p ) );
1115    pIter( p );
1116  }
1117  return count;
1118}
1119
[48884f2]1120/*2
1121* returns the length of a (numbers of monomials)
1122* respect syzComp
1123*/
1124poly pLast(poly a, int &l)
1125{
1126  if (a == NULL)
1127  {
1128    l = 0;
1129    return NULL;
1130  }
1131  l = 1;
1132  if (! rIsSyzIndexRing(currRing))
1133  {
1134    while (pNext(a)!=NULL)
1135    {
1136      pIter(a);
1137      l++;
1138    }
1139  }
1140  else
1141  {
1142    int curr_limit = rGetCurrSyzLimit(currRing);
1143    poly pp = a;
1144    while ((a=pNext(a))!=NULL)
1145    {
1146      if (pGetComp(a)<=curr_limit/*syzComp*/)
1147        l++;
1148      else break;
1149      pp = a;
1150    }
1151    a=pp;
1152  }
1153  return a;
1154}
1155
Note: See TracBrowser for help on using the repository browser.