source: git/kernel/p_polys.cc @ bcd38b

spielwiese
Last change on this file since bcd38b was 8c5988, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: const ring git-svn-id: file:///usr/local/Singular/svn/trunk@7851 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 17.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    p_polys.cc
6 *  Purpose: implementation of currRing independent poly procedures
7 *  Author:  obachman (Olaf Bachmann)
8 *  Created: 8/00
9 *  Version: $Id: p_polys.cc,v 1.2 2005-04-20 17:25:51 Singular Exp $
10 *******************************************************************/
11
12#include "mod2.h"
13#include "structs.h"
14#include "p_polys.h"
15#include "ring.h"
16#include "febase.h"
17
18/***************************************************************
19 *
20 * Completing what needs to be set for the monomial
21 *
22 ***************************************************************/
23// this is special for the syz stuff
24static int* _Components = NULL;
25static long* _ShiftedComponents = NULL;
26static int _ExternalComponents = 0;
27
28void p_Setm_General(poly p, ring r)
29{
30  p_LmCheckPolyRing(p, r);
31  int pos=0;
32  if (r->typ!=NULL)
33  {
34    loop
35    {
36      long ord=0;
37      sro_ord* o=&(r->typ[pos]);
38      switch(o->ord_typ)
39      {
40        case ro_dp:
41        {
42          int a,e;
43          a=o->data.dp.start;
44          e=o->data.dp.end;
45          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r);
46          p->exp[o->data.dp.place]=ord;
47          break;
48        }
49        case ro_wp_neg:
50          ord=POLY_NEGWEIGHT_OFFSET;
51          // no break;
52        case ro_wp:
53        {
54          int a,e;
55          a=o->data.wp.start;
56          e=o->data.wp.end;
57          int *w=o->data.wp.weights;
58          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r)*w[i-a];
59          p->exp[o->data.wp.place]=ord;
60          break;
61        }
62        case ro_cp:
63        {
64          int a,e;
65          a=o->data.cp.start;
66          e=o->data.cp.end;
67          int pl=o->data.cp.place;
68          for(int i=a;i<=e;i++) { p->exp[pl]=p_GetExp(p,i,r); pl++; }
69          break;
70        }
71        case ro_syzcomp:
72        {
73          int c=p_GetComp(p,r);
74          long sc = c;
75          int* Components = (_ExternalComponents ? _Components :
76                             o->data.syzcomp.Components);
77          long* ShiftedComponents = (_ExternalComponents ? _ShiftedComponents:
78                                     o->data.syzcomp.ShiftedComponents);
79          if (ShiftedComponents != NULL)
80          {
81            assume(Components != NULL);
82            assume(c == 0 || Components[c] != 0);
83            sc = ShiftedComponents[Components[c]];
84            assume(c == 0 || sc != 0);
85          }
86          p->exp[o->data.syzcomp.place]=sc;
87          break;
88        }
89        case ro_syz:
90        {
91          int c=p_GetComp(p, r);
92          if (c > o->data.syz.limit)
93            p->exp[o->data.syz.place] = o->data.syz.curr_index;
94          else if (c > 0)
95            p->exp[o->data.syz.place]= o->data.syz.syz_index[c];
96          else
97          {
98            assume(c == 0);
99            p->exp[o->data.syz.place]= 0;
100          }
101          break;
102        }
103        default:
104          dReportError("wrong ord in rSetm:%d\n",o->ord_typ);
105          return;
106      }
107      pos++;
108      if (pos == r->OrdSize) return;
109    }
110  }
111}
112
113void p_Setm_Syz(poly p, ring r, int* Components, long* ShiftedComponents)
114{
115  _Components = Components;
116  _ShiftedComponents = ShiftedComponents;
117  _ExternalComponents = 1;
118  p_Setm_General(p, r);
119  _ExternalComponents = 0;
120}
121
122// dummy for lp, ls, etc
123void p_Setm_Dummy(poly p, ring r)
124{
125  p_LmCheckPolyRing(p, r);
126}
127
128// for dp, Dp, ds, etc
129void p_Setm_TotalDegree(poly p, ring r)
130{
131  p_LmCheckPolyRing(p, r);
132  p->exp[r->pOrdIndex] = p_ExpVectorQuerSum(p, r);
133}
134
135// for wp, Wp, ws, etc
136void p_Setm_WFirstTotalDegree(poly p, ring r)
137{
138  p_LmCheckPolyRing(p, r);
139  p->exp[r->pOrdIndex] = pWFirstTotalDegree(p, r);
140}
141
142p_SetmProc p_GetSetmProc(ring r)
143{
144  // covers lp, rp, ls,
145  if (r->typ == NULL) return p_Setm_Dummy;
146
147  if (r->OrdSize == 1)
148  {
149    if (r->typ[0].ord_typ == ro_dp &&
150        r->typ[0].data.dp.start == 1 &&
151        r->typ[0].data.dp.end == r->N &&
152        r->typ[0].data.dp.place == r->pOrdIndex)
153      return p_Setm_TotalDegree;
154    if (r->typ[0].ord_typ == ro_wp &&
155        r->typ[0].data.wp.start == 1 &&
156        r->typ[0].data.wp.end == r->N &&
157        r->typ[0].data.wp.place == r->pOrdIndex &&
158        r->typ[0].data.wp.weights == r->firstwv)
159      return p_Setm_WFirstTotalDegree;
160  }
161  return p_Setm_General;
162}
163
164
165/* -------------------------------------------------------------------*/
166/* several possibilities for pFDeg: the degree of the head term       */
167/*2
168* compute the degree of the leading monomial of p
169* the ordering is compatible with degree, use a->order
170*/
171inline long _pDeg(poly a, const ring r)
172{
173  p_LmCheckPolyRing(a, r);
174  assume(p_GetOrder(a, r) == pWTotaldegree(a, r));
175  return p_GetOrder(a, r);
176}
177
178long pDeg(poly a, const ring r)
179{
180  return _pDeg(a, r);
181}
182
183/*2
184* compute the degree of the leading monomial of p
185* with respect to weigths 1
186* (all are 1 so save multiplications or they are of different signs)
187* the ordering is not compatible with degree so do not use p->Order
188*/
189inline long _pTotaldegree(poly p, ring r)
190{
191  p_LmCheckPolyRing(p, r);
192  return (long) p_ExpVectorQuerSum(p, r);
193}
194long pTotaldegree(poly p, ring r)
195{
196  return (long) _pTotaldegree(p, r);
197}
198
199// pWTotalDegree for weighted orderings
200// whose first block covers all variables
201inline long _pWFirstTotalDegree(poly p, ring r)
202{
203  int i;
204  long sum = 0;
205
206  for (i=1; i<= r->firstBlockEnds; i++)
207  {
208    sum += p_GetExp(p, i, r)*r->firstwv[i-1];
209  }
210  return sum;
211}
212
213long pWFirstTotalDegree(poly p, ring r)
214{
215  return (long) _pWFirstTotalDegree(p, r);
216}
217
218/*2
219* compute the degree of the leading monomial of p
220* with respect to weigths from the ordering
221* the ordering is not compatible with degree so do not use p->Order
222*/
223long pWTotaldegree(poly p, ring r)
224{
225  p_LmCheckPolyRing(p, r);
226  int i, k;
227  long j =0;
228
229  // iterate through each block:
230  for (i=0;r->order[i]!=0;i++)
231  {
232    switch(r->order[i])
233    {
234      case ringorder_wp:
235      case ringorder_ws:
236      case ringorder_Wp:
237      case ringorder_Ws:
238        for (k=r->block0[i];k<=r->block1[i];k++)
239        { // in jedem block:
240          j+= p_GetExp(p,k,r)*r->wvhdl[i][k - r->block0[i]];
241        }
242        break;
243      case ringorder_M:
244      case ringorder_lp:
245      case ringorder_ls:
246      case ringorder_dp:
247      case ringorder_ds:
248      case ringorder_Dp:
249      case ringorder_Ds:
250      case ringorder_rp:
251        for (k=r->block0[i];k<=r->block1[i];k++)
252        {
253          j+= p_GetExp(p,k,r);
254        }
255        break;
256      case ringorder_c:
257      case ringorder_C:
258      case ringorder_S:
259      case ringorder_s:
260      case ringorder_aa:
261        break;
262      case ringorder_a:
263        for (k=r->block0[i];k<=r->block1[i];k++)
264        { // only one line
265          j+= p_GetExp(p,k, r)*r->wvhdl[i][ k- r->block0[i]];
266        }
267        return j;
268#ifndef NDEBUG
269      default:
270        Print("missing order %d in pWTotaldegree\n",r->order[i]);
271        break;
272#endif
273    }
274  }
275  return  j;
276}
277
278int pWeight(int i, const ring r)
279{
280  if ((r->firstwv==NULL) || (i>r->firstBlockEnds))
281  {
282    return 1;
283  }
284  return r->firstwv[i-1];
285}
286
287long pWDegree(poly p, ring r)
288{
289  if (r->firstwv==NULL) return pTotaldegree(p, r);
290  p_LmCheckPolyRing(p, r);
291  int i, k;
292  long j =0;
293
294  for(i=1;i<=r->firstBlockEnds;i++)
295    j+=p_GetExp(p, i, r)*r->firstwv[i-1];
296
297  for (;i<=r->N;i++)
298    j+=p_GetExp(p,i, r)*pWeight(i, r);
299
300  return j;
301}
302
303
304/* ---------------------------------------------------------------------*/
305/* several possibilities for pLDeg: the maximal degree of a monomial in p*/
306/*  compute in l also the pLength of p                                   */
307
308/*2
309* compute the length of a polynomial (in l)
310* and the degree of the monomial with maximal degree: the last one
311*/
312long pLDeg0(poly p,int *l, ring r)
313{
314  p_CheckPolyRing(p, r);
315  Exponent_t k= p_GetComp(p, r);
316  int ll=1;
317
318  if (k > 0)
319  {
320    while ((pNext(p)!=NULL) && (p_GetComp(pNext(p), r)==k))
321    {
322      pIter(p);
323      ll++;
324    }
325  }
326  else
327  {
328     while (pNext(p)!=NULL)
329     {
330       pIter(p);
331       ll++;
332     }
333  }
334  *l=ll;
335  return r->pFDeg(p, r);
336}
337
338/*2
339* compute the length of a polynomial (in l)
340* and the degree of the monomial with maximal degree: the last one
341* but search in all components before syzcomp
342*/
343long pLDeg0c(poly p,int *l, ring r)
344{
345  assume(p!=NULL);
346#ifdef PDEBUG
347  _p_Test(p,r,PDEBUG);
348#endif
349  p_CheckPolyRing(p, r);
350  long o;
351  int ll=1;
352
353  if (! rIsSyzIndexRing(r))
354  {
355    while (pNext(p) != NULL)
356    {
357      pIter(p);
358      ll++;
359    }
360    o = r->pFDeg(p, r);
361  }
362  else
363  {
364    int curr_limit = rGetCurrSyzLimit(r);
365    poly pp = p;
366    while ((p=pNext(p))!=NULL)
367    {
368      if (p_GetComp(p, r)<=curr_limit/*syzComp*/)
369        ll++;
370      else break;
371      pp = p;
372    }
373#ifdef PDEBUG
374    _p_Test(pp,r,PDEBUG);
375#endif
376    o = r->pFDeg(pp, r);
377  }
378  *l=ll;
379  return o;
380}
381
382/*2
383* compute the length of a polynomial (in l)
384* and the degree of the monomial with maximal degree: the first one
385* this works for the polynomial case with degree orderings
386* (both c,dp and dp,c)
387*/
388long pLDegb(poly p,int *l, ring r)
389{
390  p_CheckPolyRing(p, r);
391  Exponent_t k= p_GetComp(p, r);
392  long o = r->pFDeg(p, r);
393  int ll=1;
394
395  if (k != 0)
396  {
397    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
398    {
399      ll++;
400    }
401  }
402  else
403  {
404    while ((p=pNext(p)) !=NULL)
405    {
406      ll++;
407    }
408  }
409  *l=ll;
410  return o;
411}
412
413/*2
414* compute the length of a polynomial (in l)
415* and the degree of the monomial with maximal degree:
416* this is NOT the last one, we have to look for it
417*/
418long pLDeg1(poly p,int *l, ring r)
419{
420  p_CheckPolyRing(p, r);
421  Exponent_t k= p_GetComp(p, r);
422  int ll=1;
423  long  t,max;
424
425  max=r->pFDeg(p, r);
426  if (k > 0)
427  {
428    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
429    {
430      t=r->pFDeg(p, r);
431      if (t>max) max=t;
432      ll++;
433    }
434  }
435  else
436  {
437    while ((p=pNext(p))!=NULL)
438    {
439      t=r->pFDeg(p, r);
440      if (t>max) max=t;
441      ll++;
442    }
443  }
444  *l=ll;
445  return max;
446}
447
448/*2
449* compute the length of a polynomial (in l)
450* and the degree of the monomial with maximal degree:
451* this is NOT the last one, we have to look for it
452* in all components
453*/
454long pLDeg1c(poly p,int *l, ring r)
455{
456  p_CheckPolyRing(p, r);
457  int ll=1;
458  long  t,max;
459
460  max=r->pFDeg(p, r);
461  if (rIsSyzIndexRing(r))
462  {
463    long limit = rGetCurrSyzLimit(r);
464    while ((p=pNext(p))!=NULL)
465    {
466      if (p_GetComp(p, r)<=limit)
467      {
468        if ((t=r->pFDeg(p, r))>max) max=t;
469        ll++;
470      }
471      else break;
472    }
473  }
474  else
475  {
476    while ((p=pNext(p))!=NULL)
477    {
478      if ((t=r->pFDeg(p, r))>max) max=t;
479      ll++;
480    }
481  }
482  *l=ll;
483  return max;
484}
485
486// like pLDeg1, only pFDeg == pDeg
487long pLDeg1_Deg(poly p,int *l, ring r)
488{
489  assume(r->pFDeg == pDeg);
490  p_CheckPolyRing(p, r);
491  Exponent_t k= p_GetComp(p, r);
492  int ll=1;
493  long  t,max;
494
495  max=_pDeg(p, r);
496  if (k > 0)
497  {
498    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
499    {
500      t=_pDeg(p, r);
501      if (t>max) max=t;
502      ll++;
503    }
504  }
505  else
506  {
507    while ((p=pNext(p))!=NULL)
508    {
509      t=_pDeg(p, r);
510      if (t>max) max=t;
511      ll++;
512    }
513  }
514  *l=ll;
515  return max;
516}
517
518long pLDeg1c_Deg(poly p,int *l, ring r)
519{
520  assume(r->pFDeg == pDeg);
521  p_CheckPolyRing(p, r);
522  int ll=1;
523  long  t,max;
524
525  max=_pDeg(p, r);
526  if (rIsSyzIndexRing(r))
527  {
528    long limit = rGetCurrSyzLimit(r);
529    while ((p=pNext(p))!=NULL)
530    {
531      if (p_GetComp(p, r)<=limit)
532      {
533        if ((t=_pDeg(p, r))>max) max=t;
534        ll++;
535      }
536      else break;
537    }
538  }
539  else
540  {
541    while ((p=pNext(p))!=NULL)
542    {
543      if ((t=_pDeg(p, r))>max) max=t;
544      ll++;
545    }
546  }
547  *l=ll;
548  return max;
549}
550
551// like pLDeg1, only pFDeg == pTotoalDegree
552long pLDeg1_Totaldegree(poly p,int *l, ring r)
553{
554  p_CheckPolyRing(p, r);
555  Exponent_t k= p_GetComp(p, r);
556  int ll=1;
557  long  t,max;
558
559  max=_pTotaldegree(p, r);
560  if (k > 0)
561  {
562    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
563    {
564      t=_pTotaldegree(p, r);
565      if (t>max) max=t;
566      ll++;
567    }
568  }
569  else
570  {
571    while ((p=pNext(p))!=NULL)
572    {
573      t=_pTotaldegree(p, r);
574      if (t>max) max=t;
575      ll++;
576    }
577  }
578  *l=ll;
579  return max;
580}
581
582long pLDeg1c_Totaldegree(poly p,int *l, ring r)
583{
584  p_CheckPolyRing(p, r);
585  int ll=1;
586  long  t,max;
587
588  max=_pTotaldegree(p, r);
589  if (rIsSyzIndexRing(r))
590  {
591    long limit = rGetCurrSyzLimit(r);
592    while ((p=pNext(p))!=NULL)
593    {
594      if (p_GetComp(p, r)<=limit)
595      {
596        if ((t=_pTotaldegree(p, r))>max) max=t;
597        ll++;
598      }
599      else break;
600    }
601  }
602  else
603  {
604    while ((p=pNext(p))!=NULL)
605    {
606      if ((t=_pTotaldegree(p, r))>max) max=t;
607      ll++;
608    }
609  }
610  *l=ll;
611  return max;
612}
613
614// like pLDeg1, only pFDeg == pWFirstTotalDegree
615long pLDeg1_WFirstTotalDegree(poly p,int *l, ring r)
616{
617  p_CheckPolyRing(p, r);
618  Exponent_t k= p_GetComp(p, r);
619  int ll=1;
620  long  t,max;
621
622  max=_pWFirstTotalDegree(p, r);
623  if (k > 0)
624  {
625    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
626    {
627      t=_pWFirstTotalDegree(p, r);
628      if (t>max) max=t;
629      ll++;
630    }
631  }
632  else
633  {
634    while ((p=pNext(p))!=NULL)
635    {
636      t=_pWFirstTotalDegree(p, r);
637      if (t>max) max=t;
638      ll++;
639    }
640  }
641  *l=ll;
642  return max;
643}
644
645long pLDeg1c_WFirstTotalDegree(poly p,int *l, ring r)
646{
647  p_CheckPolyRing(p, r);
648  int ll=1;
649  long  t,max;
650
651  max=_pWFirstTotalDegree(p, r);
652  if (rIsSyzIndexRing(r))
653  {
654    long limit = rGetCurrSyzLimit(r);
655    while ((p=pNext(p))!=NULL)
656    {
657      if (p_GetComp(p, r)<=limit)
658      {
659        if ((t=_pTotaldegree(p, r))>max) max=t;
660        ll++;
661      }
662      else break;
663    }
664  }
665  else
666  {
667    while ((p=pNext(p))!=NULL)
668    {
669      if ((t=_pTotaldegree(p, r))>max) max=t;
670      ll++;
671    }
672  }
673  *l=ll;
674  return max;
675}
676
677/***************************************************************
678 *
679 * Maximal Exponent business
680 *
681 ***************************************************************/
682
683static inline unsigned long
684p_GetMaxExpL2(unsigned long l1, unsigned long l2, ring r,
685              unsigned long number_of_exp)
686{
687  const unsigned long bitmask = r->bitmask;
688  unsigned long ml1 = l1 & bitmask;
689  unsigned long ml2 = l2 & bitmask;
690  unsigned long max = (ml1 > ml2 ? ml1 : ml2);
691  unsigned long j = number_of_exp - 1;
692
693  if (j > 0)
694  {
695    unsigned long mask = bitmask << r->BitsPerExp;
696    while (1)
697    {
698      ml1 = l1 & mask;
699      ml2 = l2 & mask;
700      max |= ((ml1 > ml2 ? ml1 : ml2) & mask);
701      j--;
702      if (j == 0) break;
703      mask = mask << r->BitsPerExp;
704    }
705  }
706  return max;
707}
708
709static inline unsigned long
710p_GetMaxExpL2(unsigned long l1, unsigned long l2, ring r)
711{
712  return p_GetMaxExpL2(l1, l2, r, r->ExpPerLong);
713}
714
715poly p_GetMaxExpP(poly p, ring r)
716{
717  p_CheckPolyRing(p, r);
718  if (p == NULL) return p_Init(r);
719  poly max = p_LmInit(p, r);
720  pIter(p);
721  if (p == NULL) return max;
722  int i, offset;
723  unsigned long l_p, l_max;
724  unsigned long divmask = r->divmask;
725
726  do
727  {
728    offset = r->VarL_Offset[0];
729    l_p = p->exp[offset];
730    l_max = max->exp[offset];
731    // do the divisibility trick to find out whether l has an exponent
732    if (l_p > l_max ||
733        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
734      max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
735
736    for (i=1; i<r->VarL_Size; i++)
737    {
738      offset = r->VarL_Offset[i];
739      l_p = p->exp[offset];
740      l_max = max->exp[offset];
741      // do the divisibility trick to find out whether l has an exponent
742      if (l_p > l_max ||
743          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
744        max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
745    }
746    pIter(p);
747  }
748  while (p != NULL);
749  return max;
750}
751
752unsigned long p_GetMaxExpL(poly p, ring r, unsigned long l_max)
753{
754  unsigned long l_p, divmask = r->divmask;
755  int i;
756
757  while (p != NULL)
758  {
759    l_p = p->exp[r->VarL_Offset[0]];
760    if (l_p > l_max ||
761        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
762      l_max = p_GetMaxExpL2(l_max, l_p, r);
763    for (i=1; i<r->VarL_Size; i++)
764    {
765      l_p = p->exp[r->VarL_Offset[i]];
766      // do the divisibility trick to find out whether l has an exponent
767      if (l_p > l_max ||
768          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
769        l_max = p_GetMaxExpL2(l_max, l_p, r);
770    }
771    pIter(p);
772  }
773  return l_max;
774}
775
776/***************************************************************
777 *
778 * Misc things
779 *
780 ***************************************************************/
781// returns TRUE, if all monoms have the same component
782BOOLEAN p_OneComp(poly p, ring r)
783{
784  if(p!=NULL)
785  {
786    long i = p_GetComp(p, r);
787    while (pNext(p)!=NULL)
788    {
789      pIter(p);
790      if(i != p_GetComp(p, r)) return FALSE;
791    }
792  }
793  return TRUE;
794}
795
796/*2
797*test if a monomial /head term is a pure power
798*/
799int p_IsPurePower(const poly p, const ring r)
800{
801  int i,k=0;
802
803  for (i=r->N;i;i--)
804  {
805    if (p_GetExp(p,i, r)!=0)
806    {
807      if(k!=0) return 0;
808      k=i;
809    }
810  }
811  return k;
812}
813
814/*2
815* returns a polynomial representing the integer i
816*/
817poly p_ISet(int i, ring r)
818{
819  poly rc = NULL;
820  if (i!=0)
821  {
822    rc = p_Init(r);
823    pSetCoeff0(rc,r->cf->nInit(i));
824    if (r->cf->nIsZero(p_GetCoeff(rc,r)))
825      p_DeleteLm(&rc,r);
826  }
827  return rc;
828}
829
830/*2
831* returns a polynomial representing the number n
832* destroys n
833*/
834poly p_NSet(number n, ring r)
835{
836  if (r->cf->nIsZero(n))
837  {
838    r->cf->cfDelete(&n, r);
839    return NULL;
840  }
841  else
842  {
843    poly rc = p_Init(r);
844    pSetCoeff0(rc,n);
845    return rc;
846  }
847}
848
849/***************************************************************
850 *
851 * p_ShallowDelete
852 *
853 ***************************************************************/
854#undef LINKAGE
855#define LINKAGE
856#undef p_Delete
857#define p_Delete p_ShallowDelete
858#undef n_Delete
859#define n_Delete(n, r) ((void)0)
860
861#include "p_Delete__T.cc"
862
Note: See TracBrowser for help on using the repository browser.