source: git/kernel/p_polys.cc @ 2132395

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