My Project
spectrum.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // spectrum.cc
3 // begin of file
4 // Stephan Endrass, endrass@mathematik.uni-mainz.de
5 // 23.7.99
6 // ----------------------------------------------------------------------------
7 
8 #define SPECTRUM_CC
9 
10 
11 
12 
13 #include "kernel/mod2.h"
14 
15 #ifdef HAVE_SPECTRUM
16 
17 #ifdef SPECTRUM_PRINT
18 #include <iostream.h>
19 #ifndef SPECTRUM_IOSTREAM
20 #include <stdio.h>
21 #endif
22 #endif
23 
24 #include "misc/mylimits.h"
25 
26 #include "coeffs/numbers.h"
28 #include "polys/simpleideals.h"
29 #include "misc/intvec.h"
30 #include "polys/monomials/ring.h"
31 //extern ring currRing;
32 
33 #include "kernel/GBEngine/kstd1.h"
36 #include "kernel/spectrum/GMPrat.h"
39 #include "kernel/spectrum/splist.h"
40 #include "kernel/spectrum/semic.h"
41 
42 // ----------------------------------------------------------------------------
43 // test if the polynomial h has a term of total degree d
44 // ----------------------------------------------------------------------------
45 
46 BOOLEAN hasTermOfDegree( poly h, int d, const ring r )
47 {
48  do
49  {
50  if( p_Totaldegree( h,r )== d )
51  return TRUE;
52  pIter(h);
53  }
54  while( h!=NULL );
55 
56  return FALSE;
57 }
58 
59 // ----------------------------------------------------------------------------
60 // test if the polynomial h has a constant term
61 // ----------------------------------------------------------------------------
62 
63 static BOOLEAN inline hasConstTerm( poly h, const ring r )
64 {
65  return hasTermOfDegree(h,0,r);
66 }
67 
68 // ----------------------------------------------------------------------------
69 // test if the polynomial h has a linear term
70 // ----------------------------------------------------------------------------
71 
72 static BOOLEAN inline hasLinearTerm( poly h, const ring r )
73 {
74  return hasTermOfDegree(h,1,r);
75 }
76 
77 // ----------------------------------------------------------------------------
78 // test if the ideal J has a lead monomial on the axis number k
79 // ----------------------------------------------------------------------------
80 
81 BOOLEAN hasAxis( ideal J,int k, const ring r )
82 {
83  int i;
84 
85  for( i=0; i<IDELEMS(J); i++ )
86  {
87  if (p_IsPurePower( J->m[i],r ) == k) return TRUE;
88  }
89  return FALSE;
90 }
91 
92 // ----------------------------------------------------------------------------
93 // test if the ideal J has 1 as a lead monomial
94 // ----------------------------------------------------------------------------
95 
96 int hasOne( ideal J, const ring r )
97 {
98  int i;
99 
100  for( i=0; i<IDELEMS(J); i++ )
101  {
102  if (p_IsConstant(J->m[i],r)) return TRUE;
103  }
104  return FALSE;
105 }
106 // ----------------------------------------------------------------------------
107 // test if m is a multiple of one of the monomials of f
108 // ----------------------------------------------------------------------------
109 
110 int isMultiple( poly f,poly m, const ring r )
111 {
112  while( f!=NULL )
113  {
114  // ---------------------------------------------------
115  // for a local order f|m is only possible if f>=m
116  // ---------------------------------------------------
117 
118  if( p_LmCmp( f,m,r )>=0 )
119  {
120  if( p_LmDivisibleByNoComp( f,m,r ) )
121  {
122  return TRUE;
123  }
124  else
125  {
126  pIter( f );
127  }
128  }
129  else
130  {
131  return FALSE;
132  }
133  }
134 
135  return FALSE;
136 }
137 
138 // ----------------------------------------------------------------------------
139 // compute the minimal monomial of minimmal weight>=max_weight
140 // ----------------------------------------------------------------------------
141 
142 poly computeWC( const newtonPolygon &np,Rational max_weight, const ring r )
143 {
144  poly m = p_One(r);
145  poly wc = NULL;
146  int mdegree;
147 
148  for( int i=1; i<=r->N; i++ )
149  {
150  mdegree = 1;
151  p_SetExp( m,i,mdegree,r );
152  // pSetm( m );
153  // np.weight_shift does not need pSetm( m ), postpone it
154 
155  while( np.weight_shift( m,r )<max_weight )
156  {
157  mdegree++;
158  p_SetExp( m,i,mdegree,r );
159  // pSetm( m );
160  // np.weight_shift does not need pSetm( m ), postpone it
161  }
162  p_Setm( m,r );
163 
164  if( i==1 || p_Cmp( m,wc,r )<0 )
165  {
166  p_Delete( &wc,r );
167  wc = p_Head( m,r );
168  }
169 
170  p_SetExp( m,i,0,r );
171  }
172 
173  p_Delete( &m,r );
174 
175  return wc;
176 }
177 
178 // ----------------------------------------------------------------------------
179 // deletes all monomials of f which are less than hc
180 // ----------------------------------------------------------------------------
181 
182 static inline poly normalFormHC( poly f,poly hc, const ring r )
183 {
184  poly *ptr = &f;
185 
186  while( (*ptr)!=NULL )
187  {
188  if( p_LmCmp( *ptr,hc,r )>=0 )
189  {
190  ptr = &(pNext( *ptr ));
191  }
192  else
193  {
194  p_Delete( ptr,r );
195  return f;
196  }
197  }
198 
199  return f;
200 }
201 
202 // ----------------------------------------------------------------------------
203 // deletes all monomials of f which are multiples of monomials of Z
204 // ----------------------------------------------------------------------------
205 
206 static inline poly normalFormZ( poly f,poly Z, const ring r )
207 {
208  poly *ptr = &f;
209 
210  while( (*ptr)!=NULL )
211  {
212  if( !isMultiple( Z,*ptr,r ) )
213  {
214  ptr = &(pNext( *ptr ));
215  }
216  else
217  {
218  p_LmDelete(ptr,r);
219  }
220  }
221 
222  return f;
223 }
224 
225 // ?? HS:
226 // Looks for the shortest polynomial f in stdJ which is divisible
227 // by the monimial m, return its index in stdJ
228 // ----------------------------------------------------------------------------
229 // Looks for the first polynomial f in stdJ which satisfies m=LT(f)*eta
230 // for a monomial eta. The return eta*f-m and cancel all monomials
231 // which are smaller than the highest corner hc
232 // ----------------------------------------------------------------------------
233 
234 static inline int isLeadMonomial( poly m,ideal stdJ, const ring r )
235 {
236  int length = MAX_INT_VAL;
237  int result = -1;
238 
239  for( int i=0; i<IDELEMS(stdJ); i++ )
240  {
241  if( p_Cmp( stdJ->m[i],m,r )>=0 && p_DivisibleBy( stdJ->m[i],m,r ) )
242  {
243  int tmp = pLength( stdJ->m[i] );
244 
245  if( tmp < length )
246  {
247  length = tmp;
248  result = i;
249  }
250  }
251  }
252 
253  return result;
254 }
255 
256 // ----------------------------------------------------------------------------
257 // set the exponent of a monomial t an integer array
258 // ----------------------------------------------------------------------------
259 
260 static void setExp( poly m,int *r, const ring s )
261 {
262  for( int i=s->N; i>0; i-- )
263  {
264  p_SetExp( m,i,r[i-1],s );
265  }
266  p_Setm( m,s );
267 }
268 
269 // ----------------------------------------------------------------------------
270 // check if the ordering is a reverse wellordering, i.e. every monomial
271 // is smaller than only finitely monomials
272 // ----------------------------------------------------------------------------
273 
274 static BOOLEAN isWell( const ring r )
275 {
276  int b = rBlocks( r );
277 
278  if( b==3 &&
279  ( r->order[0] == ringorder_ds ||
280  r->order[0] == ringorder_Ds ||
281  r->order[0] == ringorder_ws ||
282  r->order[0] == ringorder_Ws ) )
283  {
284  return TRUE;
285  }
286  else if( b>=3
287  && (( r->order[0] ==ringorder_a
288  && r->block1[0]==r->N )
289  || (r->order[0]==ringorder_M
290  && r->block1[0]==r->N*r->N )))
291  {
292  for( int i=r->N-1; i>=0; i-- )
293  {
294  if( r->wvhdl[0][i]>=0 )
295  {
296  return FALSE;
297  }
298  }
299  return TRUE;
300  }
301 
302  return FALSE;
303 }
304 
305 // ----------------------------------------------------------------------------
306 // compute all monomials not in stdJ and their normal forms
307 // ----------------------------------------------------------------------------
308 
309 void computeNF( ideal stdJ,poly hc,poly wc,spectrumPolyList *NF, const ring r )
310 {
311  int carry,k;
312  multiCnt C( r->N,0 );
313  poly Z = NULL;
314 
315  int well = isWell(r);
316 
317  do
318  {
319  poly m = p_One(r);
320  setExp( m,C.cnt,r );
321 
322  carry = FALSE;
323 
324  k = isLeadMonomial( m,stdJ,r );
325 
326  if( k < 0 )
327  {
328  // ---------------------------
329  // m is not a lead monomial
330  // ---------------------------
331 
332  NF->insert_node( m,NULL,r );
333  }
334  else if( isMultiple( Z,m,r ) )
335  {
336  // ------------------------------------
337  // m is trivially in the ideal stdJ
338  // ------------------------------------
339 
340  p_Delete( &m,r );
341  carry = TRUE;
342  }
343  else if( p_Cmp( m,hc,r ) < 0 || p_Cmp( m,wc,r ) < 0 )
344  {
345  // -------------------
346  // we do not need m
347  // -------------------
348 
349  p_Delete( &m,r );
350  carry = TRUE;
351  }
352  else
353  {
354  // --------------------------
355  // compute lazy normal form
356  // --------------------------
357 
358  poly multiplicant = p_MDivide( m,stdJ->m[k],r );
359  pGetCoeff( multiplicant ) = n_Init(1,r->cf);
360 
361  poly nf = p_Mult_mm( p_Copy( stdJ->m[k],r ), multiplicant,r );
362 
363  p_Delete( &multiplicant,r );
364 
365  nf = normalFormHC( nf,hc,r );
366 
367  if( pNext( nf )==NULL )
368  {
369  // ----------------------------------
370  // normal form of m is m itself
371  // ----------------------------------
372 
373  p_Delete( &nf,r );
374  NF->delete_monomial( m,r );
375  Z = p_Add_q( Z,m,r );
376  carry = TRUE;
377  }
378  else
379  {
380  nf = normalFormZ( nf,Z,r );
381 
382  if( pNext( nf )==NULL )
383  {
384  // ----------------------------------
385  // normal form of m is m itself
386  // ----------------------------------
387 
388  p_Delete( &nf,r );
389  NF->delete_monomial( m,r );
390  Z = p_Add_q( Z,m,r );
391  carry = TRUE;
392  }
393  else
394  {
395  // ------------------------------------
396  // normal form of m is a polynomial
397  // ------------------------------------
398 
399  p_Norm( nf,r );
400  if( well==TRUE )
401  {
402  NF->insert_node( m,nf,r );
403  }
404  else
405  {
406  poly nfhard = kNF( stdJ,(ideal)NULL,pNext( nf ),0,0 );
407  nfhard = normalFormHC( nfhard,hc,r );
408  nfhard = normalFormZ ( nfhard,Z,r );
409 
410  if( nfhard==NULL )
411  {
412  NF->delete_monomial( m,r );
413  Z = p_Add_q( Z,m,r );
414  carry = TRUE;
415  }
416  else
417  {
418  p_Delete( &pNext( nf ),r );
419  pNext( nf ) = nfhard;
420  NF->insert_node( m,nf,r );
421  }
422  }
423  }
424  }
425  }
426  }
427  while( C.inc( carry ) );
428 
429  // delete single monomials
430 
431  BOOLEAN not_finished;
432 
433  do
434  {
435  not_finished = FALSE;
436 
437  spectrumPolyNode *node = NF->root;
438 
439  while( node!=(spectrumPolyNode*)NULL )
440  {
441  if( node->nf!=NULL && pNext( node->nf )==NULL )
442  {
443  NF->delete_monomial( node->nf,r );
444  not_finished = TRUE;
445  node = (spectrumPolyNode*)NULL;
446  }
447  else
448  {
449  node = node->next;
450  }
451  }
452  } while( not_finished );
453 
454  p_Delete( &Z,r );
455 }
456 
457 // ----------------------------------------------------------------------------
458 // check if currRing is local
459 // ----------------------------------------------------------------------------
460 
461 BOOLEAN ringIsLocal( const ring r )
462 {
463  poly m = p_One(r);
464  poly one = p_One(r);
465  BOOLEAN res=TRUE;
466 
467  for( int i=r->N; i>0; i-- )
468  {
469  p_SetExp( m,i,1,r );
470  p_Setm( m,r );
471 
472  if( p_Cmp( m,one,r )>0 )
473  {
474  res=FALSE;
475  break;
476  }
477  p_SetExp( m,i,0,r );
478  }
479 
480  p_Delete( &m,r );
481  p_Delete( &one,r );
482 
483  return res;
484 }
485 
486 // ----------------------------------------------------------------------------
487 // print error message corresponding to spectrumState state:
488 // ----------------------------------------------------------------------------
489 /*void spectrumPrintError(spectrumState state)
490 {
491  switch( state )
492  {
493  case spectrumZero:
494  WerrorS( "polynomial is zero" );
495  break;
496  case spectrumBadPoly:
497  WerrorS( "polynomial has constant term" );
498  break;
499  case spectrumNoSingularity:
500  WerrorS( "not a singularity" );
501  break;
502  case spectrumNotIsolated:
503  WerrorS( "the singularity is not isolated" );
504  break;
505  case spectrumNoHC:
506  WerrorS( "highest corner cannot be computed" );
507  break;
508  case spectrumDegenerate:
509  WerrorS( "principal part is degenerate" );
510  break;
511  case spectrumOK:
512  break;
513 
514  default:
515  WerrorS( "unknown error occurred" );
516  break;
517  }
518 }*/
519 #endif /* HAVE_SPECTRUM */
520 // ----------------------------------------------------------------------------
521 // spectrum.cc
522 // end of file
523 // ----------------------------------------------------------------------------
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
CanonicalForm b
Definition: cfModGcd.cc:4103
FILE * f
Definition: checklibs.c:9
int * cnt
Definition: multicnt.h:21
void inc(void)
Definition: multicnt.cc:154
Rational weight_shift(poly, const ring r) const
Definition: npolygon.cc:585
void delete_monomial(poly, const ring)
Definition: splist.cc:269
spectrumPolyNode * root
Definition: splist.h:60
void insert_node(poly, poly, const ring)
Definition: splist.cc:184
spectrumPolyNode * next
Definition: splist.h:39
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:535
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
CanonicalForm res
Definition: facAbsFact.cc:60
static BOOLEAN length(leftv result, leftv arg)
Definition: interval.cc:257
STATIC_VAR Poly * h
Definition: janet.cc:971
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:3182
#define pIter(p)
Definition: monomials.h:37
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
const int MAX_INT_VAL
Definition: mylimits.h:12
#define NULL
Definition: omList.c:12
int p_IsPurePower(const poly p, const ring r)
return i, if head depends only on var(i)
Definition: p_polys.cc:1226
poly p_MDivide(poly a, poly b, const ring r)
Definition: p_polys.cc:1488
void p_Norm(poly p1, const ring r)
Definition: p_polys.cc:3715
poly p_One(const ring r)
Definition: p_polys.cc:1313
static int pLength(poly a)
Definition: p_polys.h:188
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:934
static void p_LmDelete(poly p, const ring r)
Definition: p_polys.h:721
static int p_Cmp(poly p1, poly p2, ring r)
Definition: p_polys.h:1725
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent @Note: VarOffset encodes the position in p->exp
Definition: p_polys.h:486
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:231
static poly p_Head(const poly p, const ring r)
copy the (leading) term of p
Definition: p_polys.h:858
static int p_LmCmp(poly p, poly q, const ring r)
Definition: p_polys.h:1578
static BOOLEAN p_LmDivisibleByNoComp(poly a, poly b, const ring r)
Definition: p_polys.h:1875
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition: p_polys.h:1962
static BOOLEAN p_DivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1898
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:899
static poly p_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:1049
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:844
static long p_Totaldegree(poly p, const ring r)
Definition: p_polys.h:1505
static int rBlocks(const ring r)
Definition: ring.h:568
@ ringorder_a
Definition: ring.h:70
@ ringorder_ds
Definition: ring.h:84
@ ringorder_Ds
Definition: ring.h:85
@ ringorder_ws
Definition: ring.h:86
@ ringorder_Ws
Definition: ring.h:87
@ ringorder_M
Definition: ring.h:74
#define IDELEMS(i)
Definition: simpleideals.h:23
BOOLEAN hasAxis(ideal J, int k, const ring r)
Definition: spectrum.cc:81
int hasOne(ideal J, const ring r)
Definition: spectrum.cc:96
BOOLEAN ringIsLocal(const ring r)
Definition: spectrum.cc:461
static poly normalFormHC(poly f, poly hc, const ring r)
Definition: spectrum.cc:182
static BOOLEAN hasConstTerm(poly h, const ring r)
Definition: spectrum.cc:63
static void setExp(poly m, int *r, const ring s)
Definition: spectrum.cc:260
poly computeWC(const newtonPolygon &np, Rational max_weight, const ring r)
Definition: spectrum.cc:142
static BOOLEAN isWell(const ring r)
Definition: spectrum.cc:274
int isMultiple(poly f, poly m, const ring r)
Definition: spectrum.cc:110
static poly normalFormZ(poly f, poly Z, const ring r)
Definition: spectrum.cc:206
static int isLeadMonomial(poly m, ideal stdJ, const ring r)
Definition: spectrum.cc:234
static BOOLEAN hasLinearTerm(poly h, const ring r)
Definition: spectrum.cc:72
BOOLEAN hasTermOfDegree(poly h, int d, const ring r)
Definition: spectrum.cc:46
void computeNF(ideal stdJ, poly hc, poly wc, spectrumPolyList *NF, const ring r)
Definition: spectrum.cc:309
Definition: gnumpfl.cc:25