Context Navigation

source: git/Singular/polys.cc @ 954622

spielwiese

Last change on this file since 954622 was 954622, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: fixes for "TEST_MAC_ORDER" and some optimizations for binomials (binom.cc binom.h claptmpl.cc ipid.h kstd1.cc kstd2.cc kutil.cc polys-impl.h polys.cc polys.h polys0.cc polys1.cc ring.h spolys.cc spolys0.h) git-svn-id: file:///usr/local/Singular/svn/trunk@1005 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 68.9 KB

Line
1
2	/****************************************
3	* Computer Algebra System SINGULAR *
4	****************************************/
5	/* $Id: polys.cc,v 1.12 1998-01-05 16:39:26 Singular Exp $ */
6
7	/*
8	* ABSTRACT - all basic methods to manipulate polynomials
9	*/
10
11	/* includes */
12	#include <stdio.h>
13	#include <string.h>
14	#include <ctype.h>
15	#include "mod2.h"
16	#include "tok.h"
17	#include "mmemory.h"
18	#include "febase.h"
19	#include "numbers.h"
20	#include "polys.h"
21	#include "ring.h"
22	#include "binom.h"
23
24	#ifdef COMP_FAST
25	#include "polys-comp.h"
26	#endif
27
28	/* ----------- global variables, set by pChangeRing --------------------- */
29	/* initializes the internal data from the exp vector */
30	pSetmProc pSetm;
31	/* computes length and maximal degree of a POLYnomial */
32	pLDegProc pLDeg;
33	/* computes the degree of the initial term, used for std */
34	pFDegProc pFDeg;
35	/* the monomial ordering of the head monomials a and b */
36	pCompProc pComp0;
37	/* returns -1 if a comes before b, 0 if a=b, 1 otherwise */
38
39	/* the number of variables */
40	/* 1 for polynomial ring, -1 otherwise */
41	int pOrdSgn;
42	/* TRUE for momomial output as x2y, FALSE for x^2y /
43	int pShortOut = (int)TRUE;
44	// it is of type int, not BOOLEAN because it is also in ip
45	/* TRUE if the monomial ordering is not compatible with pFDeg */
46	BOOLEAN pLexOrder;
47	/* TRUE if the monomial ordering has polynomial and power series blocks */
48	BOOLEAN pMixedOrder;
49
50	#if defined(TEST_MAC_ORDER) \|\| defined(COMP_FAST)
51	int pComponentOrder;
52	#else
53	static int pComponentOrder;
54	#endif
55
56	#ifdef DRING
57	int p2;
58	BOOLEAN pDRING=FALSE;
59	#endif
60
61	#ifdef SRING
62	int pAltVars;
63	BOOLEAN pSRING=FALSE;
64	#endif
65
66	#ifdef SDRING
67	BOOLEAN pSDRING=FALSE;
68	#include "polys.inc"
69	#endif
70
71	/* ----------- global variables, set by procedures from hecke/kstd1 ----- */
72	/* the highest monomial below pHEdge */
73	poly ppNoether = NULL;
74
75	/* -------------- static variables --------------------------------------- */
76	/is the basic comparing procedure during a computation of syzygies/
77	static pCompProc pCompOld;
78	/for grouping module indecees during computations/
79	#ifndef COMP_FAST
80	static int maxBound = 0;
81	#else
82	int maxBound = 0;
83	#endif
84
85	/contains the headterms for the Schreyer orderings/
86	static int* SchreyerOrd;
87	static int maxSchreyer=0;
88	static int indexShift=0;
89	static pLDegProc pLDegOld;
90
91	typedef int (bcompProc)(poly p1, poly p2, int i1, int i2, short w);
92	static bcompProc bcomph[20];
93	static short** polys_wv;
94	static short * firstwv;
95	static int * block0;
96	static int * block1;
97	static int firstBlockEnds;
98	static int * order;
99
100	/0 implementation/
101	/-------- the several possibilities for pSetm:-----------------------/
102
103	/* Remark: These could be made more efficient by avoiding using pGetExp */
104
105	/*2
106	* define the order of p with respect to lex. ordering, N=1
107	*/
108	static void setlex1(poly p)
109	{
110	p->Order = (Order_t)pGetExp(p,1);
111	}
112
113	/*2
114	* define the order of p with respect to lex. ordering, N>1
115	*/
116	static void setlex2(poly p)
117	{
118	p->Order = (((Order_t)pGetExp(p,1))<<(sizeof(short)*8))
119	+ (Order_t)pGetExp(p,2);
120	}
121
122	/*2
123	* define the order of p with respect to a degree ordering
124	*/
125	static void setdeg1(poly p)
126	{
127	p->Order = pExpQuerSum1(p, firstBlockEnds);
128	#ifdef COMP_DEBUG
129	int i, j = pGetExp(p,1);
130
131	for (i = firstBlockEnds; i>1; i--) j += pGetExp(p,i);
132	if (j != p->Order)
133	{
134	Print("Error in setdeg1");
135	pExpQuerSum1(p, firstBlockEnds);
136	}
137	#endif
138	}
139
140	/*2
141	* define the order of p with respect to a degree ordering
142	* with weigthts
143	*/
144	static void setdeg1w(poly p)
145	{
146	Order_t i, j = 0;
147
148	for (i = firstBlockEnds; i>0; i--)
149	{
150	j += ((Order_t) pGetExp(p,i))*firstwv[i-1];
151	}
152	p->Order = j;
153	}
154
155
156	/-------- IMPLEMENTATION OF MONOMIAL COMPARISONS ---------------------/
157
158
159	/***************************************************************
160	***************************************************************
161	*
162	* MONOMIAL COMPARIONS:
163	* They are influenced by the following macros:
164	* COMP_TRADITIONAL -- (un)set in polys-impl.h
165	Keeps the traditional comparison routines
166	defined -- needed as long as their might be comparisons with
167	negativ components.
168	All the traditional routines are prefixed by t_
169
170	* COMP_FAST -- (un)set in polys-impl.h
171	Implements monomial comparisons using the fast
172	techniques. Undefine in case there are problems. All the fast
173	routines are prefixed by f_
174
175	* COMP_STATISTIC -- (un)set in polys-impl.h
176	Provides several routines for accumulating statistics on monomial
177	comparisons
178
179	* COMP_DEBUG -- (un)set in polys-impl.h
180	Turns on debugging of COMP_FAST by comparing the resutsl of fast
181	comparison with traditional comparison
182
183	*
184	*
185	***************************************************************
186	***************************************************************/
187
188	#ifdef COMP_DEBUG
189	static int debug_comp(poly p1, poly p2);
190	#endif // COMP_DEBUG
191
192	#ifdef COMP_STATISTICS
193	// Initializes, resets and outputs Comp statistics
194	void InitCompStatistics(int nvars);
195	void ResetCompStatistics();
196	void OutputCompStatistics();
197	unsigned long MonomCountTotal = 0;
198	unsigned long MonomCountOrderS = 0;
199	unsigned long MonomCountOrderG = 0;
200	unsigned long* MonomCountExp = NULL;
201	unsigned long MonomCountExpS = 0;
202	unsigned long MonomCountExpG = 0;
203	unsigned long MonomCountCompS = 0;
204	unsigned long MonomCountCompG = 0;
205	unsigned long MonomCountEqual = 0;
206	#endif // COMP_STATISTICS
207
208
209	#define NonZeroR(l, actionG, actionS) \
210	do \
211	{ \
212	long _l = l; \
213	if (_l) \
214	{ \
215	if (_l > 0) actionG; \
216	actionS; \
217	} \
218	} \
219	while(0)
220
221	/***************************************************************
222	*
223	* FAST COMPARISONS
224	*
225	***************************************************************/
226	#ifdef COMP_FAST
227	static pCompProc f_pComp0 = NULL;
228	static int f_comp_1(poly p1, poly p2);
229	static int f_comp_c_1(poly p1, poly p2);
230	static int f_comp_1_c(poly p1, poly p2);
231	static int f_comp_2(poly p1, poly p2);
232	static int f_comp_c_2(poly p1, poly p2);
233	static int f_comp_2_c(poly p1, poly p2);
234	static int f_comp_2i(poly p1, poly p2);
235	static int f_comp_c_2i(poly p1, poly p2);
236	static int f_comp_2i_c(poly p1, poly p2);
237	static int f_comp_2i_1(poly p1, poly p2);
238	static int f_comp_c_2i_1(poly p1, poly p2);
239	static int f_comp_2i_1_c(poly p1, poly p2);
240
241	#define Mreturn(d, multiplier) \
242	{ \
243	if (d > 0) return multiplier; \
244	return -multiplier; \
245	}
246
247	// comp_1 is used if pVariables1W == 1 and component is compatible with ordering
248	static int f_comp_1(poly p1, poly p2)
249	{
250	register long d = pGetOrder(p1) - pGetOrder(p2);
251
252	if (d) Mreturn(d, pOrdSgn);
253	_pMonCmp_1(p1, p2, d, goto NotEqual, return 0);
254
255	NotEqual:
256	Mreturn(d, pLexSgn);
257	}
258
259	// comp_1_c is used if pVariables1W == 1, priority is given to exponents,
260	// component is incompatible with ordering
261	static int f_comp_1_c(poly p1, poly p2)
262	{
263	register long d = pGetOrder(p1) - pGetOrder(p2);
264
265	if (d) Mreturn(d, pOrdSgn);
266	_pMonCmp_1_c(p1, p2, d, goto NotEqual , return 0);
267
268	NotEqual:
269	Mreturn(d, pLexSgn)
270	}
271
272	// comp_1_c is used if pVariables1W == 1, priority is given to component,
273	// component is incompatible with ordering
274	static int f_comp_c_1(poly p1, poly p2)
275	{
276	register long d = pGetComp(p2) - pGetComp(p1);
277	if (d) Mreturn(d, pComponentOrder);
278	d = pGetOrder(p1) - pGetOrder(p2);
279	if (d) Mreturn(d, pOrdSgn);
280	_pMonCmp_1(p1, p2, d, goto NotEqual , return 0);
281
282	NotEqual:
283	Mreturn(d, pLexSgn);
284	}
285
286	// comp_2 is used if pVariables1W == 2 and component is compatible with ordering
287	static int f_comp_2(poly p1, poly p2)
288	{
289	register long d = pGetOrder(p1) - pGetOrder(p2);
290
291	if (d) Mreturn(d, pOrdSgn);
292	_pMonCmp_2(p1, p2, d, goto NotEqual , return 0);
293
294	NotEqual:
295	Mreturn(d, pLexSgn);
296	}
297
298	// comp_2_c is used if pVariables1W == 2, priority is given to exponents,
299	// component is incompatible with ordering
300	static int f_comp_2_c(poly p1, poly p2)
301	{
302	register long d = pGetOrder(p1) - pGetOrder(p2);
303
304	if (d) Mreturn(d, pOrdSgn);
305	_pMonCmp_2_c(p1, p2, d, goto NotEqual, return 0);
306
307	NotEqual:
308	Mreturn(d, pLexSgn);
309	}
310
311	// comp_2_c is used if pVariables1W == 2, priority is given to component,
312	// component is incompatible with ordering
313	static int f_comp_c_2(poly p1, poly p2)
314	{
315	register long d = pGetComp(p2) - pGetComp(p1);
316	if (d) Mreturn(d, pComponentOrder);
317	d = pGetOrder(p1) - pGetOrder(p2);
318	if (d) Mreturn(d, pOrdSgn);
319	_pMonCmp_2(p1, p2, d, goto NotEqual , return 0);
320
321	NotEqual:
322	Mreturn(d, pLexSgn);
323	}
324
325	// comp_2i is used if pVariables1W == 2*i and component is compatible
326	// with ordering
327	static int f_comp_2i(poly p1, poly p2)
328	{
329	register long d = pGetOrder(p1) - pGetOrder(p2);
330
331	if (d) Mreturn(d, pOrdSgn);
332	_pMonCmp_2i(p1, p2, pVariables1W, d, goto NotEqual , return 0);
333
334	NotEqual:
335	Mreturn(d, pLexSgn);
336	}
337
338	// comp_2i_c is used if pVariables1W == 2*i, priority is given to exponents,
339	// component is incompatible with ordering
340	static int f_comp_2i_c(poly p1, poly p2)
341	{
342	register long d = pGetOrder(p1) - pGetOrder(p2);
343
344	if (d) Mreturn(d, pOrdSgn);
345	_pMonCmp_2i_c(p1, p2, pVariables1W, d, goto NotEqual , return 0);
346
347	NotEqual:
348	Mreturn(d, pLexSgn);
349	}
350
351	// comp_2i_c is used if pVariables1W == 2*i, priority is given to component,
352	// component is incompatible with ordering
353	static int f_comp_c_2i(poly p1, poly p2)
354	{
355	register long d = pGetComp(p2) - pGetComp(p1);
356	if (d) Mreturn(d, pComponentOrder);
357	d = pGetOrder(p1) - pGetOrder(p2);
358	if (d) Mreturn(d, pOrdSgn);
359	_pMonCmp_2i(p1, p2, pVariablesW, d, goto NotEqual, return 0);
360
361	NotEqual:
362	Mreturn(d, pLexSgn);
363	}
364
365	// comp_2i_1 is used if pVariables1W == 2*i and component is compatible
366	// with ordering
367	static int f_comp_2i_1(poly p1, poly p2)
368	{
369	register long d = pGetOrder(p1) - pGetOrder(p2);
370
371	if (d) Mreturn(d, pOrdSgn);
372	_pMonCmp_2i_1(p1, p2, pVariables1W, d, goto NotEqual, return 0);
373
374	NotEqual:
375	Mreturn(d, pLexSgn);
376	}
377
378	// comp_2i_1_c is used if pVariables1W == 2*i, priority is given to exponents,
379	// component is incompatible with ordering
380	static int f_comp_2i_1_c(poly p1, poly p2)
381	{
382	register long d = pGetOrder(p1) - pGetOrder(p2);
383	if (d)
384	{
385	if (d > 0) return pOrdSgn;
386	return -pOrdSgn;
387	}
388
389	_pMonCmp_2i_1_c(p1, p2, pVariables1W, d, goto NotEqual , return 0);
390
391	NotEqual:
392	Mreturn(d, pLexSgn);
393	}
394
395	// comp_c_2i_1 is used if pVariables1W == 2*i, priority is given to component,
396	// component is incompatible with ordering
397	static int f_comp_c_2i_1(poly p1, poly p2)
398	{
399	register long d = pGetComp(p2) - pGetComp(p1);
400	if (d) Mreturn(d, pComponentOrder);
401	d = pGetOrder(p1) - pGetOrder(p2);
402	if (d) Mreturn(d, pOrdSgn);
403	_pMonCmp_2i_1(p1, p2, pVariablesW, d, goto NotEqual, return 0);
404
405	NotEqual:
406	Mreturn(d, pLexSgn);
407	}
408	#endif // COMP_FAST
409
410	/***************************************************************
411	*
412	* FAST COMPARISONS
413	*
414	***************************************************************/
415	#ifdef COMP_TRADITIONAL
416
417	pCompProc t_pComp0;
418
419	static int t_comp_c_1(poly p1, poly p2);
420	static int t_comp_1_c(poly p1, poly p2);
421	static int t_comp_lex_c_i(poly p1, poly p2);
422	static int t_comp_c_lex_i(poly p1, poly p2);
423	static int t_comp_revlex_c_i(poly p1, poly p2);
424	static int t_comp_c_revlex_i(poly p1, poly p2);
425
426	inline long LexComp(poly p1, poly p2)
427	{
428	long d, i;
429	for (i=1; i<pVariables; i++)
430	{
431	d = pGetExpDiff(p1, p2, i);
432	if (d) return d;
433	}
434	return pGetExpDiff(p1, p2, pVariables);
435	}
436
437	inline long RevLexComp(poly p1, poly p2)
438	{
439	long d, i;
440	for (i=pVariables; i>1; i--)
441	{
442	d = pGetExpDiff(p1, p2, i);
443	if (d) return d;
444	}
445	return pGetExpDiff(p1, p2, 1);
446	}
447
448	static int t_comp_c_1(poly p1, poly p2)
449	{
450	NonZeroR(pGetComp(p1) - pGetComp(p2),
451	return -pComponentOrder, return pComponentOrder);
452	NonZeroR(pGetOrder(p1) - pGetOrder(p2), return pOrdSgn, return -pOrdSgn);
453	return 0;
454	}
455
456	static int t_comp_1_c(poly p1, poly p2)
457	{
458	NonZeroR(pGetOrder(p1) - pGetOrder(p2), return pOrdSgn, return -pOrdSgn);
459	NonZeroR(pGetComp(p1) - pGetComp(p2),
460	return -pComponentOrder, return pComponentOrder);
461	return 0;
462	}
463
464	static int t_comp_lex_c_i(poly p1, poly p2)
465	{
466	NonZeroR(pGetComp(p1) - pGetComp(p2),
467	return -pComponentOrder, return pComponentOrder);
468	NonZeroR(pGetOrder(p1) - pGetOrder(p2), return pOrdSgn, return -pOrdSgn);
469	NonZeroR(LexComp(p1, p2), return pLexSgn, return -pLexSgn);
470	return 0;
471	}
472
473	static int t_comp_lex_i_c(poly p1, poly p2)
474	{
475	NonZeroR(pGetOrder(p1) - pGetOrder(p2), return pOrdSgn, return -pOrdSgn);
476	NonZeroR(LexComp(p1, p2), return pLexSgn, return -pLexSgn);
477	NonZeroR(pGetComp(p1) - pGetComp(p2),
478	return -pComponentOrder, return pComponentOrder);
479	return 0;
480	}
481
482	static int t_comp_revlex_c_i(poly p1, poly p2)
483	{
484	NonZeroR(pGetComp(p1) - pGetComp(p2),
485	return -pComponentOrder, return pComponentOrder);
486	NonZeroR(pGetOrder(p1) - pGetOrder(p2), return pOrdSgn, return -pOrdSgn);
487	NonZeroR(RevLexComp(p1, p2), return pLexSgn, return -pLexSgn);
488	return 0;
489	}
490
491	static int t_comp_revlex_i_c(poly p1, poly p2)
492	{
493	NonZeroR(pGetOrder(p1) - pGetOrder(p2), return pOrdSgn, return -pOrdSgn);
494	NonZeroR(RevLexComp(p1, p2), return pLexSgn, return -pLexSgn);
495	NonZeroR(pGetComp(p1) - pGetComp(p2),
496	return -pComponentOrder, return pComponentOrder);
497	return 0;
498	}
499
500	#endif // COMP_TRADITIONAL
501
502	/*2
503	* compare the head monomial of p1 and p2 with weight vector
504	*/
505	static int comp1a ( poly p1, poly p2, int f, int l, short * w )
506	{
507	int d= p1->Order - p2->Order;
508	if ( d > 0 /p1->Order > p2->Order/ )
509	return 1;
510	else if ( d < 0 /p1->Order < p2->Order/ )
511	return -1;
512	return 0;
513	}
514
515
516	/---------------------------------------------------/
517
518	/* These functions could be made faster if you use pointers to the
519	* exponent vectors and pointer arithmetic instead of using the
520	* macro pGetExp !!!
521	*/
522
523	/*2
524	* compare the head monomial of p1 and p2 with lexicographic ordering
525	*/
526	static int comp_lp ( poly p1, poly p2, int f, int l, short * w )
527	{
528	int i = f;
529	while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
530	i++;
531	if ( i > l )
532	return 0;
533	if ( pGetExp(p1,i) > pGetExp(p2,i) )
534	return 1;
535	return -1;
536	}
537
538	/*2
539	* compare the head monomial of p1 and p2 with degree reverse lexicographic
540	* ordering
541	*/
542	static int comp_dp ( poly p1, poly p2, int f, int l, short * w )
543	{
544	int i, s1 = 0, s2 = 0;
545
546	for ( i = f; i <= l; i++ )
547	{
548	s1 += pGetExp(p1,i);
549	s2 += pGetExp(p2,i);
550	}
551	if ( s1 == s2 )
552	{
553	i = l;
554	while ( (i >= f ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
555	i--;
556	if ( i < f )
557	return 0;
558	if ( pGetExp(p1,i) > pGetExp(p2,i) )
559	return -1;
560	return 1;
561	}
562	if ( s1 > s2 )
563	return 1;
564	return -1;
565	}
566
567	/*2
568	* compare the head monomial of p1 and p2 with degree lexicographic ordering
569	*/
570	static int comp_Dp ( poly p1, poly p2, int f, int l, short * w )
571	{
572	int i, s1 = 0, s2 = 0;
573
574	for ( i = f; i <= l; i++ )
575	{
576	s1 += pGetExp(p1,i);
577	s2 += pGetExp(p2,i);
578	}
579	if ( s1 == s2 )
580	{
581	i = f;
582	while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
583	i++;
584	if ( i > l )
585	return 0;
586	if ( pGetExp(p1,i) > pGetExp(p2,i) )
587	return 1;
588	return -1;
589	}
590	if ( s1 > s2 )
591	return 1;
592	return -1;
593	}
594
595	/*2
596	* compare the head monomial of p1 and p2 with weighted degree reverse
597	* lexicographic ordering
598	*/
599	static int comp_wp ( poly p1, poly p2, int f, int l, short * w )
600	{
601	int i, s1 = 0, s2 = 0;
602
603	for ( i = f; i <= l; i++, w++ )
604	{
605	s1 += (int)pGetExp(p1,i)(w);
606	s2 += (int)pGetExp(p2,i)(w);
607	}
608	if ( s1 == s2 )
609	{
610	i = l;
611	while ( ( i >= f ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
612	i--;
613	if ( i < f )
614	return 0;
615	if ( pGetExp(p1,i) > pGetExp(p2,i) )
616	return -1;
617	return 1;
618	}
619	if ( s1 > s2 )
620	return 1;
621	return -1;
622	}
623
624	/*2
625	* compare the head monomial of p1 and p2 with weighted degree lexicographic
626	* ordering
627	*/
628	static int comp_Wp ( poly p1, poly p2, int f, int l, short * w )
629	{
630	int i, s1 = 0, s2 = 0;
631
632	for ( i = f; i <= l; i++, w++ )
633	{
634	s1 += (int)pGetExp(p1,i)(w);
635	s2 += (int)pGetExp(p2,i)(w);
636	}
637	if ( s1 == s2 )
638	{
639	i = f;
640	while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
641	i++;
642	if ( i > l )
643	return 0;
644	if ( pGetExp(p1,i) > pGetExp(p2,i) )
645	return 1;
646	return -1;
647	}
648	if ( s1 > s2 )
649	return 1;
650	return -1;
651	}
652
653	/*2
654	* compare the head monomial of p1 and p2 with lexicographic ordering
655	* (power series case)
656	*/
657	static int comp_ls ( poly p1, poly p2, int f, int l, short * w )
658	{
659	int i;
660
661	i = f;
662	while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
663	i++;
664	if ( i > l )
665	return 0;
666	if ( pGetExp(p1,i) < pGetExp(p2,i) )
667	return 1;
668	return -1;
669	}
670
671	/*2
672	* compare the head monomial of p1 and p2 with degree reverse lexicographic
673	* ordering (power series case)
674	*/
675	static int comp_ds ( poly p1, poly p2, int f, int l, short * w )
676	{
677	int i, s1 = 0, s2 = 0;
678
679	for ( i = f; i <= l; i++ )
680	{
681	s1 += pGetExp(p1,i);
682	s2 += pGetExp(p2,i);
683	}
684	if ( s1 == s2 )
685	{
686	i = l;
687	while ( ( i >= f ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
688	i--;
689	if ( i < f )
690	return 0;
691	if ( pGetExp(p1,i) < pGetExp(p2,i) )
692	return 1;
693	return -1;
694	}
695	if ( s1 < s2 )
696	return 1;
697	return -1;
698	}
699
700	/*2
701	* compare the head monomial of p1 and p2 with degree lexicographic ordering
702	* (power series case)
703	*/
704	static int comp_Ds ( poly p1, poly p2, int f, int l, short * w )
705	{
706	int i, s1 = 0, s2 = 0;
707
708	for ( i = f; i <= l; i++ )
709	{
710	s1 += pGetExp(p1,i);
711	s2 += pGetExp(p2,i);
712	}
713	if ( s1 == s2 )
714	{
715	i = f;
716	while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
717	i++;
718	if ( i > l )
719	return 0;
720	if ( pGetExp(p1,i) < pGetExp(p2,i) )
721	return -1;
722	return 1;
723	}
724	if ( s1 < s2 )
725	return 1;
726	return -1;
727	}
728
729	/*2
730	* compare the head monomial of p1 and p2 with weighted degree reverse
731	* lexicographic ordering (power series case)
732	*/
733	static int comp_ws ( poly p1, poly p2, int f, int l, short * w )
734	{
735	int i, s1 = 0, s2 = 0;
736
737	for ( i = f; i <= l; i++, w++ )
738	{
739	s1 += (int)pGetExp(p1,i)(w);
740	s2 += (int)pGetExp(p2,i)(w);
741	}
742	if ( s1 == s2 )
743	{
744	i = l;
745	while ( ( i >= f ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
746	i--;
747	if ( i < f )
748	return 0;
749	if ( pGetExp(p1,i) < pGetExp(p2,i) )
750	return 1;
751	return -1;
752	}
753	if ( s1 < s2 )
754	return 1;
755	return -1;
756	}
757
758	/*2
759	* compare the head monomial of p1 and p2 with weighted degree lexicographic
760	* ordering (power series case)
761	*/
762	static int comp_Ws ( poly p1, poly p2, int f, int l, short * w )
763	{
764	int i, s1 = 0, s2 = 0;
765
766	for ( i = f; i <= l; i++, w++ )
767	{
768	s1 += (int)pGetExp(p1,i)(w);
769	s2 += (int)pGetExp(p2,i)(w);
770	}
771	if ( s1 == s2 )
772	{
773	i = f;
774	while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
775	i++;
776	if ( i > l )
777	return 0;
778	if ( pGetExp(p1,i) < pGetExp(p2,i) )
779	return -1;
780	return 1;
781	}
782	if ( s1 < s2 )
783	return 1;
784	return -1;
785	}
786
787	/*2
788	* compare the head monomial of p1 and p2 with matrix order
789	* w contains a series of l-f+1 lines
790	*/
791	static int comp_M ( poly p1, poly p2, int f, int l, short * w )
792	{
793	int i, j, s1, s2;
794
795	for ( i = f; i <= l; i++ )
796	{
797	s1 = s2 = 0;
798	for ( j = f; j <= l; j++, w++ )
799	{
800	s1 += (int)pGetExp(p1,j)(int)(w);
801	s2 += (int)pGetExp(p2,j)(int)(w);
802	}
803	if ( s1 < s2 )
804	return -1;
805	if ( s1 > s2 )
806	return 1;
807	/* now w points to the last element of the current row, the next w++ */
808	/* moves on to the first element of the next row ! */
809	}
810	return 0;
811	}
812
813	/*2
814	* compare the head monomial of p1 and p2 with weight vector
815	*/
816	static int comp_a ( poly p1, poly p2, int f, int l, short * w )
817	{
818	int i, s1 = 0, s2 = 0;
819
820	for ( i = f; i <= l; i++, w++ )
821	{
822	s1 += (int)pGetExp(p1,i)(w);
823	s2 += (int)pGetExp(p2,i)(w);
824	}
825	if ( s1 > s2 )
826	return 1;
827	if ( s1 < s2 )
828	return -1;
829	return 0;
830	}
831
832	/*2
833	* compare the head monomial of p1 and p2 with module component
834	*/
835	static int comp_c ( poly p1, poly p2, int f, int l, short * w )
836	{
837	if ( pGetComp(p1) > pGetComp(p2) )
838	return -pComponentOrder;
839	if ( pGetComp(p1) < pGetComp(p2) )
840	return pComponentOrder;
841	return 0;
842	}
843
844	/---------------------------------------------------------------/
845
846	/*2
847	* compare p1 and p2 by a block ordering
848	* uses (*bcomph[])() to do the real work
849	*/
850	static int BlockComp(poly p1, poly p2)
851	{
852	int res, i, e, a;
853
854	/4 compare in all blocks,
855	* each block has var numbers a(=block0[i]) to e (=block1[i])*
856	* the block number starts with 0*/
857	e = 0;
858	i = 0;
859	loop
860	{
861	a = block0[i];
862	e = block1[i];
863	res = (*bcomph[i])(p1, p2, a, e , polys_wv[i]);
864	if (res)
865	return res;
866	i++;
867	if (order[i]==0)
868	break;
869	}
870	return 0;
871	}
872
873	int pComp(poly p1, poly p2)
874	{
875	if (p2==NULL)
876	return 1;
877	if (p1==NULL)
878	return -1;
879	return pComp0(p1,p2);
880	}
881
882	/----------pComp handling for syzygies---------------------/
883	static void newHeadsB(polyset actHeads,int length)
884	{
885	int i;
886	int* newOrder=(int)Alloc(lengthsizeof(int));
887
888	for (i=0;i<length;i++)
889	{
890	if (actHeads[i])
891	{
892	newOrder[i] = SchreyerOrd[pGetComp(actHeads[i])-1];
893	}
894	else
895	{
896	newOrder[i]=0;
897	}
898	}
899	Free((ADDRESS)SchreyerOrd,maxSchreyer*sizeof(int));
900	SchreyerOrd = newOrder;
901	maxSchreyer = length;
902	/*
903	*for (i=0;i<maxSchreyer;i++); Print("%d ",SchreyerOrd[i]);
904	*PrintLn();
905	*/
906	}
907
908	int mcompSchrB(poly p1,poly p2)
909	{
910	int CompP1=pGetComp(p1),CompP2=pGetComp(p2),result,
911	cP1=SchreyerOrd[CompP1-1],cP2=SchreyerOrd[CompP2-1];
912
913	if (CompP1==CompP2) return pCompOld(p1,p2);
914	pSetComp(p1,cP1);
915	pSetComp(p2,cP2);
916	result = pCompOld(p1,p2);
917	pSetComp(p1,CompP1);
918	pSetComp(p2,CompP2);
919	if (!result)
920	{
921	if (CompP1>CompP2)
922	return -1;
923	else if (CompP1<CompP2)
924	return 1;
925	}
926	return result;
927	}
928
929
930	static void newHeadsM(polyset actHeads,int length)
931	{
932	int i;
933	int* newOrder=
934	(int)Alloc((length+maxSchreyer-indexShift)sizeof(int));
935
936	for (i=0;i<length+maxSchreyer-indexShift;i++)
937	newOrder[i]=0;
938	for (i=indexShift;i<maxSchreyer;i++)
939	{
940	newOrder[i-indexShift] = SchreyerOrd[i];
941	SchreyerOrd[i] = 0;
942	}
943	for (i=maxSchreyer-indexShift;i<length+maxSchreyer-indexShift;i++)
944	newOrder[i] = newOrder[pGetComp(actHeads[i-maxSchreyer+indexShift])-1];
945	Free((ADDRESS)SchreyerOrd,maxSchreyer*sizeof(int));
946	SchreyerOrd = newOrder;
947	indexShift = maxSchreyer-indexShift;
948	maxSchreyer = length+indexShift;
949	}
950
951	/*2
952	* compute the length of a polynomial (in l)
953	* and the degree of the monomial with maximal degree:
954	* this is NOT the last one and the module component
955	* has to be <= indexShift
956	*/
957	static int ldegSchrM(poly p,int *l)
958	{
959	int t,max;
960
961	(*l)=1;
962	max=pFDeg(p);
963	while ((pNext(p)!=NULL) && (pGetComp(pNext(p))<=indexShift))
964	{
965	pIter(p);
966	t=pFDeg(p);
967	if (t>max) max=t;
968	(*l)++;
969	}
970	return max;
971	}
972
973	int mcompSchrM(poly p1,poly p2)
974	{
975	if ( pGetComp(p1)<=indexShift)
976	{
977	if ( pGetComp(p2)>indexShift) return 1;
978	}
979	else if ( pGetComp(p2)<=indexShift) return -1;
980	return mcompSchrB(p1,p2);
981	}
982
983	void pSetSchreyerOrdM(polyset nextOrder, int length,int comps)
984	{
985	int i;
986
987	if (length!=0)
988	{
989	if (maxSchreyer!=0)
990	newHeadsM(nextOrder, length);
991	else
992	{
993	indexShift = comps;
994	if (!indexShift) indexShift = 1;
995	SchreyerOrd = (int)Alloc((indexShift+length)sizeof(int));
996	maxSchreyer = length+indexShift;
997	for (i=0;i<indexShift;i++)
998	SchreyerOrd[i] = i;
999	for (i=indexShift;i<maxSchreyer;i++)
1000	SchreyerOrd[i] = pGetComp(nextOrder[i-indexShift]);
1001	#ifdef COMP_DEBUG
1002	pCompOld = t_pComp0;
1003	#else
1004	pCompOld = pComp0;
1005	#endif
1006	pComp0 = mcompSchrM;
1007	t_pComp0 = mcompSchrM;
1008	pLDegOld = pLDeg;
1009	pLDeg = ldegSchrM;
1010	}
1011	}
1012	else
1013	{
1014	if (maxSchreyer!=0)
1015	{
1016	Free((ADDRESS)SchreyerOrd,maxSchreyer*sizeof(int));
1017	maxSchreyer = 0;
1018	indexShift = 0;
1019	#ifdef COMP_DEBUG
1020	pComp0 = debug_comp;
1021	#else
1022	pComp0 = pCompOld;
1023	#endif
1024	t_pComp0 = pCompOld;
1025	pLDeg = pLDegOld;
1026	}
1027	}
1028	}
1029
1030	/*2
1031	*the pComp0 for normal syzygies
1032	*compares monomials in the usual ring order (pCompOld)
1033	*but groups module indecees according indexBounds befor
1034	*/
1035	static int mcompSyz(poly p1,poly p2)
1036	{
1037	if (pGetComp(p1)<=maxBound)
1038	{
1039	if (pGetComp(p2)>maxBound) return 1;
1040	}
1041	else if (pGetComp(p2)<=maxBound)
1042	{
1043	return -1;
1044	}
1045	return pCompOld(p1,p2);
1046	}
1047
1048	void pSetSyzComp(int k)
1049	{
1050	if (k!=0)
1051	{
1052	if (maxBound==0)
1053	{
1054	#ifdef COMP_DEBUG
1055	pCompOld = t_pComp0;
1056	pComp0 = mcompSyz;
1057	t_pComp0 = mcompSyz;
1058	#else
1059	pCompOld = pComp0;
1060	pComp0 = mcompSyz;
1061	t_pComp0 = mcompSyz;
1062	#endif
1063	}
1064	maxBound = k;
1065	}
1066	else
1067	{
1068	if (maxBound!=0)
1069	{
1070	#ifdef COMP_DEBUG
1071	pComp0 = debug_comp;
1072	#else
1073	pComp0 = pCompOld;
1074	#endif
1075	t_pComp0 = pCompOld;
1076	maxBound = 0;
1077	}
1078	}
1079	}
1080
1081	/*2
1082	* the type of the module ordering: C: -1, c: 1
1083	*/
1084	int pModuleOrder()
1085	{
1086	return pComponentOrder;
1087	}
1088
1089	/* -------------------------------------------------------------------*/
1090	/* several possibilities for pFDeg: the degree of the head term */
1091	/*2
1092	* compute the degree of the leading monomial of p
1093	* the ordering is compatible with degree, use a->order
1094	*/
1095	int pDeg(poly a)
1096	{
1097	return ((a!=NULL) ? (a->Order) : (-1));
1098	}
1099
1100	/*2
1101	* compute the degree of the leading monomial of p
1102	* with respect to weigths 1
1103	* (all are 1 so save multiplications or they are of different signs)
1104	* the ordering is not compatible with degree so do not use p->Order
1105	*/
1106	int pTotaldegree(poly p)
1107	{
1108	return pExpQuerSum(p);
1109	}
1110
1111	/*2
1112	* compute the degree of the leading monomial of p
1113	* with respect to weigths from the ordering
1114	* the ordering is not compatible with degree so do not use p->Order
1115	*/
1116	int pWTotaldegree(poly p)
1117	{
1118	int i, k;
1119	int j =0;
1120
1121	// iterate through each block:
1122	for (i=0;order[i]!=0;i++)
1123	{
1124	switch(order[i])
1125	{
1126	case ringorder_wp:
1127	case ringorder_ws:
1128	case ringorder_Wp:
1129	case ringorder_Ws:
1130	for (k=block0[i];k<=block1[i];k++)
1131	{ // in jedem block:
1132	j+= pGetExp(p,k)*polys_wv[i][k-block0[i]];
1133	}
1134	break;
1135	case ringorder_M:
1136	case ringorder_lp:
1137	case ringorder_dp:
1138	case ringorder_ds:
1139	case ringorder_Dp:
1140	case ringorder_Ds:
1141	for (k=block0[i];k<=block1[i];k++)
1142	{
1143	j+= pGetExp(p,k);
1144	}
1145	break;
1146	case ringorder_c:
1147	case ringorder_C:
1148	break;
1149	case ringorder_a:
1150	for (k=block0[i];k<=block1[i];k++)
1151	{ // only one line
1152	j+= pGetExp(p,k)*polys_wv[i][k-block0[i]];
1153	}
1154	return j;
1155	}
1156	}
1157	return j;
1158	}
1159	/* ---------------------------------------------------------------------*/
1160	/* several possibilities for pLDeg: the maximal degree of a monomial in p*/
1161	/* compute in l also the pLength of p */
1162
1163	/*2
1164	* compute the length of a polynomial (in l)
1165	* and the degree of the monomial with maximal degree: the last one
1166	*/
1167	static int ldeg0(poly p,int *l)
1168	{
1169	Exponent_t k= pGetComp(p);
1170	int ll=1;
1171
1172	while ((pNext(p)!=NULL) && (pGetComp(pNext(p))==k))
1173	{
1174	pIter(p);
1175	ll++;
1176	}
1177	*l=ll;
1178	return (p->Order);
1179	}
1180
1181	/*2
1182	* compute the length of a polynomial (in l)
1183	* and the degree of the monomial with maximal degree: the last one
1184	* but search in all components before syzcomp
1185	*/
1186	static int ldeg0c(poly p,int *l)
1187	{
1188	int o=pFDeg(p);
1189	int ll=1;
1190
1191	if (maxBound/syzComp/==0)
1192	{
1193	while ((p=pNext(p))!=NULL)
1194	{
1195	o=pFDeg(p);
1196	ll++;
1197	}
1198	}
1199	else
1200	{
1201	while ((p=pNext(p))!=NULL)
1202	{
1203	if (pGetComp(p)<=maxBound/syzComp/)
1204	{
1205	o=pFDeg(p);
1206	ll++;
1207	}
1208	else break;
1209	}
1210	}
1211	*l=ll;
1212	return o;
1213	}
1214
1215	/*2
1216	* compute the length of a polynomial (in l)
1217	* and the degree of the monomial with maximal degree: the first one
1218	* this works for the polynomial case with degree orderings
1219	* (both c,dp and dp,c)
1220	*/
1221	static int ldegb(poly p,int *l)
1222	{
1223	Exponent_t k= pGetComp(p);
1224	int o = p->Order;
1225	int ll=1;
1226
1227	while (((p=pNext(p))!=NULL) && (pGetComp(p)==k))
1228	{
1229	ll++;
1230	}
1231	*l=ll;
1232	return o;
1233	}
1234
1235	/*2
1236	* compute the length of a polynomial (in l)
1237	* and the degree of the monomial with maximal degree:
1238	* this is NOT the last one, we have to look for it
1239	*/
1240	static int ldeg1(poly p,int *l)
1241	{
1242	Exponent_t k= pGetComp(p);
1243	int ll=1;
1244	int t,max;
1245
1246	max=pFDeg(p);
1247	while (((p=pNext(p))!=NULL) && (pGetComp(p)==k))
1248	{
1249	t=pFDeg(p);
1250	if (t>max) max=t;
1251	ll++;
1252	}
1253	*l=ll;
1254	return max;
1255	}
1256
1257	/*2
1258	* compute the length of a polynomial (in l)
1259	* and the degree of the monomial with maximal degree:
1260	* this is NOT the last one, we have to look for it
1261	* in all components
1262	*/
1263	static int ldeg1c(poly p,int *l)
1264	{
1265	int ll=1;
1266	int t,max;
1267
1268	max=pFDeg(p);
1269	while ((p=pNext(p))!=NULL)
1270	{
1271	if ((maxBound/syzComp/==0) \|\| (pGetComp(p)<=maxBound/syzComp/))
1272	{
1273	if ((t=pFDeg(p))>max) max=t;
1274	ll++;
1275	}
1276	else break;
1277	}
1278	*l=ll;
1279	return max;
1280	}
1281
1282	/* -------------------------------------------------------- */
1283	/* set the variables for a choosen ordering */
1284
1285
1286	/*2
1287	* sets the comparision routine for monomials: for the first block
1288	* of variables (or is the number of the ordering)
1289	*/
1290	#ifdef COMP_FAST
1291	static void SimpleChoose(int or, int comp_order, pCompProc *p)
1292	#else
1293	static void SimpleChoose(int or, pCompProc *p)
1294	#endif
1295	{
1296	#ifdef COMP_TRADITIONAL
1297	if (pVariables <= 1)
1298	{
1299	t_pComp0 = t_comp_1_c;
1300	}
1301	else
1302	{
1303	switch(or)
1304	{
1305	case ringorder_lp:
1306	case ringorder_Dp:
1307	case ringorder_Wp:
1308	case ringorder_Ds:
1309	case ringorder_Ws:
1310	case ringorder_ls:
1311	t_pComp0 = t_comp_lex_i_c;
1312	pLexSgn = 1;
1313	break;
1314
1315	#ifdef PDEBUG
1316	case ringorder_unspec:
1317	case ringorder_dp:
1318	case ringorder_wp:
1319	case ringorder_ds:
1320	case ringorder_ws:
1321	#else
1322	default:
1323	#endif
1324	t_pComp0 = t_comp_revlex_i_c;
1325	pLexSgn = -1;
1326	break;
1327	#ifdef PDEBUG
1328	default:
1329	Werror("wrong internal ordering:%d at %s, l:%d\n",or,__FILE__,__LINE__);
1330	#endif
1331	}
1332	}
1333	if (or == ringorder_lp \|\| or == ringorder_ls)
1334	{
1335	pLexOrder=TRUE;
1336	pFDeg = pTotaldegree;
1337	pLDeg = ldeg1c;
1338	if (or == ringorder_ls) pLexSgn = -1;
1339	}
1340
1341	*p = t_pComp0;
1342	#endif
1343	#ifdef COMP_FAST
1344	switch(or)
1345	{
1346	case ringorder_dp:
1347	case ringorder_wp:
1348	case ringorder_ds:
1349	case ringorder_ws:
1350	case ringorder_ls:
1351	case ringorder_unspec:
1352	pSetVarIndicies_RevLex(pVariables);
1353	pLexSgn = -1;
1354	if (comp_order == ringorder_C \|\| or == ringorder_unspec)
1355	{
1356	if (pVariables1W == 1)
1357	f_pComp0 = f_comp_1_c;
1358	else if (pVariables1W == 2)
1359	f_pComp0 = f_comp_2_c;
1360	else if (pVariables1W & 1)
1361	f_pComp0 = f_comp_2i_1_c;
1362	else
1363	f_pComp0 = f_comp_2i_c;
1364	}
1365	else
1366	{
1367	if (pVariables1W == 1)
1368	f_pComp0 = f_comp_1;
1369	else if (pVariables1W == 2)
1370	f_pComp0 = f_comp_2;
1371	else if (pVariables1W & 1)
1372	f_pComp0 = f_comp_2i_1;
1373	else
1374	f_pComp0 = f_comp_2i;
1375	}
1376	break;
1377
1378	#ifdef PDEBUG
1379	case ringorder_lp:
1380	case ringorder_Dp:
1381	case ringorder_Wp:
1382	case ringorder_Ds:
1383	case ringorder_Ws:
1384	#else
1385	default:
1386	#endif
1387	pSetVarIndicies_Lex(pVariables);
1388	pLexSgn = 1;
1389	if (comp_order == ringorder_c)
1390	{
1391	if (pVariables1W == 1)
1392	f_pComp0 = f_comp_1_c;
1393	else if (pVariables1W == 2)
1394	f_pComp0 = f_comp_2_c;
1395	else if (pVariables1W & 1)
1396	f_pComp0 = f_comp_2i_1_c;
1397	else
1398	f_pComp0 = f_comp_2i_c;
1399	}
1400	else
1401	{
1402	if (pVariables1W == 1)
1403	f_pComp0 = f_comp_1;
1404	else if (pVariables1W == 2)
1405	f_pComp0 = f_comp_2;
1406	else if (pVariables1W & 1)
1407	f_pComp0 = f_comp_2i_1;
1408	else
1409	f_pComp0 = f_comp_2i;
1410	}
1411	#ifdef PDEBUG
1412	break;
1413	default:
1414	Werror("wrong internal ordering:%d at %s, l:%d\n",or,__FILE__,__LINE__);
1415	#endif
1416	}
1417
1418	if (or == ringorder_lp \|\| or == ringorder_ls)
1419	{
1420	pLexOrder=TRUE;
1421	pFDeg = pTotaldegree;
1422	pLDeg = ldeg1c;
1423	if (or == ringorder_ls)
1424	pSetVarIndicies_Lex(pVariables);
1425	}
1426	*p = f_pComp0;
1427	#endif // COMP_FAST
1428
1429	#ifdef COMP_DEBUG
1430	*p = debug_comp;
1431	#endif
1432	}
1433
1434	/*2
1435	* sets pSetm
1436	* (according or = order of first block)
1437	*/
1438	static void SetpSetm(int or, int ip)
1439	{
1440	switch(or)
1441	{
1442	case ringorder_lp:
1443	case ringorder_ls:
1444	if (pVariables>1)
1445	pSetm= setlex2;
1446	else
1447	pSetm= setlex1;
1448	break;
1449	case ringorder_dp:
1450	case ringorder_Dp:
1451	case ringorder_ds:
1452	case ringorder_Ds:
1453	case ringorder_unspec:
1454	pSetm= setdeg1;
1455	break;
1456	case ringorder_a:
1457	case ringorder_wp:
1458	case ringorder_Wp:
1459	case ringorder_ws:
1460	case ringorder_Ws:
1461	case ringorder_M:
1462	pSetm= setdeg1w;
1463	firstwv=polys_wv[ip];
1464	break;
1465	case ringorder_c:
1466	case ringorder_C:
1467	return;
1468	/do not set firstBlockEnds for this orderings/
1469	#ifdef TEST
1470	default:
1471	Werror("wrong internal ordering:%d at %s, l:%d\n",or,__FILE__,__LINE__);
1472	#endif
1473	}
1474	firstBlockEnds=block1[ip];
1475	}
1476
1477	/*2
1478	* sets the comparision routine for monomials: for the first block
1479	* of variables (or is the number of the ordering)
1480	*/
1481	#ifdef COMP_FAST
1482	static void SimpleChooseC(int or, int comp_order, pCompProc *p)
1483	#else
1484	static void SimpleChooseC(int or, pCompProc *p)
1485	#endif
1486	{
1487	#ifdef COMP_TRADITIONAL
1488	if (pVariables <= 1)
1489	{
1490	t_pComp0 = t_comp_c_1;
1491	}
1492	else
1493	{
1494	switch(or)
1495	{
1496	case ringorder_lp:
1497	case ringorder_Dp:
1498	case ringorder_Wp:
1499	case ringorder_Ds:
1500	case ringorder_Ws:
1501	case ringorder_ls:
1502	t_pComp0 = t_comp_lex_c_i;
1503	pLexSgn = 1;
1504	break;
1505
1506	#ifdef PDEBUG
1507	case ringorder_unspec:
1508	case ringorder_dp:
1509	case ringorder_wp:
1510	case ringorder_ds:
1511	case ringorder_ws:
1512	#else
1513	default:
1514	#endif
1515	t_pComp0 = t_comp_revlex_c_i;
1516	pLexSgn = -1;
1517	break;
1518	#ifdef PDEBUG
1519	default:
1520	Werror("wrong internal ordering:%d at %s, l:%d\n",or,__FILE__,__LINE__);
1521	#endif
1522	}
1523	}
1524	if (or == ringorder_lp \|\| or == ringorder_ls)
1525	{
1526	pLexOrder=TRUE;
1527	pFDeg = pTotaldegree;
1528	pLDeg = ldeg1c;
1529	if (or == ringorder_ls) pLexSgn = -1;
1530	}
1531
1532	*p = t_pComp0;
1533	#endif
1534	#ifdef COMP_FAST
1535	switch(or)
1536	{
1537	case ringorder_dp:
1538	case ringorder_wp:
1539	case ringorder_ds:
1540	case ringorder_ls:
1541	case ringorder_ws:
1542	pSetVarIndicies_RevLex(pVariables);
1543	pLexSgn = -1;
1544	if (pVariablesW == 1)
1545	f_pComp0 = f_comp_c_1;
1546	else if (pVariablesW == 2)
1547	f_pComp0 = f_comp_c_2;
1548	else if (pVariablesW & 1)
1549	f_pComp0 = f_comp_c_2i_1;
1550	else
1551	f_pComp0 = f_comp_c_2i;
1552	break;
1553
1554	#ifdef PDEBUG
1555	case ringorder_lp:
1556	case ringorder_Dp:
1557	case ringorder_Wp:
1558	case ringorder_Ds:
1559	case ringorder_Ws:
1560	#else
1561	default:
1562	#endif
1563	pSetVarIndicies_Lex(pVariables);
1564	pLexSgn = 1;
1565	if (pVariablesW == 1)
1566	f_pComp0 = f_comp_c_1;
1567	else if (pVariablesW == 2)
1568	f_pComp0 = f_comp_c_2;
1569	else if (pVariablesW & 1)
1570	f_pComp0 = f_comp_c_2i_1;
1571	else
1572	f_pComp0 = f_comp_c_2i;
1573	#ifdef PDEBUG
1574	break;
1575	default:
1576	Werror("wrong internal ordering:%d at %s, l:%d\n",or,__FILE__,__LINE__);
1577	#endif
1578	}
1579	if (or == ringorder_lp \|\| or == ringorder_ls)
1580	{
1581	pLexOrder=TRUE;
1582	pFDeg = pTotaldegree;
1583	pLDeg = ldeg1c;
1584	if (or == ringorder_ls)
1585	pSetVarIndicies_Lex(pVariables);
1586	}
1587	*p = f_pComp0;
1588	#endif // COMP_FAST
1589
1590	#ifdef COMP_DEBUG
1591	*p = debug_comp;
1592	#endif
1593	}
1594
1595	/*2
1596	* sets the comparision routine for monomials: for all but the first
1597	* block of variables (ip is the block number, or the number of the ordering)
1598	*/
1599	static void HighSet(int ip, int or)
1600	{
1601	switch(or)
1602	{
1603	case ringorder_lp:
1604	bcomph[ip]=comp_lp;
1605	if (pOrdSgn==-1) pMixedOrder=TRUE;
1606	break;
1607	case ringorder_dp:
1608	bcomph[ip]=comp_dp;
1609	if (pOrdSgn==-1) pMixedOrder=TRUE;
1610	break;
1611	case ringorder_Dp:
1612	bcomph[ip]=comp_Dp;
1613	if (pOrdSgn==-1) pMixedOrder=TRUE;
1614	break;
1615	case ringorder_wp:
1616	bcomph[ip]=comp_wp;
1617	if (pOrdSgn==-1) pMixedOrder=TRUE;
1618	break;
1619	case ringorder_Wp:
1620	bcomph[ip]=comp_Wp;
1621	if (pOrdSgn==-1) pMixedOrder=TRUE;
1622	break;
1623	case ringorder_ls:
1624	bcomph[ip]=comp_ls;
1625	break;
1626	case ringorder_ds:
1627	bcomph[ip]=comp_ds;
1628	break;
1629	case ringorder_Ds:
1630	bcomph[ip]=comp_Ds;
1631	break;
1632	case ringorder_ws:
1633	bcomph[ip]=comp_ws;
1634	break;
1635	case ringorder_Ws:
1636	bcomph[ip]=comp_Ws;
1637	break;
1638	case ringorder_c:
1639	pComponentOrder=1;
1640	bcomph[ip]=comp_c;
1641	break;
1642	case ringorder_C:
1643	pComponentOrder=-1;
1644	bcomph[ip]=comp_c;
1645	break;
1646	case ringorder_M:
1647	bcomph[ip]=comp_M;
1648	pMixedOrder=TRUE;
1649	break;
1650	case ringorder_a:
1651	if (pOrdSgn==-1) pMixedOrder=TRUE;
1652	if (ip==0)
1653	bcomph[0]=comp1a;
1654	else
1655	bcomph[ip]=comp_a;
1656	break;
1657	#ifdef TEST
1658	default:
1659	Werror("wrong internal ordering:%d at %s, l:%d\n",or,__FILE__,__LINE__);
1660	#endif
1661	}
1662	}
1663
1664	#ifdef DIV_COUNT
1665	struct div_triple
1666	{
1667	unsigned long greater;
1668	unsigned long equal;
1669	unsigned long smaller;
1670	};
1671
1672	static div_triple* DivTriples = NULL;
1673	static int NumberOfDivTriples = 0;
1674	static unsigned long DivCountTotal = 0;
1675	static unsigned long DivCountDiv = 0;
1676	static unsigned long LexDivCount = 0;
1677	static unsigned long RevLexDivCount = 0;
1678	static unsigned long LexDivCount2 = 0;
1679	static unsigned long RevLexDivCount2 = 0;
1680	static int pVariables2;
1681	static int pVariablesEven;
1682
1683	void InitDivCount(int nvars)
1684	{
1685	if (nvars & 1) pVariablesEven = 0;
1686	else pVariablesEven = 1;
1687
1688	pVariables2 = (pVariables -1) / 2;
1689
1690	if (DivTriples == NULL)
1691	{
1692	DivTriples = (div_triple) Alloc0(nvarssizeof(div_triple));
1693	NumberOfDivTriples = nvars;
1694	}
1695	else if (nvars > NumberOfDivTriples)
1696	{
1697	DivTriples = (div_triple*) ReAlloc(DivTriples,
1698	NumberOfDivTriples*sizeof(div_triple),
1699	nvars*sizeof(div_triple));
1700	for (int i=NumberOfDivTriples; i<nvars; i++)
1701	{
1702	DivTriples[i].greater = 0;
1703	DivTriples[i].smaller = 0;
1704	DivTriples[i].equal = 0;
1705	}
1706
1707	NumberOfDivTriples = nvars;
1708	}
1709	}
1710
1711	void ResetDivCount()
1712	{
1713	int i;
1714	for (i=0; i<NumberOfDivTriples;i++)
1715	{
1716	DivTriples[i].greater = 0;
1717	DivTriples[i].equal = 0;
1718	DivTriples[i].smaller = 0;
1719	}
1720
1721	DivCountTotal = 0;
1722	DivCountDiv = 0;
1723	LexDivCount = 0;
1724	RevLexDivCount = 0;
1725	LexDivCount2 = 0;
1726	RevLexDivCount2 = 0;
1727	}
1728
1729	void OutputDivCount()
1730	{
1731	printf("Total : %10u\n", DivCountTotal);
1732	if (DivCountTotal == 0) DivCountTotal++;
1733	printf("Div : %10u \t %.4f\n", DivCountDiv,
1734	(float) DivCountDiv / (float) DivCountTotal);
1735	printf("LexCount : %10u \t %.4f\n", LexDivCount,
1736	(float) LexDivCount / (float) DivCountTotal);
1737	printf("RevLexCount: %10u \t %.4f\n", RevLexDivCount,
1738	(float) RevLexDivCount / (float) DivCountTotal);
1739	printf("LexCount2 : %10u \t %.4f\n", LexDivCount2,
1740	(float) LexDivCount2 / (float) DivCountTotal);
1741	printf("RevLexCount2: %10u \t %.4f\n", RevLexDivCount2,
1742	(float) RevLexDivCount2 / (float) DivCountTotal);
1743
1744	for (int i=0; i<NumberOfDivTriples; i++)
1745	{
1746	unsigned long total = DivTriples[i].greater +
1747	DivTriples[i].equal + DivTriples[i].smaller;
1748	// avoid total==0
1749	if (total == 0) total=1;
1750	printf(" [%d]: %10u (%.4f) %10u (%.4f) %10u (%.4f)\n", i,
1751	DivTriples[i].greater,
1752	(float) DivTriples[i].greater / (float) total ,
1753	DivTriples[i].equal,
1754	(float) DivTriples[i].equal / (float) total ,
1755	DivTriples[i].smaller,
1756	(float) DivTriples[i].smaller / (float) total );
1757	}
1758	}
1759
1760	BOOLEAN pDivisibleBy(poly a, poly b)
1761	{
1762	if ((a!=NULL)&&(( pGetComp(a)==0) \|\| ( pGetComp(a) == pGetComp(b))))
1763	{
1764	int i;
1765	Exponent_t *e1=&( pGetExp(a,1));
1766	Exponent_t *e2=&( pGetExp(b,1));
1767	BOOLEAN res = TRUE, res2 = TRUE;
1768	DivCountTotal++;
1769
1770	for (i=0; i<pVariables; i++)
1771	{
1772	if (res == TRUE) LexDivCount++;
1773	if (e1 > e2)
1774	{
1775	DivTriples[i].greater++;
1776	res = FALSE;
1777	}
1778	else if (e1 == e2) DivTriples[i].equal++;
1779	else DivTriples[i].smaller++;
1780	e1++;
1781	e2++;
1782	}
1783	if (res == TRUE) DivCountDiv++;
1784	e1 = &(pGetExp(a,pVariables));
1785	e2 = &(pGetExp(b,pVariables));
1786
1787	for (i=0; i<pVariables; i++)
1788	{
1789	RevLexDivCount++;
1790	if (e1 > e2)
1791	{
1792	res2 = FALSE;
1793	break;
1794	}
1795	e1--;
1796	e2--;
1797	}
1798	if (res != res2)
1799	fprintf(stderr, "Error in RevLexCount\n");
1800
1801	if (pVariables < 4) return res;
1802	LexTop:
1803	res2 = TRUE;
1804	int* s1 = (int*) &( pGetExp(a,2));
1805	int* s2 = (int*) &( pGetExp(b,2));
1806	for (i=0; i<pVariables2; i++, s1++, s2++)
1807	{
1808	LexDivCount2++;
1809	if (s1 > s2)
1810	{
1811	res2 = FALSE;
1812	break;
1813	}
1814	}
1815	if (res2 == TRUE)
1816	{
1817	LexDivCount2++;
1818	if ( pGetExp(a,1) > pGetExp(b,1)) res2 = FALSE;
1819	if (res2 == TRUE && pVariablesEven)
1820	{
1821	LexDivCount2++;
1822	if (pGetExp(a,pVariables) > pGetExp(b,pVariables))
1823	res2 = FALSE;
1824	}
1825	}
1826	if (res2 == TRUE)
1827	{
1828	e1 = &( pGetExp(a,2));
1829	e2 = &( pGetExp(b,2));
1830	for (i=0; i<pVariables2; i++, e1 += 2, e2 += 2)
1831	{
1832	LexDivCount++;
1833	#ifdef WORDS_BIG_ENDIAN
1834	if (e1[1] > e2[1])
1835	#else
1836	if (e1 > e2)
1837	#endif
1838	{
1839	res2 = FALSE;
1840	break;
1841	}
1842	}
1843	}
1844	if (res != res2)
1845	{
1846	fprintf(stderr, "Error in LexDivCount2\n");
1847	// goto LexTop;
1848	}
1849
1850	RevLexTop:
1851	res2 = TRUE;
1852	if (pVariablesEven)
1853	{
1854	s1 = (int*) &(pGetExp(a,pVariables-2));
1855	s2 = (int*) &(pGetExp(b,pVariables-2));
1856	}
1857	else
1858	{
1859	s1 = (int*) &(pGetExp(a,pVariables-1));
1860	s2 = (int*) &(pGetExp(b,pVariables-1));
1861	}
1862	for (i=0; i< pVariables2; i++, s1--, s2--)
1863	{
1864	RevLexDivCount2++;
1865	if (s1 > s2)
1866	{
1867	res2 = FALSE;
1868	break;
1869	}
1870	}
1871	if (res2 == TRUE)
1872	{
1873	if (pVariablesEven)
1874	{
1875	RevLexDivCount2++;
1876	if (pGetExp(a,pVariables) > pGetExp(b,pVariables))
1877	res2 = FALSE;
1878	}
1879	if (res2 == TRUE)
1880	{
1881	RevLexDivCount2++;
1882	if ( pGetExp(a,1) > pGetExp(b,1)) res2 = FALSE;
1883	}
1884	}
1885	if (res2 == TRUE)
1886	{
1887	if (pVariablesEven)
1888	{
1889	e1 = &(pGetExp(a,pVariables-2));
1890	e2 = &(pGetExp(b,pVariables-2));
1891	}
1892	else
1893	{
1894	e1 = &(pGetExp(a,pVariables-1));
1895	e2 = &(pGetExp(b,pVariables-1));
1896	}
1897	for (i=0; i<pVariables2; i++, e1 -= 2, e2 -= 2)
1898	{
1899	#ifdef WORDS_BIG_ENDIAN
1900	if (e1[1] > e2[1])
1901	#else
1902	if (e1 > e2)
1903	#endif
1904	{
1905	res2 = FALSE;
1906	break;
1907	}
1908	}
1909	}
1910	if (res != res2)
1911	{
1912	fprintf(stderr, "Error in RevLexDivCount2\n");
1913	// goto RevLexTop;
1914	}
1915	return res;
1916	}
1917	return FALSE;
1918	}
1919	#endif
1920
1921	/* -------------------------------------------------------- */
1922	/*2
1923	* change all variables to fit the description of the new ring
1924	*/
1925	void pChangeRing(int n, int Sgn, int * orders, int * b0, int * b1,
1926	short ** wv)
1927	{
1928	#ifdef TEST_MAC_ORDER
1929	bNoAdd=FALSE;
1930	#endif
1931	int i;
1932	pComponentOrder=1;
1933	if (ppNoether!=NULL) pDelete(&ppNoether);
1934	#ifdef SRING
1935	pSRING=FALSE;
1936	pAltVars=n+1;
1937	#endif
1938	t_pComp0 = NULL;
1939	#ifdef COMP_FAST
1940	f_pComp0 = NULL;
1941	#endif
1942	pVariables = n;
1943
1944	// set the various size parameters and initialize memory
1945	pMonomSize = POLYSIZE + (pVariables + 1) * sizeof(Exponent_t);
1946	if ((pMonomSize % sizeof(void*)) == 0)
1947	{
1948	pMonomSizeW = pMonomSize/sizeof(void*);
1949	}
1950	else
1951	{
1952	pMonomSizeW = pMonomSize/sizeof(void*) + 1;
1953	pMonomSize = pMonomSizeWsizeof(void);
1954	}
1955
1956	#ifdef COMP_FAST
1957	if ((((pVariables+1)sizeof(Exponent_t)) % sizeof(void)) == 0)
1958	pVariables1W = (pVariables+1)sizeof(Exponent_t) / sizeof(void);
1959	else
1960	pVariables1W = ((pVariables+1)sizeof(Exponent_t) / sizeof(void)) + 1;
1961	if ((((pVariables)sizeof(Exponent_t)) % sizeof(void)) == 0)
1962	pVariablesW = (pVariables)sizeof(Exponent_t) / sizeof(void);
1963	else
1964	pVariablesW = ((pVariables)sizeof(Exponent_t) / sizeof(void)) + 1;
1965
1966	// Set default Var Indicies
1967	pSetVarIndicies(pVariables);
1968	#endif
1969	mmSpecializeBlock(pMonomSize);
1970
1971	pOrdSgn = Sgn;
1972	pVectorOut=(orders[0]==ringorder_c);
1973	order=orders;
1974	block0=b0;
1975	block1=b1;
1976	firstwv=NULL;
1977	polys_wv=wv;
1978	/------- only one real block ----------------------/
1979	pLexOrder=FALSE;
1980	pMixedOrder=FALSE;
1981	pFDeg=pDeg;
1982	if (pOrdSgn == 1) pLDeg = ldegb;
1983	else pLDeg = ldeg0;
1984	/======== ordering type is (_,c) =========================/
1985	if ((orders[0]==ringorder_unspec)
1986	\|\|(
1987	((orders[1]==ringorder_c)\|\|(orders[1]==ringorder_C))
1988	&& (orders[0]!=ringorder_M)
1989	&& (orders[2]==0))
1990	)
1991	{
1992	if ((orders[0]!=ringorder_unspec)
1993	&& (orders[1]==ringorder_C))
1994	pComponentOrder=-1;
1995	if (pOrdSgn == -1) pLDeg = ldeg0c;
1996	#ifdef COMP_FAST
1997	SimpleChoose(orders[0],orders[1], &pComp0);
1998	#else
1999	SimpleChoose(orders[0],&pComp0);
2000	#endif
2001	SetpSetm(orders[0],0);
2002	#ifdef TEST_MAC_ORDER
2003	if (orders[0]==ringorder_dp)
2004	bBinomSet(orders);
2005	#endif
2006	}
2007	/======== ordering type is (c,_) =========================/
2008	else if (((orders[0]==ringorder_c)\|\|(orders[0]==ringorder_C))
2009	&& (orders[1]!=ringorder_M)
2010	&& (orders[2]==0))
2011	{
2012	/* pLDeg = ldeg0; is standard*/
2013	if (orders[0]==ringorder_C)
2014	pComponentOrder=-1;
2015	#ifdef COMP_FAST
2016	SimpleChooseC(orders[1],pComponentOrder, &pComp0);
2017	#else
2018	SimpleChooseC(orders[1],&pComp0);
2019	#endif
2020	SetpSetm(orders[1],1);
2021	#ifdef TEST_MAC_ORDER
2022	if (orders[1]==ringorder_dp)
2023	bBinomSet(orders);
2024	#endif
2025	}
2026	/------- more than one block ----------------------/
2027	else
2028	{
2029	//pLexOrder=TRUE;
2030	pVectorOut=orders[0]==ringorder_c;
2031	if ((pVectorOut)\|\|(orders[0]==ringorder_C))
2032	{
2033	if(block1[1]!=pVariables) pLexOrder=TRUE;
2034	}
2035	else
2036	{
2037	if(block1[0]!=pVariables) pLexOrder=TRUE;
2038	}
2039	/the number of orderings:/
2040	i = 0;
2041	while (orders[++i] != 0);
2042	do
2043	{
2044	i--;
2045	HighSet(i, orders[i]);/sets also pMixedOrder to TRUE, if.../
2046	SetpSetm(orders[i],i);
2047	}
2048	while (i != 0);
2049
2050	pComp0 = BlockComp;
2051	if ((orders[0]!=ringorder_c)&&(orders[0]!=ringorder_C))
2052	{
2053	pLDeg = ldeg1c;
2054	}
2055	else
2056	{
2057	pLDeg = ldeg1;
2058	}
2059	pFDeg = pWTotaldegree; // may be improved: pTotaldegree for lp/dp/ls/.. blocks
2060	}
2061	if ((pLexOrder) \|\| (pOrdSgn==-1))
2062	{
2063	test &= ~Sy_bit(OPT_REDTAIL); /* noredTail */
2064	}
2065	#ifdef COMP_TRADITIONAL
2066	if (t_pComp0 == NULL)
2067	t_pComp0 = pComp0;
2068	#endif
2069
2070	#ifdef COMP_FAST
2071	if (f_pComp0 == NULL)
2072	f_pComp0 = pComp0;
2073	#endif
2074
2075	#ifdef DIV_COUNT
2076	InitDivCount(pVariables);
2077	#endif
2078	}
2079
2080	/* -------------------------------------------------------- */
2081
2082	static BOOLEAN pMultT_nok;
2083	/*2
2084	* update the polynomial a by multipying it by
2085	* the (number) coefficient
2086	* and the exponent vector (of) exp (a well initialized polynomial)
2087	*/
2088	poly pMultT(poly a, poly exp )
2089	{
2090	int i;
2091	number t,x,y=pGetCoeff(exp);
2092	poly aa=a;
2093	poly prev=NULL;
2094	#ifdef SDRING
2095	poly pDRINGres=NULL;
2096	#endif
2097
2098	pMultT_nok = pGetComp(exp);
2099	#ifdef PDEBUG
2100	pTest(a);
2101	pTest(exp);
2102	#endif
2103	while (a !=NULL)
2104	{
2105	x=pGetCoeff(a);
2106	t=nMult(x/pGetCoeff(a)/,y/pGetCoeff(exp)/);
2107	nDelete(&x/pGetCoeff(a)/);
2108	pSetCoeff0(a,t);
2109	if (nIsZero(t))
2110	{
2111	if (prev==NULL) { pDelete1(&a); aa=a; }
2112	else { pDelete1(&prev->next); a=prev->next;}
2113	}
2114	else
2115	{
2116	#ifdef DRING
2117	if (pDRING)
2118	{
2119	if (pdDFlag(a)==1)
2120	{
2121	if (pdDFlag(exp)==1)
2122	{
2123	pDRINGres=pAdd(pDRINGres,pMultDD(a,exp));
2124	}
2125	else
2126	{
2127	pDRINGres=pAdd(pDRINGres,pMultDT(a,exp));
2128	}
2129	}
2130	else
2131	{
2132	if (pdDFlag(exp)==1)
2133	{
2134	pDRINGres=pAdd(pDRINGres,pMultDD(a,exp));
2135	}
2136	else
2137	{
2138	pDRINGres=pAdd(pDRINGres,pMultTT(a,exp));
2139	}
2140	}
2141	}
2142	else
2143	#endif
2144	#ifdef SRING
2145	if (pSRING)
2146	{
2147	pDRINGres=pAdd(pDRINGres,psMultM(a,exp));
2148	}
2149	else
2150	#endif
2151	{
2152	for (i=pVariables; i != 0; i--)
2153	{
2154	pAddExp(a,i, pGetExp(exp,i));
2155	}
2156	#ifdef TEST_MAC_ORDER
2157	if (bNoAdd)
2158	pSetm(a);
2159	else
2160	#endif
2161	a->Order += exp->Order;
2162	if (pMultT_nok)
2163	{
2164	if (pGetComp(a) == 0)
2165	{
2166	pSetComp(a, pGetComp(exp));
2167	}
2168	else
2169	{
2170	return NULL /FALSE/;
2171	}
2172	}
2173	}
2174	prev=a;
2175	pIter(a);
2176	}
2177	}
2178	pMultT_nok=0;
2179	#ifdef DRING
2180	if (pDRING)
2181	{
2182	pDelete(&aa);
2183	return pDRINGres;
2184	}
2185	#endif
2186	#ifdef SRING
2187	if (pSRING)
2188	{
2189	pDelete(&aa);
2190	return pDRINGres;
2191	}
2192	#endif
2193	return aa; /TRUE/
2194	}
2195
2196	/*2
2197	* multiply p1 with p2, p1 and p2 are destroyed
2198	* do not put attention on speed: the procedure is only used in the interpreter
2199	*/
2200	poly pMult(poly p1, poly p2)
2201	{
2202	poly res, r, rn, a;
2203	BOOLEAN cont;
2204
2205	if ((p1!=NULL) && (p2!=NULL))
2206	{
2207	#ifdef PDEBUG
2208	pTest(p1);
2209	pTest(p2);
2210	#endif
2211	cont = TRUE;
2212	a = p1;
2213	if (pNext(p2)!=NULL)
2214	a = pCopy(a);
2215	else
2216	cont = FALSE;
2217	res = pMultT(a, p2);
2218	if (pMultT_nok)
2219	{
2220	if (cont) pDelete(&p1);
2221	pDelete(&res);
2222	pDelete(&p2);
2223	return NULL;
2224	}
2225	pDelete1(&p2);
2226	r = res;
2227	if (r!=NULL) rn = pNext(r);
2228	else rn=NULL;
2229	while (cont)
2230	{
2231	if (pNext(p2)==NULL)
2232	{
2233	a = p1;
2234	cont = FALSE;
2235	}
2236	else
2237	{
2238	a = pCopy(p1);
2239	}
2240	a=pMultT(a, p2); //sets pMultT_nok
2241	if (pMultT_nok)
2242	{
2243	if (cont) pDelete(&p1);
2244	pDelete(&a);
2245	pDelete(&res);
2246	pDelete(&p2);
2247	return NULL;
2248	}
2249	while ((rn!=NULL) && (pComp0(rn,a)>0))
2250	{
2251	r = rn;
2252	pIter(rn);
2253	}
2254	if (r!=NULL) pNext(r) = rn = pAdd(a, rn);
2255	else res=r=a;
2256	pDelete1(&p2);
2257	}
2258	return res;
2259	}
2260	pDelete(&p1);
2261	pDelete(&p2);
2262	return NULL;
2263	}
2264
2265	/*2
2266	* update a by multiplying it with c (c will not be destroyed)
2267	*/
2268	void pMultN(poly a, number c)
2269	{
2270	number t;
2271
2272	while (a!=NULL)
2273	{
2274	t=nMult(pGetCoeff(a), c);
2275	//nNormalize(t);
2276	pSetCoeff(a,t);
2277	pIter(a);
2278	}
2279	}
2280
2281	/*2
2282	* return a copy of the poly a times the number c (a,c will not be destroyed)
2283	*/
2284	poly pMultCopyN(poly a, number c)
2285	{
2286	poly result=NULL,hp;
2287
2288	if (a != NULL)
2289	{
2290	result=pNew();
2291	pCopy2(result,a);
2292	pNext(result)=NULL;
2293	pGetCoeff(result)=nMult(pGetCoeff(a),c);
2294	pIter(a);
2295	hp=result;
2296	while (a)
2297	{
2298	hp=pNext(hp)=pNew();
2299	pCopy2(hp,a);
2300	pSetCoeff0(hp,nMult(pGetCoeff(a), c));
2301	pIter(a);
2302	}
2303	pNext(hp)=NULL;
2304	}
2305	return result;
2306	}
2307
2308	/*2
2309	* returns TRUE if the head term of b is a multiple of the head term of a
2310	*/
2311	#if defined(macintosh)
2312	BOOLEAN pDivisibleBy(poly a, poly b)
2313	{
2314	if ((a!=NULL)&&(( pGetComp(a)==0) \|\| ( pGetComp(a) == pGetComp(b))))
2315	{
2316	int i=pVariables;
2317	Exponent_t *e1=&( pGetExp(a,1));
2318	Exponent_t *e2=&( pGetExp(b,1));
2319	if ((e1) > (e2)) return FALSE;
2320	do
2321	{
2322	i--;
2323	if (i == 0) return TRUE;
2324	e1++;
2325	e2++;
2326	} while ((e1) <= (e2));
2327	}
2328	return FALSE;
2329	}
2330	#endif
2331
2332
2333	/*2
2334	* assumes that the head term of b is a multiple of the head term of a
2335	* and return the multiplicant *m
2336	*/
2337	poly pDivide(poly a, poly b)
2338	{
2339	int i;
2340	poly result=pInit();
2341
2342	for(i=(int)pVariables; i; i--)
2343	pSetExp(result,i, pGetExp(a,i)- pGetExp(b,i));
2344	pSetComp(result, pGetComp(a) - pGetComp(b));
2345	pSetm(result);
2346	return result;
2347	}
2348
2349	/*2
2350	* divides a by the monomial b, ignores monomials wihich are not divisible
2351	* assumes that b is not NULL
2352	*/
2353	poly pDivideM(poly a, poly b)
2354	{
2355	if (a==NULL) return NULL;
2356	poly result=a;
2357	poly prev=NULL;
2358	int i;
2359	number inv=nInvers(pGetCoeff(b));
2360
2361	while (a!=NULL)
2362	{
2363	if (pDivisibleBy(b,a))
2364	{
2365	for(i=(int)pVariables; i; i--)
2366	pSubExp(a,i, pGetExp(b,i));
2367	pSubComp(a, pGetComp(b));
2368	pSetm(a);
2369	prev=a;
2370	pIter(a);
2371	}
2372	else
2373	{
2374	if (prev==NULL)
2375	{
2376	pDelete1(&result);
2377	a=result;
2378	}
2379	else
2380	{
2381	pDelete1(&pNext(prev));
2382	a=pNext(prev);
2383	}
2384	}
2385	}
2386	pMultN(result,inv);
2387	nDelete(&inv);
2388	pDelete(&b);
2389	return result;
2390	}
2391
2392	/*2
2393	* returns the LCM of the head terms of a and b in *m
2394	*/
2395	void pLcm(poly a, poly b, poly m)
2396	{
2397	int i;
2398	for (i=pVariables; i; i--)
2399	{
2400	pSetExp(m,i, max( pGetExp(a,i), pGetExp(b,i)));
2401	}
2402	pSetComp(m, max(pGetComp(a), pGetComp(b)));
2403	}
2404
2405	/*2
2406	* convert monomial given as string to poly, e.g. 1x3y5z
2407	*/
2408	poly pmInit(char *st, BOOLEAN &ok)
2409	{
2410	int i,j;
2411	ok=FALSE;
2412	BOOLEAN b=FALSE;
2413	poly rc = pInit();
2414	char *s = nRead(st,&(rc->coef));
2415	if (s==st)
2416	/* i.e. it does not start with a coeff: test if it is a ringvar*/
2417	{
2418	j = rIsRingVar(s);
2419	if (j >= 0)
2420	{
2421	pIncrExp(rc,1+j);
2422	goto done;
2423	}
2424	}
2425	else
2426	b=TRUE;
2427	while (*s!='\0')
2428	{
2429	char ss[2];
2430	ss[0] = *s++;
2431	ss[1] = '\0';
2432	j = rIsRingVar(ss);
2433	if (j >= 0)
2434	{
2435	s = eati(s,&i);
2436	pAddExp(rc,1+j, (Exponent_t)i);
2437	}
2438	else
2439	{
2440	if ((s!=st)&&isdigit(st[0]))
2441	{
2442	errorreported=TRUE;
2443	}
2444	pDelete(&rc);
2445	return NULL;
2446	}
2447	}
2448	done:
2449	ok=!errorreported;
2450	if (nIsZero(pGetCoeff(rc))) pDelete1(&rc);
2451	else
2452	{
2453	#ifdef DRING
2454	if (pDRING)
2455	{
2456	for(i=1;i<=pdN;i++)
2457	{
2458	if(pGetExp(rc,pdDX(i))>0)
2459	{
2460	pdDFlag(rc)=1;
2461	break;
2462	}
2463	}
2464	}
2465	#endif
2466	pSetm(rc);
2467	}
2468	return rc;
2469	}
2470
2471	/*2
2472	*make p homgeneous by multiplying the monomials by powers of x_varnum
2473	*/
2474	poly pHomogen (poly p, int varnum)
2475	{
2476	poly q=NULL;
2477	poly res;
2478	int o,ii;
2479
2480	if (p!=NULL)
2481	{
2482	if ((varnum < 1) \|\| (varnum > pVariables))
2483	{
2484	return NULL;
2485	}
2486	o=pWTotaldegree(p);
2487	q=pNext(p);
2488	while (q != NULL)
2489	{
2490	ii=pWTotaldegree(q);
2491	if (ii>o) o=ii;
2492	pIter(q);
2493	}
2494	q = pCopy(p);
2495	res = q;
2496	while (q != NULL)
2497	{
2498	ii = o-pWTotaldegree(q);
2499	if (ii!=0)
2500	{
2501	pAddExp(q,varnum, (Exponent_t)ii);
2502	pSetm(q);
2503	}
2504	pIter(q);
2505	}
2506	q = pOrdPolyInsertSetm(res);
2507	}
2508	return q;
2509	}
2510
2511
2512	/*2
2513	*replaces the maximal powers of the leading monomial of p2 in p1 by
2514	*the same powers of n, utility for dehomogenization
2515	*/
2516	poly pDehomogen (poly p1,poly p2,number n)
2517	{
2518	polyset P;
2519	int SizeOfSet=5;
2520	int i;
2521	poly p;
2522	number nn;
2523
2524	P = (polyset)Alloc(5*sizeof(poly));
2525	for (i=0; i<5; i++)
2526	{
2527	P[i] = NULL;
2528	}
2529	pCancelPolyByMonom(p1,p2,&P,&SizeOfSet);
2530	p = P[0];
2531	//P[0] = NULL ;// for safety, may be remoeved later
2532	for (i=1; i<SizeOfSet; i++)
2533	{
2534	if (P[i] != NULL)
2535	{
2536	nPower(n,i,&nn);
2537	pMultN(P[i],nn);
2538	p = pAdd(p,P[i]);
2539	//P[i] =NULL; // for safety, may be removed later
2540	nDelete(&nn);
2541	}
2542	}
2543	Free((ADDRESS)P,SizeOfSet*sizeof(poly));
2544	return p;
2545	}
2546
2547	/*4
2548	*Returns the exponent of the maximal power of the leading monomial of
2549	*p2 in that of p1
2550	*/
2551	static int pGetMaxPower (poly p1,poly p2)
2552	{
2553	int i,k,res = 32000; /a very large integer/
2554
2555	if (p1 == NULL) return 0;
2556	for (i=1; i<=pVariables; i++)
2557	{
2558	if ( pGetExp(p2,i) != 0)
2559	{
2560	k = pGetExp(p1,i) / pGetExp(p2,i);
2561	if (k < res) res = k;
2562	}
2563	}
2564	return res;
2565	}
2566
2567	/*2
2568	*Returns as i-th entry of P the coefficient of the (i-1) power of
2569	*the leading monomial of p2 in p1
2570	*/
2571	void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet)
2572	{
2573	int maxPow;
2574	poly p,qp,Coeff;
2575
2576	if (*P == NULL)
2577	{
2578	P = (polyset) Alloc(5sizeof(poly));
2579	*SizeOfSet = 5;
2580	}
2581	p = pCopy(p1);
2582	while (p != NULL)
2583	{
2584	qp = p->next;
2585	p->next = NULL;
2586	maxPow = pGetMaxPower(p,p2);
2587	Coeff = pDivByMonom(p,p2);
2588	if (maxPow > *SizeOfSet)
2589	{
2590	pEnlargeSet(P,SizeOfSet,maxPow+1-SizeOfSet);
2591	*SizeOfSet = maxPow+1;
2592	}
2593	(P)[maxPow] = pAdd((P)[maxPow],Coeff);
2594	pDelete(&p);
2595	p = qp;
2596	}
2597	}
2598
2599	/*2
2600	*returns the leading monomial of p1 divided by the maximal power of that
2601	*of p2
2602	*/
2603	poly pDivByMonom (poly p1,poly p2)
2604	{
2605	int k, i;
2606
2607	if (p1 == NULL) return NULL;
2608	k = pGetMaxPower(p1,p2);
2609	if (k == 0)
2610	return pHead(p1);
2611	else
2612	{
2613	number n;
2614	poly p = pInit();
2615
2616	p->next = NULL;
2617	for (i=1; i<=pVariables; i++)
2618	{
2619	pSetExp(p,i, pGetExp(p1,i)-k* pGetExp(p2,i));
2620	}
2621	nPower(p2->coef,k,&n);
2622	pSetCoeff0(p,nDiv(p1->coef,n));
2623	nDelete(&n);
2624	pSetm(p);
2625	return p;
2626	}
2627	}
2628	/----------utilities for syzygies--------------/
2629	poly pTakeOutComp(poly * p, int k)
2630	{
2631	poly q = *p,qq=NULL,result = NULL;
2632
2633	if (q==NULL) return NULL;
2634	if (pGetComp(q)==k)
2635	{
2636	result = q;
2637	while ((q!=NULL) && (pGetComp(q)==k))
2638	{
2639	pSetComp(q,0);
2640	qq = q;
2641	pIter(q);
2642	}
2643	*p = q;
2644	pNext(qq) = NULL;
2645	}
2646	if (q==NULL) return result;
2647	if (pGetComp(q) > k) pDecrComp(q);
2648	poly pNext_q;
2649	while ((pNext_q=pNext(q))!=NULL)
2650	{
2651	if (pGetComp(pNext_q)==k)
2652	{
2653	if (result==NULL)
2654	{
2655	result = pNext_q;
2656	qq = result;
2657	}
2658	else
2659	{
2660	pNext(qq) = pNext_q;
2661	pIter(qq);
2662	}
2663	pNext(q) = pNext(pNext_q);
2664	pNext(qq) =NULL;
2665	pSetComp(qq,0);
2666	}
2667	else
2668	{
2669	/pIter(q);/ q=pNext_q;
2670	if (pGetComp(q) > k) pDecrComp(q);
2671	}
2672	}
2673	return result;
2674	}
2675
2676	poly pTakeOutComp1(poly * p, int k)
2677	{
2678	poly q = *p;
2679
2680	if (q==NULL) return NULL;
2681
2682	poly qq=NULL,result = NULL;
2683
2684	if (pGetComp(q)==k)
2685	{
2686	result = q; /* p /
2687	while ((q!=NULL) && (pGetComp(q)==k))
2688	{
2689	pSetComp(q,0);
2690	qq = q;
2691	pIter(q);
2692	}
2693	*p = q;
2694	pNext(qq) = NULL;
2695	}
2696	if (q==NULL) return result;
2697	// if (pGetComp(q) > k) pGetComp(q)--;
2698	while (pNext(q)!=NULL)
2699	{
2700	if (pGetComp(pNext(q))==k)
2701	{
2702	if (result==NULL)
2703	{
2704	result = pNext(q);
2705	qq = result;
2706	}
2707	else
2708	{
2709	pNext(qq) = pNext(q);
2710	pIter(qq);
2711	}
2712	pNext(q) = pNext(pNext(q));
2713	pNext(qq) =NULL;
2714	pSetComp(qq,0);
2715	}
2716	else
2717	{
2718	pIter(q);
2719	// if (pGetComp(q) > k) pGetComp(q)--;
2720	}
2721	}
2722	return result;
2723	}
2724
2725	void pDeleteComp(poly * p,int k)
2726	{
2727	poly q;
2728
2729	while ((p!=NULL) && (pGetComp(p)==k)) pDelete1(p);
2730	if (*p==NULL) return;
2731	q = *p;
2732	if (pGetComp(q)>k) pDecrComp(q);
2733	while (pNext(q)!=NULL)
2734	{
2735	if (pGetComp(pNext(q))==k)
2736	pDelete1(&(pNext(q)));
2737	else
2738	{
2739	pIter(q);
2740	if (pGetComp(q)>k) pDecrComp(q);
2741	}
2742	}
2743	}
2744	/----------end of utilities for syzygies--------------/
2745
2746	/*2
2747	* pair has no common factor ? or is no polynomial
2748	*/
2749	#ifdef COMP_FAST
2750	BOOLEAN pHasNotCF(poly p1, poly p2)
2751	{
2752	#ifdef SRING
2753	if (pSRING)
2754	return FALSE;
2755	#endif
2756
2757	if (pGetComp(p1) > 0 \|\| pGetComp(p2) > 0)
2758	return FALSE;
2759	Exponent_t * m1 = &(p1->exp[pVarLowIndex]);
2760	Exponent_t * m2 = &(p2->exp[pVarLowIndex]);
2761	int i = 1;
2762	loop
2763	{
2764	if (((m1) > 0) && ((m2) > 0))
2765	return FALSE;
2766	if (i == pVariables)
2767	return TRUE;
2768	m1++;m2++;
2769	i++;
2770	}
2771	}
2772	#else
2773	BOOLEAN pHasNotCF(poly p1, poly p2)
2774	{
2775	#ifdef SRING
2776	if (pSRING)
2777	return FALSE;
2778	#endif
2779	short * m1 = p1->exp;
2780	short * m2 = p2->exp;
2781
2782	if (((m1) > 0)\|\|((m2) > 0))
2783	return FALSE;
2784	int i = 1;
2785	loop
2786	{
2787	m1++;m2++;
2788	if (((m1) > 0) && ((m2) > 0))
2789	return FALSE;
2790	if (i == pVariables)
2791	return TRUE;
2792	i++;
2793	}
2794	}
2795	#endif
2796	/*
2797	*void pSFactors(poly f, poly g, poly a, poly b)
2798	*{
2799	* int i,d;
2800	*
2801	* for (i=pVariables;i>0;i--)
2802	* {
2803	* d = pGetExp(f,i)- pGetExp(g,i);
2804	* if (d >= 0)
2805	* {
2806	* pGetExp(a,i) = 0;
2807	* pGetExp(b,i) = d;
2808	* }
2809	* else
2810	* {
2811	* pGetExp(a,i) = -d;
2812	* pGetExp(b,i) = 0;
2813	* }
2814	* }
2815	* pGetComp(a) = 0;
2816	* pGetComp(b) = 0;
2817	* pSetm(a);
2818	* pSetm(b);
2819	*}
2820	*/
2821
2822	/*
2823	*void pSDiv(poly f, poly g, poly b)
2824	*{
2825	* int i,d;
2826	*
2827	* for (i=pVariables;i>0;i--)
2828	* {
2829	* d = pGetExp(f,i)- pGetExp(g,i);
2830	* pGetExp(b,i) = d;
2831	* }
2832	* pGetComp(b) = 0;
2833	* pSetm(b);
2834	*}
2835	*/
2836
2837	/*2
2838	* update the initial term of a polynomial a by multipying it by
2839	* the (number) coefficient
2840	* and the exponent vector (of) exp (a well initialized polynomial)
2841	*/
2842	/*
2843	*void pSMultBy(poly f, poly m)
2844	*{
2845	* number t;
2846	* int i;
2847	* // short notok;
2848	*
2849	* t=nMult(f->coef, m->coef);
2850	* nDelete(&(f->coef));
2851	* f->coef = t;
2852	* f->Order += m->Order;
2853	* for (i=pVariables; i; i--)
2854	* pGetExp(f,i) += pGetExp(m,i);
2855	* // if (notok)
2856	* // {
2857	* if (!( pGetComp(f)))
2858	* {
2859	* pGetComp(f) = pGetComp(m);
2860	* }
2861	* // else
2862	* // {
2863	* // HALT;
2864	* // }
2865	* // }
2866	*}
2867	*/
2868
2869
2870	/*2
2871	*should return 1 if p divides q and p<q,
2872	* -1 if q divides p and q<p
2873	* 0 otherwise
2874	*/
2875	#ifdef COMP_FAST
2876	int pDivComp(poly p, poly q)
2877	{
2878	if (pGetComp(p) == pGetComp(q))
2879	{
2880	Exponent_t * mp = &(p->exp[pVarLowIndex]);
2881	Exponent_t * mq = &(q->exp[pVarLowIndex]);
2882	int i=pVariables;
2883	BOOLEAN a=FALSE, b=FALSE;
2884	for (; i>0; i--)
2885	{
2886	if (mp<mq)
2887	{
2888	if (b) return 0;
2889	a =TRUE;
2890	}
2891	else if (mp>mq)
2892	{
2893	if (a) return 0;
2894	b = TRUE;
2895	}
2896	mp++;mq++;
2897	}
2898	if (a) return 1;
2899	else if (b) return -1;
2900	}
2901	return 0;
2902	}
2903	#else
2904	int pDivComp(poly p, poly q)
2905	{
2906	short * mp = p->exp;
2907	short * mq = q->exp;
2908
2909	if (mp==mq)
2910	{
2911	int i=pVariables;
2912	BOOLEAN a=FALSE, b=FALSE;
2913	for (; i>0; i--)
2914	{
2915	mp++;mq++;
2916	if (mp<mq)
2917	{
2918	if (b) return 0;
2919	a =TRUE;
2920	}
2921	else if (mp>mq)
2922	{
2923	if (a) return 0;
2924	b = TRUE;
2925	}
2926	}
2927	if (a) return 1;
2928	else if (b) return -1;
2929	}
2930	return 0;
2931	}
2932	#endif
2933	/*2
2934	*divides p1 by its leading monomial
2935	*/
2936	void pNorm(poly p1)
2937	{
2938	poly h;
2939	number k, c;
2940
2941	if (p1!=NULL)
2942	{
2943	if (!nIsOne(pGetCoeff(p1)))
2944	{
2945	nNormalize(pGetCoeff(p1));
2946	k=pGetCoeff(p1);
2947	c = nInit(1);
2948	pSetCoeff0(p1,c);
2949	h = pNext(p1);
2950	while (h!=NULL)
2951	{
2952	c=nDiv(pGetCoeff(h),k);
2953	if (!nIsOne(c)) nNormalize(c);
2954	pSetCoeff(h,c);
2955	pIter(h);
2956	}
2957	nDelete(&k);
2958	}
2959	else
2960	{
2961	h = pNext(p1);
2962	while (h!=NULL)
2963	{
2964	nNormalize(pGetCoeff(h));
2965	pIter(h);
2966	}
2967	}
2968	}
2969	}
2970
2971	/*2
2972	*normalize all coeffizients
2973	*/
2974	void pNormalize(poly p)
2975	{
2976	while (p!=NULL)
2977	{
2978	nTest(pGetCoeff(p));
2979	nNormalize(pGetCoeff(p));
2980	pIter(p);
2981	}
2982	}
2983
2984	/*3
2985	* substitute the n-th variable by 1 in p
2986	* destroy p
2987	*/
2988	static poly pSubst1 (poly p,int n)
2989	{
2990	poly qq,result = NULL;
2991
2992	while (p != NULL)
2993	{
2994	qq = p;
2995	pIter(p);
2996	qq->next = NULL;
2997	pSetExp(qq,n,0);
2998	pSetm(qq);
2999	result = pAdd(result,qq);
3000	}
3001	return result;
3002	}
3003
3004	/*2
3005	* substitute the n-th variable by e in p
3006	* destroy p
3007	*/
3008	poly pSubst(poly p, int n, poly e)
3009	{
3010	if ((e!=NULL)&&(pIsConstant(e))&&(nIsOne(pGetCoeff(e))))
3011	return pSubst1(p,n);
3012
3013	int exponent,i;
3014	poly h, res, m;
3015	Exponent_t me,ee;
3016	number nu,nu1;
3017
3018	me=(Exponent_t )Alloc((pVariables+1)sizeof(Exponent_t));
3019	ee=(Exponent_t )Alloc((pVariables+1)sizeof(Exponent_t));
3020	if (e!=NULL) pGetExpV(e,ee);
3021	res=NULL;
3022	h=p;
3023	while (h!=NULL)
3024	{
3025	if ((e!=NULL) \|\| (pGetExp(h,n)==0))
3026	{
3027	m=pHead(h);
3028	pGetExpV(m,me);
3029	exponent=me[n];
3030	me[n]=0;
3031	for(i=1;i<=pVariables;i++)
3032	me[i]+=exponent*ee[i];
3033	pSetExpV(m,me);
3034	if (e!=NULL)
3035	{
3036	nPower(pGetCoeff(e),exponent,&nu);
3037	nu1=nMult(pGetCoeff(m),nu);
3038	nDelete(&nu);
3039	pSetCoeff(m,nu1);
3040	}
3041	res=pAdd(res,m);
3042	}
3043	pDelete1(&h);
3044	}
3045	Free((ADDRESS)me,(pVariables+1)*sizeof(Exponent_t));
3046	Free((ADDRESS)ee,(pVariables+1)*sizeof(Exponent_t));
3047	return res;
3048	}
3049
3050	BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
3051	{
3052	int k, j;
3053
3054	if (lcm==NULL) return FALSE;
3055
3056	for (j=pVariables; j; j--)
3057	if ( pGetExp(p,j) > pGetExp(lcm,j)) return FALSE;
3058	if ( pGetComp(p) != pGetComp(lcm)) return FALSE;
3059	for (j=pVariables; j; j--)
3060	{
3061	if (pGetExp(p1,j)!=pGetExp(lcm,j))
3062	{
3063	if (pGetExp(p,j)!=pGetExp(lcm,j))
3064	{
3065	for (k=pVariables; k>j; k--)
3066	{
3067	if ((pGetExp(p,k)!=pGetExp(lcm,k))
3068	&& (pGetExp(p2,k)!=pGetExp(lcm,k)))
3069	return TRUE;
3070	}
3071	for (k=j-1; k; k--)
3072	{
3073	if ((pGetExp(p,k)!=pGetExp(lcm,k))
3074	&& (pGetExp(p2,k)!=pGetExp(lcm,k)))
3075	return TRUE;
3076	}
3077	return FALSE;
3078	}
3079	}
3080	else if (pGetExp(p2,j)!=pGetExp(lcm,j))
3081	{
3082	if (pGetExp(p,j)!=pGetExp(lcm,j))
3083	{
3084	for (k=pVariables; k>j; k--)
3085	{
3086	if ((pGetExp(p,k)!=pGetExp(lcm,k))
3087	&& (pGetExp(p1,k)!=pGetExp(lcm,k)))
3088	return TRUE;
3089	}
3090	for (k=j-1; k; k--)
3091	{
3092	if ((pGetExp(p,k)!=pGetExp(lcm,k))
3093	&& (pGetExp(p1,k)!=pGetExp(lcm,k)))
3094	return TRUE;
3095	}
3096	return FALSE;
3097	}
3098	}
3099	}
3100	return FALSE;
3101	}
3102
3103	int pWeight(int i)
3104	{
3105	if ((firstwv==NULL) \|\| (i>firstBlockEnds))
3106	{
3107	return 1;
3108	}
3109	return firstwv[i-1];
3110	}
3111
3112	#ifdef COMP_DEBUG
3113	static int debug_comp(poly p1, poly p2)
3114	{
3115	int t_d = t_pComp0(p1, p2);
3116	int f_d = f_pComp0(p1, p2);
3117
3118	if (t_d != f_d)
3119	{
3120	fprintf(stderr, "Error in comp1lpc\n");
3121	t_pComp0(p1, p2);
3122	f_pComp0(p1, p2);
3123	}
3124	return t_d;
3125	}
3126	#endif
3127
3128	#ifdef COMP_STATISTICS
3129	static int s_comp_lp_c_1(poly p1, poly p2)
3130	{
3131	MonomCountTotal++;
3132
3133	OrderCmp(p2, p1,
3134	{MonomCountOrderS++; return -1;},
3135	{MonomCountOrderG++; return 1;});
3136
3137	CompCmp(p2, p1,
3138	{MonomCountCompS++; return -pComponentOrder},
3139	{MonomCountCompG++; return pComponentOrder});
3140
3141	MonomCountEqual++;
3142	return 0;
3143	}
3144	#endif
3145
3146	#ifdef COMP_STATISTICS
3147	static int s_comp_lp_c_i(poly p1, poly p2)
3148	{
3149	MonomCountTotal++;
3150
3151	OrderCmp(p2, p1,
3152	{MonomCountOrderS++; return -1;},
3153	{MonomCountOrderG++; return 1;});
3154
3155	int i = SIZEOF_ORDER / SIZEOF_EXPONENT + 1;
3156	Exponent_t d;
3157
3158
3159	for (;;)
3160	{
3161	d = pGetExp(p1, i) - pGetExp(p2, i);
3162	if (d)
3163	{
3164	MonomCountExp[i]++;
3165	if (d < 0)
3166	{
3167	return -1;
3168	MonomCountExpS++;
3169	}
3170	MonomCountExpG++;
3171	return 1;
3172	}
3173	i++;
3174	if (i == pVariables) break;
3175	}
3176
3177	CompCmp(p2, p1,
3178	{MonomCountCompS++; return -pComponentOrder},
3179	{MonomCountCompG++; return pComponentOrder});
3180
3181	MonomCountEqual++;
3182	return 0;
3183	}
3184	#endif
3185
3186	#ifdef MONOM_COUNT
3187	static void InitMonomCount(unsigned int nvars)
3188	{
3189	if (gMonomCount == NULL)
3190	{
3191	gMonomCount = (unsigned long ) Alloc0(nvarssizeof(int));
3192	gMonomCountLength = nvars;
3193	}
3194	else if (gMonomCountLength < nvars)
3195	{
3196	int i;
3197	gMonomCount = (unsigned long *) ReAlloc(gMonomCount,
3198	gMonomCountLength*sizeof(int),
3199	nvars*sizeof(int));
3200	for (i = gMonomCountLength; i< nvars; i++)
3201	gMonomCount[i] = 0;
3202	gMonomCountLength = nvars;
3203	}
3204	}
3205
3206	void OutputMonomCount()
3207	{
3208	unsigned int i;
3209	unsigned long max = gMonomCount0;
3210	float fmax;
3211	unsigned long sum = 0;
3212
3213
3214	// check that everything went ok
3215	if (gMonomCountS !=
3216	gMonomCount0 + gMonomCountOrder + gMonomCountL + gMonomCountG)
3217	{
3218	printf("MonomCountTotalError: %10lu %10lu\n", gMonomCountS,
3219	gMonomCount0 + gMonomCountOrder + gMonomCountL + gMonomCountG);
3220	}
3221	for (i=0; i<gMonomCountLength; i++)
3222	sum += gMonomCount[i];
3223	if (sum != gMonomCountL + gMonomCountG)
3224	{
3225	printf("MonomCountNotEqualError: %10lu %10lu\n",
3226	gMonomCountL + gMonomCountG, sum);
3227	}
3228
3229	printf("Total : %10lu\n", gMonomCountS);
3230	if (gMonomCountS == 0) gMonomCountS = 1;
3231
3232	printf("Order : %10lu \t %.4f\n",gMonomCountOrder,
3233	(float) gMonomCountOrder / (float) gMonomCountS);
3234	printf("Equal : %10lu \t %.4f\n",gMonomCount0,
3235	(float) gMonomCount0 / (float) gMonomCountS);
3236	printf("NotEqual: %10lu \t %.4f\n",gMonomCountL + gMonomCountG,
3237	((float) gMonomCountL + gMonomCountG) / (float) gMonomCountS);
3238	printf("\tGreater: %10lu \t %.4f\n",gMonomCountG,
3239	(float) gMonomCountG / (float) gMonomCountS);
3240	printf("\tLess : %10lu \t %.4f\n",gMonomCountL,
3241	(float) gMonomCountL / (float) gMonomCountS);
3242
3243	for (i=0; i<gMonomCountLength; i++)
3244	{
3245	printf("\t\t[%d] : %10lu \t %.4f\n", i, gMonomCount[i],
3246	(float) gMonomCount[i] / (float) gMonomCountS);
3247	}
3248	printf("E1Neg : %10lu \t %.4f\n",gMonomCountN1,
3249	(float) gMonomCountN1 / (float) gMonomCountS);
3250	printf("E2Neg : %10lu \t %.4f\n",gMonomCountN2,
3251	(float) gMonomCountN2 / (float) gMonomCountS);
3252	printf("BothNeg : %10lu \t %.4f\n",gMonomCountNN,
3253	(float) gMonomCountNN / (float) gMonomCountS);
3254	printf("Mcount2 : %u\n", gMcount2);
3255	}
3256
3257	void ResetMonomCount()
3258	{
3259	int i;
3260
3261	for (i=0; i<gMonomCountLength; i++)
3262	gMonomCount[i] = 0;
3263	gMonomCountOrder = 0;
3264	gMonomCount0 = 0;
3265	gMonomCountS = 0;
3266	gMonomCountG = 0;
3267	gMonomCountL = 0;
3268	gMonomCountN1 = 0;
3269	gMonomCountN2 = 0;
3270	gMonomCountNN = 0;
3271	gMcount2 = 0;
3272	}
3273
3274	#endif
3275	#ifdef MONOM_COUNT
3276	unsigned long gMonomCount0 = 0;
3277	unsigned long gMonomCountLength = 0;
3278	unsigned long gMonomCountOrder = 0;
3279	unsigned long gMonomCountG = 0;
3280	unsigned long gMonomCountL = 0;
3281	unsigned long gMonomCountS = 0;
3282	unsigned long gMonomCountN1 = 0;
3283	unsigned long gMonomCountN2 = 0;
3284	unsigned long gMonomCountNN = 0;
3285	unsigned long gMcount2 = 0;
3286	unsigned long* gMonomCount = NULL;
3287	#endif

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

Download in other formats: