source: git/Singular/polys-impl.cc @ 939ac5f

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