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"
27 #include "polys/monomials/p_polys.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"
34 #include "kernel/combinatorics/stairc.h"
35 #include "kernel/spectrum/multicnt.h"
36 #include "kernel/spectrum/GMPrat.h"
37 #include "kernel/spectrum/kmatrix.h"
38 #include "kernel/spectrum/npolygon.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
hasTermOfDegree(poly h,int d,const ring r)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
hasConstTerm(poly h,const ring r)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
hasLinearTerm(poly h,const ring r)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
hasAxis(ideal J,int k,const ring r)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
hasOne(ideal J,const ring r)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
isMultiple(poly f,poly m,const ring r)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
computeWC(const newtonPolygon & np,Rational max_weight,const ring r)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
normalFormHC(poly f,poly hc,const ring r)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
normalFormZ(poly f,poly Z,const ring r)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
isLeadMonomial(poly m,ideal stdJ,const ring r)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
setExp(poly m,int * r,const ring s)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
isWell(const ring r)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
computeNF(ideal stdJ,poly hc,poly wc,spectrumPolyList * NF,const ring r)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
ringIsLocal(const ring r)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 // ----------------------------------------------------------------------------
524