source: git/Singular/polys-impl.h @ c232af

spielwiese
Last change on this file since c232af was c232af, checked in by Olaf Bachmann <obachman@…>, 24 years ago
* omalloc stuff git-svn-id: file:///usr/local/Singular/svn/trunk@4524 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 16.9 KB
Line 
1#ifndef POLYS_IMPL_H
2#define POLYS_IMPL_H
3/****************************************
4*  Computer Algebra System SINGULAR     *
5****************************************/
6/* $Id: polys-impl.h,v 1.46 2000-08-14 12:56:44 obachman Exp $ */
7
8/***************************************************************
9 *
10 * File:       polys-impl.h
11 * Purpose:    low-level definition and declarations for polys
12 *
13 * If you touch anything here, you better know what you are doing.
14 * What is here should not be used directly from other routines -- the
15 * encapsulations in polys.h should be used, instead.
16 *
17 ***************************************************************/
18#include "tok.h"
19#include "structs.h"
20#include <omalloc.h>
21
22/***************************************************************
23 *
24 * definition of the poly structure and its fields
25 *
26 ***************************************************************/
27
28#ifndef NDEBUG
29#define VARS (10)   /*max. number of variables as constant*/
30#else
31#define VARS (0)
32#endif
33union s_exp
34{
35#ifdef HAVE_SHIFTED_EXPONENTS
36   long           e[VARS +1];
37   unsigned long  l[VARS +1];
38#else
39   Exponent_t     e[VARS +1];
40   unsigned long  l[(VARS +1)/(sizeof(long)/sizeof(Exponent_t))];
41#endif
42};
43
44typedef s_exp  monomial;
45typedef Exponent_t* Exponent_pt;
46
47typedef unsigned long Order_t;
48struct  spolyrec
49{
50  poly      next; // next needs to be the first field
51  number    coef; // and coef the second --- do not change this !!!
52  monomial  exp; // make sure that exp is aligned
53};
54
55/***************************************************************
56 *
57 * variables/defines used for managment of monomials
58 *
59 ***************************************************************/
60
61#define POLYSIZE (sizeof(poly) + sizeof(number))
62#define POLYSIZEW (POLYSIZE / sizeof(long))
63#define MAX_EXPONENT_NUMBER ((MAX_BLOCK_SIZE - POLYSIZE) / SIZEOF_EXPONENT)
64
65// number of Variables
66extern int pVariables;
67// size of a monom in bytes - always a multiple of sizeof(void*)
68extern int pMonomSize;
69// size of a monom in units of sizeof(void*) -- i.e. in words
70extern int pMonomSizeW;
71// Ceiling((pVariables+1) / sizeof(void*)) == length of exp-vector in words
72// extern int pVariables1W;
73// Ceiling((pVariables) / sizeof(void*))
74// extern int pVariablesW;
75extern int *pVarOffset;
76// extern int pVarLowIndex;
77// extern int pVarHighIndex;
78// extern int pVarCompIndex;
79extern omBin currPolyBin;
80
81/***************************************************************
82 *
83 * Primitives for determening/setting  the way exponents are arranged
84 *
85 ***************************************************************/
86#define _pExpIndex(i) (currRing->VarOffset[(i)])
87#define _pRingExpIndex(r, i)  (r)->VarOffset[(i)]
88
89#define _pCompIndex        (currRing->pCompIndex)
90#define _pRingCompIndex(r)  ((r)->pCompIndex)
91
92/***************************************************************
93 *
94 * Primitives for accessing and setting fields of a poly
95 *
96 ***************************************************************/
97#define _pNext(p)           ((p)->next)
98#define _pIter(p)           ((p) = (p)->next)
99
100#define _pGetCoeff(p)       ((p)->coef)
101#define _pSetCoeff(p,n)     {nDelete(&((p)->coef));(p)->coef=n;}
102#define _pSetCoeff0(p,n)    (p)->coef=n
103
104#define _pGetOrder(p)       ((p)->exp.l[currRing->pOrdIndex])
105
106#if defined(PDEBUG) && PDEBUG > 1
107extern Exponent_t pPDSetExp(poly p, int v, Exponent_t e, char* f, int l);
108extern Exponent_t pPDGetExp(poly p, int v, char* f, int l);
109extern Exponent_t pPDIncrExp(poly p, int v, char* f, int l);
110extern Exponent_t pPDDecrExp(poly p, int v, char* f, int l);
111extern Exponent_t pPDAddExp(poly p, int v, Exponent_t e, char* f, int l);
112extern Exponent_t pPDMultExp(poly p, int v, Exponent_t e, char* f, int l);
113extern Exponent_t pPDSubExp(poly p, int v, Exponent_t e, char* f, int l);
114
115extern Exponent_t pPDRingSetExp(ring r,poly p,int v,Exponent_t e,char* f,int l);
116extern Exponent_t pPDRingGetExp(ring r,poly p, int v, char* f, int l);
117
118extern Exponent_t pDBSetComp(poly p, Exponent_t k, int l, char* f, int l);
119extern Exponent_t pDBDecrComp(poly p, char* f, int l);
120extern Exponent_t pDBAddComp(poly p, Exponent_t k, int l, char* f, int l);
121extern Exponent_t pDBSubComp(poly p, Exponent_t k, char* f, int l);
122extern Exponent_t pDBRingSetComp(ring r, poly p, Exponent_t k, char* f, int l);
123
124
125#define _pSetExp(p,v,e)     pPDSetExp(p,v,e,__FILE__,__LINE__)
126#define _pGetExp(p,v)       pPDGetExp(p,v,__FILE__,__LINE__)
127#define _pIncrExp(p,v)      pPDIncrExp(p,v,__FILE__,__LINE__)
128#define _pDecrExp(p,v)      pPDDecrExp(p,v,__FILE__,__LINE__)
129#define _pAddExp(p,i,v)     pPDAddExp(p,i,v,__FILE__,__LINE__)
130#define _pSubExp(p,i,v)     pPDSubExp(p,i,v,__FILE__,__LINE__)
131#define _pMultExp(p,i,v)    pPDMultExp(p,i,v,__FILE__,__LINE__)
132
133#define _pRingSetExp(r,p,v,e)     pPDRingSetExp(r,p,v,e,__FILE__,__LINE__)
134#define _pRingGetExp(r,p,v)       pPDRingGetExp(r,p,v,__FILE__,__LINE__)
135
136#define _pSetComp(p,k)      pDBSetComp(p, k, 0, __FILE__, __LINE__)
137#define _pDecrComp(p)       pDBDecrComp(p, __FILE__, __LINE__)
138#define _pAddComp(p,v)      pDBAddComp(p,v, 0, __FILE__, __LINE__)
139#define _pSubComp(p,v)      pDBSubComp(p,v, __FILE__, __LINE__)
140#define _pRingSetComp(r,p,k)  pDBRingSetComp(r, p, k, __FILE__, __LINE__)
141
142#define pSetCompS(p, k, l) pDBSetComp(p, k, l, __FILE__, __LINE__)
143
144#else  // ! (defined(PDEBUG) && PDEBUG > 1)
145#ifndef HAVE_SHIFTED_EXPONENTS
146
147#define _pSetExp(p,v,E)     (p)->exp.e[_pExpIndex(v)]=(E)
148#define _pGetExp(p,v)       (p)->exp.e[_pExpIndex(v)]
149#define _pIncrExp(p,v)      ((p)->exp.e[_pExpIndex(v)])++
150#define _pDecrExp(p,v)      ((p)->exp.e[_pExpIndex(v)])--
151#define _pAddExp(p,i,v)     ((p)->exp.e[_pExpIndex(i)]) += (v)
152#define _pSubExp(p,i,v)     ((p)->exp.e[_pExpIndex(i)]) -= (v)
153#define _pMultExp(p,i,v)    ((p)->exp.e[_pExpIndex(i)]) *= (v)
154
155#define _pRingSetExp(r,p,e_i,e_v)     (p)->exp.e[_pRingExpIndex(r,e_i)]=(e_v)
156#define _pRingGetExp(r,p,v)       (p)->exp.e[_pRingExpIndex(r,v)]
157#else
158inline Exponent_t pSGetExp(poly p, int v, ring r)
159{
160  assume(v>0);
161  assume(v<=r->N);
162  return (p->exp.l[(r->VarOffset[v] & 0xffffff)] >> (r->VarOffset[v] >> 24))
163          & r->bitmask;
164}
165
166inline Exponent_t pSSetExp(poly p, int v, int e, ring r)
167{
168  assume(e>=0);
169  assume(v>0);
170  assume(v<=r->N);
171  assume(e<=((int)r->bitmask));
172
173  // shift e to the left:
174  register int shift = r->VarOffset[v] >> 24;
175  unsigned long ee = ((unsigned long)e) << shift /*(r->VarOffset[v] >> 24)*/;
176  // clear the bits in the exponent vector:
177  register int offset = (r->VarOffset[v] & 0xffffff);
178  p->exp.l[offset]  &=
179    ~( r->bitmask << shift );
180  // insert e with |
181  p->exp.l[ offset ] |= ee;
182  return e;
183}
184
185inline Exponent_t pSIncrExp(poly p, int v, ring r)
186{
187  assume(v>0);
188  assume(v<=r->N);
189
190  Exponent_t e=pSGetExp(p,v,r);
191  e++;
192  return pSSetExp(p,v,e,r);
193}
194
195inline Exponent_t pSDecrExp(poly p, int v, ring r)
196{
197  assume(v>0);
198  assume(v<=r->N);
199
200  Exponent_t e=pSGetExp(p,v,r);
201  e--;
202  return pSSetExp(p,v,e,r);
203}
204
205inline Exponent_t pSAddExp(poly p, int v, Exponent_t e, ring r)
206{
207  assume(v>0);
208  assume(v<=r->N);
209
210  Exponent_t ee=pSGetExp(p,v,r);
211  ee+=e;
212  return pSSetExp(p,v,ee,r);
213}
214
215inline Exponent_t pSSubExp(poly p, int v, Exponent_t e, ring r)
216{
217  assume(v>0);
218  assume(v<=r->N);
219
220  Exponent_t ee=pSGetExp(p,v,r);
221  ee-=e;
222  return pSSetExp(p,v,ee,r);
223}
224
225inline Exponent_t pSMultExp(poly p, int v, Exponent_t e, ring r)
226{
227  assume(v>0);
228  assume(v<=r->N);
229
230  Exponent_t ee=pSGetExp(p,v,r);
231  ee*=e;
232  return pSSetExp(p,v,ee,r);
233}
234
235#define _pSetExp(p,v,E)     pSSetExp(p,v,E,currRing)
236#define _pGetExp(p,v)       pSGetExp(p,v,currRing)
237#define _pIncrExp(p,v)      pSIncrExp(p,v,currRing)
238#define _pDecrExp(p,v)      pSDecrExp(p,v,currRing)
239#define _pAddExp(p,i,v)     pSAddExp(p,i,v,currRing)
240#define _pSubExp(p,i,v)     pSSubExp(p,i,v,currRing)
241#define _pMultExp(p,i,v)    pSMultExp(p,i,v,currRing)
242
243#define _pRingSetExp(r,p,e_i,e_v)  pSSetExp(p,e_i,e_v,r)
244#define _pRingGetExp(r,p,v)        pSGetExp(p,v,r)
245#endif
246
247#define _pSetComp(p,k)      _pGetComp(p) = (k)
248#define _pDecrComp(p)       _pGetComp(p)--
249#define _pAddComp(p,v)      _pGetComp(p) += (v)
250#define _pSubComp(p,v)      _pGetComp(p) -= (v)
251#define _pRingSetComp(r,p,k)      (_pRingGetComp(r, p) = (k))
252#define pSetCompS(p, k,l)     _pSetComp(p, k)
253
254#endif // defined(PDEBUG) && PDEBUG > 1
255
256#define _pGetComp(p)        ((p)->exp.e[_pCompIndex])
257#define _pIncrComp(p)       _pGetComp(p)++
258#define _pRingGetComp(r,p)        ((p)->exp.e[_pRingCompIndex(r)])
259
260#ifndef HAVE_SHIFTED_EXPONENTS
261inline Exponent_t _pGetExpSum(poly p1, poly p2, int i)
262{
263  int index = _pExpIndex(i);
264  return p1->exp.e[index] + p2->exp.e[index];
265}
266inline Exponent_t _pGetExpDiff(poly p1, poly p2, int i)
267{
268  int index = _pExpIndex(i);
269  return p1->exp.e[index] - p2->exp.e[index];
270}
271#else
272inline Exponent_t _pGetExpSum(poly p1, poly p2, int i)
273{
274  return _pGetExp(p1,i) + _pGetExp(p2,i);
275}
276inline Exponent_t _pGetExpDiff(poly p1, poly p2, int i)
277{
278  return _pGetExp(p1,i) - _pGetExp(p2,i);
279}
280#endif
281
282inline void _pGetExpV(poly p, Exponent_t *ev)
283{
284  for (int j = pVariables; j; j--)
285      ev[j] = _pGetExp(p, j);
286
287  ev[0] = _pGetComp(p);
288}
289
290extern pSetmProc pSetm;
291inline void _pSetExpV(poly p, Exponent_t *ev)
292{
293  for (int j = pVariables; j; j--)
294      _pSetExp(p, j, ev[j]);
295
296  _pSetComp(p, ev[0]);
297  pSetm(p);
298}
299
300/***************************************************************
301 *
302 * Storage Managament Routines
303 *
304 ***************************************************************/
305#define prAllocMonom(p, r)  omTypeAllocBin(poly, p, r->PolyBin)
306#define _pNew(h)            (poly) omAllocBin(h)
307#define _pInit(h)           (poly) omAlloc0Bin(h)
308#define _pFree1(a, h)       omFreeBin((void*) a, h)
309#define _pRingFree1(r, a)   omFreeBin((void*) a, r->PolyBin)
310
311extern void    _pDelete(poly * a, omBin h);
312extern void    _pDelete1(poly * a, omBin h);
313
314extern poly    _pCopy(poly a);
315extern poly    _pCopy(omBin h, poly a);
316extern poly    _pCopy1(poly a);
317extern poly    _pHead(omBin h, poly a);
318extern poly    _pHead0(poly a);
319extern poly    _pShallowCopyDeleteHead(omBin d_h, poly *s_p, omBin s_h);
320extern poly    _pShallowCopyDelete(omBin d_h, poly *s_p, omBin s_h);
321#define _pCopy2(p1, p2)     omMemcpyW(p1, p2, pMonomSizeW)
322
323
324/***************************************************************
325 *
326 * Routines which work on vectors instead of single exponents
327 *
328 ***************************************************************/
329// Here is a handy Macro which disables inlining when run with
330// profiling and enables it otherwise
331
332#ifdef DO_DEEP_PROFILE
333
334#ifndef POLYS_IMPL_CC
335
336#define DECLARE(type, arglist) type arglist; \
337   static type dummy_##arglist
338#else
339#define DECLARE(type, arglist) type arglist
340#endif // POLYS_IMPL_CC
341
342#else //! DO_DEEP_PROFILE
343
344#define DECLARE(type, arglist ) inline type arglist
345
346#endif // DO_DEEP_PROFILE
347
348
349#if defined(PDEBUG) && PDEBUG > 1
350#define _pMonAddOn(p1, p2)  pDBMonAddOn(p1, p2, __FILE__, __LINE__)
351extern  void pDBMonAddOn(poly p1, poly p2, char* f, int l);
352inline void __pMonAddOn(poly p1, poly p2)
353#else
354  DECLARE(void, _pMonAddOn(poly p1, poly p2))
355#endif // defined(PDEBUG) && PDEBUG > 1
356{
357  int i = currRing->ExpLSize;
358  unsigned long* s1 = &(p1->exp.l[0]);
359  const unsigned long* s2 = &(p2->exp.l[0]);
360  for (;;)
361  {
362    *s1 += *s2;
363    i--;
364    if (i==0) break;
365    s1++;
366    s2++;
367  }
368}
369
370#if defined(PDEBUG) && PDEBUG > 1
371#define _pMonSubFrom(p1, p2)  pDBMonSubFrom(p1, p2, __FILE__, __LINE__)
372extern  void pDBMonSubFrom(poly p1, poly p2, char* f, int l);
373inline void __pMonSubFrom(poly p1, poly p2)
374#else
375  DECLARE(void, _pMonSubFrom(poly p1, poly p2))
376#endif // defined(PDEBUG) && PDEBUG > 1
377{
378  int i = currRing->ExpLSize;
379  unsigned long* s1 = &(p1->exp.l[0]);
380  const unsigned long* s2 = &(p2->exp.l[0]);
381
382  for (;;)
383  {
384    *s1 -= *s2;
385    i--;
386    if (i==0) break;
387    s1++;
388    s2++;
389  }
390}
391
392// Makes p1 a copy of p2 and adds on exponents of p3
393#if defined(PDEBUG) && PDEBUG > 1
394#define _prMonAdd(p1, p2, p3)  pDBMonAdd(p1, p2, p3, __FILE__, __LINE__)
395extern  void prDBMonAdd(poly p1, poly p2, poly p3, char* f, int l);
396inline void __prMonAdd(poly p1, poly p2, poly p3)
397#else
398  DECLARE(void, _prMonAdd(poly p1, poly p2, poly p3, ring r))
399#endif // defined(PDEBUG) && PDEBUG > 1
400{
401  unsigned long* s1 = &(p1->exp.l[0]);
402  const unsigned long* s2 = &(p2->exp.l[0]);
403  const unsigned long* s3 = &(p3->exp.l[0]);
404  const unsigned long* const ub = s3 + r->ExpLSize;
405
406  for (;;)
407  {
408    *s1 = *s2 + *s3;
409    s3++;
410    if (s3 == ub) break;
411    s1++;
412    s2++;
413  }
414}
415
416
417#if SIZEOF_LONG == 4
418
419#if SIZEOF_EXPONENT == 1
420#define P_DIV_MASK 0x80808080
421#define EXPONENT_MAX     0x7f
422#else // SIZEOF_EXPONENT == 2
423#define P_DIV_MASK 0x80008000
424#define EXPONENT_MAX   0x7fff
425#endif
426
427#else // SIZEOF_LONG == 8
428
429#if SIZEOF_EXPONENT == 1
430#define P_DIV_MASK 0x8080808080808080
431#define EXPONENT_MAX             0x7f
432#elif  SIZEOF_EXPONENT == 2
433#define P_DIV_MASK 0x8000800080008000
434#define EXPONENT_MAX           0x7fff
435#else // SIZEOF_EXPONENT == 4
436#define P_DIV_MASK 0x8000000080000000
437#define EXPONENT_MAX       0x7fffffff
438#endif
439
440#endif
441
442// #define LONG_MONOMS
443
444#ifdef LONG_MONOMS
445DECLARE(BOOLEAN, __pDivisibleBy(poly a, poly b))
446{
447  const unsigned long* const lb = (unsigned long*) &(a->exp.l[currRing->pDivLow]);
448  const unsigned long* s1 = (unsigned long*) &(a->exp.l[currRing->pDivHigh]);
449  const unsigned long* s2 = (unsigned long*) &(b->exp.l[currRing->pDivHigh]);
450
451  for (;;)
452  {
453    // Yes, the following is correct, provided that the exponents do
454    // not have their first bit set
455    if ((*s2 - *s1) & P_DIV_MASK) return FALSE;
456    if (s1 == lb) return TRUE;
457    s1--;
458    s2--;
459  }
460}
461#else
462DECLARE(BOOLEAN, __pDivisibleBy(poly a, poly b))
463{
464  int i=pVariables; // assume i>0
465
466  for (;;)
467  {
468    if (_pGetExp(a,i) > _pGetExp(b,i))
469      return FALSE;
470    i--;
471    if (i==0)
472      return TRUE;
473  }
474}
475#endif
476
477#if defined(PDEBUG) && PDEBUG > 1
478#define _pDivisibleBy(a,b)   pDBDivisibleBy(a, b, __FILE__, __LINE__)
479extern  BOOLEAN pDBDivisibleBy(poly p1, poly p2, char* f, int l);
480inline BOOLEAN _pDivisibleBy_orig(poly a, poly b)
481#else
482inline BOOLEAN _pDivisibleBy(poly a, poly b)
483#endif // defined(PDEBUG) && PDEBUG > 1
484{
485  if ((a!=NULL)&&((_pGetComp(a)==0) || (_pGetComp(a) == _pGetComp(b))))
486  {
487    return __pDivisibleBy(a,b);
488  }
489  return FALSE;
490}
491
492#if defined(PDEBUG) && PDEBUG > 1
493#define _pDivisibleBy1(a,b)   pDBDivisibleBy1(a, b, __FILE__, __LINE__)
494extern  BOOLEAN pDBDivisibleBy1(poly p1, poly p2, char* f, int l);
495inline BOOLEAN _pDivisibleBy1_orig(poly a, poly b)
496#else
497inline BOOLEAN _pDivisibleBy1(poly a, poly b)
498#endif // defined(PDEBUG) && PDEBUG > 1
499{
500  if (_pGetComp(a) == 0 || _pGetComp(a) == _pGetComp(b))
501    return __pDivisibleBy(a,b);
502  return FALSE;
503}
504
505#if defined(PDEBUG) && PDEBUG > 1
506#define _pDivisibleBy2(a,b)   pDBDivisibleBy2(a, b, __FILE__, __LINE__)
507extern  BOOLEAN pDBDivisibleBy2(poly p1, poly p2, char* f, int l);
508#else
509#define _pDivisibleBy2(a,b) __pDivisibleBy(a,b)
510#endif // defined(PDEBUG) && PDEBUG > 1
511
512#ifdef PDEBUG
513#define _pEqual(p, q)       mmDBEqual(p, q, __FILE__, __LINE__)
514BOOLEAN mmDBEqual(poly p, poly q, char* file, int line);
515#else
516#define _pEqual __pEqual
517#endif
518
519DECLARE(BOOLEAN, __pEqual(poly p1, poly p2))
520{
521  const unsigned long *s1 = (unsigned long*) &(p1->exp.l[0]);
522  const unsigned long *s2 = (unsigned long*) &(p2->exp.l[0]);
523  const unsigned long* const lb = s1 + currRing->ExpLSize;
524
525  for(;;)
526  {
527    if (*s1 != *s2) return FALSE;
528    s1++;
529    if (s1 == lb) return TRUE;
530    s2++;
531  }
532}
533
534/***************************************************************
535 *
536 * Misc. things
537 *
538 *
539 ***************************************************************/
540// Divisiblity tests based on Short Exponent Vectors
541// define to enable debugging of this
542// #define PDIV_DEBUG
543#if defined(PDEBUG) && ! defined(PDIV_DEBUG)
544#define PDIV_DEBUG
545#endif
546
547#ifdef PDIV_DEBUG
548#define _pShortDivisibleBy(a, sev_a, b, not_sev_b) \
549  pDBShortDivisibleBy(a, sev_a, b, not_sev_b, __FILE__, __LINE__)
550BOOLEAN pDBShortDivisibleBy(poly p1, unsigned long sev_1,
551                            poly p2, unsigned long not_sev_2,
552                            char* f, int l);
553#else
554#define _pShortDivisibleBy(a, sev_a, b, not_sev_b) \
555  ( ! (sev_a & not_sev_b) && pDivisibleBy(a, b))
556#endif
557
558
559/***************************************************************
560 *
561 * Routines which implement low-level manipulations/operations
562 * on exponents and "are allowed" to access single exponetns
563 *
564 ***************************************************************/
565
566#ifdef LONG_MONOMS
567DECLARE(int, __pExpQuerSum2(poly p, int from, int to))
568{
569  int j = 0;
570  int i = from ;
571
572  for(;;)
573  {
574    if (i > to) break;
575    j += p->exp.e[i];
576    i++;
577  }
578  if (from <= _pCompIndex && to >= _pCompIndex)
579    return j - _pGetComp(p);
580  return j;
581}
582
583#define _pExpQuerSum(p)  __pExpQuerSum2(p, currRing->pVarLowIndex, currRing->pVarHighIndex)
584
585inline int _pExpQuerSum2(poly p,int from,int to)
586{
587  int ei_to = _pExpIndex(to);
588  int ei_from = _pExpIndex(from);
589
590  if (ei_from > ei_to)
591    return __pExpQuerSum2(p, ei_to, ei_from);
592  else
593    return __pExpQuerSum2(p, ei_from, ei_to);
594}
595
596#else
597DECLARE(int, _pExpQuerSum(poly p))
598{
599  int s = 0;
600  int i = pVariables;
601  for (;;)
602  {
603    s += _pGetExp(p, i);
604    i--;
605    if (i==0) return s;
606  }
607}
608#endif
609
610#endif // POLYS_IMPL_H
Note: See TracBrowser for help on using the repository browser.