source: git/Singular/p_polys.cc @ 7d51c4

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