source: git/kernel/spectrum.cc @ 6ce030f

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