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

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