source: git/Singular/polys-impl.h @ 5e4527

spielwiese
Last change on this file since 5e4527 was 5e4527, checked in by Hans Schönemann <hannes@…>, 25 years ago
*hannes: C++-fixes git-svn-id: file:///usr/local/Singular/svn/trunk@3102 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 16.1 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.30 1999-06-08 07:50:58 Singular 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#include "mmheap.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
33
34typedef Exponent_t  monomial[VARS + 1];
35typedef Exponent_t* Exponent_pt;
36
37typedef long Order_t;
38struct  spolyrec
39{
40  poly      next; // next needs to be the first field
41  number    coef; // and coef the second --- do not change this !!!
42  Order_t   Order;
43  monomial  exp; // make sure that exp is aligned
44};
45
46/***************************************************************
47 *
48 * variables/defines used for managment of monomials
49 *
50 ***************************************************************/
51
52#define POLYSIZE (sizeof(poly) + sizeof(number) + sizeof(Order_t))
53#define POLYSIZEW (POLYSIZE / sizeof(long))
54#define MAX_EXPONENT_NUMBER ((MAX_BLOCK_SIZE - POLYSIZE) / SIZEOF_EXPONENT)
55
56// number of Variables
57extern int pVariables;
58// size of a monom in bytes - always a multiple of sizeof(void*)
59extern int pMonomSize;
60// size of a monom in units of sizeof(void*) -- i.e. in words
61extern int pMonomSizeW;
62// Ceiling((pVariables+1) / sizeof(void*)) == length of exp-vector in words
63extern int pVariables1W;
64// Ceiling((pVariables) / sizeof(void*))
65extern int pVariablesW;
66extern int *pVarOffset;
67extern int pVarLowIndex;
68extern int pVarHighIndex;
69extern int pVarCompIndex;
70
71/***************************************************************
72 *
73 * Primitives for determening/setting  the way exponents are arranged
74 *
75 ***************************************************************/
76#define _pExpIndex(i) pVarOffset[(i)]
77#define _pRingExpIndex(r, i)  (r)->VarOffset[(i)]
78
79#define _pCompIndex        pVarOffset[0]
80#define _pRingCompIndex(r)  ((r)->VarOffset[0])
81
82// for simple, lex orderings
83extern void pGetVarIndicies_Lex(int nvars, int* VarOffset, int &VarCompIndex,
84                                int &VarLowIndex, int &VarHighIndex);
85// for simple, revlex orderings
86extern void pGetVarIndicies_RevLex(int nvars,int *VarOffset,int &VarCompIndex,
87                                   int &VarLowIndex, int &VarHighIndex);
88// for all non-simple orderings
89extern void pGetVarIndicies(int nvars, int *VarOffset, int &VarCompIndex,
90                            int &VarLowIndex, int &VarHighIndex);
91// gets var indicies w.r.t. the ring r --
92// determines which one of three pGetVarIndicies((int nvars, ...) to use
93// based on properties of r
94extern void pGetVarIndicies(ring r, int *VarOffset, int &VarCompIndex,
95                            int &VarLowIndex, int &VarHighIndex);
96
97/***************************************************************
98 *
99 * Primitives for accessing and setting fields of a poly
100 *
101 ***************************************************************/
102#define _pNext(p)           ((p)->next)
103#define _pIter(p)           ((p) = (p)->next)
104
105#define _pGetCoeff(p)       ((p)->coef)
106#define _pSetCoeff(p,n)     {nDelete(&((p)->coef));(p)->coef=n;}
107#define _pSetCoeff0(p,n)    (p)->coef=n
108
109#define _pGetOrder(p)       ((p)->Order)
110
111#if defined(PDEBUG) && PDEBUG != 0
112extern Exponent_t pPDSetExp(poly p, int v, Exponent_t e, char* f, int l);
113extern Exponent_t pPDGetExp(poly p, int v, char* f, int l);
114extern Exponent_t pPDIncrExp(poly p, int v, char* f, int l);
115extern Exponent_t pPDDecrExp(poly p, int v, char* f, int l);
116extern Exponent_t pPDAddExp(poly p, int v, Exponent_t e, char* f, int l);
117extern Exponent_t pPDMultExp(poly p, int v, Exponent_t e, char* f, int l);
118extern Exponent_t pPDSubExp(poly p, int v, Exponent_t e, char* f, int l);
119
120extern Exponent_t pPDRingSetExp(ring r,poly p,int v,Exponent_t e,char* f,int l);
121extern Exponent_t pPDRingGetExp(ring r,poly p, int v, char* f, int l);
122
123#define _pSetExp(p,v,e)     pPDSetExp(p,v,e,__FILE__,__LINE__)
124#define _pGetExp(p,v)       pPDGetExp(p,v,__FILE__,__LINE__)
125#define _pIncrExp(p,v)      pPDIncrExp(p,v,__FILE__,__LINE__)
126#define _pDecrExp(p,v)      pPDDecrExp(p,v,__FILE__,__LINE__)
127#define _pAddExp(p,i,v)     pPDAddExp(p,i,v,__FILE__,__LINE__)
128#define _pSubExp(p,i,v)     pPDSubExp(p,i,v,__FILE__,__LINE__)
129#define _pMultExp(p,i,v)    pPDMultExp(p,i,v,__FILE__,__LINE__)
130
131#define _pRingSetExp(r,p,v,e)     pPDRingSetExp(r,p,v,e,__FILE__,__LINE__)
132#define _pRingGetExp(r,p,v)       pPDRingGetExp(r,p,v,__FILE__,__LINE__)
133
134#else  // ! (defined(PDEBUG) && PDEBUG != 0)
135
136#define _pSetExp(p,v,e)     (p)->exp[_pExpIndex(v)]=(e)
137#define _pGetExp(p,v)       (p)->exp[_pExpIndex(v)]
138#define _pIncrExp(p,v)      ((p)->exp[_pExpIndex(v)])++
139#define _pDecrExp(p,v)      ((p)->exp[_pExpIndex(v)])--
140#define _pAddExp(p,i,v)     ((p)->exp[_pExpIndex(i)]) += (v)
141#define _pSubExp(p,i,v)     ((p)->exp[_pExpIndex(i)]) -= (v)
142#define _pMultExp(p,i,v)    ((p)->exp[_pExpIndex(i)]) *= (v)
143
144#define _pRingSetExp(r,p,v,e)     (p)->exp[_pRingExpIndex(r,v)]=(e)
145#define _pRingGetExp(r,p,v)       (p)->exp[_pRingExpIndex(r,v)]
146
147#endif // defined(PDEBUG) && PDEBUG != 0
148
149inline Exponent_t _pGetExpSum(poly p1, poly p2, int i)
150{
151  int index = _pExpIndex(i);
152  return p1->exp[index] + p2->exp[index];
153}
154inline Exponent_t _pGetExpDiff(poly p1, poly p2, int i)
155{
156  int index = _pExpIndex(i);
157  return p1->exp[index] - p2->exp[index];
158}
159
160#define _pGetComp(p)        ((p)->exp[_pCompIndex])
161#define _pSetComp(p,k)      _pGetComp(p) = (k)
162#define _pIncrComp(p)       _pGetComp(p)++
163#define _pDecrComp(p)       _pGetComp(p)--
164#define _pAddComp(p,v)      _pGetComp(p) += (v)
165#define _pSubComp(p,v)      _pGetComp(p) -= (v)
166
167#define _pRingGetComp(r,p)        ((p)->exp[_pRingCompIndex(r)])
168#define _pRingSetComp(r,p,k)      (_pRingGetComp(p) = (k))
169
170inline void _pGetExpV(poly p, Exponent_t *ev)
171{
172  for (int j = pVariables; j; j--)
173      ev[j] = _pGetExp(p, j);
174
175  ev[0] = _pGetComp(p);
176}
177
178extern pSetmProc pSetm;
179inline void _pSetExpV(poly p, Exponent_t *ev)
180{
181  for (int j = pVariables; j; j--)
182      _pSetExp(p, j, ev[j]);
183
184  _pSetComp(p, ev[0]);
185  pSetm(p);
186}
187
188/***************************************************************
189 *
190 * Storage Managament Routines
191 *
192 ***************************************************************/
193#ifdef MDEBUG
194
195poly    pDBInit(char *f, int l);
196poly    pDBCopy(poly a, char *f, int l);
197poly    pDBCopy1(poly a, char *f, int l);
198poly    pDBHead(poly a, char *f, int l);
199poly    pDBHead0(poly a, char *f, int l);
200poly    pDBFetchCopy(ring r, poly a, char *f, int l);
201
202void    pDBDelete(poly * a, char * f, int l);
203void    pDBDelete1(poly * a, char * f, int l);
204
205#define pDBNew(f,l)  (poly) mmDBAllocHeap(mm_specHeap, f,l)
206#define _pNew()      (poly) mmDBAllocHeap(mm_specHeap, __FILE__, __LINE__)
207#define _pInit()     (poly) pDBInit(__FILE__,__LINE__)
208
209#define pDBFree1(a,f,l) mmDBFreeHeap((void*)a, mm_specHeap, f, l)
210
211#define _pDelete(a)     pDBDelete((a),__FILE__,__LINE__)
212#define _pDelete1(a)    pDBDelete1((a),__FILE__,__LINE__)
213
214#define _pCopy(A)       pDBCopy(A,__FILE__,__LINE__)
215#define _pCopy1(A)      pDBCopy1(A, __FILE__,__LINE__)
216#define _pHead(A)       pDBHead(A,__FILE__,__LINE__)
217#define _pHead0(A)      pDBHead0(A, __FILE__,__LINE__)
218#define _pFetchCopy(r,A)    pDBFetchCopy(r, A,__FILE__,__LINE__)
219
220#else // ! MDEBUG
221
222inline poly _pNew()
223{
224  void * p;
225  AllocHeap(p, mm_specHeap);
226  return (poly)p;
227}
228
229inline poly    _pInit(void)
230{
231  poly p = _pNew();
232  memsetW((long*)p,0, pMonomSizeW);
233  return p;
234}
235
236extern void    _pDelete(poly * a);
237extern void    _pDelete1(poly * a);
238
239extern poly    _pCopy(poly a);
240extern poly    _pCopy1(poly a);
241extern poly    _pHead(poly a);
242extern poly    _pHead0(poly a);
243extern poly    _pFetchCopy(ring r,poly a);
244#endif // MDEBUG
245
246#define _pCopy2(p1, p2)     memcpyW(p1, p2, pMonomSizeW)
247#define _pFree1(a)          FreeHeap(a, mm_specHeap)
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 _pMonAddFast(p1, p2)  pDBMonAddFast(p1, p2, __FILE__, __LINE__)
276extern  void pDBMonAddFast(poly p1, poly p2, char* f, int l);
277inline void __pMonAddFast(poly p1, poly p2)
278#else
279  DECLARE(void, _pMonAddFast(poly p1, poly p2))
280#endif // defined(PDEBUG) && PDEBUG == 1
281{
282  // OK -- this might be the only place where we are kind of quick and
283  // dirty: the following only works correctly if all exponents are
284  // positive and the sum of two exponents does not exceed
285  // EXPONENT_MAX
286  Exponent_t c2 = _pGetComp(p2);
287  int i = pVariables1W;
288  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
289  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
290  // set comp of p2 temporarily to 0, so that nothing is added to comp of p1
291  _pSetComp(p2, 0);
292
293  for (;;)
294  {
295    *s1 += *s2;
296    i--;
297    if (i==0) break;
298    s1++;
299    s2++;
300  }
301  // reset comp of p2
302  _pSetComp(p2, c2);
303  _pGetOrder(p1) += _pGetOrder(p2);
304}
305
306// Makes p1 a copy of p2 and adds on exponets of p3
307#if defined(PDEBUG) && PDEBUG == 1
308#define _pCopyAddFast(p1, p2, p3)  pDBCopyAddFast(p1, p2, p3, __FILE__, __LINE__)
309extern  void pDBCopyAddFast(poly p1, poly p2, poly p3, char* f, int l);
310inline void __pCopyAddFast(poly p1, poly p2, poly p3)
311#else
312  DECLARE(void, _pCopyAddFast(poly p1, poly p2, poly p3))
313#endif // defined(PDEBUG) && PDEBUG == 1
314{
315  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
316  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
317  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
318  const unsigned long* const ub = s3 + pVariables1W;
319
320  p1->next = p2->next;
321  p1->coef = p2->coef;
322
323  for (;;)
324  {
325    *s1 = *s2 + *s3;
326    s3++;
327    if (s3 == ub) break;
328    s1++;
329    s2++;
330  }
331  // we first are supposed to do a copy from p2 to p1 -- therefore,
332  // component of p1 is set to comp of p2
333  _pSetComp(p1, _pGetComp(p2));
334  _pGetOrder(p1) = _pGetOrder(p2) + _pGetOrder(p3);
335}
336
337// Similar to pCopyAddFast, except that we do not care about the "next" field
338#if defined(PDEBUG) && PDEBUG == 1
339#define _pCopyAddFast0(p1, p2, p3)  pDBCopyAddFast(p1, p2, p3, __FILE__, __LINE__)
340extern  void pDBCopyAddFast(poly p1, poly p2, poly p3, char* f, int l);
341inline void __pCopyAddFast0(poly p1, poly p2, poly p3)
342#else
343  DECLARE(void, _pCopyAddFast0(poly p1, poly p2, poly p3))
344#endif // defined(PDEBUG) && PDEBUG == 1
345{
346  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
347  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
348  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
349  const unsigned long* const ub = s3 + pVariables1W;
350
351  p1->coef = p2->coef;
352
353  for (;;)
354  {
355    *s1 = *s2 + *s3;
356    s3++;
357    if (s3 == ub) break;
358    s1++;
359    s2++;
360  }
361  _pSetComp(p1, _pGetComp(p2));
362  _pGetOrder(p1) = _pGetOrder(p2) + _pGetOrder(p3);
363}
364
365// Similar to pCopyAddFast0, except that we do not recompute the Order,
366// but assume that it is the sum of the Order of p2 and p3
367#if defined(PDEBUG) && PDEBUG == 1
368#define _pCopyAddFastHomog(p1, p2, p3, Order)  \
369   pDBCopyAddFastHomog(p1, p2, p3, Order, __FILE__, __LINE__)
370extern  void pDBCopyAddFastHomog(poly p1, poly p2, poly p3, Order_t Order,
371                                 char* f, int l);
372inline void __pCopyAddFastHomog(poly p1, poly p2, poly p3, Order_t Order)
373#else
374  DECLARE(void, _pCopyAddFastHomog(poly p1, poly p2, poly p3, Order_t Order))
375#endif // defined(PDEBUG) && PDEBUG == 1
376{
377  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
378  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
379  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
380  const unsigned long* const ub = s3 + pVariables1W;
381
382  p1->coef = p2->coef;
383  p1->Order = Order;
384
385  for (;;)
386  {
387    *s1 = *s2 + *s3;
388    s3++;
389    if (s3 == ub) break;
390    s1++;
391    s2++;
392  }
393  _pSetComp(p1, _pGetComp(p2));
394}
395
396#if SIZEOF_LONG == 4
397
398#if SIZEOF_EXPONENT == 1
399#define P_DIV_MASK 0x80808080
400#else // SIZEOF_EXPONENT == 2
401#define P_DIV_MASK 0x80008000
402#endif
403
404#else // SIZEOF_LONG == 8
405
406#if SIZEOF_EXPONENT == 1
407#define P_DIV_MASK 0x8080808080808080
408#elif  SIZEOF_EXPONENT == 2
409#define P_DIV_MASK 0x8000800080008000
410#else // SIZEOF_EXPONENT == 4
411#define P_DIV_MASK 0x8000000080000000
412#endif
413
414#endif
415
416DECLARE(BOOLEAN, __pDivisibleBy(poly a, poly b))
417{
418#ifdef WORDS_BIGENDIAN
419  const unsigned long* const lb = (unsigned long*) &(a->exp[0]);;
420  const unsigned long* s1 = ((unsigned long*) a) + pMonomSizeW -1;
421  const unsigned long* s2 = ((unsigned long*) b) + pMonomSizeW -1;
422#else
423  const unsigned long* const lb = ((unsigned long*) a) + pMonomSizeW;
424  const unsigned long* s1 = (unsigned long*) &(a->exp[0]);
425  const unsigned long* s2 = (unsigned long*) &(b->exp[0]);
426#endif
427
428  for (;;)
429  {
430    // Yes, the following is correct, provided that the exponents do
431    // not have their first bit set
432    if ((*s2 - *s1) & P_DIV_MASK) return FALSE;
433#ifdef WORDS_BIGENDIAN
434    if (s1 == lb) return TRUE;
435    s1--;
436    s2--;
437#else
438    s1++;
439    if (s1 == lb) return TRUE;
440    s2++;
441#endif
442  }
443}
444
445#if defined(PDEBUG) && PDEBUG == 1
446#define _pDivisibleBy(a,b)   pDBDivisibleBy(a, b, __FILE__, __LINE__)
447extern  BOOLEAN pDBDivisibleBy(poly p1, poly p2, char* f, int l);
448inline BOOLEAN _pDivisibleBy_orig(poly a, poly b)
449#else
450inline BOOLEAN _pDivisibleBy(poly a, poly b)
451#endif // defined(PDEBUG) && PDEBUG == 1
452{
453  if ((a!=NULL)&&((_pGetComp(a)==0) || (_pGetComp(a) == _pGetComp(b))))
454  {
455    return __pDivisibleBy(a,b);
456  }
457  return FALSE;
458}
459
460#if defined(PDEBUG) && PDEBUG == 1
461#define _pDivisibleBy1(a,b)   pDBDivisibleBy1(a, b, __FILE__, __LINE__)
462extern  BOOLEAN pDBDivisibleBy1(poly p1, poly p2, char* f, int l);
463inline BOOLEAN _pDivisibleBy1_orig(poly a, poly b)
464#else
465inline BOOLEAN _pDivisibleBy1(poly a, poly b)
466#endif // defined(PDEBUG) && PDEBUG == 1
467{
468  if (_pGetComp(a) == 0 || _pGetComp(a) == _pGetComp(b))
469    return __pDivisibleBy(a,b);
470  return FALSE;
471}
472
473#if defined(PDEBUG) && PDEBUG == 1
474#define _pDivisibleBy2(a,b)   pDBDivisibleBy2(a, b, __FILE__, __LINE__)
475extern  BOOLEAN pDBDivisibleBy2(poly p1, poly p2, char* f, int l);
476#else
477#define _pDivisibleBy2(a,b) __pDivisibleBy(a,b)
478#endif // defined(PDEBUG) && PDEBUG == 1
479
480
481DECLARE(BOOLEAN, _pEqual(poly p1, poly p2))
482{
483  const long *s1 = (long*) &(p1->exp[0]);
484  const long *s2 = (long*) &(p2->exp[0]);
485  const long* const lb = s1 + pVariables1W;
486
487  for(;;)
488  {
489    if (*s1 != *s2) return FALSE;
490    s1++;
491    if (s1 == lb) return TRUE;
492    s2++;
493  }
494}
495
496/***************************************************************
497 *
498 * Routines which implement low-level manipulations/operations
499 * on exponents and "are allowed" to access single exponetns
500 *
501 ***************************************************************/
502
503DECLARE(int, __pExpQuerSum2(poly p, int from, int to))
504{
505  int j = p->exp[from];
506  int i = from + 1;
507
508  for(;;)
509  {
510    if (i > to) return j;
511    j += p->exp[i];
512    i++;
513  }
514}
515
516#define _pExpQuerSum(p)  __pExpQuerSum2(p, pVarLowIndex, pVarHighIndex)
517
518inline int _pExpQuerSum1(poly p, int to)
519{
520  int ei_to = _pExpIndex(to);
521  int ei_1 = _pExpIndex(1);
522
523  if (ei_1 > ei_to)
524    return __pExpQuerSum2(p, ei_to, ei_1);
525  else
526    return __pExpQuerSum2(p, ei_1, ei_to);
527}
528
529
530inline int _pExpQuerSum2(poly p,int from,int to)
531{
532  int ei_to = _pExpIndex(to);
533  int ei_from = _pExpIndex(from);
534
535  if (ei_from > ei_to)
536    return __pExpQuerSum2(p, ei_to, ei_from);
537  else
538    return __pExpQuerSum2(p, ei_from, ei_to);
539}
540
541#endif // POLYS_IMPL_H
Note: See TracBrowser for help on using the repository browser.