source: git/Singular/polys-impl.cc @ 64eef3

fieker-DuValspielwiese
Last change on this file since 64eef3 was ec7aac, checked in by Olaf Bachmann <obachman@…>, 24 years ago
* replaced prProcs by fast poly procs * fixed various memory leaks * added dError stuff git-svn-id: file:///usr/local/Singular/svn/trunk@4565 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 20.6 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: polys-impl.cc,v 1.49 2000-09-04 13:39:03 obachman Exp $ */
5
6/***************************************************************
7 *
8 * File:       polys-impl.cc
9 * Purpose:    low-level procuders for polys.
10 *
11 * If you touch anything here, you better know what you are doing.
12 * What is here should not be used directly from other routines -- the
13 * encapsulations in polys.[h,cc] should be used, instead.
14 *
15 ***************************************************************/
16#ifndef POLYS_IMPL_CC
17#define POLYS_IMPL_CC
18
19#include <stdio.h>
20#include <string.h>
21#include "mod2.h"
22#include <omalloc.h>
23#include "tok.h"
24#include "structs.h"
25#include "febase.h"
26#include "numbers.h"
27#include "polys.h"
28#include "ring.h"
29#include "polys-impl.h"
30#include "p_Procs.h"
31#include "dError.h"
32
33#ifdef HAVE_SHIFTED_EXPONENTS
34#ifdef PDEBUG
35int pDBsyzComp=0;
36#endif
37#endif
38
39/***************************************************************
40 *
41 * Storage Managament Routines
42 *
43 ***************************************************************/
44
45
46/*2
47* delete a poly, resets pointer
48* put the monomials in the freelist
49*/
50void _pDelete(poly* p, omBin heap)
51{
52  poly h = *p;
53  poly pp;
54
55  while (h!=NULL)
56  {
57    nDelete(&(h->coef));
58    pp=h;
59    pIter(h);
60    _pFree1((ADDRESS)pp, heap);
61  }
62  *p = NULL;
63}
64
65/*2
66* remove first monom
67*/
68void _pDelete1(poly* p, omBin heap)
69{
70  poly h = *p;
71
72  if (h==NULL) return;
73  nDelete(&(h->coef));
74  *p = pNext(h);
75  _pFree1((ADDRESS)h, heap);
76}
77
78
79void ppDelete(poly* p, ring rg)
80{
81  ring origRing = currRing;
82  rChangeCurrRing(rg, FALSE);
83  pDelete(p);
84  rChangeCurrRing(origRing, FALSE);
85}
86
87/*2
88* creates a copy of p
89*/
90poly _pCopy(omBin d_h, poly s_p)
91{
92  spolyrec dp;
93  poly d_p = &dp;
94
95  assume(d_h != NULL && (d_h == currPolyBin) ||
96         d_h->sizeW == currPolyBin->sizeW);
97  pTest(s_p);
98  while (s_p != NULL)
99  {
100    omTypeAllocBin( poly,d_p->next, d_h);
101    d_p = d_p->next;
102    pSetCoeff0(d_p, nCopy(pGetCoeff(s_p)));
103    omMemcpyW(&(d_p->exp.l[0]), &(s_p->exp.l[0]), currRing->ExpLSize);
104    pIter(s_p);
105  }
106  pNext(d_p) = NULL;
107  pHeapTest(dp.next, d_h);
108  return dp.next;
109}
110
111poly _pCopy(poly s_p)
112{
113  return _pCopy(currPolyBin, s_p);
114}
115
116
117poly _pShallowCopyDelete(omBin d_h, poly *p, omBin s_h)
118{
119  spolyrec dp;
120  poly d_p = &dp, tmp;
121  poly s_p = *p;
122
123  assume(d_h != NULL && s_h != NULL &&
124         d_h->sizeW == s_h->sizeW);
125
126  if (currRing->ExpLSize <= 2)
127  {
128    if (currRing->ExpLSize == 1)
129    {
130      while (s_p != NULL)
131      {
132        omTypeAllocBin( poly,d_p->next, d_h);
133        d_p = d_p->next;
134
135        d_p->coef = s_p->coef;
136        d_p->exp.l[0] = s_p->exp.l[0];
137
138        tmp = pNext(s_p);
139        omFreeBin(s_p, s_h);
140        s_p = tmp;
141      }
142    }
143    else
144    {
145      while (s_p != NULL)
146      {
147        omTypeAllocBin( poly,d_p->next, d_h);
148        d_p = d_p->next;
149
150        d_p->coef = s_p->coef;
151        d_p->exp.l[0] = s_p->exp.l[0];
152        d_p->exp.l[1] = s_p->exp.l[1];
153
154        tmp = pNext(s_p);
155        omFreeBin(s_p, s_h);
156        s_p = tmp;
157      }
158    }
159  }
160  else
161  {
162    if (currRing->ExpLSize & 1)
163    {
164      while (s_p != NULL)
165      {
166
167        omTypeAllocBin( poly,d_p->next, d_h);
168        d_p = d_p->next;
169
170        d_p->coef = s_p->coef;
171        omMemcpy_nwODD(&(d_p->exp.l[0]), &(s_p->exp.l[1]), currRing->ExpLSize);
172
173        tmp = pNext(s_p);
174        omFreeBin(s_p, s_h);
175        s_p = tmp;
176      }
177    }
178    else
179    {
180      while (s_p != NULL)
181      {
182
183        omTypeAllocBin(poly,d_p->next, d_h);
184        d_p = d_p->next;
185
186        d_p->coef = s_p->coef;
187        omMemcpy_nwEVEN(&(d_p->exp.l[0]), &(s_p->exp.l[1]), currRing->ExpLSize);
188
189        tmp = pNext(s_p);
190        omFreeBin(s_p, s_h);
191        s_p = tmp;
192      }
193    }
194  }
195  pNext(d_p) = NULL;
196  pHeapTest(dp.next, d_h);
197  *p = NULL;
198  return pNext(dp.next);
199}
200
201
202/*2
203* creates a copy of the initial monomial of p
204* sets the coeff of the copy to a defined value
205*/
206poly _pCopy1(poly p)
207{
208  poly w;
209  w = pNew();
210  pCopy2(w,p);
211  nNew(&(w->coef));
212  pNext(w) = NULL;
213  return w;
214}
215
216/*2
217* returns (a copy of) the head term of a
218*/
219poly _pHead(omBin heap, poly p)
220{
221  poly w=NULL;
222
223  if (p!=NULL)
224  {
225    assume(heap != NULL && (heap == currPolyBin) ||
226           heap->sizeW == currPolyBin->sizeW);
227
228    omTypeAllocBin( poly,w, heap);
229    omMemcpyW(&(w->exp.l[0]), &(p->exp.l[0]), currRing->ExpLSize);
230    pSetCoeff0(w,nCopy(pGetCoeff(p)));
231    pNext(w) = NULL;
232  }
233  return w;
234}
235
236poly _pShallowCopyDeleteHead(omBin d_h, poly *s_p, omBin s_h)
237{
238  poly w = NULL;
239  poly p = *s_p;
240
241  if (p!=NULL)
242  {
243    assume(d_h != NULL && s_h != NULL &&
244           d_h->sizeW == s_h->sizeW);
245
246    omTypeAllocBin( poly,w, d_h);
247    omMemcpyW(&(w->exp.l[0]), &(p->exp.l[0]), currRing->ExpLSize);
248    pSetCoeff0(w,pGetCoeff(p));
249    pNext(w) = NULL;
250
251    *s_p = pNext(p);
252    omFreeBin(p, s_h);
253  }
254  return w;
255}
256
257
258
259poly pHeadProc(poly p)
260{
261  return pHead(p);
262}
263
264/*2
265* returns (a copy of) the head term of a without the coef
266*/
267poly _pHead0(poly p)
268{
269  poly w=NULL;
270
271  if (p!=NULL)
272  {
273    w = pNew();
274    pCopy2(w,p);
275    pSetCoeff0(w,NULL);
276    pNext(w) = NULL;
277  }
278  return w;
279}
280
281
282/***************************************************************
283 *
284 * Routines for turned on debugging
285 *
286 ***************************************************************/
287
288#if defined(PDEBUG) && PDEBUG > 1
289Exponent_t pPDSetExp(poly p, int v, Exponent_t e, char* f, int l)
290{
291  if (v == 0)
292  {
293    dReportError("zero index to exponent in %s:%d\n", f, l);
294  }
295  if (v > pVariables)
296  {
297    dReportError("index %d to exponent too large in %s:%d\n", v, f, l);
298  }
299  return (p)->exp.e[_pExpIndex(v)]=(e);
300}
301
302Exponent_t pPDGetExp(poly p, int v, char* f, int l)
303{
304  if (v == 0)
305  {
306    dReportError("zero index to exponent in %s:%d\n", f, l);
307  }
308  if (v > pVariables)
309  {
310    dReportError("index %d to exponent too large in %s:%d\n", v, f, l);
311  }
312  return (p)->exp.e[_pExpIndex(v)];
313}
314
315Exponent_t pPDRingSetExp(ring r, poly p, int v, Exponent_t e, char* f, int l)
316{
317  if (v == 0)
318  {
319    dReportError("zero index to exponent in %s:%d\n", f, l);
320  }
321  if (v > r->N)
322  {
323    dReportError("index %d to exponent too large in %s:%d\n", v, f, l);
324  }
325  return (p)->exp.e[_pRingExpIndex(r, v)]=(e);
326}
327
328Exponent_t pPDRingGetExp(ring r, poly p, int v, char* f, int l)
329{
330  if (v == 0)
331  {
332    dReportError("zero index to exponent in %s:%d\n", f, l);
333  }
334  if (v > r->N)
335  {
336    dReportError("index %d to exponent too large in %s:%d\n", v, f, l);
337  }
338  return (p)->exp.e[_pRingExpIndex(r,v)];
339}
340
341Exponent_t pPDIncrExp(poly p, int v, char* f, int l)
342{
343  if (v == 0)
344  {
345    dReportError("zero index to exponent in %s:%d\n", f, l);
346  }
347  if (v > pVariables)
348  {
349    dReportError("index %d to exponent too large in %s:%d\n", v, f, l);
350  }
351  return ((p)->exp.e[_pExpIndex(v)])++;
352}
353
354Exponent_t pPDDecrExp(poly p, int v, char* f, int l)
355{
356  if (v == 0)
357  {
358    dReportError("zero index to exponent in %s:%d\n", f, l);
359  }
360  if (v > pVariables)
361  {
362    dReportError("index %d to exponent too large in %s:%d\n", v, f, l);
363  }
364  return ((p)->exp.e[_pExpIndex(v)])--;
365}
366
367Exponent_t pPDAddExp(poly p, int v, Exponent_t e, char* f, int l)
368{
369  if (v == 0)
370  {
371    dReportError("zero index to exponent in %s:%d\n", f, l);
372  }
373  if (v > pVariables)
374  {
375    dReportError("index %d to exponent too large in %s:%d\n", v, f, l);
376  }
377  return ((p)->exp.e[_pExpIndex(v)]) += (e);
378}
379
380Exponent_t pPDSubExp(poly p, int v, Exponent_t e, char* f, int l)
381{
382  if (v == 0)
383  {
384    dReportError("zero index to exponent in %s:%d\n", f, l);
385  }
386  if (v > pVariables)
387  {
388    dReportError("index %d to exponent too large in %s:%d\n", v, f, l);
389  }
390  return ((p)->exp.e[_pExpIndex(v)]) -= (e);
391}
392
393Exponent_t pPDMultExp(poly p, int v, Exponent_t e, char* f, int l)
394{
395  if (v == 0)
396  {
397    dReportError("zero index to exponent in %s:%d\n", f, l);
398  }
399  if (v > pVariables)
400  {
401    dReportError("index %d to exponent too large in %s:%d\n", v, f, l);
402  }
403  return ((p)->exp.e[_pExpIndex(v)]) *= (e);
404}
405
406// routines on components
407Exponent_t pDBSetComp(poly p, Exponent_t k, int le, char* f, int l)
408{
409  if (k < 0)
410  {
411    dReportError("set negative component %d in %s:%d", k, f, l);
412  }
413  if (currRing->order[1] == ringorder_S)
414  {
415    if (le <= 0) le = currRing->typ[1].data.syzcomp.length;
416    if (k > l)
417    {
418      dReportError("component %d larger then max %d in %s:%d",
419            k, le, f, l);
420    }
421  }
422  return _pGetComp(p) = (k);
423}
424
425Exponent_t pDBDecrComp(poly p, char* f, int l)
426{
427  if (_pGetComp(p) < 1)
428  {
429    dReportError("decrement to negative component %d in %s:%d\n", _pGetComp(p), f, l);
430  }
431  return _pGetComp(p)--;
432}
433
434Exponent_t pDBAddComp(poly p, Exponent_t k, int le, char* f, int l)
435{
436  if (_pGetComp(p) + k < 0)
437  {
438    dReportError("add to negative component %d + %d = %d in %s:%d\n", _pGetComp(p),
439          k, _pGetComp(p) + k, f, l);
440  }
441  _pGetComp(p) += (k);
442
443  if (currRing->order[1] == ringorder_S)
444  {
445    if (le <= 0) le = currRing->typ[1].data.syzcomp.length;
446    if (_pGetComp(p) > le)
447    {
448      dReportError("sum of components %d larger then max %d in %s:%d\n",
449            _pGetComp(p), le, f, l);
450      assume(0);
451    }
452  }
453  return _pGetComp(p);
454}
455
456Exponent_t pDBSubComp(poly p, Exponent_t k, char* f, int l)
457{
458  if (_pGetComp(p) - k < 0)
459  {
460    dReportError("sub to negative component %d - %d = %d in %s:%d\n", _pGetComp(p),
461          k, _pGetComp(p) - k, f, l);
462  }
463  return _pGetComp(p) -= (k);
464}
465
466Exponent_t pDBRingSetComp(ring r, poly p, Exponent_t k, char* f, int l)
467{
468  if (k < 0)
469  {
470    dReportError("set negative component %d in %s:%d\n", k, f, l);
471  }
472  return _pRingGetComp(r, p) = (k);
473}
474
475// checks whether fast monom add did not overflow
476void pDBMonAddOn(poly p1, poly p2, char* f, int l)
477{
478  poly ptemp = pHead0(p1);
479
480  if (pGetComp(p1) != 0 && pGetComp(p2) != 0)
481  {
482    dReportError("Error in pMonAddOn: both components %d:%d !=0 in %s:%d",
483         pGetComp(p1), pGetComp(p2), f, l);
484  }
485
486  __pMonAddOn(p1, p2);
487
488  for (int i=1; i<=pVariables; i++)
489  {
490    pAddExp(ptemp, i, pGetExp(p2, i));
491  }
492  pAddComp(ptemp, pGetComp(p2));
493  pSetm(ptemp);
494
495  if (! pEqual(ptemp, p1))
496  {
497    dReportError("Error in pMonAddOn in %s:%d\n", f, l);
498  }
499
500  pFree1(ptemp);
501}
502
503void pDBMonSubFrom(poly p1, poly p2, char* f, int l)
504{
505  poly ptemp = pNew();
506  pCopy2(ptemp, p1);
507
508  if ((pGetComp(p1) != pGetComp(p2)) && (pGetComp(p2)!=0))
509  {
510    dReportError("Error in pMonSubFrom: components are different %d:%d in %s:%d",
511         pGetComp(p1), pGetComp(p2), f, l);
512  }
513
514  __pMonSubFrom(p1, p2);
515
516  for (int i=1; i<=pVariables; i++)
517  {
518    if (pGetExp(ptemp, i) < pGetExp(p2, i))
519    {
520      dReportError("Error in pMonSubFrom: %dth exponent %d of p1 smaller than %d of p2", i, pGetExp(ptemp, i), pGetExp(p2, i));
521    }
522    pSubExp(ptemp, i, pGetExp(p2, i));
523  }
524  pSetComp(ptemp, pGetComp(ptemp)-pGetComp(p2));
525  pSetm(ptemp);
526
527  if (! pEqual(ptemp, p1))
528  {
529    dReportError("Error in pMonSubFrom in %s:%d", f, l);
530  }
531
532  pFree1(ptemp);
533}
534
535void prDBMonAdd(poly p1, poly p2, poly p3, ring r, char* f, int l)
536{
537  if (r == currRing)
538  {
539    if (pGetComp(p3) != 0 && pGetComp(p2) != 0)
540    {
541      dReportError("Error in pMonAdd: both components %d:%d !=0 in %s:%d",
542           pGetComp(p3), pGetComp(p2), f, l);
543    }
544    if (p2 == p1 || p3 == p1)
545    {
546      dReportError("Error in pMonAdd: Destination equals source in %s:%d", f, l);
547    }
548  }
549  __pMonAdd(p1, p2, p3, r);
550
551  if (r == currRing)
552  {
553    poly ptemp = pInit();
554    for (int i=1; i<=pVariables; i++)
555    {
556      pSetExp(ptemp, i, pGetExp(p2, i) + pGetExp(p3, i));
557      if (pGetExp(ptemp, i) != pGetExp(p1, i))
558      {
559        dReportError("Error in pMonAdd: %th exponent: %d != (%d == %d + %d)",
560             i, pGetExp(p1, i), pGetExp(ptemp, i), pGetExp(p2, i),
561             pGetExp(p3, i));
562      }
563    }
564    pSetComp(ptemp, pGetComp(p2) + pGetComp(p3));
565    pSetm(ptemp);
566
567    if (! pEqual(ptemp, p1))
568      dReportError("Error in pMonAdd in %s:%d", f, l);
569    pFree1(ptemp);
570  }
571}
572
573static BOOLEAN OldpDivisibleBy(poly a, poly b)
574{
575  if ((a!=NULL)&&((pGetComp(a)==0) || (pGetComp(a) == pGetComp(b))))
576  {
577    for (int i = 1; i<=pVariables; i++)
578      if (pGetExp(a, i) > pGetExp(b,i)) return FALSE;
579    return TRUE;
580  }
581  return FALSE;
582}
583
584
585BOOLEAN pDBDivisibleBy(poly a, poly b, char* f, int l)
586{
587  BOOLEAN istrue = OldpDivisibleBy(a,b);
588  BOOLEAN f_istrue = _pDivisibleBy_orig(a, b);
589
590  if (istrue != f_istrue)
591  {
592    dReportError("Error in pDivisibleBy in %s:%d\n", f, l);
593    _pDivisibleBy_orig(a, b);
594  }
595  return istrue;
596}
597
598BOOLEAN pDBDivisibleBy1(poly a, poly b, char* f, int l)
599{
600  BOOLEAN istrue = OldpDivisibleBy(a,b);
601  BOOLEAN f_istrue = _pDivisibleBy1_orig(a, b);
602
603  if (istrue != f_istrue)
604  {
605    dReportError("Error in pDivisibleBy1 in %s:%d\n", f, l);
606    _pDivisibleBy1_orig(a, b);
607  }
608  return istrue;
609}
610
611BOOLEAN pDBDivisibleBy2(poly a, poly b, char* f, int l)
612{
613  BOOLEAN istrue = OldpDivisibleBy(a,b);
614  BOOLEAN f_istrue = __pDivisibleBy(a, b);
615
616  if (istrue != f_istrue)
617  {
618    dReportError("Error in pDivisibleBy2 in %s:%d\n", f, l);
619    __pDivisibleBy(a, b);
620  }
621  return f_istrue;
622}
623
624#endif //  defined(PDEBUG) && PDEBUG > 1
625
626#ifdef PDEBUG
627BOOLEAN mmDBEqual(poly p, poly q, char *f, int l)
628{
629  int i;
630
631  for (i = 1; i<=pVariables; i++)
632  {
633    if (pGetExp(p,i) != pGetExp(q, i)) return FALSE;
634  }
635  if (pGetComp(p) != pGetComp(q)) return FALSE;
636  if (__pEqual(p, q) != TRUE)
637  {
638    dReportError("Error in pEqual: exp/comp same, bug monoms different in %s:%d",
639         f, l);
640  }
641  return TRUE;
642}
643
644BOOLEAN prDBTest(poly p, ring r, char* f, int l)
645{
646  ring old_ring = NULL;
647  BOOLEAN res;
648
649  if (r != currRing)
650  {
651    old_ring = currRing;
652    rChangeCurrRing(r);
653  }
654  res = pDBTest(p, currRing->PolyBin, f, l);
655  if (old_ring != NULL)
656  {
657    rChangeCurrRing(old_ring);
658  }
659  return res;
660}
661
662
663BOOLEAN pDBTest(poly p, char *f, int l)
664{
665  return pDBTest(p, currPolyBin, f,l);
666}
667
668BOOLEAN pDBTest(poly p, omBin heap, char *f, int l)
669{
670  poly old=NULL;
671  BOOLEAN ismod=FALSE;
672  if (heap == NULL) heap = currPolyBin;
673
674  while (p!=NULL)
675  {
676    omCheckIf(omCheckAddrBin(p, heap), return FALSE);
677#ifdef LDEBUG
678    if (!nDBTest(p->coef,f,l))
679      return FALSE;
680#endif
681    if ((p->coef==NULL)&&(nGetChar()<2))
682    {
683      dReportError("NULL coef in poly in %s:%d\n",f,l);
684      return FALSE;
685    }
686    if (nIsZero(p->coef))
687    {
688      dReportError("zero coef in poly in %s:%d\n",f,l);
689      return FALSE;
690    }
691    int i=pVariables;
692#ifndef HAVE_SHIFTED_EXPONENTS
693    // can not hapen for SHIFTED_EXPONENTS
694    for(;i;i--)
695    {
696      if (pGetExp(p,i)<0)
697      {
698        dReportError("neg. Exponent %d of x(%d) in %s:%d\n",pGetExp(p,i),i,f,l);
699        return FALSE;
700      }
701    }
702#endif
703    if (pGetComp(p)<0)
704    {
705      dReportError("neg Component in %s:%d\n",f,l);
706      return FALSE;
707    }
708    if (ismod==0)
709    {
710      if (pGetComp(p)>0) ismod=1;
711      else               ismod=2;
712    }
713    else if (ismod==1)
714    {
715      if (pGetComp(p)==0)
716      {
717        dReportError("mix vec./poly in %s:%d\n",f,l);
718        return FALSE;
719      }
720    }
721    else if (ismod==2)
722    {
723      if (pGetComp(p)!=0)
724      {
725        dReportError("mix poly/vec. in %s:%d\n",f,l);
726        return FALSE;
727      }
728    }
729    if (currRing->order[1] == ringorder_S)
730    {
731      long c1, cc1, ccc1, ec1;
732      sro_ord* o = &(currRing->typ[1]);
733
734      c1 = pGetComp(p);
735      cc1 = o->data.syzcomp.Components[c1];
736      ccc1 = o->data.syzcomp.ShiftedComponents[cc1];
737      if (! (c1 == 0 || cc1 != 0))
738      {
739        dReportError("Component <-> TrueComponent zero mismatch\n", f, l);
740        return FALSE;
741      }
742      if (! (c1 == 0 || ccc1 != 0))
743      {
744        dReportError("Component <-> ShiftedComponent zero mismatch\n", f, l);
745        return FALSE;
746      }
747      ec1 = p->exp.l[currRing->typ[1].data.syzcomp.place];
748      if (ec1 != ccc1)
749      {
750        dReportError("Shifted comp out of sync. should %d, is %d", ccc1, ec1);
751        return FALSE;
752      }
753    }
754    if (currRing->order[0] == ringorder_s)
755    {
756      unsigned long syzindex = p->exp.l[currRing->typ[0].data.syz.place];
757      pSetm(p);
758      if (p->exp.l[currRing->typ[0].data.syz.place] != syzindex)
759      {
760        dReportError("Syzindex wrong: Was %dl but should be %d in %s:%d\n",
761             syzindex, p->exp.l[currRing->typ[0].data.syz.place], f, l);
762      }
763    }
764    old=p;
765    pIter(p);
766    if (pComp(old,p)!=1)
767    {
768      dReportError("wrong order (");
769      wrp(old);
770      Print(") in %s:%d (pComp=%d)\n",f,l,pComp(old,p));
771      return FALSE;
772    }
773    if (p != NULL)
774    {
775      if (pGetComp(old) == pGetComp(p))
776      {
777        i = pVariables;
778        for (i=pVariables;i; i--)
779        {
780          if (pGetExp(old, i) != pGetExp(p, i)) break;
781        }
782        if (i == 0)
783        {
784          dReportError("different Compare, but same exponent vector for");
785          wrp(old);
786          Print(" in %s%d\n", f, l);
787          return FALSE;
788        }
789      }
790    }
791  }
792  return TRUE;
793}
794
795#endif // PDEBUG
796
797static unsigned long GetBitFields(Exponent_t e,
798                                  unsigned int s, unsigned int n)
799{
800  unsigned int i = 0, ev = 0;
801  assume(n > 0 && s < BIT_SIZEOF_LONG);
802  do
803  {
804    assume(s+i < BIT_SIZEOF_LONG);
805    if (e > (Exponent_t) i) ev |= Sy_bit(s+i);
806    else break;
807    i++;
808  }
809  while (i < n);
810  return ev;
811}
812
813// Short Exponent Vectors are used for fast divisibility tests
814// ShortExpVectors "squeeze" an exponent vector into one word as follows:
815// Let n = BIT_SIZEOF_LONG / pVariables.
816// If n == 0 (i.e. pVariables > BIT_SIZE_OF_LONG), let m == the number
817// of non-zero exponents. If (m>BIT_SIZEOF_LONG), then sev = ~0, else
818// first m bits of sev are set to 1.
819// Otherwise (i.e. pVariables <= BIT_SIZE_OF_LONG)
820// represented by a bit-field of length n (resp. n+1 for some
821// exponents). If the value of an exponent is greater or equal to n, then
822// all of its respective n bits are set to 1. If the value of an exponent
823// is smaller than n, say m, then only the first m bits of the respective
824// n bits are set to 1, the others are set to 0.
825// This way, we have:
826// exp1 / exp2 ==> (ev1 & ~ev2) == 0, i.e.,
827// if (ev1 & ~ev2) then exp1 does not divide exp2
828unsigned long pGetShortExpVector(poly p)
829{
830  assume(p != NULL);
831  if (p == NULL) return 0;
832  unsigned long ev = 0; // short exponent vector
833  unsigned int n = BIT_SIZEOF_LONG / pVariables; // number of bits per exp
834  unsigned int m1; // highest bit which is filled with (n+1)
835  unsigned int i = 0, j=1;
836
837  if (n == 0)
838  {
839    for (; j<=(unsigned long) pVariables; j++)
840    {
841      if (pGetExp(p,j) > 0) i++;
842      if (i == BIT_SIZEOF_LONG) break;
843    }
844    ev = (unsigned long) ~0 >> ((unsigned long) (BIT_SIZEOF_LONG - i));
845    return ev;
846  }
847  else
848  {
849    m1 = (n+1)*(BIT_SIZEOF_LONG - n*pVariables);
850  }
851
852  n++;
853  while (i<m1)
854  {
855    ev |= GetBitFields(pGetExp(p, j), i, n);
856    i += n;
857    j++;
858  }
859
860  n--;
861  while (i<BIT_SIZEOF_LONG)
862  {
863    ev |= GetBitFields(pGetExp(p, j), i, n);
864    i += n;
865    j++;
866  }
867  return ev;
868}
869
870#ifdef PDIV_DEBUG
871static int pDivisibleBy_number = 1;
872static int pDivisibleBy_FALSE = 1;
873static int pDivisibleBy_ShortFalse = 1;
874static int pDivisibleBy_Null = 1;
875BOOLEAN pDBShortDivisibleBy(poly p1, unsigned long sev_1,
876                            poly p2, unsigned long not_sev_2,
877                            char* f, int l)
878{
879  if (sev_1 != 0 && pGetShortExpVector(p1) != sev_1)
880  {
881    dReportError("sev1 is %o but should be %o in %s:%d\n", sev_1,
882          pGetShortExpVector(p1), f, l);
883    assume(0);
884  }
885  if (~ pGetShortExpVector(p2) != not_sev_2)
886  {
887    dReportError("not_sev2 is %o but should be %o in %s:%d\n", not_sev_2,
888          ~ pGetShortExpVector(p2), f, l);
889    assume(0);
890  }
891  if (sev_1 == 0) pDivisibleBy_Null++;
892  pDivisibleBy_number++;
893  BOOLEAN ret = pDivisibleBy(p1, p2);
894  if (! ret) pDivisibleBy_FALSE++;
895  if (sev_1 & not_sev_2)
896  {
897    pDivisibleBy_ShortFalse++;
898    if (ret)
899    {
900      dReportError("p1 divides p2, but sev's are wrong in %s:%d\n", f, l);
901      assume(0);
902    }
903  }
904  return ret;
905}
906
907void pPrintDivisbleByStat()
908{
909  Print("#Tests: %d; #FALSE %d(%d); #SHORT %d(%d) #NULL:%d(%d)\n",
910        pDivisibleBy_number,
911        pDivisibleBy_FALSE, pDivisibleBy_FALSE*100/pDivisibleBy_number,
912        pDivisibleBy_ShortFalse, pDivisibleBy_ShortFalse*100/pDivisibleBy_FALSE,
913        pDivisibleBy_Null, pDivisibleBy_Null*100/pDivisibleBy_number);
914
915}
916#endif
917
918#ifdef HAVE_SHIFTED_EXPONENTS
919int rComp0_Func(poly p1,poly p2)
920{
921  int i;
922  for(i=0; i<=currRing->pCompHighIndex;i++)
923  {
924    if (p1->exp.l[i] != p2->exp.l[i])
925    {
926      if (p1->exp.l[i] > p2->exp.l[i])
927        return currRing->ordsgn[i];
928      else
929        return -currRing->ordsgn[i];
930    }
931  }
932  return 0;
933}
934
935#ifdef PDEBUG
936int rComp0(poly p1,poly p2)
937{
938  int i;
939  for(i=0; i<=currRing->pCompHighIndex;i++)
940  {
941    if (p1->exp.l[i] != p2->exp.l[i])
942    {
943      if (p1->exp.l[i] > p2->exp.l[i])
944        return currRing->ordsgn[i];
945      else
946        return -currRing->ordsgn[i];
947    }
948  }
949  return 0;
950}
951#endif
952#endif
953#endif // POLYS_IMPL_CC
Note: See TracBrowser for help on using the repository browser.