Context Navigation

source: git/Singular/polys.cc @ a29c12

spielwiese

Last change on this file since a29c12 was a29c12, checked in by Wilfred Pohl <pohl@…>, 26 years ago
or -> o_r (mac) git-svn-id: file:///usr/local/Singular/svn/trunk@1094 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 69.2 KB

Line
1	/****************************************
2	* Computer Algebra System SINGULAR *
3	****************************************/
4	/* $Id: polys.cc,v 1.15 1998-01-27 16:20:42 pohl Exp $ */
5
6	/*
7	* ABSTRACT - all basic methods to manipulate polynomials
8	*/
9
10	/* includes */
11	#include <stdio.h>
12	#include <string.h>
13	#include <ctype.h>
14	#include "mod2.h"
15	#include "tok.h"
16	#include "mmemory.h"
17	#include "febase.h"
18	#include "numbers.h"
19	#include "polys.h"
20	#include "ring.h"
21	#include "binom.h"
22	#include "ipid.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(Exponent_t)*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 (o_r is the number of the ordering)
1289	*/
1290	#ifdef COMP_FAST
1291	static void SimpleChoose(int o_r, int comp_order, pCompProc *p)
1292	#else
1293	static void SimpleChoose(int o_r, 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(o_r)
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",o_r,__FILE__,__LINE__);
1330	#endif
1331	}
1332	}
1333	if (o_r == ringorder_lp \|\| o_r == ringorder_ls)
1334	{
1335	pLexOrder=TRUE;
1336	pFDeg = pTotaldegree;
1337	pLDeg = ldeg1c;
1338	if (o_r == ringorder_ls) pLexSgn = -1;
1339	}
1340
1341	*p = t_pComp0;
1342	#endif
1343	#ifdef COMP_FAST
1344	switch(o_r)
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 \|\| o_r == 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",o_r,__FILE__,__LINE__);
1415	#endif
1416	}
1417
1418	if (o_r == ringorder_lp \|\| o_r == ringorder_ls)
1419	{
1420	pLexOrder=TRUE;
1421	pFDeg = pTotaldegree;
1422	pLDeg = ldeg1c;
1423	if (o_r == 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 o_r = order of first block)
1437	*/
1438	static void SetpSetm(int o_r, int ip)
1439	{
1440	switch(o_r)
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",o_r,__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 (o_r is the number of the ordering)
1480	*/
1481	#ifdef COMP_FAST
1482	static void SimpleChooseC(int o_r, int comp_order, pCompProc *p)
1483	#else
1484	static void SimpleChooseC(int o_r, 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(o_r)
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",o_r,__FILE__,__LINE__);
1521	#endif
1522	}
1523	}
1524	if (o_r == ringorder_lp \|\| o_r == ringorder_ls)
1525	{
1526	pLexOrder=TRUE;
1527	pFDeg = pTotaldegree;
1528	pLDeg = ldeg1c;
1529	if (o_r == ringorder_ls) pLexSgn = -1;
1530	}
1531
1532	*p = t_pComp0;
1533	#endif
1534	#ifdef COMP_FAST
1535	switch(o_r)
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",o_r,__FILE__,__LINE__);
1577	#endif
1578	}
1579	if (o_r == ringorder_lp \|\| o_r == ringorder_ls)
1580	{
1581	pLexOrder=TRUE;
1582	pFDeg = pTotaldegree;
1583	pLDeg = ldeg1c;
1584	if (o_r == 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, o_r the number of the ordering)
1598	*/
1599	static void HighSet(int ip, int o_r)
1600	{
1601	switch(o_r)
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",o_r,__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
1926	void pChangeRing(int n, int Sgn, int * orders, int * b0, int * b1,
1927	short ** wv)
1928	{
1929	sip_sring tmpR;
1930	memset(&tmpR, 0, sizeof(sip_sring));
1931	tmpR.N = n;
1932	tmpR.OrdSgn = Sgn;
1933	tmpR.order = orders;
1934	tmpR.block0 = b0;
1935	tmpR.block1 = b1;
1936	tmpR.wvhdl = wv;
1937	pSetGlobals(&tmpR);
1938	}
1939
1940	void pSetGlobals(ring r, BOOLEAN complete)
1941	{
1942	#ifdef TEST_MAC_ORDER
1943	bNoAdd=FALSE;
1944	#endif
1945	int i;
1946	pComponentOrder=1;
1947	if (ppNoether!=NULL) pDelete(&ppNoether);
1948	#ifdef SRING
1949	pSRING=FALSE;
1950	pAltVars=r->N+1;
1951	#endif
1952	t_pComp0 = NULL;
1953	#ifdef COMP_FAST
1954	f_pComp0 = NULL;
1955	#endif
1956	pVariables = r->N;
1957
1958	// set the various size parameters and initialize memory
1959	pMonomSize = POLYSIZE + (pVariables + 1) * sizeof(Exponent_t);
1960	if ((pMonomSize % sizeof(void*)) == 0)
1961	{
1962	pMonomSizeW = pMonomSize/sizeof(void*);
1963	}
1964	else
1965	{
1966	pMonomSizeW = pMonomSize/sizeof(void*) + 1;
1967	pMonomSize = pMonomSizeWsizeof(void);
1968	}
1969
1970	#ifdef COMP_FAST
1971	if ((((pVariables+1)sizeof(Exponent_t)) % sizeof(void)) == 0)
1972	pVariables1W = (pVariables+1)sizeof(Exponent_t) / sizeof(void);
1973	else
1974	pVariables1W = ((pVariables+1)sizeof(Exponent_t) / sizeof(void)) + 1;
1975	if ((((pVariables)sizeof(Exponent_t)) % sizeof(void)) == 0)
1976	pVariablesW = (pVariables)sizeof(Exponent_t) / sizeof(void);
1977	else
1978	pVariablesW = ((pVariables)sizeof(Exponent_t) / sizeof(void)) + 1;
1979
1980	// Set default Var Indicies
1981	pSetVarIndicies(pVariables);
1982	#endif
1983	mmSpecializeBlock(pMonomSize);
1984
1985	pOrdSgn = r->OrdSgn;
1986	pVectorOut=(r->order[0]==ringorder_c);
1987	order=r->order;
1988	block0=r->block0;
1989	block1=r->block1;
1990	firstwv=NULL;
1991	polys_wv=r->wvhdl;
1992	/------- only one real block ----------------------/
1993	pLexOrder=FALSE;
1994	pMixedOrder=FALSE;
1995	pFDeg=pDeg;
1996	if (pOrdSgn == 1) pLDeg = ldegb;
1997	else pLDeg = ldeg0;
1998	/======== ordering type is (_,c) =========================/
1999	if ((order[0]==ringorder_unspec)
2000	\|\|(
2001	((order[1]==ringorder_c)\|\|(order[1]==ringorder_C))
2002	&& (order[0]!=ringorder_M)
2003	&& (order[2]==0))
2004	)
2005	{
2006	if ((order[0]!=ringorder_unspec)
2007	&& (order[1]==ringorder_C))
2008	pComponentOrder=-1;
2009	if (pOrdSgn == -1) pLDeg = ldeg0c;
2010	#ifdef COMP_FAST
2011	SimpleChoose(order[0],order[1], &pComp0);
2012	#else
2013	SimpleChoose(order[0],&pComp0);
2014	#endif
2015	SetpSetm(order[0],0);
2016	#ifdef TEST_MAC_ORDER
2017	if (order[0]==ringorder_dp)
2018	bBinomSet(order);
2019	#endif
2020	}
2021	/======== ordering type is (c,_) =========================/
2022	else if (((order[0]==ringorder_c)\|\|(order[0]==ringorder_C))
2023	&& (order[1]!=ringorder_M)
2024	&& (order[2]==0))
2025	{
2026	/* pLDeg = ldeg0; is standard*/
2027	if (order[0]==ringorder_C)
2028	pComponentOrder=-1;
2029	#ifdef COMP_FAST
2030	SimpleChooseC(order[1],pComponentOrder, &pComp0);
2031	#else
2032	SimpleChooseC(order[1],&pComp0);
2033	#endif
2034	SetpSetm(order[1],1);
2035	#ifdef TEST_MAC_ORDER
2036	if (order[1]==ringorder_dp)
2037	bBinomSet(order);
2038	#endif
2039	}
2040	/------- more than one block ----------------------/
2041	else
2042	{
2043	//pLexOrder=TRUE;
2044	pVectorOut=order[0]==ringorder_c;
2045	if ((pVectorOut)\|\|(order[0]==ringorder_C))
2046	{
2047	if(block1[1]!=pVariables) pLexOrder=TRUE;
2048	}
2049	else
2050	{
2051	if(block1[0]!=pVariables) pLexOrder=TRUE;
2052	}
2053	/the number of orderings:/
2054	i = 0;
2055	while (order[++i] != 0);
2056	do
2057	{
2058	i--;
2059	HighSet(i, order[i]);/sets also pMixedOrder to TRUE, if.../
2060	SetpSetm(order[i],i);
2061	}
2062	while (i != 0);
2063
2064	pComp0 = BlockComp;
2065	if ((order[0]!=ringorder_c)&&(order[0]!=ringorder_C))
2066	{
2067	pLDeg = ldeg1c;
2068	}
2069	else
2070	{
2071	pLDeg = ldeg1;
2072	}
2073	pFDeg = pWTotaldegree; // may be improved: pTotaldegree for lp/dp/ls/.. blocks
2074	}
2075	if ((pLexOrder) \|\| (pOrdSgn==-1))
2076	{
2077	test &= ~Sy_bit(OPT_REDTAIL); /* noredTail */
2078	}
2079	#ifdef COMP_TRADITIONAL
2080	if (t_pComp0 == NULL)
2081	t_pComp0 = pComp0;
2082	#endif
2083
2084	#ifdef COMP_FAST
2085	if (f_pComp0 == NULL)
2086	f_pComp0 = pComp0;
2087	#endif
2088
2089	#ifdef DIV_COUNT
2090	InitDivCount(pVariables);
2091	#endif
2092	}
2093
2094	/* -------------------------------------------------------- */
2095
2096	static BOOLEAN pMultT_nok;
2097	/*2
2098	* update the polynomial a by multipying it by
2099	* the (number) coefficient
2100	* and the exponent vector (of) exp (a well initialized polynomial)
2101	*/
2102	poly pMultT(poly a, poly exp )
2103	{
2104	int i;
2105	number t,x,y=pGetCoeff(exp);
2106	poly aa=a;
2107	poly prev=NULL;
2108	#ifdef SDRING
2109	poly pDRINGres=NULL;
2110	#endif
2111
2112	pMultT_nok = pGetComp(exp);
2113	#ifdef PDEBUG
2114	pTest(a);
2115	pTest(exp);
2116	#endif
2117	while (a !=NULL)
2118	{
2119	x=pGetCoeff(a);
2120	t=nMult(x/pGetCoeff(a)/,y/pGetCoeff(exp)/);
2121	nDelete(&x/pGetCoeff(a)/);
2122	pSetCoeff0(a,t);
2123	if (nIsZero(t))
2124	{
2125	if (prev==NULL) { pDelete1(&a); aa=a; }
2126	else { pDelete1(&prev->next); a=prev->next;}
2127	}
2128	else
2129	{
2130	#ifdef DRING
2131	if (pDRING)
2132	{
2133	if (pdDFlag(a)==1)
2134	{
2135	if (pdDFlag(exp)==1)
2136	{
2137	pDRINGres=pAdd(pDRINGres,pMultDD(a,exp));
2138	}
2139	else
2140	{
2141	pDRINGres=pAdd(pDRINGres,pMultDT(a,exp));
2142	}
2143	}
2144	else
2145	{
2146	if (pdDFlag(exp)==1)
2147	{
2148	pDRINGres=pAdd(pDRINGres,pMultDD(a,exp));
2149	}
2150	else
2151	{
2152	pDRINGres=pAdd(pDRINGres,pMultTT(a,exp));
2153	}
2154	}
2155	}
2156	else
2157	#endif
2158	#ifdef SRING
2159	if (pSRING)
2160	{
2161	pDRINGres=pAdd(pDRINGres,psMultM(a,exp));
2162	}
2163	else
2164	#endif
2165	{
2166	for (i=pVariables; i != 0; i--)
2167	{
2168	pAddExp(a,i, pGetExp(exp,i));
2169	}
2170	#ifdef TEST_MAC_ORDER
2171	if (bNoAdd)
2172	pSetm(a);
2173	else
2174	#endif
2175	a->Order += exp->Order;
2176	if (pMultT_nok)
2177	{
2178	if (pGetComp(a) == 0)
2179	{
2180	pSetComp(a, pGetComp(exp));
2181	}
2182	else
2183	{
2184	return NULL /FALSE/;
2185	}
2186	}
2187	}
2188	prev=a;
2189	pIter(a);
2190	}
2191	}
2192	pMultT_nok=0;
2193	#ifdef DRING
2194	if (pDRING)
2195	{
2196	pDelete(&aa);
2197	return pDRINGres;
2198	}
2199	#endif
2200	#ifdef SRING
2201	if (pSRING)
2202	{
2203	pDelete(&aa);
2204	return pDRINGres;
2205	}
2206	#endif
2207	return aa; /TRUE/
2208	}
2209
2210	/*2
2211	* multiply p1 with p2, p1 and p2 are destroyed
2212	* do not put attention on speed: the procedure is only used in the interpreter
2213	*/
2214	poly pMult(poly p1, poly p2)
2215	{
2216	poly res, r, rn, a;
2217	BOOLEAN cont;
2218
2219	if ((p1!=NULL) && (p2!=NULL))
2220	{
2221	#ifdef PDEBUG
2222	pTest(p1);
2223	pTest(p2);
2224	#endif
2225	cont = TRUE;
2226	a = p1;
2227	if (pNext(p2)!=NULL)
2228	a = pCopy(a);
2229	else
2230	cont = FALSE;
2231	res = pMultT(a, p2);
2232	if (pMultT_nok)
2233	{
2234	if (cont) pDelete(&p1);
2235	pDelete(&res);
2236	pDelete(&p2);
2237	return NULL;
2238	}
2239	pDelete1(&p2);
2240	r = res;
2241	if (r!=NULL) rn = pNext(r);
2242	else rn=NULL;
2243	while (cont)
2244	{
2245	if (pNext(p2)==NULL)
2246	{
2247	a = p1;
2248	cont = FALSE;
2249	}
2250	else
2251	{
2252	a = pCopy(p1);
2253	}
2254	a=pMultT(a, p2); //sets pMultT_nok
2255	if (pMultT_nok)
2256	{
2257	if (cont) pDelete(&p1);
2258	pDelete(&a);
2259	pDelete(&res);
2260	pDelete(&p2);
2261	return NULL;
2262	}
2263	while ((rn!=NULL) && (pComp0(rn,a)>0))
2264	{
2265	r = rn;
2266	pIter(rn);
2267	}
2268	if (r!=NULL) pNext(r) = rn = pAdd(a, rn);
2269	else res=r=a;
2270	pDelete1(&p2);
2271	}
2272	return res;
2273	}
2274	pDelete(&p1);
2275	pDelete(&p2);
2276	return NULL;
2277	}
2278
2279	/*2
2280	* update a by multiplying it with c (c will not be destroyed)
2281	*/
2282	void pMultN(poly a, number c)
2283	{
2284	number t;
2285
2286	while (a!=NULL)
2287	{
2288	t=nMult(pGetCoeff(a), c);
2289	//nNormalize(t);
2290	pSetCoeff(a,t);
2291	pIter(a);
2292	}
2293	}
2294
2295	/*2
2296	* return a copy of the poly a times the number c (a,c will not be destroyed)
2297	*/
2298	poly pMultCopyN(poly a, number c)
2299	{
2300	poly result=NULL,hp;
2301
2302	if (a != NULL)
2303	{
2304	result=pNew();
2305	pCopy2(result,a);
2306	pNext(result)=NULL;
2307	pGetCoeff(result)=nMult(pGetCoeff(a),c);
2308	pIter(a);
2309	hp=result;
2310	while (a)
2311	{
2312	hp=pNext(hp)=pNew();
2313	pCopy2(hp,a);
2314	pSetCoeff0(hp,nMult(pGetCoeff(a), c));
2315	pIter(a);
2316	}
2317	pNext(hp)=NULL;
2318	}
2319	return result;
2320	}
2321
2322	/*2
2323	* returns TRUE if the head term of b is a multiple of the head term of a
2324	*/
2325	#if defined(macintosh)
2326	BOOLEAN pDivisibleBy(poly a, poly b)
2327	{
2328	if ((a!=NULL)&&(( pGetComp(a)==0) \|\| ( pGetComp(a) == pGetComp(b))))
2329	{
2330	int i=pVariables;
2331	Exponent_t *e1=&( pGetExp(a,1));
2332	Exponent_t *e2=&( pGetExp(b,1));
2333	if ((e1) > (e2)) return FALSE;
2334	do
2335	{
2336	i--;
2337	if (i == 0) return TRUE;
2338	e1++;
2339	e2++;
2340	} while ((e1) <= (e2));
2341	}
2342	return FALSE;
2343	}
2344	#endif
2345
2346
2347	/*2
2348	* assumes that the head term of b is a multiple of the head term of a
2349	* and return the multiplicant *m
2350	*/
2351	poly pDivide(poly a, poly b)
2352	{
2353	int i;
2354	poly result=pInit();
2355
2356	for(i=(int)pVariables; i; i--)
2357	pSetExp(result,i, pGetExp(a,i)- pGetExp(b,i));
2358	pSetComp(result, pGetComp(a) - pGetComp(b));
2359	pSetm(result);
2360	return result;
2361	}
2362
2363	/*2
2364	* divides a by the monomial b, ignores monomials wihich are not divisible
2365	* assumes that b is not NULL
2366	*/
2367	poly pDivideM(poly a, poly b)
2368	{
2369	if (a==NULL) return NULL;
2370	poly result=a;
2371	poly prev=NULL;
2372	int i;
2373	number inv=nInvers(pGetCoeff(b));
2374
2375	while (a!=NULL)
2376	{
2377	if (pDivisibleBy(b,a))
2378	{
2379	for(i=(int)pVariables; i; i--)
2380	pSubExp(a,i, pGetExp(b,i));
2381	pSubComp(a, pGetComp(b));
2382	pSetm(a);
2383	prev=a;
2384	pIter(a);
2385	}
2386	else
2387	{
2388	if (prev==NULL)
2389	{
2390	pDelete1(&result);
2391	a=result;
2392	}
2393	else
2394	{
2395	pDelete1(&pNext(prev));
2396	a=pNext(prev);
2397	}
2398	}
2399	}
2400	pMultN(result,inv);
2401	nDelete(&inv);
2402	pDelete(&b);
2403	return result;
2404	}
2405
2406	/*2
2407	* returns the LCM of the head terms of a and b in *m
2408	*/
2409	void pLcm(poly a, poly b, poly m)
2410	{
2411	int i;
2412	for (i=pVariables; i; i--)
2413	{
2414	pSetExp(m,i, max( pGetExp(a,i), pGetExp(b,i)));
2415	}
2416	pSetComp(m, max(pGetComp(a), pGetComp(b)));
2417	}
2418
2419	/*2
2420	* convert monomial given as string to poly, e.g. 1x3y5z
2421	*/
2422	poly pmInit(char *st, BOOLEAN &ok)
2423	{
2424	int i,j;
2425	ok=FALSE;
2426	BOOLEAN b=FALSE;
2427	poly rc = pInit();
2428	char *s = nRead(st,&(rc->coef));
2429	if (s==st)
2430	/* i.e. it does not start with a coeff: test if it is a ringvar*/
2431	{
2432	j = rIsRingVar(s);
2433	if (j >= 0)
2434	{
2435	pIncrExp(rc,1+j);
2436	goto done;
2437	}
2438	}
2439	else
2440	b=TRUE;
2441	while (*s!='\0')
2442	{
2443	char ss[2];
2444	ss[0] = *s++;
2445	ss[1] = '\0';
2446	j = rIsRingVar(ss);
2447	if (j >= 0)
2448	{
2449	s = eati(s,&i);
2450	pAddExp(rc,1+j, (Exponent_t)i);
2451	}
2452	else
2453	{
2454	if ((s!=st)&&isdigit(st[0]))
2455	{
2456	errorreported=TRUE;
2457	}
2458	pDelete(&rc);
2459	return NULL;
2460	}
2461	}
2462	done:
2463	ok=!errorreported;
2464	if (nIsZero(pGetCoeff(rc))) pDelete1(&rc);
2465	else
2466	{
2467	#ifdef DRING
2468	if (pDRING)
2469	{
2470	for(i=1;i<=pdN;i++)
2471	{
2472	if(pGetExp(rc,pdDX(i))>0)
2473	{
2474	pdDFlag(rc)=1;
2475	break;
2476	}
2477	}
2478	}
2479	#endif
2480	pSetm(rc);
2481	}
2482	return rc;
2483	}
2484
2485	/*2
2486	*make p homgeneous by multiplying the monomials by powers of x_varnum
2487	*/
2488	poly pHomogen (poly p, int varnum)
2489	{
2490	poly q=NULL;
2491	poly res;
2492	int o,ii;
2493
2494	if (p!=NULL)
2495	{
2496	if ((varnum < 1) \|\| (varnum > pVariables))
2497	{
2498	return NULL;
2499	}
2500	o=pWTotaldegree(p);
2501	q=pNext(p);
2502	while (q != NULL)
2503	{
2504	ii=pWTotaldegree(q);
2505	if (ii>o) o=ii;
2506	pIter(q);
2507	}
2508	q = pCopy(p);
2509	res = q;
2510	while (q != NULL)
2511	{
2512	ii = o-pWTotaldegree(q);
2513	if (ii!=0)
2514	{
2515	pAddExp(q,varnum, (Exponent_t)ii);
2516	pSetm(q);
2517	}
2518	pIter(q);
2519	}
2520	q = pOrdPolyInsertSetm(res);
2521	}
2522	return q;
2523	}
2524
2525
2526	/*2
2527	*replaces the maximal powers of the leading monomial of p2 in p1 by
2528	*the same powers of n, utility for dehomogenization
2529	*/
2530	poly pDehomogen (poly p1,poly p2,number n)
2531	{
2532	polyset P;
2533	int SizeOfSet=5;
2534	int i;
2535	poly p;
2536	number nn;
2537
2538	P = (polyset)Alloc(5*sizeof(poly));
2539	for (i=0; i<5; i++)
2540	{
2541	P[i] = NULL;
2542	}
2543	pCancelPolyByMonom(p1,p2,&P,&SizeOfSet);
2544	p = P[0];
2545	//P[0] = NULL ;// for safety, may be remoeved later
2546	for (i=1; i<SizeOfSet; i++)
2547	{
2548	if (P[i] != NULL)
2549	{
2550	nPower(n,i,&nn);
2551	pMultN(P[i],nn);
2552	p = pAdd(p,P[i]);
2553	//P[i] =NULL; // for safety, may be removed later
2554	nDelete(&nn);
2555	}
2556	}
2557	Free((ADDRESS)P,SizeOfSet*sizeof(poly));
2558	return p;
2559	}
2560
2561	/*4
2562	*Returns the exponent of the maximal power of the leading monomial of
2563	*p2 in that of p1
2564	*/
2565	static int pGetMaxPower (poly p1,poly p2)
2566	{
2567	int i,k,res = 32000; /a very large integer/
2568
2569	if (p1 == NULL) return 0;
2570	for (i=1; i<=pVariables; i++)
2571	{
2572	if ( pGetExp(p2,i) != 0)
2573	{
2574	k = pGetExp(p1,i) / pGetExp(p2,i);
2575	if (k < res) res = k;
2576	}
2577	}
2578	return res;
2579	}
2580
2581	/*2
2582	*Returns as i-th entry of P the coefficient of the (i-1) power of
2583	*the leading monomial of p2 in p1
2584	*/
2585	void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet)
2586	{
2587	int maxPow;
2588	poly p,qp,Coeff;
2589
2590	if (*P == NULL)
2591	{
2592	P = (polyset) Alloc(5sizeof(poly));
2593	*SizeOfSet = 5;
2594	}
2595	p = pCopy(p1);
2596	while (p != NULL)
2597	{
2598	qp = p->next;
2599	p->next = NULL;
2600	maxPow = pGetMaxPower(p,p2);
2601	Coeff = pDivByMonom(p,p2);
2602	if (maxPow > *SizeOfSet)
2603	{
2604	pEnlargeSet(P,SizeOfSet,maxPow+1-SizeOfSet);
2605	*SizeOfSet = maxPow+1;
2606	}
2607	(P)[maxPow] = pAdd((P)[maxPow],Coeff);
2608	pDelete(&p);
2609	p = qp;
2610	}
2611	}
2612
2613	/*2
2614	*returns the leading monomial of p1 divided by the maximal power of that
2615	*of p2
2616	*/
2617	poly pDivByMonom (poly p1,poly p2)
2618	{
2619	int k, i;
2620
2621	if (p1 == NULL) return NULL;
2622	k = pGetMaxPower(p1,p2);
2623	if (k == 0)
2624	return pHead(p1);
2625	else
2626	{
2627	number n;
2628	poly p = pInit();
2629
2630	p->next = NULL;
2631	for (i=1; i<=pVariables; i++)
2632	{
2633	pSetExp(p,i, pGetExp(p1,i)-k* pGetExp(p2,i));
2634	}
2635	nPower(p2->coef,k,&n);
2636	pSetCoeff0(p,nDiv(p1->coef,n));
2637	nDelete(&n);
2638	pSetm(p);
2639	return p;
2640	}
2641	}
2642	/----------utilities for syzygies--------------/
2643	poly pTakeOutComp(poly * p, int k)
2644	{
2645	poly q = *p,qq=NULL,result = NULL;
2646
2647	if (q==NULL) return NULL;
2648	if (pGetComp(q)==k)
2649	{
2650	result = q;
2651	while ((q!=NULL) && (pGetComp(q)==k))
2652	{
2653	pSetComp(q,0);
2654	qq = q;
2655	pIter(q);
2656	}
2657	*p = q;
2658	pNext(qq) = NULL;
2659	}
2660	if (q==NULL) return result;
2661	if (pGetComp(q) > k) pDecrComp(q);
2662	poly pNext_q;
2663	while ((pNext_q=pNext(q))!=NULL)
2664	{
2665	if (pGetComp(pNext_q)==k)
2666	{
2667	if (result==NULL)
2668	{
2669	result = pNext_q;
2670	qq = result;
2671	}
2672	else
2673	{
2674	pNext(qq) = pNext_q;
2675	pIter(qq);
2676	}
2677	pNext(q) = pNext(pNext_q);
2678	pNext(qq) =NULL;
2679	pSetComp(qq,0);
2680	}
2681	else
2682	{
2683	/pIter(q);/ q=pNext_q;
2684	if (pGetComp(q) > k) pDecrComp(q);
2685	}
2686	}
2687	return result;
2688	}
2689
2690	poly pTakeOutComp1(poly * p, int k)
2691	{
2692	poly q = *p;
2693
2694	if (q==NULL) return NULL;
2695
2696	poly qq=NULL,result = NULL;
2697
2698	if (pGetComp(q)==k)
2699	{
2700	result = q; /* p /
2701	while ((q!=NULL) && (pGetComp(q)==k))
2702	{
2703	pSetComp(q,0);
2704	qq = q;
2705	pIter(q);
2706	}
2707	*p = q;
2708	pNext(qq) = NULL;
2709	}
2710	if (q==NULL) return result;
2711	// if (pGetComp(q) > k) pGetComp(q)--;
2712	while (pNext(q)!=NULL)
2713	{
2714	if (pGetComp(pNext(q))==k)
2715	{
2716	if (result==NULL)
2717	{
2718	result = pNext(q);
2719	qq = result;
2720	}
2721	else
2722	{
2723	pNext(qq) = pNext(q);
2724	pIter(qq);
2725	}
2726	pNext(q) = pNext(pNext(q));
2727	pNext(qq) =NULL;
2728	pSetComp(qq,0);
2729	}
2730	else
2731	{
2732	pIter(q);
2733	// if (pGetComp(q) > k) pGetComp(q)--;
2734	}
2735	}
2736	return result;
2737	}
2738
2739	void pDeleteComp(poly * p,int k)
2740	{
2741	poly q;
2742
2743	while ((p!=NULL) && (pGetComp(p)==k)) pDelete1(p);
2744	if (*p==NULL) return;
2745	q = *p;
2746	if (pGetComp(q)>k) pDecrComp(q);
2747	while (pNext(q)!=NULL)
2748	{
2749	if (pGetComp(pNext(q))==k)
2750	pDelete1(&(pNext(q)));
2751	else
2752	{
2753	pIter(q);
2754	if (pGetComp(q)>k) pDecrComp(q);
2755	}
2756	}
2757	}
2758	/----------end of utilities for syzygies--------------/
2759
2760	/*2
2761	* pair has no common factor ? or is no polynomial
2762	*/
2763	#ifdef COMP_FAST
2764	BOOLEAN pHasNotCF(poly p1, poly p2)
2765	{
2766	#ifdef SRING
2767	if (pSRING)
2768	return FALSE;
2769	#endif
2770
2771	if (pGetComp(p1) > 0 \|\| pGetComp(p2) > 0)
2772	return FALSE;
2773	Exponent_t * m1 = &(p1->exp[pVarLowIndex]);
2774	Exponent_t * m2 = &(p2->exp[pVarLowIndex]);
2775	int i = 1;
2776	loop
2777	{
2778	if (((m1) > 0) && ((m2) > 0))
2779	return FALSE;
2780	if (i == pVariables)
2781	return TRUE;
2782	m1++;m2++;
2783	i++;
2784	}
2785	}
2786	#else
2787	BOOLEAN pHasNotCF(poly p1, poly p2)
2788	{
2789	#ifdef SRING
2790	if (pSRING)
2791	return FALSE;
2792	#endif
2793	short * m1 = p1->exp;
2794	short * m2 = p2->exp;
2795
2796	if (((m1) > 0)\|\|((m2) > 0))
2797	return FALSE;
2798	int i = 1;
2799	loop
2800	{
2801	m1++;m2++;
2802	if (((m1) > 0) && ((m2) > 0))
2803	return FALSE;
2804	if (i == pVariables)
2805	return TRUE;
2806	i++;
2807	}
2808	}
2809	#endif
2810	/*
2811	*void pSFactors(poly f, poly g, poly a, poly b)
2812	*{
2813	* int i,d;
2814	*
2815	* for (i=pVariables;i>0;i--)
2816	* {
2817	* d = pGetExp(f,i)- pGetExp(g,i);
2818	* if (d >= 0)
2819	* {
2820	* pGetExp(a,i) = 0;
2821	* pGetExp(b,i) = d;
2822	* }
2823	* else
2824	* {
2825	* pGetExp(a,i) = -d;
2826	* pGetExp(b,i) = 0;
2827	* }
2828	* }
2829	* pGetComp(a) = 0;
2830	* pGetComp(b) = 0;
2831	* pSetm(a);
2832	* pSetm(b);
2833	*}
2834	*/
2835
2836	/*
2837	*void pSDiv(poly f, poly g, poly b)
2838	*{
2839	* int i,d;
2840	*
2841	* for (i=pVariables;i>0;i--)
2842	* {
2843	* d = pGetExp(f,i)- pGetExp(g,i);
2844	* pGetExp(b,i) = d;
2845	* }
2846	* pGetComp(b) = 0;
2847	* pSetm(b);
2848	*}
2849	*/
2850
2851	/*2
2852	* update the initial term of a polynomial a by multipying it by
2853	* the (number) coefficient
2854	* and the exponent vector (of) exp (a well initialized polynomial)
2855	*/
2856	/*
2857	*void pSMultBy(poly f, poly m)
2858	*{
2859	* number t;
2860	* int i;
2861	* // short notok;
2862	*
2863	* t=nMult(f->coef, m->coef);
2864	* nDelete(&(f->coef));
2865	* f->coef = t;
2866	* f->Order += m->Order;
2867	* for (i=pVariables; i; i--)
2868	* pGetExp(f,i) += pGetExp(m,i);
2869	* // if (notok)
2870	* // {
2871	* if (!( pGetComp(f)))
2872	* {
2873	* pGetComp(f) = pGetComp(m);
2874	* }
2875	* // else
2876	* // {
2877	* // HALT;
2878	* // }
2879	* // }
2880	*}
2881	*/
2882
2883
2884	/*2
2885	*should return 1 if p divides q and p<q,
2886	* -1 if q divides p and q<p
2887	* 0 otherwise
2888	*/
2889	#ifdef COMP_FAST
2890	int pDivComp(poly p, poly q)
2891	{
2892	if (pGetComp(p) == pGetComp(q))
2893	{
2894	Exponent_t * mp = &(p->exp[pVarLowIndex]);
2895	Exponent_t * mq = &(q->exp[pVarLowIndex]);
2896	int i=pVariables;
2897	BOOLEAN a=FALSE, b=FALSE;
2898	for (; i>0; i--)
2899	{
2900	if (mp<mq)
2901	{
2902	if (b) return 0;
2903	a =TRUE;
2904	}
2905	else if (mp>mq)
2906	{
2907	if (a) return 0;
2908	b = TRUE;
2909	}
2910	mp++;mq++;
2911	}
2912	if (a) return 1;
2913	else if (b) return -1;
2914	}
2915	return 0;
2916	}
2917	#else
2918	int pDivComp(poly p, poly q)
2919	{
2920	short * mp = p->exp;
2921	short * mq = q->exp;
2922
2923	if (mp==mq)
2924	{
2925	int i=pVariables;
2926	BOOLEAN a=FALSE, b=FALSE;
2927	for (; i>0; i--)
2928	{
2929	mp++;mq++;
2930	if (mp<mq)
2931	{
2932	if (b) return 0;
2933	a =TRUE;
2934	}
2935	else if (mp>mq)
2936	{
2937	if (a) return 0;
2938	b = TRUE;
2939	}
2940	}
2941	if (a) return 1;
2942	else if (b) return -1;
2943	}
2944	return 0;
2945	}
2946	#endif
2947	/*2
2948	*divides p1 by its leading monomial
2949	*/
2950	void pNorm(poly p1)
2951	{
2952	poly h;
2953	number k, c;
2954
2955	if (p1!=NULL)
2956	{
2957	if (!nIsOne(pGetCoeff(p1)))
2958	{
2959	nNormalize(pGetCoeff(p1));
2960	k=pGetCoeff(p1);
2961	c = nInit(1);
2962	pSetCoeff0(p1,c);
2963	h = pNext(p1);
2964	while (h!=NULL)
2965	{
2966	c=nDiv(pGetCoeff(h),k);
2967	if (!nIsOne(c)) nNormalize(c);
2968	pSetCoeff(h,c);
2969	pIter(h);
2970	}
2971	nDelete(&k);
2972	}
2973	else
2974	{
2975	h = pNext(p1);
2976	while (h!=NULL)
2977	{
2978	nNormalize(pGetCoeff(h));
2979	pIter(h);
2980	}
2981	}
2982	}
2983	}
2984
2985	/*2
2986	*normalize all coeffizients
2987	*/
2988	void pNormalize(poly p)
2989	{
2990	while (p!=NULL)
2991	{
2992	nTest(pGetCoeff(p));
2993	nNormalize(pGetCoeff(p));
2994	pIter(p);
2995	}
2996	}
2997
2998	/*3
2999	* substitute the n-th variable by 1 in p
3000	* destroy p
3001	*/
3002	static poly pSubst1 (poly p,int n)
3003	{
3004	poly qq,result = NULL;
3005
3006	while (p != NULL)
3007	{
3008	qq = p;
3009	pIter(p);
3010	qq->next = NULL;
3011	pSetExp(qq,n,0);
3012	pSetm(qq);
3013	result = pAdd(result,qq);
3014	}
3015	return result;
3016	}
3017
3018	/*2
3019	* substitute the n-th variable by e in p
3020	* destroy p
3021	*/
3022	poly pSubst(poly p, int n, poly e)
3023	{
3024	if ((e!=NULL)&&(pIsConstant(e))&&(nIsOne(pGetCoeff(e))))
3025	return pSubst1(p,n);
3026
3027	int exponent,i;
3028	poly h, res, m;
3029	Exponent_t me,ee;
3030	number nu,nu1;
3031
3032	me=(Exponent_t )Alloc((pVariables+1)sizeof(Exponent_t));
3033	ee=(Exponent_t )Alloc((pVariables+1)sizeof(Exponent_t));
3034	if (e!=NULL) pGetExpV(e,ee);
3035	res=NULL;
3036	h=p;
3037	while (h!=NULL)
3038	{
3039	if ((e!=NULL) \|\| (pGetExp(h,n)==0))
3040	{
3041	m=pHead(h);
3042	pGetExpV(m,me);
3043	exponent=me[n];
3044	me[n]=0;
3045	for(i=1;i<=pVariables;i++)
3046	me[i]+=exponent*ee[i];
3047	pSetExpV(m,me);
3048	if (e!=NULL)
3049	{
3050	nPower(pGetCoeff(e),exponent,&nu);
3051	nu1=nMult(pGetCoeff(m),nu);
3052	nDelete(&nu);
3053	pSetCoeff(m,nu1);
3054	}
3055	res=pAdd(res,m);
3056	}
3057	pDelete1(&h);
3058	}
3059	Free((ADDRESS)me,(pVariables+1)*sizeof(Exponent_t));
3060	Free((ADDRESS)ee,(pVariables+1)*sizeof(Exponent_t));
3061	return res;
3062	}
3063
3064	BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
3065	{
3066	int k, j;
3067
3068	if (lcm==NULL) return FALSE;
3069
3070	for (j=pVariables; j; j--)
3071	if ( pGetExp(p,j) > pGetExp(lcm,j)) return FALSE;
3072	if ( pGetComp(p) != pGetComp(lcm)) return FALSE;
3073	for (j=pVariables; j; j--)
3074	{
3075	if (pGetExp(p1,j)!=pGetExp(lcm,j))
3076	{
3077	if (pGetExp(p,j)!=pGetExp(lcm,j))
3078	{
3079	for (k=pVariables; k>j; k--)
3080	{
3081	if ((pGetExp(p,k)!=pGetExp(lcm,k))
3082	&& (pGetExp(p2,k)!=pGetExp(lcm,k)))
3083	return TRUE;
3084	}
3085	for (k=j-1; k; k--)
3086	{
3087	if ((pGetExp(p,k)!=pGetExp(lcm,k))
3088	&& (pGetExp(p2,k)!=pGetExp(lcm,k)))
3089	return TRUE;
3090	}
3091	return FALSE;
3092	}
3093	}
3094	else if (pGetExp(p2,j)!=pGetExp(lcm,j))
3095	{
3096	if (pGetExp(p,j)!=pGetExp(lcm,j))
3097	{
3098	for (k=pVariables; k>j; k--)
3099	{
3100	if ((pGetExp(p,k)!=pGetExp(lcm,k))
3101	&& (pGetExp(p1,k)!=pGetExp(lcm,k)))
3102	return TRUE;
3103	}
3104	for (k=j-1; k; k--)
3105	{
3106	if ((pGetExp(p,k)!=pGetExp(lcm,k))
3107	&& (pGetExp(p1,k)!=pGetExp(lcm,k)))
3108	return TRUE;
3109	}
3110	return FALSE;
3111	}
3112	}
3113	}
3114	return FALSE;
3115	}
3116
3117	int pWeight(int i)
3118	{
3119	if ((firstwv==NULL) \|\| (i>firstBlockEnds))
3120	{
3121	return 1;
3122	}
3123	return firstwv[i-1];
3124	}
3125
3126	#ifdef COMP_DEBUG
3127	static int debug_comp(poly p1, poly p2)
3128	{
3129	int t_d = t_pComp0(p1, p2);
3130	int f_d = f_pComp0(p1, p2);
3131
3132	if (t_d != f_d)
3133	{
3134	fprintf(stderr, "Error in comp1lpc\n");
3135	t_pComp0(p1, p2);
3136	f_pComp0(p1, p2);
3137	}
3138	return t_d;
3139	}
3140	#endif
3141
3142	#ifdef COMP_STATISTICS
3143	static int s_comp_lp_c_1(poly p1, poly p2)
3144	{
3145	MonomCountTotal++;
3146
3147	OrderCmp(p2, p1,
3148	{MonomCountOrderS++; return -1;},
3149	{MonomCountOrderG++; return 1;});
3150
3151	CompCmp(p2, p1,
3152	{MonomCountCompS++; return -pComponentOrder},
3153	{MonomCountCompG++; return pComponentOrder});
3154
3155	MonomCountEqual++;
3156	return 0;
3157	}
3158	#endif
3159
3160	#ifdef COMP_STATISTICS
3161	static int s_comp_lp_c_i(poly p1, poly p2)
3162	{
3163	MonomCountTotal++;
3164
3165	OrderCmp(p2, p1,
3166	{MonomCountOrderS++; return -1;},
3167	{MonomCountOrderG++; return 1;});
3168
3169	int i = SIZEOF_ORDER / SIZEOF_EXPONENT + 1;
3170	Exponent_t d;
3171
3172
3173	for (;;)
3174	{
3175	d = pGetExp(p1, i) - pGetExp(p2, i);
3176	if (d)
3177	{
3178	MonomCountExp[i]++;
3179	if (d < 0)
3180	{
3181	return -1;
3182	MonomCountExpS++;
3183	}
3184	MonomCountExpG++;
3185	return 1;
3186	}
3187	i++;
3188	if (i == pVariables) break;
3189	}
3190
3191	CompCmp(p2, p1,
3192	{MonomCountCompS++; return -pComponentOrder},
3193	{MonomCountCompG++; return pComponentOrder});
3194
3195	MonomCountEqual++;
3196	return 0;
3197	}
3198	#endif
3199
3200	#ifdef MONOM_COUNT
3201	static void InitMonomCount(unsigned int nvars)
3202	{
3203	if (gMonomCount == NULL)
3204	{
3205	gMonomCount = (unsigned long ) Alloc0(nvarssizeof(int));
3206	gMonomCountLength = nvars;
3207	}
3208	else if (gMonomCountLength < nvars)
3209	{
3210	int i;
3211	gMonomCount = (unsigned long *) ReAlloc(gMonomCount,
3212	gMonomCountLength*sizeof(int),
3213	nvars*sizeof(int));
3214	for (i = gMonomCountLength; i< nvars; i++)
3215	gMonomCount[i] = 0;
3216	gMonomCountLength = nvars;
3217	}
3218	}
3219
3220	void OutputMonomCount()
3221	{
3222	unsigned int i;
3223	unsigned long max = gMonomCount0;
3224	float fmax;
3225	unsigned long sum = 0;
3226
3227
3228	// check that everything went ok
3229	if (gMonomCountS !=
3230	gMonomCount0 + gMonomCountOrder + gMonomCountL + gMonomCountG)
3231	{
3232	printf("MonomCountTotalError: %10lu %10lu\n", gMonomCountS,
3233	gMonomCount0 + gMonomCountOrder + gMonomCountL + gMonomCountG);
3234	}
3235	for (i=0; i<gMonomCountLength; i++)
3236	sum += gMonomCount[i];
3237	if (sum != gMonomCountL + gMonomCountG)
3238	{
3239	printf("MonomCountNotEqualError: %10lu %10lu\n",
3240	gMonomCountL + gMonomCountG, sum);
3241	}
3242
3243	printf("Total : %10lu\n", gMonomCountS);
3244	if (gMonomCountS == 0) gMonomCountS = 1;
3245
3246	printf("Order : %10lu \t %.4f\n",gMonomCountOrder,
3247	(float) gMonomCountOrder / (float) gMonomCountS);
3248	printf("Equal : %10lu \t %.4f\n",gMonomCount0,
3249	(float) gMonomCount0 / (float) gMonomCountS);
3250	printf("NotEqual: %10lu \t %.4f\n",gMonomCountL + gMonomCountG,
3251	((float) gMonomCountL + gMonomCountG) / (float) gMonomCountS);
3252	printf("\tGreater: %10lu \t %.4f\n",gMonomCountG,
3253	(float) gMonomCountG / (float) gMonomCountS);
3254	printf("\tLess : %10lu \t %.4f\n",gMonomCountL,
3255	(float) gMonomCountL / (float) gMonomCountS);
3256
3257	for (i=0; i<gMonomCountLength; i++)
3258	{
3259	printf("\t\t[%d] : %10lu \t %.4f\n", i, gMonomCount[i],
3260	(float) gMonomCount[i] / (float) gMonomCountS);
3261	}
3262	printf("E1Neg : %10lu \t %.4f\n",gMonomCountN1,
3263	(float) gMonomCountN1 / (float) gMonomCountS);
3264	printf("E2Neg : %10lu \t %.4f\n",gMonomCountN2,
3265	(float) gMonomCountN2 / (float) gMonomCountS);
3266	printf("BothNeg : %10lu \t %.4f\n",gMonomCountNN,
3267	(float) gMonomCountNN / (float) gMonomCountS);
3268	printf("Mcount2 : %u\n", gMcount2);
3269	}
3270
3271	void ResetMonomCount()
3272	{
3273	int i;
3274
3275	for (i=0; i<gMonomCountLength; i++)
3276	gMonomCount[i] = 0;
3277	gMonomCountOrder = 0;
3278	gMonomCount0 = 0;
3279	gMonomCountS = 0;
3280	gMonomCountG = 0;
3281	gMonomCountL = 0;
3282	gMonomCountN1 = 0;
3283	gMonomCountN2 = 0;
3284	gMonomCountNN = 0;
3285	gMcount2 = 0;
3286	}
3287
3288	#endif
3289	#ifdef MONOM_COUNT
3290	unsigned long gMonomCount0 = 0;
3291	unsigned long gMonomCountLength = 0;
3292	unsigned long gMonomCountOrder = 0;
3293	unsigned long gMonomCountG = 0;
3294	unsigned long gMonomCountL = 0;
3295	unsigned long gMonomCountS = 0;
3296	unsigned long gMonomCountN1 = 0;
3297	unsigned long gMonomCountN2 = 0;
3298	unsigned long gMonomCountNN = 0;
3299	unsigned long gMcount2 = 0;
3300	unsigned long* gMonomCount = NULL;
3301	#endif

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

Download in other formats: