source: git/Singular/polys.cc @ 94fcf4

fieker-DuValspielwiese
Last change on this file since 94fcf4 was 311499, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: DRING-changes git-svn-id: file:///usr/local/Singular/svn/trunk@2007 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 49.7 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: polys.cc,v 1.28 1998-05-28 16:50:50 Singular 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#include "polys-comp.h"
24
25/* ----------- global variables, set by pChangeRing --------------------- */
26/* initializes the internal data from the exp vector */
27pSetmProc pSetm;
28/* computes length and maximal degree of a POLYnomial */
29pLDegProc pLDeg;
30/* computes the degree of the initial term, used for std */
31pFDegProc pFDeg;
32/* the monomial ordering of the head monomials a and b */
33/* returns -1 if a comes before b, 0 if a=b, 1 otherwise */
34pCompProc pComp0;
35
36int pVariables;     // number of variables
37int pVariablesW;    // number of words of pVariables exponents
38int pVariables1W;   // number of words of (pVariables+1) exponents
39int pMonomSize;     // size of monom (in bytes)
40int pMonomSizeW;    // size of monom (in words)
41int pLexSgn;        // 1, for lex monom comps; -1 otherwise (exception: ls)
42int pVarOffset;     // controls the way exponents are stored in a vector
43int pVarLowIndex;   // lowest exponent index
44int pVarHighIndex;  // highest exponent index
45int pVarCompIndex;  // Location of component in exponent vector
46
47/* 1 for polynomial ring, -1 otherwise */
48int     pOrdSgn;
49/* TRUE for momomial output as x2y, FALSE for x^2*y */
50int pShortOut = (int)TRUE;
51// it is of type int, not BOOLEAN because it is also in ip
52/* TRUE if the monomial ordering is not compatible with pFDeg */
53BOOLEAN pLexOrder;
54/* TRUE if the monomial ordering has polynomial and power series blocks */
55BOOLEAN pMixedOrder;
56/* 1 for c ordering, -1 otherwise (i.e. for C ordering) */
57int  pComponentOrder;
58
59#ifdef DRING
60int      p2;
61BOOLEAN  pDRING=FALSE;
62#endif
63
64#ifdef SRING
65int      pAltVars;
66BOOLEAN  pSRING=FALSE;
67#endif
68
69#ifdef SDRING
70BOOLEAN  pSDRING=FALSE;
71#include "polys.inc"
72#endif
73
74/* ----------- global variables, set by procedures from hecke/kstd1 ----- */
75/* the highest monomial below pHEdge */
76poly      ppNoether = NULL;
77
78/* -------------- static variables --------------------------------------- */
79/*is the basic comparing procedure during a computation of syzygies*/
80static pCompProc pCompOld;
81/*for grouping module indicies during computations*/
82int pMaxBound = 0;
83
84/*contains the headterms for the Schreyer orderings*/
85static int* SchreyerOrd;
86static int maxSchreyer=0;
87static int indexShift=0;
88static pLDegProc pLDegOld;
89
90typedef int (*bcompProc)(poly p1, poly p2, int i1, int i2, short * w);
91static bcompProc bcomph[20];
92static short**   polys_wv;
93static short *   firstwv;
94static int * block0;
95static int * block1;
96static int firstBlockEnds;
97static int * order;
98
99/*0 implementation*/
100/*-------- the several possibilities for pSetm:-----------------------*/
101
102/* Remark: These could be made more efficient by avoiding using pGetExp */
103
104/*2
105* define the order of p with respect to lex. ordering, N=1
106*/
107static void setlex1(poly p)
108{
109  p->Order = (Order_t)pGetExp(p,1);
110}
111
112/*2
113* define the order of p with respect to lex. ordering, N>1
114*/
115static void setlex2(poly p)
116{
117  p->Order = (((Order_t)pGetExp(p,1))<<(sizeof(Exponent_t)*8))
118    + (Order_t)pGetExp(p,2);
119}
120
121/*2
122* define the order of p with respect to a degree ordering
123*/
124static void setdeg1(poly p)
125{
126  p->Order = pExpQuerSum1(p, firstBlockEnds);
127}
128
129/*2
130* define the order of p with respect to a degree ordering
131* with weigthts
132*/
133static void setdeg1w(poly p)
134{
135  Order_t i, j = 0;
136
137  for (i = firstBlockEnds; i>0; i--)
138  {
139     j += ((Order_t) pGetExp(p,i))*firstwv[i-1];
140  }
141  p->Order = j;
142}
143
144
145/*-------- IMPLEMENTATION OF MONOMIAL COMPARISONS ---------------------*/
146
147
148#define NonZeroR(l, actionG, actionS)           \
149do                                              \
150{                                               \
151  long _l = l;                                  \
152  if (_l)                                       \
153  {                                             \
154    if (_l > 0) actionG;                        \
155    actionS;                                    \
156  }                                             \
157}                                               \
158while(0)
159
160#define Mreturn(d, multiplier)                      \
161{                                                   \
162  if (d > 0) return multiplier;                     \
163  return -multiplier;                               \
164}
165
166static int pComp_otEXP_nwONE(poly p1, poly p2);
167static int pComp_otCOMPEXP_nwONE(poly p1, poly p2);
168static int pComp_otEXPCOMP_nwONE(poly p1, poly p2);
169static int pComp_otEXP_nwTWO(poly p1, poly p2);
170static int pComp_otCOMPEXP_nwTWO(poly p1, poly p2);
171static int pComp_otEXPCOMP_nwTWO(poly p1, poly p2);
172static int pComp_otEXP_nwEVEN(poly p1, poly p2);
173static int pComp_otCOMPEXP_nwEVEN(poly p1, poly p2);
174static int pComp_otEXPCOMP_nwEVEN(poly p1, poly p2);
175static int pComp_otEXP_nwODD(poly p1, poly p2);
176static int pComp_otCOMPEXP_nwODD(poly p1, poly p2);
177static int pComp_otEXPCOMP_nwODD(poly p1, poly p2);
178
179
180// comp_nwONE is used if pVariables1W == 1 and component is compatible
181// with ordering
182static int pComp_otEXP_nwONE(poly p1, poly p2)
183{
184  register long d = pGetOrder(p1) - pGetOrder(p2);
185
186  if (d) Mreturn(d, pOrdSgn);
187  _pMonComp_otEXP_nwONE(p1, p2, d, goto NotEqual, return 0);
188
189  NotEqual:
190  Mreturn(d, pLexSgn);
191}
192
193// comp_otEXPCOMP_nwONE :  pVariables1W == 1, priority is
194// given to exponents, component is incompatible with ordering
195static int pComp_otEXPCOMP_nwONE(poly p1, poly p2)
196{
197  register long d = pGetOrder(p1) - pGetOrder(p2);
198
199  if (d) Mreturn(d, pOrdSgn);
200  _pMonComp_otEXPCOMP_nwONE(p1, p2, d, goto NotEqual , return 0);
201
202  NotEqual:
203  Mreturn(d, pLexSgn)
204}
205
206// comp_otEXPCOMP_nwONE :  pVariables1W == 1, priority is given to component,
207// component is incompatible with ordering
208static int pComp_otCOMPEXP_nwONE(poly p1, poly p2)
209{
210  register long d = pGetComp(p2) - pGetComp(p1);
211  if (d) Mreturn(d, pComponentOrder);
212  d = pGetOrder(p1) - pGetOrder(p2);
213  if (d) Mreturn(d, pOrdSgn);
214  _pMonComp_otEXP_nwONE(p1, p2, d, goto NotEqual , return 0);
215
216  NotEqual:
217  Mreturn(d, pLexSgn);
218}
219
220// comp_nwTWO :  pVariables1W == 2 and component is compatible with ordering
221static int pComp_otEXP_nwTWO(poly p1, poly p2)
222{
223  register long d = pGetOrder(p1) - pGetOrder(p2);
224
225  if (d) Mreturn(d, pOrdSgn);
226  _pMonComp_otEXP_nwTWO(p1, p2, d, goto NotEqual , return 0);
227
228  NotEqual:
229  Mreturn(d, pLexSgn);
230}
231
232// comp_otEXPCOMP_nwTWO :  pVariables1W == 2, priority is given to exponents,
233// component is incompatible with ordering
234static int pComp_otEXPCOMP_nwTWO(poly p1, poly p2)
235{
236  register long d = pGetOrder(p1) - pGetOrder(p2);
237
238  if (d) Mreturn(d, pOrdSgn);
239  _pMonComp_otEXPCOMP_nwTWO(p1, p2, d, goto NotEqual, return 0);
240
241  NotEqual:
242  Mreturn(d, pLexSgn);
243}
244
245// comp_otEXPCOMP_nwTWO :  pVariables1W == 2, priority is given to component,
246// component is incompatible with ordering
247static int pComp_otCOMPEXP_nwTWO(poly p1, poly p2)
248{
249  register long d = pGetComp(p2) - pGetComp(p1);
250  if (d) Mreturn(d, pComponentOrder);
251  d = pGetOrder(p1) - pGetOrder(p2);
252  if (d) Mreturn(d, pOrdSgn);
253  _pMonComp_otEXP_nwTWO(p1, p2, d, goto NotEqual , return 0);
254
255  NotEqual:
256  Mreturn(d, pLexSgn);
257}
258
259// comp_nwEVEN :  pVariables1W == 2*i and component is compatible
260// with ordering
261static int pComp_otEXP_nwEVEN(poly p1, poly p2)
262{
263  register long d = pGetOrder(p1) - pGetOrder(p2);
264
265  if (d) Mreturn(d, pOrdSgn);
266  _pMonComp_otEXP_nwEVEN(p1, p2, pVariables1W, d, goto NotEqual , return 0);
267
268  NotEqual:
269  Mreturn(d, pLexSgn);
270}
271
272// comp_otEXPCOMP_nwEVEN : pVariables1W == 2*i, priority is given to exponents,
273// component is incompatible with ordering
274static int pComp_otEXPCOMP_nwEVEN(poly p1, poly p2)
275{
276  register long d = pGetOrder(p1) - pGetOrder(p2);
277
278  if (d) Mreturn(d, pOrdSgn);
279  _pMonComp_otEXPCOMP_nwEVEN(p1, p2, pVariables1W, d, goto NotEqual , return 0);
280
281  NotEqual:
282  Mreturn(d, pLexSgn);
283}
284
285// comp_otEXPCOMP_nwEVEN : pVariables1W == 2*i, priority is given to component,
286// component is incompatible with ordering
287static int pComp_otCOMPEXP_nwEVEN(poly p1, poly p2)
288{
289  register long d = pGetComp(p2) - pGetComp(p1);
290  if (d) Mreturn(d, pComponentOrder);
291  d = pGetOrder(p1) - pGetOrder(p2);
292  if (d) Mreturn(d, pOrdSgn);
293  _pMonComp_otEXP_nwEVEN(p1, p2, pVariablesW, d, goto NotEqual, return 0);
294
295  NotEqual:
296  Mreturn(d, pLexSgn);
297}
298
299// comp_nwODD : pVariables1W == 2*i and component is compatible
300// with ordering
301static int pComp_otEXP_nwODD(poly p1, poly p2)
302{
303  register long d = pGetOrder(p1) - pGetOrder(p2);
304
305  if (d) Mreturn(d, pOrdSgn);
306  _pMonComp_otEXP_nwODD(p1, p2, pVariables1W, d, goto NotEqual, return 0);
307
308  NotEqual:
309  Mreturn(d, pLexSgn);
310}
311
312// comp_otEXPCOMP_nwODD : pVariables1W == 2*i, priority is given to exponents,
313// component is incompatible with ordering
314static int pComp_otEXPCOMP_nwODD(poly p1, poly p2)
315{
316  register long d = pGetOrder(p1) - pGetOrder(p2);
317  if (d)
318  {
319    if (d > 0) return pOrdSgn;
320    return -pOrdSgn;
321  }
322
323  _pMonComp_otEXPCOMP_nwODD(p1, p2, pVariables1W, d, goto NotEqual , return 0);
324
325  NotEqual:
326  Mreturn(d, pLexSgn);
327}
328
329// comp_otCOMPEXP_nwODD : pVariables1W == 2*i, priority is given to component,
330// component is incompatible with ordering
331static int pComp_otCOMPEXP_nwODD(poly p1, poly p2)
332{
333  register long d = pGetComp(p2) - pGetComp(p1);
334  if (d) Mreturn(d, pComponentOrder);
335  d = pGetOrder(p1) - pGetOrder(p2);
336  if (d) Mreturn(d, pOrdSgn);
337  _pMonComp_otEXP_nwODD(p1, p2, pVariablesW, d, goto NotEqual, return 0);
338
339  NotEqual:
340  Mreturn(d, pLexSgn);
341}
342/*2
343* compare the head monomial of p1 and p2 with weight vector
344*/
345static int comp1a  ( poly p1, poly p2, int f, int l, short * w )
346{
347  int d= pGetOrder(p1) - pGetOrder(p2);
348  if ( d > 0 /*p1->Order > p2->Order*/ )
349    return 1;
350  else if ( d < 0 /*p1->Order < p2->Order*/ )
351    return -1;
352  return 0;
353}
354
355/*---------------------------------------------------*/
356
357/* These functions could be made faster if you use pointers to the
358* exponent vectors and pointer arithmetic instead of using the
359* macro pGetExp !!!
360*/
361
362/*2
363* compare the head monomial of p1 and p2 with lexicographic ordering
364*/
365static int comp_lp ( poly p1, poly p2, int f, int l, short * w )
366{
367  int i = f;
368  while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
369    i++;
370  if ( i > l )
371    return 0;
372  if ( pGetExp(p1,i) > pGetExp(p2,i) )
373    return 1;
374  return -1;
375}
376
377/*2
378* compare the head monomial of p1 and p2 with degree reverse lexicographic
379* ordering
380*/
381static int comp_dp ( poly p1, poly p2, int f, int l, short * w )
382{
383  int i, s1 = 0, s2 = 0;
384
385  for ( i = f; i <= l; i++ )
386  {
387    s1 += pGetExp(p1,i);
388    s2 += pGetExp(p2,i);
389  }
390  if ( s1 == s2 )
391  {
392    i = l;
393    while ( (i >= f ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
394      i--;
395    if ( i < f )
396      return 0;
397    if ( pGetExp(p1,i) > pGetExp(p2,i) )
398      return -1;
399    return 1;
400  }
401  if ( s1 > s2 )
402    return 1;
403  return -1;
404}
405
406/*2
407* compare the head monomial of p1 and p2 with degree lexicographic ordering
408*/
409static int comp_Dp ( poly p1, poly p2, int f, int l, short * w )
410{
411  int i, s1 = 0, s2 = 0;
412
413  for ( i = f; i <= l; i++ )
414  {
415    s1 += pGetExp(p1,i);
416    s2 += pGetExp(p2,i);
417  }
418  if ( s1 == s2 )
419  {
420    i = f;
421    while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
422      i++;
423    if ( i > l )
424      return 0;
425    if ( pGetExp(p1,i) > pGetExp(p2,i) )
426      return 1;
427    return -1;
428  }
429  if ( s1 > s2 )
430    return 1;
431  return -1;
432}
433
434/*2
435* compare the head monomial of p1 and p2 with weighted degree reverse
436* lexicographic ordering
437*/
438static int comp_wp ( poly p1, poly p2, int f, int l, short * w )
439{
440  int i, s1 = 0, s2 = 0;
441
442  for ( i = f; i <= l; i++, w++ )
443  {
444    s1 += (int)pGetExp(p1,i)*(*w);
445    s2 += (int)pGetExp(p2,i)*(*w);
446  }
447  if ( s1 == s2 )
448  {
449    i = l;
450    while ( ( i >= f ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
451      i--;
452    if ( i < f )
453      return 0;
454    if ( pGetExp(p1,i) > pGetExp(p2,i) )
455      return -1;
456    return 1;
457  }
458  if ( s1 > s2 )
459    return 1;
460  return -1;
461}
462
463/*2
464* compare the head monomial of p1 and p2 with weighted degree lexicographic
465* ordering
466*/
467static int comp_Wp ( poly p1, poly p2, int f, int l, short * w )
468{
469  int i, s1 = 0, s2 = 0;
470
471  for ( i = f; i <= l; i++, w++ )
472  {
473    s1 += (int)pGetExp(p1,i)*(*w);
474    s2 += (int)pGetExp(p2,i)*(*w);
475  }
476  if ( s1 == s2 )
477  {
478    i = f;
479    while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
480      i++;
481    if ( i > l )
482      return 0;
483    if ( pGetExp(p1,i) > pGetExp(p2,i) )
484      return 1;
485    return -1;
486  }
487  if ( s1 > s2 )
488    return 1;
489  return -1;
490}
491
492/*2
493* compare the head monomial of p1 and p2 with lexicographic ordering
494* (power series case)
495*/
496static int comp_ls ( poly p1, poly p2, int f, int l, short * w )
497{
498  int i;
499
500  i = f;
501  while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
502    i++;
503  if ( i > l )
504    return 0;
505  if ( pGetExp(p1,i) < pGetExp(p2,i) )
506    return 1;
507  return -1;
508}
509
510/*2
511* compare the head monomial of p1 and p2 with degree reverse lexicographic
512* ordering (power series case)
513*/
514static int comp_ds ( poly p1, poly p2, int f, int l, short * w )
515{
516  int i, s1 = 0, s2 = 0;
517
518  for ( i = f; i <= l; i++ )
519  {
520    s1 += pGetExp(p1,i);
521    s2 += pGetExp(p2,i);
522  }
523  if ( s1 == s2 )
524  {
525    i = l;
526    while ( ( i >= f ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
527      i--;
528    if ( i < f )
529      return 0;
530    if ( pGetExp(p1,i) < pGetExp(p2,i) )
531      return 1;
532    return -1;
533  }
534  if ( s1 < s2 )
535    return 1;
536  return -1;
537}
538
539/*2
540* compare the head monomial of p1 and p2 with degree lexicographic ordering
541* (power series case)
542*/
543static int comp_Ds ( 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 = f;
555    while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
556      i++;
557    if ( i > l )
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 weighted degree reverse
570* lexicographic ordering (power series case)
571*/
572static int comp_ws ( poly p1, poly p2, int f, int l, short * w )
573{
574  int i, s1 = 0, s2 = 0;
575
576  for ( i = f; i <= l; i++, w++ )
577  {
578    s1 += (int)pGetExp(p1,i)*(*w);
579    s2 += (int)pGetExp(p2,i)*(*w);
580  }
581  if ( s1 == s2 )
582  {
583    i = l;
584    while ( ( i >= f ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
585      i--;
586    if ( i < f )
587      return 0;
588    if ( pGetExp(p1,i) < pGetExp(p2,i) )
589      return 1;
590    return -1;
591  }
592  if ( s1 < s2 )
593    return 1;
594  return -1;
595}
596
597/*2
598* compare the head monomial of p1 and p2 with weighted degree lexicographic
599* ordering (power series case)
600*/
601static int comp_Ws ( poly p1, poly p2, int f, int l, short * w )
602{
603  int i, s1 = 0, s2 = 0;
604
605  for ( i = f; i <= l; i++, w++ )
606  {
607    s1 += (int)pGetExp(p1,i)*(*w);
608    s2 += (int)pGetExp(p2,i)*(*w);
609  }
610  if ( s1 == s2 )
611  {
612    i = f;
613    while ( ( i <= l ) && ( pGetExp(p1,i) == pGetExp(p2,i) ) )
614      i++;
615    if ( i > l )
616      return 0;
617    if ( pGetExp(p1,i) < pGetExp(p2,i) )
618      return -1;
619    return 1;
620  }
621  if ( s1 < s2 )
622    return 1;
623  return -1;
624}
625
626/*2
627* compare the head monomial of p1 and p2 with matrix order
628* w contains a series of l-f+1 lines
629*/
630static int comp_M ( poly p1, poly p2, int f, int l, short * w )
631{
632  int i, j, s1, s2;
633
634  for ( i = f; i <= l; i++ )
635  {
636    s1 = s2 = 0;
637    for ( j = f; j <= l; j++, w++ )
638    {
639      s1 += (int)pGetExp(p1,j)*(int)(*w);
640      s2 += (int)pGetExp(p2,j)*(int)(*w);
641    }
642    if ( s1 < s2 )
643      return -1;
644    if ( s1 > s2 )
645      return 1;
646    /* now w points to the last element of the current row, the next w++ */
647    /* moves on to the first element of the next row ! */
648  }
649  return 0;
650}
651
652/*2
653* compare the head monomial of p1 and p2 with weight vector
654*/
655static int comp_a  ( poly p1, poly p2, int f, int l, short * w )
656{
657  int i, s1 = 0, s2 = 0;
658
659  for ( i = f; i <= l; i++, w++ )
660  {
661    s1 += (int)pGetExp(p1,i)*(*w);
662    s2 += (int)pGetExp(p2,i)*(*w);
663  }
664  if ( s1 > s2 )
665    return 1;
666  if ( s1 < s2 )
667    return -1;
668  return 0;
669}
670
671/*2
672* compare the head monomial of p1 and p2 with module component
673*/
674static int comp_c  ( poly p1, poly p2, int f, int l, short * w )
675{
676  if ( pGetComp(p1) > pGetComp(p2) )
677    return -pComponentOrder;
678  if ( pGetComp(p1) < pGetComp(p2) )
679    return pComponentOrder;
680  return 0;
681}
682
683/*---------------------------------------------------------------*/
684
685/*2
686* compare p1 and p2 by a block ordering
687* uses (*bcomph[])() to do the real work
688*/
689static int BlockComp(poly p1, poly p2)
690{
691  int res, i, e, a;
692
693  /*4 compare in all blocks,*
694  * each block has var numbers a(=block0[i]) to e (=block1[i])*
695  * the block number starts with 0*/
696  e = 0;
697  i = 0;
698  loop
699  {
700    a = block0[i];
701    e = block1[i];
702    res = (*bcomph[i])(p1, p2, a, e , polys_wv[i]);
703    if (res)
704      return res;
705    i++;
706    if (order[i]==0)
707      break;
708  }
709  return 0;
710}
711
712int pComp(poly p1, poly p2)
713{
714  if (p2==NULL)
715    return 1;
716  if (p1==NULL)
717    return -1;
718  return pComp0(p1,p2);
719}
720
721/*----------pComp handling for syzygies---------------------*/
722static void newHeadsB(polyset actHeads,int length)
723{
724  int i;
725  int* newOrder=(int*)Alloc(length*sizeof(int));
726
727  for (i=0;i<length;i++)
728  {
729    if (actHeads[i]!=NULL)
730    {
731      newOrder[i] = SchreyerOrd[pGetComp(actHeads[i])-1];
732    }
733    else
734    {
735      newOrder[i]=0;
736    }
737  }
738  Free((ADDRESS)SchreyerOrd,maxSchreyer*sizeof(int));
739  SchreyerOrd = newOrder;
740  maxSchreyer = length;
741/*
742*for (i=0;i<maxSchreyer;i++); Print("%d  ",SchreyerOrd[i]);
743*PrintLn();
744*/
745}
746
747int mcompSchrB(poly p1,poly p2)
748{
749  int CompP1=pGetComp(p1),CompP2=pGetComp(p2),result,
750      cP1=SchreyerOrd[CompP1-1],cP2=SchreyerOrd[CompP2-1];
751
752  if (CompP1==CompP2) return pCompOld(p1,p2);
753  pSetComp(p1,cP1);
754  pSetComp(p2,cP2);
755  result = pCompOld(p1,p2);
756  pSetComp(p1,CompP1);
757  pSetComp(p2,CompP2);
758  if (!result)
759  {
760    if (CompP1>CompP2)
761      return -1;
762    else if (CompP1<CompP2)
763      return 1;
764  }
765  return result;
766}
767
768
769static void newHeadsM(polyset actHeads,int length)
770{
771  int i;
772  int* newOrder=
773    (int*)Alloc0((length+maxSchreyer-indexShift)*sizeof(int));
774
775  //for (i=0;i<length+maxSchreyer-indexShift;i++)
776  //  newOrder[i]=0;
777  for (i=indexShift;i<maxSchreyer;i++)
778  {
779    newOrder[i-indexShift] = SchreyerOrd[i];
780    SchreyerOrd[i] = 0;
781  }
782  for (i=maxSchreyer-indexShift;i<length+maxSchreyer-indexShift;i++)
783    newOrder[i] = newOrder[pGetComp(actHeads[i-maxSchreyer+indexShift])-1];
784  Free((ADDRESS)SchreyerOrd,maxSchreyer*sizeof(int));
785  SchreyerOrd = newOrder;
786  indexShift = maxSchreyer-indexShift;
787  maxSchreyer = length+indexShift;
788}
789
790/*2
791* compute the length of a polynomial (in l)
792* and the degree of the monomial with maximal degree:
793* this is NOT the last one and the module component
794* has to be <= indexShift
795*/
796static int ldegSchrM(poly p,int *l)
797{
798  int  t,max;
799
800  (*l)=1;
801  max=pFDeg(p);
802  while ((pNext(p)!=NULL) && (pGetComp(pNext(p))<=indexShift))
803  {
804    pIter(p);
805    t=pFDeg(p);
806    if (t>max) max=t;
807    (*l)++;
808  }
809  return max;
810}
811
812int mcompSchrM(poly p1,poly p2)
813{
814  if ( pGetComp(p1)<=indexShift)
815  {
816    if ( pGetComp(p2)>indexShift) return 1;
817  }
818  else if ( pGetComp(p2)<=indexShift)  return -1;
819  return mcompSchrB(p1,p2);
820}
821
822void pSetSchreyerOrdM(polyset nextOrder, int length,int comps)
823{
824  int i;
825
826  if (length!=0)
827  {
828    if (maxSchreyer!=0)
829      newHeadsM(nextOrder, length);
830    else
831    {
832      indexShift = comps;
833      if (indexShift==0) indexShift = 1;
834      SchreyerOrd = (int*)Alloc((indexShift+length)*sizeof(int));
835      maxSchreyer = length+indexShift;
836      for (i=0;i<indexShift;i++)
837        SchreyerOrd[i] = i;
838      for (i=indexShift;i<maxSchreyer;i++)
839        SchreyerOrd[i] = pGetComp(nextOrder[i-indexShift]);
840      pCompOld = pComp0;
841      pComp0 = mcompSchrM;
842      pLDegOld = pLDeg;
843      pLDeg = ldegSchrM;
844    }
845  }
846  else
847  {
848    if (maxSchreyer!=0)
849    {
850      Free((ADDRESS)SchreyerOrd,maxSchreyer*sizeof(int));
851      maxSchreyer = 0;
852      indexShift = 0;
853      pComp0 = pCompOld;
854      pLDeg = pLDegOld;
855    }
856  }
857}
858
859/*2
860*the pComp0 for normal syzygies
861*compares monomials in the usual ring order (pCompOld)
862*but groups module indecees according indexBounds befor
863*/
864static int mcompSyz(poly p1,poly p2)
865{
866  if (pGetComp(p1)<=pMaxBound)
867  {
868    if (pGetComp(p2)>pMaxBound) return 1;
869  }
870  else if (pGetComp(p2)<=pMaxBound)
871  {
872    return -1;
873  }
874  return pCompOld(p1,p2);
875}
876
877void pSetSyzComp(int k)
878{
879  if (k!=0)
880  {
881    if (pMaxBound==0)
882    {
883      pCompOld = pComp0;
884      pComp0 = mcompSyz;
885    }
886    pMaxBound = k;
887  }
888  else
889  {
890    if (pMaxBound!=0)
891    {
892      pComp0 = pCompOld;
893      pMaxBound = 0;
894    }
895  }
896}
897
898/*2
899* the type of the module ordering: C: -1, c: 1
900*/
901int pModuleOrder()
902{
903  return pComponentOrder;
904}
905
906/* -------------------------------------------------------------------*/
907/* several possibilities for pFDeg: the degree of the head term       */
908/*2
909* compute the degree of the leading monomial of p
910* the ordering is compatible with degree, use a->order
911*/
912int pDeg(poly a)
913{
914  return ((a!=NULL) ? (a->Order) : (-1));
915}
916
917/*2
918* compute the degree of the leading monomial of p
919* with respect to weigths 1
920* (all are 1 so save multiplications or they are of different signs)
921* the ordering is not compatible with degree so do not use p->Order
922*/
923int pTotaldegree(poly p)
924{
925  return pExpQuerSum(p);
926}
927
928/*2
929* compute the degree of the leading monomial of p
930* with respect to weigths from the ordering
931* the ordering is not compatible with degree so do not use p->Order
932*/
933int pWTotaldegree(poly p)
934{
935  int i, k;
936  int j =0;
937
938  // iterate through each block:
939  for (i=0;order[i]!=0;i++)
940  {
941    switch(order[i])
942    {
943      case ringorder_wp:
944      case ringorder_ws:
945      case ringorder_Wp:
946      case ringorder_Ws:
947        for (k=block0[i];k<=block1[i];k++)
948        { // in jedem block:
949          j+= pGetExp(p,k)*polys_wv[i][k-block0[i]];
950        }
951        break;
952      case ringorder_M:
953      case ringorder_lp:
954      case ringorder_dp:
955      case ringorder_ds:
956      case ringorder_Dp:
957      case ringorder_Ds:
958        for (k=block0[i];k<=block1[i];k++)
959        {
960          j+= pGetExp(p,k);
961        }
962        break;
963      case ringorder_c:
964      case ringorder_C:
965        break;
966      case ringorder_a:
967        for (k=block0[i];k<=block1[i];k++)
968        { // only one line
969          j+= pGetExp(p,k)*polys_wv[i][k-block0[i]];
970        }
971        return j;
972    }
973  }
974  return  j;
975}
976/* ---------------------------------------------------------------------*/
977/* several possibilities for pLDeg: the maximal degree of a monomial in p*/
978/*  compute in l also the pLength of p                                   */
979
980/*2
981* compute the length of a polynomial (in l)
982* and the degree of the monomial with maximal degree: the last one
983*/
984static int ldeg0(poly p,int *l)
985{
986  Exponent_t k= pGetComp(p);
987  int ll=1;
988
989  while ((pNext(p)!=NULL) && (pGetComp(pNext(p))==k))
990  {
991    pIter(p);
992    ll++;
993  }
994  *l=ll;
995  return (p->Order);
996}
997
998/*2
999* compute the length of a polynomial (in l)
1000* and the degree of the monomial with maximal degree: the last one
1001* but search in all components before syzcomp
1002*/
1003static int ldeg0c(poly p,int *l)
1004{
1005  int o=pFDeg(p);
1006  int ll=1;
1007
1008  if (pMaxBound/*syzComp*/==0)
1009  {
1010    while ((p=pNext(p))!=NULL)
1011    {
1012      o=pFDeg(p);
1013      ll++;
1014    }
1015  }
1016  else
1017  {
1018    while ((p=pNext(p))!=NULL)
1019    {
1020      if (pGetComp(p)<=pMaxBound/*syzComp*/)
1021      {
1022        o=pFDeg(p);
1023        ll++;
1024      }
1025      else break;
1026    }
1027  }
1028  *l=ll;
1029  return o;
1030}
1031
1032/*2
1033* compute the length of a polynomial (in l)
1034* and the degree of the monomial with maximal degree: the first one
1035* this works for the polynomial case with degree orderings
1036* (both c,dp and dp,c)
1037*/
1038static int ldegb(poly p,int *l)
1039{
1040  Exponent_t k= pGetComp(p);
1041  int o = p->Order;
1042  int ll=1;
1043
1044  while (((p=pNext(p))!=NULL) && (pGetComp(p)==k))
1045  {
1046    ll++;
1047  }
1048  *l=ll;
1049  return o;
1050}
1051
1052/*2
1053* compute the length of a polynomial (in l)
1054* and the degree of the monomial with maximal degree:
1055* this is NOT the last one, we have to look for it
1056*/
1057static int ldeg1(poly p,int *l)
1058{
1059  Exponent_t k= pGetComp(p);
1060  int ll=1;
1061  int  t,max;
1062
1063  max=pFDeg(p);
1064  while (((p=pNext(p))!=NULL) && (pGetComp(p)==k))
1065  {
1066    t=pFDeg(p);
1067    if (t>max) max=t;
1068    ll++;
1069  }
1070  *l=ll;
1071  return max;
1072}
1073
1074/*2
1075* compute the length of a polynomial (in l)
1076* and the degree of the monomial with maximal degree:
1077* this is NOT the last one, we have to look for it
1078* in all components
1079*/
1080static int ldeg1c(poly p,int *l)
1081{
1082  int ll=1;
1083  int  t,max;
1084
1085  max=pFDeg(p);
1086  while ((p=pNext(p))!=NULL)
1087  {
1088    if ((pMaxBound/*syzComp*/==0) || (pGetComp(p)<=pMaxBound/*syzComp*/))
1089    {
1090       if ((t=pFDeg(p))>max) max=t;
1091       ll++;
1092    }
1093    else break;
1094  }
1095  *l=ll;
1096  return max;
1097}
1098
1099/* -------------------------------------------------------- */
1100/* set the variables for a choosen ordering                 */
1101
1102
1103/*
1104* sets the comparision routine for monomials: for simple monomial orderings
1105* Priority is given to exponent vector
1106*/
1107static void SimpleChoose(int o_r, int comp_order, pCompProc *p)
1108{
1109  switch(o_r)
1110  {
1111      case ringorder_dp:
1112      case ringorder_wp:
1113      case ringorder_ds:
1114      case ringorder_ws:
1115      case ringorder_ls:
1116      case ringorder_unspec:
1117        pSetVarIndicies_RevLex(pVariables);
1118        pLexSgn = -1;
1119        if (comp_order == ringorder_C || o_r == ringorder_unspec)
1120        {
1121          if (pVariables1W == 1)        *p = pComp_otEXPCOMP_nwONE;
1122          else if (pVariables1W == 2)   *p = pComp_otEXPCOMP_nwTWO;
1123          else if (pVariables1W & 1)    *p = pComp_otEXPCOMP_nwODD;
1124          else                          *p = pComp_otEXPCOMP_nwEVEN;
1125        }
1126        else
1127        {
1128          // component is compatible with exponent vector
1129          if (pVariables1W == 1)        *p = pComp_otEXP_nwONE;
1130          else if (pVariables1W == 2)   *p = pComp_otEXP_nwTWO;
1131          else if (pVariables1W & 1)    *p = pComp_otEXP_nwODD;
1132          else                          *p = pComp_otEXP_nwEVEN;
1133        }
1134        break;
1135
1136#ifdef PDEBUG
1137      case ringorder_lp:
1138      case ringorder_Dp:
1139      case ringorder_Wp:
1140      case ringorder_Ds:
1141      case ringorder_Ws:
1142#else
1143      default:
1144#endif
1145        pSetVarIndicies_Lex(pVariables);
1146        pLexSgn = 1;
1147        if (comp_order == ringorder_c)
1148        {
1149          if (pVariables1W == 1)        *p = pComp_otEXPCOMP_nwONE;
1150          else if (pVariables1W == 2)   *p = pComp_otEXPCOMP_nwTWO;
1151          else if (pVariables1W & 1)    *p = pComp_otEXPCOMP_nwODD;
1152          else                          *p = pComp_otEXPCOMP_nwEVEN;
1153        }
1154        else
1155        {
1156          // component is compatible with exponent vector
1157          if (pVariables1W == 1)        *p = pComp_otEXP_nwONE;
1158          else if (pVariables1W == 2)   *p = pComp_otEXP_nwTWO;
1159          else if (pVariables1W & 1)    *p = pComp_otEXP_nwODD;
1160          else                          *p = pComp_otEXP_nwEVEN;
1161        }
1162#ifdef PDEBUG
1163        break;
1164      default:
1165        Werror("wrong internal ordering:%d at %s, l:%d\n",o_r,__FILE__,__LINE__);
1166#endif
1167  }
1168
1169  if (o_r == ringorder_lp || o_r == ringorder_ls)
1170  {
1171    pLexOrder=TRUE;
1172    pFDeg = pTotaldegree;
1173    pLDeg = ldeg1c;
1174    if (o_r == ringorder_ls)
1175      pSetVarIndicies_Lex(pVariables);
1176  }
1177}
1178
1179/*
1180* sets the comparision routine for monomials: for simple monomial orderings
1181* Priority is given to component
1182*/
1183static void SimpleChooseC(int o_r, pCompProc *p)
1184{
1185  switch(o_r)
1186  {
1187      case ringorder_dp:
1188      case ringorder_wp:
1189      case ringorder_ds:
1190      case ringorder_ls:
1191      case ringorder_ws:
1192        pSetVarIndicies_RevLex(pVariables);
1193        pLexSgn = -1;
1194        if (pVariablesW == 1)
1195          *p = pComp_otCOMPEXP_nwONE;
1196        else if (pVariablesW == 2)
1197          *p = pComp_otCOMPEXP_nwTWO;
1198        else if (pVariablesW & 1)
1199          *p = pComp_otCOMPEXP_nwODD;
1200        else
1201          *p = pComp_otCOMPEXP_nwEVEN;
1202        break;
1203
1204#ifdef PDEBUG
1205      case ringorder_lp:
1206      case ringorder_Dp:
1207      case ringorder_Wp:
1208      case ringorder_Ds:
1209      case ringorder_Ws:
1210#else
1211      default:
1212#endif
1213        pSetVarIndicies_Lex(pVariables);
1214        pLexSgn = 1;
1215        if (pVariablesW == 1)
1216          *p = pComp_otCOMPEXP_nwONE;
1217        else if (pVariablesW == 2)
1218          *p = pComp_otCOMPEXP_nwTWO;
1219        else if (pVariablesW & 1)
1220          *p = pComp_otCOMPEXP_nwODD;
1221        else
1222          *p = pComp_otCOMPEXP_nwEVEN;
1223#ifdef PDEBUG
1224        break;
1225      default:
1226        Werror("wrong internal ordering:%d at %s, l:%d\n",o_r,__FILE__,__LINE__);
1227#endif
1228  }
1229  if (o_r == ringorder_lp || o_r == ringorder_ls)
1230  {
1231    pLexOrder=TRUE;
1232    pFDeg = pTotaldegree;
1233    pLDeg = ldeg1c;
1234    if (o_r == ringorder_ls)
1235     pSetVarIndicies_Lex(pVariables);
1236  }
1237}
1238
1239/*2
1240* sets pSetm
1241* (according o_r = order of first block)
1242*/
1243static void SetpSetm(int o_r, int ip)
1244{
1245  switch(o_r)
1246  {
1247    case ringorder_lp:
1248    case ringorder_ls:
1249      if (pVariables>1)
1250        pSetm= setlex2;
1251      else
1252        pSetm= setlex1;
1253      break;
1254    case ringorder_dp:
1255    case ringorder_Dp:
1256    case ringorder_ds:
1257    case ringorder_Ds:
1258    case ringorder_unspec:
1259      pSetm= setdeg1;
1260      break;
1261    case ringorder_a:
1262    case ringorder_wp:
1263    case ringorder_Wp:
1264    case ringorder_ws:
1265    case ringorder_Ws:
1266    case ringorder_M:
1267      pSetm= setdeg1w;
1268      firstwv=polys_wv[ip];
1269      break;
1270    case ringorder_c:
1271    case ringorder_C:
1272      return;
1273      /*do not set firstBlockEnds for this orderings*/
1274#ifdef TEST
1275    default:
1276      Werror("wrong internal ordering:%d at %s, l:%d\n",o_r,__FILE__,__LINE__);
1277#endif
1278  }
1279  firstBlockEnds=block1[ip];
1280}
1281
1282
1283/*2
1284* sets the comparision routine for monomials: for all but the first
1285* block of variables (ip is the block number, o_r the number of the ordering)
1286*/
1287static void HighSet(int ip, int o_r)
1288{
1289  switch(o_r)
1290  {
1291    case ringorder_lp:
1292      bcomph[ip]=comp_lp;
1293      if (pOrdSgn==-1) pMixedOrder=TRUE;
1294      break;
1295    case ringorder_dp:
1296      bcomph[ip]=comp_dp;
1297      if (pOrdSgn==-1) pMixedOrder=TRUE;
1298      break;
1299    case ringorder_Dp:
1300      bcomph[ip]=comp_Dp;
1301      if (pOrdSgn==-1) pMixedOrder=TRUE;
1302      break;
1303    case ringorder_wp:
1304      bcomph[ip]=comp_wp;
1305      if (pOrdSgn==-1) pMixedOrder=TRUE;
1306      break;
1307    case ringorder_Wp:
1308      bcomph[ip]=comp_Wp;
1309      if (pOrdSgn==-1) pMixedOrder=TRUE;
1310      break;
1311    case ringorder_ls:
1312      bcomph[ip]=comp_ls;
1313      break;
1314    case ringorder_ds:
1315      bcomph[ip]=comp_ds;
1316      break;
1317    case ringorder_Ds:
1318      bcomph[ip]=comp_Ds;
1319      break;
1320    case ringorder_ws:
1321      bcomph[ip]=comp_ws;
1322      break;
1323    case ringorder_Ws:
1324      bcomph[ip]=comp_Ws;
1325      break;
1326    case ringorder_c:
1327      pComponentOrder=1;
1328      bcomph[ip]=comp_c;
1329      break;
1330    case ringorder_C:
1331      pComponentOrder=-1;
1332      bcomph[ip]=comp_c;
1333      break;
1334    case ringorder_M:
1335      bcomph[ip]=comp_M;
1336      pMixedOrder=TRUE;
1337      break;
1338    case ringorder_a:
1339      if (pOrdSgn==-1) pMixedOrder=TRUE;
1340      if (ip==0)
1341        bcomph[0]=comp1a;
1342      else
1343        bcomph[ip]=comp_a;
1344      break;
1345#ifdef TEST
1346    default:
1347      Werror("wrong internal ordering:%d at %s, l:%d\n",o_r,__FILE__,__LINE__);
1348#endif
1349  }
1350}
1351
1352/* -------------------------------------------------------- */
1353/*2
1354* change all variables to fit the description of the new ring
1355*/
1356
1357void pChangeRing(int n, int Sgn, int * orders, int * b0, int * b1,
1358         short ** wv)
1359{
1360  sip_sring tmpR;
1361  memset(&tmpR, 0, sizeof(sip_sring));
1362  tmpR.N = n;
1363  tmpR.OrdSgn = Sgn;
1364  tmpR.order = orders;
1365  tmpR.block0 = b0;
1366  tmpR.block1 = b1;
1367  tmpR.wvhdl = wv;
1368  pSetGlobals(&tmpR);
1369}
1370
1371void pSetGlobals(ring r, BOOLEAN complete)
1372{
1373  int i;
1374  pComponentOrder=1;
1375  if (ppNoether!=NULL) pDelete(&ppNoether);
1376#ifdef SRING
1377  pSRING=FALSE;
1378  pAltVars=r->N+1;
1379#endif
1380  pVariables = r->N;
1381
1382  // set the various size parameters and initialize memory
1383  if ((((pVariables+1)*sizeof(Exponent_t)) % sizeof(void*)) == 0)
1384    pVariables1W = (pVariables+1)*sizeof(Exponent_t) / sizeof(void*);
1385  else
1386    pVariables1W = ((pVariables+1)*sizeof(Exponent_t) / sizeof(void*)) + 1;
1387
1388  if ((((pVariables)*sizeof(Exponent_t)) % sizeof(void*)) == 0)
1389    pVariablesW = (pVariables)*sizeof(Exponent_t) / sizeof(void*);
1390  else
1391    pVariablesW = ((pVariables)*sizeof(Exponent_t) / sizeof(void*)) + 1;
1392
1393  pMonomSize = POLYSIZE + (pVariables + 1) * sizeof(Exponent_t);
1394  if ((pMonomSize % sizeof(void*)) == 0)
1395  {
1396    pMonomSizeW = pMonomSize/sizeof(void*);
1397  }
1398  else
1399  {
1400    pMonomSizeW = pMonomSize/sizeof(void*) + 1;
1401    pMonomSize = pMonomSizeW*sizeof(void*);
1402  }
1403
1404  // Set default Var Indicies
1405  pSetVarIndicies(pVariables);
1406
1407  // Initialize memory management
1408  mmSpecializeBlock(pMonomSize);
1409
1410  pOrdSgn = r->OrdSgn;
1411  pVectorOut=(r->order[0]==ringorder_c);
1412  order=r->order;
1413  block0=r->block0;
1414  block1=r->block1;
1415  firstwv=NULL;
1416  polys_wv=r->wvhdl;
1417  /*------- only one real block ----------------------*/
1418  pLexOrder=FALSE;
1419  pMixedOrder=FALSE;
1420  pFDeg=pDeg;
1421  if (pOrdSgn == 1) pLDeg = ldegb;
1422  else              pLDeg = ldeg0;
1423  /*======== ordering type is (_,c) =========================*/
1424  if ((order[0]==ringorder_unspec)
1425  ||(
1426    ((order[1]==ringorder_c)||(order[1]==ringorder_C))
1427    && (order[0]!=ringorder_M)
1428    && (order[2]==0))
1429  )
1430  {
1431    if ((order[0]!=ringorder_unspec)
1432    && (order[1]==ringorder_C))
1433      pComponentOrder=-1;
1434    if (pOrdSgn == -1) pLDeg = ldeg0c;
1435    SimpleChoose(order[0],order[1], &pComp0);
1436    SetpSetm(order[0],0);
1437  }
1438  /*======== ordering type is (c,_) =========================*/
1439  else if (((order[0]==ringorder_c)||(order[0]==ringorder_C))
1440  && (order[1]!=ringorder_M)
1441  &&  (order[2]==0))
1442  {
1443    /* pLDeg = ldeg0; is standard*/
1444    if (order[0]==ringorder_C)
1445      pComponentOrder=-1;
1446    SimpleChooseC(order[1], &pComp0);
1447    SetpSetm(order[1],1);
1448  }
1449  /*------- more than one block ----------------------*/
1450  else
1451  {
1452    //pLexOrder=TRUE;
1453    pVectorOut=order[0]==ringorder_c;
1454    if ((pVectorOut)||(order[0]==ringorder_C))
1455    {
1456      if(block1[1]!=pVariables) pLexOrder=TRUE;
1457    }
1458    else
1459    {
1460      if(block1[0]!=pVariables) pLexOrder=TRUE;
1461    }
1462    /*the number of orderings:*/
1463    i = 0;
1464    while (order[++i] != 0);
1465    do
1466    {
1467      i--;
1468      HighSet(i, order[i]);/*sets also pMixedOrder to TRUE, if...*/
1469      SetpSetm(order[i],i);
1470    }
1471    while (i != 0);
1472
1473    pComp0 = BlockComp;
1474    if ((order[0]!=ringorder_c)&&(order[0]!=ringorder_C))
1475    {
1476      pLDeg = ldeg1c;
1477    }
1478    else
1479    {
1480      pLDeg = ldeg1;
1481    }
1482    pFDeg = pWTotaldegree; // may be improved: pTotaldegree for lp/dp/ls/.. blocks
1483  }
1484  if (complete)
1485  {
1486    if ((pLexOrder) || (pOrdSgn==-1))
1487    {
1488      test &= ~Sy_bit(OPT_REDTAIL); /* noredTail */
1489    }
1490  }
1491}
1492
1493/* -------------------------------------------------------- */
1494
1495static BOOLEAN pMultT_nok;
1496/*2
1497* update the polynomial a by multipying it by
1498* the (number) coefficient
1499* and the exponent vector (of) exp (a well initialized polynomial)
1500*/
1501poly pMultT(poly a, poly exp )
1502{
1503  int i;
1504  number t,x,y=pGetCoeff(exp);
1505  poly aa=a;
1506  poly prev=NULL;
1507#ifdef SDRING
1508  poly pDRINGres=NULL;
1509#endif
1510
1511  pMultT_nok = pGetComp(exp);
1512#ifdef PDEBUG
1513  pTest(a);
1514  pTest(exp);
1515#endif
1516  while (a !=NULL)
1517  {
1518    x=pGetCoeff(a);
1519    t=nMult(x/*pGetCoeff(a)*/,y/*pGetCoeff(exp)*/);
1520    nDelete(&x/*pGetCoeff(a)*/);
1521    pSetCoeff0(a,t);
1522    if (nIsZero(t))
1523    {
1524      if (prev==NULL) { pDelete1(&a); aa=a; }
1525      else            { pDelete1(&prev->next); a=prev->next;}
1526    }
1527    else
1528    {
1529#ifdef DRING
1530      if (pDRING)
1531      {
1532         if (pdDFlag(a)==1)
1533         {
1534           if (pdDFlag(exp)==1)
1535           {
1536             pDRINGres=pAdd(pDRINGres,pMultDD(a,exp));
1537           }
1538           else
1539           {
1540             pDRINGres=pAdd(pDRINGres,pMultDT(a,exp));
1541           }
1542         }
1543         else
1544         {
1545           if (pdDFlag(exp)==1)
1546           {
1547             pDRINGres=pAdd(pDRINGres,pMultDD(a,exp));
1548           }
1549           else
1550           {
1551             pDRINGres=pAdd(pDRINGres,pMultTT(a,exp));
1552           }
1553         }
1554      }
1555      else
1556#endif
1557#ifdef SRING
1558      if (pSRING)
1559      {
1560        pDRINGres=pAdd(pDRINGres,psMultM(a,exp));
1561      }
1562      else
1563#endif
1564      {
1565        for (i=pVariables; i != 0; i--)
1566        {
1567           pAddExp(a,i, pGetExp(exp,i));
1568        }
1569        #ifdef TEST_MAC_ORDER
1570        if (bNoAdd)
1571          pSetm(a);
1572        else
1573        #endif
1574          a->Order += exp->Order;
1575        if (pMultT_nok)
1576        {
1577          if (pGetComp(a) == 0)
1578          {
1579             pSetComp(a, pGetComp(exp));
1580          }
1581          else
1582          {
1583            return NULL /*FALSE*/;
1584          }
1585        }
1586      }
1587      prev=a;
1588      pIter(a);
1589    }
1590  }
1591  pMultT_nok=0;
1592#ifdef DRING
1593   if (pDRING)
1594   {
1595     pDelete(&aa);
1596     return pDRINGres;
1597   }
1598#endif
1599#ifdef SRING
1600   if (pSRING)
1601   {
1602     pDelete(&aa);
1603     return pDRINGres;
1604   }
1605#endif
1606   return aa; /*TRUE*/
1607}
1608
1609/*2
1610* multiply p1 with p2, p1 and p2 are destroyed
1611* do not put attention on speed: the procedure is only used in the interpreter
1612*/
1613poly pMult(poly p1, poly p2)
1614{
1615  poly res, r, rn, a;
1616  BOOLEAN cont;
1617
1618  if ((p1!=NULL) && (p2!=NULL))
1619  {
1620#ifdef PDEBUG
1621    pTest(p1);
1622    pTest(p2);
1623#endif
1624    cont = TRUE;
1625    a = p1;
1626    if (pNext(p2)!=NULL)
1627      a = pCopy(a);
1628    else
1629      cont = FALSE;
1630    res = pMultT(a, p2);
1631    if (pMultT_nok)
1632    {
1633      if (cont) pDelete(&p1);
1634      pDelete(&res);
1635      pDelete(&p2);
1636      return NULL;
1637    }
1638    pDelete1(&p2);
1639    r = res;
1640    if (r!=NULL) rn = pNext(r);
1641    else         rn=NULL;
1642    while (cont)
1643    {
1644      if (pNext(p2)==NULL)
1645      {
1646        a = p1;
1647        cont = FALSE;
1648      }
1649      else
1650      {
1651        a = pCopy(p1);
1652      }
1653      a=pMultT(a, p2); //sets pMultT_nok
1654      if (pMultT_nok)
1655      {
1656        if (cont) pDelete(&p1);
1657        pDelete(&a);
1658        pDelete(&res);
1659        pDelete(&p2);
1660        return NULL;
1661      }
1662      while ((rn!=NULL) && (pComp0(rn,a)>0))
1663      {
1664        r = rn;
1665        pIter(rn);
1666      }
1667      if (r!=NULL) pNext(r) = rn = pAdd(a, rn);
1668      else         res=r=a;
1669      pDelete1(&p2);
1670    }
1671    return res;
1672  }
1673  pDelete(&p1);
1674  pDelete(&p2);
1675  return NULL;
1676}
1677
1678/*2
1679* update a by multiplying it with c (c will not be destroyed)
1680*/
1681void pMultN(poly a, number c)
1682{
1683  number t;
1684
1685  while (a!=NULL)
1686  {
1687    t=nMult(pGetCoeff(a), c);
1688    //nNormalize(t);
1689    pSetCoeff(a,t);
1690    pIter(a);
1691  }
1692}
1693
1694/*2
1695* return a copy of the poly a times the number c (a,c will not be destroyed)
1696*/
1697poly pMultCopyN(poly a, number c)
1698{
1699  poly result=NULL,hp;
1700
1701  if (a != NULL)
1702  {
1703    result=pNew();
1704    pCopy2(result,a);
1705    pNext(result)=NULL;
1706    pGetCoeff(result)=nMult(pGetCoeff(a),c);
1707    pIter(a);
1708    hp=result;
1709    while (a)
1710    {
1711      hp=pNext(hp)=pNew();
1712      pCopy2(hp,a);
1713      pSetCoeff0(hp,nMult(pGetCoeff(a), c));
1714      pIter(a);
1715    }
1716    pNext(hp)=NULL;
1717  }
1718  return result;
1719}
1720
1721/*2
1722* assumes that the head term of b is a multiple of the head term of a
1723* and return the multiplicant *m
1724*/
1725poly pDivide(poly a, poly b)
1726{
1727  int i;
1728  poly result=pInit();
1729
1730  for(i=(int)pVariables; i; i--)
1731    pSetExp(result,i, pGetExp(a,i)- pGetExp(b,i));
1732  pSetComp(result, pGetComp(a) - pGetComp(b));
1733  pSetm(result);
1734  return result;
1735}
1736
1737/*2
1738* divides a by the monomial b, ignores monomials wihich are not divisible
1739* assumes that b is not NULL
1740*/
1741poly pDivideM(poly a, poly b)
1742{
1743  if (a==NULL) return NULL;
1744  poly result=a;
1745  poly prev=NULL;
1746  int i;
1747  number inv=nInvers(pGetCoeff(b));
1748
1749  while (a!=NULL)
1750  {
1751    if (pDivisibleBy(b,a))
1752    {
1753      for(i=(int)pVariables; i; i--)
1754         pSubExp(a,i, pGetExp(b,i));
1755      pSubComp(a, pGetComp(b));
1756      pSetm(a);
1757      prev=a;
1758      pIter(a);
1759    }
1760    else
1761    {
1762      if (prev==NULL)
1763      {
1764        pDelete1(&result);
1765        a=result;
1766      }
1767      else
1768      {
1769        pDelete1(&pNext(prev));
1770        a=pNext(prev);
1771      }
1772    }
1773  }
1774  pMultN(result,inv);
1775  nDelete(&inv);
1776  pDelete(&b);
1777  return result;
1778}
1779
1780/*2
1781* returns the LCM of the head terms of a and b in *m
1782*/
1783void pLcm(poly a, poly b, poly m)
1784{
1785  int i;
1786  for (i=pVariables; i; i--)
1787  {
1788    pSetExp(m,i, max( pGetExp(a,i), pGetExp(b,i)));
1789  }
1790  pSetComp(m, max(pGetComp(a), pGetComp(b)));
1791}
1792
1793/*2
1794* convert monomial given as string to poly, e.g. 1x3y5z
1795*/
1796poly pmInit(char *st, BOOLEAN &ok)
1797{
1798  int i,j;
1799  ok=FALSE;
1800  BOOLEAN b=FALSE;
1801  poly rc = pInit();
1802  char *s = nRead(st,&(rc->coef));
1803  if (s==st)
1804  /* i.e. it does not start with a coeff: test if it is a ringvar*/
1805  {
1806    j = rIsRingVar(s);
1807    if (j >= 0)
1808    {
1809      pIncrExp(rc,1+j);
1810      goto done;
1811    }
1812  }
1813  else
1814    b=TRUE;
1815  while (*s!='\0')
1816  {
1817    char ss[2];
1818    ss[0] = *s++;
1819    ss[1] = '\0';
1820    j = rIsRingVar(ss);
1821    if (j >= 0)
1822    {
1823      s = eati(s,&i);
1824      pAddExp(rc,1+j, (Exponent_t)i);
1825    }
1826    else
1827    {
1828      if ((s!=st)&&isdigit(st[0]))
1829      {
1830        errorreported=TRUE;
1831      }
1832      pDelete(&rc);
1833      return NULL;
1834    }
1835  }
1836done:
1837  ok=!errorreported;
1838  if (nIsZero(pGetCoeff(rc))) pDelete1(&rc);
1839  else
1840  {
1841#ifdef DRING
1842    if (pDRING)
1843    {
1844      for(i=1;i<=pdN;i++)
1845      {
1846        if(pGetExp(rc,pdDX(i))>0)
1847        {
1848          pdSetDFlag(rc,1);
1849          break;
1850        }
1851      }
1852    }
1853#endif
1854    pSetm(rc);
1855  }
1856  return rc;
1857}
1858
1859/*2
1860*make p homgeneous by multiplying the monomials by powers of x_varnum
1861*/
1862poly pHomogen (poly p, int varnum)
1863{
1864  poly q=NULL;
1865  poly res;
1866  int  o,ii;
1867
1868  if (p!=NULL)
1869  {
1870    if ((varnum < 1) || (varnum > pVariables))
1871    {
1872      return NULL;
1873    }
1874    o=pWTotaldegree(p);
1875    q=pNext(p);
1876    while (q != NULL)
1877    {
1878      ii=pWTotaldegree(q);
1879      if (ii>o) o=ii;
1880      pIter(q);
1881    }
1882    q = pCopy(p);
1883    res = q;
1884    while (q != NULL)
1885    {
1886      ii = o-pWTotaldegree(q);
1887      if (ii!=0)
1888      {
1889        pAddExp(q,varnum, (Exponent_t)ii);
1890        pSetm(q);
1891      }
1892      pIter(q);
1893    }
1894    q = pOrdPolyInsertSetm(res);
1895  }
1896  return q;
1897}
1898
1899
1900/*2
1901*replaces the maximal powers of the leading monomial of p2 in p1 by
1902*the same powers of n, utility for dehomogenization
1903*/
1904poly pDehomogen (poly p1,poly p2,number n)
1905{
1906  polyset P;
1907  int     SizeOfSet=5;
1908  int     i;
1909  poly    p;
1910  number  nn;
1911
1912  P = (polyset)Alloc0(5*sizeof(poly));
1913  //for (i=0; i<5; i++)
1914  //{
1915  //  P[i] = NULL;
1916  //}
1917  pCancelPolyByMonom(p1,p2,&P,&SizeOfSet);
1918  p = P[0];
1919  //P[0] = NULL ;// for safety, may be remoeved later
1920  for (i=1; i<SizeOfSet; i++)
1921  {
1922    if (P[i] != NULL)
1923    {
1924      nPower(n,i,&nn);
1925      pMultN(P[i],nn);
1926      p = pAdd(p,P[i]);
1927      //P[i] =NULL; // for safety, may be removed later
1928      nDelete(&nn);
1929    }
1930  }
1931  Free((ADDRESS)P,SizeOfSet*sizeof(poly));
1932  return p;
1933}
1934
1935/*4
1936*Returns the exponent of the maximal power of the leading monomial of
1937*p2 in that of p1
1938*/
1939static int pGetMaxPower (poly p1,poly p2)
1940{
1941  int     i,k,res = 32000; /*a very large integer*/
1942
1943  if (p1 == NULL) return 0;
1944  for (i=1; i<=pVariables; i++)
1945  {
1946    if ( pGetExp(p2,i) != 0)
1947    {
1948      k =  pGetExp(p1,i) /  pGetExp(p2,i);
1949      if (k < res) res = k;
1950    }
1951  }
1952  return res;
1953}
1954
1955/*2
1956*Returns as i-th entry of P the coefficient of the (i-1) power of
1957*the leading monomial of p2 in p1
1958*/
1959void pCancelPolyByMonom (poly p1,poly p2,polyset * P,int * SizeOfSet)
1960{
1961  int   maxPow;
1962  poly  p,qp,Coeff;
1963
1964  if (*P == NULL)
1965  {
1966    *P = (polyset) Alloc(5*sizeof(poly));
1967    *SizeOfSet = 5;
1968  }
1969  p = pCopy(p1);
1970  while (p != NULL)
1971  {
1972    qp = p->next;
1973    p->next = NULL;
1974    maxPow = pGetMaxPower(p,p2);
1975    Coeff = pDivByMonom(p,p2);
1976    if (maxPow > *SizeOfSet)
1977    {
1978      pEnlargeSet(P,*SizeOfSet,maxPow+1-*SizeOfSet);
1979      *SizeOfSet = maxPow+1;
1980    }
1981    (*P)[maxPow] = pAdd((*P)[maxPow],Coeff);
1982    pDelete(&p);
1983    p = qp;
1984  }
1985}
1986
1987/*2
1988*returns the leading monomial of p1 divided by the maximal power of that
1989*of p2
1990*/
1991poly pDivByMonom (poly p1,poly p2)
1992{
1993  int     k, i;
1994
1995  if (p1 == NULL) return NULL;
1996  k = pGetMaxPower(p1,p2);
1997  if (k == 0)
1998    return pHead(p1);
1999  else
2000  {
2001    number n;
2002    poly p = pInit();
2003
2004    p->next = NULL;
2005    for (i=1; i<=pVariables; i++)
2006    {
2007       pSetExp(p,i, pGetExp(p1,i)-k* pGetExp(p2,i));
2008    }
2009    nPower(p2->coef,k,&n);
2010    pSetCoeff0(p,nDiv(p1->coef,n));
2011    nDelete(&n);
2012    pSetm(p);
2013    return p;
2014  }
2015}
2016/*----------utilities for syzygies--------------*/
2017poly pTakeOutComp(poly * p, int k)
2018{
2019  poly q = *p,qq=NULL,result = NULL;
2020
2021  if (q==NULL) return NULL;
2022  if (pGetComp(q)==k)
2023  {
2024    result = q;
2025    while ((q!=NULL) && (pGetComp(q)==k))
2026    {
2027      pSetComp(q,0);
2028      qq = q;
2029      pIter(q);
2030    }
2031    *p = q;
2032    pNext(qq) = NULL;
2033  }
2034  if (q==NULL) return result;
2035  if (pGetComp(q) > k) pDecrComp(q);
2036  poly pNext_q;
2037  while ((pNext_q=pNext(q))!=NULL)
2038  {
2039    if (pGetComp(pNext_q)==k)
2040    {
2041      if (result==NULL)
2042      {
2043        result = pNext_q;
2044        qq = result;
2045      }
2046      else
2047      {
2048        pNext(qq) = pNext_q;
2049        pIter(qq);
2050      }
2051      pNext(q) = pNext(pNext_q);
2052      pNext(qq) =NULL;
2053      pSetComp(qq,0);
2054    }
2055    else
2056    {
2057      /*pIter(q);*/ q=pNext_q;
2058      if (pGetComp(q) > k) pDecrComp(q);
2059    }
2060  }
2061  return result;
2062}
2063
2064poly pTakeOutComp1(poly * p, int k)
2065{
2066  poly q = *p;
2067
2068  if (q==NULL) return NULL;
2069
2070  poly qq=NULL,result = NULL;
2071
2072  if (pGetComp(q)==k)
2073  {
2074    result = q; /* *p */
2075    while ((q!=NULL) && (pGetComp(q)==k))
2076    {
2077      pSetComp(q,0);
2078      qq = q;
2079      pIter(q);
2080    }
2081    *p = q;
2082    pNext(qq) = NULL;
2083  }
2084  if (q==NULL) return result;
2085//  if (pGetComp(q) > k) pGetComp(q)--;
2086  while (pNext(q)!=NULL)
2087  {
2088    if (pGetComp(pNext(q))==k)
2089    {
2090      if (result==NULL)
2091      {
2092        result = pNext(q);
2093        qq = result;
2094      }
2095      else
2096      {
2097        pNext(qq) = pNext(q);
2098        pIter(qq);
2099      }
2100      pNext(q) = pNext(pNext(q));
2101      pNext(qq) =NULL;
2102      pSetComp(qq,0);
2103    }
2104    else
2105    {
2106      pIter(q);
2107//      if (pGetComp(q) > k) pGetComp(q)--;
2108    }
2109  }
2110  return result;
2111}
2112
2113void pDeleteComp(poly * p,int k)
2114{
2115  poly q;
2116
2117  while ((*p!=NULL) && (pGetComp(*p)==k)) pDelete1(p);
2118  if (*p==NULL) return;
2119  q = *p;
2120  if (pGetComp(q)>k) pDecrComp(q);
2121  while (pNext(q)!=NULL)
2122  {
2123    if (pGetComp(pNext(q))==k)
2124      pDelete1(&(pNext(q)));
2125    else
2126    {
2127      pIter(q);
2128      if (pGetComp(q)>k) pDecrComp(q);
2129    }
2130  }
2131}
2132/*----------end of utilities for syzygies--------------*/
2133
2134/*2
2135* pair has no common factor ? or is no polynomial
2136*/
2137BOOLEAN pHasNotCF(poly p1, poly p2)
2138{
2139#ifdef SRING
2140  if (pSRING)
2141    return FALSE;
2142#endif
2143
2144  if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
2145    return FALSE;
2146  int i = 1;
2147  loop
2148  {
2149    if ((pGetExp(p1, i) > 0) && (pGetExp(p2, i) > 0))   return FALSE;
2150    if (i == pVariables)                                return TRUE;
2151    i++;
2152  }
2153}
2154
2155
2156/*2
2157*should return 1 if p divides q and p<q,
2158*             -1 if q divides p and q<p
2159*              0 otherwise
2160*/
2161int     pDivComp(poly p, poly q)
2162{
2163  if (pGetComp(p) == pGetComp(q))
2164  {
2165    int i=pVariables;
2166    long d;
2167    BOOLEAN a=FALSE, b=FALSE;
2168    for (; i>0; i--)
2169    {
2170      d = pGetExpDiff(p, q, i);
2171      if (d)
2172      {
2173        if (d < 0)
2174        {
2175          if (b) return 0;
2176          a =TRUE;
2177        }
2178        else
2179        {
2180          if (a) return 0;
2181          b = TRUE;
2182        }
2183      }
2184    }
2185    if (a) return 1;
2186    else if (b)  return -1;
2187  }
2188  return 0;
2189}
2190/*2
2191*divides p1 by its leading monomial
2192*/
2193void pNorm(poly p1)
2194{
2195  poly h;
2196  number k, c;
2197
2198  if (p1!=NULL)
2199  {
2200    if (!nIsOne(pGetCoeff(p1)))
2201    {
2202      nNormalize(pGetCoeff(p1));
2203      k=pGetCoeff(p1);
2204      c = nInit(1);
2205      pSetCoeff0(p1,c);
2206      h = pNext(p1);
2207      while (h!=NULL)
2208      {
2209        c=nDiv(pGetCoeff(h),k);
2210        if (!nIsOne(c)) nNormalize(c);
2211        pSetCoeff(h,c);
2212        pIter(h);
2213      }
2214      nDelete(&k);
2215    }
2216    else
2217    {
2218      h = pNext(p1);
2219      while (h!=NULL)
2220      {
2221        nNormalize(pGetCoeff(h));
2222        pIter(h);
2223      }
2224    }
2225  }
2226}
2227
2228/*2
2229*normalize all coeffizients
2230*/
2231void pNormalize(poly p)
2232{
2233  while (p!=NULL)
2234  {
2235    nTest(pGetCoeff(p));
2236    nNormalize(pGetCoeff(p));
2237    pIter(p);
2238  }
2239}
2240
2241/*3
2242* substitute the n-th variable by 1 in p
2243* destroy p
2244*/
2245static poly pSubst1 (poly p,int n)
2246{
2247  poly qq,result = NULL;
2248
2249  while (p != NULL)
2250  {
2251    qq = p;
2252    pIter(p);
2253    qq->next = NULL;
2254    pSetExp(qq,n,0);
2255    pSetm(qq);
2256    result = pAdd(result,qq);
2257  }
2258  return result;
2259}
2260
2261/*2
2262* substitute the n-th variable by e in p
2263* destroy p
2264*/
2265poly pSubst(poly p, int n, poly e)
2266{
2267  if ((e!=NULL)&&(pIsConstant(e))&&(nIsOne(pGetCoeff(e))))
2268    return pSubst1(p,n);
2269
2270  int exponent,i;
2271  poly h, res, m;
2272  Exponent_t *me,*ee;
2273  number nu,nu1;
2274
2275  me=(Exponent_t *)Alloc((pVariables+1)*sizeof(Exponent_t));
2276  ee=(Exponent_t *)Alloc((pVariables+1)*sizeof(Exponent_t));
2277  if (e!=NULL) pGetExpV(e,ee);
2278  res=NULL;
2279  h=p;
2280  while (h!=NULL)
2281  {
2282    if ((e!=NULL) || (pGetExp(h,n)==0))
2283    {
2284      m=pHead(h);
2285      pGetExpV(m,me);
2286      exponent=me[n];
2287      me[n]=0;
2288      for(i=pVariables;i>0;i--)
2289        me[i]+=exponent*ee[i];
2290      pSetExpV(m,me);
2291      if (e!=NULL)
2292      {
2293        nPower(pGetCoeff(e),exponent,&nu);
2294        nu1=nMult(pGetCoeff(m),nu);
2295        nDelete(&nu);
2296        pSetCoeff(m,nu1);
2297      }
2298      res=pAdd(res,m);
2299    }
2300    pDelete1(&h);
2301  }
2302  Free((ADDRESS)me,(pVariables+1)*sizeof(Exponent_t));
2303  Free((ADDRESS)ee,(pVariables+1)*sizeof(Exponent_t));
2304  return res;
2305}
2306
2307BOOLEAN pCompareChain (poly p,poly p1,poly p2,poly lcm)
2308{
2309  int k, j;
2310
2311  if (lcm==NULL) return FALSE;
2312
2313  for (j=pVariables; j; j--)
2314    if ( pGetExp(p,j) >  pGetExp(lcm,j)) return FALSE;
2315  if ( pGetComp(p) !=  pGetComp(lcm)) return FALSE;
2316  for (j=pVariables; j; j--)
2317  {
2318    if (pGetExp(p1,j)!=pGetExp(lcm,j))
2319    {
2320      if (pGetExp(p,j)!=pGetExp(lcm,j))
2321      {
2322        for (k=pVariables; k>j; k--)
2323        {
2324          if ((pGetExp(p,k)!=pGetExp(lcm,k))
2325          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
2326            return TRUE;
2327        }
2328        for (k=j-1; k; k--)
2329        {
2330          if ((pGetExp(p,k)!=pGetExp(lcm,k))
2331          && (pGetExp(p2,k)!=pGetExp(lcm,k)))
2332            return TRUE;
2333        }
2334        return FALSE;
2335      }
2336    }
2337    else if (pGetExp(p2,j)!=pGetExp(lcm,j))
2338    {
2339      if (pGetExp(p,j)!=pGetExp(lcm,j))
2340      {
2341        for (k=pVariables; k>j; k--)
2342        {
2343          if ((pGetExp(p,k)!=pGetExp(lcm,k))
2344          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
2345            return TRUE;
2346        }
2347        for (k=j-1; k!=0 ; k--)
2348        {
2349          if ((pGetExp(p,k)!=pGetExp(lcm,k))
2350          && (pGetExp(p1,k)!=pGetExp(lcm,k)))
2351            return TRUE;
2352        }
2353        return FALSE;
2354      }
2355    }
2356  }
2357  return FALSE;
2358}
2359
2360int pWeight(int i)
2361{
2362  if ((firstwv==NULL) || (i>firstBlockEnds))
2363  {
2364    return 1;
2365  }
2366  return firstwv[i-1];
2367}
Note: See TracBrowser for help on using the repository browser.