source: git/kernel/spectrum.cc @ fbc7cb

jengelh-datetimespielwiese
Last change on this file since fbc7cb was ba5e9e, checked in by Oleksandr Motsak <motsak@…>, 9 years ago
Changed configure-scripts to generate individual public config files for each package: resources, libpolys, singular (main) fix: sources should include correct corresponding config headers.
  • Property mode set to 100644
File size: 12.4 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#ifdef HAVE_CONFIG_H
11#include "singularconfig.h"
12#endif /* HAVE_CONFIG_H */
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/kstd1.h>
34#include <kernel/stairc.h>
35#include <kernel/multicnt.h>
36#include <kernel/GMPrat.h>
37#include <kernel/kmatrix.h>
38#include <kernel/npolygon.h>
39#include <kernel/splist.h>
40#include <kernel/semic.h>
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_Divide( 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);
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// ----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.