source: git/Singular/polys.cc @ 76352f9

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