Changeset 337919 in git
- Timestamp:
- Jan 11, 2001, 6:13:00 PM (23 years ago)
- Branches:
- (u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')
- Children:
- 167e413179ab275fe34d1eb07187736261e4ec66
- Parents:
- 2dd32e6b8afefafb14da9f8d5f5e44b849c103a4
- Location:
- Singular/LIB
- Files:
-
- 4 edited
-
gaussman.lib (modified) (4 diffs)
-
linalg.lib (modified) (71 diffs)
-
mondromy.lib (modified) (4 diffs)
-
rinvar.lib (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/gaussman.lib
r2dd32e r337919 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: gaussman.lib,v 1.2 4 2001-01-10 17:42:26 mschulzeExp $";2 version="$Id: gaussman.lib,v 1.25 2001-01-11 17:12:58 Singular Exp $"; 3 3 category="Singularities"; 4 4 … … 18 18 gamma4(...); Hertling's gamma4 invariant 19 19 20 SEE ALSO: mondromy .lib, spectrum.lib, jordan.lib20 SEE ALSO: mondromy_lib, spectrum_lib, jordan_lib 21 21 22 22 KEYWORDS: singularities; Gauss-Manin connection; monodromy; spectrum; … … 250 250 intvec b; 251 251 for(i=ncols(e);i>=1;i--) 252 252 { 253 253 b[i]=sum(l[2][i]); 254 254 } … … 774 774 int l[2][i]: multiplicity of spectral number l[1][i] 775 775 list l[3]: 776 module l[3][i]: vector space basis of the l[1][i]-th graded part 776 module l[3][i]: vector space basis of the l[1][i]-th graded part 777 777 of the V-filtration on the Jacobian algebra 778 778 in terms of l[4] -
Singular/LIB/linalg.lib
r2dd32e r337919 1 //GMG last modified: 04/25/2000 2 ////////////////////////////////////////////////////////////////////////////// 3 version="$Id: linalg.lib,v 1. 9 2001-01-07 02:06:39 greuelExp $";1 //GMG last modified: 04/25/2000 2 ////////////////////////////////////////////////////////////////////////////// 3 version="$Id: linalg.lib,v 1.10 2001-01-11 17:12:59 Singular Exp $"; 4 4 category="Linear Algebra"; 5 5 info=" 6 LIBRARY: linalg.lib 6 LIBRARY: linalg.lib Algorithmic Linear Algebra 7 7 AUTHORS: Ivor Saynisch (ivs@math.tu-cottbus.de) 8 8 @* Mathias Schulze (mschulze@mathematik.uni-kl.de) 9 9 10 10 PROCEDURES: 11 inverse(A); 11 inverse(A); matrix, the inverse of A 12 12 inverse_B(A); list(matrix Inv,poly p),Inv*A=p*En ( using busadj(A) ) 13 13 inverse_L(A); list(matrix Inv,poly p),Inv*A=p*En ( using lift ) … … 16 16 diag_test(A); test whether A can be diagnolized 17 17 busadj(A); coefficients of Adj(E*t-A) and coefficients of det(E*t-A) 18 charpoly(A,v); characteristic polynomial of A ( using busadj(A) ) 18 charpoly(A,v); characteristic polynomial of A ( using busadj(A) ) 19 19 adjoint(A); adjoint of A ( using busadj(A) ) 20 20 det_B(A); determinant of A ( using busadj(A) ) … … 61 61 kill @R; 62 62 setring BR; 63 return(1); 63 return(1); 64 64 } 65 65 ////////////////////////////////////////////////////////////////////////////// … … 70 70 RETURN: matrix 71 71 " 72 { 72 { 73 73 module m=module(A); 74 74 75 75 if(A[i,i]==0){ 76 76 m[i]=m[i]+m[j]; … … 79 79 m=module(transpose(matrix(m))); 80 80 } 81 81 82 82 A=matrix(m); 83 83 m[j]=m[j]-(A[i,j]/A[i,i])*m[i]; … … 85 85 m[j]=m[j]-(A[i,j]/A[i,i])*m[i]; 86 86 m=module(transpose(matrix(m))); 87 87 88 88 return(matrix(m)); 89 89 } … … 94 94 { 95 95 int k; 96 if (nrows(v2)>nrows(v1)) { k=nrows(v2); } else { k=nrows(v1); } 96 if (nrows(v2)>nrows(v1)) { k=nrows(v2); } else { k=nrows(v1); } 97 97 return ((transpose(matrix(v1,k,1))*matrix(v2,k,1))[1,1]); 98 98 } … … 104 104 proc inverse(matrix A, list #) 105 105 "USAGE: inverse(A [,opt]); A a square matrix, opt integer 106 RETURN: 107 @format 106 RETURN: 107 @format 108 108 a matrix: 109 109 - the inverse matrix of A, if A is invertible; 110 110 - the 1x1 0-matrix if A is not invertible (in the polynomial ring!). 111 111 There are the following options: 112 - opt=0 or not given: heuristically best option from below 112 - opt=0 or not given: heuristically best option from below 113 113 - opt=1 : apply std to (transpose(E,A)), ordering (C,dp). 114 114 - opt=2 : apply interred (transpose(E,A)), ordering (C,dp). 115 - opt=3 : apply lift(A,E), ordering (C,dp). 115 - opt=3 : apply lift(A,E), ordering (C,dp). 116 116 @end format 117 117 NOTE: parameters and minpoly are allowed; opt=2 is only correct for … … 123 123 //--------------------------- initialization and check ------------------------ 124 124 int ii,u,i,opt; 125 matrix invA; 125 matrix invA; 126 126 int db = printlevel-voice+3; //used for comments 127 127 def R=basering; … … 149 149 {opt = 1;} 150 150 else 151 {opt = 1;} 151 {opt = 1;} 152 152 } 153 153 else {opt = 1;} … … 177 177 if( opt==1 or opt==2 ) 178 178 { 179 option(redSB); 179 option(redSB); 180 180 module B = freemodule(n),M; 181 181 if(opt == 2) … … 186 186 if(opt == 1) 187 187 { 188 module D = std(transpose(B)); 188 module D = std(transpose(B)); 189 189 D = transpose(simplify(D,1)); 190 190 } … … 192 192 module D1 = D[n+1..2*n]; 193 193 //----------------------- check if matrix is invertible ---------------------- 194 for (ii=1; ii<=n; ii++) 195 { 194 for (ii=1; ii<=n; ii++) 195 { 196 196 if ( D1[ii] != gen(ii) ) 197 197 { … … 225 225 print(inverse(A));""; 226 226 matrix B[3][3]= 227 y+1, x+y, y, 228 z, z+1, z, 227 y+1, x+y, y, 228 z, z+1, z, 229 229 y+z+2,x+y+z+2,y+z+1; 230 230 print(inverse(B)); … … 244 244 "// ** input is not a square matrix";; 245 245 return(A); 246 } 247 246 } 247 248 248 if(!const_mat(A)){ 249 249 "// ** input is not a constant matrix"; 250 250 return(A); 251 251 } 252 253 if(deg(std(A-transpose(A))[1])!=-1){ 252 253 if(deg(std(A-transpose(A))[1])!=-1){ 254 254 "// ** input is not a symmetric matrix"; 255 255 return(A); … … 261 261 } 262 262 } 263 263 264 264 return(A); 265 265 } … … 273 273 274 274 ////////////////////////////////////////////////////////////////////////////// 275 proc orthogonalize(matrix A) 275 proc orthogonalize(matrix A) 276 276 "USAGE: orthogonalize(A); A = constant matrix 277 277 RETURN: matrix, orthogonal basis of the colum space of A … … 286 286 matrix B; 287 287 return(B); 288 } 288 } 289 289 290 290 module B=module(interred(A)); 291 291 292 292 for(i=1;i<=n;i=i+1) { 293 293 for(j=1;j<i;j=j+1) { … … 297 297 } 298 298 } 299 299 300 300 return(matrix(B)); 301 301 } … … 310 310 //////////////////////////////////////////////////////////////////////////// 311 311 proc diag_test(matrix A) 312 "USAGE: 312 "USAGE: diag_test(A); A = const square matrix 313 313 RETURN: int, 1 if A is diagonalisable, 0 if not 314 314 -1 no statement is possible, since A does not split. … … 319 319 { 320 320 int i,j; 321 int n = nrows(A); 321 int n = nrows(A); 322 322 string mp = string(minpoly); 323 323 string cs = charstr(basering); 324 324 int np=0; 325 325 326 if(ncols(A) != n) { 326 if(ncols(A) != n) { 327 327 "// input is not a square matrix"; 328 328 return(-1); 329 } 329 } 330 330 331 331 if(!const_mat(A)){ 332 332 "// input is not a constant matrix"; 333 333 return(-1); 334 } 334 } 335 335 336 336 //Parameterring wegen factorize nicht erlaubt … … 339 339 } 340 340 if(np>0){ 341 342 343 } 344 345 //speichern des aktuellen Rings 341 "// rings with parameters not allowed"; 342 return(-1); 343 } 344 345 //speichern des aktuellen Rings 346 346 def BR=basering; 347 347 //setze R[t] 348 348 execute("ring rt=("+charstr(basering)+"),(@t,"+varstr(basering)+"),lp;"); 349 349 execute("minpoly="+mp+";"); 350 matrix A=imap(BR,A); 350 matrix A=imap(BR,A); 351 351 352 352 intvec z; 353 353 intvec s; 354 poly X; //characteristisches Polynom 354 poly X; //characteristisches Polynom 355 355 poly dXdt; //Ableitung von X nach t 356 ideal g; 357 poly b; 356 ideal g; //ggT(X,dXdt) 357 poly b; //Komponente der Busadjunkten-Matrix 358 358 matrix E[n][n]; //Einheitsmatrix 359 359 360 360 E=E+1; 361 361 A=E*@t-A; 362 362 X=det(A); 363 363 364 matrix Xfactors=matrix(factorize(X,1)); 364 matrix Xfactors=matrix(factorize(X,1)); //zerfaellt die Matrtix ? 365 365 int nf=ncols(Xfactors); 366 366 … … 371 371 return(-1); 372 372 } 373 } 373 } 374 374 375 375 dXdt=diff(X,@t); … … 377 377 378 378 //Busadjunkte 379 z=2..n; 379 z=2..n; 380 380 for(i=1;i<=n;i++){ 381 s=2..n; 381 s=2..n; 382 382 for(j=1;j<=n;j++){ 383 383 b=det(submat(A,z,s)); 384 384 385 385 if(0!=reduce(b,g)){ 386 387 388 389 } 390 386 //" matrix not diagonalizable"; 387 setring BR; 388 return(0); 389 } 390 391 391 s[j]=j; 392 392 } 393 393 z[i]=i; 394 394 } 395 396 //"Die Matrix ist diagonalisierbar"; 395 396 //"Die Matrix ist diagonalisierbar"; 397 397 setring BR; 398 398 return(1); … … 410 410 "USAGE: busadj(A); A = square matrix (of size nxn) 411 411 RETURN: list L: 412 @format 413 L[1] contains the (n+1) coefficients of the characteristic 412 @format 413 L[1] contains the (n+1) coefficients of the characteristic 414 414 polynomial X of A, i.e. 415 415 X = L[1][1]+..+L[1][k]*t^(k-1)+..+(L[1][n+1])*t^n 416 416 L[2] contains the n (nxn)-matrices Hk which are the coefficients of 417 417 the busadjoint bA = adjoint(E*t-A) of A, i.e. 418 bA = (Hn-1)*t^(n-1)+...+Hk*t^k+...+H0, ( Hk=L[2][k+1] ) 418 bA = (Hn-1)*t^(n-1)+...+Hk*t^k+...+H0, ( Hk=L[2][k+1] ) 419 419 @end format 420 420 EXAMPLE: example busadj; shows an example" … … 428 428 poly a; 429 429 430 if(ncols(A) != n) { 430 if(ncols(A) != n) { 431 431 "input is not a square matrix"; 432 432 return(L); 433 } 433 } 434 434 435 435 bA = E; … … 457 457 list L = busadj(A); 458 458 poly X = L[1][1]+L[1][2]*t+L[1][3]*t2; X; 459 matrix bA[2][2] = L[2][1]+L[2][2]*t; 459 matrix bA[2][2] = L[2][1]+L[2][2]*t; 460 460 print(bA); //the busadjoint of A; 461 print(bA*(t*unitmat(2)-A)); 461 print(bA*(t*unitmat(2)-A)); 462 462 } 463 463 … … 465 465 proc charpoly(matrix A, list #) 466 466 "USAGE: charpoly(A[,v]); A square matrix, v string, name of a variable 467 RETURN: poly, the characteristic polynomial det(E*v-A) 467 RETURN: poly, the characteristic polynomial det(E*v-A) 468 468 (default: v=name of last variable) 469 469 NOTE: A must be independent of the variable v. The computation uses det. … … 477 477 poly X; 478 478 if(ncols(A) != n) 479 { 479 { 480 480 "// input is not a square matrix"; 481 481 return(X); 482 } 482 } 483 483 //---------------------- test for correct variable ------------------------- 484 if( size(#)==0 ){ 485 #[1] = varstr(z); 484 if( size(#)==0 ){ 485 #[1] = varstr(z); 486 486 } 487 487 if( typeof(#[1]) == "string") { v = #[1]; } … … 490 490 "// 2nd argument must be a name of a variable not contained in the matrix"; 491 491 return(X); 492 } 492 } 493 493 j=-1; 494 494 for(i=1; i<=z; i++) … … 536 536 string mp = string(minpoly); 537 537 538 if(ncols(A) != n){ 538 if(ncols(A) != n){ 539 539 "// input is not a square matrix"; 540 540 return(X); 541 } 541 } 542 542 543 543 //test for correct variable 544 if( size(#)==0 ){ 545 #[1] = varstr(nvars(BR)); 544 if( size(#)==0 ){ 545 #[1] = varstr(nvars(BR)); 546 546 } 547 547 if( typeof(#[1]) == "string"){ … … 551 551 "// 2nd argument must be a name of a variable not contained in the matrix"; 552 552 return(X); 553 } 554 553 } 554 555 555 j=-1; 556 556 for(i=1; i<=nvars(BR); i++){ … … 561 561 return(X); 562 562 } 563 563 564 564 //var can not be in A 565 565 s="Wp("; … … 590 590 execute("X=X+L[1][i]*"+v+"^"+string(i-1)+";"); 591 591 } 592 593 return(X); 592 593 return(X); 594 594 } 595 595 example … … 612 612 list L; 613 613 614 if(ncols(A) != n) { 614 if(ncols(A) != n) { 615 615 "// input is not a square matrix"; 616 616 return(Adj); 617 } 618 617 } 618 619 619 L = busadj(A); 620 620 Adj= (-1)^(n-1)*L[2][1]; 621 621 return(Adj); 622 622 623 623 } 624 624 example … … 636 636 "USAGE: inverse_B(A); A = square matrix 637 637 RETURN: list Inv with 638 - Inv[1] = matrix I and 639 - Inv[2] = poly p 640 such that I*A = unitmat(n)*p; 638 - Inv[1] = matrix I and 639 - Inv[2] = poly p 640 such that I*A = unitmat(n)*p; 641 641 NOTE: p=1 if 1/det(A) is computable and p=det(A) if not; 642 642 the computation uses busadj. … … 650 650 list L; 651 651 list Inv; 652 653 if(ncols(A) != n) { 652 653 if(ncols(A) != n) { 654 654 "input is not a square matrix"; 655 655 return(I); 656 } 657 656 } 657 658 658 L=busadj(A); 659 I=module(-L[2][1]); //+-Adj(A) 660 661 if(reduce(1,std(L[1][1]))==0){ 662 I=I*lift(L[1][1],1)[1][1]; 659 I=module(-L[2][1]); //+-Adj(A) 660 661 if(reduce(1,std(L[1][1]))==0){ 662 I=I*lift(L[1][1],1)[1][1]; 663 663 factor=1; 664 664 } … … 682 682 ////////////////////////////////////////////////////////////////////////////// 683 683 proc det_B(matrix A) 684 "USAGE: det_B(A); A any matrix 684 "USAGE: det_B(A); A any matrix 685 685 RETURN: returns the determinant of A 686 686 NOTE: the computation uses the busadj algorithm 687 687 EXAMPLE: example det_B; shows an example" 688 688 { 689 int n=nrows(A); 689 int n=nrows(A); 690 690 list L; 691 691 … … 696 696 } 697 697 example 698 { "EXAMPLE"; echo=2; 698 { "EXAMPLE"; echo=2; 699 699 ring r=0,(x),dp; 700 700 matrix A[10][10]=random(2,10,10)+unitmat(10)*x; 701 701 print(A); 702 det_B(A); 702 det_B(A); 703 703 } 704 704 … … 707 707 "USAGE: inverse_L(A); A = square matrix 708 708 RETURN: list Inv representing a left inverse of A, i.e 709 - Inv[1] = matrix I and 709 - Inv[1] = matrix I and 710 710 - Inv[2] = poly p 711 such that I*A = unitmat(n)*p; 711 such that I*A = unitmat(n)*p; 712 712 NOTE: p=1 if 1/det(A) is computable and p=det(A) if not; 713 713 the computation computes first det(A) and then uses lift 714 SEE ALSO: inverse, inverse_B _714 SEE ALSO: inverse, inverse_B 715 715 EXAMPLE: example inverse_L; shows an example" 716 716 { … … 735 735 // test if 1/det(A) exists 736 736 if(reduce(1,std(d))!=0){ E=E*d;} 737 737 738 738 I=lift(A,E); 739 if(I==unitmat(n)-unitmat(n)){ //catch error in lift 739 if(I==unitmat(n)-unitmat(n)){ //catch error in lift 740 740 "// matrix is not invertible"; 741 741 return(Inv); … … 745 745 Inv=insert(Inv,factor); 746 746 Inv=insert(Inv,I); 747 747 748 748 return(Inv); 749 749 } … … 761 761 ////////////////////////////////////////////////////////////////////////////// 762 762 proc gaussred(matrix A) 763 "USAGE: gaussred(A); A any constant matrix 763 "USAGE: gaussred(A); A any constant matrix 764 764 RETURN: list Z: Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A) 765 765 gives a row reduced matrix S, a permutation matrix P and a 766 normalized lower triangular matrix U, with P*A=U*S 766 normalized lower triangular matrix U, with P*A=U*S 767 767 NOTE: This procedure is designed for teaching purposes mainly. 768 The straight forward implementation in the interpreted library 769 is not very efficient (no standard basis computation). 768 The straight forward implementation in the interpreted library 769 is not very efficient (no standard basis computation). 770 770 EXAMPLE: example gaussred; shows an example" 771 771 { … … 787 787 788 788 for(i=1;i<=mr;i=i+1){ 789 if((i+k)>m){break}; 790 789 if((i+k)>m){break}; 790 791 791 //Test: Diagonalelement=0 792 792 if(A[i,i+k]==0){ 793 793 jp=i;pivo=0; 794 794 for(j=i+1;j<=n;j=j+1){ 795 796 795 c=abs(A[j,i+k]); 796 if(pivo<c){ pivo=c;jp=j;} 797 797 } 798 798 if(jp != i){ //Zeilentausch 799 800 c=A[i,j]; 801 802 803 804 805 c=P[i,j]; 806 807 P[jp,j]=c; 808 799 for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter) 800 c=A[i,j]; 801 A[i,j]=A[jp,j]; 802 A[jp,j]=c; 803 } 804 for(j=1;j<=n;j=j+1){ //Zeilentausch in P 805 c=P[i,j]; 806 P[i,j]=P[jp,j]; 807 P[jp,j]=c; 808 } 809 809 } 810 810 if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe ! 811 811 } //i sollte im naechsten Lauf nicht erhoeht sein 812 813 //Eliminationsschritt 814 for(j=i+1;j<=n;j=j+1){ 812 813 //Eliminationsschritt 814 for(j=i+1;j<=n;j=j+1){ 815 815 c=A[j,i+k]/A[i,i+k]; 816 816 for(l=i+k+1;l<=m;l=l+1){ 817 817 A[j,l]=A[j,l]-A[i,l]*c; 818 818 } 819 819 A[j,i+k]=0; // nur wichtig falls k>0 ist 820 820 A[j,i]=c; // bildet U 821 821 } 822 rang=i; 823 } 824 822 rang=i; 823 } 824 825 825 for(i=1;i<=mr;i=i+1){ 826 826 for(j=i+1;j<=n;j=j+1){ … … 829 829 } 830 830 } 831 831 832 832 Z=insert(Z,rang); 833 833 Z=insert(Z,A); 834 834 Z=insert(Z,U); 835 835 Z=insert(Z,P); 836 836 837 837 return(Z); 838 838 } … … 857 857 gives n row reduced matrix S, a permutation matrix P and a 858 858 normalized lower triangular matrix U, with P*A=U*S 859 NOTE: with row pivoting 859 NOTE: with row pivoting 860 860 EXAMPLE: example gaussred_pivot; shows an example" 861 861 { … … 877 877 878 878 for(i=1;i<=mr;i=i+1){ 879 if((i+k)>m){break}; 880 879 if((i+k)>m){break}; 880 881 881 //Pivotisierung 882 882 pivo=abs(A[i,i+k]);jp=i; … … 887 887 if(jp != i){ //Zeilentausch 888 888 for(j=1;j<=m;j=j+1){ //Zeilentausch in A (und U) (i-te mit jp-ter) 889 c=A[i,j]; 890 891 889 c=A[i,j]; 890 A[i,j]=A[jp,j]; 891 A[jp,j]=c; 892 892 } 893 893 for(j=1;j<=n;j=j+1){ //Zeilentausch in P 894 c=P[i,j]; 895 896 P[jp,j]=c; 894 c=P[i,j]; 895 P[i,j]=P[jp,j]; 896 P[jp,j]=c; 897 897 } 898 898 } 899 899 if(pivo==0){k++;continue;} //eine von selbst auftauchende Stufe ! 900 900 //i sollte im naechsten Lauf nicht erhoeht sein 901 //Eliminationsschritt 902 for(j=i+1;j<=n;j=j+1){ 901 //Eliminationsschritt 902 for(j=i+1;j<=n;j=j+1){ 903 903 c=A[j,i+k]/A[i,i+k]; 904 904 for(l=i+k+1;l<=m;l=l+1){ 905 905 A[j,l]=A[j,l]-A[i,l]*c; 906 906 } 907 907 A[j,i+k]=0; // nur wichtig falls k>0 ist 908 908 A[j,i]=c; // bildet U 909 909 } 910 rang=i; 911 } 912 910 rang=i; 911 } 912 913 913 for(i=1;i<=mr;i=i+1){ 914 914 for(j=i+1;j<=n;j=j+1){ … … 917 917 } 918 918 } 919 919 920 920 Z=insert(Z,rang); 921 921 Z=insert(Z,A); 922 922 Z=insert(Z,U); 923 923 Z=insert(Z,P); 924 924 925 925 return(Z); 926 926 } … … 966 966 proc mat_rk(matrix A) 967 967 "USAGE: mat_rk(A); A any constant matrix 968 RETURN: int, rank of A 968 RETURN: int, rank of A 969 969 EXAMPLE: example mat_rk; shows an example" 970 970 { … … 988 988 "USAGE: U_D_O(A); constant invertible matrix A 989 989 RETURN: list Z: Z[1]=P , Z[2]=U , Z[3]=D , Z[4]=O 990 gives a permutation matrix P, 990 gives a permutation matrix P, 991 991 a normalized lower triangular matrix U , 992 a diagonal matrix D, and 992 a diagonal matrix D, and 993 993 a normalized upper triangular matrix O 994 994 with P*A=U*D*O … … 1010 1010 return(Z); 1011 1011 } 1012 1012 1013 1013 L=gaussred(A); 1014 1014 1015 1015 if(L[4]!=n){ 1016 1016 "// input is not an invertible matrix"; 1017 Z=insert(Z,-1); //hint for calling procedures 1017 Z=insert(Z,-1); //hint for calling procedures 1018 1018 return(Z); 1019 1019 } … … 1037 1037 { "EXAMPLE";echo=2; 1038 1038 ring r = 0,(x),dp; 1039 matrix A[5][5] = 10, 4, 0, -9, 8, 1040 -3, 6, -6, -4, 9, 1039 matrix A[5][5] = 10, 4, 0, -9, 8, 1040 -3, 6, -6, -4, 9, 1041 1041 0, 3, -1, -9, -8, 1042 1042 -4,-2, -6, -10,10, 1043 1043 -9, 5, -1, -6, 5; 1044 list Z = U_D_O(A); //construct P,U,D,O s.t. P*A=U*D*O 1044 list Z = U_D_O(A); //construct P,U,D,O s.t. P*A=U*D*O 1045 1045 print(Z[1]); //P 1046 1046 print(Z[2]); //U … … 1053 1053 ////////////////////////////////////////////////////////////////////////////// 1054 1054 proc pos_def(matrix A) 1055 "USAGE: pos_def(A); A = constant, symmetric square matrix 1056 RETURN: int: 1057 1 if A is positive definit , 1058 0 if not, 1059 -1 if unknown 1055 "USAGE: pos_def(A); A = constant, symmetric square matrix 1056 RETURN: int: 1057 1 if A is positive definit , 1058 0 if not, 1059 -1 if unknown 1060 1060 EXAMPLE: example pos_def; shows an example" 1061 1061 { … … 1072 1072 "// input is not a constant matrix"; 1073 1073 return(-1); 1074 } 1075 if(deg(std(A-transpose(A))[1])!=-1){ 1074 } 1075 if(deg(std(A-transpose(A))[1])!=-1){ 1076 1076 "// input is not a hermitian (symmetric) matrix"; 1077 1077 return(-1); 1078 1078 } 1079 1079 1080 1080 Z=U_D_O(A); 1081 1081 … … 1085 1085 1086 1086 H=Z[1]; 1087 //es fand Zeilentausch statt: also nicht positiv definit 1087 //es fand Zeilentausch statt: also nicht positiv definit 1088 1088 if(deg(std(H-unitmat(n))[1])!=-1){ 1089 1089 return(0); 1090 1090 } 1091 1091 1092 1092 H=Z[3]; 1093 1093 1094 1094 for(j=1;j<=n;j=j+1){ 1095 if(H[j,j]<=0){ 1095 if(H[j,j]<=0){ 1096 1096 return(0); 1097 1097 } //eigenvalue<=0, not pos.definit … … 1105 1105 matrix A[5][5] = 20, 4, 0, -9, 8, 1106 1106 4, 12, -6, -4, 9, 1107 0, -6, -2, -9, -8, 1108 -9, -4, -9, -20, 10, 1107 0, -6, -2, -9, -8, 1108 -9, -4, -9, -20, 10, 1109 1109 8, 9, -8, 10, 10; 1110 1110 pos_def(A); … … 1118 1118 proc linsolve(matrix A, matrix b) 1119 1119 "USAGE: linsolve(A,b); A a constant nxm-matrix, b a constant nx1-matrix 1120 RETURN: a 1xm matrix X, solution of inhomogeneous linear system A*X = b 1120 RETURN: a 1xm matrix X, solution of inhomogeneous linear system A*X = b 1121 1121 return the 0-matrix if system is not solvable 1122 1122 NOTE: uses gaussred … … 1131 1131 matrix Ab[n][m+1]; 1132 1132 matrix X[m][1]; 1133 1133 1134 1134 if(ncols(b)!=1){ 1135 1135 "// right hand side b is not a nx1 matrix"; … … 1140 1140 "// input hand is not a constant matrix"; 1141 1141 return(X); 1142 } 1143 1142 } 1143 1144 1144 if(n_b>n){ 1145 1145 for(i=n; i<=n_b; i++){ 1146 1146 if(b[i,1]!=0){ 1147 1148 1149 } 1150 } 1151 } 1152 1153 if(n_b<n){ 1147 "// right hand side b not in Image(A)"; 1148 return X; 1149 } 1150 } 1151 } 1152 1153 if(n_b<n){ 1154 1154 matrix copy[n_b][1]=b; 1155 1155 matrix b[n][1]=0; … … 1158 1158 } 1159 1159 } 1160 1160 1161 1161 r=mat_rk(A); 1162 1162 1163 1163 //1. b constant vector 1164 if(const_mat(b)){ 1164 if(const_mat(b)){ 1165 1165 //extend A with b 1166 1166 for(i=1; i<=n; i++){ … … 1170 1170 Ab[i,m+1]=b[i,1]; 1171 1171 } 1172 1172 1173 1173 //Gauss reduction 1174 1174 Z = gaussred(Ab); 1175 1175 Ab = Z[3]; //normal form 1176 rc = Z[4]; //rank(Ab) 1176 rc = Z[4]; //rank(Ab) 1177 1177 //print(Ab); 1178 1178 1179 1179 if(r<rc){ 1180 1181 return(X); 1182 } 1183 k=m; 1180 "// no solution"; 1181 return(X); 1182 } 1183 k=m; 1184 1184 for(i=r;i>=1;i=i-1){ 1185 1186 j=1; 1187 while(Ab[i,j]==0){j=j+1;}// suche Ecke 1188 1185 1186 j=1; 1187 while(Ab[i,j]==0){j=j+1;}// suche Ecke 1188 1189 1189 for(;k>j;k=k-1){ X[k]=0;}//springe zur Ecke 1190 1190 1191 1191 1192 1192 c=Ab[i,m+1]; //i-te Komponene von b 1193 1193 for(j=m;j>k;j=j-1){ 1194 1194 c=c-X[j,1]*Ab[i,j]; 1195 1195 } 1196 1196 if(Ab[i,k]==0){ 1197 1198 } 1199 else{ 1200 1197 X[k,1]=1; //willkuerlich 1198 } 1199 else{ 1200 X[k,1]=c/Ab[i,k]; 1201 1201 } 1202 1202 k=k-1; 1203 1203 if(k==0){break;} 1204 1204 } 1205 1206 1205 1206 1207 1207 }//endif (const b) 1208 1208 else{ //b not constant 1209 1209 "// !not implemented!"; 1210 1210 1211 1211 } 1212 1212 … … 1217 1217 ring r=0,(x),dp; 1218 1218 matrix A[3][2] = -4,-6, 1219 2, 3, 1219 2, 3, 1220 1220 -5, 7; 1221 1221 matrix b[3][1] = 10, … … 1277 1277 ASSUME: The eigenvalues of M are in the coefficient field. 1278 1278 RETURN: The procedure returns a list jd with 3 entries: 1279 @format 1279 @format 1280 1280 - jd[1] ideal, eigenvalues of M, 1281 - jd[2] list of intvecs, jd[2][i][j] size of j-th Jordan block with 1281 - jd[2] list of intvecs, jd[2][i][j] size of j-th Jordan block with 1282 1282 eigenvalue jd[1][i], and 1283 1283 - jd[3] a matrix, jd[3]^(-1)*M*jd[3] in Jordan normal form. … … 1745 1745 system("--random", 12345678); 1746 1746 int n = 70; 1747 matrix m = sparsemat(n,n,50,100); 1747 matrix m = sparsemat(n,n,50,100); 1748 1748 option(prot,mem); 1749 1749 1750 1750 int t=timer; 1751 matrix im = inverse(m,1)[1]; 1751 matrix im = inverse(m,1)[1]; 1752 1752 timer-t; 1753 1753 print(im*m); 1754 //list l0 = watchdog(100,"inverse("+"m"+",3)"); 1754 //list l0 = watchdog(100,"inverse("+"m"+",3)"); 1755 1755 //bricht bei 100 sec ab und gibt l0[1]: string Killed zurueck 1756 1756 1757 //inverse(m,1): std 5sec 5,5 MB 1757 //inverse(m,1): std 5sec 5,5 MB 1758 1758 //inverse(m,2): interred 12sec 1759 1759 //inverse(m,2): lift nach 180 sec 13MB abgebrochen 1760 1760 //n=60: linalgorig: 3 linalg: 5 1761 //n=70: linalgorig: 6,7 linalg: 11,12 1762 // aber linalgorig rechnet falsch! 1761 //n=70: linalgorig: 6,7 linalg: 11,12 1762 // aber linalgorig rechnet falsch! 1763 1763 1764 1764 2. Sparse poly Matrizen … … 1771 1771 1772 1772 int t=timer; 1773 matrix im = inverse(m); 1773 matrix im = inverse(m); 1774 1774 timer-t; 1775 1775 print(im*m); … … 1790 1790 1791 1791 int t=timer; 1792 matrix im = inverse(m); 1792 matrix im = inverse(m); 1793 1793 timer-t; 1794 1794 print(im*m); … … 1807 1807 1808 1808 int t=timer; 1809 matrix im = inverse(m,3); 1809 matrix im = inverse(m,3); 1810 1810 timer-t; 1811 1811 print(im*m); … … 1818 1818 int n =3; 1819 1819 ring r= 0,(x,y,z),(C,dp); 1820 matrix A=triagmatrix(n,n,1,0,0,50,2); 1821 intmat B=sparsetriag(n,n,20,1); 1820 matrix A=triagmatrix(n,n,1,0,0,50,2); 1821 intmat B=sparsetriag(n,n,20,1); 1822 1822 matrix M = A*transpose(B); 1823 1823 M=M*transpose(M); … … 1835 1835 //eifacheres Gegenbeispiel: 1836 1836 matrix M = 1837 9yz+3y+3z+2, 9y2+6y+1, 1837 9yz+3y+3z+2, 9y2+6y+1, 1838 1838 9xyz+3xy+3xz-9z2+2x-6z-1,9xy2+6xy-9yz+x-3y-3z 1839 1839 //det M=1, inverse(M,2); ->// ** matrix is not invertible … … 1842 1842 5. charpoly: 1843 1843 ----------- 1844 //ring rp=(0,A,B,C),(x),dp; 1845 ring r=0,(A,B,C,x),dp; 1844 //ring rp=(0,A,B,C),(x),dp; 1845 ring r=0,(A,B,C,x),dp; 1846 1846 matrix m[12][12]= 1847 AC,BC,-3BC,0,-A2+B2,-3AC+1,B2, B2, 1, 0, -C2+1,0, 1848 1, 1, 2C, 0,0, B, -A, -4C, 2A+1,0, 0, 0, 1849 0, 0, 0, 1,0, 2C+1, -4C+1,-A, B+1, 0, B+1, 3B, 1847 AC,BC,-3BC,0,-A2+B2,-3AC+1,B2, B2, 1, 0, -C2+1,0, 1848 1, 1, 2C, 0,0, B, -A, -4C, 2A+1,0, 0, 0, 1849 0, 0, 0, 1,0, 2C+1, -4C+1,-A, B+1, 0, B+1, 3B, 1850 1850 AB,B2,0, 1,0, 1, 0, 1, A, 0, 1, B+1, 1851 1, 0, 1, 0,0, 1, 0, -C2, 0, 1, 0, 1, 1852 0, 0, 2, 1,2A, 1, 0, 0, 0, 0, 1, 1, 1853 0, 1, 0, 1,1, 2, A, 3B+1,1, B2,1, 1, 1854 0, 1, 0, 1,1, 1, 1, 1, 2, 0, 0, 0, 1855 1, 0, 1, 0,0, 0, 1, 0, 1, 1, 0, 3, 1856 1, 3B,B2+1,0,0, 1, 0, 1, 0, 0, 1, 0, 1857 0, 0, 1, 0,0, 0, 0, 1, 0, 0, 0, 0, 1858 0, 1, 0, 1,1, 3, 3B+1, 0, 1, 1, 1, 0; 1851 1, 0, 1, 0,0, 1, 0, -C2, 0, 1, 0, 1, 1852 0, 0, 2, 1,2A, 1, 0, 0, 0, 0, 1, 1, 1853 0, 1, 0, 1,1, 2, A, 3B+1,1, B2,1, 1, 1854 0, 1, 0, 1,1, 1, 1, 1, 2, 0, 0, 0, 1855 1, 0, 1, 0,0, 0, 1, 0, 1, 1, 0, 3, 1856 1, 3B,B2+1,0,0, 1, 0, 1, 0, 0, 1, 0, 1857 0, 0, 1, 0,0, 0, 0, 1, 0, 0, 0, 0, 1858 0, 1, 0, 1,1, 3, 3B+1, 0, 1, 1, 1, 0; 1859 1859 option(prot,mem); 1860 1860 -
Singular/LIB/mondromy.lib
r2dd32e r337919 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: mondromy.lib,v 1. 19 2000-12-22 15:30:39 greuelExp $";2 version="$Id: mondromy.lib,v 1.20 2001-01-11 17:12:59 Singular Exp $"; 3 3 category="Singularities"; 4 4 info=" … … 6 6 AUTHOR: Mathias Schulze, email: mschulze@mathematik.uni-kl.de 7 7 8 OVERVIEW: 9 A library to compute the monodromy of an isolated hypersurface singularity. 10 It uses an algorithm by Brieskorn (manuscripta math. 2 (1970), 103-161) to 11 compute a connection matrix of the meromorphic Gauss-Manin connection up to 12 arbitrarily high order, and an algorithm of Gerard and Levelt (Ann. Inst. 8 OVERVIEW: 9 A library to compute the monodromy of an isolated hypersurface singularity. 10 It uses an algorithm by Brieskorn (manuscripta math. 2 (1970), 103-161) to 11 compute a connection matrix of the meromorphic Gauss-Manin connection up to 12 arbitrarily high order, and an algorithm of Gerard and Levelt (Ann. Inst. 13 13 Fourier, Grenoble 23,1 (1973), pp. 157-195) to transform it to a simple pole. 14 14 … … 23 23 Brieskorn lattice 24 24 25 SEE ALSO: gaussman .lib25 SEE ALSO: gaussman_lib 26 26 "; 27 27 … … 179 179 v=v+ui; 180 180 } 181 v=jet(v,n)/u0; 181 v=jet(v,n)/u0; 182 182 dbprint(printlevel-voice+2,"//...inverse computed ["+string(timer-t)+ 183 183 " secs, "+string((memory(1)+1023)/1024)+" K]"); -
Singular/LIB/rinvar.lib
r2dd32e r337919 1 1 // Last change 10.12.2000 (TB) 2 2 /////////////////////////////////////////////////////////////////////////////// 3 version="$Id: rinvar.lib,v 1. 5 2000-12-22 14:43:04 greuelExp $";3 version="$Id: rinvar.lib,v 1.6 2001-01-11 17:13:00 Singular Exp $"; 4 4 category="Invariant theory"; 5 5 info=" 6 6 LIBRARY: rinvar.lib Invariant Rings of Reductive Groups 7 AUTHOR: Thomas Bayer, tbayer@in.tum.de 7 AUTHOR: Thomas Bayer, tbayer@in.tum.de 8 8 http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/ 9 9 Current Adress: Institut fuer Informatik, TU Muenchen 10 OVERVIEW: 11 Implementation based on Derksen's algorithm. Written in the frame of the 12 diploma thesis (advisor: Prof. Gert-Martin Greuel) 'Computations of moduli 10 OVERVIEW: 11 Implementation based on Derksen's algorithm. Written in the frame of the 12 diploma thesis (advisor: Prof. Gert-Martin Greuel) 'Computations of moduli 13 13 spaces of semiquasihomogenous singularities and an implementation in Singular' 14 14 … … 30 30 TransferIdeal(R,name,nA); transfer the ideal 'name' from R to basering 31 31 32 SEE ALSO: qhmoduli .lib, zeroset.lib32 SEE ALSO: qhmoduli_lib, zeroset_lib 33 33 "; 34 34
Note: See TracChangeset
for help on using the changeset viewer.
