Changeset 61549b in git
 Timestamp:
 Nov 5, 2001, 10:17:05 AM (22 years ago)
 Branches:
 (u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
 Children:
 1f8111cd2acad6f5a631540b24272cf5aca44000
 Parents:
 5f403d5fffa49541688757a3ff48dc8d052361f2
 Location:
 Singular/LIB
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/gaussman.lib
r5f403d5 r61549b 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: gaussman.lib,v 1.5 4 20010831 18:18:08 mschulze Exp $";2 version="$Id: gaussman.lib,v 1.55 20011105 09:16:48 mschulze Exp $"; 3 3 category="Singularities"; 4 4 … … 17 17 monodromy(t); Jordan data of monodromy of t 18 18 spectrum(t); spectrum of t 19 sppnormalize(a,w[,m]); normalize spectral pairs20 19 sppairs(t); spectral pairs of t 20 sppnf(a,w[,m][,vV]); normalize spectral pairs 21 vfilt(t); Vfiltration of t on H''/H' 22 vwfilt(t); weight refined Vfiltration of t on H''/H' 21 23 saito(t); matrix A0+A1*s of t on H'' 22 vfilt(t[,opt]); Vfiltration on H''/H' or spectrum of t 23 endfilt(t,V); endomorphism filtration of Vfiltration V 24 endvfilt(V); endomorphism Vfiltration 24 25 spprint(Sp); print spectrum Sp 25 26 sppprint(Spp); print spectral pairs Spp … … 413 414 dbprint(printlevelvoice+2,"// compute basis of saturation of H''"); 414 415 H=std(H0); 415 int d0=maxdeg1(H);416 416 417 417 dbprint(printlevelvoice+2,"// transform H'' to saturation of H''"); … … 419 419 H0=l[1]; 420 420 421 return(A0,r,H,H0 );422 } 423 /////////////////////////////////////////////////////////////////////////////// 424 425 static proc tmatrix(matrix A0,ideal r,module H,int K,int K0)421 return(A0,r,H,H0,k); 422 } 423 /////////////////////////////////////////////////////////////////////////////// 424 425 static proc tmatrix(matrix A0,ideal r,module H,int k0,int K,int K0) 426 426 { 427 427 dbprint(printlevelvoice+2,"// compute matrix A of t"); 428 int d0=maxdeg1(H); 429 dbprint(printlevelvoice+2,"// k="+string(K+d0+1)); 430 list l=gmscoeffs(r,K+d0+1,K0+d0+1); 428 dbprint(printlevelvoice+2,"// k="+string(K+k0+1)); 429 list l=gmscoeffs(r,K+k0+1,K0+k0+1); 431 430 matrix C; 432 431 C,r=l[1..2]; … … 441 440 /////////////////////////////////////////////////////////////////////////////// 442 441 443 static proc eigenvals(matrix A0,ideal r,module H,int K0)442 static proc eigenvals(matrix A0,ideal r,module H,int k0,int K0) 444 443 { 445 444 dbprint(printlevelvoice+2, 446 445 "// compute eigenvalues e with multiplicities m of A"); 447 446 matrix A; 448 A,A0,r=tmatrix(A0,r,H, 0,K0);447 A,A0,r=tmatrix(A0,r,H,k0,0,K0); 449 448 list l=eigenvalues(A); 450 449 def e,m=l[1..2]; … … 452 451 dbprint(printlevelvoice+2,"// m="+string(m)); 453 452 454 return(e,m,A0,r); 455 } 456 /////////////////////////////////////////////////////////////////////////////// 457 458 static proc transform(ideal e,intvec m,matrix A,matrix A0,ideal r,module H,module H0,int K,int K0) 459 { 453 return(e,m,A0,r,int(max(e)min(e))); 454 } 455 /////////////////////////////////////////////////////////////////////////////// 456 457 static proc transform(matrix A,matrix A0,ideal r,module H,module H0,ideal e,intvec m,int k0,int k1,int K,int K0) 458 { 459 int mu=ncols(gmsbasis); 460 460 461 dbprint(printlevelvoice+2,"// compute minimum e0 and maximum e1 of e"); 461 462 number e0,e1=min(e),max(e); 462 463 dbprint(printlevelvoice+2,"// e0="+string(e0)); 463 464 dbprint(printlevelvoice+2,"// e1="+string(e1)); 464 int d1=int(e1e0);465 A,A0,r=tmatrix(A0,r,H,K+d1,K0+d1);465 A,A0,r=tmatrix(A0,r,H,k0,K+k1,K0+k1); 466 module U0=s^k0*freemodule(mu); 466 467 467 468 if(e1>=e0+1) … … 482 483 A=V*A*U; 483 484 H0=V*H0; 485 U0=U0*U; 484 486 485 487 dbprint(printlevelvoice+2,"// transform to reduce e1 by 1"); … … 514 516 A[i,i]=A[i,i]1; 515 517 H0[i]=H0[i]*s; 518 U0[i]=U0[i]/s; 516 519 } 517 520 e[i0]=e[i0]1; … … 524 527 H0=transpose(H0); 525 528 526 l=spn ormalize(e,m);529 l=spnf(e,m); 527 530 e,m=l[1..2]; 528 531 … … 534 537 } 535 538 536 return( e,m,A,A0,r,H0);539 return(A,A0,r,H0,U0,e,m); 537 540 } 538 541 /////////////////////////////////////////////////////////////////////////////// … … 553 556 setring(G); 554 557 555 int mu=ncols(gmsbasis);556 558 matrix A; 559 module U0; 557 560 ideal e; 558 561 intvec m; 559 560 def A0,r,H,H0=saturate(n); 561 e,m,A0,r=eigenvals(A0,r,H,n); 562 e,m,A,A0,r,H0=transform(e,m,A,A0,r,H,H0,0,0); 562 int k1; 563 564 def A0,r,H,H0,k0=saturate(n); 565 e,m,A0,r,k1=eigenvals(A0,r,H,k0,n); 566 A,A0,r,H0,U0,e,m=transform(A,A0,r,H,H0,e,m,k0,k1,0,0); 563 567 564 568 setring(R); 565 matrix A=imap(G,A); 566 ideal e=imap(G,e); 567 568 return(jordan(A)); 569 return(jordan(imap(G,A),imap(G,e),m)); 569 570 } 570 571 example … … 583 584 @format 584 585 list Sp: spectrum of t 585 ideal Sp[1]: spectral numbers in increasing order 586 intvec Sp[2]: 587 int Sp[2][i]: multiplicity of spectral number Sp[1][i] 586 ideal Sp[1]: Vfiltration indices in increasing order 587 intvec Sp[2]: weight filtration indices in decreasing order 588 intvec Sp[3]: 589 int Sp[3][i]: multiplicity of spectral number Sp[1][i] 588 590 @end format 589 591 SEE ALSO: spectrum_lib 590 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; spectrum 591 EXAMPLE: example spnumbers; shows examples 592 " 593 { 594 list l=sppairs(t); 595 return(spnormalize(l[1],l[3])); 592 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 593 mixed Hodge structure; spectrum 594 EXAMPLE: example spectrum; shows examples 595 " 596 { 597 list l=vwfilt(t); 598 return(spnf(l[1],l[3])); 596 599 } 597 600 example … … 600 603 poly t=x5+x2y2+y5; 601 604 spprint(spectrum(t)); 602 }603 ///////////////////////////////////////////////////////////////////////////////604 605 static proc sppappend(list l,number a,int w,int m)606 {607 if(size(l)==0)608 {609 l=list(ideal(a),intvec(w),intvec(m));610 }611 else612 {613 int n=ncols(l[1]);614 l[1][n+1]=a;615 l[2][n+1]=w;616 l[3][n+1]=m;617 }618 return(l);619 }620 ///////////////////////////////////////////////////////////////////////////////621 622 proc sppnormalize(ideal a,intvec w,list #)623 "USAGE: sppnormalize(a,w[,m]);624 RETURN:625 @format626 list Spp: normalized spectral pairs (a,w,m)627 ideal Spp[1]: numbers in a in increasing order628 intvec Spp[2]: integers in w in decreasing order629 intvec Spp[3]:630 int Spp[3][i]: multiplicity of pair (Spp[1][i],Spp[2][i]) in a,w631 @end format632 EXAMPLE: example sppnormalize; shows examples633 "634 {635 intvec m;636 int i,j;637 if(size(#)==0)638 {639 for(i=ncols(a);i>=1;i)640 {641 m[i]=1;642 }643 }644 else645 {646 m=#[1];647 }648 649 list l;650 number a0;651 int w0,m0;652 for(i=ncols(a);i>=1;i)653 {654 if(m[i]!=0)655 {656 for(j=i1;j>=1;j)657 {658 if(m[j]!=0)659 {660 if(number(a[i])>number(a[j])661 (number(a[i])==number(a[j])&&w[i]<w[j]))662 {663 a0=number(a[i]);664 a[i]=a[j];665 a[j]=a0;666 w0=w[i];667 w[i]=w[j];668 w[j]=w0;669 m0=m[i];670 m[i]=m[j];671 m[j]=m0;672 }673 if(number(a[i])==number(a[j])&&w[i]==w[j])674 {675 m[i]=m[i]+m[j];676 m[j]=0;677 }678 }679 }680 l=sppappend(l,number(a[i]),w[i],m[i]);681 }682 }683 684 return(l);685 }686 example687 { "EXAMPLE:"; echo=2;688 605 } 689 606 /////////////////////////////////////////////////////////////////////////////// … … 695 612 RETURN: 696 613 @format 697 list Spp: 698 ideal Spp[1]: spectral numbers in increasing order 699 intvec Spp[2]: weight filtration indices in decreasing order 614 list Spp: spectrum of t 615 ideal Spp[1],intvec Spp[2]: spectral pairs in in/decreasing lex. order 700 616 intvec Spp[3]: 701 617 int Spp[3][i]: multiplicity of spectral pair (Spp[1][i],Spp[2][i]) … … 703 619 SEE ALSO: spectrum_lib 704 620 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 621 mixed Hodge structure; spectrum; spectral pairs 622 EXAMPLE: example sppairs; shows examples 623 " 624 { 625 list l=vwfilt(t); 626 return(list(l[1],l[2],l[3])); 627 } 628 example 629 { "EXAMPLE:"; echo=2; 630 ring R=0,(x,y),ds; 631 poly t=x5+x2y2+y5; 632 sppprint(sppairs(t)); 633 } 634 /////////////////////////////////////////////////////////////////////////////// 635 636 proc sppnf(ideal a,intvec w,list #) 637 "USAGE: sppnf(a,w[,m][,vV]); ideal a, intvec w, intvec m, module v, list V 638 RETURN: 639 @format 640 list Spp: normalized spectral pairs (a,w,m[,V]) 641 ideal Spp[1]: numbers in a in increasing order 642 intvec Spp[2]: integers in w in decreasing order 643 intvec Spp[3]: 644 int Spp[3][i]: multiplicity of pair (Spp[1][i],Spp[2][i]) in (a,w) 645 list Spp[4]: 646 module Spp[4][i]: module associated to pair (Spp[1][i],Spp[2][i]) 647 @end format 648 EXAMPLE: example sppnorm; shows examples 649 " 650 { 651 int n=ncols(a); 652 intvec m; 653 module v; 654 list V; 655 int i,j; 656 while(i<size(#)) 657 { 658 i++; 659 if(typeof(#[i])=="intvec") 660 { 661 m=#[i]; 662 } 663 if(typeof(#[i])=="module") 664 { 665 v=#[i]; 666 for(j=n;j>=1;j) 667 { 668 V[j]=module(v[j]); 669 } 670 } 671 if(typeof(#[i])=="list") 672 { 673 V=#[i]; 674 } 675 } 676 if(m==0) 677 { 678 for(i=n;i>=1;i) 679 { 680 m[i]=1; 681 } 682 } 683 684 int k; 685 ideal a0; 686 intvec w0,m0; 687 list V0; 688 number a1; 689 int w1,m1; 690 for(i=n;i>=1;i) 691 { 692 if(m[i]!=0) 693 { 694 for(j=i1;j>=1;j) 695 { 696 if(m[j]!=0) 697 { 698 if(number(a[i])>number(a[j]) 699 (number(a[i])==number(a[j])&&w[i]<w[j])) 700 { 701 a1=number(a[i]); 702 a[i]=a[j]; 703 a[j]=a1; 704 w1=w[i]; 705 w[i]=w[j]; 706 w[j]=w1; 707 m1=m[i]; 708 m[i]=m[j]; 709 m[j]=m1; 710 if(size(V)>0) 711 { 712 v=V[i]; 713 V[i]=V[j]; 714 V[j]=v; 715 } 716 } 717 if(number(a[i])==number(a[j])&&w[i]==w[j]) 718 { 719 m[i]=m[i]+m[j]; 720 m[j]=0; 721 if(size(V)>0) 722 { 723 V[i]=V[i]+V[j]; 724 } 725 } 726 } 727 } 728 k++; 729 a0[k]=a[i]; 730 w0[k]=w[i]; 731 m0[k]=m[i]; 732 if(size(V)>0) 733 { 734 V0[k]=V[i]; 735 } 736 } 737 } 738 739 if(size(V0)>0) 740 { 741 n=size(V0); 742 module U=std(V0[n]); 743 for(i=n1;i>=1;i) 744 { 745 V0[i]=simplify(reduce(V0[i],U),1); 746 if(i>=2) 747 { 748 U=std(U+V0[i]); 749 } 750 } 751 } 752 753 list l; 754 if(k>0) 755 { 756 l=a0,w0,m0; 757 if(size(V0)>0) 758 { 759 l[4]=V0; 760 } 761 } 762 return(l); 763 } 764 example 765 { "EXAMPLE:"; echo=2; 766 } 767 /////////////////////////////////////////////////////////////////////////////// 768 769 proc vfilt(poly t) 770 "USAGE: vfilt(t); poly t 771 ASSUME: basering with characteristic 0 and local degree ordering, 772 t with isolated citical point 0 773 RETURN: 774 @format 775 list V: Vfiltration on H''/H' 776 ideal V[1]: spectral numbers in increasing order 777 intvec V[2]: 778 int V[2][i]: multiplicity of spectral number V[1][i] 779 list V[4]: 780 module V[4][i]: vector space basis of V[1][i]th graded part 781 in terms of V[4] 782 ideal V[4]: monomial vector space basis of H''/H' 783 ideal V[5]: standard basis of Jacobian ideal 784 @end format 785 SEE ALSO: spectrum_lib 786 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 787 mixed Hodge structure; Vfiltration; spectrum 788 EXAMPLE: example vfilt; shows examples 789 " 790 { 791 list l=vwfilt(t); 792 return(spnf(l[1],l[3],l[4])+list(l[5],l[6])); 793 } 794 example 795 { "EXAMPLE:"; echo=2; 796 ring R=0,(x,y),ds; 797 poly t=x5+x2y2+y5; 798 vfilt(t); 799 } 800 /////////////////////////////////////////////////////////////////////////////// 801 802 proc vwfilt(poly t) 803 "USAGE: vwfilt(t); poly t 804 ASSUME: basering with characteristic 0 and local degree ordering, 805 t with isolated citical point 0 806 RETURN: 807 @format 808 list VW: weight refined Vfiltration on H''/H' 809 ideal VW[1]: spectral numbers in increasing order 810 intvec VW[2]: weights in decreasing order 811 intvec VW[3]: 812 int VW[3][i]: multiplicity of spectral pair (VW[1][i],VW[2][i]) 813 list VW[4]: 814 module VW[4][i]: vector space basis of (VW[1][i],VW[2][i])th graded part 815 in terms of VW[5] 816 ideal VW[5]: monomial vector space basis of H''/H' 817 ideal VW[6]: standard basis of Jacobian ideal 818 @end format 819 SEE ALSO: spectrum_lib 820 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 821 mixed Hodge structure; Vfiltration; weight filtration; 705 822 spectrum; spectral pairs 706 EXAMPLE: example sppairs; shows examples823 EXAMPLE: example vwfilt; shows examples 707 824 " 708 825 { … … 714 831 int mu=ncols(gmsbasis); 715 832 matrix A; 833 module U0; 716 834 ideal e; 717 835 intvec m; 718 719 def A0,r,H,H0=saturate(n); 720 e,m,A0,r=eigenvals(A0,r,H,n); 721 e,m,A,A0,r,H0=transform(e,m,A,A0,r,H,H0,0,0); 836 int k1; 837 838 def A0,r,H,H0,k0=saturate(n); 839 e,m,A0,r,k1=eigenvals(A0,r,H,k0,n); 840 A,A0,r,H0,U0,e,m=transform(A,A0,r,H,H0,e,m,k0,k1,0,0); 722 841 723 842 dbprint(printlevelvoice+2,"// compute weight filtration basis"); 724 843 list l=jordanbasis(A,e,m); 725 844 def U,v=l[1..2]; 845 kill l; 726 846 vector u0; 727 847 int v0; 728 int i,j,k, k0;729 for(k, k0=1,1;k0<=ncols(e);k,k0=k+m[k0],k0+1)730 { 731 for(i=k+m[ k0]1;i>=k+1;i)848 int i,j,k,l; 849 for(k,l=1,1;l<=ncols(e);k,l=k+m[l],l+1) 850 { 851 for(i=k+m[l]1;i>=k+1;i) 732 852 { 733 853 for(j=i1;j>=k;j) … … 747 867 dbprint(printlevelvoice+2,"// compute normal form of H''"); 748 868 H0=std(V*H0); 869 U0=U0*U; 749 870 750 871 dbprint(printlevelvoice+2,"// compute spectral pairs"); … … 754 875 { 755 876 j=leadexp(H0[i])[nvars(basering)+1]; 756 a[i]=A[j,j]+ deg(lead(H0[i]))/deg(s)1;877 a[i]=A[j,j]+ord(H0[i])/deg(s)1; 757 878 w[i]=v[j]+n; 758 879 } 880 kill v; 881 module v=simplify(grmat(H*U0*H0,2*k0),1); 759 882 760 883 setring(R); 761 762 return(sppnormalize(imap(G,a),w)); 884 ideal g=imap(G,gmsstd); 885 attrib(g,"isSB",1); 886 return(sppnf(imap(G,a),w,imap(G,v))+list(imap(G,gmsbasis),g)); 763 887 } 764 888 example … … 766 890 ring R=0,(x,y),ds; 767 891 poly t=x5+x2y2+y5; 768 sppprint(sppairs(t)); 769 } 770 /////////////////////////////////////////////////////////////////////////////// 771 772 static proc grmat(matrix A,int k) 773 { 774 int i,j; 775 for(i=1;i<=ncols(A);i++) 776 { 777 for(j=1;j<=nrows(A);j++) 778 { 779 A[i,j]=jet(A[i,j]/var(1)^k,0); 780 } 781 } 782 return(A); 783 } 784 /////////////////////////////////////////////////////////////////////////////// 785 786 proc saito(poly t,list #) 787 "USAGE: saito(t); poly t 788 ASSUME: basering with characteristic 0 and local degree ordering, 789 t with isolated citical point 0 790 RETURN: list A: matrix A[1]+A[2]*s of t on H'' 791 SEE ALSO: spectrum_lib 792 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 793 monodromy; spectrum; spectral pairs 794 EXAMPLE: example saito; shows examples 795 " 796 { 797 def R=basering; 798 int n=nvars(R)1; 799 def G=gmsring(t,"s"); 800 setring(G); 801 802 int mu=ncols(gmsbasis); 803 matrix A; 804 ideal e; 805 intvec m; 806 807 def A0,r,H,H0=saturate(2*n+mu+1); 808 e,m,A0,r=eigenvals(A0,r,H,n+1); 809 int d0,d1=maxdeg1(H),int(max(e)min(e)); 810 e,m,A,A0,r,H0=transform(e,m,A,A0,r,H,H0,d0+d1+1,d0+d1+1); 811 d0=maxdeg1(H0); 812 813 dbprint(printlevelvoice+2,"// transform to Jordan basis"); 814 module U=jordanbasis(A,e,m)[1]; 815 module V=inverse(U); 816 A=V*A*U; 817 H=V*H0; 818 819 dbprint(printlevelvoice+2,"// compute splitting of V filtration"); 820 int i,j,k; 821 U=freemodule(mu); 822 matrix U0[mu][mu]; 823 matrix u0[mu^2][1]; 824 A0=commutator(jet(A,0)); 825 for(k=1;k<=d0+1;k++) 826 { 827 V=0; 828 for(j=0;j<k;j++) 829 { 830 U0=matrix(U0)grmat(A,kj)*grmat(U,j); 831 } 832 u0=U0[1..mu,1..mu]; 833 u0=inverse(A0+k)*u0; 834 U0=u0[1..mu^2,1]; 835 U=matrix(U)+s^k*U0; 836 } 837 838 dbprint(printlevelvoice+2,"// transform to splitting basis"); 839 A=jet(A,0); 840 H=std(jet(division(U,H)[1],d0)); 841 842 dbprint(printlevelvoice+2,"// compute N"); 843 matrix N=A; 844 for(i=1;i<=ncols(N);i++) 845 { 846 N[i,i]=0; 847 } 848 849 dbprint(printlevelvoice+2,"// compute N splitting of Hodge filtration"); 850 module F; 851 matrix K[mu][1]; 852 U=0; 853 i=1; 854 for(k=1;k<=ncols(e);k++) 855 { 856 j=0; 857 K=matrix(0,mu,0); 858 while(j<m[k]) 859 { 860 K=K+gen(i); 861 i++; 862 j++; 863 } 864 j=0; 865 F=grmat(H,0); 866 V=intersect(F,K); 867 while(size(V)==0) 868 { 869 j++; 870 V=V+grmat(H,j); 871 V=intersect(F,K) 872 } 873 while(size(V)>0) 874 { 875 U=U+V; 876 V=N*V; 877 } 878 } 879 880 dbprint(printlevelvoice+2,"// transform to N splitting basis"); 881 V=inverse(U); 882 A=V*A*U; 883 H=V*H*U; 884 885 dbprint(printlevelvoice+2,"// compute matrix A0+sA1 of t"); 886 // degBound=d0; 887 // option("redSB"); 888 // H=std(H); 889 H=simplify(lead(std(H)),1); 890 // degBound=0; 891 list l=division(H,s*A*H+s^2*diff(matrix(H),s)); 892 A=jet(l[1],l[2],d0+1); 893 A0=grmat(A,0); 894 A=grmat(A,1); 895 896 setring(R); 897 898 return(list(imap(G,A0),imap(G,A))); 899 } 900 example 901 { "EXAMPLE:"; echo=2; 902 ring R=0,(x,y),ds; 903 poly t=x5+x2y2+y5; 904 list A=saito(t); 905 print(A[1]); 906 print(A[2]); 907 } 908 /////////////////////////////////////////////////////////////////////////////// 909 910 proc vfilt(poly t,list #) 911 "USAGE: vfilt(t[,opt]); poly t, int opt 892 vwfilt(t); 893 } 894 /////////////////////////////////////////////////////////////////////////////// 895 896 proc vfilt1(poly t,list #) 897 "USAGE: vfilt1(t[,opt]); poly t, int opt 912 898 ASSUME: basering with characteristic 0 and local degree ordering, 913 899 t with isolated citical point 0 … … 930 916 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 931 917 Hodge filtration; Vfiltration; spectrum 932 EXAMPLE: example vfilt ; shows examples918 EXAMPLE: example vfilt1; shows examples 933 919 " 934 920 { … … 1190 1176 ring R=0,(x,y),ds; 1191 1177 poly t=x5+x2y2+y5; 1192 vfilt(t); 1193 } 1194 /////////////////////////////////////////////////////////////////////////////// 1195 1196 proc endfilt(list V) 1197 "USAGE: endfilt(V); list V 1178 vfilt1(t); 1179 } 1180 /////////////////////////////////////////////////////////////////////////////// 1181 1182 static proc grmat(matrix A,int k) 1183 { 1184 int i,j; 1185 for(i=1;i<=ncols(A);i++) 1186 { 1187 for(j=1;j<=nrows(A);j++) 1188 { 1189 A[i,j]=jet(A[i,j]/var(1)^k,0); 1190 } 1191 } 1192 return(A); 1193 } 1194 /////////////////////////////////////////////////////////////////////////////// 1195 1196 proc saito(poly t,list #) 1197 "USAGE: saito(t); poly t 1198 ASSUME: basering with characteristic 0 and local degree ordering, 1199 t with isolated citical point 0 1200 RETURN: list A: matrix A[1]+A[2]*s of t on H'' 1201 SEE ALSO: spectrum_lib 1202 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 1203 mixed Hodge structure; Vfiltration; weight filtration; 1204 spectrum; spectral pairs 1205 EXAMPLE: example saito; shows examples 1206 " 1207 { 1208 def R=basering; 1209 int n=nvars(R)1; 1210 def G=gmsring(t,"s"); 1211 setring(G); 1212 1213 int mu=ncols(gmsbasis); 1214 matrix A; 1215 module U0; 1216 ideal e; 1217 intvec m; 1218 int k1; 1219 1220 def A0,r,H,H0,k0=saturate(2*n+mu1); 1221 e,m,A0,r,k1=eigenvals(A0,r,H,k0,n); 1222 A,A0,r,H0,U0,e,m=transform(A,A0,r,H,H0,e,m,k0,k1,k0+k1,k0+k1); 1223 1224 dbprint(printlevelvoice+2,"// transform to Jordan basis"); 1225 module U=jordanbasis(A,e,m)[1]; 1226 matrix V=inverse(U); 1227 A=V*A*U; 1228 H=V*H0; 1229 1230 dbprint(printlevelvoice+2,"// compute splitting of Vfiltration"); 1231 int i,j,k; 1232 U=freemodule(mu); 1233 V=matrix(0,mu,mu); 1234 matrix v[mu^2][1]; 1235 A0=commutator(jet(A,0)); 1236 for(k=1;k<=k0+k1;k++) 1237 { 1238 for(j=0;j<k;j++) 1239 { 1240 V=matrix(V)grmat(A,kj)*grmat(U,j); 1241 } 1242 v=V[1..mu,1..mu]; 1243 v=inverse(A0+k)*v; 1244 V=v[1..mu^2,1]; 1245 U=matrix(U)+s^k*V; 1246 } 1247 1248 dbprint(printlevelvoice+2,"// transform to Vsplitting basis"); 1249 A=jet(A,0); 1250 H=jet(division(U,H)[1],k0+k1); 1251 H=std(H); 1252 1253 dbprint(printlevelvoice+2,"// compute Vleading terms of H''"); 1254 int i0,j0; 1255 module H1=H; 1256 for(k=ncols(H1);k>=1;k) 1257 { 1258 i0=leadexp(H1[k])[nvars(basering)+1]; 1259 j0=ord(H1[k])/deg(s); 1260 H0[k]=lead(H1[k]); 1261 H1[k]=H1[k]lead(H1[k]); 1262 if(H1[k]!=0) 1263 { 1264 i=leadexp(H1[k])[nvars(basering)+1]; 1265 j=ord(H1[k])/deg(s); 1266 while(A[i,i]+j==A[i0,i0]+j0) 1267 { 1268 H0[k]=H0[k]+lead(H1[k]); 1269 H1[k]=H1[k]lead(H1[k]); 1270 i=leadexp(H1[k])[nvars(basering)+1]; 1271 j=ord(H1[k])/deg(s); 1272 } 1273 } 1274 } 1275 H0=simplify(H0,1); 1276 1277 dbprint(printlevelvoice+2,"// compute N"); 1278 matrix N=A; 1279 for(i=1;i<=ncols(N);i++) 1280 { 1281 N[i,i]=0; 1282 } 1283 1284 dbprint(printlevelvoice+2,"// compute splitting of Hodge filtration"); 1285 U=0; 1286 module U1; 1287 module C; 1288 list F,I; 1289 module F0,I0; 1290 for(i0,j0=1,1;i0<=ncols(e);i0++) 1291 { 1292 C=matrix(0,mu,1); 1293 for(j=m[i0];j>=1;j,j0=j1,j0+1) 1294 { 1295 C=C+gen(j0); 1296 } 1297 F0=intersect(C,H0); 1298 F=list(); 1299 j=0; 1300 while(size(F0)>0) 1301 { 1302 j++; 1303 F[j]=matrix(0,mu,1); 1304 if(size(jet(F0,0))>0) 1305 { 1306 for(i=ncols(F0);i>=1;i) 1307 { 1308 if(ord(F0[i])==0) 1309 { 1310 F[j]=F[j]+F0[i]; 1311 } 1312 } 1313 } 1314 for(i=ncols(F0);i>=1;i) 1315 { 1316 F0[i]=F0[i]/s; 1317 } 1318 } 1319 1320 I=list(); 1321 I0=module(); 1322 U0=std(0); 1323 for(i=size(F);i>=1;i) 1324 { 1325 I[i]=module(); 1326 } 1327 for(i=1;i<=size(F);i++) 1328 { 1329 I0=reduce(F[i],U0); 1330 j=i; 1331 while(size(I0)>0) 1332 { 1333 U0=std(U0+I0); 1334 I[j]=I[j]+I0; 1335 I0=reduce(N*I0,U0); 1336 j++; 1337 } 1338 } 1339 1340 for(i=1;i<=size(I);i++) 1341 { 1342 U=U+I[i]; 1343 } 1344 } 1345 1346 dbprint(printlevelvoice+2,"// transform to Hodge splitting basis"); 1347 V=inverse(U); 1348 A=V*A*U; 1349 H=V*H; 1350 1351 dbprint(printlevelvoice+2,"// compute reduced standard basis of H''"); 1352 ring S=0,s,ds; 1353 module H=imap(G,H); 1354 degBound=k0+k1+1; 1355 option("redSB"); 1356 H=std(H); 1357 degBound=0; 1358 H=simplify(jet(H,k0+k1),1); 1359 setring(G); 1360 H=imap(S,H); 1361 1362 dbprint(printlevelvoice+2,"// compute matrix A0+sA1 of t"); 1363 list l=division(H,s*A*H+s^2*diff(matrix(H),s)); 1364 A=jet(l[1],l[2],k0+k1+1); 1365 A0=grmat(A,0); 1366 A=grmat(A,1); 1367 1368 setring(R); 1369 return(list(imap(G,A0),imap(G,A))); 1370 } 1371 example 1372 { "EXAMPLE:"; echo=2; 1373 ring R=0,(x,y),ds; 1374 poly t=x5+x2y2+y5; 1375 list A=saito(t); 1376 print(A[1]); 1377 print(A[2]); 1378 } 1379 /////////////////////////////////////////////////////////////////////////////// 1380 1381 proc endvfilt(list V) 1382 "USAGE: endvwfilt(V); list V 1198 1383 ASSUME: V computed by vfilt 1199 1384 RETURN: 1200 1385 @format 1201 list V1: endomorphim filtration of V on the Jacobian algebra 1202 ideal V1[1]: spectral numbers in increasing order 1203 intvec V1[2]: 1204 int V1[2][i]: multiplicity of spectral number V1[1][i] 1205 list V1[3]: 1206 module V1[3][i]: vector space basis of the V1[1][i]th graded part 1207 in terms of V1[4] 1208 ideal V1[4]: monomial vector space basis 1386 list EV: endomorphism Vfiltration on the Jacobian algebra 1387 ideal EV[1]: spectral numbers in increasing order 1388 intvec EV[2]: 1389 int EV[2][i]: multiplicity of spectral pair (EV[1][i],EV[2][i]) 1390 list EV[3]: 1391 module EV[3][i]: vector space basis of the (EV[1][i],EV[2][i])th 1392 graded part in terms of EV[4] 1393 ideal EV[4]: monomial vector space basis 1394 ideal EV[5]: standard basis of Jacobian ideal 1209 1395 @end format 1210 1396 SEE ALSO: spectrum_lib 1211 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; spectrum; 1212 Hodge filtration; Vfiltration 1213 EXAMPLE: example endfilt; shows examples 1397 KEYWORDS: singularities; GaussManin connection; Brieskorn lattice; 1398 mixed Hodge structure; Vfiltration; weight filtration; 1399 endomorphism filtration 1400 EXAMPLE: example endvwfilt; shows examples 1214 1401 " 1215 1402 { … … 1269 1456 } 1270 1457 1271 dbprint(printlevelvoice+2,"// collect conditions for V1["+string(b0)+"]");1458 dbprint(printlevelvoice+2,"// collect conditions for EV["+string(b0)+"]"); 1272 1459 j=ncols(a); 1273 1460 j0=mu; … … 1338 1525 i++; 1339 1526 } 1340 return(list(a1,d1,v1,m ));1527 return(list(a1,d1,v1,m,g)); 1341 1528 } 1342 1529 example … … 1344 1531 ring R=0,(x,y),ds; 1345 1532 poly t=x5+x2y2+y5; 1346 end filt(vfilt(t));1533 endvfilt(vfilt(t)); 1347 1534 } 1348 1535 /////////////////////////////////////////////////////////////////////////////// 
Singular/LIB/linalg.lib
r5f403d5 r61549b 1 1 //GMG last modified: 04/25/2000 2 2 ////////////////////////////////////////////////////////////////////////////// 3 version="$Id: linalg.lib,v 1. 19 20010825 14:56:41mschulze Exp $";3 version="$Id: linalg.lib,v 1.20 20011105 09:17:05 mschulze Exp $"; 4 4 category="Linear Algebra"; 5 5 info=" … … 26 26 pos_def(A,i); test symmetric matrix for positive definiteness 27 27 hessenberg(M); transforms M to Hessenberg form 28 spn ormalize(a[,m]);normalize spectrum29 eigenval (M);eigenvalues of M with multiplicities28 spnf(a[,m][,vV]); normalize spectrum 29 eigenvalues(M); eigenvalues of M with multiplicities 30 30 jordan(M); Jordan data of constant square matrix M 31 31 jordanbasis(M); Jordan basis of constant square matrix M … … 1335 1335 /////////////////////////////////////////////////////////////////////////////// 1336 1336 1337 static proc spappend(list l,number e,int m) 1338 { 1339 if(size(l)==0) 1340 { 1341 l=list(ideal(e),intvec(m)); 1342 } 1343 else 1344 { 1345 int n=ncols(l[1]); 1346 l[1][n+1]=e; 1347 l[2][n+1]=m; 1348 } 1349 return(l); 1350 } 1351 /////////////////////////////////////////////////////////////////////////////// 1352 1353 proc spnormalize(ideal e,list #) 1354 "USAGE: spnormalize(a,w[,m]); 1337 proc spnf(ideal e,list #) 1338 "USAGE: spnf(e[,m][,vV]); ideal e, intvec m, module v, list V 1355 1339 RETURN: 1356 1340 @format 1357 list Sp: normalized spectrum ( a,m)1341 list Sp: normalized spectrum (e,m[,V]) 1358 1342 ideal Sp[1]: numbers in a in increasing order 1359 1343 intvec Sp[2]: 1360 int Sp[2][i]: multiplicity of number Sp[1][i] in a 1344 int Sp[2][i]: multiplicity of number Sp[1][i] in e 1345 list Spp[3]: 1346 module Spp[3][i]: module associated to number Spp[1][i] 1361 1347 @end format 1362 EXAMPLE: example spnormalize; shows examples1348 EXAMPLE: example spnf; shows examples 1363 1349 " 1364 1350 { 1351 int n=ncols(e); 1365 1352 intvec m; 1353 module v; 1354 list V; 1366 1355 int i,j; 1367 if(size(#)==0) 1368 { 1369 for(i=ncols(e);i>=1;i) 1356 while(i<size(#)) 1357 { 1358 i++; 1359 if(typeof(#[i])=="intvec") 1360 { 1361 m=#[i]; 1362 } 1363 if(typeof(#[i])=="module") 1364 { 1365 v=#[i]; 1366 for(j=n;j>=1;j) 1367 { 1368 V[j]=module(v[j]); 1369 } 1370 } 1371 if(typeof(#[i])=="list") 1372 { 1373 V=#[i]; 1374 } 1375 } 1376 if(m==0) 1377 { 1378 for(i=n;i>=1;i) 1370 1379 { 1371 1380 m[i]=1; 1372 1381 } 1373 1382 } 1374 else 1375 { 1376 m=#[1]; 1377 } 1378 1379 list l; 1380 number e0; 1381 int m0; 1382 for(i=ncols(e);i>=1;i) 1383 1384 int k; 1385 ideal e0; 1386 intvec m0; 1387 list V0; 1388 number e1; 1389 int m1; 1390 for(i=n;i>=1;i) 1383 1391 { 1384 1392 if(m[i]!=0) … … 1390 1398 if(number(e[i])>number(e[j])) 1391 1399 { 1392 e 0=number(e[i]);1400 e1=number(e[i]); 1393 1401 e[i]=e[j]; 1394 e[j]=e 0;1395 m 0=m[i];1402 e[j]=e1; 1403 m1=m[i]; 1396 1404 m[i]=m[j]; 1397 m[j]=m0; 1405 m[j]=m1; 1406 if(size(V)>0) 1407 { 1408 v=V[i]; 1409 V[i]=V[j]; 1410 V[j]=v; 1411 } 1398 1412 } 1399 1413 if(number(e[i])==number(e[j])) … … 1401 1415 m[i]=m[i]+m[j]; 1402 1416 m[j]=0; 1417 if(size(V)>0) 1418 { 1419 V[i]=V[i]+V[j]; 1420 } 1403 1421 } 1404 1422 } 1405 1423 } 1406 l=spappend(l,number(e[i]),m[i]); 1407 } 1408 } 1409 1424 k++; 1425 e0[k]=e[i]; 1426 m0[k]=m[i]; 1427 if(size(V)>0) 1428 { 1429 V0[k]=V[i]; 1430 } 1431 } 1432 } 1433 1434 if(size(V0)>0) 1435 { 1436 n=size(V0); 1437 module U=std(V0[n]); 1438 for(i=n1;i>=1;i) 1439 { 1440 V0[i]=simplify(reduce(V0[i],U),1); 1441 if(i>=2) 1442 { 1443 U=std(U+V0[i]); 1444 } 1445 } 1446 } 1447 1448 list l; 1449 if(k>0) 1450 { 1451 l=e0,m0; 1452 if(size(V0)>0) 1453 { 1454 l[3]=V0; 1455 } 1456 } 1410 1457 return(l); 1411 1458 } … … 1425 1472 int l[2][i]: multiplicity of eigenvalue l[1][i] 1426 1473 @end format 1427 EXAMPLE: example eigenval ; shows examples1474 EXAMPLE: example eigenvalues; shows examples 1428 1475 " 1429 1476 { 1430 1477 M=jet(hessenberg(M),0); 1431 1478 int n=ncols(M); 1479 int k; 1480 ideal e; 1481 intvec m; 1432 1482 number e0; 1433 1483 intvec v; 1434 list l ,f;1435 int i,j ,k;1484 list l; 1485 int i,j; 1436 1486 j=1; 1437 1487 while(j<=n) … … 1454 1504 if(size(v)==1) 1455 1505 { 1456 l=spappend(l,number(M[v,v]),1); 1506 k++; 1507 e[k]=M[v,v]; 1508 m[k]=1; 1457 1509 } 1458 1510 else 1459 1511 { 1460 f=factorize(det(submat(M,v,v)var(1)));1461 for(i=size( f[1]);i>=1;i)1512 l=factorize(det(submat(M,v,v)var(1))); 1513 for(i=size(l[1]);i>=1;i) 1462 1514 { 1463 e0=number(jet( f[1][i]/var(1),0));1515 e0=number(jet(l[1][i]/var(1),0)); 1464 1516 if(e0!=0) 1465 1517 { 1466 l=spappend(l,number((e0*var(1)f[1][i])/e0),f[2][i]); 1518 k++; 1519 e[k]=(e0*var(1)l[1][i])/e0; 1520 m[k]=l[2][i]; 1467 1521 } 1468 1522 } 1469 1523 } 1470 1524 } 1471 return(spn ormalize(l[1],l[2]));1525 return(spnf(e,m)); 1472 1526 } 1473 1527 example … … 1600 1654 def e,m=#[1..2]; 1601 1655 1602 int i; 1603 for(i=1;i<=ncols(e);i++) 1656 for(int i=1;i<=ncols(e);i++) 1604 1657 { 1605 1658 if(deg(e[i]>0)) … … 1610 1663 } 1611 1664 1612 int j,k,l ;1665 int j,k,l,n; 1613 1666 matrix N0,N1; 1614 1667 module K0,K1; … … 1618 1671 intvec w; 1619 1672 1620 for(i= ncols(e);i>=1;i)1673 for(i=1;i<=ncols(e);i++) 1621 1674 { 1622 1675 N0=Me[i]*freemodule(ncols(M)); … … 1648 1701 for(j=l;j>=1;j) 1649 1702 { 1650 U=module(u)+U; 1651 w=2*jl1,w; 1703 U=U+module(u); 1704 n++; 1705 w[n]=2*jl1; 1652 1706 u=N0*u; 1653 1707 } … … 1655 1709 } 1656 1710 } 1657 w=w[1..size(w)1]; 1711 1658 1712 return(list(U,w)); 1659 1713 } … … 1696 1750 for(i=s[k];i>=2;i) 1697 1751 { 1698 J[j ,j+1]=1;1752 J[j+1,j]=1; 1699 1753 j++; 1700 1754 J[j,j]=e[k];
Note: See TracChangeset
for help on using the changeset viewer.