source: git/kernel/spectrum.cc @ f599a46

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