source: git/Singular/spectrum.cc @ 5d32fd

spielwiese
Last change on this file since 5d32fd was 5d32fd, checked in by Hans Schönemann <hannes@…>, 24 years ago
* hannes: move spectrum-stuff to a library and extra.cc, step I git-svn-id: file:///usr/local/Singular/svn/trunk@3917 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 31.6 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 "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 <limits.h>
22
23#include "numbers.h"
24#include "polys.h"
25#include "ipid.h"
26#include "ideals.h"
27#include "kstd1.h"
28#include "stairc.h"
29#include "lists.h"
30#include "intvec.h"
31#include "ring.h"
32
33#include "multicnt.h"
34#include "GMPrat.h"
35#include "kmatrix.h"
36#include "npolygon.h"
37#include "splist.h"
38#include "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 )
45{
46  do
47  {
48    if( pTotaldegree( h )== d )
49      return  TRUE;
50    pIter(h);
51  }
52  while( h!=(poly)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 )
62{
63  return  hasTermOfDegree(h,0);
64}
65
66// ----------------------------------------------------------------------------
67//  test if the polynomial  h  has a linear term
68// ----------------------------------------------------------------------------
69
70static BOOLEAN inline hasLinearTerm( poly h )
71{
72  return  hasTermOfDegree(h,1);
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 )
80{
81  int i;
82
83  for( i=0; i<IDELEMS(J); i++ )
84  {
85    if (pIsPurePower( J->m[i] ) == 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 )
95{
96  int i;
97
98  for( i=0; i<IDELEMS(J); i++ )
99  {
100    if (pIsConstant(J->m[i])) 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 )
109{
110  while( f!=(poly)NULL )
111  {
112    // ---------------------------------------------------
113    //  for a local order  f|m  is only possible if  f>=m
114    // ---------------------------------------------------
115
116    if( pComp0( f,m )>=0 )
117    {
118      if( pDivisibleBy2( f,m ) )
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 )
141{
142  poly m  = pOne();
143  poly wc = (poly)NULL;
144  int  mdegree;
145
146  for( int i=1; i<=pVariables; i++ )
147  {
148    mdegree = 1;
149    pSetExp( m,i,mdegree );
150    // pSetm( m );
151    // np.weight_shift does not need pSetm( m ), postpone it
152
153    while(  np.weight_shift( m )<max_weight )
154    {
155      mdegree++;
156      pSetExp( m,i,mdegree );
157      // pSetm( m );
158      // np.weight_shift does not need pSetm( m ), postpone it
159    }
160    pSetm( m );
161
162    if( i==1 || pComp( m,wc )<0 )
163    {
164      pDelete( &wc );
165      wc = pHead( m );
166    }
167
168    pSetExp( m,i,0 );
169  }
170
171  pDelete( &m );
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 )
181{
182  poly    *ptr = &f;
183
184  while( (*ptr)!=(poly)NULL )
185  {
186    if( pComp0( *ptr,hc )>=0 )
187    {
188      ptr = &(pNext( *ptr ));
189    }
190    else
191    {
192      pDelete( ptr );
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 )
205{
206  poly    *ptr = &f;
207
208  while( (*ptr)!=(poly)NULL )
209  {
210    if( !isMultiple( Z,*ptr ) )
211    {
212      ptr = &(pNext( *ptr ));
213    }
214    else
215    {
216      pDelete1( ptr );
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 )
233{
234  int     length = INT_MAX;
235  int     result = -1;
236
237  for( int i=0; i<IDELEMS(stdJ); i++ )
238  {
239    if( pComp( stdJ->m[i],m )>=0 && pDivisibleBy( stdJ->m[i],m ) )
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 )
259{
260  for( int i=pVariables; i>0; i-- )
261  {
262    pSetExp( m,i,r[i-1] );
263  }
264  pSetm( m );
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( void )
273{
274  int b = rBlocks( currRing );
275
276  if( b==3 &&
277      ( currRing->order[0] == ringorder_ds ||
278        currRing->order[0] == ringorder_Ds ||
279        currRing->order[0] == ringorder_ws ||
280        currRing->order[0] == ringorder_Ws ) )
281  {
282    return  TRUE;
283  }
284  else if( b>=3
285  && (( currRing->order[0] ==ringorder_a
286        && currRing->block1[0]==pVariables  )
287    || (currRing->order[0]==ringorder_M
288        && currRing->block1[0]==pVariables*pVariables )))
289  {
290    for( int i=pVariables-1; i>=0; i-- )
291    {
292      if( currRing->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
307static void    computeNF( ideal stdJ,poly hc,poly wc,spectrumPolyList *NF )
308{
309  int         carry,k;
310  multiCnt    C( pVariables,0 );
311  poly        Z = (poly)NULL;
312
313  int         well = isWell( );
314
315  do
316  {
317    poly    m = pOne();
318    setExp( m,C.cnt );
319
320    carry = FALSE;
321
322    k = isLeadMonomial( m,stdJ );
323
324    if( k < 0 )
325    {
326      // ---------------------------
327      //  m  is not a lead monomial
328      // ---------------------------
329
330      NF->insert_node( m,(poly)NULL );
331    }
332    else if( isMultiple( Z,m ) )
333    {
334      // ------------------------------------
335      //  m  is trivially in the ideal  stdJ
336      // ------------------------------------
337
338      pDelete( &m );
339      carry = TRUE;
340    }
341    else if( pComp( m,hc ) < 0 || pComp( m,wc ) < 0 )
342    {
343      // -------------------
344      //  we do not need  m
345      // -------------------
346
347      pDelete( &m );
348      carry = TRUE;
349    }
350    else
351    {
352      // --------------------------
353      //  compute lazy normal form
354      // --------------------------
355
356      poly    multiplicant = pDivide( m,stdJ->m[k] );
357      pGetCoeff( multiplicant ) = nInit(1);
358
359      poly    nf = pMultT( pCopy( stdJ->m[k] ), multiplicant );
360
361      pDelete( &multiplicant );
362
363      nf = normalFormHC( nf,hc );
364
365      if( pNext( nf )==(poly)NULL )
366      {
367        // ----------------------------------
368        //  normal form of  m  is  m  itself
369        // ----------------------------------
370
371        pDelete( &nf );
372        NF->delete_monomial( m );
373        Z = pAdd( Z,m );
374        carry = TRUE;
375      }
376      else
377      {
378        nf = normalFormZ( nf,Z );
379
380        if( pNext( nf )==(poly)NULL )
381        {
382          // ----------------------------------
383          //  normal form of  m  is  m  itself
384          // ----------------------------------
385
386          pDelete( &nf );
387          NF->delete_monomial( m );
388          Z = pAdd( Z,m );
389          carry = TRUE;
390        }
391        else
392        {
393          // ------------------------------------
394          //  normal form of  m  is a polynomial
395          // ------------------------------------
396
397          pNorm( nf );
398          if( well==TRUE )
399          {
400            NF->insert_node( m,nf );
401          }
402          else
403          {
404            poly    nfhard = kNF( stdJ,(ideal)NULL,pNext( nf ),0,0 );
405            nfhard = normalFormHC( nfhard,hc );
406            nfhard = normalFormZ ( nfhard,Z );
407
408            if( nfhard==(poly)NULL )
409            {
410              NF->delete_monomial( m );
411              Z = pAdd( Z,m );
412              carry = TRUE;
413            }
414            else
415            {
416              pDelete( &pNext( nf ) );
417              pNext( nf ) = nfhard;
418              NF->insert_node( m,nf );
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!=(poly)NULL && pNext( node->nf )==(poly)NULL )
440      {
441        NF->delete_monomial( node->nf );
442        not_finished = TRUE;
443        node = (spectrumPolyNode*)NULL;
444      }
445      else
446      {
447        node = node->next;
448      }
449    }
450  } while( not_finished );
451
452  pDelete( &Z );
453}
454
455// ----------------------------------------------------------------------------
456//  this is the main spectrum computation function
457// ----------------------------------------------------------------------------
458
459static spectrumState   spectrumCompute( poly h,lists *L,int fast )
460{
461  int i,j;
462
463  #ifdef SPECTRUM_DEBUG
464  #ifdef SPECTRUM_PRINT
465  #ifdef SPECTRUM_IOSTREAM
466    cout << "spectrumCompute\n";
467    if( fast==0 ) cout << "    no optimization" << endl;
468    if( fast==1 ) cout << "    weight optimization" << endl;
469    if( fast==2 ) cout << "    symmetry optimization" << endl;
470  #else
471    fprintf( stdout,"spectrumCompute\n" );
472    if( fast==0 ) fprintf( stdout,"    no optimization\n" );
473    if( fast==1 ) fprintf( stdout,"    weight optimization\n" );
474    if( fast==2 ) fprintf( stdout,"    symmetry optimization\n" );
475  #endif
476  #endif
477  #endif
478
479  // ----------------------
480  //  check if  h  is zero
481  // ----------------------
482
483  if( h==(poly)NULL )
484  {
485    return  spectrumZero;
486  }
487
488  // ----------------------------------
489  //  check if  h  has a constant term
490  // ----------------------------------
491
492  if( hasConstTerm( h ) )
493  {
494    return  spectrumBadPoly;
495  }
496
497  // --------------------------------
498  //  check if  h  has a linear term
499  // --------------------------------
500
501  if( hasLinearTerm( h ) )
502  {
503    *L = (lists)AllocSizeOf( slists );
504    (*L)->Init( 1 );
505    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
506    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
507
508    return  spectrumNoSingularity;
509  }
510
511  // ----------------------------------
512  //  compute the jacobi ideal of  (h)
513  // ----------------------------------
514
515  ideal J = NULL;
516  J = idInit( pVariables,1 );
517
518  #ifdef SPECTRUM_DEBUG
519  #ifdef SPECTRUM_PRINT
520  #ifdef SPECTRUM_IOSTREAM
521    cout << "\n   computing the Jacobi ideal...\n";
522  #else
523    fprintf( stdout,"\n   computing the Jacobi ideal...\n" );
524  #endif
525  #endif
526  #endif
527
528  for( i=0; i<pVariables; i++ )
529  {
530    J->m[i] = pDiff( h,i+1); //j );
531
532    #ifdef SPECTRUM_DEBUG
533    #ifdef SPECTRUM_PRINT
534    #ifdef SPECTRUM_IOSTREAM
535      cout << "        ";
536    #else
537      fprintf( stdout,"        " );
538    #endif
539      pWrite( J->m[i] );
540    #endif
541    #endif
542  }
543
544  // --------------------------------------------
545  //  compute a standard basis  stdJ  of  jac(h)
546  // --------------------------------------------
547
548  #ifdef SPECTRUM_DEBUG
549  #ifdef SPECTRUM_PRINT
550  #ifdef SPECTRUM_IOSTREAM
551    cout << endl;
552    cout << "    computing a standard basis..." << endl;
553  #else
554    fprintf( stdout,"\n" );
555    fprintf( stdout,"    computing a standard basis...\n" );
556  #endif
557  #endif
558  #endif
559
560  ideal stdJ = kStd(J,currQuotient,isNotHomog,NULL);
561  idSkipZeroes( stdJ );
562
563  #ifdef SPECTRUM_DEBUG
564  #ifdef SPECTRUM_PRINT
565    for( i=0; i<IDELEMS(stdJ); i++ )
566    {
567      #ifdef SPECTRUM_IOSTREAM
568        cout << "        ";
569      #else
570        fprintf( stdout,"        " );
571      #endif
572
573      pWrite( stdJ->m[i] );
574    }
575  #endif
576  #endif
577
578  idDelete( &J );
579
580  // ------------------------------------------
581  //  check if the  h  has a singularity
582  // ------------------------------------------
583
584  if( hasOne( stdJ ) )
585  {
586    // -------------------------------
587    //  h is smooth in the origin
588    //  return only the Milnor number
589    // -------------------------------
590
591    *L = (lists)AllocSizeOf( slists );
592    (*L)->Init( 1 );
593    (*L)->m[0].rtyp = INT_CMD;    //  milnor number
594    /* (*L)->m[0].data = (void*)0;a  -- done by Init */
595
596    return  spectrumNoSingularity;
597  }
598
599  // ------------------------------------------
600  //  check if the singularity  h  is isolated
601  // ------------------------------------------
602
603  for( i=pVariables; i>0; i-- )
604  {
605    if( hasAxis( stdJ,i )==FALSE )
606    {
607      return  spectrumNotIsolated;
608    }
609  }
610
611  // ------------------------------------------
612  //  compute the highest corner  hc  of  stdJ
613  // ------------------------------------------
614
615  #ifdef SPECTRUM_DEBUG
616  #ifdef SPECTRUM_PRINT
617  #ifdef SPECTRUM_IOSTREAM
618    cout << "\n    computing the highest corner...\n";
619  #else
620    fprintf( stdout,"\n    computing the highest corner...\n" );
621  #endif
622  #endif
623  #endif
624
625  poly hc = (poly)NULL;
626
627  scComputeHC( stdJ,0,hc );
628
629  if( hc!=(poly)NULL )
630  {
631    pGetCoeff(hc) = nInit(1);
632
633    for( i=pVariables; i>0; i-- )
634    {
635      if( pGetExp( hc,i )>0 ) pDecrExp( hc,i );
636    }
637    pSetm( hc );
638  }
639  else
640  {
641    return  spectrumNoHC;
642  }
643
644  #ifdef SPECTRUM_DEBUG
645  #ifdef SPECTRUM_PRINT
646  #ifdef SPECTRUM_IOSTREAM
647    cout << "       ";
648  #else
649    fprintf( stdout,"       " );
650  #endif
651    pWrite( hc );
652  #endif
653  #endif
654
655  // ----------------------------------------
656  //  compute the Newton polygon  nph  of  h
657  // ----------------------------------------
658
659  #ifdef SPECTRUM_DEBUG
660  #ifdef SPECTRUM_PRINT
661  #ifdef SPECTRUM_IOSTREAM
662    cout << "\n    computing the newton polygon...\n";
663  #else
664    fprintf( stdout,"\n    computing the newton polygon...\n" );
665  #endif
666  #endif
667  #endif
668
669  newtonPolygon nph( h );
670
671  #ifdef SPECTRUM_DEBUG
672  #ifdef SPECTRUM_PRINT
673    cout << nph;
674  #endif
675  #endif
676
677  // -----------------------------------------------
678  //  compute the weight corner  wc  of  (stdj,nph)
679  // -----------------------------------------------
680
681  #ifdef SPECTRUM_DEBUG
682  #ifdef SPECTRUM_PRINT
683  #ifdef SPECTRUM_IOSTREAM
684    cout << "\n    computing the weight corner...\n";
685  #else
686    fprintf( stdout,"\n    computing the weight corner...\n" );
687  #endif
688  #endif
689  #endif
690
691  poly    wc = ( fast==0 ? pCopy( hc ) :
692               ( fast==1 ? computeWC( nph,(Rational)pVariables ) :
693              /* fast==2 */computeWC( nph,((Rational)pVariables)/(Rational)2 ) ) );
694
695  #ifdef SPECTRUM_DEBUG
696  #ifdef SPECTRUM_PRINT
697  #ifdef SPECTRUM_IOSTREAM
698    cout << "        ";
699  #else
700    fprintf( stdout,"        " );
701  #endif
702    pWrite( wc );
703  #endif
704  #endif
705
706  // -------------
707  //  compute  NF
708  // -------------
709
710  #ifdef SPECTRUM_DEBUG
711  #ifdef SPECTRUM_PRINT
712  #ifdef SPECTRUM_IOSTREAM
713    cout << "\n    computing NF...\n" << endl;
714  #else
715    fprintf( stdout,"\n    computing NF...\n" );
716  #endif
717  #endif
718  #endif
719
720  spectrumPolyList NF( &nph );
721
722  computeNF( stdJ,hc,wc,&NF );
723
724  #ifdef SPECTRUM_DEBUG
725  #ifdef SPECTRUM_PRINT
726    cout << NF;
727  #ifdef SPECTRUM_IOSTREAM
728    cout << endl;
729  #else
730    fprintf( stdout,"\n" );
731  #endif
732  #endif
733  #endif
734
735  // ----------------------------
736  //  compute the spectrum of  h
737  // ----------------------------
738
739  return  NF.spectrum( L,fast );
740}
741
742// ----------------------------------------------------------------------------
743//  check if  currRing is local
744// ----------------------------------------------------------------------------
745
746static BOOLEAN ringIsLocal( void )
747{
748  poly    m   = pOne();
749  poly    one = pOne();
750
751  for( int i=pVariables; i>0; i-- )
752  {
753    pSetExp( m,i,1 );
754    pSetm( m );
755
756    if( pComp( m,one )>0 )
757    {
758      return  FALSE;
759    }
760    pSetExp( m,i,0 );
761  }
762
763  pDelete( &m );
764  pDelete( &one );
765
766  return  TRUE;
767}
768
769// ----------------------------------------------------------------------------
770// print error message corresponding to spectrumState state:
771// ----------------------------------------------------------------------------
772static void spectrumPrintError(spectrumState state)
773{
774  switch( state )
775  {
776    case spectrumZero:
777      WerrorS( "polynomial is zero" );
778      break;
779    case spectrumBadPoly:
780      WerrorS( "polynomial has constant term" );
781      break;
782    case spectrumNoSingularity:
783      WerrorS( "not a singularity" );
784      break;
785    case spectrumNotIsolated:
786      WerrorS( "the singularity is not isolated" );
787      break;
788    case spectrumNoHC:
789      WerrorS( "highest corner cannot be computed" );
790      break;
791    case spectrumDegenerate:
792      WerrorS( "principal part is degenerate" );
793      break;
794    case spectrumOK:
795      break;
796
797    default:
798      WerrorS( "unknown error occurred" );
799      break;
800  }
801}
802// ----------------------------------------------------------------------------
803//  this procedure is called from the interpreter
804// ----------------------------------------------------------------------------
805//  first  = polynomial
806//  result = list of spectrum numbers
807// ----------------------------------------------------------------------------
808
809BOOLEAN spectrumProc( leftv result,leftv first )
810{
811  spectrumState state = spectrumOK;
812
813  // -------------------
814  //  check consistency
815  // -------------------
816
817  //  check for a local ring
818
819  if( !ringIsLocal( ) )
820  {
821    WerrorS( "only works for local orderings" );
822    state = spectrumWrongRing;
823  }
824
825  //  no quotient rings are allowed
826
827  else if( currRing->qideal != NULL )
828  {
829    WerrorS( "does not work in quotient rings" );
830    state = spectrumWrongRing;
831  }
832  else
833  {
834    lists   L    = (lists)NULL;
835    int     flag = 1; // weight corner optimization is safe
836
837    state = spectrumCompute( (poly)first->Data( ),&L,flag );
838
839    if( state==spectrumOK )
840    {
841      result->rtyp = LIST_CMD;
842      result->data = (char*)L;
843    }
844    else
845    {
846      spectrumPrintError(state);
847    }
848  }
849
850  return  (state!=spectrumOK);
851}
852
853// ----------------------------------------------------------------------------
854//  this procedure is called from the interpreter
855// ----------------------------------------------------------------------------
856//  first  = polynomial
857//  result = list of spectrum numbers
858// ----------------------------------------------------------------------------
859
860BOOLEAN spectrumfProc( leftv result,leftv first )
861{
862  spectrumState state = spectrumOK;
863
864  // -------------------
865  //  check consistency
866  // -------------------
867
868  //  check for a local polynomial ring
869
870  if( currRing->OrdSgn != -1 )
871  // ?? HS: the test above is also true for k[x][[y]], k[[x]][y]
872  // or should we use:
873  //if( !ringIsLocal( ) )
874  {
875    WerrorS( "only works for local orderings" );
876    state = spectrumWrongRing;
877  }
878  else if( currRing->qideal != NULL )
879  {
880    WerrorS( "does not work in quotient rings" );
881    state = spectrumWrongRing;
882  }
883  else
884  {
885    lists   L    = (lists)NULL;
886    int     flag = 2; // symmetric optimization
887
888    state = spectrumCompute( (poly)first->Data( ),&L,flag );
889
890    if( state==spectrumOK )
891    {
892      result->rtyp = LIST_CMD;
893      result->data = (char*)L;
894    }
895    else
896    {
897      spectrumPrintError(state);
898    }
899  }
900
901  return  (state!=spectrumOK);
902}
903
904// ----------------------------------------------------------------------------
905//  check if a list is a spectrum
906//  check for:
907//      list has 6 elements
908//      1st element is int (mu=Milnor number)
909//      2nd element is int (pg=geometrical genus)
910//      3rd element is int (n =number of different spectrum numbers)
911//      4th element is intvec (num=numerators)
912//      5th element is intvec (den=denomiantors)
913//      6th element is intvec (mul=multiplicities)
914//      exactly n numerators
915//      exactly n denominators
916//      exactly n multiplicities
917//      mu>0
918//      pg>=0
919//      n>0
920//      num>0
921//      den>0
922//      mul>0
923//      symmetriy with respect to numberofvariables/2
924//      monotony
925//      mu = sum of all multiplicities
926//      pg = sum of all multiplicities where num/den<=1
927// ----------------------------------------------------------------------------
928
929semicState  list_is_spectrum( lists l )
930{
931    // -------------------
932    //  check list length
933    // -------------------
934
935    if( l->nr < 5 )
936    {
937        return  semicListTooShort;
938    }
939    else if( l->nr > 5 )
940    {
941        return  semicListTooLong;
942    }
943
944    // -------------
945    //  check types
946    // -------------
947
948    if( l->m[0].rtyp != INT_CMD )
949    {
950        return  semicListFirstElementWrongType;
951    }
952    else if( l->m[1].rtyp != INT_CMD )
953    {
954        return  semicListSecondElementWrongType;
955    }
956    else if( l->m[2].rtyp != INT_CMD )
957    {
958        return  semicListThirdElementWrongType;
959    }
960    else if( l->m[3].rtyp != INTVEC_CMD )
961    {
962        return  semicListFourthElementWrongType;
963    }
964    else if( l->m[4].rtyp != INTVEC_CMD )
965    {
966        return  semicListFifthElementWrongType;
967    }
968    else if( l->m[5].rtyp != INTVEC_CMD )
969    {
970        return  semicListSixthElementWrongType;
971    }
972
973    // -------------------------
974    //  check number of entries
975    // -------------------------
976
977    int     mu = (int)(l->m[0].Data( ));
978    int     pg = (int)(l->m[1].Data( ));
979    int     n  = (int)(l->m[2].Data( ));
980
981    if( n <= 0 )
982    {
983        return  semicListNNegative;
984    }
985
986    intvec  *num = (intvec*)l->m[3].Data( );
987    intvec  *den = (intvec*)l->m[4].Data( );
988    intvec  *mul = (intvec*)l->m[5].Data( );
989
990    if( n != num->length( ) )
991    {
992        return  semicListWrongNumberOfNumerators;
993    }
994    else if( n != den->length( ) )
995    {
996        return  semicListWrongNumberOfDenominators;
997    }
998    else if( n != mul->length( ) )
999    {
1000        return  semicListWrongNumberOfMultiplicities;
1001    }
1002
1003    // --------
1004    //  values
1005    // --------
1006
1007    if( mu <= 0 )
1008    {
1009        return  semicListMuNegative;
1010    }
1011    if( pg < 0 )
1012    {
1013        return  semicListPgNegative;
1014    }
1015
1016    int i;
1017
1018    for( i=0; i<n; i++ )
1019    {
1020        if( (*num)[i] <= 0 )
1021        {
1022            return  semicListNumNegative;
1023        }
1024        if( (*den)[i] <= 0 )
1025        {
1026            return  semicListDenNegative;
1027        }
1028        if( (*mul)[i] <= 0 )
1029        {
1030            return  semicListMulNegative;
1031        }
1032    }
1033
1034    // ----------------
1035    //  check symmetry
1036    // ----------------
1037
1038    int     j;
1039
1040    for( i=0, j=n-1; i<=j; i++,j-- )
1041    {
1042        if( (*num)[i] != pVariables*((*den)[i]) - (*num)[j] ||
1043            (*den)[i] != (*den)[j] ||
1044            (*mul)[i] != (*mul)[j] )
1045        {
1046            return  semicListNotSymmetric;
1047        }
1048    }
1049
1050    // ----------------
1051    //  check monotony
1052    // ----------------
1053
1054    for( i=0, j=1; i<n/2; i++,j++ )
1055    {
1056        if( (*num)[i]*(*den)[j] >= (*num)[j]*(*den)[i] )
1057        {
1058            return  semicListNotMonotonous;
1059        }
1060    }
1061
1062    // ---------------------
1063    //  check Milnor number
1064    // ---------------------
1065
1066    for( mu=0, i=0; i<n; i++ )
1067    {
1068        mu += (*mul)[i];
1069    }
1070
1071    if( mu != (int)(l->m[0].Data( )) )
1072    {
1073        return  semicListMilnorWrong;
1074    }
1075
1076    // -------------------------
1077    //  check geometrical genus
1078    // -------------------------
1079
1080    for( pg=0, i=0; i<n; i++ )
1081    {
1082        if( (*num)[i]<=(*den)[i] )
1083        {
1084            pg += (*mul)[i];
1085        }
1086    }
1087
1088    if( pg != (int)(l->m[1].Data( )) )
1089    {
1090        return  semicListPGWrong;
1091    }
1092
1093    return  semicOK;
1094}
1095
1096// ----------------------------------------------------------------------------
1097//  print out an error message for a spectrum list
1098// ----------------------------------------------------------------------------
1099
1100void    list_error( semicState state )
1101{
1102    switch( state )
1103    {
1104        case semicListTooShort:
1105            WerrorS( "the list is too short" );
1106            break;
1107        case semicListTooLong:
1108            WerrorS( "the list is too long" );
1109            break;
1110
1111        case semicListFirstElementWrongType:
1112            WerrorS( "first element of the list should be int" );
1113            break;
1114        case semicListSecondElementWrongType:
1115            WerrorS( "second element of the list should be int" );
1116            break;
1117        case semicListThirdElementWrongType:
1118            WerrorS( "third element of the list should be int" );
1119            break;
1120        case semicListFourthElementWrongType:
1121            WerrorS( "fourth element of the list should be intvec" );
1122            break;
1123        case semicListFifthElementWrongType:
1124            WerrorS( "fifth element of the list should be intvec" );
1125            break;
1126        case semicListSixthElementWrongType:
1127            WerrorS( "sixth element of the list should be intvec" );
1128            break;
1129
1130        case semicListNNegative:
1131            WerrorS( "first element of the list should be positive" );
1132            break;
1133        case semicListWrongNumberOfNumerators:
1134            WerrorS( "wrong number of numerators" );
1135            break;
1136        case semicListWrongNumberOfDenominators:
1137            WerrorS( "wrong number of denominators" );
1138            break;
1139        case semicListWrongNumberOfMultiplicities:
1140            WerrorS( "wrong number of multiplicities" );
1141            break;
1142
1143        case semicListMuNegative:
1144            WerrorS( "the Milnor number should be positive" );
1145            break;
1146        case semicListPgNegative:
1147            WerrorS( "the geometrical genus should be nonnegative" );
1148            break;
1149        case semicListNumNegative:
1150            WerrorS( "all numerators should be positive" );
1151            break;
1152        case semicListDenNegative:
1153            WerrorS( "all denominators should be positive" );
1154            break;
1155        case semicListMulNegative:
1156            WerrorS( "all multiplicities should be positive" );
1157            break;
1158
1159        case semicListNotSymmetric:
1160            WerrorS( "it is not symmetric" );
1161            break;
1162        case semicListNotMonotonous:
1163            WerrorS( "it is not monotonous" );
1164            break;
1165
1166        case semicListMilnorWrong:
1167            WerrorS( "the Milnor number is wrong" );
1168            break;
1169        case semicListPGWrong:
1170            WerrorS( "the geometrical genus is wrong" );
1171            break;
1172
1173        default:
1174            WerrorS( "unspecific error" );
1175            break;
1176    }
1177}
1178
1179// ----------------------------------------------------------------------------
1180//  this procedure is called from the interpreter
1181// ----------------------------------------------------------------------------
1182//  first  = list of spectrum numbers
1183//  second = list of spectrum numbers
1184//  result = sum of the two lists
1185// ----------------------------------------------------------------------------
1186
1187BOOLEAN spaddProc( leftv result,leftv first,leftv second )
1188{
1189    semicState  state;
1190
1191    // -----------------
1192    //  check arguments
1193    // -----------------
1194
1195    lists l1 = (lists)first->Data( );
1196    lists l2 = (lists)second->Data( );
1197
1198    if( (state=list_is_spectrum( l1 )) != semicOK )
1199    {
1200        WerrorS( "first argument is not a spectrum:" );
1201        list_error( state );
1202    }
1203    else if( (state=list_is_spectrum( l2 )) != semicOK )
1204    {
1205        WerrorS( "second argument is not a spectrum:" );
1206        list_error( state );
1207    }
1208    else
1209    {
1210        spectrum s1( l1 );
1211        spectrum s2( l2 );
1212        spectrum sum( s1+s2 );
1213
1214        result->rtyp = LIST_CMD;
1215        result->data = (char*)(sum.thelist( ));
1216    }
1217
1218    return  (state!=semicOK);
1219}
1220
1221// ----------------------------------------------------------------------------
1222//  this procedure is called from the interpreter
1223// ----------------------------------------------------------------------------
1224//  first  = list of spectrum numbers
1225//  second = integer
1226//  result = the multiple of the first list by the second factor
1227// ----------------------------------------------------------------------------
1228
1229BOOLEAN spmulProc( leftv result,leftv first,leftv second )
1230{
1231    semicState  state;
1232
1233    // -----------------
1234    //  check arguments
1235    // -----------------
1236
1237    lists   l = (lists)first->Data( );
1238    int     k = (int)second->Data( );
1239
1240    if( (state=list_is_spectrum( l ))!=semicOK )
1241    {
1242        WerrorS( "first argument is not a spectrum" );
1243        list_error( state );
1244    }
1245    else if( k < 0 )
1246    {
1247        WerrorS( "second argument should be positive" );
1248        state = semicMulNegative;
1249    }
1250    else
1251    {
1252        spectrum s( l );
1253        spectrum product( k*s );
1254
1255        result->rtyp = LIST_CMD;
1256        result->data = (char*)product.thelist( );
1257    }
1258
1259    return  (state!=semicOK);
1260}
1261
1262BOOLEAN    spectrumOp3  ( leftv res, leftv u, leftv v, leftv w )
1263{
1264  char v_str=(char *)v->Data();
1265  if (strcmp(v_str, "+")==0)
1266    return spaddProc(res,u,w);
1267  else if (strcmp(v_str, "*")==0)
1268    return spmulProc(res,u,w);
1269  Werror("unknown operation '%s' for spectrum",v_str);
1270  return TRUE;
1271}
1272
1273// ----------------------------------------------------------------------------
1274//  this procedure is called from the interpreter
1275// ----------------------------------------------------------------------------
1276//  first  = list of spectrum numbers
1277//  second = list of spectrum numbers
1278//  result = semicontinuity index
1279// ----------------------------------------------------------------------------
1280
1281BOOLEAN    semicProc3   ( leftv res,leftv u,leftv v,leftv w )
1282{
1283  semicState  state;
1284  BOOLEAN qh=(((int)w->Data())==1);
1285
1286  // -----------------
1287  //  check arguments
1288  // -----------------
1289
1290  lists l1 = (lists)u->Data( );
1291  lists l2 = (lists)v->Data( );
1292
1293  if( (state=list_is_spectrum( l1 ))!=semicOK )
1294  {
1295    WerrorS( "first argument is not a spectrum" );
1296    list_error( state );
1297  }
1298  else if( (state=list_is_spectrum( l2 ))!=semicOK )
1299  {
1300    WerrorS( "second argument is not a spectrum" );
1301    list_error( state );
1302  }
1303  else
1304  {
1305    spectrum s1( l1 );
1306    spectrum s2( l2 );
1307
1308    res->rtyp = INT_CMD;
1309    if (qh)
1310      res->data = (void*)(s1.mult_spectrumh( s2 ));
1311    else
1312      res->data = (void*)(s1.mult_spectrum( s2 ));
1313  }
1314
1315  // -----------------
1316  //  check status
1317  // -----------------
1318
1319  return  (state!=semicOK);
1320}
1321BOOLEAN    semicProc   ( leftv res,leftv u,leftv v )
1322{
1323  sleftv tmp;
1324  memset(&tmp,0,sizeof(tmp));
1325  tmp.rtyp=INT_CMD;
1326  /* tmp.data = (void *)0;  -- done by memset */
1327
1328  return  semicProc3(res,u,v,&tmp);
1329}
1330#endif /* HAVE_SPECTRUM */
1331// ----------------------------------------------------------------------------
1332//  spectrum.cc
1333//  end of file
1334// ----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.