source: git/kernel/p_polys.cc @ 99bdcf

jengelh-datetimespielwiese
Last change on this file since 99bdcf was 99bdcf, checked in by Hans Schoenemann <hannes@…>, 13 years ago
pTotaldegree/p_Totaldegree git-svn-id: file:///usr/local/Singular/svn/trunk@13084 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 24.8 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$
10 *******************************************************************/
11
12
13#include <kernel/mod2.h>
14
15#ifndef NDEBUG
16# define MYTEST 0
17#else /* ifndef NDEBUG */
18# define MYTEST 0
19#endif /* ifndef NDEBUG */
20
21#include <kernel/structs.h>
22#include <kernel/structs.h>
23#include <kernel/p_polys.h>
24#include <kernel/ring.h>
25#include <kernel/ideals.h>
26#include <kernel/int64vec.h>
27#ifndef NDEBUG
28#include <kernel/febase.h>
29#endif
30
31/***************************************************************
32 *
33 * Completing what needs to be set for the monomial
34 *
35 ***************************************************************/
36// this is special for the syz stuff
37static int* _Components = NULL;
38static long* _ShiftedComponents = NULL;
39static int _ExternalComponents = 0;
40
41BOOLEAN pSetm_error=0;
42
43void p_Setm_General(poly p, const ring r)
44{
45  p_LmCheckPolyRing(p, r);
46  int pos=0;
47  if (r->typ!=NULL)
48  {
49    loop
50    {
51      long ord=0;
52      sro_ord* o=&(r->typ[pos]);
53      switch(o->ord_typ)
54      {
55        case ro_dp:
56        {
57          int a,e;
58          a=o->data.dp.start;
59          e=o->data.dp.end;
60          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r);
61          p->exp[o->data.dp.place]=ord;
62          break;
63        }
64        case ro_wp_neg:
65          ord=POLY_NEGWEIGHT_OFFSET;
66          // no break;
67        case ro_wp:
68        {
69          int a,e;
70          a=o->data.wp.start;
71          e=o->data.wp.end;
72          int *w=o->data.wp.weights;
73#if 1
74          for(int i=a;i<=e;i++) ord+=p_GetExp(p,i,r)*w[i-a];
75#else
76          long ai;
77          int ei,wi;
78          for(int i=a;i<=e;i++)
79          {
80             ei=p_GetExp(p,i,r);
81             wi=w[i-a];
82             ai=ei*wi;
83             if (ai/ei!=wi) pSetm_error=TRUE;
84             ord+=ai;
85             if (ord<ai) pSetm_error=TRUE;
86          }
87#endif
88          p->exp[o->data.wp.place]=ord;
89          break;
90        }
91      case ro_wp64:
92        {
93          int64 ord=0;
94          int a,e;
95          a=o->data.wp64.start;
96          e=o->data.wp64.end;
97          int64 *w=o->data.wp64.weights64;
98          int64 ei,wi,ai;
99          for(int i=a;i<=e;i++)
100          {
101            //Print("exp %d w %d \n",p_GetExp(p,i,r),(int)w[i-a]);
102            //ord+=((int64)p_GetExp(p,i,r))*w[i-a];
103            ei=(int64)p_GetExp(p,i,r);
104            wi=w[i-a];
105            ai=ei*wi;
106            if(ei!=0 && ai/ei!=wi)
107            {
108              pSetm_error=TRUE;
109              #if SIZEOF_LONG == 4
110              Print("ai %lld, wi %lld\n",ai,wi);
111              #else
112              Print("ai %ld, wi %ld\n",ai,wi);
113              #endif
114            }
115            ord+=ai;
116            if (ord<ai)
117            {
118              pSetm_error=TRUE;
119              #if SIZEOF_LONG == 4
120              Print("ai %lld, ord %lld\n",ai,ord);
121              #else
122              Print("ai %ld, ord %ld\n",ai,ord);
123              #endif
124            }
125          }
126          int64 mask=(int64)0x7fffffff;
127          long a_0=(long)(ord&mask); //2^31
128          long a_1=(long)(ord >>31 ); /*(ord/(mask+1));*/
129
130          //Print("mask: %x,  ord: %d,  a_0: %d,  a_1: %d\n"
131          //,(int)mask,(int)ord,(int)a_0,(int)a_1);
132                    //Print("mask: %d",mask);
133
134          p->exp[o->data.wp64.place]=a_1;
135          p->exp[o->data.wp64.place+1]=a_0;
136//            if(p_Setm_error) Print("***************************\n
137//                                    ***************************\n
138//                                    **WARNING: overflow error**\n
139//                                    ***************************\n
140//                                    ***************************\n");
141          break;
142        }
143        case ro_cp:
144        {
145          int a,e;
146          a=o->data.cp.start;
147          e=o->data.cp.end;
148          int pl=o->data.cp.place;
149          for(int i=a;i<=e;i++) { p->exp[pl]=p_GetExp(p,i,r); pl++; }
150          break;
151        }
152        case ro_syzcomp:
153        {
154          int c=p_GetComp(p,r);
155          long sc = c;
156          int* Components = (_ExternalComponents ? _Components :
157                             o->data.syzcomp.Components);
158          long* ShiftedComponents = (_ExternalComponents ? _ShiftedComponents:
159                                     o->data.syzcomp.ShiftedComponents);
160          if (ShiftedComponents != NULL)
161          {
162            assume(Components != NULL);
163            assume(c == 0 || Components[c] != 0);
164            sc = ShiftedComponents[Components[c]];
165            assume(c == 0 || sc != 0);
166          }
167          p->exp[o->data.syzcomp.place]=sc;
168          break;
169        }
170        case ro_syz:
171        {
172          int c=p_GetComp(p, r);
173          if (c > o->data.syz.limit)
174            p->exp[o->data.syz.place] = o->data.syz.curr_index;
175          else if (c > 0)
176            p->exp[o->data.syz.place]= o->data.syz.syz_index[c];
177          else
178          {
179            assume(c == 0);
180            p->exp[o->data.syz.place]= 0;
181          }
182          break;
183        }
184        // Prefix for Induced Schreyer ordering
185        case ro_isTemp: // Do nothing?? (to be removed into suffix later on...?)
186        {
187          assume(p != NULL);
188
189#ifndef NDEBUG
190#if MYTEST
191          Print("isTemp ord in rSetm: pos: %d, p: ", pos);  p_DebugPrint(p, r, r, 0);
192#endif
193#endif
194          int c = p_GetComp(p, r);
195
196          assume( c >= 0 );
197          const int limit = o->data.is.limit;
198
199          assume( limit >= 0 );
200
201          // Let's simulate case ro_syz above....
202          // Should accumulate (by Suffix) and be a level indicator
203          const int* const pVarOffset = o->data.isTemp.pVarOffset;
204
205          assume( pVarOffset != NULL );
206
207          if( c > limit )
208            p->exp[o->data.isTemp.start] = 1;
209          else
210          {
211            p->exp[o->data.isTemp.start] = 0;
212          }
213
214          // TODO: Can this be done in the suffix???
215          for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
216          {
217            const int vo = pVarOffset[i];
218            if( vo != -1) // TODO: optimize: can be done once!
219            {
220              p_SetExp(p, p_GetExp(p, i, r), vo, r); // copy put them verbatim
221              assume( p_GetExp(p, vo, r) == p_GetExp(p, i, r) ); // copy put them verbatim
222            }
223          }
224
225       
226
227#ifndef NDEBUG
228          for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
229          {
230            const int vo = pVarOffset[i];
231            if( vo != -1) // TODO: optimize: can be done once!
232            {
233              assume( p_GetExp(p, vo, r) == p_GetExp(p, i, r) ); // copy put them verbatim
234            }
235          }
236
237#if MYTEST
238//          if( p->exp[o->data.isTemp.start] > 0 )
239//          {
240//            PrintS("Initial Value: "); p_DebugPrint(p, r, r, 1);
241//          }
242#endif
243#endif
244
245          break;
246        }
247
248        // Suffix for Induced Schreyer ordering
249        case ro_is:
250        {
251          assume(p != NULL);
252
253          int c = p_GetComp(p, r);
254
255          assume( c >= 0 );
256          const ideal F = o->data.is.F;
257          const int limit = o->data.is.limit;
258
259          if( F != NULL && c > limit )
260          {
261#ifndef NDEBUG
262#if MYTEST
263            Print("is ord in rSetm: pos: %d, c: %d, limit: %d\n", c, pos, limit); // p_DebugPrint(p, r, r, 1);
264#endif
265#endif
266
267            c -= limit;
268            assume( c > 0 );
269            c--;
270
271            assume( c < IDELEMS(F) ); // What about others???
272
273            const poly pp = F->m[c]; // get reference monomial!!!
274
275
276#ifndef NDEBUG
277#if MYTEST
278            Print("Respective F[c - %d: %d] pp: ", limit, c); 
279            p_DebugPrint(pp, r, r, 1);
280#endif
281#endif
282
283
284            if(pp == NULL) break;
285
286            const int start = o->data.is.start;
287            const int end = o->data.is.end;
288
289            assume(start <= end);
290            assume(pp != NULL);
291
292            for( int i = start; i <= end; i++) // v[0] may be here...
293              p->exp[i] += pp->exp[i]; // !!!!!!!! ADD corresponding LT(F)
294
295#ifndef NDEBUG
296            const int* const pVarOffset = o->data.is.pVarOffset;
297
298            assume( pVarOffset != NULL );
299
300            for( int i = 1; i <= r->N; i++ ) // No v[0] here!!!
301            {
302              const int vo = pVarOffset[i];
303              if( vo != -1) // TODO: optimize: can be done once!
304                assume( p_GetExp(p, vo, r) == (p_GetExp(p, i, r) + p_GetExp(pp, vo, r)) );
305            }
306            // TODO: how to check this for computed values???
307#endif
308#ifndef NDEBUG
309#if MYTEST
310            PrintS("IS::Suffix::Result: "); // p_Write(p, r, r);
311            p_DebugPrint(p, r, r, 1);
312#endif
313#endif
314
315          } else
316          {
317            const int* const pVarOffset = o->data.is.pVarOffset;
318
319            // What about v[0] - component: it will be added later by
320            // suffix!!!
321            // TODO: Test it!
322            const int vo = pVarOffset[0];
323            if( vo != -1 )
324              p->exp[vo] = c; // initial component v[0]!
325          }
326
327          break;
328        }
329        default:
330          dReportError("wrong ord in rSetm:%d\n",o->ord_typ);
331          return;
332      }
333      pos++;
334      if (pos == r->OrdSize) return;
335    }
336  }
337}
338
339void p_Setm_Syz(poly p, ring r, int* Components, long* ShiftedComponents)
340{
341  _Components = Components;
342  _ShiftedComponents = ShiftedComponents;
343  _ExternalComponents = 1;
344  p_Setm_General(p, r);
345  _ExternalComponents = 0;
346}
347
348// dummy for lp, ls, etc
349void p_Setm_Dummy(poly p, const ring r)
350{
351  p_LmCheckPolyRing(p, r);
352}
353
354// for dp, Dp, ds, etc
355void p_Setm_TotalDegree(poly p, const ring r)
356{
357  p_LmCheckPolyRing(p, r);
358  p->exp[r->pOrdIndex] = p_Totaldegree(p, r);
359}
360
361// for wp, Wp, ws, etc
362void p_Setm_WFirstTotalDegree(poly p, const ring r)
363{
364  p_LmCheckPolyRing(p, r);
365  p->exp[r->pOrdIndex] = pWFirstTotalDegree(p, r);
366}
367
368p_SetmProc p_GetSetmProc(ring r)
369{
370  // covers lp, rp, ls,
371  if (r->typ == NULL) return p_Setm_Dummy;
372
373  if (r->OrdSize == 1)
374  {
375    if (r->typ[0].ord_typ == ro_dp &&
376        r->typ[0].data.dp.start == 1 &&
377        r->typ[0].data.dp.end == r->N &&
378        r->typ[0].data.dp.place == r->pOrdIndex)
379      return p_Setm_TotalDegree;
380    if (r->typ[0].ord_typ == ro_wp &&
381        r->typ[0].data.wp.start == 1 &&
382        r->typ[0].data.wp.end == r->N &&
383        r->typ[0].data.wp.place == r->pOrdIndex &&
384        r->typ[0].data.wp.weights == r->firstwv)
385      return p_Setm_WFirstTotalDegree;
386  }
387  return p_Setm_General;
388}
389
390
391/* -------------------------------------------------------------------*/
392/* several possibilities for pFDeg: the degree of the head term       */
393
394/* comptible with ordering */
395long pDeg(poly a, const ring r)
396{
397  p_LmCheckPolyRing(a, r);
398  assume(p_GetOrder(a, r) == pWTotaldegree(a, r));
399  return p_GetOrder(a, r);
400}
401
402// pWTotalDegree for weighted orderings
403// whose first block covers all variables
404static inline long _pWFirstTotalDegree(poly p, const ring r)
405{
406  int i;
407  long sum = 0;
408
409  for (i=1; i<= r->firstBlockEnds; i++)
410  {
411    sum += p_GetExp(p, i, r)*r->firstwv[i-1];
412  }
413  return sum;
414}
415
416long pWFirstTotalDegree(poly p, const ring r)
417{
418  return (long) _pWFirstTotalDegree(p, r);
419}
420
421/*2
422* compute the degree of the leading monomial of p
423* with respect to weigths from the ordering
424* the ordering is not compatible with degree so do not use p->Order
425*/
426long pWTotaldegree(poly p, const ring r)
427{
428  p_LmCheckPolyRing(p, r);
429  int i, k;
430  long j =0;
431
432  // iterate through each block:
433  for (i=0;r->order[i]!=0;i++)
434  {
435    int b0=r->block0[i];
436    int b1=r->block1[i];
437    switch(r->order[i])
438    {
439      case ringorder_M:
440        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
441        { // in jedem block:
442          j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/]*r->OrdSgn;
443        }
444        break;
445      case ringorder_wp:
446      case ringorder_ws:
447      case ringorder_Wp:
448      case ringorder_Ws:
449        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
450        { // in jedem block:
451          j+= p_GetExp(p,k,r)*r->wvhdl[i][k - b0 /*r->block0[i]*/];
452        }
453        break;
454      case ringorder_lp:
455      case ringorder_ls:
456      case ringorder_rs:
457      case ringorder_dp:
458      case ringorder_ds:
459      case ringorder_Dp:
460      case ringorder_Ds:
461      case ringorder_rp:
462        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
463        {
464          j+= p_GetExp(p,k,r);
465        }
466        break;
467      case ringorder_a64:
468        {
469          int64* w=(int64*)r->wvhdl[i];
470          for (k=0;k<=(b1 /*r->block1[i]*/ - b0 /*r->block0[i]*/);k++)
471          {
472            //there should be added a line which checks if w[k]>2^31
473            j+= p_GetExp(p,k+1, r)*(long)w[k];
474          }
475          //break;
476          return j;
477        }
478      case ringorder_c:
479      case ringorder_C:
480      case ringorder_S:
481      case ringorder_s:
482      case ringorder_IS:
483      case ringorder_aa:
484        break;
485      case ringorder_a:
486        for (k=b0 /*r->block0[i]*/;k<=b1 /*r->block1[i]*/;k++)
487        { // only one line
488          j+= p_GetExp(p,k, r)*r->wvhdl[i][ k- b0 /*r->block0[i]*/];
489        }
490        //break;
491        return j;
492
493#ifndef NDEBUG
494      default:
495        Print("missing order %d in pWTotaldegree\n",r->order[i]);
496        break;
497#endif
498    }
499  }
500  return  j;
501}
502
503int pWeight(int i, const ring r)
504{
505  if ((r->firstwv==NULL) || (i>r->firstBlockEnds))
506  {
507    return 1;
508  }
509  return r->firstwv[i-1];
510}
511
512long pWDegree(poly p, const ring r)
513{
514  if (r->firstwv==NULL) return p_Totaldegree(p, r);
515  p_LmCheckPolyRing(p, r);
516  int i, k;
517  long j =0;
518
519  for(i=1;i<=r->firstBlockEnds;i++)
520    j+=p_GetExp(p, i, r)*r->firstwv[i-1];
521
522  for (;i<=r->N;i++)
523    j+=p_GetExp(p,i, r)*pWeight(i, r);
524
525  return j;
526}
527
528
529/* ---------------------------------------------------------------------*/
530/* several possibilities for pLDeg: the maximal degree of a monomial in p*/
531/*  compute in l also the pLength of p                                   */
532
533/*2
534* compute the length of a polynomial (in l)
535* and the degree of the monomial with maximal degree: the last one
536*/
537long pLDeg0(poly p,int *l, const ring r)
538{
539  p_CheckPolyRing(p, r);
540  long k= p_GetComp(p, r);
541  int ll=1;
542
543  if (k > 0)
544  {
545    while ((pNext(p)!=NULL) && (p_GetComp(pNext(p), r)==k))
546    {
547      pIter(p);
548      ll++;
549    }
550  }
551  else
552  {
553     while (pNext(p)!=NULL)
554     {
555       pIter(p);
556       ll++;
557     }
558  }
559  *l=ll;
560  return r->pFDeg(p, r);
561}
562
563/*2
564* compute the length of a polynomial (in l)
565* and the degree of the monomial with maximal degree: the last one
566* but search in all components before syzcomp
567*/
568long pLDeg0c(poly p,int *l, const ring r)
569{
570  assume(p!=NULL);
571#ifdef PDEBUG
572  _p_Test(p,r,PDEBUG);
573#endif
574  p_CheckPolyRing(p, r);
575  long o;
576  int ll=1;
577
578  if (! rIsSyzIndexRing(r))
579  {
580    while (pNext(p) != NULL)
581    {
582      pIter(p);
583      ll++;
584    }
585    o = r->pFDeg(p, r);
586  }
587  else
588  {
589    int curr_limit = rGetCurrSyzLimit(r);
590    poly pp = p;
591    while ((p=pNext(p))!=NULL)
592    {
593      if (p_GetComp(p, r)<=curr_limit/*syzComp*/)
594        ll++;
595      else break;
596      pp = p;
597    }
598#ifdef PDEBUG
599    _p_Test(pp,r,PDEBUG);
600#endif
601    o = r->pFDeg(pp, r);
602  }
603  *l=ll;
604  return o;
605}
606
607/*2
608* compute the length of a polynomial (in l)
609* and the degree of the monomial with maximal degree: the first one
610* this works for the polynomial case with degree orderings
611* (both c,dp and dp,c)
612*/
613long pLDegb(poly p,int *l, const ring r)
614{
615  p_CheckPolyRing(p, r);
616  long k= p_GetComp(p, r);
617  long o = r->pFDeg(p, r);
618  int ll=1;
619
620  if (k != 0)
621  {
622    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
623    {
624      ll++;
625    }
626  }
627  else
628  {
629    while ((p=pNext(p)) !=NULL)
630    {
631      ll++;
632    }
633  }
634  *l=ll;
635  return o;
636}
637
638/*2
639* compute the length of a polynomial (in l)
640* and the degree of the monomial with maximal degree:
641* this is NOT the last one, we have to look for it
642*/
643long pLDeg1(poly p,int *l, const ring r)
644{
645  p_CheckPolyRing(p, r);
646  long k= p_GetComp(p, r);
647  int ll=1;
648  long  t,max;
649
650  max=r->pFDeg(p, r);
651  if (k > 0)
652  {
653    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
654    {
655      t=r->pFDeg(p, r);
656      if (t>max) max=t;
657      ll++;
658    }
659  }
660  else
661  {
662    while ((p=pNext(p))!=NULL)
663    {
664      t=r->pFDeg(p, r);
665      if (t>max) max=t;
666      ll++;
667    }
668  }
669  *l=ll;
670  return max;
671}
672
673/*2
674* compute the length of a polynomial (in l)
675* and the degree of the monomial with maximal degree:
676* this is NOT the last one, we have to look for it
677* in all components
678*/
679long pLDeg1c(poly p,int *l, const ring r)
680{
681  p_CheckPolyRing(p, r);
682  int ll=1;
683  long  t,max;
684
685  max=r->pFDeg(p, r);
686  if (rIsSyzIndexRing(r))
687  {
688    long limit = rGetCurrSyzLimit(r);
689    while ((p=pNext(p))!=NULL)
690    {
691      if (p_GetComp(p, r)<=limit)
692      {
693        if ((t=r->pFDeg(p, r))>max) max=t;
694        ll++;
695      }
696      else break;
697    }
698  }
699  else
700  {
701    while ((p=pNext(p))!=NULL)
702    {
703      if ((t=r->pFDeg(p, r))>max) max=t;
704      ll++;
705    }
706  }
707  *l=ll;
708  return max;
709}
710
711// like pLDeg1, only pFDeg == pDeg
712long pLDeg1_Deg(poly p,int *l, const ring r)
713{
714  assume(r->pFDeg == pDeg);
715  p_CheckPolyRing(p, r);
716  long k= p_GetComp(p, r);
717  int ll=1;
718  long  t,max;
719
720  max=p_GetOrder(p, r);
721  if (k > 0)
722  {
723    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
724    {
725      t=p_GetOrder(p, r);
726      if (t>max) max=t;
727      ll++;
728    }
729  }
730  else
731  {
732    while ((p=pNext(p))!=NULL)
733    {
734      t=p_GetOrder(p, r);
735      if (t>max) max=t;
736      ll++;
737    }
738  }
739  *l=ll;
740  return max;
741}
742
743long pLDeg1c_Deg(poly p,int *l, const ring r)
744{
745  assume(r->pFDeg == pDeg);
746  p_CheckPolyRing(p, r);
747  int ll=1;
748  long  t,max;
749
750  max=p_GetOrder(p, r);
751  if (rIsSyzIndexRing(r))
752  {
753    long limit = rGetCurrSyzLimit(r);
754    while ((p=pNext(p))!=NULL)
755    {
756      if (p_GetComp(p, r)<=limit)
757      {
758        if ((t=p_GetOrder(p, r))>max) max=t;
759        ll++;
760      }
761      else break;
762    }
763  }
764  else
765  {
766    while ((p=pNext(p))!=NULL)
767    {
768      if ((t=p_GetOrder(p, r))>max) max=t;
769      ll++;
770    }
771  }
772  *l=ll;
773  return max;
774}
775
776// like pLDeg1, only pFDeg == pTotoalDegree
777long pLDeg1_Totaldegree(poly p,int *l, const ring r)
778{
779  p_CheckPolyRing(p, r);
780  long k= p_GetComp(p, r);
781  int ll=1;
782  long  t,max;
783
784  max=p_Totaldegree(p, r);
785  if (k > 0)
786  {
787    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
788    {
789      t=p_Totaldegree(p, r);
790      if (t>max) max=t;
791      ll++;
792    }
793  }
794  else
795  {
796    while ((p=pNext(p))!=NULL)
797    {
798      t=p_Totaldegree(p, r);
799      if (t>max) max=t;
800      ll++;
801    }
802  }
803  *l=ll;
804  return max;
805}
806
807long pLDeg1c_Totaldegree(poly p,int *l, const ring r)
808{
809  p_CheckPolyRing(p, r);
810  int ll=1;
811  long  t,max;
812
813  max=p_Totaldegree(p, r);
814  if (rIsSyzIndexRing(r))
815  {
816    long limit = rGetCurrSyzLimit(r);
817    while ((p=pNext(p))!=NULL)
818    {
819      if (p_GetComp(p, r)<=limit)
820      {
821        if ((t=p_Totaldegree(p, r))>max) max=t;
822        ll++;
823      }
824      else break;
825    }
826  }
827  else
828  {
829    while ((p=pNext(p))!=NULL)
830    {
831      if ((t=p_Totaldegree(p, r))>max) max=t;
832      ll++;
833    }
834  }
835  *l=ll;
836  return max;
837}
838
839// like pLDeg1, only pFDeg == pWFirstTotalDegree
840long pLDeg1_WFirstTotalDegree(poly p,int *l, const ring r)
841{
842  p_CheckPolyRing(p, r);
843  long k= p_GetComp(p, r);
844  int ll=1;
845  long  t,max;
846
847  max=_pWFirstTotalDegree(p, r);
848  if (k > 0)
849  {
850    while (((p=pNext(p))!=NULL) && (p_GetComp(p, r)==k))
851    {
852      t=_pWFirstTotalDegree(p, r);
853      if (t>max) max=t;
854      ll++;
855    }
856  }
857  else
858  {
859    while ((p=pNext(p))!=NULL)
860    {
861      t=_pWFirstTotalDegree(p, r);
862      if (t>max) max=t;
863      ll++;
864    }
865  }
866  *l=ll;
867  return max;
868}
869
870long pLDeg1c_WFirstTotalDegree(poly p,int *l, const ring r)
871{
872  p_CheckPolyRing(p, r);
873  int ll=1;
874  long  t,max;
875
876  max=_pWFirstTotalDegree(p, r);
877  if (rIsSyzIndexRing(r))
878  {
879    long limit = rGetCurrSyzLimit(r);
880    while ((p=pNext(p))!=NULL)
881    {
882      if (p_GetComp(p, r)<=limit)
883      {
884        if ((t=p_Totaldegree(p, r))>max) max=t;
885        ll++;
886      }
887      else break;
888    }
889  }
890  else
891  {
892    while ((p=pNext(p))!=NULL)
893    {
894      if ((t=p_Totaldegree(p, r))>max) max=t;
895      ll++;
896    }
897  }
898  *l=ll;
899  return max;
900}
901
902/***************************************************************
903 *
904 * Maximal Exponent business
905 *
906 ***************************************************************/
907
908static inline unsigned long
909p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r,
910              unsigned long number_of_exp)
911{
912  const unsigned long bitmask = r->bitmask;
913  unsigned long ml1 = l1 & bitmask;
914  unsigned long ml2 = l2 & bitmask;
915  unsigned long max = (ml1 > ml2 ? ml1 : ml2);
916  unsigned long j = number_of_exp - 1;
917
918  if (j > 0)
919  {
920    unsigned long mask = bitmask << r->BitsPerExp;
921    while (1)
922    {
923      ml1 = l1 & mask;
924      ml2 = l2 & mask;
925      max |= ((ml1 > ml2 ? ml1 : ml2) & mask);
926      j--;
927      if (j == 0) break;
928      mask = mask << r->BitsPerExp;
929    }
930  }
931  return max;
932}
933
934static inline unsigned long
935p_GetMaxExpL2(unsigned long l1, unsigned long l2, const ring r)
936{
937  return p_GetMaxExpL2(l1, l2, r, r->ExpPerLong);
938}
939
940poly p_GetMaxExpP(poly p, const ring r)
941{
942  p_CheckPolyRing(p, r);
943  if (p == NULL) return p_Init(r);
944  poly max = p_LmInit(p, r);
945  pIter(p);
946  if (p == NULL) return max;
947  int i, offset;
948  unsigned long l_p, l_max;
949  unsigned long divmask = r->divmask;
950
951  do
952  {
953    offset = r->VarL_Offset[0];
954    l_p = p->exp[offset];
955    l_max = max->exp[offset];
956    // do the divisibility trick to find out whether l has an exponent
957    if (l_p > l_max ||
958        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
959      max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
960
961    for (i=1; i<r->VarL_Size; i++)
962    {
963      offset = r->VarL_Offset[i];
964      l_p = p->exp[offset];
965      l_max = max->exp[offset];
966      // do the divisibility trick to find out whether l has an exponent
967      if (l_p > l_max ||
968          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
969        max->exp[offset] = p_GetMaxExpL2(l_max, l_p, r);
970    }
971    pIter(p);
972  }
973  while (p != NULL);
974  return max;
975}
976
977unsigned long p_GetMaxExpL(poly p, const ring r, unsigned long l_max)
978{
979  unsigned long l_p, divmask = r->divmask;
980  int i;
981
982  while (p != NULL)
983  {
984    l_p = p->exp[r->VarL_Offset[0]];
985    if (l_p > l_max ||
986        (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
987      l_max = p_GetMaxExpL2(l_max, l_p, r);
988    for (i=1; i<r->VarL_Size; i++)
989    {
990      l_p = p->exp[r->VarL_Offset[i]];
991      // do the divisibility trick to find out whether l has an exponent
992      if (l_p > l_max ||
993          (((l_max & divmask) ^ (l_p & divmask)) != ((l_max-l_p) & divmask)))
994        l_max = p_GetMaxExpL2(l_max, l_p, r);
995    }
996    pIter(p);
997  }
998  return l_max;
999}
1000
1001
1002
1003
1004/***************************************************************
1005 *
1006 * Misc things
1007 *
1008 ***************************************************************/
1009// returns TRUE, if all monoms have the same component
1010BOOLEAN p_OneComp(poly p, const ring r)
1011{
1012  if(p!=NULL)
1013  {
1014    long i = p_GetComp(p, r);
1015    while (pNext(p)!=NULL)
1016    {
1017      pIter(p);
1018      if(i != p_GetComp(p, r)) return FALSE;
1019    }
1020  }
1021  return TRUE;
1022}
1023
1024/*2
1025*test if a monomial /head term is a pure power
1026*/
1027int p_IsPurePower(const poly p, const ring r)
1028{
1029  int i,k=0;
1030
1031  for (i=r->N;i;i--)
1032  {
1033    if (p_GetExp(p,i, r)!=0)
1034    {
1035      if(k!=0) return 0;
1036      k=i;
1037    }
1038  }
1039  return k;
1040}
1041
1042/*2
1043*test if a polynomial is univariate
1044* return -1 for constant,
1045* 0 for not univariate,s
1046* i if dep. on var(i)
1047*/
1048int p_IsUnivariate(poly p, const ring r)
1049{
1050  int i,k=-1;
1051
1052  while (p!=NULL)
1053  {
1054    for (i=r->N;i;i--)
1055    {
1056      if (p_GetExp(p,i, r)!=0)
1057      {
1058        if((k!=-1)&&(k!=i)) return 0;
1059        k=i;
1060      }
1061    }
1062    pIter(p);
1063  }
1064  return k;
1065}
1066
1067// set entry e[i] to 1 if var(i) occurs in p, ignore var(j) if e[j]>0
1068int  p_GetVariables(poly p, int * e, const ring r)
1069{
1070  int i;
1071  int n=0;
1072  while(p!=NULL)
1073  {
1074    n=0;
1075    for(i=r->N; i>0; i--)
1076    {
1077      if(e[i]==0)
1078      {
1079        if (p_GetExp(p,i,r)>0)
1080        {
1081          e[i]=1;
1082          n++;
1083        }
1084      }
1085      else
1086        n++;
1087    }
1088    if (n==r->N) break;
1089    pIter(p);
1090  }
1091  return n;
1092}
1093
1094
1095/*2
1096* returns a polynomial representing the integer i
1097*/
1098poly p_ISet(int i, const ring r)
1099{
1100  poly rc = NULL;
1101  if (i!=0)
1102  {
1103    rc = p_Init(r);
1104    pSetCoeff0(rc,n_Init(i,r));
1105    if (r->cf->nIsZero(p_GetCoeff(rc,r)))
1106      p_LmDelete(&rc,r);
1107  }
1108  return rc;
1109}
1110
1111/*2
1112* an optimized version of p_ISet for the special case 1
1113*/
1114poly p_One(const ring r)
1115{
1116  poly rc = p_Init(r);
1117  pSetCoeff0(rc,n_Init(1,r));
1118  return rc;
1119}
1120
1121/*2
1122* returns a polynomial representing the number n
1123* destroys n
1124*/
1125poly p_NSet(number n, const ring r)
1126{
1127  if (r->cf->nIsZero(n))
1128  {
1129    r->cf->cfDelete(&n, r);
1130    return NULL;
1131  }
1132  else
1133  {
1134    poly rc = p_Init(r);
1135    pSetCoeff0(rc,n);
1136    return rc;
1137  }
1138}
1139
1140/***************************************************************
1141 *
1142 * p_ShallowDelete
1143 *
1144 ***************************************************************/
1145#undef LINKAGE
1146#define LINKAGE
1147#undef p_Delete
1148#define p_Delete p_ShallowDelete
1149#undef n_Delete
1150#define n_Delete(n, r) ((void)0)
1151
1152#include <kernel/p_Delete__T.cc>
1153
Note: See TracBrowser for help on using the repository browser.