source: git/Singular/spectrum.cc @ 8cfee1c

spielwiese
Last change on this file since 8cfee1c was db0e5f, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes: minor optim. git-svn-id: file:///usr/local/Singular/svn/trunk@5225 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 31.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 "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( pLmCmp( f,m )>=0 )
117    {
118      if( pLmDivisibleByNoComp( 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 || pCmp( 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( pLmCmp( *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      pDeleteLm(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( pCmp( 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( pCmp( m,hc ) < 0 || pCmp( 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 = pMult_mm( 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)omAllocBin( slists_bin);
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)omAllocBin( slists_bin);
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,currQuotient, 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  BOOLEAN res=TRUE;
751
752  for( int i=pVariables; i>0; i-- )
753  {
754    pSetExp( m,i,1 );
755    pSetm( m );
756
757    if( pCmp( m,one )>0 )
758    {
759      res=FALSE;
760      break;
761    }
762    pSetExp( m,i,0 );
763  }
764
765  pDelete( &m );
766  pDelete( &one );
767
768  return  res;
769}
770
771// ----------------------------------------------------------------------------
772// print error message corresponding to spectrumState state:
773// ----------------------------------------------------------------------------
774static void spectrumPrintError(spectrumState state)
775{
776  switch( state )
777  {
778    case spectrumZero:
779      WerrorS( "polynomial is zero" );
780      break;
781    case spectrumBadPoly:
782      WerrorS( "polynomial has constant term" );
783      break;
784    case spectrumNoSingularity:
785      WerrorS( "not a singularity" );
786      break;
787    case spectrumNotIsolated:
788      WerrorS( "the singularity is not isolated" );
789      break;
790    case spectrumNoHC:
791      WerrorS( "highest corner cannot be computed" );
792      break;
793    case spectrumDegenerate:
794      WerrorS( "principal part is degenerate" );
795      break;
796    case spectrumOK:
797      break;
798
799    default:
800      WerrorS( "unknown error occurred" );
801      break;
802  }
803}
804// ----------------------------------------------------------------------------
805//  this procedure is called from the interpreter
806// ----------------------------------------------------------------------------
807//  first  = polynomial
808//  result = list of spectrum numbers
809// ----------------------------------------------------------------------------
810
811BOOLEAN spectrumProc( leftv result,leftv first )
812{
813  spectrumState state = spectrumOK;
814
815  // -------------------
816  //  check consistency
817  // -------------------
818
819  //  check for a local ring
820
821  if( !ringIsLocal( ) )
822  {
823    WerrorS( "only works for local orderings" );
824    state = spectrumWrongRing;
825  }
826
827  //  no quotient rings are allowed
828
829  else if( currRing->qideal != NULL )
830  {
831    WerrorS( "does not work in quotient rings" );
832    state = spectrumWrongRing;
833  }
834  else
835  {
836    lists   L    = (lists)NULL;
837    int     flag = 1; // weight corner optimization is safe
838
839    state = spectrumCompute( (poly)first->Data( ),&L,flag );
840
841    if( state==spectrumOK )
842    {
843      result->rtyp = LIST_CMD;
844      result->data = (char*)L;
845    }
846    else
847    {
848      spectrumPrintError(state);
849    }
850  }
851
852  return  (state!=spectrumOK);
853}
854
855// ----------------------------------------------------------------------------
856//  this procedure is called from the interpreter
857// ----------------------------------------------------------------------------
858//  first  = polynomial
859//  result = list of spectrum numbers
860// ----------------------------------------------------------------------------
861
862BOOLEAN spectrumfProc( leftv result,leftv first )
863{
864  spectrumState state = spectrumOK;
865
866  // -------------------
867  //  check consistency
868  // -------------------
869
870  //  check for a local polynomial ring
871
872  if( currRing->OrdSgn != -1 )
873  // ?? HS: the test above is also true for k[x][[y]], k[[x]][y]
874  // or should we use:
875  //if( !ringIsLocal( ) )
876  {
877    WerrorS( "only works for local orderings" );
878    state = spectrumWrongRing;
879  }
880  else if( currRing->qideal != NULL )
881  {
882    WerrorS( "does not work in quotient rings" );
883    state = spectrumWrongRing;
884  }
885  else
886  {
887    lists   L    = (lists)NULL;
888    int     flag = 2; // symmetric optimization
889
890    state = spectrumCompute( (poly)first->Data( ),&L,flag );
891
892    if( state==spectrumOK )
893    {
894      result->rtyp = LIST_CMD;
895      result->data = (char*)L;
896    }
897    else
898    {
899      spectrumPrintError(state);
900    }
901  }
902
903  return  (state!=spectrumOK);
904}
905
906// ----------------------------------------------------------------------------
907//  check if a list is a spectrum
908//  check for:
909//      list has 6 elements
910//      1st element is int (mu=Milnor number)
911//      2nd element is int (pg=geometrical genus)
912//      3rd element is int (n =number of different spectrum numbers)
913//      4th element is intvec (num=numerators)
914//      5th element is intvec (den=denomiantors)
915//      6th element is intvec (mul=multiplicities)
916//      exactly n numerators
917//      exactly n denominators
918//      exactly n multiplicities
919//      mu>0
920//      pg>=0
921//      n>0
922//      num>0
923//      den>0
924//      mul>0
925//      symmetriy with respect to numberofvariables/2
926//      monotony
927//      mu = sum of all multiplicities
928//      pg = sum of all multiplicities where num/den<=1
929// ----------------------------------------------------------------------------
930
931semicState  list_is_spectrum( lists l )
932{
933    // -------------------
934    //  check list length
935    // -------------------
936
937    if( l->nr < 5 )
938    {
939        return  semicListTooShort;
940    }
941    else if( l->nr > 5 )
942    {
943        return  semicListTooLong;
944    }
945
946    // -------------
947    //  check types
948    // -------------
949
950    if( l->m[0].rtyp != INT_CMD )
951    {
952        return  semicListFirstElementWrongType;
953    }
954    else if( l->m[1].rtyp != INT_CMD )
955    {
956        return  semicListSecondElementWrongType;
957    }
958    else if( l->m[2].rtyp != INT_CMD )
959    {
960        return  semicListThirdElementWrongType;
961    }
962    else if( l->m[3].rtyp != INTVEC_CMD )
963    {
964        return  semicListFourthElementWrongType;
965    }
966    else if( l->m[4].rtyp != INTVEC_CMD )
967    {
968        return  semicListFifthElementWrongType;
969    }
970    else if( l->m[5].rtyp != INTVEC_CMD )
971    {
972        return  semicListSixthElementWrongType;
973    }
974
975    // -------------------------
976    //  check number of entries
977    // -------------------------
978
979    int     mu = (int)(l->m[0].Data( ));
980    int     pg = (int)(l->m[1].Data( ));
981    int     n  = (int)(l->m[2].Data( ));
982
983    if( n <= 0 )
984    {
985        return  semicListNNegative;
986    }
987
988    intvec  *num = (intvec*)l->m[3].Data( );
989    intvec  *den = (intvec*)l->m[4].Data( );
990    intvec  *mul = (intvec*)l->m[5].Data( );
991
992    if( n != num->length( ) )
993    {
994        return  semicListWrongNumberOfNumerators;
995    }
996    else if( n != den->length( ) )
997    {
998        return  semicListWrongNumberOfDenominators;
999    }
1000    else if( n != mul->length( ) )
1001    {
1002        return  semicListWrongNumberOfMultiplicities;
1003    }
1004
1005    // --------
1006    //  values
1007    // --------
1008
1009    if( mu <= 0 )
1010    {
1011        return  semicListMuNegative;
1012    }
1013    if( pg < 0 )
1014    {
1015        return  semicListPgNegative;
1016    }
1017
1018    int i;
1019
1020    for( i=0; i<n; i++ )
1021    {
1022        if( (*num)[i] <= 0 )
1023        {
1024            return  semicListNumNegative;
1025        }
1026        if( (*den)[i] <= 0 )
1027        {
1028            return  semicListDenNegative;
1029        }
1030        if( (*mul)[i] <= 0 )
1031        {
1032            return  semicListMulNegative;
1033        }
1034    }
1035
1036    // ----------------
1037    //  check symmetry
1038    // ----------------
1039
1040    int     j;
1041
1042    for( i=0, j=n-1; i<=j; i++,j-- )
1043    {
1044        if( (*num)[i] != pVariables*((*den)[i]) - (*num)[j] ||
1045            (*den)[i] != (*den)[j] ||
1046            (*mul)[i] != (*mul)[j] )
1047        {
1048            return  semicListNotSymmetric;
1049        }
1050    }
1051
1052    // ----------------
1053    //  check monotony
1054    // ----------------
1055
1056    for( i=0, j=1; i<n/2; i++,j++ )
1057    {
1058        if( (*num)[i]*(*den)[j] >= (*num)[j]*(*den)[i] )
1059        {
1060            return  semicListNotMonotonous;
1061        }
1062    }
1063
1064    // ---------------------
1065    //  check Milnor number
1066    // ---------------------
1067
1068    for( mu=0, i=0; i<n; i++ )
1069    {
1070        mu += (*mul)[i];
1071    }
1072
1073    if( mu != (int)(l->m[0].Data( )) )
1074    {
1075        return  semicListMilnorWrong;
1076    }
1077
1078    // -------------------------
1079    //  check geometrical genus
1080    // -------------------------
1081
1082    for( pg=0, i=0; i<n; i++ )
1083    {
1084        if( (*num)[i]<=(*den)[i] )
1085        {
1086            pg += (*mul)[i];
1087        }
1088    }
1089
1090    if( pg != (int)(l->m[1].Data( )) )
1091    {
1092        return  semicListPGWrong;
1093    }
1094
1095    return  semicOK;
1096}
1097
1098// ----------------------------------------------------------------------------
1099//  print out an error message for a spectrum list
1100// ----------------------------------------------------------------------------
1101
1102void    list_error( semicState state )
1103{
1104    switch( state )
1105    {
1106        case semicListTooShort:
1107            WerrorS( "the list is too short" );
1108            break;
1109        case semicListTooLong:
1110            WerrorS( "the list is too long" );
1111            break;
1112
1113        case semicListFirstElementWrongType:
1114            WerrorS( "first element of the list should be int" );
1115            break;
1116        case semicListSecondElementWrongType:
1117            WerrorS( "second element of the list should be int" );
1118            break;
1119        case semicListThirdElementWrongType:
1120            WerrorS( "third element of the list should be int" );
1121            break;
1122        case semicListFourthElementWrongType:
1123            WerrorS( "fourth element of the list should be intvec" );
1124            break;
1125        case semicListFifthElementWrongType:
1126            WerrorS( "fifth element of the list should be intvec" );
1127            break;
1128        case semicListSixthElementWrongType:
1129            WerrorS( "sixth element of the list should be intvec" );
1130            break;
1131
1132        case semicListNNegative:
1133            WerrorS( "first element of the list should be positive" );
1134            break;
1135        case semicListWrongNumberOfNumerators:
1136            WerrorS( "wrong number of numerators" );
1137            break;
1138        case semicListWrongNumberOfDenominators:
1139            WerrorS( "wrong number of denominators" );
1140            break;
1141        case semicListWrongNumberOfMultiplicities:
1142            WerrorS( "wrong number of multiplicities" );
1143            break;
1144
1145        case semicListMuNegative:
1146            WerrorS( "the Milnor number should be positive" );
1147            break;
1148        case semicListPgNegative:
1149            WerrorS( "the geometrical genus should be nonnegative" );
1150            break;
1151        case semicListNumNegative:
1152            WerrorS( "all numerators should be positive" );
1153            break;
1154        case semicListDenNegative:
1155            WerrorS( "all denominators should be positive" );
1156            break;
1157        case semicListMulNegative:
1158            WerrorS( "all multiplicities should be positive" );
1159            break;
1160
1161        case semicListNotSymmetric:
1162            WerrorS( "it is not symmetric" );
1163            break;
1164        case semicListNotMonotonous:
1165            WerrorS( "it is not monotonous" );
1166            break;
1167
1168        case semicListMilnorWrong:
1169            WerrorS( "the Milnor number is wrong" );
1170            break;
1171        case semicListPGWrong:
1172            WerrorS( "the geometrical genus is wrong" );
1173            break;
1174
1175        default:
1176            WerrorS( "unspecific error" );
1177            break;
1178    }
1179}
1180
1181// ----------------------------------------------------------------------------
1182//  this procedure is called from the interpreter
1183// ----------------------------------------------------------------------------
1184//  first  = list of spectrum numbers
1185//  second = list of spectrum numbers
1186//  result = sum of the two lists
1187// ----------------------------------------------------------------------------
1188
1189BOOLEAN spaddProc( leftv result,leftv first,leftv second )
1190{
1191    semicState  state;
1192
1193    // -----------------
1194    //  check arguments
1195    // -----------------
1196
1197    lists l1 = (lists)first->Data( );
1198    lists l2 = (lists)second->Data( );
1199
1200    if( (state=list_is_spectrum( l1 )) != semicOK )
1201    {
1202        WerrorS( "first argument is not a spectrum:" );
1203        list_error( state );
1204    }
1205    else if( (state=list_is_spectrum( l2 )) != semicOK )
1206    {
1207        WerrorS( "second argument is not a spectrum:" );
1208        list_error( state );
1209    }
1210    else
1211    {
1212        spectrum s1( l1 );
1213        spectrum s2( l2 );
1214        spectrum sum( s1+s2 );
1215
1216        result->rtyp = LIST_CMD;
1217        result->data = (char*)(sum.thelist( ));
1218    }
1219
1220    return  (state!=semicOK);
1221}
1222
1223// ----------------------------------------------------------------------------
1224//  this procedure is called from the interpreter
1225// ----------------------------------------------------------------------------
1226//  first  = list of spectrum numbers
1227//  second = integer
1228//  result = the multiple of the first list by the second factor
1229// ----------------------------------------------------------------------------
1230
1231BOOLEAN spmulProc( leftv result,leftv first,leftv second )
1232{
1233    semicState  state;
1234
1235    // -----------------
1236    //  check arguments
1237    // -----------------
1238
1239    lists   l = (lists)first->Data( );
1240    int     k = (int)second->Data( );
1241
1242    if( (state=list_is_spectrum( l ))!=semicOK )
1243    {
1244        WerrorS( "first argument is not a spectrum" );
1245        list_error( state );
1246    }
1247    else if( k < 0 )
1248    {
1249        WerrorS( "second argument should be positive" );
1250        state = semicMulNegative;
1251    }
1252    else
1253    {
1254        spectrum s( l );
1255        spectrum product( k*s );
1256
1257        result->rtyp = LIST_CMD;
1258        result->data = (char*)product.thelist( );
1259    }
1260
1261    return  (state!=semicOK);
1262}
1263
1264// ----------------------------------------------------------------------------
1265//  this procedure is called from the interpreter
1266// ----------------------------------------------------------------------------
1267//  first  = list of spectrum numbers
1268//  second = list of spectrum numbers
1269//  result = semicontinuity index
1270// ----------------------------------------------------------------------------
1271
1272BOOLEAN    semicProc3   ( leftv res,leftv u,leftv v,leftv w )
1273{
1274  semicState  state;
1275  BOOLEAN qh=(((int)w->Data())==1);
1276
1277  // -----------------
1278  //  check arguments
1279  // -----------------
1280
1281  lists l1 = (lists)u->Data( );
1282  lists l2 = (lists)v->Data( );
1283
1284  if( (state=list_is_spectrum( l1 ))!=semicOK )
1285  {
1286    WerrorS( "first argument is not a spectrum" );
1287    list_error( state );
1288  }
1289  else if( (state=list_is_spectrum( l2 ))!=semicOK )
1290  {
1291    WerrorS( "second argument is not a spectrum" );
1292    list_error( state );
1293  }
1294  else
1295  {
1296    spectrum s1( l1 );
1297    spectrum s2( l2 );
1298
1299    res->rtyp = INT_CMD;
1300    if (qh)
1301      res->data = (void*)(s1.mult_spectrumh( s2 ));
1302    else
1303      res->data = (void*)(s1.mult_spectrum( s2 ));
1304  }
1305
1306  // -----------------
1307  //  check status
1308  // -----------------
1309
1310  return  (state!=semicOK);
1311}
1312BOOLEAN    semicProc   ( leftv res,leftv u,leftv v )
1313{
1314  sleftv tmp;
1315  memset(&tmp,0,sizeof(tmp));
1316  tmp.rtyp=INT_CMD;
1317  /* tmp.data = (void *)0;  -- done by memset */
1318
1319  return  semicProc3(res,u,v,&tmp);
1320}
1321#endif /* HAVE_SPECTRUM */
1322// ----------------------------------------------------------------------------
1323//  spectrum.cc
1324//  end of file
1325// ----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.