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

spielwiese
Last change on this file since c2b529 was c2b529, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: bug fixes to TEST_MAC_ORDER (binom.* polys-impl.* spSpolyLoop.cc) git-svn-id: file:///usr/local/Singular/svn/trunk@1096 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 32.4 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.16 1998-01-27 18:51:21 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 "tok.h"
19#include "structs.h"
20#include "mmemory.h"
21#include "binom.h"
22
23/***************************************************************
24 *
25 * definition of the poly structure and its fields
26 *
27 ***************************************************************/
28
29// EXPONENT_TYPE is determined by configure und defined in mod2.h
30typedef EXPONENT_TYPE Exponent_t;
31
32#define VARS (100)   /*max. number of variables as constant*/
33
34typedef Exponent_t  monomial[VARS + 1];
35typedef Exponent_t* Exponent_pt;
36
37typedef long Order_t;
38// make sure that exp is aligned
39struct  spolyrec
40{
41  poly      next;
42  number    coef;
43  Order_t   Order;
44  #ifdef TEST_MAC_DEBUG
45  Order_t   MOrder;
46  #endif
47  monomial  exp;
48};
49
50/***************************************************************
51 * MACROS CONTROLING MONOMIAL COMPARIONS:
52
53 * COMP_TRADITIONAL
54     Keeps the traditional comparison routines
55     defined -- needed as long as their might be comparisons with
56     negativ components.
57     All the traditional routines are prefixed by t_
58
59 * COMP_FAST
60     Implements monomial operations using the fast vector
61     techniques and several other extensions which go along with that.
62     Undefine in case there are problems.
63
64 * COMP_STATISTIC
65     Provides several routines for accumulating statistics on monomial
66     comparisons and divisibility tests
67
68 * COMP_DEBUG
69     Turns on debugging of COMP_FAST by comparing the results of fast
70     comparison with traditional comparison
71
72 * COMP_NO_EXP_VECTOR_OPS
73    Like COMP_FAST, except that it turns off "vector techniques" of
74    monomial operations, i.e. does everything exponent-wise.
75 ***************************************************************/
76// #define COMP_FAST
77// #define COMP_DEBUG
78// #define COMP_NO_EXP_VECTOR_OPS
79#define COMP_TRADITIONAL
80
81#if defined(COMP_NO_EXP_VECTOR_OPS) && ! defined(COMP_FAST)
82#define COMP_FAST
83#endif
84
85#if defined(COMP_FAST) && ! defined(NDEBUG)
86#define COMP_DEBUG
87#endif
88
89// some relations between these flags
90#ifdef COMP_DEBUG
91#define COMP_TRADITIONAL
92#define COMP_FAST
93#undef  COMP_PROFILE
94#undef  COMP_STATISTICS
95#endif // COMP_DEBUG
96
97#ifdef COMP_STATISTICS
98#undef COMP_FAST
99#endif // COMP_STATISTICS
100
101// for the time being COMP_TRADITIONAL always has to be defined, since
102// traditional routines are needed in spolys.cc -- monomials with
103// negative exponents are compared there!
104#define COMP_TRADITIONAL
105
106/***************************************************************
107 *
108 * variables used for storage management and monomial traversions
109 *
110 ***************************************************************/
111
112// size of poly without exponents
113#ifdef TEST_MAC_DEBUG
114#define POLYSIZE (sizeof(poly) + sizeof(number) + 2*sizeof(Order_t))
115#else
116#define POLYSIZE (sizeof(poly) + sizeof(number) + sizeof(Order_t))
117#endif
118#define POLYSIZEW (POLYSIZE / sizeof(long))
119// number of Variables
120extern int pVariables;
121// size of a monom in bytes - always a multiple of sizeof(void*)
122extern int pMonomSize;
123// size of a monom in units of sizeof(void*) -- i.e. in words
124extern int pMonomSizeW;
125#ifdef COMP_FAST
126// Ceiling((pVariables+1) / sizeof(void*)) == length of exp-vector in words
127extern int pVariables1W;
128// Ceiling((pVariables) / sizeof(void*))
129extern int pVariablesW;
130extern int pCompIndex;
131extern int pVarOffset;
132extern int pVarLowIndex;
133extern int pVarHighIndex;
134#else
135#define pCompIndex 0
136#endif
137
138/***************************************************************
139 *
140 * Primitives for determening/setting  the way exponents are arranged
141 *
142 ***************************************************************/
143// And here is how we determine the way exponents are stored:
144// There are the following four possibilities:
145//
146//
147// BIGENDIAN -- lex order
148// e_1, e_2, ... , e_n,..,comp : pVarOffset = -1,
149//                               pCompIndex = pVariables + #(..)
150//                               pVarLowIndex = 0,
151//                               pVarHighIndex = pVariables-1
152// BIGENDIAN -- rev lex order
153// e_n, ... , e_2, e_1,..,comp : pVarOffset = pVariables,
154//                               pCompIndex = pVariables + #(..)
155//                               pVarLowIndex = 0,
156//                               pVarHighIndex = pVariables-1
157// LITTLEENDIAN -- rev lex order
158// comp,.., e_1, e_2, ... , e_n : pVarOffset = #(..),
159//                                pCompIndex = 0,
160//                                pVarLowIndex = 1 + #(..),
161//                                pVarHighIndex = #(..) + pVariables
162// LITTLEENDIAN -- lex order
163// comp,..,e_n, .... , e_2, e_1 : pVarOffset = pVariables + 1 + #(..)
164//                                pCompIndex = 0
165//                                pVarLowIndex = 1 + #(..)
166//                                pVarHighIndex = #(..) + pVariables
167//
168// Furthermore, the size of the exponent vector is always a multiple
169// of the word size -- "empty exponents" (exactly #(..) ones) are
170// filled in between comp and first/last exponent -- i.e. comp and
171// first/last exponent might not be next to each other
172
173#ifdef COMP_FAST
174
175#ifdef WORDS_BIGENDIAN
176
177#define _pHasReverseExp    (pVarOffset != -1)
178#define _pExpIndex(i)                           \
179  (pVarOffset == -1 ? (i) - 1 : pVarOffset - (i))
180#define _pRingExpIndex(r, i)                           \
181  ((r)->VarOffset == -1 ? (i) - 1 : (r)->VarOffset - (i))
182
183#else // ! WORDS_BIGENDIAN
184
185#define _pHasReverseExp    (pVarOffset > (SIZEOF_LONG / SIZEOF_EXPONENT) - 1)
186#define _pExpIndex(i)                                   \
187  (pVarOffset > (SIZEOF_LONG / SIZEOF_EXPONENT) - 1 ?   \
188   pVarOffset - (i) : pVarOffset + (i))
189#define _pRingExpIndex(r, i)                                \
190  ((r)->VarOffset > (SIZEOF_LONG / SIZEOF_EXPONENT) - 1 ?   \
191   (r)->VarOffset - (i) : (r)->VarOffset + (i))
192
193#endif // WORDS_BIGENDIAN
194
195inline void pGetVarIndicies_Lex(int nvars,
196                                int &VarOffset, int &VarCompIndex,
197                                int &VarLowIndex, int &VarHighIndex)
198{
199  long temp = (nvars+1)*sizeof(Exponent_t);
200  if ((temp % sizeof(long)) == 0)
201    temp = temp / sizeof(long);
202  else
203    temp = (temp / sizeof(long)) + 1; // now temp == nvars1W
204#ifdef WORDS_BIGENDIAN
205  VarCompIndex = temp * sizeof(long)/sizeof(Exponent_t) - 1;
206  VarOffset    = -1;
207  VarLowIndex  = 0;
208  VarHighIndex = nvars - 1;
209#else //  ! WORDS_BIGENDIAN
210  VarHighIndex = temp * sizeof(long)/sizeof(Exponent_t) - 1;
211  VarCompIndex = 0;
212  VarOffset    = VarHighIndex + 1;
213  VarLowIndex  = VarOffset - nvars;
214#endif // WORDS_BIGENDIAN
215}
216#define pSetVarIndicies_Lex(nvars) \
217  pGetVarIndicies_Lex(nvars,pVarOffset,pCompIndex,pVarLowIndex,pVarHighIndex)
218
219inline void pGetVarIndicies_RevLex(int nvars,
220                                   int &VarOffset, int &VarCompIndex,
221                                   int &VarLowIndex, int &VarHighIndex)
222{
223  long temp = (nvars+1)*sizeof(Exponent_t);
224  if ((temp % sizeof(long)) == 0)
225    temp = temp / sizeof(long);
226  else
227    temp = (temp / sizeof(long)) + 1;
228#ifdef WORDS_BIGENDIAN
229  VarCompIndex = temp * sizeof(long)/sizeof(Exponent_t) - 1;
230  VarOffset    = nvars;
231  VarLowIndex  = 0;
232  VarHighIndex = nvars-1;
233#else //  ! WORDS_BIGENDIAN
234  // comp, ..., e_1, e_2, ... , e_n
235  VarHighIndex = temp * sizeof(long)/sizeof(Exponent_t) - 1;
236  VarCompIndex = 0;
237  VarLowIndex  = VarHighIndex - nvars + 1;
238  VarOffset    = VarLowIndex - 1;
239#endif // WORDS_BIGENDIAN
240}
241#define pSetVarIndicies_RevLex(nvars) \
242 pGetVarIndicies_RevLex(nvars,pVarOffset,pCompIndex,pVarLowIndex,pVarHighIndex)
243
244// The default settings:
245inline void pGetVarIndicies(int nvars,
246                            int &VarOffset, int &VarCompIndex,
247                            int &VarLowIndex, int &VarHighIndex)
248{
249  pGetVarIndicies_Lex(nvars,VarOffset,VarCompIndex,VarLowIndex,VarHighIndex);
250}
251
252// gets var indicies w.r.t. the ring r
253extern void pGetVarIndicies(ring r,
254                            int &VarOffset, int &VarCompIndex,
255                            int &VarLowIndex, int &VarHighIndex);
256
257#define pSetVarIndicies(nvars) \
258  pGetVarIndicies(nvars, pVarOffset, pCompIndex, pVarLowIndex, pVarHighIndex)
259
260
261#else  // ! COMP_FAST
262#define _pExpIndex(i)       (i)
263#define _pRingExpIndex(r,i) (i)
264#endif // COMP_FAST
265
266/***************************************************************
267 *
268 * Primitives for accessing and seeting fields of a poly
269 *
270 ***************************************************************/
271#define _pNext(p)           ((p)->next)
272#define _pIter(p)           ((p) = (p)->next)
273
274#define _pGetCoeff(p)       ((p)->coef)
275#define _pSetCoeff(p,n)     {nDelete(&((p)->coef));(p)->coef=n;}
276#define _pSetCoeff0(p,n)    (p)->coef=n
277
278#ifndef TEST_MAC_DEBUG
279#define _pGetOrder(p)       ((p)->Order)
280#else
281#define _pGetOrder(p)       ((p)->MOrder)
282#endif
283
284#if defined(PDEBUG) && PDEBUG != 0
285extern Exponent_t pPDSetExp(poly p, int v, Exponent_t e, char* f, int l);
286extern Exponent_t pPDGetExp(poly p, int v, char* f, int l);
287extern Exponent_t pPDIncrExp(poly p, int v, char* f, int l);
288extern Exponent_t pPDDecrExp(poly p, int v, char* f, int l);
289extern Exponent_t pPDAddExp(poly p, int v, Exponent_t e, char* f, int l);
290extern Exponent_t pPDMultExp(poly p, int v, Exponent_t e, char* f, int l);
291extern Exponent_t pPDSubExp(poly p, int v, Exponent_t e, char* f, int l);
292
293extern Exponent_t pPDRingSetExp(ring r,poly p,int v,Exponent_t e,char* f,int l);
294extern Exponent_t pPDRingGetExp(ring r,poly p, int v, char* f, int l);
295
296#define _pSetExp(p,v,e)     pPDSetExp(p,v,e,__FILE__,__LINE__)
297#define _pGetExp(p,v)       pPDGetExp(p,v,__FILE__,__LINE__)
298#define _pIncrExp(p,v)      pPDIncrExp(p,v,__FILE__,__LINE__)
299#define _pDecrExp(p,v)      pPDDecrExp(p,v,__FILE__,__LINE__)
300#define _pAddExp(p,i,v)     pPDAddExp(p,i,v,__FILE__,__LINE__)
301#define _pSubExp(p,i,v)     pPDSubExp(p,i,v,__FILE__,__LINE__)
302#define _pMultExp(p,i,v)    pPDMultExp(p,i,v,__FILE__,__LINE__)
303
304#define _pRingSetExp(r,p,v,e)     pPDRingSetExp(r,p,v,e,__FILE__,__LINE__)
305#define _pRingGetExp(r,p,v)       pPDRingGetExp(r,p,v,__FILE__,__LINE__)
306
307#else  // ! (defined(PDEBUG) && PDEBUG != 0)
308
309#define _pSetExp(p,v,e)     (p)->exp[_pExpIndex(v)]=(e)
310#define _pGetExp(p,v)       (p)->exp[_pExpIndex(v)]
311#define _pIncrExp(p,v)      ((p)->exp[_pExpIndex(v)])++
312#define _pDecrExp(p,v)      ((p)->exp[_pExpIndex(v)])--
313#define _pAddExp(p,i,v)     ((p)->exp[_pExpIndex(i)]) += (v)
314#define _pSubExp(p,i,v)     ((p)->exp[_pExpIndex(i)]) -= (v)
315#define _pMultExp(p,i,v)    ((p)->exp[_pExpIndex(i)]) *= (v)
316
317#define _pRingSetExp(r,p,v,e)     (p)->exp[_pRingExpIndex(r,v)]=(e)
318#define _pRingGetExp(r,p,v)       (p)->exp[_pRingExpIndex(r,v)]
319
320#endif // defined(PDEBUG) && PDEBUG != 0
321
322inline Exponent_t _pGetExpSum(poly p1, poly p2, int i)
323{
324  int index = _pExpIndex(i);
325  return p1->exp[index] + p2->exp[index];
326}
327inline Exponent_t _pGetExpDiff(poly p1, poly p2, int i)
328{
329  int index = _pExpIndex(i);
330  return p1->exp[index] - p2->exp[index];
331}
332
333#define _pSetComp(p,k)      (p)->exp[pCompIndex] = (k)
334#define _pGetComp(p)        (p)->exp[pCompIndex]
335#define _pIncrComp(p)       (p)->exp[pCompIndex]++
336#define _pDecrComp(p)       (p)->exp[pCompIndex]--
337#define _pAddComp(p,v)      (p)->exp[pCompIndex] += (v)
338#define _pSubComp(p,v)      (p)->exp[pCompIndex] -= (v)
339
340#ifdef COMP_FAST
341#define _pRingSetComp(r,p,k)      (p)->exp[(r)->CompIndex] = (k)
342#define _pRingGetComp(r,p)        (p)->exp[(r)->CompIndex]
343#else
344#define _pRingSetComp(r,p,k)      _pSetComp(p,k)
345#define _pRingGetComp(r,p)        _pGetComp(p)
346#endif
347
348
349/***************************************************************
350 *
351 * Storage Managament Routines
352 *
353 ***************************************************************/
354#ifdef PDEBUG
355
356poly    pDBNew(char *f, int l);
357poly    pDBInit(char * f,int l);
358
359void    pDBDelete(poly * a, char * f, int l);
360void    pDBDelete1(poly * a, char * f, int l);
361void    pDBFree1(poly a, char * f, int l);
362
363poly    pDBCopy(poly a, char *f, int l);
364poly    pDBCopy1(poly a, char *f, int l);
365poly    pDBHead(poly a, char *f, int l);
366poly    pDBHead0(poly a, char *f, int l);
367poly    pDBFetchCopy(ring r, poly a, char *f, int l);
368
369void    pDBDelete(poly * a, char * f, int l);
370void    pDBDelete1(poly * a, char * f, int l);
371void    pDBFree1(poly a, char * f, int l);
372
373#define _pNew()         pDBNew(__FILE__,__LINE__)
374#define _pInit()        pDBInit(__FILE__,__LINE__)
375
376#define _pDelete(a)     pDBDelete((a),__FILE__,__LINE__)
377#define _pDelete1(a)    pDBDelete1((a),__FILE__,__LINE__)
378#define _pFree1(a)                              \
379do                                              \
380{                                               \
381  pDBFree1((a),__FILE__,__LINE__);              \
382  (a)=NULL;                                     \
383}                                               \
384while(0)
385
386#define _pCopy(A)       pDBCopy(A,__FILE__,__LINE__)
387#define _pCopy1(A)      pDBCopy1(A, __FILE__,__LINE__)
388#define _pHead(A)       pDBHead(A,__FILE__,__LINE__)
389#define _pHead0(A)      pDBHead0(A, __FILE__,__LINE__)
390#ifdef COMP_FAST
391#define _pFetchCopy(r,A)    pDBFetchCopy(r, A,__FILE__,__LINE__)
392#else
393#define _pFetchCopy(r,p)    pOrdPolyInsertSetm(pCopy(p))
394#endif
395
396#else // ! PDEBUG
397
398#ifdef MDEBUG
399#define _pNew()         (poly) mmDBAllocSpecialized(__FILE__,__LINE__)
400#else
401#define _pNew()         (poly) mmAllocSpecialized()
402#endif
403// #define _pNew() _pInit()
404
405#include <string.h>
406
407inline poly    _pInit(void)
408{
409#ifdef MDEBUG
410  poly p=(poly)mmDBAllocSpecialized(__FILE__,__LINE__);
411#else
412  poly p=(poly)mmAllocSpecialized();
413#endif
414  memset(p,0, pMonomSize);
415  return p;
416}
417
418extern void    _pDelete(poly * a);
419extern void    _pDelete1(poly * a);
420#ifdef MDEBUG
421#define _pFree1(a)       mmDBFreeSpecialized((ADDRESS)a,__FILE__,__LINE__)
422#else
423#define _pFree1(a)       mmFreeSpecialized((ADDRESS)a)
424#endif
425
426extern poly    _pCopy(poly a);
427extern poly    _pCopy1(poly a);
428extern poly    _pHead(poly a);
429extern poly    _pHead0(poly a);
430#ifdef COMP_FAST
431extern poly    _pFetchCopy(ring r,poly a);
432#else
433#define _pFetchCopy(r,p)  pOrdPolyInsertSetm(pCopy(p))
434#endif
435
436#endif // PDEBUG
437
438#define _pCopy2(p1, p2)     memcpyW(p1, p2, pMonomSizeW)
439
440// Here is a handy Macro which disables inlining when run with
441// profiling and enables it otherwise
442
443#ifdef DO_PROFILE
444
445#ifndef POLYS_IMPL_CC
446#define DECLARE(type, arglist) type arglist; \
447   static type dummy_##arglist
448#else
449#define DECLARE(type, arglist) type arglist
450#endif // POLYS_IMPL_CC
451
452#else //! DO_PROFILE
453
454#define DECLARE(type, arglist ) inline type arglist
455
456#endif // DO_PROFILE
457
458
459/***************************************************************
460 *
461 * Routines which work on vectors instead of single exponents
462 *
463 ***************************************************************/
464
465#ifdef COMP_FAST
466
467// nice declaration isn't it ??
468#if defined(PDEBUG) && PDEBUG == 1
469#define pMonAddFast(p1, p2)  pDBMonAddFast(p1, p2, __FILE__, __LINE__)
470extern  void pDBMonAddFast(poly p1, poly p2, char* f, int l);
471inline void _pMonAddFast(poly p1, poly p2)
472#else
473  DECLARE(void, pMonAddFast(poly p1, poly p2))
474#endif // defined(PDEBUG) && PDEBUG == 1
475{
476  // OK -- this might be the only place where we are kind of quick and
477  // dirty: the following only works correctly if all exponents are
478  // positive and the sum of two exponents does not exceed
479  // EXPONENT_MAX
480#ifndef COMP_NO_EXP_VECTOR_OPS
481  Exponent_t c2 = _pGetComp(p2);
482  int i = pVariables1W;
483  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
484  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
485  // set comp of p2 temporarily to 0, so that nothing is added to comp of p1
486  _pSetComp(p2, 0);
487#else
488  int i = pVariables;
489  Exponent_pt s1 = &(p1->exp[pVarLowIndex]);
490  Exponent_pt s2 = &(p2->exp[pVarLowIndex]);
491#endif
492
493  for (;;)
494  {
495    *s1 += *s2;
496    i--;
497    if (i==0) break;
498    s1++;
499    s2++;
500  }
501#ifndef COMP_NO_EXP_VECTOR_OPS
502  // reset comp of p2
503  _pSetComp(p2, c2);
504#endif
505#ifdef TEST_MAC_ORDER
506  if (bNoAdd) bSetm(p1);else
507#endif
508  _pGetOrder(p1) += _pGetOrder(p2);
509}
510
511// Makes p1 a copy of p2 and adds on exponets of p3
512#if defined(PDEBUG) && PDEBUG == 1
513#define _pCopyAddFast(p1, p2, p3)  pDBCopyAddFast(p1, p2, p3, __FILE__, __LINE__)
514extern  void pDBCopyAddFast(poly p1, poly p2, poly p3, char* f, int l);
515inline void __pCopyAddFast(poly p1, poly p2, poly p3)
516#else
517  DECLARE(void, _pCopyAddFast(poly p1, poly p2, poly p3))
518#endif // defined(PDEBUG) && PDEBUG == 1
519{
520  p1->next = p2->next;
521  p1->coef = p2->coef;
522//  memset(p1, 0, pMonomSize);
523
524  //memset(p1,0,pMonomSize);
525#ifndef COMP_NO_EXP_VECTOR_OPS
526  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
527  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
528  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
529  const unsigned long* const ub = s3 + pVariables1W;
530#else
531  Exponent_t* s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
532  const Exponent_t* s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
533  const Exponent_t* s3 = (Exponent_t*) &(p3->exp[pVarLowIndex]);
534  const Exponent_t* const ub = s3 + pVariables;
535// need to zero the "fill in" slots (i.e., empty exponents)
536#ifdef  WORDS_BIGENDIAN
537  *((unsigned long *) ((unsigned long*) p1) + pMonomSizeW -1) = 0;
538#else
539  *((unsigned long *) p1->exp) = 0;
540#endif
541#endif
542
543  for (;;)
544  {
545    *s1 = *s2 + *s3;
546    s3++;
547    if (s3 == ub) break;
548    s1++;
549    s2++;
550  }
551  // we first are supposed to do a copy from p2 to p1 -- therefore,
552  // component of p1 is set to comp of p2
553  _pSetComp(p1, _pGetComp(p2));
554  #ifdef TEST_MAC_ORDER
555  if (bNoAdd) bSetm(p1);else
556  #endif
557  _pGetOrder(p1) = _pGetOrder(p2) + _pGetOrder(p3);
558}
559
560// Similar to pCopyAddFast, except that we assume that the component
561// of p2 and p3 is zero component
562#if defined(PDEBUG) && PDEBUG == 1
563#define _pCopyAddFast1(p1, p2, p3)  pDBCopyAddFast(p1, p2, p3, __FILE__, __LINE__)
564extern  void pDBCopyAddFast(poly p1, poly p2, poly p3, char* f, int l);
565inline void __pCopyAddFast1(poly p1, poly p2, poly p3)
566#else
567  DECLARE(void, _pCopyAddFast1(poly p1, poly p2, poly p3))
568#endif // defined(PDEBUG) && PDEBUG == 1
569{
570#ifndef COMP_NO_EXP_VECTOR_OPS
571  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
572  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
573  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
574  const unsigned long* const ub = s3 + pVariables1W;
575#else
576  Exponent_t* s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
577  const Exponent_t* s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
578  const Exponent_t* s3 = (Exponent_t*) &(p3->exp[pVarLowIndex]);
579  const Exponent_t* const ub = s3 + pVariables;
580#ifdef  WORDS_BIGENDIAN
581  *((unsigned long *) ((unsigned long*) p1) + pMonomSizeW -1) = 0;
582#else
583  *((unsigned long *) p1->exp) = 0;
584#endif
585#endif
586
587  for (;;)
588  {
589    *s1 = *s2 + *s3;
590    s3++;
591    if (s3 == ub) break;
592    s1++;
593    s2++;
594  }
595  #ifdef TEST_MAC_ORDER
596  if (bNoAdd) bSetm(p1);else
597  #endif
598  _pGetOrder(p1) = _pGetOrder(p2) + _pGetOrder(p3);
599}
600
601
602#ifndef COMP_NO_EXP_VECTOR_OPS
603
604#if SIZEOF_LONG == 4
605
606#if SIZEOF_EXPONENT == 1
607#define P_DIV_MASK 0x80808080
608#else // SIZEOF_EXPONENT == 2
609#define P_DIV_MASK 0x80008000
610#endif
611
612#else // SIZEOF_LONG == 8
613
614#if SIZEOF_EXPONENT == 1
615#define P_DIV_MASK 0x8080808080808080
616#elif  SIZEOF_EXPONENT == 2
617#define P_DIV_MASK 0x8000800080008000
618#else // SIZEOF_EXPONENT == 4
619#define P_DIV_MASK 0x8000000080000000
620#endif
621
622#endif
623
624DECLARE(BOOLEAN, __pDivisibleBy(poly a, poly b))
625{
626  const unsigned long* s1;
627  const unsigned long* s2;
628  const unsigned long* lb;
629
630#ifdef WORDS_BIGENDIAN
631  lb = (unsigned long*) &(a->exp[0]);
632  if (pVariables & ((SIZEOF_LONG / SIZEOF_EXPONENT) - 1))
633  {
634    // now pVariables == pVariables1W, i.e. there are exponents in the
635    // "first" word of exponentvector
636    s1 = ((unsigned long*) a) + pMonomSizeW -1;
637    s2 = ((unsigned long*) b) + pMonomSizeW -1;
638  }
639  else
640  {
641    // first exponent word has only component as significant field --
642    // Hence, do not bother
643    s1 = ((unsigned long*) a) + pMonomSizeW -2;
644    s2 = ((unsigned long*) b) + pMonomSizeW -2;
645  }
646#else // !WORDS_BIGENDIAN
647  lb = ((unsigned long*) a) + pMonomSizeW;
648  if (pVariables & ((SIZEOF_LONG / SIZEOF_EXPONENT) - 1))
649  {
650    s1 = (unsigned long*) &(a->exp[0]);
651    s2 = (unsigned long*) &(b->exp[0]);
652  }
653  else
654  {
655    s1 = (unsigned long*) &(a->exp[0]) + 1;
656    s2 = (unsigned long*) &(b->exp[0]) + 1;
657  }
658#endif
659  for (;;)
660  {
661// O.K. -- and now comes a bit of magic. The following _really_
662// works. Think about it! If you can prove it, please tell me, for I
663// did not bother to prove it formally (Hint: We can assume that our
664// exponents are always positive).
665    if ((*s2 - *s1) & P_DIV_MASK) return FALSE;
666#ifdef WORDS_BIGENDIAN
667    if (s1 == lb) return TRUE;
668    s1--;
669    s2--;
670#else
671    s1++;
672    if (s1 == lb) return TRUE;
673    s2++;
674#endif
675  }
676}
677
678#else //  ! COMP_NO_EXP_VECTOR_OPS
679
680DECLARE(BOOLEAN, __pDivisibleBy(poly a, poly b))
681{
682#ifdef WORDS_BIGENDIAN
683  const Exponent_t* s1 = &(a->exp[pVarHighIndex]);
684  const Exponent_t* s2 = &(b->exp[pVarHighIndex]);
685  const Exponent_t* lb = s1 - pVariables;
686
687  for (;;)
688  {
689    if (*s1 > *s2) return FALSE;
690    s1--;
691    if (s1 == lb) return TRUE;
692    s2--;
693  }
694
695#else // !WORDS_BIGENDIAN
696  const Exponent_t* s1 = &(a->exp[pVarLowIndex]);
697  const Exponent_t* s2 = &(b->exp[pVarLowIndex]);
698  const Exponent_t* lb = s1 + pVariables;
699
700  for (;;)
701  {
702   if (*s1 > *s2) return FALSE;
703   s1++;
704   if (s1 == lb) return TRUE;
705   s2++;
706  }
707#endif  // WORDS_BIGENDIAN
708}
709
710#endif //  COMP_NO_EXP_VECTOR_OPS
711
712#if defined(PDEBUG) && PDEBUG == 1
713#define _pDivisibleBy(a,b)   pDBDivisibleBy(a, b, __FILE__, __LINE__)
714extern  BOOLEAN pDBDivisibleBy(poly p1, poly p2, char* f, int l);
715inline BOOLEAN _pDivisibleBy_orig(poly a, poly b)
716#else
717inline BOOLEAN _pDivisibleBy(poly a, poly b)
718#endif // defined(PDEBUG) && PDEBUG == 1
719{
720  if ((a!=NULL)&&((_pGetComp(a)==0) || (_pGetComp(a) == _pGetComp(b))))
721  {
722    return __pDivisibleBy(a,b);
723  }
724  return FALSE;
725}
726
727#if defined(PDEBUG) && PDEBUG == 1
728#define _pDivisibleBy1(a,b)   pDBDivisibleBy1(a, b, __FILE__, __LINE__)
729extern  BOOLEAN pDBDivisibleBy1(poly p1, poly p2, char* f, int l);
730inline BOOLEAN _pDivisibleBy1_orig(poly a, poly b)
731#else
732inline BOOLEAN _pDivisibleBy1(poly a, poly b)
733#endif // defined(PDEBUG) && PDEBUG == 1
734{
735  if (_pGetComp(a) == 0 || _pGetComp(a) == _pGetComp(b))
736    return __pDivisibleBy(a,b);
737  return FALSE;
738}
739
740#if defined(PDEBUG) && PDEBUG == 1
741#define _pDivisibleBy2(a,b)   pDBDivisibleBy2(a, b, __FILE__, __LINE__)
742extern  BOOLEAN pDBDivisibleBy2(poly p1, poly p2, char* f, int l);
743#else
744#define _pDivisibleBy2(a,b) __pDivisibleBy(a,b)
745#endif // defined(PDEBUG) && PDEBUG == 1
746
747
748DECLARE(BOOLEAN, _pEqual(poly p1, poly p2))
749{
750#ifndef COMP_NO_EXP_VECTOR_OPS
751  const long *s1 = (long*) &(p1->exp[0]);
752  const long *s2 = (long*) &(p2->exp[0]);
753  const long* const lb = s1 + pVariables1W;
754#else
755  const Exponent_t *s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
756  const Exponent_t *s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
757  const Exponent_t* const lb = s1 + pVariables;
758  if (_pGetComp(p1) != _pGetComp(p2)) return FALSE;
759#endif
760
761  for(;;)
762  {
763    if (*s1 != *s2) return FALSE;
764    s1++;
765    if (s1 == lb) return TRUE;
766    s2++;
767  }
768}
769
770inline void _pGetExpV(poly p, Exponent_t *ev)
771{
772  if (_pHasReverseExp)
773  {
774    for (int i = pVarLowIndex, j = pVariables; j; i++, j--)
775      ev[j] = p->exp[i];
776  }
777  else
778    memcpy(&ev[1], &(p->exp[pVarLowIndex]), pVariables*sizeof(Exponent_t));
779  ev[0] = _pGetComp(p);
780}
781
782extern pSetmProc pSetm;
783inline void _pSetExpV(poly p, Exponent_t *ev)
784{
785  if (_pHasReverseExp)
786  {
787    for (int i = pVarLowIndex, j = pVariables; j; i++, j--)
788      p->exp[i] = ev[j];
789  }
790  else
791    memcpy(&(p->exp[pVarLowIndex]), &ev[1], pVariables*sizeof(Exponent_t));
792  _pSetComp(p, ev[0]);
793  pSetm(p);
794}
795
796#else // ! COMP_FAST
797
798DECLARE(BOOLEAN, _pDivisibleBy(poly a, poly b))
799{
800  if ((a!=NULL)&&((a->exp[0]==0) || (a->exp[0] == b->exp[0])))
801  {
802    int i=pVariables;
803    short *e1=&(a->exp[1]);
804    short *e2=&(b->exp[1]);
805    if ((*e1) > (*e2)) return FALSE;
806    do
807    {
808      i--;
809      if (i == 0) return TRUE;
810      e1++;
811      e2++;
812    } while ((*e1) <= (*e2));
813  }
814  return FALSE;
815}
816
817#define _pDivisibleBy1(a,b) _pDivisibleBy(a,b)
818#define _pDivisibleBy2(a,b) _pDivisibleBy(a,b)
819
820#ifdef TEST_MAC_ORDER
821DECLARE(void, pMonAddFast(poly a, poly m))
822{
823  for(int ii =pVariables; ii; ii--) (a)->exp[ii] += (m)->exp[ii];\
824  if (bNoAdd) bSetm(a); else
825  _pGetOrder(a) += _pGetOrder(m);
826}
827#else
828DECLARE(void, pMonAddFast(poly a, poly m))
829{
830  for(int ii =pVariables; ii; ii--) (a)->exp[ii] += (m)->exp[ii];\
831  _pGetOrder(a) += _pGetOrder(m);
832}
833#endif
834
835DECLARE(BOOLEAN, _pEqual(poly p1, poly p2))
836{
837  int i;
838  short *e1=p1->exp;
839  short *e2=p2->exp;
840
841  if (p1->Order != p2->Order) return FALSE;
842  for (i=pVariables; i>=0; i--,e1++,e2++)
843    if (*e1 != *e2) return FALSE;
844  return TRUE;
845}
846
847#define _pGetExpV(p,e)    memcpy((e),(p)->exp,(pVariables+1)*sizeof(short));
848//void    pSetExpV(poly p, short * exp);
849#define _pSetExpV(p,e) {memcpy((p)->exp,(e),(pVariables+1)*sizeof(short));pSetm(p);}
850
851#endif // COMP_FAST
852
853/***************************************************************
854 *
855 * Routines which implement low-level manipulations/operations on exponents
856 *
857 ***************************************************************/
858
859DECLARE(int, __pExpQuerSum2(poly p, int from, int to))
860{
861  int j = p->exp[from];
862  int i = from + 1;
863
864  for(;;)
865  {
866    if (i > to) return j;
867    j += p->exp[i];
868    i++;
869  }
870}
871
872#ifdef COMP_FAST
873#define _pExpQuerSum(p)  __pExpQuerSum2(p, pVarLowIndex, pVarHighIndex)
874
875#define _pExpQuerSum1(p,to) \
876 (_pHasReverseExp ? \
877    __pExpQuerSum2(p, _pExpIndex(to), _pExpIndex(1)) : \
878    __pExpQuerSum2(p, _pExpIndex(1), _pExpIndex(to)))
879
880#define _pExpQuerSum2(p,from,to) \
881 (_pHasReverseExp ? \
882    __pExpQuerSum2(p, _pExpIndex(to), _pExpIndex(from)) : \
883    __pExpQuerSum2(p, _pExpIndex(from), _pExpIndex(to)))
884#else
885
886#define _pExpQuerSum(p)             __pExpQuerSum2(p, 1, pVariables)
887#define _pExpQuerSum1(p, to)        __pExpQuerSum2(p, 1, to)
888#define _pExpQuerSum2(p, from, to)  __pExpQuerSum2(p, from, to)
889
890#endif
891
892/***************************************************************
893 *
894 * Routines which implement macaulay ordering routines
895 *
896 ***************************************************************/
897#ifdef TEST_MAC_ORDER
898
899DECLARE(void, _bSetm0(poly p))
900{
901
902  int i=1;
903  int ord = -INT_MAX;
904  Exponent_t *ep;
905
906  if(_pHasReverseExp)
907  {
908    ep=&(p->exp[pVarHighIndex]);
909    int *ip=bBinomials+(*ep); /*_pGetExp(p,1);*/
910    loop
911    {
912      ord += (*ip);
913      if (i==pVariables) break;
914      i++;
915      //ip+=bHighdeg_1+_pGetExp(p,i);
916      ep--;
917      ip+=bHighdeg_1+(*ep);
918    }
919#if 0
920    int *ip=bBinomials+_pGetExp(p,1);
921    loop
922    {
923      ord += (*ip);
924      if (i==pVariables) break;
925      i++;
926      ip+=bHighdeg_1+_pGetExp(p,i);
927    }
928#endif
929  }
930  else
931  {
932    ep=&(p->exp[pVarLowIndex]);
933    int *ip=bBinomials+(*ep); /*_pGetExp(p,1);*/
934    loop
935    {
936      ord += (*ip);
937      if (i==pVariables) break;
938      i++;
939      //ip+=bHighdeg_1+_pGetExp(p,i);
940      ep++;
941      ip+=bHighdeg_1+(*ep);
942    }
943#if 0
944    int *ip=bBinomials+_pGetExp(p,1);
945    loop
946    {
947      ord += (*ip);
948      if (i==pVariables) break;
949      i++;
950      ip+=bHighdeg_1+_pGetExp(p,i);
951    }
952#endif
953  }
954  p->Order=ord;
955}
956
957DECLARE(void, _bSetm(poly p))
958{
959  int ord = _pExpQuerSum(p);
960
961  if (ord<bHighdeg)
962    _bSetm0(p);
963  else
964    p->Order=ord;
965}
966
967// ordering dp,c or c,dp, general case
968#if defined(PDEBUG) && PDEBUG == 1
969#define pbMonAddFast(p1, p2)  pDBMonAddFast(p1, p2, __FILE__, __LINE__)
970extern  void pbDBMonAddFast(poly p1, poly p2, char* f, int l);
971inline void _pbMonAddFast(poly p1, poly p2)
972#else
973#define pbMonAddFast(p1, p2)  _pbMonAddFast(p1, p2)
974DECLARE(void, _pbMonAddFast(poly p1, poly p2))
975#endif // defined(PDEBUG) && PDEBUG == 1
976{
977  // OK -- this might be the only place where we are kind of quick and
978  // dirty: the following only works correctly if all exponents are
979  // positive and the sum of two exponents does not exceed
980  // EXPONENT_MAX
981#ifndef COMP_NO_EXP_VECTOR_OPS
982  Exponent_t c2 = _pGetComp(p2);
983  int i = pVariables1W;
984  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
985  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
986  // set comp of p2 temporarily to 0, so that nothing is added to comp of p1
987  _pSetComp(p2, 0);
988#else
989  int i = pVariables;
990  Exponent_pt s1 = &(p1->exp[pVarLowIndex]);
991  Exponent_pt s2 = &(p2->exp[pVarLowIndex]);
992#endif
993
994  for (;;)
995  {
996    *s1 += *s2;
997    i--;
998    if (i==0) break;
999    s1++;
1000    s2++;
1001  }
1002#ifndef COMP_NO_EXP_VECTOR_OPS
1003  // reset comp of p2
1004  _pSetComp(p2, c2);
1005#endif
1006  if ((_pGetOrder(p1)|_pGetOrder(p2))>0) // i.e. overflow of mac order
1007    _pGetOrder(p1) += _pGetOrder(p2);
1008  else
1009    _bSetm(p1);
1010}
1011
1012// ordering dp,c or c,dp, below degree limit
1013#if defined(PDEBUG) && PDEBUG == 1
1014#define pbMonAddFast0(p1, p2)  pbDBMonAddFast0(p1, p2, __FILE__, __LINE__)
1015extern  void pbDBMonAddFast0(poly p1, poly p2, char* f, int l);
1016inline void _pbMonAddFast0(poly p1, poly p2)
1017#else
1018  DECLARE(void, pbMonAddFast0(poly p1, poly p2))
1019#endif // defined(PDEBUG) && PDEBUG == 1
1020{
1021  // OK -- this might be the only place where we are kind of quick and
1022  // dirty: the following only works correctly if all exponents are
1023  // positive and the sum of two exponents does not exceed
1024  // EXPONENT_MAX
1025#ifndef COMP_NO_EXP_VECTOR_OPS
1026  Exponent_t c2 = _pGetComp(p2);
1027  int i = pVariables1W;
1028  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
1029  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
1030  // set comp of p2 temporarily to 0, so that nothing is added to comp of p1
1031  _pSetComp(p2, 0);
1032#else
1033  int i = pVariables;
1034  Exponent_pt s1 = &(p1->exp[pVarLowIndex]);
1035  Exponent_pt s2 = &(p2->exp[pVarLowIndex]);
1036#endif
1037
1038  for (;;)
1039  {
1040    *s1 += *s2;
1041    i--;
1042    if (i==0) break;
1043    s1++;
1044    s2++;
1045  }
1046#ifndef COMP_NO_EXP_VECTOR_OPS
1047  // reset comp of p2
1048  _pSetComp(p2, c2);
1049#endif
1050  _bSetm0(p1);
1051}
1052
1053// ordering dp,c or c,dp, below degree limit
1054// Makes p1 a copy of p2 and adds on exponets of p3
1055#if defined(PDEBUG) && PDEBUG == 1
1056#define _pbCopyAddFast0(p1, p2, p3)  pDBCopyAddFast(p1, p2, p3, __FILE__, __LINE__)
1057inline void __pbCopyAddFast0(poly p1, poly p2, poly p3)
1058#else
1059  DECLARE(void, _pbCopyAddFast0(poly p1, poly p2, poly p3))
1060#endif // defined(PDEBUG) && PDEBUG == 1
1061{
1062  p1->next = p2->next;
1063  p1->coef = p2->coef;
1064
1065#ifndef COMP_NO_EXP_VECTOR_OPS
1066  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
1067  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
1068  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
1069  const unsigned long* const ub = s3 + pVariables1W;
1070#else
1071  Exponent_t* s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
1072  const Exponent_t* s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
1073  const Exponent_t* s3 = (Exponent_t*) &(p3->exp[pVarLowIndex]);
1074  const Exponent_t* const ub = s3 + pVariables;
1075// need to zero the "fill in" slots (i.e., empty exponents)
1076#ifdef  WORDS_BIGENDIAN
1077  *((unsigned long*) p1 + pMonomSizeW -1) = 0;
1078#else
1079  *((unsigned long *) p1->exp) = 0;
1080#endif
1081#endif
1082
1083  for (;;)
1084  {
1085    *s1 = *s2 + *s3;
1086    s3++;
1087    if (s3 == ub) break;
1088    s1++;
1089    s2++;
1090  }
1091  // we first are supposed to do a copy from p2 to p1 -- therefore,
1092  // component of p1 is set to comp of p2
1093  _pSetComp(p1, _pGetComp(p2));
1094  _bSetm0(p1);
1095}
1096
1097// Similar to pCopyAddFast, except that we assume that the component
1098// of p2 and p3 is zero component
1099#if defined(PDEBUG) && PDEBUG == 1
1100#define _pbCopyAddFast10(p1, p2, p3)  pbDBCopyAddFast0(p1, p2, p3, __FILE__, __LINE__)
1101extern  void pbDBCopyAddFast0(poly p1, poly p2, poly p3, char* f, int l);
1102inline void __pbCopyAddFast10(poly p1, poly p2, poly p3)
1103#else
1104  DECLARE(void, _pbCopyAddFast10(poly p1, poly p2, poly p3))
1105#endif // defined(PDEBUG) && PDEBUG == 1
1106{
1107  p1->next = p2->next;
1108  p1->coef = p2->coef;
1109
1110#ifndef COMP_NO_EXP_VECTOR_OPS
1111  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
1112  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
1113  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
1114  const unsigned long* const ub = s3 + pVariables1W;
1115#else
1116  Exponent_t* s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
1117  const Exponent_t* s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
1118  const Exponent_t* s3 = (Exponent_t*) &(p3->exp[pVarLowIndex]);
1119  const Exponent_t* const ub = s3 + pVariables;
1120#ifdef  WORDS_BIGENDIAN
1121  *((unsigned long*) p1 + pMonomSizeW -1) = 0;
1122#else
1123  *((unsigned long *) p1->exp) = 0;
1124#endif
1125#endif
1126
1127  for (;;)
1128  {
1129    *s1 = *s2 + *s3;
1130    s3++;
1131    if (s3 == ub) break;
1132    s1++;
1133    s2++;
1134  }
1135  _bSetm0(p1);
1136}
1137
1138#endif
1139
1140#endif // POLYS_IMPL_H
1141
Note: See TracBrowser for help on using the repository browser.