source: git/libpolys/polys/polys.cc @ 014b65

spielwiese
Last change on this file since 014b65 was 014b65, checked in by Mohamed Barakat <mohamed.barakat@…>, 13 years ago
- moved misc,reporter,resources,coeffs,polys -> (new) libpolys (Hans agreed) - migrated to automake in coeffs, misc status: everything builds (except polys) todo: . migrate resources and reporter to automake . create autoconf macros for omalloc, factory, and libpolys
  • Property mode set to 100644
File size: 13.6 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
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 <auxiliary.h>
15#include "options.h"
16#include "omalloc.h"
17#include "reporter.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/* 1 for polynomial ring, -1 otherwise */
28int     pOrdSgn;
29// it is of type int, not BOOLEAN because it is also in ip
30/* TRUE if the monomial ordering is not compatible with pFDeg */
31BOOLEAN pLexOrder;
32
33/* ----------- global variables, set by procedures from hecke/kstd1 ----- */
34/* the highest monomial below pHEdge */
35poly      ppNoether = NULL;
36
37#ifdef HAVE_RINGS   //TODO Oliver
38#define pDiv_nn(p, n)              p_Div_nn(p, n, currRing)
39
40poly p_Div_nn(poly p, const number n, const ring r)
41{
42  pAssume(!n_IsZero(n,r));
43  p_Test(p, r);
44
45  poly q = p;
46  while (p != NULL)
47  {
48    number nc = pGetCoeff(p);
49    pSetCoeff0(p, n_Div(nc, n, r));
50    n_Delete(&nc, r);
51    pIter(p);
52  }
53  p_Test(q, r);
54  return q;
55}
56#endif
57
58#ifdef HAVE_RINGS
59/* TRUE iff LT(f) | LT(g) */
60BOOLEAN pDivisibleByRingCase(poly f, poly g)
61{
62  int exponent;
63  for(int i = (int)pVariables; i; i--)
64  {
65    exponent = pGetExp(g, i) - pGetExp(f, i);
66    if (exponent < 0) return FALSE;
67  }
68  return nDivBy(pGetCoeff(g), pGetCoeff(f));
69}
70#endif
71
72BOOLEAN _p_Test(poly p, ring r, int level);
73
74/*2
75*make p homogeneous by multiplying the monomials by powers of x_varnum
76*assume: deg(var(varnum))==1
77*/
78poly pHomogen (poly p, int varnum)
79{
80  pFDegProc deg;
81  if (pLexOrder && (currRing->order[0]==ringorder_lp))
82    deg=p_Totaldegree;
83  else
84    deg=pFDeg;
85
86  poly q=NULL, qn;
87  int  o,ii;
88  sBucket_pt bp;
89
90  if (p!=NULL)
91  {
92    if ((varnum < 1) || (varnum > pVariables))
93    {
94      return NULL;
95    }
96    o=deg(p,currRing);
97    q=pNext(p);
98    while (q != NULL)
99    {
100      ii=deg(q,currRing);
101      if (ii>o) o=ii;
102      pIter(q);
103    }
104    q = pCopy(p);
105    bp = sBucketCreate(currRing);
106    while (q != NULL)
107    {
108      ii = o-deg(q,currRing);
109      if (ii!=0)
110      {
111        pAddExp(q,varnum, (long)ii);
112        pSetm(q);
113      }
114      qn = pNext(q);
115      pNext(q) = NULL;
116      sBucket_Add_p(bp, q, 1);
117      q = qn;
118    }
119    sBucketDestroyAdd(bp, &q, &ii);
120  }
121  return q;
122}
123
124/*----------utilities for syzygies--------------*/
125poly pTakeOutComp(poly * p, int k)
126{
127  poly q = *p,qq=NULL,result = NULL;
128
129  if (q==NULL) return NULL;
130  BOOLEAN use_setmcomp=rOrd_SetCompRequiresSetm(currRing);
131  if (pGetComp(q)==k)
132  {
133    result = q;
134    do
135    {
136      pSetComp(q,0);
137      if (use_setmcomp) pSetmComp(q);
138      qq = q;
139      pIter(q);
140    }
141    while ((q!=NULL) && (pGetComp(q)==k));
142    *p = q;
143    pNext(qq) = NULL;
144  }
145  if (q==NULL) return result;
146  if (pGetComp(q) > k)
147  {
148    pSubComp(q,1);
149    if (use_setmcomp) pSetmComp(q);
150  }
151  poly pNext_q;
152  while ((pNext_q=pNext(q))!=NULL)
153  {
154    if (pGetComp(pNext_q)==k)
155    {
156      if (result==NULL)
157      {
158        result = pNext_q;
159        qq = result;
160      }
161      else
162      {
163        pNext(qq) = pNext_q;
164        pIter(qq);
165      }
166      pNext(q) = pNext(pNext_q);
167      pNext(qq) =NULL;
168      pSetComp(qq,0);
169      if (use_setmcomp) pSetmComp(qq);
170    }
171    else
172    {
173      /*pIter(q);*/ q=pNext_q;
174      if (pGetComp(q) > k)
175      {
176        pSubComp(q,1);
177        if (use_setmcomp) pSetmComp(q);
178      }
179    }
180  }
181  return result;
182}
183
184// Splits *p into two polys: *q which consists of all monoms with
185// component == comp and *p of all other monoms *lq == pLength(*q)
186void pTakeOutComp(poly *r_p, long comp, poly *r_q, int *lq)
187{
188  spolyrec pp, qq;
189  poly p, q, p_prev;
190  int l = 0;
191
192#ifdef HAVE_ASSUME
193  int lp = pLength(*r_p);
194#endif
195
196  pNext(&pp) = *r_p;
197  p = *r_p;
198  p_prev = &pp;
199  q = &qq;
200
201  while(p != NULL)
202  {
203    while (pGetComp(p) == comp)
204    {
205      pNext(q) = p;
206      pIter(q);
207      pSetComp(p, 0);
208      pSetmComp(p);
209      pIter(p);
210      l++;
211      if (p == NULL)
212      {
213        pNext(p_prev) = NULL;
214        goto Finish;
215      }
216    }
217    pNext(p_prev) = p;
218    p_prev = p;
219    pIter(p);
220  }
221
222  Finish:
223  pNext(q) = NULL;
224  *r_p = pNext(&pp);
225  *r_q = pNext(&qq);
226  *lq = l;
227#ifdef HAVE_ASSUME
228  assume(pLength(*r_p) + pLength(*r_q) == lp);
229#endif
230  pTest(*r_p);
231  pTest(*r_q);
232}
233
234#if 1
235poly pTakeOutComp1(poly * p, int k)
236{
237  poly q = *p;
238
239  if (q==NULL) return NULL;
240
241  poly qq=NULL,result = NULL;
242
243  if (pGetComp(q)==k)
244  {
245    result = q; /* *p */
246    while ((q!=NULL) && (pGetComp(q)==k))
247    {
248      pSetComp(q,0);
249      pSetmComp(q);
250      qq = q;
251      pIter(q);
252    }
253    *p = q;
254    pNext(qq) = NULL;
255  }
256  if (q==NULL) return result;
257//  if (pGetComp(q) > k) pGetComp(q)--;
258  while (pNext(q)!=NULL)
259  {
260    if (pGetComp(pNext(q))==k)
261    {
262      if (result==NULL)
263      {
264        result = pNext(q);
265        qq = result;
266      }
267      else
268      {
269        pNext(qq) = pNext(q);
270        pIter(qq);
271      }
272      pNext(q) = pNext(pNext(q));
273      pNext(qq) =NULL;
274      pSetComp(qq,0);
275      pSetmComp(qq);
276    }
277    else
278    {
279      pIter(q);
280//      if (pGetComp(q) > k) pGetComp(q)--;
281    }
282  }
283  return result;
284}
285#endif
286
287void pDeleteComp(poly * p,int k)
288{
289  poly q;
290
291  while ((*p!=NULL) && (pGetComp(*p)==k)) pLmDelete(p);
292  if (*p==NULL) return;
293  q = *p;
294  if (pGetComp(q)>k)
295  {
296    pSubComp(q,1);
297    pSetmComp(q);
298  }
299  while (pNext(q)!=NULL)
300  {
301    if (pGetComp(pNext(q))==k)
302      pLmDelete(&(pNext(q)));
303    else
304    {
305      pIter(q);
306      if (pGetComp(q)>k)
307      {
308        pSubComp(q,1);
309        pSetmComp(q);
310      }
311    }
312  }
313}
314/*----------end of utilities for syzygies--------------*/
315
316/*2
317* divides p1 by its leading coefficient if it is a unit
318* (this will always be true over fields; but not over coefficient rings)
319*/
320void pNorm(poly p1)
321{
322#ifdef HAVE_RINGS
323  if (rField_is_Ring(currRing))
324  {
325    if (!nIsUnit(pGetCoeff(p1))) return;
326  }
327#endif
328  if (p1!=NULL)
329  {
330    if (pNext(p1)==NULL)
331    {
332      pSetCoeff(p1,nInit(1));
333      return;
334    }
335    poly h;
336    if (!nIsOne(pGetCoeff(p1)))
337    {
338      number k, c;
339      nNormalize(pGetCoeff(p1));
340      k = pGetCoeff(p1);
341      c = nInit(1);
342      pSetCoeff0(p1,c);
343      h = pNext(p1);
344      while (h!=NULL)
345      {
346        c=nDiv(pGetCoeff(h),k);
347        // no need to normalize: Z/p, R
348        // normalize already in nDiv: Q_a, Z/p_a
349        // remains: Q
350        if (rField_is_Q() && (!nIsOne(c))) nNormalize(c);
351        pSetCoeff(h,c);
352        pIter(h);
353      }
354      nDelete(&k);
355    }
356    else
357    {
358      if (nNormalize != nDummy2)
359      {
360        h = pNext(p1);
361        while (h!=NULL)
362        {
363          nNormalize(pGetCoeff(h));
364          pIter(h);
365        }
366      }
367    }
368  }
369}
370
371/*2
372*normalize all coefficients
373*/
374void p_Normalize(poly p,const ring r)
375{
376  if (rField_has_simple_inverse(r)) return; /* Z/p, GF(p,n), R, long R/C */
377  while (p!=NULL)
378  {
379#ifdef LDEBUG
380    if (currRing==r) {nTest(pGetCoeff(p));}
381#endif
382    n_Normalize(pGetCoeff(p),r);
383    pIter(p);
384  }
385}
386
387// splits p into polys with Exp(n) == 0 and Exp(n) != 0
388// Poly with Exp(n) != 0 is reversed
389static void pSplitAndReversePoly(poly p, int n, poly *non_zero, poly *zero)
390{
391  if (p == NULL)
392  {
393    *non_zero = NULL;
394    *zero = NULL;
395    return;
396  }
397  spolyrec sz;
398  poly z, n_z, next;
399  z = &sz;
400  n_z = NULL;
401
402  while(p != NULL)
403  {
404    next = pNext(p);
405    if (pGetExp(p, n) == 0)
406    {
407      pNext(z) = p;
408      pIter(z);
409    }
410    else
411    {
412      pNext(p) = n_z;
413      n_z = p;
414    }
415    p = next;
416  }
417  pNext(z) = NULL;
418  *zero = pNext(&sz);
419  *non_zero = n_z;
420  return;
421}
422
423/*3
424* substitute the n-th variable by 1 in p
425* destroy p
426*/
427static poly pSubst1 (poly p,int n)
428{
429  poly qq=NULL, result = NULL;
430  poly zero=NULL, non_zero=NULL;
431
432  // reverse, so that add is likely to be linear
433  pSplitAndReversePoly(p, n, &non_zero, &zero);
434
435  while (non_zero != NULL)
436  {
437    assume(pGetExp(non_zero, n) != 0);
438    qq = non_zero;
439    pIter(non_zero);
440    qq->next = NULL;
441    pSetExp(qq,n,0);
442    pSetm(qq);
443    result = pAdd(result,qq);
444  }
445  p = pAdd(result, zero);
446  pTest(p);
447  return p;
448}
449
450/*3
451* substitute the n-th variable by number e in p
452* destroy p
453*/
454static poly pSubst2 (poly p,int n, number e)
455{
456  assume( ! nIsZero(e) );
457  poly qq,result = NULL;
458  number nn, nm;
459  poly zero, non_zero;
460
461  // reverse, so that add is likely to be linear
462  pSplitAndReversePoly(p, n, &non_zero, &zero);
463
464  while (non_zero != NULL)
465  {
466    assume(pGetExp(non_zero, n) != 0);
467    qq = non_zero;
468    pIter(non_zero);
469    qq->next = NULL;
470    nPower(e, pGetExp(qq, n), &nn);
471    nm = nMult(nn, pGetCoeff(qq));
472#ifdef HAVE_RINGS
473    if (nIsZero(nm))
474    {
475      pLmFree(&qq);
476      nDelete(&nm);
477    }
478    else
479#endif
480    {
481      pSetCoeff(qq, nm);
482      pSetExp(qq, n, 0);
483      pSetm(qq);
484      result = pAdd(result,qq);
485    }
486    nDelete(&nn);
487  }
488  p = pAdd(result, zero);
489  pTest(p);
490  return p;
491}
492
493
494/* delete monoms whose n-th exponent is different from zero */
495poly pSubst0(poly p, int n)
496{
497  spolyrec res;
498  poly h = &res;
499  pNext(h) = p;
500
501  while (pNext(h)!=NULL)
502  {
503    if (pGetExp(pNext(h),n)!=0)
504    {
505      pLmDelete(&pNext(h));
506    }
507    else
508    {
509      pIter(h);
510    }
511  }
512  pTest(pNext(&res));
513  return pNext(&res);
514}
515
516/*2
517* substitute the n-th variable by e in p
518* destroy p
519*/
520poly pSubst(poly p, int n, poly e)
521{
522  if (e == NULL) return pSubst0(p, n);
523
524  if (pIsConstant(e))
525  {
526    if (nIsOne(pGetCoeff(e))) return pSubst1(p,n);
527    else return pSubst2(p, n, pGetCoeff(e));
528  }
529
530#ifdef HAVE_PLURAL
531  if (rIsPluralRing(currRing))
532  {
533    return nc_pSubst(p,n,e);
534  }
535#endif
536
537  int exponent,i;
538  poly h, res, m;
539  int *me,*ee;
540  number nu,nu1;
541
542  me=(int *)omAlloc((pVariables+1)*sizeof(int));
543  ee=(int *)omAlloc((pVariables+1)*sizeof(int));
544  if (e!=NULL) pGetExpV(e,ee);
545  res=NULL;
546  h=p;
547  while (h!=NULL)
548  {
549    if ((e!=NULL) || (pGetExp(h,n)==0))
550    {
551      m=pHead(h);
552      pGetExpV(m,me);
553      exponent=me[n];
554      me[n]=0;
555      for(i=pVariables;i>0;i--)
556        me[i]+=exponent*ee[i];
557      pSetExpV(m,me);
558      if (e!=NULL)
559      {
560        nPower(pGetCoeff(e),exponent,&nu);
561        nu1=nMult(pGetCoeff(m),nu);
562        nDelete(&nu);
563        pSetCoeff(m,nu1);
564      }
565      res=pAdd(res,m);
566    }
567    pLmDelete(&h);
568  }
569  omFreeSize((ADDRESS)me,(pVariables+1)*sizeof(int));
570  omFreeSize((ADDRESS)ee,(pVariables+1)*sizeof(int));
571  return res;
572}
573
574/* Returns TRUE if
575     * LM(p) | LM(lcm)
576     * LC(p) | LC(lcm) only if ring
577     * Exists i, j:
578         * LE(p, i)  != LE(lcm, i)
579         * LE(p1, i) != LE(lcm, i)   ==> LCM(p1, p) != lcm
580         * LE(p, j)  != LE(lcm, j)
581         * LE(p2, j) != LE(lcm, j)   ==> LCM(p2, p) != lcm
582*/
583BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
584{
585  int k, j;
586
587  if (lcm==NULL) return FALSE;
588
589  for (j=pVariables; j; j--)
590    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
591  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
592  for (j=pVariables; j; j--)
593  {
594    if (pGetExp(p1,j)!=pGetExp(lcm,j))
595    {
596      if (pGetExp(p,j)!=pGetExp(lcm,j))
597      {
598        for (k=pVariables; k>j; k--)
599        {
600          if ((pGetExp(p,k)!=pGetExp(lcm,k))
601          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
602            return TRUE;
603        }
604        for (k=j-1; k; k--)
605        {
606          if ((pGetExp(p,k)!=pGetExp(lcm,k))
607          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
608            return TRUE;
609        }
610        return FALSE;
611      }
612    }
613    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
614    {
615      if (pGetExp(p,j)!=pGetExp(lcm,j))
616      {
617        for (k=pVariables; k>j; k--)
618        {
619          if ((pGetExp(p,k)!=pGetExp(lcm,k))
620          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
621            return TRUE;
622        }
623        for (k=j-1; k!=0 ; k--)
624        {
625          if ((pGetExp(p,k)!=pGetExp(lcm,k))
626          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
627            return TRUE;
628        }
629        return FALSE;
630      }
631    }
632  }
633  return FALSE;
634}
635#ifdef HAVE_RATGRING
636BOOLEAN pCompareChainPart (poly p,poly p1,poly p2,poly lcm)
637{
638  int k, j;
639
640  if (lcm==NULL) return FALSE;
641
642  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
643    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
644  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
645  for (j=currRing->real_var_end; j>=currRing->real_var_start; j--)
646  {
647    if (pGetExp(p1,j)!=pGetExp(lcm,j))
648    {
649      if (pGetExp(p,j)!=pGetExp(lcm,j))
650      {
651        for (k=pVariables; k>j; k--)
652        for (k=currRing->real_var_end; k>j; k--)
653        {
654          if ((pGetExp(p,k)!=pGetExp(lcm,k))
655          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
656            return TRUE;
657        }
658        for (k=j-1; k>=currRing->real_var_start; k--)
659        {
660          if ((pGetExp(p,k)!=pGetExp(lcm,k))
661          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
662            return TRUE;
663        }
664        return FALSE;
665      }
666    }
667    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
668    {
669      if (pGetExp(p,j)!=pGetExp(lcm,j))
670      {
671        for (k=currRing->real_var_end; k>j; k--)
672        {
673          if ((pGetExp(p,k)!=pGetExp(lcm,k))
674          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
675            return TRUE;
676        }
677        for (k=j-1; k>=currRing->real_var_start; k--)
678        {
679          if ((pGetExp(p,k)!=pGetExp(lcm,k))
680          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
681            return TRUE;
682        }
683        return FALSE;
684      }
685    }
686  }
687  return FALSE;
688}
689#endif
690
691/*2
692* returns the length of a (numbers of monomials)
693* respect syzComp
694*/
695poly pLast(poly a, int &l)
696{
697  if (a == NULL)
698  {
699    l = 0;
700    return NULL;
701  }
702  l = 1;
703  if (! rIsSyzIndexRing(currRing))
704  {
705    while (pNext(a)!=NULL)
706    {
707      pIter(a);
708      l++;
709    }
710  }
711  else
712  {
713    int curr_limit = rGetCurrSyzLimit(currRing);
714    poly pp = a;
715    while ((a=pNext(a))!=NULL)
716    {
717      if (pGetComp(a)<=curr_limit/*syzComp*/)
718        l++;
719      else break;
720      pp = a;
721    }
722    a=pp;
723  }
724  return a;
725}
726
Note: See TracBrowser for help on using the repository browser.