Changeset 86c1f0 in git for Singular/LIB/gaussman.lib
 Timestamp:
 Aug 9, 2001, 10:39:54 AM (22 years ago)
 Branches:
 (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
 Children:
 85140660b24ae4a0f2f58450c04b514020b5d8f1
 Parents:
 1ac173caaa52482213be39480214765a2e70eb5f
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/gaussman.lib
r1ac173 r86c1f0 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: gaussman.lib,v 1.4 8 20010801 13:21:36mschulze Exp $";2 version="$Id: gaussman.lib,v 1.49 20010809 08:39:54 mschulze Exp $"; 3 3 category="Singularities"; 4 4 … … 12 12 13 13 PROCEDURES: 14 gmsring( f,s); Brieskorn lattice in the GaussManin system14 gmsring(t,s); Brieskorn lattice in GaussManin system of t 15 15 gmsnf(p,K[,Kmax]); GaussManin system normal form 16 16 gmscoeffs(p,K[,Kmax]); GaussManin system basis representation 17 monodromy( f[,opt]); monodromy matrix or spectrum of monodromy of f18 vfilt(f[,opt]); Vfiltration on H''/H' or spectrum of f19 endfilt(poly f,list V); endomorphism filtration of filtration V20 spectrum(f); spectrum of f21 sppairs(f[,opt]); spectral pairs or spectrum of f17 monodromy(t[,opt]); Jordan data or eigenvalues of monodromy of t 18 spectrum(t); spectrum of t 19 sppairs(t[,opt]); spectral pairs or spectrum of t 20 vfilt(t[,opt]); Vfiltration on H''/H' or spectrum of t 21 endfilt(t,V); endomorphism filtration of Vfiltration V 22 22 spgen(a); generate spectrum defined by a 23 23 sppgen(a,w); generate spectral pairs defined by a and w … … 106 106 107 107 proc gmsring(poly t,string s) 108 "USAGE: gmsring( f,s); poly f, string s;108 "USAGE: gmsring(t,s); poly t, string s; 109 109 ASSUME: basering with characteristic 0 and local degree ordering, 110 fwith isolated singularity at 0110 t with isolated singularity at 0 111 111 RETURN: 112 112 @format 113 ring G: 114 G=C{{s}}*basering, 115 s is the inverse of the GaussManin connection of f, 116 G contains: 117 poly gmspoly: image of f 118 ideal gmsjacob: image of Jacobian ideal 119 ideal gmsstd: image of standard basis of Jacobian ideal 120 matrix gmsmatrix: matrix(gmsjacob)*gmsmatrix=matrix(gmsstd) 121 ideal gmsbasis: image of monomial vector space basis of Jacobian algebra 122 int gmsmaxweight: maximal weight of variables of basering 123 G projects to H''=C{{s}}*gmsbasis 113 ring G: C{{s}}*basering, 114 poly gmspoly: image of t 115 ideal gmsjacob: image of Jacobian ideal 116 ideal gmsstd: image of standard basis of Jacobian ideal 117 matrix gmsmatrix: matrix(gmsjacob)*gmsmatrix=matrix(gmsstd) 118 ideal gmsbasis: image of monomial vector space basis of Jacobian algebra 119 int gmsmaxweight: maximal weight of variables of basering 124 120 @end format 125 121 NOTE: do not modify gms variables if you want to use gms procedures … … 203 199 { "EXAMPLE:"; echo=2; 204 200 ring R=0,(x,y),ds; 205 poly f=x5+x2y2+y5;206 def G=gmsring( f,"s");201 poly t=x5+x2y2+y5; 202 def G=gmsring(t,"s"); 207 203 setring(G); 208 204 gmspoly; … … 210 206 print(gmsstd); 211 207 print(gmsmatrix); 208 print(gmsbasis); 212 209 gmsmaxweight; 213 210 } … … 216 213 proc gmsnf(ideal p,int K,list #) 217 214 "USAGE: gmsnf(p,K[,Kmax]); poly p, int K[, int Kmax]; 218 ASSUME: basering was constructed by gms, K<=Kmax215 ASSUME: basering constructed by gmsring, K<=Kmax 219 216 RETURN: 220 217 @format … … 305 302 { "EXAMPLE:"; echo=2; 306 303 ring R=0,(x,y),ds; 307 poly f=x5+x2y2+y5;308 def G=gmsring( f,"s");304 poly t=x5+x2y2+y5; 305 def G=gmsring(t,"s"); 309 306 setring(G); 310 307 list l0=gmsnf(gmspoly,0); … … 319 316 proc gmscoeffs(ideal p,int K,list #) 320 317 "USAGE: gmscoeffs(p,K[,Kmax]); poly p, int K[, int Kmax]; 321 ASSUME: basering was constructed by gms, K<=Kmax318 ASSUME: basering constructed by gmsring, K<=Kmax 322 319 RETURN: 323 320 @format … … 345 342 { "EXAMPLE:"; echo=2; 346 343 ring R=0,(x,y),ds; 347 poly f=x5+x2y2+y5;348 def G=gmsring( f,"s");344 poly t=x5+x2y2+y5; 345 def G=gmsring(t,"s"); 349 346 setring(G); 350 347 list l0=gmscoeffs(gmspoly,0); … … 359 356 static proc maxintdif(ideal e) 360 357 { 361 int i,j, id;362 int mid=0;358 int i,j,d; 359 int d0=0; 363 360 for(i=ncols(e);i>=1;i) 364 361 { 365 362 for(j=i1;j>=1;j) 366 363 { 367 id=int(e[i]e[j]);368 if( id<0)369 { 370 id=id;371 } 372 if( id>mid)373 { 374 mid=id;375 } 376 } 377 } 378 return( mid);364 d=int(e[i]e[j]); 365 if(d<0) 366 { 367 d=d; 368 } 369 if(d>d0) 370 { 371 d0=d; 372 } 373 } 374 } 375 return(d0); 379 376 } 380 377 /////////////////////////////////////////////////////////////////////////////// … … 406 403 /////////////////////////////////////////////////////////////////////////////// 407 404 408 proc monodromy(poly f,list #)409 "USAGE: monodromy( f[,opt]); poly f, int opt405 proc monodromy(poly t,list #) 406 "USAGE: monodromy(t[,opt]); poly t[, int opt] 410 407 ASSUME: basering with characteristic 0 and local degree ordering, 411 fwith isolated singularity at 0408 t with isolated singularity at 0 412 409 RETURN: 413 410 @format … … 435 432 def R=basering; 436 433 int n=nvars(R)1; 437 def G=gmsring( f,"s");434 def G=gmsring(t,"s"); 438 435 setring G; 439 436 440 int mu ,modm=ncols(gmsbasis),maxorddif(gmsbasis);437 int mu=ncols(gmsbasis); 441 438 ideal r=gmspoly*gmsbasis; 442 439 list l; 443 440 matrix U=freemodule(mu); 444 441 matrix A0[mu][mu],C; 445 module H, dH=freemodule(mu),freemodule(mu);442 module H,H1=freemodule(mu),freemodule(mu); 446 443 module H0; 447 444 int k=1; … … 464 461 dbprint(printlevelvoice+2,"// compute saturation of H''"); 465 462 H0=H; 466 dH=jet(module(A0*dH+s^2*diff(matrix(dH),s)),k);467 H=H*s+ dH;463 H1=jet(module(A0*H1+s^2*diff(matrix(H1),s)),k); 464 H=H*s+H1; 468 465 } 469 466 A0=A0k*s; … … 471 468 dbprint(printlevelvoice+2,"// compute basis of saturation of H''"); 472 469 H=std(H0); 473 int modH=maxorddif(H)/deg(s);474 dbprint(printlevelvoice+2,"// k="+string( modH+1));470 int d0=maxorddif(H)/deg(s); 471 dbprint(printlevelvoice+2,"// k="+string(d0+1)); 475 472 dbprint(printlevelvoice+2,"// compute matrix A of t"); 476 473 if(opt==0) 477 474 { 478 l=gmscoeffs(r, modH+1,modH+1);475 l=gmscoeffs(r,d0+1,d0+1); 479 476 } 480 477 else 481 478 { 482 l=gmscoeffs(r, modH+1,modH+1+n);479 l=gmscoeffs(r,d0+1,d0+1+n); 483 480 } 484 481 C,r=l[1..2]; … … 502 499 } 503 500 504 dbprint(printlevelvoice+2,"// compute maximal integer difference of e");505 int mide=maxintdif(e);506 dbprint(printlevelvoice+2,"// mide="+string(mide));507 508 if( mide>0)509 { 510 dbprint(printlevelvoice+2,"// k="+string( modH+1+mide));501 dbprint(printlevelvoice+2,"// compute maximal integer difference d1 of e"); 502 int d1=maxintdif(e); 503 dbprint(printlevelvoice+2,"// d1="+string(d1)); 504 505 if(d1>0) 506 { 507 dbprint(printlevelvoice+2,"// k="+string(d0+d1+1)); 511 508 dbprint(printlevelvoice+2,"// compute matrix A of t"); 512 l=gmscoeffs(r, modH+1+mide,modH+1+mide);509 l=gmscoeffs(r,d0+d1+1,d0+d1+1); 513 510 C,r=l[1..2]; 514 511 A0=A0+C; 515 512 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 516 513 l=division(H*s,A0*H+s^2*diff(matrix(H),s)); 517 A=jet(l[1],l[2], mide);518 519 intvec ide;520 ide[mu]=0;514 A=jet(l[1],l[2],d1); 515 516 intvec d; 517 d[mu]=0; 521 518 int i,j; 522 519 for(i=ncols(e);i>=1;i) … … 525 522 { 526 523 k=int(e[j]e[i]); 527 if(k> ide[j])528 { 529 ide[j]=k;530 } 531 if(k> ide[i])532 { 533 ide[i]=k;524 if(k>d[j]) 525 { 526 d[j]=k; 527 } 528 if(k>d[i]) 529 { 530 d[i]=k; 534 531 } 535 532 } … … 539 536 for(i=m[j];i>=1;i) 540 537 { 541 ide[k]=ide[j];538 d[k]=d[j]; 542 539 k; 543 540 } 544 541 } 545 542 546 while( mide>0)547 { 548 dbprint(printlevelvoice+2,"// transform basis to reduce mideby 1");543 while(d1>0) 544 { 545 dbprint(printlevelvoice+2,"// transform basis to reduce d1 by 1"); 549 546 550 547 A0=jet(A,0); … … 560 557 for(j=mu;j>=1;j) 561 558 { 562 if( ide[i]==0&&ide[j]>0)559 if(d[i]==0&&d[j]>0) 563 560 { 564 561 A[i,j]=A[i,j]/s; … … 566 563 else 567 564 { 568 if( ide[i]>0&&ide[j]==0)565 if(d[i]>0&&d[j]==0) 569 566 { 570 567 A[i,j]=A[i,j]*s; … … 575 572 for(i=mu;i>=1;i) 576 573 { 577 if( ide[i]>0)574 if(d[i]>0) 578 575 { 579 576 A[i,i]=A[i,i]1; 580 ide[i]=ide[i]1;581 } 582 } 583 584 mide;585 dbprint(printlevelvoice+2,"// mide="+string(mide));577 d[i]=d[i]1; 578 } 579 } 580 581 d1; 582 dbprint(printlevelvoice+2,"// d1="+string(d1)); 586 583 } 587 584 … … 601 598 /////////////////////////////////////////////////////////////////////////////// 602 599 603 proc vfilt(poly f,list #)604 "USAGE: vfilt(f[,opt]); poly f, int opt600 proc spectrum(poly t) 601 "USAGE: spectrum(t); poly t 605 602 ASSUME: basering with characteristic 0 and local degree ordering, 606 fwith isolated singularity at 0603 t with isolated singularity at 0 607 604 RETURN: 608 605 @format 609 list V: Vfiltration of f on H''/H' 610 if opt=0 or opt=1: 611 intvec V[1]: numerators of spectral numbers 612 intvec V[2]: denominators of spectral numbers 613 intvec V[3]: 614 int V[3][i]: multiplicity of spectral number V[1][i]/V[2][i] 615 if opt=1: 606 list Sp: spectrum of t 607 ideal Sp[1]: spectral numbers in increasing order 608 intvec Sp[2]: 609 int Sp[2][i]: multiplicity of spectral number Sp[1][i] 610 @end format 611 SEE ALSO: spectrum_lib 612 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; spectrum 613 EXAMPLE: example spnumbers; shows examples 614 " 615 { 616 return(sppairs(t,0)); 617 } 618 example 619 { "EXAMPLE:"; echo=2; 620 ring R=0,(x,y),ds; 621 poly t=x5+x2y2+y5; 622 spprint(spectrum(t)); 623 } 624 /////////////////////////////////////////////////////////////////////////////// 625 626 proc sppairs(poly t,list #) 627 "USAGE: sppairs(t[,opt]); poly t[, int opt] 628 ASSUME: basering with characteristic 0 and local degree ordering, 629 t with isolated singularity at 0 630 RETURN: list l: 631 @format 632 ideal l[1]: spectral numbers in increasing order 633 if opt<=0: 634 intvec l[2]: 635 int l[2][i]: multiplicity of spectral number l[1][i] 636 if opt>=1: 637 intvec l[2]: weight filtration indices in decreasing order 638 intvec l[3]: 639 int l[3][i]: multiplicity of spectral pair (l[1][i],l[2][i]) 640 default: opt=1 641 @end format 642 SEE ALSO: spectrum_lib 643 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 644 spectrum; spectral pairs 645 EXAMPLE: example sppairs; shows examples 646 " 647 { 648 int opt=1; 649 if(size(#)>0) 650 { 651 if(typeof(#[1])=="int") 652 { 653 opt=#[1]; 654 } 655 } 656 657 def R=basering; 658 int n=nvars(R)1; 659 def G=gmsring(t,"s"); 660 setring(G); 661 662 int mu=ncols(gmsbasis); 663 ideal r=gmspoly*gmsbasis; 664 list l; 665 matrix U=freemodule(mu); 666 matrix A0[mu][mu],C; 667 module H0; 668 module H,H1=freemodule(mu),freemodule(mu); 669 int k=1; 670 while(size(reduce(H,std(H0*s)))>0) 671 { 672 k++; 673 dbprint(printlevelvoice+2,"// k="+string(k)); 674 dbprint(printlevelvoice+2,"// compute matrix A of t"); 675 l=gmscoeffs(r,k,mu+n); 676 C,r=l[1..2]; 677 A0=A0+C; 678 679 dbprint(printlevelvoice+2,"// compute saturation of H''"); 680 H0=H; 681 H1=jet(module(A0*H1+s^2*diff(matrix(H1),s)),k); 682 H=H*s+H1; 683 } 684 A0=A0k*s; 685 686 dbprint(printlevelvoice+2,"// compute basis of saturation of H''"); 687 H=std(H0); 688 int d0=maxorddif(H)/deg(s); 689 dbprint(printlevelvoice+2,"// transform H'' to saturation of H''"); 690 l=division(H,freemodule(mu)*s^k); 691 H0=l[1]; 692 693 dbprint(printlevelvoice+2,"// k="+string(d0+1)); 694 dbprint(printlevelvoice+2,"// compute matrix A of t"); 695 l=gmscoeffs(r,d0+1,d0+1+n); 696 C,r=l[1..2]; 697 A0=A0+C; 698 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 699 l=division(H*s,A0*H+s^2*diff(matrix(H),s)); 700 matrix A=jet(l[1],l[2],0); 701 702 int i,j; 703 matrix V; 704 if(opt<=0) 705 { 706 dbprint(printlevelvoice+2,"// transform to Jordan basis"); 707 U=jordanbasis(A)[1]; 708 V=lift(U,freemodule(mu)); 709 A=V*A*U; 710 dbprint(printlevelvoice+2,"// compute normal form of H''"); 711 H0=std(V*H0); 712 713 dbprint(printlevelvoice+2,"// compute spectral numbers"); 714 ideal a; 715 for(i=1;i<=mu;i++) 716 { 717 j=leadexp(H0[i])[nvars(basering)+1]; 718 a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)1; 719 } 720 721 setring(R); 722 return(spgen(imap(G,a))); 723 } 724 725 dbprint(printlevelvoice+2, 726 "// compute eigenvalues e with multiplicities m of A"); 727 l=eigenval(A); 728 def e,m=l[1..2]; 729 dbprint(printlevelvoice+2,"// e="+string(e)); 730 dbprint(printlevelvoice+2,"// m="+string(m)); 731 dbprint(printlevelvoice+2,"// compute maximal integer difference d1 of e"); 732 int d1=maxintdif(e); 733 dbprint(printlevelvoice+2,"// d1="+string(d1)); 734 735 if(d1>0) 736 { 737 dbprint(printlevelvoice+2,"// k="+string(d0+d1+1)); 738 dbprint(printlevelvoice+2,"// compute matrix A of t"); 739 l=gmscoeffs(r,d0+d1+1,d0+d1+1); 740 C,r=l[1..2]; 741 A0=A0+C; 742 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 743 l=division(H*s,A0*H+s^2*diff(matrix(H),s)); 744 A=jet(l[1],l[2],d1); 745 746 intvec d; 747 d[mu]=0; 748 for(i=ncols(e);i>=1;i) 749 { 750 for(j=i1;j>=1;j) 751 { 752 k=int(e[j]e[i]); 753 if(k>d[j]) 754 { 755 d[j]=k; 756 } 757 if(k>d[i]) 758 { 759 d[i]=k; 760 } 761 } 762 } 763 for(j,k=ncols(e),mu;j>=1;j) 764 { 765 for(i=m[j];i>=1;i) 766 { 767 d[k]=d[j]; 768 k; 769 } 770 } 771 772 while(d1>0) 773 { 774 dbprint(printlevelvoice+2,"// transform to separate eigenvalues"); 775 A0=jet(A,0); 776 U=0; 777 for(i=ncols(e);i>=1;i) 778 { 779 U=syz(power(A0e[i],n+1))+U; 780 } 781 V=lift(U,freemodule(mu)); 782 A=V*A*U; 783 H0=V*H0; 784 785 dbprint(printlevelvoice+2,"// transform to reduce d1 by 1"); 786 for(i=mu;i>=1;i) 787 { 788 for(j=mu;j>=1;j) 789 { 790 if(d[i]==0&&d[j]>0) 791 { 792 A[i,j]=A[i,j]/s; 793 } 794 else 795 { 796 if(d[i]>0&&d[j]==0) 797 { 798 A[i,j]=A[i,j]*s; 799 } 800 } 801 } 802 } 803 H0=transpose(H0); 804 for(i=mu;i>=1;i) 805 { 806 if(d[i]>0) 807 { 808 A[i,i]=A[i,i]1; 809 d[i]=d[i]1; 810 H0[i]=H0[i]*s; 811 } 812 } 813 H0=transpose(H0); 814 815 d1; 816 dbprint(printlevelvoice+2,"// d1="+string(d1)); 817 } 818 819 A=jet(A,0); 820 } 821 822 dbprint(printlevelvoice+2,"// transform to Jordan basis"); 823 intvec w0; 824 l=jordanbasis(A); 825 U,w0=l[1..2]; 826 V=lift(U,freemodule(mu)); 827 A0=V*A*U; 828 829 dbprint(printlevelvoice+2,"// compute weight filtration basis"); 830 vector u; 831 i=1; 832 while(i<ncols(A)) 833 { 834 j=i+1; 835 while(j<ncols(A)&&A[i,i]==A[j,j]) 836 { 837 if(w0[i]<w0[j]) 838 { 839 k=w0[i]; 840 w0[i]=w0[j]; 841 w0[i]=k; 842 u=U[i]; 843 U[i]=U[j]; 844 U[j]=u; 845 } 846 j++; 847 } 848 if(j==ncols(A)&&A[i,i]==A[j,j]&&w0[i]<w0[j]) 849 { 850 k=w0[i]; 851 w0[i]=w0[j]; 852 w0[i]=k; 853 u=U[i]; 854 U[i]=U[j]; 855 U[j]=u; 856 } 857 i=j; 858 } 859 860 dbprint(printlevelvoice+2,"// transform to weight filtration basis"); 861 V=lift(U,freemodule(mu)); 862 A=V*A*U; 863 dbprint(printlevelvoice+2,"// compute normal form of H''"); 864 H0=std(V*H0); 865 866 dbprint(printlevelvoice+2,"// compute spectral pairs"); 867 ideal a; 868 intvec w; 869 for(i=1;i<=mu;i++) 870 { 871 j=leadexp(H0[i])[nvars(basering)+1]; 872 a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)1; 873 w[i]=w0[j]+n; 874 } 875 setring(R); 876 return(sppgen(imap(G,a),w)); 877 } 878 example 879 { "EXAMPLE:"; echo=2; 880 ring R=0,(x,y),ds; 881 poly t=x5+x2y2+y5; 882 spprint(sppairs(t)); 883 } 884 /////////////////////////////////////////////////////////////////////////////// 885 886 proc vfilt(poly t,list #) 887 "USAGE: vfilt(t[,opt]); poly t, int opt 888 ASSUME: basering with characteristic 0 and local degree ordering, 889 t with isolated singularity at 0 890 RETURN: 891 @format 892 list V: Vfiltration of t on H''/H' 893 intvec V[1]: spectral numbers in increasing order 894 intvec V[2]: 895 int V[2][i]: multiplicity of spectral number V[1][i]/V[2][i] 896 if opt>=1: 616 897 list V[4]: 617 module V[ 4][i]: vector space basis of V[1][i]/V[2][i]th graded part898 module V[3][i]: vector space basis of V[1][i]/V[2][i]th graded part 618 899 in terms of V[5] 619 ideal V[5]: monomial vector space basis 900 ideal V[4]: monomial vector space basis of H''/H' 901 ideal V[5]: standard basis of Jacobian ideal 620 902 default: opt=1 621 903 @end format … … 638 920 def R=basering; 639 921 int n=nvars(R)1; 640 def G=gmsring( f,"s");922 def G=gmsring(t,"s"); 641 923 setring G; 642 924 643 int mu ,modm=ncols(gmsbasis),maxorddif(gmsbasis);925 int mu=ncols(gmsbasis); 644 926 ideal r=gmspoly*gmsbasis; 645 927 list l; 646 928 matrix U=freemodule(mu); 647 929 matrix A[mu][mu],C; 648 module H, dH=freemodule(mu),freemodule(mu);930 module H,H1=freemodule(mu),freemodule(mu); 649 931 module H0; 650 932 int k=1; … … 662 944 dbprint(printlevelvoice+2,"// compute saturation of H''"); 663 945 H0=H; 664 dH=jet(module(A*dH+s^2*diff(matrix(dH),s)),k);665 H=H*s+ dH;946 H1=jet(module(A*H1+s^2*diff(matrix(H1),s)),k); 947 H=H*s+H1; 666 948 } 667 949 A=Ak*s; … … 669 951 dbprint(printlevelvoice+2,"// compute basis of saturation of H''"); 670 952 H=std(H0); 671 int modH=maxorddif(H)/deg(s);672 dbprint(printlevelvoice+2,"// k="+string( N+modH));953 int d0=maxorddif(H)/deg(s); 954 dbprint(printlevelvoice+2,"// k="+string(d0+N)); 673 955 dbprint(printlevelvoice+2,"// compute matrix A of t"); 674 l=gmscoeffs(r, N+modH,N+modH);956 l=gmscoeffs(r,d0+N,d0+N); 675 957 C,r=l[1..2]; 676 958 A=A+C; … … 775 1057 V[ncols(eM)+1]=interred(V1); 776 1058 intvec d; 777 if(opt ==0)1059 if(opt<=0) 778 1060 { 779 1061 for(i=ncols(eM);number(eM[i])1>number(n1)/2;i) … … 876 1158 list v=imap(G,v); 877 1159 ideal m=imap(G,gmsbasis); 878 return(list(a,d,v,m)); 1160 ideal g=imap(G,gmsstd); 1161 attrib(g,"isSB",1); 1162 return(list(a,d,v,m,g)); 879 1163 } 880 1164 } … … 882 1166 { "EXAMPLE:"; echo=2; 883 1167 ring R=0,(x,y),ds; 884 poly f=x5+x2y2+y5; 885 vfilt(f); 886 } 887 /////////////////////////////////////////////////////////////////////////////// 888 889 proc endfilt(poly f,list V) 890 "USAGE: endfilt(f,V); poly f, list V 891 ASSUME: basering with characteristic 0 and local degree ordering, 892 f with isolated singularity at 0 1168 poly t=x5+x2y2+y5; 1169 vfilt(t); 1170 } 1171 /////////////////////////////////////////////////////////////////////////////// 1172 1173 proc endfilt(list V) 1174 "USAGE: endfilt(V); list V 1175 ASSUME: V computed by vfilt 893 1176 RETURN: 894 1177 @format 895 list V1: endomorphim filtration of V on the Jacobian algebra of f1178 list V1: endomorphim filtration of V on the Jacobian algebra 896 1179 ideal V1[1]: spectral numbers in increasing order 897 1180 intvec V1[2]: … … 908 1191 " 909 1192 { 910 if(charstr(basering)!="0") 911 { 912 ERROR("characteristic 0 expected"); 913 } 914 int n=nvars(basering)1; 915 for(int i=n+1;i>=1;i) 916 { 917 if(var(i)>1) 918 { 919 ERROR("local ordering expected"); 920 } 921 } 922 ideal sJ=std(jacob(f)); 923 if(vdim(sJ)<=0) 924 { 925 if(vdim(sJ)==0) 926 { 927 ERROR("singularity at 0 expected"); 928 } 929 else 930 { 931 ERROR("isolated singularity at 0 expected"); 932 } 933 } 934 935 def a,d,v,m=V[1..4]; 1193 def a,d,v,m,g=V[1..5]; 936 1194 int mu=ncols(m); 937 1195 938 1196 module V0=v[1]; 939 for(i =2;i<=size(v);i++)1197 for(int i=2;i<=size(v);i++) 940 1198 { 941 1199 V0=V0,v[i]; … … 947 1205 for(i=ncols(m);i>=1;i) 948 1206 { 949 M[i]=lift(V0,coeffs(reduce(m[i]*m,U, sJ),m)*V0);1207 M[i]=lift(V0,coeffs(reduce(m[i]*m,U,g),m)*V0); 950 1208 } 951 1209 … … 1062 1320 { "EXAMPLE:"; echo=2; 1063 1321 ring R=0,(x,y),ds; 1064 poly f=x5+x2y2+y5; 1065 endfilt(f,vfilt(f)); 1066 } 1067 /////////////////////////////////////////////////////////////////////////////// 1068 1069 proc spectrum(poly f) 1070 "USAGE: spectrum(f); poly f 1071 ASSUME: basering with characteristic 0 and local degree ordering, 1072 f with isolated singularity at 0 1073 RETURN: 1074 @format 1075 list Sp: spectrum of f 1076 ideal Sp[1]: spectral numbers in increasing order 1077 intvec Sp[2]: 1078 int Sp[2][i]: multiplicity of spectral number Sp[1][i] 1079 @end format 1080 SEE ALSO: spectrum_lib 1081 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; spectrum 1082 EXAMPLE: example spnumbers; shows examples 1083 " 1084 { 1085 return(sppairs(f,0)); 1086 } 1087 example 1088 { "EXAMPLE:"; echo=2; 1089 ring R=0,(x,y),ds; 1090 poly f=x5+x2y2+y5; 1091 spprint(spectrum(f)); 1092 } 1093 /////////////////////////////////////////////////////////////////////////////// 1094 1095 proc sppairs(poly f,list #) 1096 "USAGE: sppairs(f[,opt]); poly f, int opt 1097 ASSUME: basering with characteristic 0 and local degree ordering, 1098 f with isolated singularity at 0 1099 RETURN: 1100 @format 1101 if opt=0 1102 list Sp: spectrum of f 1103 ideal Sp[1]: spectral numbers in increasing order 1104 intvec Sp[2]: corresponding multiplicities of spectral numbers 1105 if opt=1: 1106 list Spp: spectral pairs of f 1107 ideal Spp[1]: spectral numbers in increasing order 1108 intvec Spp[2]: corresponding weight filtration indices in increasing order 1109 intvec Spp[3]: corresponding multiplicities of spectral pairs 1110 default: opt=1 1111 @end format 1112 SEE ALSO: spectrum_lib 1113 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 1114 spectrum; spectral pairs 1115 EXAMPLE: example sppairs; shows examples 1116 " 1117 { 1118 int opt=1; 1119 if(size(#)>0) 1120 { 1121 if(typeof(#[1])=="int") 1122 { 1123 opt=#[1]; 1124 } 1125 } 1126 1127 def R=basering; 1128 int n=nvars(R)1; 1129 def G=gmsring(f,"s"); 1130 setring(G); 1131 1132 int mu,modm=ncols(gmsbasis),maxorddif(gmsbasis); 1133 ideal r=gmspoly*gmsbasis; 1134 list l; 1135 matrix U=freemodule(mu); 1136 matrix A0[mu][mu],C; 1137 module H0; 1138 module H,dH=freemodule(mu),freemodule(mu); 1139 int k=1; 1140 while(size(reduce(H,std(H0*s)))>0) 1141 { 1142 k++; 1143 dbprint(printlevelvoice+2,"// k="+string(k)); 1144 dbprint(printlevelvoice+2,"// compute matrix A of t"); 1145 l=gmscoeffs(r,k,mu+n); 1146 C,r=l[1..2]; 1147 A0=A0+C; 1148 1149 dbprint(printlevelvoice+2,"// compute saturation of H''"); 1150 H0=H; 1151 dH=jet(module(A0*dH+s^2*diff(matrix(dH),s)),k); 1152 H=H*s+dH; 1153 } 1154 A0=A0k*s; 1155 1156 dbprint(printlevelvoice+2,"// compute basis of saturation of H''"); 1157 H=std(H0); 1158 int modH=maxorddif(H)/deg(s); 1159 dbprint(printlevelvoice+2,"// transform H'' to saturation of H''"); 1160 l=division(H,freemodule(mu)*s^k); 1161 H0=l[1]; 1162 1163 dbprint(printlevelvoice+2,"// k="+string(modH+1)); 1164 dbprint(printlevelvoice+2,"// compute matrix A of t"); 1165 l=gmscoeffs(r,modH+1,modH+1+n); 1166 C,r=l[1..2]; 1167 A0=A0+C; 1168 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 1169 l=division(H*s,A0*H+s^2*diff(matrix(H),s)); 1170 matrix A=jet(l[1],l[2],0); 1171 1172 int i,j; 1173 matrix V; 1174 if(!opt) 1175 { 1176 dbprint(printlevelvoice+2,"// transform to Jordan basis"); 1177 U=jordanbasis(A)[1]; 1178 V=lift(U,freemodule(mu)); 1179 A=V*A*U; 1180 dbprint(printlevelvoice+2,"// compute normal form of H''"); 1181 H0=std(V*H0); 1182 1183 dbprint(printlevelvoice+2,"// compute spectral numbers"); 1184 ideal a; 1185 for(i=1;i<=mu;i++) 1186 { 1187 j=leadexp(H0[i])[nvars(basering)+1]; 1188 a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)1; 1189 } 1190 1191 setring(R); 1192 return(spgen(imap(G,a))); 1193 } 1194 1195 dbprint(printlevelvoice+2, 1196 "// compute eigenvalues e with multiplicities m of A"); 1197 l=eigenval(A); 1198 def e,m=l[1..2]; 1199 dbprint(printlevelvoice+2,"// e="+string(e)); 1200 dbprint(printlevelvoice+2,"// m="+string(m)); 1201 dbprint(printlevelvoice+2,"// compute maximal integer difference of e"); 1202 int mide=maxintdif(e); 1203 dbprint(printlevelvoice+2,"// mide="+string(mide)); 1204 1205 if(mide>0) 1206 { 1207 dbprint(printlevelvoice+2,"// k="+string(modH+1+mide)); 1208 dbprint(printlevelvoice+2,"// compute matrix A of t"); 1209 l=gmscoeffs(r,modH+1+mide,modH+1+mide); 1210 C,r=l[1..2]; 1211 A0=A0+C; 1212 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 1213 l=division(H*s,A0*H+s^2*diff(matrix(H),s)); 1214 A=jet(l[1],l[2],mide); 1215 1216 intvec ide; 1217 ide[mu]=0; 1218 for(i=ncols(e);i>=1;i) 1219 { 1220 for(j=i1;j>=1;j) 1221 { 1222 k=int(e[j]e[i]); 1223 if(k>ide[j]) 1224 { 1225 ide[j]=k; 1226 } 1227 if(k>ide[i]) 1228 { 1229 ide[i]=k; 1230 } 1231 } 1232 } 1233 for(j,k=ncols(e),mu;j>=1;j) 1234 { 1235 for(i=m[j];i>=1;i) 1236 { 1237 ide[k]=ide[j]; 1238 k; 1239 } 1240 } 1241 1242 while(mide>0) 1243 { 1244 dbprint(printlevelvoice+2,"// transform to separate eigenvalues"); 1245 A0=jet(A,0); 1246 U=0; 1247 for(i=ncols(e);i>=1;i) 1248 { 1249 U=syz(power(A0e[i],n+1))+U; 1250 } 1251 V=lift(U,freemodule(mu)); 1252 A=V*A*U; 1253 H0=V*H0; 1254 1255 dbprint(printlevelvoice+2,"// transform to reduce mide by 1"); 1256 for(i=mu;i>=1;i) 1257 { 1258 for(j=mu;j>=1;j) 1259 { 1260 if(ide[i]==0&&ide[j]>0) 1261 { 1262 A[i,j]=A[i,j]/s; 1263 } 1264 else 1265 { 1266 if(ide[i]>0&&ide[j]==0) 1267 { 1268 A[i,j]=A[i,j]*s; 1269 } 1270 } 1271 } 1272 } 1273 H0=transpose(H0); 1274 for(i=mu;i>=1;i) 1275 { 1276 if(ide[i]>0) 1277 { 1278 A[i,i]=A[i,i]1; 1279 ide[i]=ide[i]1; 1280 H0[i]=H0[i]*s; 1281 } 1282 } 1283 H0=transpose(H0); 1284 1285 mide; 1286 dbprint(printlevelvoice+2,"// mide="+string(mide)); 1287 } 1288 1289 A=jet(A,0); 1290 } 1291 1292 dbprint(printlevelvoice+2,"// transform to Jordan basis"); 1293 intvec w0; 1294 l=jordanbasis(A); 1295 U,w0=l[1..2]; 1296 V=lift(U,freemodule(mu)); 1297 A0=V*A*U; 1298 1299 dbprint(printlevelvoice+2,"// compute weight filtration basis"); 1300 vector u; 1301 i=1; 1302 while(i<ncols(A)) 1303 { 1304 j=i+1; 1305 while(j<ncols(A)&&A[i,i]==A[j,j]) 1306 { 1307 if(w0[i]<w0[j]) 1308 { 1309 k=w0[i]; 1310 w0[i]=w0[j]; 1311 w0[i]=k; 1312 u=U[i]; 1313 U[i]=U[j]; 1314 U[j]=u; 1315 } 1316 j++; 1317 } 1318 if(j==ncols(A)&&A[i,i]==A[j,j]&&w0[i]<w0[j]) 1319 { 1320 k=w0[i]; 1321 w0[i]=w0[j]; 1322 w0[i]=k; 1323 u=U[i]; 1324 U[i]=U[j]; 1325 U[j]=u; 1326 } 1327 i=j; 1328 } 1329 1330 dbprint(printlevelvoice+2,"// transform to weight filtration basis"); 1331 V=lift(U,freemodule(mu)); 1332 A=V*A*U; 1333 dbprint(printlevelvoice+2,"// compute normal form of H''"); 1334 H0=std(V*H0); 1335 1336 dbprint(printlevelvoice+2,"// compute spectral pairs"); 1337 ideal a; 1338 intvec w; 1339 for(i=1;i<=mu;i++) 1340 { 1341 j=leadexp(H0[i])[nvars(basering)+1]; 1342 a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)1; 1343 w[i]=w0[j]+n; 1344 } 1345 setring(R); 1346 return(sppgen(imap(G,a),w)); 1347 } 1348 example 1349 { "EXAMPLE:"; echo=2; 1350 ring R=0,(x,y),ds; 1351 poly f=x5+x2y2+y5; 1352 spprint(sppairs(f)); 1322 poly t=x5+x2y2+y5; 1323 endfilt(vfilt(t)); 1353 1324 } 1354 1325 /////////////////////////////////////////////////////////////////////////////// … … 1360 1331 list Sp: numbers in a with multiplicities 1361 1332 ideal Sp[1]: numbers in a in increasing order 1362 intvec Sp[2]: corresponding multiplicities 1333 intvec Sp[2]: 1334 int Sp[2][i]: multiplicity of number Sp[1][i] in a 1363 1335 @end format 1364 1336 EXAMPLE: example spgen; shows examples … … 1407 1379 1408 1380 proc sppgen(ideal a,intvec w) 1409 "USAGE: sppgen(a ); ideal a1381 "USAGE: sppgen(a,w); ideal a, intvec w 1410 1382 RETURN: 1411 1383 @format 1412 1384 list Spp: pairs in a and w with multiplicities 1413 ideal Spp[1]: numbers in a 1414 intvec Spp[2]: corresponding integers in w 1415 intvec Spp[3]: corresponding multiplicities 1385 ideal Spp[1]: numbers in a in increasing order 1386 intvec Spp[2]: integers in w in decreasing order 1387 intvec Spp[3]: 1388 int Spp[3][i]: multiplicity of pair (Spp[1][i],Spp[2][i]) in a,w 1416 1389 @end format 1417 1390 EXAMPLE: example sppgen; shows examples
Note: See TracChangeset
for help on using the changeset viewer.