Changeset 1418c4 in git
 Timestamp:
 Aug 1, 2001, 2:33:41 PM (22 years ago)
 Branches:
 (u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
 Children:
 4b6c751b1fadaa54d2622f7071de1b4113148c5f
 Parents:
 34a9eb1eb751113c0633d266eb927061823c9864
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/gaussman.lib
r34a9eb1 r1418c4 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: gaussman.lib,v 1.4 6 20010731 17:37:43mschulze Exp $";2 version="$Id: gaussman.lib,v 1.47 20010801 12:33:41 mschulze Exp $"; 3 3 category="Singularities"; 4 4 … … 359 359 static proc maxintdif(ideal e) 360 360 { 361 dbprint(printlevelvoice+2,"//gaussman::maxintdif");362 361 int i,j,id; 363 362 int mid=0; … … 383 382 static proc maxorddif(matrix H) 384 383 { 385 dbprint(printlevelvoice+2,"//gaussman::maxorddif");386 384 int i,j,d; 387 385 int d0,d1=1,1; … … 439 437 def G=gmsring(f,"s"); 440 438 setring G; 439 441 440 int mu,modm=ncols(gmsbasis),maxorddif(gmsbasis); 442 441 ideal r=gmspoly*gmsbasis; … … 450 449 { 451 450 k++; 452 dbprint(printlevelvoice+2,"// gaussman::monodromy:k="+string(k));453 dbprint(printlevelvoice+2,"// gaussman::monodromy: compute C");451 dbprint(printlevelvoice+2,"// k="+string(k)); 452 dbprint(printlevelvoice+2,"// compute matrix A of t"); 454 453 if(opt==0) 455 454 { … … 463 462 A0=A0+C; 464 463 464 dbprint(printlevelvoice+2,"// compute saturation of H''"); 465 465 H0=H; 466 dbprint(printlevelvoice+2,"//gaussman::monodromy: compute dH");467 466 dH=jet(module(A0*dH+s^2*diff(matrix(dH),s)),k); 468 467 H=H*s+dH; 469 470 dbprint(printlevelvoice+2,"//gaussman::monodromy: test dH==0");471 468 } 472 469 A0=A0k*s; 473 470 474 dbprint(printlevelvoice+2, 475 "//gaussman::monodromy: compute basis of saturation"); 471 dbprint(printlevelvoice+2,"// compute basis of saturation of H''"); 476 472 H=std(H0); 477 473 int modH=maxorddif(H)/deg(s); 478 dbprint(printlevelvoice+2,"// gaussman::monodromy:k="+string(modH+1));479 dbprint(printlevelvoice+2,"// gaussman::monodromy: compute C");474 dbprint(printlevelvoice+2,"// k="+string(modH+1)); 475 dbprint(printlevelvoice+2,"// compute matrix A of t"); 480 476 if(opt==0) 481 477 { … … 488 484 C,r=l[1..2]; 489 485 A0=A0+C; 490 dbprint(printlevelvoice+2, 491 "//gaussman::monodromy: compute A on saturation"); 486 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 492 487 l=division(H*s,A0*H+s^2*diff(matrix(H),s)); 493 488 matrix A=jet(l[1],l[2],0); 494 489 495 490 dbprint(printlevelvoice+2, 496 "//gaussman::monodromy: compute eigenvalues e and "+ 497 "multiplicities b of A"); 491 "// compute eigenvalues e with multiplicities m of A"); 498 492 l=eigenval(jet(A,0)); 499 ideal e=l[1]; 500 intvec b=l[2]; 501 dbprint(printlevelvoice+2,"//gaussman::monodromy: e="+string(e)); 502 dbprint(printlevelvoice+2,"//gaussman::monodromy: b="+string(b)); 493 def e,m=l[1..2]; 494 dbprint(printlevelvoice+2,"// e="+string(e)); 495 dbprint(printlevelvoice+2,"// m="+string(m)); 503 496 504 497 if(opt==0) … … 509 502 } 510 503 504 dbprint(printlevelvoice+2,"// compute maximal integer difference of e"); 511 505 int mide=maxintdif(e); 506 dbprint(printlevelvoice+2,"// mide="+string(mide)); 507 512 508 if(mide>0) 513 509 { 514 dbprint(printlevelvoice+2, 515 "//gaussman::monodromy: k="+string(modH+1+mide)); 516 dbprint(printlevelvoice+2,"//gaussman::monodromy: compute C"); 510 dbprint(printlevelvoice+2,"// k="+string(modH+1+mide)); 511 dbprint(printlevelvoice+2,"// compute matrix A of t"); 517 512 l=gmscoeffs(r,modH+1+mide,modH+1+mide); 518 513 C,r=l[1..2]; 519 514 A0=A0+C; 515 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 520 516 l=division(H*s,A0*H+s^2*diff(matrix(H),s)); 521 517 A=jet(l[1],l[2],mide); … … 541 537 for(j,k=ncols(e),mu;j>=1;j) 542 538 { 543 for(i= b[j];i>=1;i)539 for(i=m[j];i>=1;i) 544 540 { 545 541 ide[k]=ide[j]; … … 550 546 while(mide>0) 551 547 { 552 dbprint(printlevelvoice+2,"// gaussman::monodromy: mide="+string(mide));548 dbprint(printlevelvoice+2,"// transform basis to reduce mide by 1"); 553 549 554 550 A0=jet(A,0); … … 585 581 } 586 582 } 583 587 584 mide; 585 dbprint(printlevelvoice+2,"// mide="+string(mide)); 588 586 } 589 587 … … 639 637 640 638 def R=basering; 639 int n=nvars(R)1; 641 640 def G=gmsring(f,"s"); 642 641 setring G; 643 642 644 int n=nvars(R)1;645 643 int mu,modm=ncols(gmsbasis),maxorddif(gmsbasis); 646 644 ideal r=gmspoly*gmsbasis; … … 650 648 module H,dH=freemodule(mu),freemodule(mu); 651 649 module H0; 652 int sdH=1;653 650 int k=1; 654 651 int N=n+1; 655 652 656 while(s dH>0)653 while(size(reduce(H,std(H0*s)))>0) 657 654 { 658 655 k++; 659 dbprint(printlevelvoice+2,"// gaussman::vfilt:k="+string(k));660 dbprint(printlevelvoice+2,"// gaussman::vfilt: compute C");656 dbprint(printlevelvoice+2,"// k="+string(k)); 657 dbprint(printlevelvoice+2,"// compute matrix A of t"); 661 658 l=gmscoeffs(r,k); 662 659 C,r=l[1..2]; 663 660 A=A+C; 664 661 662 dbprint(printlevelvoice+2,"// compute saturation of H''"); 665 663 H0=H; 666 dbprint(printlevelvoice+2,"//gaussman::vfilt: compute dH");667 664 dH=jet(module(A*dH+s^2*diff(matrix(dH),s)),k); 668 665 H=H*s+dH; 669 670 dbprint(printlevelvoice+2,"//gaussman::vfilt: test dH==0");671 sdH=size(reduce(H,std(H0*s)));672 666 } 673 667 A=Ak*s; 674 668 675 dbprint(printlevelvoice+2,"// gaussman::vfilt: compute basis of saturation");676 H= minbase(H0);669 dbprint(printlevelvoice+2,"// compute basis of saturation of H''"); 670 H=std(H0); 677 671 int modH=maxorddif(H)/deg(s); 678 dbprint(printlevelvoice+2,"// gaussman::monodromy:k="+string(N+modH));679 dbprint(printlevelvoice+2,"// gaussman::monodromy: compute C");672 dbprint(printlevelvoice+2,"// k="+string(N+modH)); 673 dbprint(printlevelvoice+2,"// compute matrix A of t"); 680 674 l=gmscoeffs(r,N+modH,N+modH); 681 675 C,r=l[1..2]; 682 676 A=A+C; 683 677 684 dbprint(printlevelvoice+2,"// gaussman::vfilt: transform H0 to saturation");678 dbprint(printlevelvoice+2,"// transform H'' to saturation of H''"); 685 679 l=division(H,freemodule(mu)*s^k); 686 680 H0=jet(l[1],l[2],N1); 687 681 688 dbprint(printlevelvoice+2, 689 "//gaussman::vfilt: compute H0 as vector space V0"); 690 dbprint(printlevelvoice+2, 691 "//gaussman::vfilt: compute H1 as vector space V1"); 682 dbprint(printlevelvoice+2,"// compute vector spaces"); 692 683 poly p; 693 684 int i0,j0,i1,j1; … … 716 707 } 717 708 718 dbprint(printlevelvoice+2,"// gaussman::vfilt: compute A on saturation");709 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 719 710 l=division(H*s,A*H+s^2*diff(matrix(H),s)); 720 711 A=jet(l[1],l[2],N1); 721 712 722 dbprint(printlevelvoice+2,"// gaussman::vfilt:compute matrix M of A");713 dbprint(printlevelvoice+2,"// compute matrix M of A"); 723 714 matrix M[mu*N][mu*N]; 724 715 for(i0=mu;i0>=1;i0) … … 746 737 } 747 738 748 dbprint(printlevelvoice+2,"// gaussman::vfilt:compute eigenvalues eA of A");739 dbprint(printlevelvoice+2,"// compute eigenvalues eA of A"); 749 740 ideal eA=eigenval(jet(A,0))[1]; 750 dbprint(printlevelvoice+2,"// gaussman::vfilt:eA="+string(eA));751 752 dbprint(printlevelvoice+2,"// gaussman::vfilt:compute eigenvalues eM of M");741 dbprint(printlevelvoice+2,"// eA="+string(eA)); 742 743 dbprint(printlevelvoice+2,"// compute eigenvalues eM of M"); 753 744 ideal eM; 754 745 k=0; … … 776 767 } 777 768 } 778 dbprint(printlevelvoice+2,"//gaussman::vfilt: eM="+string(eM)); 779 780 dbprint(printlevelvoice+2, 781 "//gaussman::vfilt: compute Vfiltration on H0/H1"); 769 dbprint(printlevelvoice+2,"// eM="+string(eM)); 770 771 dbprint(printlevelvoice+2,"// compute Vfiltration on H''/sH''"); 782 772 ideal a; 783 773 k=0; … … 789 779 for(i=ncols(eM);number(eM[i])1>number(n1)/2;i) 790 780 { 791 dbprint(printlevelvoice+2, 792 "//gaussman::vfilt: compute V["+string(i)+"]"); 781 dbprint(printlevelvoice+2,"// compute V["+string(i)+"]"); 793 782 V1=module(V1)+syz(power(MeM[i],n+1)); 794 783 V[i]=interred(intersect(V1,V0)); … … 802 791 } 803 792 804 dbprint(printlevelvoice+2, 805 "//gaussman::vfilt: symmetry index found"); 793 dbprint(printlevelvoice+2,"// symmetry index found"); 806 794 int j=k; 807 795 808 796 if(number(eM[i])1==number(n1)/2) 809 797 { 810 dbprint(printlevelvoice+2, 811 "//gaussman::vfilt: compute V["+string(i)+"]"); 798 dbprint(printlevelvoice+2,"// compute V["+string(i)+"]"); 812 799 V1=module(V1)+syz(power(MeM[i],n+1)); 813 800 V[i]=interred(intersect(V1,V0)); … … 821 808 } 822 809 823 dbprint(printlevelvoice+2,"// gaussman::vfilt:apply symmetry");810 dbprint(printlevelvoice+2,"// apply symmetry"); 824 811 while(j>=1) 825 812 { … … 841 828 for(i=ncols(eM);i>=1;i) 842 829 { 843 dbprint(printlevelvoice+2, 844 "//gaussman::vfilt: compute V["+string(i)+"]"); 830 dbprint(printlevelvoice+2,"// compute V["+string(i)+"]"); 845 831 V1=module(V1)+syz(power(MeM[i],n+1)); 846 832 V[i]=interred(intersect(V1,V0)); … … 852 838 k++; 853 839 a[k]=eM[i]1; 854 dbprint(printlevelvoice+2,855 "//gaussman::vfilt: transform to V0");856 840 v[k]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1]; 857 841 } … … 873 857 a[j]=eM[i]1; 874 858 v[k]=v[j]; 875 dbprint(printlevelvoice+2,876 "//gaussman::vfilt: transform to V0");877 859 v[j]=matrix(freemodule(ncols(V[i])),mu,mu*N)*division(V0,V[i])[1]; 878 860 j; … … 881 863 } 882 864 883 dbprint(printlevelvoice+2,"// gaussman::vfilt:compute graded parts");865 dbprint(printlevelvoice+2,"// compute graded parts"); 884 866 for(k=1;k<size(v);k++) 885 867 { … … 960 942 } 961 943 962 dbprint(printlevelvoice+2, 963 "//gaussman::endfilt: compute multiplication in Jacobian algebra"); 944 dbprint(printlevelvoice+2,"// compute multiplication in Jacobian algebra"); 964 945 list M; 965 946 matrix U=freemodule(ncols(m)); … … 979 960 while(b0<number(a[ncols(a)]a[1])) 980 961 { 981 dbprint(printlevelvoice+2, 982 "//gaussman::endfilt: find next possible index"); 962 dbprint(printlevelvoice+2,"// find next possible index"); 983 963 b1=number(a[ncols(a)]a[1]); 984 964 for(j=ncols(a);j>=1;j) … … 1008 988 } 1009 989 1010 dbprint(printlevelvoice+2, 1011 "//gaussman::endfilt: collect conditions for V1["+string(b0)+"]"); 990 dbprint(printlevelvoice+2,"// collect conditions for V1["+string(b0)+"]"); 1012 991 j=ncols(a); 1013 992 j0=mu; … … 1038 1017 } 1039 1018 1040 dbprint(printlevelvoice+2, 1041 "//gaussman::endfilt: compose condition matrix"); 1019 dbprint(printlevelvoice+2,"// compose condition matrix"); 1042 1020 L=transpose(module(l[1])); 1043 1021 for(k=2;k<=ncols(m);k++) … … 1046 1024 } 1047 1025 1048 dbprint(printlevelvoice+2, 1049 "//gaussman::endfilt: compute kernel of condition matrix"); 1026 dbprint(printlevelvoice+2,"// compute kernel of condition matrix"); 1050 1027 v0=v0+list(syz(L)); 1051 1028 a0=a0,b0; 1052 1029 } 1053 1030 1054 dbprint(printlevelvoice+2,"// gaussman::endfilt:compute graded parts");1031 dbprint(printlevelvoice+2,"// compute graded parts"); 1055 1032 option(redSB); 1056 1033 for(i=1;i<size(v0);i++) … … 1060 1037 } 1061 1038 1062 dbprint(printlevelvoice+2, 1063 "//gaussman::endfilt: remove trivial graded parts"); 1039 dbprint(printlevelvoice+2,"// remove trivial graded parts"); 1064 1040 i=1; 1065 1041 while(size(v0[i])==0) … … 1130 1106 list Spp: spectral pairs of f 1131 1107 ideal Spp[1]: spectral numbers in increasing order 1132 intvec Spp[2]: corresponding nilpotencyindices in increasing order1108 intvec Spp[2]: corresponding weight filtration indices in increasing order 1133 1109 intvec Spp[3]: corresponding multiplicities of spectral pairs 1134 1110 default: opt=1 … … 1153 1129 def G=gmsring(f,"s"); 1154 1130 setring(G); 1131 1155 1132 int mu,modm=ncols(gmsbasis),maxorddif(gmsbasis); 1156 1133 ideal r=gmspoly*gmsbasis; … … 1164 1141 { 1165 1142 k++; 1143 dbprint(printlevelvoice+2,"// k="+string(k)); 1144 dbprint(printlevelvoice+2,"// compute matrix A of t"); 1166 1145 l=gmscoeffs(r,k,mu+n); 1167 1146 C,r=l[1..2]; 1168 1147 A0=A0+C; 1148 1149 dbprint(printlevelvoice+2,"// compute saturation of H''"); 1169 1150 H0=H; 1170 1151 dH=jet(module(A0*dH+s^2*diff(matrix(dH),s)),k); … … 1173 1154 A0=A0k*s; 1174 1155 1156 dbprint(printlevelvoice+2,"// compute basis of saturation of H''"); 1175 1157 H=std(H0); 1176 1158 int modH=maxorddif(H)/deg(s); 1159 dbprint(printlevelvoice+2,"// transform H'' to saturation of H''"); 1177 1160 l=division(H,freemodule(mu)*s^k); 1178 1161 H0=l[1]; 1162 1163 dbprint(printlevelvoice+2,"// k="+string(modH+1)); 1164 dbprint(printlevelvoice+2,"// compute matrix A of t"); 1179 1165 l=gmscoeffs(r,modH+1,modH+1+n); 1180 1166 C,r=l[1..2]; 1181 1167 A0=A0+C; 1168 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 1182 1169 l=division(H*s,A0*H+s^2*diff(matrix(H),s)); 1183 1170 matrix A=jet(l[1],l[2],0); … … 1187 1174 if(!opt) 1188 1175 { 1189 U=jordanbasis(A); 1176 dbprint(printlevelvoice+2,"// transform to Jordan basis"); 1177 U=jordanbasis(A)[1]; 1190 1178 V=lift(U,freemodule(mu)); 1191 1179 A=V*A*U; 1180 dbprint(printlevelvoice+2,"// compute normal form of H''"); 1192 1181 H0=std(V*H0); 1182 1183 dbprint(printlevelvoice+2,"// compute spectral numbers"); 1193 1184 ideal a; 1194 1185 for(i=1;i<=mu;i++) … … 1197 1188 a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)1; 1198 1189 } 1190 1199 1191 setring(R); 1200 1192 return(spgen(imap(G,a))); 1201 1193 } 1202 1194 1195 dbprint(printlevelvoice+2, 1196 "// compute eigenvalues e with multiplicities m of A"); 1203 1197 l=eigenval(A); 1204 def e,b=l[1..2]; 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"); 1205 1202 int mide=maxintdif(e); 1203 dbprint(printlevelvoice+2,"// mide="+string(mide)); 1204 1206 1205 if(mide>0) 1207 1206 { 1207 dbprint(printlevelvoice+2,"// k="+string(modH+1+mide)); 1208 dbprint(printlevelvoice+2,"// compute matrix A of t"); 1208 1209 l=gmscoeffs(r,modH+1+mide,modH+1+mide); 1209 1210 C,r=l[1..2]; 1210 1211 A0=A0+C; 1212 dbprint(printlevelvoice+2,"// transform A to saturation of H''"); 1211 1213 l=division(H*s,A0*H+s^2*diff(matrix(H),s)); 1212 1214 A=jet(l[1],l[2],mide); … … 1231 1233 for(j,k=ncols(e),mu;j>=1;j) 1232 1234 { 1233 for(i= b[j];i>=1;i)1235 for(i=m[j];i>=1;i) 1234 1236 { 1235 1237 ide[k]=ide[j]; … … 1240 1242 while(mide>0) 1241 1243 { 1244 dbprint(printlevelvoice+2,"// transform to separate eigenvalues"); 1242 1245 A0=jet(A,0); 1243 1246 U=0; … … 1250 1253 H0=V*H0; 1251 1254 1255 dbprint(printlevelvoice+2,"// transform to reduce mide by 1"); 1252 1256 for(i=mu;i>=1;i) 1253 1257 { … … 1278 1282 } 1279 1283 H0=transpose(H0); 1284 1280 1285 mide; 1286 dbprint(printlevelvoice+2,"// mide="+string(mide)); 1281 1287 } 1282 1288 … … 1284 1290 } 1285 1291 1286 U=jordanbasis(A); 1292 dbprint(printlevelvoice+2,"// transform to Jordan basis"); 1293 intvec w0; 1294 l=jordanbasis(A); 1295 U,w0=l[1..2]; 1287 1296 V=lift(U,freemodule(mu)); 1288 1297 A0=V*A*U; 1289 1298 1290 intvec w; 1291 w[mu]=1; 1292 j=1; 1293 for(i=mu1;i>=1;i) 1294 { 1295 j++; 1296 if(A0[i,i+1]==0) 1297 { 1298 j=1; 1299 } 1300 w[i]=j; 1301 } 1302 1299 dbprint(printlevelvoice+2,"// compute weight filtration basis"); 1303 1300 vector u; 1304 for(i=1;i<ncols(A);i++) 1301 i=1; 1302 while(i<ncols(A)) 1305 1303 { 1306 1304 j=i+1; 1307 1305 while(j<ncols(A)&&A[i,i]==A[j,j]) 1308 1306 { 1309 if(w [i]<w[j])1310 { 1311 k=w [i];1312 w [i]=w[j];1313 w [i]=k;1307 if(w0[i]>w0[j]) 1308 { 1309 k=w0[i]; 1310 w0[i]=w0[j]; 1311 w0[i]=k; 1314 1312 u=U[i]; 1315 1313 U[i]=U[j]; … … 1318 1316 j++; 1319 1317 } 1320 } 1321 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"); 1322 1331 V=lift(U,freemodule(mu)); 1323 1332 A=V*A*U; 1333 dbprint(printlevelvoice+2,"// compute normal form of H''"); 1324 1334 H0=std(V*H0); 1335 1336 dbprint(printlevelvoice+2,"// compute spectral pairs"); 1325 1337 ideal a; 1338 intvec w; 1326 1339 for(i=1;i<=mu;i++) 1327 1340 { 1328 1341 j=leadexp(H0[i])[nvars(basering)+1]; 1329 1342 a[i]=A[j,j]+deg(lead(H0[i]))/deg(s)1; 1343 w[i]=w0[j]+n; 1330 1344 } 1331 1345 setring(R); … … 1353 1367 ideal a0=jet(a,0); 1354 1368 int i,j; 1355 number n;1369 poly p; 1356 1370 for(i=1;i<=ncols(a0);i++) 1357 1371 { … … 1360 1374 if(number(a0[i])>number(a0[j])) 1361 1375 { 1362 n=a0[i];1376 p=a0[i]; 1363 1377 a0[i]=a0[j]; 1364 a0[j]= n;1378 a0[j]=p; 1365 1379 } 1366 1380 } … … 1407 1421 intvec w0=w; 1408 1422 int i,j,k; 1409 number n;1423 poly p; 1410 1424 for(i=1;i<=ncols(a0);i++) 1411 1425 { … … 1414 1428 if(number(a0[i])>number(a0[j])a0[i]==a0[j]&&w0[i]>w0[j]) 1415 1429 { 1416 n=a0[i];1430 p=a0[i]; 1417 1431 a0[i]=a0[j]; 1418 a0[j]= n;1432 a0[j]=p; 1419 1433 k=w0[i]; 1420 1434 w0[i]=w0[j];
Note: See TracChangeset
for help on using the changeset viewer.