Context Navigation

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

spielwiese

Last change on this file since c860e9 was c860e9, checked in by Hans Schönemann <hannes@…>, 25 years ago
* hannes: introduced EXPONENT_MAX in polys-impl.h fixed exponent overflow in pPower * GMG: fixes in general.lib git-svn-id: file:///usr/local/Singular/svn/trunk@3528 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 16.3 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.31 1999-08-23 13:15: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
34	typedef Exponent_t monomial[VARS + 1];
35	typedef Exponent_t* Exponent_pt;
36
37	typedef long Order_t;
38	struct 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
57	extern int pVariables;
58	// size of a monom in bytes - always a multiple of sizeof(void*)
59	extern int pMonomSize;
60	// size of a monom in units of sizeof(void*) -- i.e. in words
61	extern int pMonomSizeW;
62	// Ceiling((pVariables+1) / sizeof(void*)) == length of exp-vector in words
63	extern int pVariables1W;
64	// Ceiling((pVariables) / sizeof(void*))
65	extern int pVariablesW;
66	extern int *pVarOffset;
67	extern int pVarLowIndex;
68	extern int pVarHighIndex;
69	extern 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
83	extern void pGetVarIndicies_Lex(int nvars, int* VarOffset, int &VarCompIndex,
84	int &VarLowIndex, int &VarHighIndex);
85	// for simple, revlex orderings
86	extern void pGetVarIndicies_RevLex(int nvars,int *VarOffset,int &VarCompIndex,
87	int &VarLowIndex, int &VarHighIndex);
88	// for all non-simple orderings
89	extern 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
94	extern 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
112	extern Exponent_t pPDSetExp(poly p, int v, Exponent_t e, char* f, int l);
113	extern Exponent_t pPDGetExp(poly p, int v, char* f, int l);
114	extern Exponent_t pPDIncrExp(poly p, int v, char* f, int l);
115	extern Exponent_t pPDDecrExp(poly p, int v, char* f, int l);
116	extern Exponent_t pPDAddExp(poly p, int v, Exponent_t e, char* f, int l);
117	extern Exponent_t pPDMultExp(poly p, int v, Exponent_t e, char* f, int l);
118	extern Exponent_t pPDSubExp(poly p, int v, Exponent_t e, char* f, int l);
119
120	extern Exponent_t pPDRingSetExp(ring r,poly p,int v,Exponent_t e,char* f,int l);
121	extern 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
149	inline Exponent_t _pGetExpSum(poly p1, poly p2, int i)
150	{
151	int index = _pExpIndex(i);
152	return p1->exp[index] + p2->exp[index];
153	}
154	inline 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
170	inline 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
178	extern pSetmProc pSetm;
179	inline 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
195	poly pDBInit(char *f, int l);
196	poly pDBCopy(poly a, char *f, int l);
197	poly pDBCopy1(poly a, char *f, int l);
198	poly pDBHead(poly a, char *f, int l);
199	poly pDBHead0(poly a, char *f, int l);
200	poly pDBFetchCopy(ring r, poly a, char *f, int l);
201
202	void pDBDelete(poly * a, char * f, int l);
203	void 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
222	inline poly _pNew()
223	{
224	void * p;
225	AllocHeap(p, mm_specHeap);
226	return (poly)p;
227	}
228
229	inline poly _pInit(void)
230	{
231	poly p = _pNew();
232	memsetW((long*)p,0, pMonomSizeW);
233	return p;
234	}
235
236	extern void _pDelete(poly * a);
237	extern void _pDelete1(poly * a);
238
239	extern poly _pCopy(poly a);
240	extern poly _pCopy1(poly a);
241	extern poly _pHead(poly a);
242	extern poly _pHead0(poly a);
243	extern 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__)
276	extern void pDBMonAddFast(poly p1, poly p2, char* f, int l);
277	inline 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__)
309	extern void pDBCopyAddFast(poly p1, poly p2, poly p3, char* f, int l);
310	inline 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__)
340	extern void pDBCopyAddFast(poly p1, poly p2, poly p3, char* f, int l);
341	inline 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__)
370	extern void pDBCopyAddFastHomog(poly p1, poly p2, poly p3, Order_t Order,
371	char* f, int l);
372	inline 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	#define EXPONENT_MAX 0x7f
401	#else // SIZEOF_EXPONENT == 2
402	#define P_DIV_MASK 0x80008000
403	#define EXPONENT_MAX 0x7fff
404	#endif
405
406	#else // SIZEOF_LONG == 8
407
408	#if SIZEOF_EXPONENT == 1
409	#define P_DIV_MASK 0x8080808080808080
410	#define EXPONENT_MAX 0x7f
411	#elif SIZEOF_EXPONENT == 2
412	#define P_DIV_MASK 0x8000800080008000
413	#define EXPONENT_MAX 0x7fff
414	#else // SIZEOF_EXPONENT == 4
415	#define P_DIV_MASK 0x8000000080000000
416	#define EXPONENT_MAX 0x7fffffff
417	#endif
418
419	#endif
420
421	DECLARE(BOOLEAN, __pDivisibleBy(poly a, poly b))
422	{
423	#ifdef WORDS_BIGENDIAN
424	const unsigned long* const lb = (unsigned long*) &(a->exp[0]);;
425	const unsigned long* s1 = ((unsigned long*) a) + pMonomSizeW -1;
426	const unsigned long* s2 = ((unsigned long*) b) + pMonomSizeW -1;
427	#else
428	const unsigned long* const lb = ((unsigned long*) a) + pMonomSizeW;
429	const unsigned long* s1 = (unsigned long*) &(a->exp[0]);
430	const unsigned long* s2 = (unsigned long*) &(b->exp[0]);
431	#endif
432
433	for (;;)
434	{
435	// Yes, the following is correct, provided that the exponents do
436	// not have their first bit set
437	if ((s2 - s1) & P_DIV_MASK) return FALSE;
438	#ifdef WORDS_BIGENDIAN
439	if (s1 == lb) return TRUE;
440	s1--;
441	s2--;
442	#else
443	s1++;
444	if (s1 == lb) return TRUE;
445	s2++;
446	#endif
447	}
448	}
449
450	#if defined(PDEBUG) && PDEBUG == 1
451	#define _pDivisibleBy(a,b) pDBDivisibleBy(a, b, __FILE__, __LINE__)
452	extern BOOLEAN pDBDivisibleBy(poly p1, poly p2, char* f, int l);
453	inline BOOLEAN _pDivisibleBy_orig(poly a, poly b)
454	#else
455	inline BOOLEAN _pDivisibleBy(poly a, poly b)
456	#endif // defined(PDEBUG) && PDEBUG == 1
457	{
458	if ((a!=NULL)&&((_pGetComp(a)==0) \|\| (_pGetComp(a) == _pGetComp(b))))
459	{
460	return __pDivisibleBy(a,b);
461	}
462	return FALSE;
463	}
464
465	#if defined(PDEBUG) && PDEBUG == 1
466	#define _pDivisibleBy1(a,b) pDBDivisibleBy1(a, b, __FILE__, __LINE__)
467	extern BOOLEAN pDBDivisibleBy1(poly p1, poly p2, char* f, int l);
468	inline BOOLEAN _pDivisibleBy1_orig(poly a, poly b)
469	#else
470	inline BOOLEAN _pDivisibleBy1(poly a, poly b)
471	#endif // defined(PDEBUG) && PDEBUG == 1
472	{
473	if (_pGetComp(a) == 0 \|\| _pGetComp(a) == _pGetComp(b))
474	return __pDivisibleBy(a,b);
475	return FALSE;
476	}
477
478	#if defined(PDEBUG) && PDEBUG == 1
479	#define _pDivisibleBy2(a,b) pDBDivisibleBy2(a, b, __FILE__, __LINE__)
480	extern BOOLEAN pDBDivisibleBy2(poly p1, poly p2, char* f, int l);
481	#else
482	#define _pDivisibleBy2(a,b) __pDivisibleBy(a,b)
483	#endif // defined(PDEBUG) && PDEBUG == 1
484
485
486	DECLARE(BOOLEAN, _pEqual(poly p1, poly p2))
487	{
488	const long s1 = (long) &(p1->exp[0]);
489	const long s2 = (long) &(p2->exp[0]);
490	const long* const lb = s1 + pVariables1W;
491
492	for(;;)
493	{
494	if (s1 != s2) return FALSE;
495	s1++;
496	if (s1 == lb) return TRUE;
497	s2++;
498	}
499	}
500
501	/***************************************************************
502	*
503	* Routines which implement low-level manipulations/operations
504	* on exponents and "are allowed" to access single exponetns
505	*
506	***************************************************************/
507
508	DECLARE(int, __pExpQuerSum2(poly p, int from, int to))
509	{
510	int j = p->exp[from];
511	int i = from + 1;
512
513	for(;;)
514	{
515	if (i > to) return j;
516	j += p->exp[i];
517	i++;
518	}
519	}
520
521	#define _pExpQuerSum(p) __pExpQuerSum2(p, pVarLowIndex, pVarHighIndex)
522
523	inline int _pExpQuerSum1(poly p, int to)
524	{
525	int ei_to = _pExpIndex(to);
526	int ei_1 = _pExpIndex(1);
527
528	if (ei_1 > ei_to)
529	return __pExpQuerSum2(p, ei_to, ei_1);
530	else
531	return __pExpQuerSum2(p, ei_1, ei_to);
532	}
533
534
535	inline int _pExpQuerSum2(poly p,int from,int to)
536	{
537	int ei_to = _pExpIndex(to);
538	int ei_from = _pExpIndex(from);
539
540	if (ei_from > ei_to)
541	return __pExpQuerSum2(p, ei_to, ei_from);
542	else
543	return __pExpQuerSum2(p, ei_from, ei_to);
544	}
545
546	#endif // POLYS_IMPL_H

Note: See TracBrowser for help on using the repository browser.

Download in other formats: