My Project
Loading...
Searching...
No Matches
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"
31//extern ring currRing;
32
41
42// ----------------------------------------------------------------------------
43// test if the polynomial h has a term of total degree d
44// ----------------------------------------------------------------------------
45
46BOOLEAN 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
63static 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
72static 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
81BOOLEAN 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
96int 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
110int 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
142poly 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
182static 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
206static 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
234static 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
260static 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
274static 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
309void 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
461BOOLEAN ringIsLocal( const ring r )
462{
463 poly m = p_One(r);
464 poly one = p_One(r);
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:1492
void p_Norm(poly p1, const ring r)
Definition: p_polys.cc:3719
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