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

spielwiese
Last change on this file since c9567bf was c9567bf, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: bug fixes to TEST_MAC_ORDER (binom.cc binom.h polys-impl.h spSpolyLoop.cc) git-svn-id: file:///usr/local/Singular/svn/trunk@1062 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 31.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.15 1998-01-24 17:22:06 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#define _pGetOrder(p)       ((p)->Order)
279
280#if defined(PDEBUG) && PDEBUG != 0
281extern Exponent_t pPDSetExp(poly p, int v, Exponent_t e, char* f, int l);
282extern Exponent_t pPDGetExp(poly p, int v, char* f, int l);
283extern Exponent_t pPDIncrExp(poly p, int v, char* f, int l);
284extern Exponent_t pPDDecrExp(poly p, int v, char* f, int l);
285extern Exponent_t pPDAddExp(poly p, int v, Exponent_t e, char* f, int l);
286extern Exponent_t pPDMultExp(poly p, int v, Exponent_t e, char* f, int l);
287extern Exponent_t pPDSubExp(poly p, int v, Exponent_t e, char* f, int l);
288
289extern Exponent_t pPDRingSetExp(ring r,poly p,int v,Exponent_t e,char* f,int l);
290extern Exponent_t pPDRingGetExp(ring r,poly p, int v, char* f, int l);
291
292#define _pSetExp(p,v,e)     pPDSetExp(p,v,e,__FILE__,__LINE__)
293#define _pGetExp(p,v)       pPDGetExp(p,v,__FILE__,__LINE__)
294#define _pIncrExp(p,v)      pPDIncrExp(p,v,__FILE__,__LINE__)
295#define _pDecrExp(p,v)      pPDDecrExp(p,v,__FILE__,__LINE__)
296#define _pAddExp(p,i,v)     pPDAddExp(p,i,v,__FILE__,__LINE__)
297#define _pSubExp(p,i,v)     pPDSubExp(p,i,v,__FILE__,__LINE__)
298#define _pMultExp(p,i,v)    pPDMultExp(p,i,v,__FILE__,__LINE__)
299
300#define _pRingSetExp(r,p,v,e)     pPDRingSetExp(r,p,v,e,__FILE__,__LINE__)
301#define _pRingGetExp(r,p,v)       pPDRingGetExp(r,p,v,__FILE__,__LINE__)
302
303#else  // ! (defined(PDEBUG) && PDEBUG != 0)
304
305#define _pSetExp(p,v,e)     (p)->exp[_pExpIndex(v)]=(e)
306#define _pGetExp(p,v)       (p)->exp[_pExpIndex(v)]
307#define _pIncrExp(p,v)      ((p)->exp[_pExpIndex(v)])++
308#define _pDecrExp(p,v)      ((p)->exp[_pExpIndex(v)])--
309#define _pAddExp(p,i,v)     ((p)->exp[_pExpIndex(i)]) += (v)
310#define _pSubExp(p,i,v)     ((p)->exp[_pExpIndex(i)]) -= (v)
311#define _pMultExp(p,i,v)    ((p)->exp[_pExpIndex(i)]) *= (v)
312
313#define _pRingSetExp(r,p,v,e)     (p)->exp[_pRingExpIndex(r,v)]=(e)
314#define _pRingGetExp(r,p,v)       (p)->exp[_pRingExpIndex(r,v)]
315
316#endif // defined(PDEBUG) && PDEBUG != 0
317
318inline Exponent_t _pGetExpSum(poly p1, poly p2, int i)
319{
320  int index = _pExpIndex(i);
321  return p1->exp[index] + p2->exp[index];
322}
323inline Exponent_t _pGetExpDiff(poly p1, poly p2, int i)
324{
325  int index = _pExpIndex(i);
326  return p1->exp[index] - p2->exp[index];
327}
328
329#define _pSetComp(p,k)      (p)->exp[pCompIndex] = (k)
330#define _pGetComp(p)        (p)->exp[pCompIndex]
331#define _pIncrComp(p)       (p)->exp[pCompIndex]++
332#define _pDecrComp(p)       (p)->exp[pCompIndex]--
333#define _pAddComp(p,v)      (p)->exp[pCompIndex] += (v)
334#define _pSubComp(p,v)      (p)->exp[pCompIndex] -= (v)
335
336#ifdef COMP_FAST
337#define _pRingSetComp(r,p,k)      (p)->exp[(r)->CompIndex] = (k)
338#define _pRingGetComp(r,p)        (p)->exp[(r)->CompIndex]
339#else
340#define _pRingSetComp(r,p,k)      _pSetComp(p,k)
341#define _pRingGetComp(r,p)        _pGetComp(p)
342#endif
343
344
345/***************************************************************
346 *
347 * Storage Managament Routines
348 *
349 ***************************************************************/
350#ifdef MDEBUG
351
352poly    pDBNew(char *f, int l);
353poly    pDBInit(char * f,int l);
354
355void    pDBDelete(poly * a, char * f, int l);
356void    pDBDelete1(poly * a, char * f, int l);
357void    pDBFree1(poly a, char * f, int l);
358
359poly    pDBCopy(poly a, char *f, int l);
360poly    pDBCopy1(poly a, char *f, int l);
361poly    pDBHead(poly a, char *f, int l);
362poly    pDBHead0(poly a, char *f, int l);
363poly    pDBFetchCopy(ring r, poly a, char *f, int l);
364
365void    pDBDelete(poly * a, char * f, int l);
366void    pDBDelete1(poly * a, char * f, int l);
367void    pDBFree1(poly a, char * f, int l);
368
369#define _pNew()         pDBNew(__FILE__,__LINE__)
370#define _pInit()        pDBInit(__FILE__,__LINE__)
371
372#define _pDelete(a)     pDBDelete((a),__FILE__,__LINE__)
373#define _pDelete1(a)    pDBDelete1((a),__FILE__,__LINE__)
374#define _pFree1(a)                              \
375do                                              \
376{                                               \
377  pDBFree1((a),__FILE__,__LINE__);              \
378  (a)=NULL;                                     \
379}                                               \
380while(0)
381
382#define _pCopy(A)       pDBCopy(A,__FILE__,__LINE__)
383#define _pCopy1(A)      pDBCopy1(A, __FILE__,__LINE__)
384#define _pHead(A)       pDBHead(A,__FILE__,__LINE__)
385#define _pHead0(A)      pDBHead0(A, __FILE__,__LINE__)
386#ifdef COMP_FAST
387#define _pFetchCopy(r,A)    pDBFetchCopy(r, A,__FILE__,__LINE__)
388#else
389#define _pFetchCopy(r,p)    pOrdPolyInsertSetm(pCopy(p))
390#endif
391
392#else // ! MDEBUG
393
394#define _pNew()         (poly) mmAllocSpecialized()
395// #define _pNew() _pInit()
396
397#include <string.h>
398
399inline poly    _pInit(void)
400{
401  poly p=(poly)mmAllocSpecialized();
402  memset(p,0, pMonomSize);
403  return p;
404}
405
406extern void    _pDelete(poly * a);
407extern void    _pDelete1(poly * a);
408#define _pFree1(a)       mmFreeSpecialized((ADDRESS)a)
409
410extern poly    _pCopy(poly a);
411extern poly    _pCopy1(poly a);
412extern poly    _pHead(poly a);
413extern poly    _pHead0(poly a);
414#ifdef COMP_FAST
415extern poly    _pFetchCopy(ring r,poly a);
416#else
417#define _pFetchCopy(r,p)  pOrdPolyInsertSetm(pCopy(p))
418#endif
419
420#endif // MDEBUG
421
422#define _pCopy2(p1, p2)     memcpyW(p1, p2, pMonomSizeW)
423
424// Here is a handy Macro which disables inlining when run with
425// profiling and enables it otherwise
426
427#ifdef DO_PROFILE
428
429#ifndef POLYS_IMPL_CC
430#define DECLARE(type, arglist) type arglist; \
431   static type dummy_##arglist
432#else
433#define DECLARE(type, arglist) type arglist
434#endif // POLYS_IMPL_CC
435
436#else //! DO_PROFILE
437
438#define DECLARE(type, arglist ) inline type arglist
439
440#endif // DO_PROFILE
441
442
443/***************************************************************
444 *
445 * Routines which work on vectors instead of single exponents
446 *
447 ***************************************************************/
448
449#ifdef COMP_FAST
450
451// nice declaration isn't it ??
452#if defined(PDEBUG) && PDEBUG == 1
453#define pMonAddFast(p1, p2)  pDBMonAddFast(p1, p2, __FILE__, __LINE__)
454extern  void pDBMonAddFast(poly p1, poly p2, char* f, int l);
455inline void _pMonAddFast(poly p1, poly p2)
456#else
457  DECLARE(void, pMonAddFast(poly p1, poly p2))
458#endif // defined(PDEBUG) && PDEBUG == 1
459{
460  // OK -- this might be the only place where we are kind of quick and
461  // dirty: the following only works correctly if all exponents are
462  // positive and the sum of two exponents does not exceed
463  // EXPONENT_MAX
464#ifndef COMP_NO_EXP_VECTOR_OPS
465  Exponent_t c2 = _pGetComp(p2);
466  int i = pVariables1W;
467  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
468  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
469  // set comp of p2 temporarily to 0, so that nothing is added to comp of p1
470  _pSetComp(p2, 0);
471#else
472  int i = pVariables;
473  Exponent_pt s1 = &(p1->exp[pVarLowIndex]);
474  Exponent_pt s2 = &(p2->exp[pVarLowIndex]);
475#endif
476
477  for (;;)
478  {
479    *s1 += *s2;
480    i--;
481    if (i==0) break;
482    s1++;
483    s2++;
484  }
485#ifndef COMP_NO_EXP_VECTOR_OPS
486  // reset comp of p2
487  _pSetComp(p2, c2);
488#endif
489#ifdef TEST_MAC_ORDER
490  if (bNoAdd) bSetm(p1);else
491#endif
492  _pGetOrder(p1) += _pGetOrder(p2);
493}
494
495// Makes p1 a copy of p2 and adds on exponets of p3
496#if defined(PDEBUG) && PDEBUG == 1
497#define _pCopyAddFast(p1, p2, p3)  pDBCopyAddFast(p1, p2, p3, __FILE__, __LINE__)
498extern  void pDBCopyAddFast(poly p1, poly p2, poly p3, char* f, int l);
499inline void __pCopyAddFast(poly p1, poly p2, poly p3)
500#else
501  DECLARE(void, _pCopyAddFast(poly p1, poly p2, poly p3))
502#endif // defined(PDEBUG) && PDEBUG == 1
503{
504  p1->next = p2->next;
505  p1->coef = p2->coef;
506//  memset(p1, 0, pMonomSize);
507
508  //memset(p1,0,pMonomSize);
509#ifndef COMP_NO_EXP_VECTOR_OPS
510  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
511  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
512  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
513  const unsigned long* const ub = s3 + pVariables1W;
514#else
515  Exponent_t* s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
516  const Exponent_t* s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
517  const Exponent_t* s3 = (Exponent_t*) &(p3->exp[pVarLowIndex]);
518  const Exponent_t* const ub = s3 + pVariables;
519// need to zero the "fill in" slots (i.e., empty exponents) 
520#ifdef  WORDS_BIGENDIAN
521  *((unsigned long *) ((unsigned long*) p1) + pMonomSizeW -1) = 0;
522#else
523  *((unsigned long *) p1->exp) = 0;
524#endif
525#endif
526
527  for (;;)
528  {
529    *s1 = *s2 + *s3;
530    s3++;
531    if (s3 == ub) break;
532    s1++;
533    s2++;
534  }
535  // we first are supposed to do a copy from p2 to p1 -- therefore,
536  // component of p1 is set to comp of p2
537  _pSetComp(p1, _pGetComp(p2));
538  #ifdef TEST_MAC_ORDER
539  if (bNoAdd) bSetm(p1);else
540  #endif
541  _pGetOrder(p1) = _pGetOrder(p2) + _pGetOrder(p3);
542}
543
544// Similar to pCopyAddFast, except that we assume that the component
545// of p2 and p3 is zero component
546#if defined(PDEBUG) && PDEBUG == 1
547#define _pCopyAddFast1(p1, p2, p3)  pDBCopyAddFast(p1, p2, p3, __FILE__, __LINE__)
548extern  void pDBCopyAddFast(poly p1, poly p2, poly p3, char* f, int l);
549inline void __pCopyAddFast1(poly p1, poly p2, poly p3)
550#else
551  DECLARE(void, _pCopyAddFast1(poly p1, poly p2, poly p3))
552#endif // defined(PDEBUG) && PDEBUG == 1
553{
554#ifndef COMP_NO_EXP_VECTOR_OPS
555  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
556  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
557  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
558  const unsigned long* const ub = s3 + pVariables1W;
559#else
560  Exponent_t* s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
561  const Exponent_t* s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
562  const Exponent_t* s3 = (Exponent_t*) &(p3->exp[pVarLowIndex]);
563  const Exponent_t* const ub = s3 + pVariables;
564#ifdef  WORDS_BIGENDIAN
565  *((unsigned long *) ((unsigned long*) p1) + pMonomSizeW -1) = 0;
566#else
567  *((unsigned long *) p1->exp) = 0;
568#endif
569#endif
570
571  for (;;)
572  {
573    *s1 = *s2 + *s3;
574    s3++;
575    if (s3 == ub) break;
576    s1++;
577    s2++;
578  }
579  #ifdef TEST_MAC_ORDER
580  if (bNoAdd) bSetm(p1);else
581  #endif
582  _pGetOrder(p1) = _pGetOrder(p2) + _pGetOrder(p3);
583}
584
585
586#ifndef COMP_NO_EXP_VECTOR_OPS
587
588#if SIZEOF_LONG == 4
589
590#if SIZEOF_EXPONENT == 1
591#define P_DIV_MASK 0x80808080
592#else // SIZEOF_EXPONENT == 2
593#define P_DIV_MASK 0x80008000
594#endif
595
596#else // SIZEOF_LONG == 8
597
598#if SIZEOF_EXPONENT == 1
599#define P_DIV_MASK 0x8080808080808080
600#elif  SIZEOF_EXPONENT == 2
601#define P_DIV_MASK 0x8000800080008000
602#else // SIZEOF_EXPONENT == 4
603#define P_DIV_MASK 0x8000000080000000
604#endif
605
606#endif
607
608DECLARE(BOOLEAN, __pDivisibleBy(poly a, poly b))
609{
610  const unsigned long* s1;
611  const unsigned long* s2;
612  const unsigned long* lb;
613
614#ifdef WORDS_BIGENDIAN
615  lb = (unsigned long*) &(a->exp[0]);
616  if (pVariables & ((SIZEOF_LONG / SIZEOF_EXPONENT) - 1))
617  {
618    // now pVariables == pVariables1W, i.e. there are exponents in the
619    // "first" word of exponentvector
620    s1 = ((unsigned long*) a) + pMonomSizeW -1;
621    s2 = ((unsigned long*) b) + pMonomSizeW -1;
622  }
623  else
624  {
625    // first exponent word has only component as significant field --
626    // Hence, do not bother
627    s1 = ((unsigned long*) a) + pMonomSizeW -2;
628    s2 = ((unsigned long*) b) + pMonomSizeW -2;
629  }
630#else // !WORDS_BIGENDIAN
631  lb = ((unsigned long*) a) + pMonomSizeW;
632  if (pVariables & ((SIZEOF_LONG / SIZEOF_EXPONENT) - 1))
633  {
634    s1 = (unsigned long*) &(a->exp[0]);
635    s2 = (unsigned long*) &(b->exp[0]);
636  }
637  else
638  {
639    s1 = (unsigned long*) &(a->exp[0]) + 1;
640    s2 = (unsigned long*) &(b->exp[0]) + 1;
641  }
642#endif
643  for (;;)
644  {
645// O.K. -- and now comes a bit of magic. The following _really_
646// works. Think about it! If you can prove it, please tell me, for I
647// did not bother to prove it formally (Hint: We can assume that our
648// exponents are always positive).
649    if ((*s2 - *s1) & P_DIV_MASK) return FALSE;
650#ifdef WORDS_BIGENDIAN
651    if (s1 == lb) return TRUE;
652    s1--;
653    s2--;
654#else
655    s1++;
656    if (s1 == lb) return TRUE;
657    s2++;
658#endif
659  }
660}
661
662#else //  ! COMP_NO_EXP_VECTOR_OPS
663
664DECLARE(BOOLEAN, __pDivisibleBy(poly a, poly b))
665{
666#ifdef WORDS_BIGENDIAN
667  const Exponent_t* s1 = &(a->exp[pVarHighIndex]);
668  const Exponent_t* s2 = &(b->exp[pVarHighIndex]);
669  const Exponent_t* lb = s1 - pVariables;
670
671  for (;;)
672  {
673    if (*s1 > *s2) return FALSE;
674    s1--;
675    if (s1 == lb) return TRUE;
676    s2--;
677  }
678
679#else // !WORDS_BIGENDIAN
680  const Exponent_t* s1 = &(a->exp[pVarLowIndex]);
681  const Exponent_t* s2 = &(b->exp[pVarLowIndex]);
682  const Exponent_t* lb = s1 + pVariables;
683
684  for (;;)
685  {
686   if (*s1 > *s2) return FALSE;
687   s1++;
688   if (s1 == lb) return TRUE;
689   s2++;
690  }
691#endif  // WORDS_BIGENDIAN
692}
693
694#endif //  COMP_NO_EXP_VECTOR_OPS
695
696#if defined(PDEBUG) && PDEBUG == 1
697#define _pDivisibleBy(a,b)   pDBDivisibleBy(a, b, __FILE__, __LINE__)
698extern  BOOLEAN pDBDivisibleBy(poly p1, poly p2, char* f, int l);
699inline BOOLEAN _pDivisibleBy_orig(poly a, poly b)
700#else
701inline BOOLEAN _pDivisibleBy(poly a, poly b)
702#endif // defined(PDEBUG) && PDEBUG == 1
703{
704  if ((a!=NULL)&&((_pGetComp(a)==0) || (_pGetComp(a) == _pGetComp(b))))
705  {
706    return __pDivisibleBy(a,b);
707  }
708  return FALSE;
709}
710
711#if defined(PDEBUG) && PDEBUG == 1
712#define _pDivisibleBy1(a,b)   pDBDivisibleBy1(a, b, __FILE__, __LINE__)
713extern  BOOLEAN pDBDivisibleBy1(poly p1, poly p2, char* f, int l);
714inline BOOLEAN _pDivisibleBy1_orig(poly a, poly b)
715#else
716inline BOOLEAN _pDivisibleBy1(poly a, poly b)
717#endif // defined(PDEBUG) && PDEBUG == 1
718{
719  if (_pGetComp(a) == 0 || _pGetComp(a) == _pGetComp(b))
720    return __pDivisibleBy(a,b);
721  return FALSE;
722}
723
724#if defined(PDEBUG) && PDEBUG == 1
725#define _pDivisibleBy2(a,b)   pDBDivisibleBy2(a, b, __FILE__, __LINE__)
726extern  BOOLEAN pDBDivisibleBy2(poly p1, poly p2, char* f, int l);
727#else
728#define _pDivisibleBy2(a,b) __pDivisibleBy(a,b)
729#endif // defined(PDEBUG) && PDEBUG == 1
730
731
732DECLARE(BOOLEAN, _pEqual(poly p1, poly p2))
733{
734#ifndef COMP_NO_EXP_VECTOR_OPS
735  const long *s1 = (long*) &(p1->exp[0]);
736  const long *s2 = (long*) &(p2->exp[0]);
737  const long* const lb = s1 + pVariables1W;
738#else
739  const Exponent_t *s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
740  const Exponent_t *s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
741  const Exponent_t* const lb = s1 + pVariables;
742  if (_pGetComp(p1) != _pGetComp(p2)) return FALSE;
743#endif
744
745  for(;;)
746  {
747    if (*s1 != *s2) return FALSE;
748    s1++;
749    if (s1 == lb) return TRUE;
750    s2++;
751  }
752}
753
754inline void _pGetExpV(poly p, Exponent_t *ev)
755{
756  if (_pHasReverseExp)
757  {
758    for (int i = pVarLowIndex, j = pVariables; j; i++, j--)
759      ev[j] = p->exp[i];
760  }
761  else
762    memcpy(&ev[1], &(p->exp[pVarLowIndex]), pVariables*sizeof(Exponent_t));
763  ev[0] = _pGetComp(p);
764}
765
766extern pSetmProc pSetm;
767inline void _pSetExpV(poly p, Exponent_t *ev)
768{
769  if (_pHasReverseExp)
770  {
771    for (int i = pVarLowIndex, j = pVariables; j; i++, j--)
772      p->exp[i] = ev[j];
773  }
774  else
775    memcpy(&(p->exp[pVarLowIndex]), &ev[1], pVariables*sizeof(Exponent_t));
776  _pSetComp(p, ev[0]);
777  pSetm(p);
778}
779
780#else // ! COMP_FAST
781
782DECLARE(BOOLEAN, _pDivisibleBy(poly a, poly b))
783{
784  if ((a!=NULL)&&((a->exp[0]==0) || (a->exp[0] == b->exp[0])))
785  {
786    int i=pVariables;
787    short *e1=&(a->exp[1]);
788    short *e2=&(b->exp[1]);
789    if ((*e1) > (*e2)) return FALSE;
790    do
791    {
792      i--;
793      if (i == 0) return TRUE;
794      e1++;
795      e2++;
796    } while ((*e1) <= (*e2));
797  }
798  return FALSE;
799}
800
801#define _pDivisibleBy1(a,b) _pDivisibleBy(a,b)
802#define _pDivisibleBy2(a,b) _pDivisibleBy(a,b)
803
804#ifdef TEST_MAC_ORDER
805DECLARE(void, pMonAddFast(poly a, poly m))
806{
807  for(int ii =pVariables; ii; ii--) (a)->exp[ii] += (m)->exp[ii];\
808  if (bNoAdd) bSetm(a); else
809  _pGetOrder(a) += _pGetOrder(m);
810}
811#else
812DECLARE(void, pMonAddFast(poly a, poly m))
813{
814  for(int ii =pVariables; ii; ii--) (a)->exp[ii] += (m)->exp[ii];\
815  _pGetOrder(a) += _pGetOrder(m);
816}
817#endif
818
819DECLARE(BOOLEAN, _pEqual(poly p1, poly p2))
820{
821  int i;
822  short *e1=p1->exp;
823  short *e2=p2->exp;
824
825  if (p1->Order != p2->Order) return FALSE;
826  for (i=pVariables; i>=0; i--,e1++,e2++)
827    if (*e1 != *e2) return FALSE;
828  return TRUE;
829}
830
831#define _pGetExpV(p,e)    memcpy((e),(p)->exp,(pVariables+1)*sizeof(short));
832//void    pSetExpV(poly p, short * exp);
833#define _pSetExpV(p,e) {memcpy((p)->exp,(e),(pVariables+1)*sizeof(short));pSetm(p);}
834
835#endif // COMP_FAST
836
837/***************************************************************
838 *
839 * Routines which implement low-level manipulations/operations on exponents
840 *
841 ***************************************************************/
842
843DECLARE(int, __pExpQuerSum2(poly p, int from, int to))
844{
845  int j = p->exp[from];
846  int i = from + 1;
847
848  for(;;)
849  {
850    if (i > to) return j;
851    j += p->exp[i];
852    i++;
853  }
854}
855
856#ifdef COMP_FAST
857#define _pExpQuerSum(p)  __pExpQuerSum2(p, pVarLowIndex, pVarHighIndex)
858
859#define _pExpQuerSum1(p,to) \
860 (_pHasReverseExp ? \
861    __pExpQuerSum2(p, _pExpIndex(to), _pExpIndex(1)) : \
862    __pExpQuerSum2(p, _pExpIndex(1), _pExpIndex(to)))
863
864#define _pExpQuerSum2(p,from,to) \
865 (_pHasReverseExp ? \
866    __pExpQuerSum2(p, _pExpIndex(to), _pExpIndex(from)) : \
867    __pExpQuerSum2(p, _pExpIndex(from), _pExpIndex(to)))
868#else
869
870#define _pExpQuerSum(p)             __pExpQuerSum2(p, 1, pVariables)
871#define _pExpQuerSum1(p, to)        __pExpQuerSum2(p, 1, to)
872#define _pExpQuerSum2(p, from, to)  __pExpQuerSum2(p, from, to)
873
874#endif
875
876/***************************************************************
877 *
878 * Routines which implement macaulay ordering routines
879 *
880 ***************************************************************/
881#ifdef TEST_MAC_ORDER
882
883DECLARE(void, _bSetm0(poly p))
884{
885
886  int i=1;
887  int ord = -INT_MAX;
888  Exponent_t *ep;
889
890  if(_pHasReverseExp)
891  {
892    ep=&(p->exp[pVarHighIndex]);
893    int *ip=bBinomials+(*ep); /*_pGetExp(p,1);*/
894    loop
895    {
896      ord += (*ip);
897      if (i==pVariables) break;
898      i++;
899      //ip+=bHighdeg_1+_pGetExp(p,i);
900      ep--;
901      ip+=bHighdeg_1+(*ep);
902    }
903  }
904  else
905  {
906    ep=&(p->exp[pVarLowIndex]);
907    int *ip=bBinomials+(*ep); /*_pGetExp(p,1);*/
908    loop
909    {
910      ord += (*ip);
911      if (i==pVariables) break;
912      i++;
913      //ip+=bHighdeg_1+_pGetExp(p,i);
914      ep++;
915      ip+=bHighdeg_1+(*ep);
916    }
917  }
918  p->Order=ord;
919}
920
921DECLARE(void, _bSetm(poly p))
922{
923  int ord = _pExpQuerSum(p);
924
925  if (ord<bHighdeg)
926    _bSetm0(p);
927  else
928    p->Order=ord;
929}
930
931// ordering dp,c or c,dp, general case
932#if defined(PDEBUG) && PDEBUG == 1
933#define pbMonAddFast(p1, p2)  pDBMonAddFast(p1, p2, __FILE__, __LINE__)
934extern  void pbDBMonAddFast(poly p1, poly p2, char* f, int l);
935inline void _pbMonAddFast(poly p1, poly p2)
936#else
937#define pbMonAddFast(p1, p2)  _pbMonAddFast(p1, p2)
938DECLARE(void, _pbMonAddFast(poly p1, poly p2))
939#endif // defined(PDEBUG) && PDEBUG == 1
940{
941  // OK -- this might be the only place where we are kind of quick and
942  // dirty: the following only works correctly if all exponents are
943  // positive and the sum of two exponents does not exceed
944  // EXPONENT_MAX
945#ifndef COMP_NO_EXP_VECTOR_OPS
946  Exponent_t c2 = _pGetComp(p2);
947  int i = pVariables1W;
948  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
949  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
950  // set comp of p2 temporarily to 0, so that nothing is added to comp of p1
951  _pSetComp(p2, 0);
952#else
953  int i = pVariables;
954  Exponent_pt s1 = &(p1->exp[pVarLowIndex]);
955  Exponent_pt s2 = &(p2->exp[pVarLowIndex]);
956#endif
957
958  for (;;)
959  {
960    *s1 += *s2;
961    i--;
962    if (i==0) break;
963    s1++;
964    s2++;
965  }
966#ifndef COMP_NO_EXP_VECTOR_OPS
967  // reset comp of p2
968  _pSetComp(p2, c2);
969#endif
970  if ((_pGetOrder(p1)|_pGetOrder(p2))>0) // i.e. overflow of mac order
971    _pGetOrder(p1) += _pGetOrder(p2);
972  else
973    _bSetm(p1);
974}
975
976// ordering dp,c or c,dp, below degree limit
977#if defined(PDEBUG) && PDEBUG == 1
978#define pbMonAddFast0(p1, p2)  pbDBMonAddFast0(p1, p2, __FILE__, __LINE__)
979extern  void pbDBMonAddFast0(poly p1, poly p2, char* f, int l);
980inline void _pbMonAddFast0(poly p1, poly p2)
981#else
982  DECLARE(void, pbMonAddFast0(poly p1, poly p2))
983#endif // defined(PDEBUG) && PDEBUG == 1
984{
985  // OK -- this might be the only place where we are kind of quick and
986  // dirty: the following only works correctly if all exponents are
987  // positive and the sum of two exponents does not exceed
988  // EXPONENT_MAX
989#ifndef COMP_NO_EXP_VECTOR_OPS
990  Exponent_t c2 = _pGetComp(p2);
991  int i = pVariables1W;
992  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
993  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
994  // set comp of p2 temporarily to 0, so that nothing is added to comp of p1
995  _pSetComp(p2, 0);
996#else
997  int i = pVariables;
998  Exponent_pt s1 = &(p1->exp[pVarLowIndex]);
999  Exponent_pt s2 = &(p2->exp[pVarLowIndex]);
1000#endif
1001
1002  for (;;)
1003  {
1004    *s1 += *s2;
1005    i--;
1006    if (i==0) break;
1007    s1++;
1008    s2++;
1009  }
1010#ifndef COMP_NO_EXP_VECTOR_OPS
1011  // reset comp of p2
1012  _pSetComp(p2, c2);
1013#endif
1014  _bSetm0(p1);
1015}
1016
1017// ordering dp,c or c,dp, below degree limit
1018// Makes p1 a copy of p2 and adds on exponets of p3
1019#if defined(PDEBUG) && PDEBUG == 1
1020#define _pbCopyAddFast0(p1, p2, p3)  pDBCopyAddFast(p1, p2, p3, __FILE__, __LINE__)
1021inline void __pbCopyAddFast0(poly p1, poly p2, poly p3)
1022#else
1023  DECLARE(void, _pbCopyAddFast0(poly p1, poly p2, poly p3))
1024#endif // defined(PDEBUG) && PDEBUG == 1
1025{
1026  p1->next = p2->next;
1027  p1->coef = p2->coef;
1028
1029#ifndef COMP_NO_EXP_VECTOR_OPS
1030  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
1031  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
1032  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
1033  const unsigned long* const ub = s3 + pVariables1W;
1034#else
1035  Exponent_t* s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
1036  const Exponent_t* s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
1037  const Exponent_t* s3 = (Exponent_t*) &(p3->exp[pVarLowIndex]);
1038  const Exponent_t* const ub = s3 + pVariables;
1039// need to zero the "fill in" slots (i.e., empty exponents)
1040#ifdef  WORDS_BIGENDIAN
1041  *((unsigned long*) p1 + pMonomSizeW -1) = 0;
1042#else
1043  *((unsigned long *) p1->exp) = 0;
1044#endif
1045#endif
1046
1047  for (;;)
1048  {
1049    *s1 = *s2 + *s3;
1050    s3++;
1051    if (s3 == ub) break;
1052    s1++;
1053    s2++;
1054  }
1055  // we first are supposed to do a copy from p2 to p1 -- therefore,
1056  // component of p1 is set to comp of p2
1057  _pSetComp(p1, _pGetComp(p2));
1058  _bSetm0(p1);
1059}
1060
1061// Similar to pCopyAddFast, except that we assume that the component
1062// of p2 and p3 is zero component
1063#if defined(PDEBUG) && PDEBUG == 1
1064#define _pbCopyAddFast10(p1, p2, p3)  pbDBCopyAddFast0(p1, p2, p3, __FILE__, __LINE__)
1065extern  void pbDBCopyAddFast0(poly p1, poly p2, poly p3, char* f, int l);
1066inline void __pbCopyAddFast10(poly p1, poly p2, poly p3)
1067#else
1068  DECLARE(void, _pbCopyAddFast10(poly p1, poly p2, poly p3))
1069#endif // defined(PDEBUG) && PDEBUG == 1
1070{
1071  p1->next = p2->next;
1072  p1->coef = p2->coef;
1073
1074#ifndef COMP_NO_EXP_VECTOR_OPS
1075  unsigned long* s1 = (unsigned long*) &(p1->exp[0]);
1076  const unsigned long* s2 = (unsigned long*) &(p2->exp[0]);
1077  const unsigned long* s3 = (unsigned long*) &(p3->exp[0]);
1078  const unsigned long* const ub = s3 + pVariables1W;
1079#else
1080  Exponent_t* s1 = (Exponent_t*) &(p1->exp[pVarLowIndex]);
1081  const Exponent_t* s2 = (Exponent_t*) &(p2->exp[pVarLowIndex]);
1082  const Exponent_t* s3 = (Exponent_t*) &(p3->exp[pVarLowIndex]);
1083  const Exponent_t* const ub = s3 + pVariables;
1084#ifdef  WORDS_BIGENDIAN
1085  *((unsigned long*) p1 + pMonomSizeW -1) = 0;
1086#else
1087  *((unsigned long *) p1->exp) = 0;
1088#endif
1089#endif
1090
1091  for (;;)
1092  {
1093    *s1 = *s2 + *s3;
1094    s3++;
1095    if (s3 == ub) break;
1096    s1++;
1097    s2++;
1098  }
1099  _bSetm0(p1);
1100}
1101
1102#endif
1103
1104#endif // POLYS_IMPL_H
1105
Note: See TracBrowser for help on using the repository browser.