Context Navigation

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

spielwiese

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

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

Download in other formats: