source: git/kernel/spectrum.cc @ 210e07

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