Context Navigation

source: git/Singular/polys.cc @ ca7a56

fieker-DuValspielwiese

Last change on this file since ca7a56 was eb17bd3, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: changes for TEST_MAC_ORDER (binom.cc binom.h polys-impl.cc polys-impl.h polys.cc spSpolyLoop.cc spolys0.h) git-svn-id: file:///usr/local/Singular/svn/trunk@1029 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 69.2 KB

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

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

Download in other formats: