Changeset 08fa62 in git
- Timestamp:
- Oct 14, 2013, 6:15:36 PM (9 years ago)
- Branches:
- (u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
- Children:
- 76d26c42ef4754edaa37b1777b4d2831f6f7e349
- Parents:
- ac665c6c03e4e4354b586d1fd770d8d5669253cf
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/derham.lib
rac665c r08fa62 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version=""; 3 category="Noncommutative"; 4 info=" 5 LIBRARY: derham.lib Computation of deRham cohomology 6 AUTHORS: Cornelia Rottner, rottner@mathematik.uni-kl.de 7 OVERVIEW: 8 PROCEDURES: 9 10 "; 11 12 13 LIB "nctools.lib"; 14 LIB "matrix.lib"; 15 LIB "qhmoduli.lib"; 16 LIB "general.lib"; 17 LIB "dmod.lib"; 18 LIB "bfun.lib"; 19 LIB "dmodapp.lib"; 20 LIB "poly.lib"; 21 22 ///////////////////////////////////////////////////////////////////////////// 23 24 static proc divdr(matrix m, matrix n) 25 { 26 m=transpose(m); 27 n=transpose(n); 28 matrix con=concat(m,n); 29 matrix s=syz(con); 30 s=submat(s,1..ncols(m),1..ncols(s)); 31 s=transpose(compress(s)); 32 return(s); 33 } 34 35 ///////////////////////////////////////////////////////////////////////////// 36 37 38 static proc matrixlift(matrix M, matrix N) 39 { 40 // option(noreturnSB); 41 matrix l=transpose(lift(transpose(M),transpose(N))); 42 return(l); 43 } 44 45 /////////////////////////////////////////////////////////////////////////////// 46 47 proc shortexactpieces(list #) 48 { 49 matrix Bnew= divdr(#[2],#[3]); 50 matrix Bold=Bnew; 51 matrix Z=divdr(Bnew,#[1]); 52 list bzh; list zcb; 53 bzh=list(list(),list(),Z,unitmat(ncols(Z)),Z); 54 zcb=(Z, Bnew, #[1], unitmat(ncols(#[1])), Bnew); 55 list sep; 56 sep[1]=(list(bzh,zcb)); 57 int i; 58 list out; 59 for (i=3; i<=size(#)-2; i=i+2) 60 { 61 out=bzhzcb(Bold, #[i-1] , #[i], #[i+1], #[i+2]); 62 sep[size(sep)+1]=out[1]; 63 Bold=out[2]; 64 } 65 bzh=(divdr(#[size(#)-2], #[size(#)-1]),#[size(#)-2], #[size(#)-1],unitmat(ncols(#[size(#)-1])),transpose(concat(transpose(#[size(#)-2]),transpose(#[size(#)-1])))); 66 zcb=(#[size(#)-1], unitmat(ncols(#[size(#)-1])), #[size(#)-1],list(),list()); 67 sep[size(sep)+1]=list(bzh,zcb); 68 return(sep); 69 } 70 71 //////////////////////////////////////////////////////////////////////////////////////// 72 73 static proc bzhzcb (matrix Bold, matrix f0, matrix C1, matrix f1,matrix C2) 74 { 75 matrix Bnew=divdr(f1,C2); 76 matrix Z= divdr(Bnew,C1); 77 matrix lift1= matrixlift(Bnew,f0); 78 list bzh=(Bold, lift1, Z, unitmat(ncols(Z)), transpose(concat(transpose(lift1),transpose(Z)))); 79 list zcb=(Z, Bnew, C1, unitmat(ncols(C1)),Bnew); 80 list out=(list(bzh, zcb), Bnew); 81 return(out); 82 } 83 84 ////////////////////////////////////////////////////////////////////////////////////// 85 86 proc VdstrictGB (matrix M, int d ,list #); 87 "USAGE:VdstrictGB(M,d[,v]); M a matrix, d an integer, v an optional intvec(shift vector) 88 RETURN:matrix M; the rows of M foem a Vd-strict Groebner basis for imM 89 ASSUME:1<=d<=nvars(basering)/2; size(v)=ncols(M) 90 " 91 { 92 if (M==matrix(0,nrows(M),ncols(M))) 93 { 94 return (matrix(0,1,ncols(M))); 95 } 96 def W =basering; 97 int ncM=ncols(M); 98 list Data=ringlist(W); 99 Data[2]=list("nhv")+Data[2]; 100 Data[3][3]=Data[3][1]; 101 Data[3][1]=Data[3][2]; 102 Data[3][2]=list("dp",intvec(1)); 103 matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2])); 104 Data[5]=re; 105 int k; int l; 106 Data[6]=transpose(concat(matrix(0,1,1),transpose(concat(matrix(0,1,1),Data[6])))); 107 def Whom=ring(Data); 108 setring Whom; 109 matrix Mnew=imap(W,M); 110 intvec v; 111 if (size(#)!=0) 112 { 113 v=#[1]; 114 } 115 if (size(v) < ncM) 116 { 117 v=v,0:(ncM-size(v)); 118 } 119 Mnew=homogenize(Mnew, d, v); 120 Mnew=transpose(Mnew); 121 Mnew=std(Mnew); 122 Mnew=subst(Mnew,nhv,1); 123 Mnew=transpose(Mnew); 124 setring W; 125 M=imap(Whom,Mnew); 126 return(M); 127 } 128 129 //////////////////////////////////////////////////////////////////////////////////// 130 131 static proc Vdnormalform(matrix F, matrix M, int d, intvec v) 132 { 133 def W =basering; 134 int c=ncols(M); 135 F=submat(F,intvec(1..nrows(F)),intvec(1..c)); 136 list Data=ringlist(W); 137 Data[2]=list("nhv")+Data[2]; 138 Data[3][3]=Data[3][1]; 139 Data[3][1]=Data[3][2]; 140 Data[3][2]=list("dp",intvec(1)); 141 matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2])); 142 Data[5]=re; 143 int k; 144 int l; 145 matrix rep[size(Data[2])][size(Data[2])]; 146 for (l=size(Data[2])-1;l>=1; l--) 147 { 148 for (k=l-1; k>=1;k--) 149 { 150 rep[k+1,l+1]=Data[6][k,l]; 151 } 152 } 153 Data[6]=rep; 154 def Whom=ring(Data); 155 setring Whom; 156 matrix Mnew=imap(W,M); 157 Mnew=(homogenize(Mnew, d, v));//doppelte Berechung unnötig->muss noch geändert werden!!! 158 matrix Fnew=imap(W,F); 159 matrix Fb; 160 for (l=1; l<=nrows(Fnew); l++) 161 { 162 Fb=homogenize(submat(Fnew,l,intvec(1..ncols(Fnew))),d,v); 163 Fb=transpose(reduce(transpose(Fb),std(transpose(Mnew))));// doppelte Berechnung unnötig, unterdrückt aber Fehler meldung 164 for (k=1; k<=ncols(Fnew);k++) 165 { 166 Fnew[l,k]=Fb[1,k]; 167 } 168 } 169 Fnew=subst(Fnew,nhv,1); 170 setring W; 171 F=imap(Whom,Fnew); 172 return(F); 173 } 174 175 176 /////////////////////////////////////////////////////////////////////////////// 177 178 static proc homogenize (matrix M, int d, intvec v) 179 { 180 int l; poly f; int s; int i; intvec vnm;int kmin; list findmin; 181 int n=(nvars(basering)-1) div 2; 182 list rempoly; 183 list remk; 184 list rem1; 185 list rem2; 186 for (int k=1; k<=nrows(M); k++) //man könnte auch paralell immer weiter homogenisieren, d.h. immer ein enues Minimum finden und das dann machen 187 { 188 for (l=1; l<=ncols (M); l++) 189 { 190 f=M[k,l]; 191 s=size(f); 192 for (i=1; i<=s; i++) 193 { 194 vnm=leadexp(f); 195 vnm=vnm[n+2..n+d+1]-vnm[2..d+1]; 196 kmin=sum(vnm)+v[l]; 197 rem1[size(rem1)+1]=lead(f); 198 rem2[size(rem2)+1]=kmin; 199 findmin=insert(findmin,kmin); 200 f=f-lead(f); 201 } 202 rempoly[l]=rem1; 203 remk[l]=rem2; 204 rem1=list(); 205 rem2=list(); 206 } 207 if (size(findmin)!=0) 208 { 209 kmin=Min(findmin); 210 } 211 for (l=1; l<=ncols(M); l++) 212 { 213 if (M[k,l]!=0) 214 { 215 M[k,l]=0; 216 for (i=1; i<=size(rempoly[l]);i++) 217 { 218 M[k,l]=M[k,l]+nhv^(remk[l][i]-kmin)*rempoly[l][i]; 219 } 220 } 221 } 222 rempoly=list(); 223 remk=list(); 224 findmin=list(); 225 } 226 return(M); 227 } 228 229 230 ////////////////////////////////////////////////////////////////////////////////////// 231 232 static proc soldr (matrix M, matrix N) 233 { 234 int n=nrows(M); 235 int q=ncols(M); 236 matrix S=concat(transpose(M),transpose(N)); 237 def W=basering; 238 list Data=ringlist(W); 239 list Save=Data[3]; 240 Data[3]=list(list("c",0),list("dp",intvec(1..nvars(W)))); 241 def Wmod=ring(Data); 242 setring Wmod; 243 matrix Smod=imap(W,S); 244 matrix E[q][1]; 245 matrix Smod2; 246 matrix Smodnew; 247 option(returnSB); 248 int i; int j; 249 for (i=1;i<=q;i++) 250 { 251 E[i,1]=1; 252 Smod2=concat(E,Smod); 253 print (Smod2); 254 Smod2=syz(Smod2); 255 E[i,1]=0; 256 for (j=1;j<=ncols(Smod2);j++) 257 { 258 if (Smod2[1,j]==1) 259 { 260 Smodnew=concat(Smodnew,(-1)*(submat(Smod2,intvec(2..n+1),j))); 261 break; 262 } 263 } 264 } 265 Smodnew=transpose(submat(Smodnew,intvec(1..n),intvec(2..q+1))); 266 setring W; 267 matrix Snew=imap(Wmod,Smodnew); 268 return (Snew); 269 } 270 271 272 ///////////////////////////////////////////////////////////////////////////// 273 274 proc toVdstrictsequence (list C,int n, intvec v) 275 276 { 277 matrix J_C=VdstrictGB(C[5],n,list(v)); 278 matrix J_A=C[1]; 279 matrix f_CB=C[4]; 280 matrix f_ACB=transpose(concat(transpose(C[2]),transpose(f_CB))); 281 matrix J_AC=divdr(f_ACB,C[3]); 282 matrix P=matrixlift(J_AC * prodr(ncols(C[1]),ncols(C[5])) ,J_C); 283 list storePi; 284 matrix Pi[1][ncols(J_AC)]; 285 int i;int j; 286 for (i=1; i<=nrows(J_C); i++) 287 { 288 for (j=1; j<=nrows(J_AC);j++) 289 { 290 Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC))); 291 } 292 storePi[i]=Pi; 293 Pi=0; 294 } 295 intvec m_a; 296 list findMin; 297 int comMin; 298 for (i=1; i<=ncols(J_A); i++) 299 { 300 for (j=1; j<=size(storePi);j++) 301 { 302 if (storePi[j][1,i]!=0) 303 { 304 comMin=Vddeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v)-Vddeg(storePi[j][1,i],n,intvec(0)); 305 findMin[size(findMin)+1]=comMin; 306 } 307 } 308 if (size(findMin)!=0) 309 { 310 m_a[i]=Min(findMin); 311 findMin=list(); 312 } 313 else 314 { 315 m_a[i]=0; 316 } 317 } 318 matrix zero[ncols(J_A)][ncols(J_C)]; 319 320 matrix g_AB=concat(unitmat(ncols(J_A)),zero); 321 matrix g_BC= transpose(concat(transpose(zero),transpose(unitmat(ncols(J_C))))); 322 intvec m_b=m_a,v; 323 324 J_A=VdstrictGB(J_A,n,m_a); 325 J_AC=transpose(storePi[1]); 326 for (i=2; i<= size(storePi); i++) 327 { 328 J_AC=concat(J_AC, transpose(storePi[i])); 329 } 330 J_AC=transpose(concat(transpose(matrix(J_A,nrows(J_A),nrows(J_AC))),J_AC)); 331 332 list Vdstrict=(list(J_A),list(g_AB),list(J_AC),list(g_BC),list(J_C),list(m_a),list(m_b),list(v)); 333 return (Vdstrict); 334 } 335 336 337 ///////////////////////////////////////////////////////////////////////// 338 339 static proc prodr (int k, int l) 340 { 341 if (k==0) 342 { 343 matrix P=unitmat(l); 344 return (P); 345 } 346 matrix O[l][k]; 347 matrix P=transpose(concat(O,unitmat(l))); 348 return (P); 349 } 350 351 ///////////////////////////////////////////////////////////////////////// 352 353 proc Vddeg(matrix M, int d, intvec v, list #)//Aternative: in WHom Leadmonom ausrechnen!! 354 "USAGE: Vddeg(M,d,v); M 1xr-matrix, d int, v intvec of size r 355 RETURN:int; the Vd-degree of M 356 " 357 { 358 int i;int j; 359 int n=nvars(basering) div 2; 360 intvec e; 361 int etoint; 362 list findmax; 363 int c=ncols(M); 364 poly l; 365 list positionpoly; 366 list positionVd; 367 for (i=1; i<=c; i++) 368 { 369 positionpoly[i]=list(); 370 positionVd[i]=list(); 371 while (M[1,i]!=0) 372 { 373 l=lead(M[1,i]); 374 positionpoly[i][size(positionpoly[i])+1]=l; 375 e=leadexp(l); 376 e=-e[1..d]+e[n+1..n+d]; 377 e=sum(e)+v[i]; 378 etoint=e[1]; 379 positionVd[i][size(positionVd[i])+1]=etoint; 380 findmax[size(findmax)+1]=etoint; 381 M[1,i]=M[1,i]-l; 382 } 383 } 384 if (size(findmax)!=0) 385 { 386 int maxVd=Max(findmax); 387 if (size(#)==0) 388 { 389 return (maxVd); 390 } 391 } 392 else // M is 0-modul 393 { 394 return(int(0)); 395 } 396 l=0; 397 for (i=c; i>=1; i--) 398 { 399 for (j=1; j<=size(positionVd[i]); j++) 400 { 401 if (positionVd[i][j]==maxVd) 402 { 403 l=l+positionpoly[i][j]; 404 } 405 } 406 if (l!=0) 407 { 408 return (list(l,i)); 409 } 410 } 411 412 } 413 414 415 /////////////////////////////////////////////////////////////////////////////// 416 417 proc toVdstrictsequences (list L,int d, intvec v) 418 "USAGE: toVdstrictsequences(L,d,v); L list, d int, v intvec, L contains two lists of short exact sequences(D,f_DA,A,f_AF,F) and (A,f_AB,B,f_BC,C), v is a shift vector on the range of C 419 RETURN: list of two lists; each lists contains Vd-strict exact sequences with corresponding shift vectors 420 " 421 { 422 matrix J_F=L[1][5]; 423 matrix J_D=L[1][1]; 424 matrix f_FA=L[1][4]; 425 matrix f_DFA=transpose(concat(transpose(L[1][2]),transpose(f_FA))); 426 matrix J_DF=divdr(f_DFA,L[1][3]); 427 matrix J_C=L[2][5]; 428 matrix f_CB=L[2][4]; 429 matrix f_DFCB=transpose(concat(transpose(f_DFA*L[2][2]),transpose(f_CB))); 430 matrix J_DFC=divdr(f_DFCB,L[2][3]); 431 matrix P=matrixlift(J_DFC*prodr(ncols(J_DF),ncols(L[2][5])),J_C); 432 list storePi; 433 matrix Pi[1][ncols(J_DFC)]; 434 int i; int j; 435 for (i=1; i<=nrows(J_C); i++) 436 { 437 438 for (j=1; j<=nrows(J_DFC);j++) 439 { 440 Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC))); 441 } 442 storePi[i]=Pi; 443 Pi=0; 444 } 445 intvec m_a; 446 list findMin; 447 list noMin; 448 int comMin; 449 for (i=1; i<=ncols(J_DF); i++) 450 { 451 for (j=1; j<=size(storePi);j++) 452 { 453 if (storePi[j][1,i]!=0) 454 { 455 comMin=Vddeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v)-Vddeg(storePi[j][1,i],d,intvec(0)); 456 findMin[size(findMin)+1]=comMin; 457 } 458 } 459 if (size(findMin)!=0) 460 { 461 m_a[i]=Min(findMin); 462 findMin=list(); 463 noMin[i]=0; 464 } 465 else 466 { 467 noMin[i]=1; 468 } 469 } 470 if (size(m_a) < ncols(J_DF)) 471 { 472 m_a[ncols(J_DF)]=0; 473 } 474 intvec m_f=m_a[ncols(J_D)+1..size(m_a)]; 475 J_F=VdstrictGB(J_F,d,m_f); 476 P=matrixlift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F);// selbe Prinzip wie oben--> evtl auslagern 477 list storePinew; 478 matrix Pidf[1][ncols(J_DF)]; 479 for (i=1; i<=nrows(J_F); i++) 480 { 481 for (j=1; j<=nrows(J_DF);j++) 482 { 483 Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF))); 484 } 485 storePinew[i]=Pidf; 486 Pidf=0; 487 } 488 intvec m_d; 489 for (i=1; i<=ncols(J_D); i++) 490 { 491 for (j=1; j<=size(storePinew);j++) 492 { 493 if (storePinew[j][1,i]!=0) 494 { 495 comMin=Vddeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f)-Vddeg(storePinew[j][1,i],d,intvec(0)); 496 findMin[size(findMin)+1]=comMin; 497 } 498 } 499 if (size(findMin)!=0) 500 { 501 if (noMin[i]==0) 502 { 503 m_d[i]=Min(insert(findMin,m_a[i])); 504 m_a[i]=m_d[i]; 505 } 506 else 507 { 508 m_d[i]=Min(findMin); 509 m_a[i]=m_d[i]; 510 } 511 } 512 else 513 { 514 m_d[i]=m_a[i]; 515 } 516 findMin=list(); 517 } 518 J_D=VdstrictGB(J_D,d,m_d); 519 J_DF=transpose(storePinew[1]); 520 for (i=2; i<=nrows(J_F); i++) 521 { 522 J_DF=concat(J_DF,transpose(storePinew[i])); 523 } 524 J_DF=transpose(concat(transpose(matrix(J_D,nrows(J_D),nrows(J_DF))),J_DF)); 525 J_DFC=transpose(storePi[1]); 526 for (i=2; i<=nrows(J_C); i++) 527 { 528 J_DFC=concat(J_DFC,transpose(storePi[i])); 529 } 530 J_DFC=transpose(concat(transpose(matrix(J_DF,nrows(J_DF),nrows(J_DFC))),J_DFC)); 531 intvec m_b=m_a,v; 532 matrix zero[ncols(J_D)][ncols(J_F)]; 533 matrix g_DA=concat(unitmat(ncols(J_D)),zero); 534 matrix g_AF=transpose(concat(transpose(zero),unitmat(ncols(J_F)))); 535 matrix zero1[ncols(J_DF)][ncols(J_C)]; 536 matrix g_AB=concat(unitmat(ncols(J_DF)),zero1); 537 matrix g_BC=transpose(concat(transpose(zero1),unitmat(ncols(J_C)))); 538 list out=(list(list(J_D),list(g_DA),list(J_DF),list(g_AF),list(J_F),list(m_d),list(m_a),list(m_f)),list(list(J_DF),list(g_AB),list(J_DFC),list(g_BC),list(J_C),list(m_a),list(m_b),list(v))); 539 return(out), 540 } 541 542 /////////////////////////////////////////////////////////////////////////////////////////// 543 544 proc shortexactpiecestoVdstrict(list C, int d,list #) 545 546 { 547 548 int s =size(C); 549 if (size(#)==0) 550 { 551 intvec v=0:ncols(C[s][1][5]); 552 } 553 else 554 { 555 intvec v=#[1]; 556 } 557 list out; 558 out[s]=list(toVdstrictsequence(C[s][1],d,v)); 559 out[s][2]=list(list(out[s][1][3][1]),list(unitmat(ncols(out[s][1][3][1]))),list(out[s][1][3][1]),list(list()),list(list())); 560 out[s][2][6]=list(out[s][1][7][1]); 561 out[s][2][7]=list(out[s][2][6][1]); 562 out[s][2][8]=list(list()); 563 int i; 564 for (i=s-1; i>=2; i--) 565 { 566 C[i][2][5]=out[i+1][1][1][1]; 567 out[i]=toVdstrictsequences(C[i],d,out[i+1][1][6][1]); 568 } 569 out[1]=list(list()); 570 out[1][2]=toVdstrictsequence(C[1][2],d,out[2][1][6][1]); 571 out[1][1][3]=list(out[1][2][1][1]); 572 out[1][1][5]=list(out[1][2][1][1]); 573 out[1][1][4]=list(unitmat(ncols(out[1][1][3][1]))); 574 out[1][1][7]=list(out[1][2][6][1]); 575 out[1][1][8]=list(out[1][2][6][1]); 576 out[1][1][1]=list(list()); 577 out[1][1][2]=list(list()); 578 out[1][1][6]=list(list()); 579 list Hi; 580 for (i=1; i<=size(out); i++) 581 { 582 Hi[i]=list(out[i][1][5][1],out[i][1][8][1]); 583 } 584 list outall; 585 outall[1]=out; 586 print (out); 587 outall[2]=Hi; 588 return(outall); 589 590 591 } 592 593 /////////////////////////////////////////////////////////////////////////////////////////// 594 595 proc toVdstrict2x3complex(list L, int d, list #) 596 { 597 matrix rem; int i; int j; 598 list J_A=list(list()); 599 list J_B=list(list()); 600 list J_C=list(list()); 601 list g_AB=list(list()); 602 list g_BC=list(list()); 603 list n_a=list(list()); 604 list n_b=list(list()); 605 list n_c=list(list()); 606 intvec n_b1; 607 if (size(L[5])!=0) 608 { 609 intvec n_c1; 610 for (i=1; i<=nrows(L[5]); i++) 611 { 612 rem=submat(L[5],i,intvec(1..ncols(L[5]))); 613 n_c1[i]=Vddeg(rem,d, L[8]); 614 } 615 n_c[1]=n_c1; 616 J_C[1]=transpose(syz(transpose(L[5]))); 617 if (J_C[1]!=matrix(0,nrows(J_C[1]),ncols(J_C[1]))) 618 { 619 J_C[1]=VdstrictGB(J_C[1],d,n_c1); 620 if (size(#[2])!=0) 621 { 622 n_a[1]=#[2]; 623 n_b1=n_a[1],n_c[1]; 624 n_b[1]=n_b1; 625 matrix zero[nrows(L[1])][nrows(L[5])]; 626 g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5]))); 627 if (size(#[1])!=0) 628 { 629 J_A=#[1]; 630 J_B=transpose(matrix(syz(transpose(L[3])))); 631 matrix P=matrixlift(J_B[1] * prodr(nrows(L[1]),nrows(L[5])) ,J_C[1]); 632 633 matrix Pi[1][ncols(J_B[1])]; 634 matrix Picombined; 635 for (i=1; i<=nrows(J_C[1]); i++) 636 { 637 for (j=1; j<=nrows(J_B[1]);j++) 638 { 639 Pi=Pi+P[i,j]*submat(J_B[1],j,intvec(1..ncols(J_B[1]))); 640 641 } 642 if (i==1) 643 { 644 Picombined=transpose(Pi); 645 } 646 else 647 { 648 Picombined=concat(Picombined,transpose(Pi)); 649 } 650 Pi=0; 651 } 652 Picombined=transpose(Picombined); 653 Picombined=concat(Vdnormalform(Picombined,J_A[1],d,n_a[1]),submat(Picombined,intvec(1..nrows(Picombined)),intvec((ncols(J_A[1])+1)..ncols(Picombined)))); 654 J_B[1]=transpose(concat(transpose(matrix(J_A[1],nrows(J_A[1]),ncols(J_B[1]))),transpose(Picombined))); 655 g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5])))); 656 } 657 else 658 { 659 J_B[1]=concat(matrix(0,nrows(J_C[1]),nrows(L[3])-nrows(L[5])),J_C[1]); 660 g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5])))); 661 } 662 } 663 else 664 { 665 n_b=n_c[1]; 666 J_B[1]=J_C[1]; 667 g_BC=unitmat(ncols(J_C[1])); 668 } 669 } 670 else 671 { 672 J_C=list(list()); 673 if (size(#[2])!=0) 674 { 675 matrix zero[nrows(L[1])][nrows(L[5])]; 676 g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5])))); 677 n_a[1]=#[2]; 678 n_b1=n_a[1],n_c[1]; 679 n_b[1]=n_b1; 680 g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5])));; 681 682 if (size(#[1])!=0) 683 { 684 J_A=#[1]; 685 J_B=concat(J_A[1],matrix(0,nrows(J_A[1]),nrows(L[3])-nrows(L[1]))); 686 } 687 } 688 else 689 { 690 n_b=n_c[1]; 691 g_BC=unitmat(ncols(L[5])); 692 } 693 694 } 695 } 696 else 697 { 698 if (size(#[2])!=0) 699 { 700 n_a[1]=#[2]; 701 n_b=n_a[1]; 702 g_AB=unitmat(size(n_b[1])); 703 if (size(#[1])!=0) 704 { 705 J_A=#[1]; 706 J_B[1]=J_A[1]; 707 } 708 } 709 } 710 list out=(J_A[1],g_AB[1],J_B[1],g_BC[1],J_C[1],n_a[1],n_b[1],n_c[1]); 711 return (out); 712 } 713 714 715 ////////////////////////////////////////////////////////////////////////// 716 717 proc Vdstrictdoublecompexes(list L, int d) 718 { 719 int i; int k; int c; int j; 720 intvec n_b; 721 matrix rem; 722 matrix J_B; 723 list store; 724 int t=size(L)+nvars(basering) div 2-2; 725 for (k=1; k<=(size(L)+nvars(basering) div 2-3); k++)// 726 { 727 L[1][1][1][k+1]=list(); 728 L[1][1][2][k+1]=list(); 729 L[1][1][6][k+1]=list(); 730 if (size(L[1][1][3][k])!=0) 731 { 732 for (i=1; i<=nrows(L[1][1][3][k]); i++) 733 { 734 rem=submat(L[1][1][3][k],i,(1..ncols(L[1][1][3][k]))); 735 n_b[i]=Vddeg(rem,d,L[1][1][7][k]); 736 } 737 J_B=transpose(syz(transpose(L[1][1][3][k]))); 738 L[1][1][7][k+1]=n_b; 739 L[1][1][8][k+1]=n_b; 740 L[1][1][4][k+1]=unitmat(nrows(L[1][1][3][k])); 741 if (J_B!=matrix(0,nrows(J_B),ncols(J_B))) 742 { 743 J_B=VdstrictGB(J_B,d,n_b); 744 L[1][1][3][k+1]=J_B; 745 L[1][1][5][k+1]=J_B; 746 } 747 else 748 { 749 L[1][1][3][k+1]=list(); 750 L[1][1][5][k+1]=list(); 751 } 752 n_b=0; 753 } 754 else 755 { 756 L[1][1][3][k+1]=list(); 757 L[1][1][5][k+1]=list(); 758 L[1][1][7][k+1]=list(); 759 L[1][1][8][k+1]=list(); 760 L[1][1][4][k+1]=list(); 761 } 762 for (i=1; i<size(L); i++) 763 { 764 store=toVdstrict2x3complex(list(L[i][2][1][k],L[i][2][2][k],L[i][2][3][k],L[i][2][4][k],L[i][2][5][k],L[i][2][6][k],L[i][2][7][k],L[i][2][8][k]),d,L[i][1][3][k+1],L[i][1][7][k+1]); 765 for (j=1; j<=8; j++) 766 { 767 L[i][2][j][k+1]=store[j]; 768 } 769 770 store=toVdstrict2x3complex(list(L[i+1][1][1][k],L[i+1][1][2][k],L[i+1][1][3][k],L[i+1][1][4][k],L[i+1][1][5][k],L[i+1][1][6][k],L[i+1][1][7][k],L[i+1][1][8][k]),d,L[i][2][5][k+1],L[i][2][8][k+1]); 771 772 for (j=1; j<=8; j++) 773 { 774 L[i+1][1][j][k+1]=store[j]; 775 } 776 } 777 if (size(L[size(L)][1][7][k+1])!=0) 778 { 779 L[size(L)][2][4][k+1]=list(); 780 L[size(L)][2][5][k+1]=list(); 781 L[size(L)][2][6][k+1]=L[size(L)][1][7][k+1]; 782 L[size(L)][2][7][k+1]=L[size(L)][1][7][k+1]; 783 L[size(L)][2][8][k+1]=list(); 784 L[size(L)][2][2][k+1]=unitmat(size(L[size(L)][1][7][k+1])); 785 786 if (size(L[size(L)][1][3][k+1])!=0) 787 { 788 L[size(L)][2][1][k+1]=L[size(L)][1][3][k+1]; 789 L[size(L)][2][3][k+1]=L[size(L)][1][3][k+1]; 790 } 791 else 792 { 793 L[size(L)][2][1][k+1]=list(); 794 L[size(L)][2][3][k+1]=list(); 795 } 796 } 797 else 798 { 799 for (j=1; j<=8; j++) 800 { 801 L[size(L)][2][j][k+1]=list(); 802 } 803 } 804 } 805 806 807 k=t; 808 intvec n_c; 809 intvec vn_b; 810 list N_b; 811 int n; 812 for (i=1; i<=size(L); i++) 813 { 814 for (n=1; n<=2; n++) 815 { 816 if (i==1 and n==1) 817 { 818 L[i][n][6][k+1]=list(); 819 } 820 else 821 { 822 if (n==1) 823 { 824 L[i][1][6][k+1]=L[i-1][2][8][k+1]; 825 } 826 else 827 { 828 L[i][2][6][k+1]=L[i][1][7][k+1]; 829 } 830 } 831 N_b[1]=L[i][n][6][k+1]; 832 if (size(L[i][n][5][k])!=0) 833 { 834 for (j=1; j<=nrows(L[i][n][5][k]); j++) 835 { 836 rem=submat(L[i][n][5][k],j,(1..ncols(L[i][n][5][k]))); 837 n_c[j]=Vddeg(rem,d,L[i][n][8][k]); 838 } 839 L[i][n][8][k+1]=n_c; 840 } 841 else 842 { 843 L[i][n][8][k+1]=list(); 844 } 845 N_b[2]=L[i][n][8][k+1]; 846 n_c=0; 847 if (size(N_b[1])!=0) 848 { 849 vn_b=N_b[1]; 850 if (size(N_b[2])!=0) 851 { 852 vn_b=vn_b,N_b[2]; 853 } 854 L[i][n][7][k+1]=vn_b; 855 } 856 else 857 { 858 if (size(N_b[2])!=0) 859 { 860 L[i][n][7][k+1]=N_b[2]; 861 } 862 else 863 { 864 L[i][n][7][k+1]=list(); 865 } 866 } 867 868 } 869 } 870 return(L); 871 } 872 873 //////////////////////////////////////////////////////////////////////////// 874 875 proc assemblingdoublecomplexes(list L) 876 { 877 list out; 878 int i; int j;int k;int l; int oldj; int newj; 879 for (i=1; i<=size(L); i++) 880 { 881 out[i]=list(list()); 882 out[i][1][1]=ncols(L[i][2][3][1]); 883 if (size(L[i][2][5][1])!=0) 884 { 885 out[i][1][4]=prodr(ncols(L[i][2][3][1])-ncols(L[i][2][5][1]),ncols(L[i][2][5][1])); 886 } 887 else 888 { 889 out[i][1][4]=matrix(0,ncols(L[i][2][3][1]),1); 890 } 891 892 oldj=newj; 893 for (j=1; j<=size(L[i][2][3]);j++) 894 { 895 out[i][j][2]=L[i][2][7][j]; 896 if (size(L[i][2][3][j])==0) 897 { 898 newj =j; 899 break; 900 } 901 out[i][j+1]=list(); 902 out[i][j+1][1]=nrows(L[i][2][3][j]); 903 out[i][j+1][3]=L[i][2][3][j]; 904 if (size(L[i][2][5][j])!=0) 905 { 906 out[i][j+1][4]=(-1)^j*prodr(nrows(L[i][2][3][j])-nrows(L[i][2][5][j]),nrows(L[i][2][5][j])); 907 } 908 else 909 { 910 out[i][j+1][4]=matrix(0,nrows(L[i][2][3][j]),1); 911 } 912 if(j==size(L[i][2][3])) 913 { 914 out[i][j+1][2]=L[i][2][7][j+1]; 915 newj=j+1; 916 } 917 } 918 if (i>1) 919 { 920 for (k=1; k<=Min(list(oldj,newj)); k++) 921 { 922 out[i-1][k][4]=matrix(out[i-1][k][4],nrows(out[i-1][k][4]),out[i][k][1]); 923 } 924 for (k=newj+1; k<=oldj; k++) 925 { 926 out[i-1][k]=delete(out[i-1][k],4); 927 } 928 } 929 } 930 return (out); 931 } 932 933 ////////////////////////////////////////////////////////////////////////////// 934 935 proc totalcomplex(list L); 936 { 937 list out;intvec rem1;intvec v; list remsize; int emp; 938 int i; int j; int c; int d; matrix M; int k; int l; 939 int n=nvars(basering) div 2; 940 list K; 941 for (i=1; i<=n; i++) 942 { 943 K[i]=list(); 944 } 945 L=K+L; 946 for (i=1; i<=size(L); i++) 947 { 948 emp=0; 949 if (size(L[i])!=0) 950 { 951 out[3*i-2]=L[i][1][1]; 952 v=L[i][1][1]; 953 rem1=L[i][1][2]; 954 emp=1; 955 } 956 else 957 { 958 out[3*i-2]=0; 959 v=0; 960 } 961 962 for (j=i+1; j<=size(L); j++) 963 { 964 if (size(L[j])>=j-i+1) 965 { 966 out[3*i-2]=out[3*i-2]+L[j][j-i+1][1]; 967 if (emp==0) 968 { 969 rem1=L[j][j-i+1][2]; 970 emp=1; 971 } 972 else 973 { 974 rem1=rem1,L[j][j-i+1][2]; 975 } 976 v[size(v)+1]=L[j][j-i+1][1]; 977 } 978 else 979 { 980 v[size(v)+1]=0; 981 } 982 } 983 out[3*i-1]=rem1; 984 v[size(v)+1]=0; 985 remsize[i]=v; 986 } 987 int o1; 988 int o2; 989 for (i=1; i<=size(L)-1; i++) 990 { 991 o1=1; 992 o2=1; 993 if (size(out[3*i-2])!=0) 994 { 995 o1=out[3*i-2]; 996 } 997 if (size(out[3*i+1])!=0) 998 { 999 o2=out[3*i+1]; 1000 } 1001 M=matrix(0,o1,o2); 1002 if (size(L[i])!=0) 1003 { 1004 if (size(L[i][1][4])!=0) 1005 { 1006 M=matrix(L[i][1][4],o1,o2); 1007 } 1008 } 1009 c=remsize[i][1]; 1010 // d=remsize[i+1][1]; 1011 for (j=i+1; j<=size(L); j++) 1012 { 1013 if (remsize[i][j-i+1]!=0) 1014 { 1015 for (k=c+1; k<=c+remsize[i][j-i+1]; k++) 1016 { 1017 for (l=d+1; l<=d+remsize[i+1][j-i];l++) 1018 { 1019 M[k,l]=L[j][j-i+1][3][(k-c),(l-d)]; 1020 } 1021 } 1022 d=d+remsize[i+1][j-i]; 1023 if (remsize[i+1][j-i+1]!=0) 1024 { 1025 for (k=c+1; k<=c+remsize[i][j-i+1]; k++) 1026 { 1027 for (l=d+1; l<=d+remsize[i+1][j-i+1];l++) 1028 { 1029 M[k,l]=L[j][j-i+1][4][k-c,l-d]; 1030 } 1031 } 1032 c=c+remsize[i][j-i+1]; 1033 } 1034 } 1035 else 1036 { 1037 d=d+remsize[i+1][j-i]; 1038 } 1039 } 1040 out[3*i]=M; 1041 d=0; c=0; 1042 } 1043 out[3*size(L)]=matrix(0,out[3*size(L)-2],1); 1044 return (out); 1045 } 1046 1047 ///////////////////////////////////////////////////////////////////////////////////// 1048 1049 proc toVdstrictfreecomplex(list L,list #) 1050 { 1051 def B=basering; 1052 int n=nvars(B) div 2+2; 1053 int d=nvars(B) div 2; 1054 intvec v; 1055 list out;list outall; 1056 int i;int j; 1057 matrix mem; 1058 int k; 1059 if (size(#)!=0) 1060 { 1061 for (i=1; i<=size(#); i++) 1062 { 1063 if (typeof(#[i])==intvec) 1064 { 1065 v=#[i]; 1066 } 1067 if (typeof(#[i])==int) 1068 { 1069 d=#[i]; 1070 } 1071 } 1072 } 1073 if (size(L)==2) 1074 { 1075 v=(0:ncols(L[1])); 1076 out[3*n-1]=v; 1077 out[3*n-2]=ncols(L[1]); 1078 out[3*n]=L[2]; 1079 out[3*n-3]=VdstrictGB(L[1],d,v); 1080 for (i=n-1; i>=1; i--) 1081 { 1082 out[3*i-2]=nrows(out[3*i]); 1083 v=0; 1084 for (j=1; j<=out[3*i-2]; j++) 1085 { 1086 mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); 1087 v[j]=Vddeg(mem,d, out[3*i+2]); 1088 } 1089 out[3*i-1]=v; 1090 if (i!=1) 1091 { 1092 out[3*i-3]=transpose(syz(transpose(out[3*i]))); 1093 if (out[3*i-3]!=matrix(0,nrows(out[3*i-3]),ncols(out[3*i-3]))) 1094 { 1095 out[3*i-3]=VdstrictGB(out[3*i-3],d,out[3*i-1]); 1096 } 1097 else 1098 { 1099 out[3*i-3]=matrix(0,1,ncols(out[3*i-3])); 1100 out[3*i-4]=intvec(0); 1101 out[3*i-5]=int(0); 1102 for (j=i-2; j>=1; j--) 1103 { 1104 out[3*j]=matrix(0,1,1); 1105 out[3*j-1]=intvec(0); 1106 out[3*j-2]=int(0); 1107 } 1108 break; 1109 } 1110 } 1111 } 1112 outall[1]=out; 1113 outall[2]=list(list(out[3*n-3],out[3*n-1])); 1114 return(outall); 1115 } 1116 out=shortexactpieces(L); 1117 list rem; 1118 if (v!=intvec(0:size(v))) 1119 { 1120 rem=shortexactpiecestoVdstrict(out,d,v); 1121 } 1122 else 1123 { 1124 rem=shortexactpiecestoVdstrict(out,d); 1125 } 1126 out=Vdstrictdoublecompexes(rem[1],d); 1127 out=assemblingdoublecomplexes(out); 1128 out=totalcomplex(out); 1129 outall[1]=out; 1130 outall[2]=rem[2]; 1131 return (outall); 1132 } 1133 1134 //////////////////////////////////////////////////////////////////////////////// 1135 1136 proc derhamcohomology(list L) 1137 { 1138 def R=basering; 1139 int n=nvars(R);int le=2*size(L)+n-1; 1140 def W=makeWeyl(n); 1141 setring W; 1142 list man=ringlist(W); 1143 if (n==1) 1144 { 1145 man[2][1]="x(1)"; 1146 man[2][2]="D(1)"; 1147 def Wi=ring(man); 1148 setring Wi; 1149 kill W; 1150 def W=Wi; 1151 setring W; 1152 list man=ringlist(W); 1153 } 1154 man[2][size(man[2])+1]="s";; 1155 man[3][3]=man[3][2]; 1156 man[3][2]=list("dp",intvec(1)); 1157 matrix N=UpOneMatrix(size(man[2])); 1158 man[5]=N; 1159 matrix M[1][1]; 1160 man[6]=transpose(concat(transpose(concat(man[6],M)),M)); 1161 def Ws=ring(man); setring R; int r=size(L); int i; int j;int k; int l; int count; list Fi; list subsets; list maxnum; list bernsteinpolys; list annideals; list minint; list diffmaps; 1162 for (i=1; i<=r; i++) 1163 { 1164 Fi[i]=list(); bernsteinpolys[i]=list(); annideals[i]=list(); subsets[i]=list(); 1165 maxnum[i]=list(); 1166 Fi[1][i]=L[i]; 1167 maxnum[1][i]=i; 1168 subsets[1][i]=intvec(i); 1169 } 1170 intvec v; 1171 for (i=2; i<=r; i++) 1172 { 1173 count=1; 1174 for (j=1; j<=size(Fi[i-1]);j++) 1175 { 1176 for (k=maxnum[i-1][j]+1; k<=r; k++) 1177 { 1178 maxnum[i][count]=k; 1179 v=subsets[i-1][j],k; 1180 subsets[i][count]=v; 1181 Fi[i][count]=lcm(Fi[i-1][j],L[k]);///////// 1182 count=count+1; 1183 } 1184 } 1185 } 1186 for (i=1; i<=r; i++) 1187 { 1188 for (j=1; j<=size(Fi[i]); j++) 1189 { 1190 bernsteinpolys[i][j]=bfct(Fi[i][j])[1]; 1191 for (k=1; k<=ncols(bernsteinpolys[i][j]); k++) 1192 { 1193 if (isInt(number(bernsteinpolys[i][j][k]))==1) 1194 { 1195 minint[size(minint)+1]=int(bernsteinpolys[i][j][k]); 1196 } 1197 } 1198 def D=Sannfs(Fi[i][j]); 1199 setring Ws; 1200 annideals[i][j]=fetch(D,LD); 1201 kill D; 1202 setring R; 1203 } 1204 } 1205 int m=Min(minint); 1206 list zw; 1207 for (i=1; i<r; i++) 1208 { 1209 diffmaps[i]=matrix(0,size(subsets[i]),size(subsets[i+1])); 1210 for (j=1; j<=size(subsets[i]); j++) 1211 { 1212 for (k=1; k<=size(subsets[i+1]); k++) 1213 { 1214 zw=mysubset(subsets[i][j],subsets[i+1][k]); 1215 diffmaps[i][j,k]=zw[2]*(L[zw[1]]/gcd(L[zw[1]],Fi[i][j]))^(-m); 1216 } 1217 } 1218 } 1219 diffmaps[r]=matrix(0,1,1); 1220 setring Ws; 1221 for (i=1; i<=r; i++) 1222 { 1223 for (j=1; j<=size(annideals[i]); j++) 1224 { 1225 annideals[i][j]=subst(annideals[i][j],s,m); 1226 } 1227 } 1228 setring W; 1229 list annideals=imap(Ws,annideals); 1230 list diffmaps=fetch(R,diffmaps); 1231 list fortoVdstrict; 1232 ideal IFourier=var(n+1); 1233 for (i=2;i<=n;i++) 1234 { 1235 IFourier=IFourier,var(n+i); 1236 } 1237 for (i=1; i<=n;i++) 1238 { 1239 IFourier=IFourier,-var(i); 1240 } 1241 map cFourier=W,IFourier; 1242 matrix sup; 1243 for (i=1; i<=r; i++) 1244 { 1245 sup=matrix(annideals[i][1]); 1246 fortoVdstrict[2*i-1]=transpose(cFourier(sup)); 1247 for (j=2; j<=size(annideals[i]); j++) 1248 { 1249 sup=matrix(annideals[i][j]); 1250 fortoVdstrict[2*i-1]=dsum(fortoVdstrict[2*i-1],transpose(cFourier(sup))); 1251 } 1252 sup=diffmaps[i]; 1253 fortoVdstrict[2*i]=cFourier(sup); 1254 } 1255 list rem=toVdstrictfreecomplex(fortoVdstrict); 1256 list newcomplex=rem[1]; 1257 list minmaxk=globalbfun(rem[2]); 1258 if (size(minmaxk)==0) 1259 { 1260 return (0); 1261 } 1262 list truncatedcomplex; list shorten; list upto; 1263 for (i=1; i<=size(newcomplex) div 3; i++) 1264 { 1265 shorten[3*i-1]=list(); 1266 for (j=1; j<=size(newcomplex[3*i-1]); j++) 1267 { 1268 shorten[3*i-1][j]=list(minmaxk[1]-newcomplex[3*i-1][j]+1,minmaxk[2]-newcomplex[3*i-1][j]+1); 1269 upto[size(upto)+1]=shorten[3*i-1][j][2]; 1270 if (shorten[3*i-1][j][2]<=0) 1271 { 1272 shorten[3*i-1][j]=list(); 1273 } 1274 else 1275 { 1276 if (shorten[3*i-1][j][1]<=0) 1277 { 1278 shorten[3*i-1][j][1]=1; 1279 } 1280 } 1281 } 1282 } 1283 int iupto=Max(upto); 1284 if (iupto<=0) 1285 { 1286 /////die Kohomologie ist dann überall 0, muss noch entsprechend ausgegeben werden 1287 } 1288 list allpolys; 1289 allpolys[1]=list(1); 1290 list minvar; 1291 minvar[1]=list(1); 1292 for (i=1; i<=iupto-1; i++) 1293 { 1294 allpolys[i+1]=list(); 1295 minvar[i+1]=list(); 1296 for (k=1; k<=size(allpolys[i]); k++) 1297 { 1298 for (j=minvar[i][k]; j<=nvars(W) div 2; j++) 1299 { 1300 allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j); 1301 minvar[i+1][size(minvar[i+1])+1]=j; 1302 } 1303 } 1304 } 1305 list keepformatrix;list sizetruncom;int stc;list fortrun; 1306 for (i=1; i<=size(newcomplex) div 3; i++) 1307 { 1308 truncatedcomplex[2*i-1]=list(); 1309 sizetruncom[2*i-1]=list(); 1310 sizetruncom[2*i]=list(); 1311 truncatedcomplex[2*i]=newcomplex[3*i]; 1312 v=0;count=0; 1313 sizetruncom[2*i][1]=0; 1314 for (j=1; j<=newcomplex[3*i-2]; j++) 1315 { 1316 if (size(shorten[3*i-1][j])!=0) 1317 { 1318 fortrun=sublist(allpolys,shorten[3*i-1][j][1],shorten[3*i-1][j][2]); 1319 truncatedcomplex[2*i-1][size(truncatedcomplex[2*i-1])+1]=fortrun[1]; 1320 count=count+fortrun[2]; 1321 sizetruncom[2*i-1][size(sizetruncom[2*i-1])+1]=list(int(shorten[3*i-1][j][1])-1,int(shorten[3*i-1][j][2])-1); 1322 sizetruncom[2*i][size(sizetruncom[2*i])+1]=count; 1323 if (v!=0) 1324 { 1325 v[size(v)+1]=j; 1326 } 1327 else 1328 { 1329 v[1]=j; 1330 } 1331 } 1332 } 1333 1334 if (v!=0) 1335 { 1336 truncatedcomplex[2*i]=submat(truncatedcomplex[2*i],v,1..ncols(truncatedcomplex[2*i])); 1337 if (i!=1) 1338 { 1339 truncatedcomplex[2*(i-1)]=submat(truncatedcomplex[2*(i-1)],1..nrows(truncatedcomplex[2*(i-1)]),v); 1340 } 1341 } 1342 else 1343 { 1344 truncatedcomplex[2*i]=matrix(0,1,ncols(truncatedcomplex[2*i])); 1345 if (i!=1) 1346 { 1347 truncatedcomplex[2*(i-1)]=matrix(0,nrows(truncatedcomplex[2*(i-1)]),1); 1348 } 1349 } 1350 } 1351 int b;int d;poly form;poly lform; poly nform;int ideg;int kplus; int lplus; 1352 for (i=1; i<size(truncatedcomplex) div 2; i++) 1353 { 1354 M=matrix(0,max(1,sizetruncom[2*i][size(sizetruncom[2*i])]),sizetruncom[2*i+2][size(sizetruncom[2*i+2])]); 1355 for (k=1; k<=size(truncatedcomplex[2*i-1]);k++) 1356 { 1357 for (l=1; l<=size(truncatedcomplex[2*(i+1)-1]); l++) 1358 { 1359 if (size(sizetruncom[2*i])!=1)//? 1360 { 1361 for (j=1; j<=size(truncatedcomplex[2*i-1][k]); j++) 1362 { 1363 for (b=1; b<=size(truncatedcomplex[2*i-1][k][j]); b++) 1364 { 1365 form=truncatedcomplex[2*i-1][k][j][b][1]*truncatedcomplex[2*i][k,l]; 1366 while (form!=0) 1367 { 1368 lform=lead(form); 1369 v=leadexp(lform); 1370 v=v[1..n]; 1371 if (v==(0:n)) 1372 { 1373 ideg=deg(lform)-sizetruncom[2*(i+1)-1][l][1]; 1374 if (ideg>=0) 1375 { 1376 for (d=1; d<=size(truncatedcomplex[2*(i+1)-1][l][ideg+1]);d++) 1377 { 1378 if (leadmonom(lform)==truncatedcomplex[2*(i+1)-1][l][ideg+1][d][1]) 1379 { 1380 M[sizetruncom[2*i][k]+truncatedcomplex[2*i-1][k][j][b][2],sizetruncom[2*(i+1)][l]+truncatedcomplex[2*(i+1)-1][l][ideg+1][d][2]]=leadcoef(lform); 1381 break; 1382 } 1383 } 1384 } 1385 } 1386 form=form-lform; 1387 } 1388 } 1389 } 1390 } 1391 } 1392 } 1393 truncatedcomplex[2*i]=M; 1394 truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])]; 1395 } 1396 truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])]; 1397 if (truncatedcomplex[2*i-1]!=0) 1398 { 1399 truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i-1],1); 1400 } 1401 setring R; 1402 list truncatedcomplex=imap(W,truncatedcomplex); 1403 list derhamhom=findhomology(truncatedcomplex,le); 1404 return (derhamhom); 1405 } 1406 1407 /////////////////////////////////// 1408 static proc sublist(list L, int m, int n) 1409 { 1410 list out; 1411 int i; int j; 1412 int count; 1413 for (i=m; i<=n; i++) 1414 { 1415 out[size(out)+1]=list(); 1416 for (j=1; j<=size(L[i]); j++) 1417 { 1418 count=count+1; 1419 out[size(out)][j]=list(L[i][j],count); 1420 } 1421 } 1422 list o=list(out,count); 1423 return(o); 1424 } 1425 1426 ////////////////////////////////////////////////////////////////////////// 1427 static proc mysubset(intvec L, intvec M) 1428 { 1429 int i; 1430 int j=1; 1431 list position=(M[size(M)],(-1)^(size(L))); 1432 for (i=1; i<=size(L); i++) 1433 { 1434 if (L[i]!=M[j]) 1435 { 1436 if (L[i]!=M[j+1] or j!=i) 1437 { 1438 return (L[i],0); 1439 } 1440 else 1441 { 1442 position=(M[i],(-1)^(i-1)); 1443 j=j+i; 1444 } 1445 } 1446 j=j+1; 1447 } 1448 return (position); 1449 } 1450 1451 1452 1453 //////////////////////////////////////////////////////////////////////////// 1454 1455 proc globalbfun(list L) 1456 { 1457 int i; int j; 1458 def W=basering; 1459 int n=nvars(W) div 2; 1460 list G0; 1461 ideal I; 1462 for (j=1; j<=size(L); j++) 1463 { 1464 G0[j]=list(); 1465 for (i=1; i<=ncols(L[j][1]); i++) 1466 { 1467 G0[j][i]=I; 1468 } 1469 } 1470 list out; 1471 for (j=1; j<=size(L); j++) 1472 { 1473 for (i=1; i<=nrows(L[j][1]); i++) 1474 { 1475 out=Vddeg(submat(L[j][1],i,(1..ncols(L[j][1]))),n,L[j][2],1); 1476 G0[j][out[2]][size(G0[j][out[2]])+1]=(out[1]); 1477 } 1478 } 1479 list Data=ringlist(W); 1480 for (i=1; i<=n; i++) 1481 { 1482 Data[2][2*n+i]=Data[2][i]; 1483 Data[2][3*n+i]=Data[2][n+i]; 1484 Data[2][i]="v("+string(i)+")"; 1485 Data[2][n+i]="w("+string(i)+")"; 1486 } 1487 Data[3][1][1]="M"; 1488 intvec mord=(0:16*n^2); 1489 mord[1..2*n]=(1:2*n); 1490 mord[6*n+1..8*n]=(1:2*n); 1491 for (i=0; i<=2*n-2; i++) 1492 { 1493 mord[(3+i)*4*n-i]=-1; 1494 mord[(2*n+2+i)*4*n-2*n-i]=-1; 1495 } 1496 Data[3][1][2]=mord;//ordering mh????????? 1497 matrix Ones=UpOneMatrix(4*n); 1498 Data[5]=Ones; 1499 matrix con[2*n][2*n]; 1500 Data[6]=transpose(concat(con,transpose(concat(con,Data[6])))); 1501 1502 def Wuv=ring(Data); 1503 setring Wuv; 1504 list G0=imap(W,G0); list G3; poly lterm;intvec lexp; 1505 list G1; list G2; intvec e; intvec f; int kapp; int k; int l; poly h; ideal I; 1506 for (l=1; l<=size(G0); l++) 1507 { 1508 G1[l]=list(); G2[l]=list(); G3[l]=list(); 1509 for (i=1; i<=size(G0[l]); i++) 1510 { 1511 for (j=1; j<=ncols(G0[l][i]);j++) 1512 { 1513 G0[l][i][j]=mhom(G0[l][i][j]); 1514 } 1515 for (j=1; j<=nvars(Wuv) div 4; j++) 1516 { 1517 G0[l][i][size(G0[l][i])+1]=1-v(j)*w(j); 1518 } 1519 G1[l][i]=std(G0[l][i]); 1520 G2[l][i]=I; 1521 G3[l][i]=list(); 1522 for (j=1; j<=ncols(G1[l][i]); j++) 1523 { 1524 e=leadexp(G1[l][i][j]); 1525 f=e[1..2*n]; 1526 if (f==intvec(0:(2*n))) 1527 { 1528 for (k=1; k<=n; k++) 1529 { 1530 kapp=-e[2*n+k]+e[3*n+k]; 1531 if (kapp>0) 1532 { 1533 G1[l][i][j]=(x(k)^kapp)*G1[l][i][j]; 1534 } 1535 if (kapp<0) 1536 { 1537 G1[l][i][j]=(D(k)^(-kapp))*G1[l][i][j]; 1538 } 1539 } 1540 G2[l][i][size(G2[l][i])+1]=G1[l][i][j]; 1541 G3[l][i][size(G3[l][i])+1]=list(); 1542 while (G1[l][i][j]!=0) 1543 { 1544 lterm=lead(G1[l][i][j]); 1545 G1[l][i][j]=G1[l][i][j]-lterm; 1546 lexp=leadexp(lterm); 1547 lexp=lexp[2*n+1..3*n]; 1548 G3[l][i][size(G3[l][i])][size(G3[l][i][size(G3[l][i])])+1]=list(lexp,leadcoef(lterm)); 1549 } 1550 1551 } 1552 } 1553 } 1554 } 1555 ring r=0,(s(1..n)),dp; 1556 ideal I; 1557 map G3forr=Wuv,I; 1558 list G3=G3forr(G3); 1559 poly fs; 1560 poly gs; 1561 int a; 1562 list G4; 1563 for (l=1; l<=size(G3); l++) 1564 { 1565 G4[l]=list(); 1566 for (i=1; i<=size(G3[l]);i++) 1567 { 1568 G4[l][i]=I; 1569 for (j=1; j<=size(G3[l][i]); j++) 1570 { 1571 fs=0; 1572 for (k=1; k<=size(G3[l][i][j]); k++) 1573 { 1574 gs=1; 1575 for (a=1; a<=n; a++) 1576 { 1577 if (G3[l][i][j][k][1][a]!=0) 1578 { 1579 gs=gs*permutevar(list(G3[l][i][j][k][1][a]),a); 1580 } 1581 } 1582 gs=gs*G3[l][i][j][k][2]; 1583 fs=fs+gs; 1584 } 1585 G4[l][i]=G4[l][i],fs; 1586 } 1587 } 1588 } 1589 if (n==1) 1590 { 1591 ring rnew=0,t,dp; 1592 } 1593 else 1594 { 1595 ring rnew=0,(t,s(2..n)),dp; 1596 } 1597 ideal Iformap; 1598 Iformap[1]=t; 1599 poly forel=1; 1600 for (i=2; i<=n; i++) 1601 { 1602 Iformap[1]=Iformap[1]-s(i); 1603 Iformap[i]=s(i); 1604 forel=forel*s(i); 1605 } 1606 map rtornew=r,Iformap; 1607 list G4=rtornew(G4); 1608 list getintvecs=fetch(W,L); 1609 ideal J; 1610 option(redSB); 1611 for (l=1; l<=size(G4); l++) 1612 { 1613 J=1; 1614 for (i=1; i<=size(G4[l]); i++) 1615 { 1616 G4[l][i]=eliminate(G4[l][i],forel); 1617 G4[l][i]=subst(G4[l][i],t,t-getintvecs[l][2][i]); 1618 J=intersect(J,G4[l][i]); 1619 } 1620 G4[l]=poly(std(J)[1]); 1621 } 1622 list minmax=minmaxintroot(G4);//besser factorize nehmen 1623 // Fall: keine Nullstelle muss noch weiter beruecksichtigt werden 1624 return(minmax); 1625 } 1626 1627 1628 1629 1630 ////////////////////////////////////////////////////////////////////////// 1631 1632 proc minmaxintroot(list L); 1633 { 1634 int i; int j; int k; int l; int sa; int s; number d; poly f; poly rest; list a0; list possk; list alldiv; intvec e; 1635 possk[1]=list(); 1636 for (i=1; i<=size(L); i++) 1637 { 1638 d=1; 1639 f=L[i]; 1640 while (f!=0) 1641 { 1642 rest=lead(f); 1643 d=d*denominator(leadcoef(rest)); 1644 f=f-rest; 1645 } 1646 e=leadexp(rest); 1647 if (e[1]!=0) 1648 { 1649 rest=rest/(t^(e[1])); 1650 possk[1][size(possk[1])+1]=i; 1651 } 1652 a0[i]=int(absValue(d*rest)); 1653 } 1654 int m=Max(a0); 1655 for (i=2; i<=m+1; i++) 1656 { 1657 possk[i]=list(); 1658 } 1659 list allprimefac; 1660 for (i=1; i<=size(L); i++) 1661 { 1662 allprimefac=primefactors(a0[i]); 1663 alldiv=1; 1664 possk[2][size(possk[2])+1]=i; 1665 1666 for (j=1; j<=size(allprimefac[1]); j++) 1667 { 1668 s=size(alldiv); 1669 for (k=1; k<=s; k++) 1670 { 1671 for (l=1; l<=allprimefac[2][j]; l++) 1672 { 1673 alldiv[size(alldiv)+1]=alldiv[k]*allprimefac[1][j]^l; 1674 possk[alldiv[size(alldiv)]+1][size(possk[alldiv[size(alldiv)]+1])+1]=i; 1675 } 1676 } 1677 } 1678 } 1679 int mink; 1680 int maxk; 1681 int indi; 1682 for (i=m+1; i>=1; i--) 1683 { 1684 if (size(possk[i])!=0) 1685 { 1686 for (j=1; j<=size(possk[i]); j++) 1687 { 1688 if (subst(L[possk[i][j]],t,(i-1))==0) 1689 { 1690 maxk=i-1; 1691 indi=1; 1692 break; 1693 } 1694 } 1695 } 1696 if (maxk!=0) 1697 { 1698 break; 1699 } 1700 } 1701 int indi2; 1702 for (i=m+1; i>=1; i--) 1703 { 1704 if (size(possk[i])!=0) 1705 { 1706 for (j=1; j<=size(possk[i]); j++) 1707 { 1708 if (subst(L[possk[i][j]],t,-(i-1))==0) 1709 { 1710 mink=-i+1; 1711 indi2=1; 1712 break; 1713 } 1714 } 1715 } 1716 if (mink!=0) 1717 { 1718 break; 1719 } 1720 } 1721 list mima=mink,maxk; 1722 if (indi==0) 1723 { 1724 if (indi2==0) 1725 { 1726 mima=list();//es gibt keine ganzzahlige NS 1727 } 1728 else 1729 { 1730 mima[2]=mima[1]; 1731 } 1732 } 1733 else 1734 { 1735 if (indi2==0) 1736 { 1737 mima[1]=mima[2]; 1738 } 1739 } 1740 return (mima); 1741 } 1742 1743 /////////////////////////////////////////////////////// 1744 1745 1746 proc findhomology(list L, int le) 1747 { 1748 int li; 1749 matrix M; matrix N; 1750 matrix N1; 1751 matrix lift1; 1752 list out; 1753 int i; 1754 option (redSB); 1755 for (i=2; i<=size(L); i=i+2) 1756 { 1757 if (L[i-1]==0) 1758 { 1759 li=0; 1760 out[i div 2]=0; 1761 } 1762 else 1763 { 1764 1765 if (li==0) 1766 { 1767 1768 1769 li=L[i-1]; 1770 N1=transpose(syz(transpose(L[i]))); 1771 out[i div 2]=matrix(transpose(syz(transpose(N1)))); 1772 out[i div 2]=transpose(matrix(std(transpose(out[i div 2])))); 1773 1774 } 1775 1776 else 1777 { 1778 1779 1780 li=L[i-1]; 1781 N1=transpose(syz(transpose(L[i]))); 1782 N=transpose(syz(transpose(N1))); 1783 lift1=matrixlift(N1,L[i-2]); 1784 out[i div 2]=transpose(concat(transpose(lift1),transpose(N))); 1785 out[i div 2]=transpose(matrix(std(transpose(out[i div 2])))); 1786 } 1787 } 1788 if (out[i div 2]!=matrix(0,1,ncols(out[i div 2]))) 1789 { 1790 out[i div 2]=ncols(out[i div 2])-nrows(out[i div 2]); 1791 } 1792 else 1793 { 1794 out[i div 2]=ncols(out[i div 2]); 1795 } 1796 } 1797 if (size(out)>le) 1798 { 1799 out=delete(out,1); 1800 } 1801 return(out); 1802 } 1803 1804 1805 1806 1807 ///////////////////////////////////////////////////////////////////// 1808 1809 static proc mhom(poly f) 1810 { 1811 poly g; 1812 poly l; 1813 poly add; 1814 intvec e; 1815 list minint; 1816 list remf; 1817 int i; 1818 int j; 1819 int n=nvars(basering) div 4; 1820 if (f==0) 1821 { 1822 return(f); 1823 } 1824 while (f!=0) 1825 { 1826 l=lead(f); 1827 e=leadexp(l); 1828 remf[size(remf)+1]=list(); 1829 remf[size(remf)][1]=l; 1830 for (i=1; i<=n; i++) 1831 { 1832 remf[size(remf)][i+1]=-e[2*n+i]+e[3*n+i]; 1833 if (size(minint)<i) 1834 { 1835 minint[i]=list(); 1836 } 1837 minint[i][size(minint[i])+1]=-e[2*n+i]+e[3*n+i]; 1838 } 1839 f=f-l; 1840 } 1841 for (i=1; i<=n; i++) 1842 { 1843 minint[i]=Min(minint[i]); 1844 } 1845 for (i=1; i<=size(remf); i++) 1846 { 1847 add=remf[i][1]; 1848 for (j=1; j<=n; j++) 1849 { 1850 add=v(j)^(remf[i][j+1]-minint[j])*add; 1851 } 1852 g=g+add; 1853 } 1854 return (g); 1855 } 1856 1857 1858 1859 1860 ////////////////////////////////////////////////////////////////////////// 1861 1862 static proc permutevar(list L,int n) 1863 { 1864 if (typeof(L[1])=="intvec") 1865 { 1866 intvec v=L[1]; 1867 } 1868 else 1869 { 1870 intvec v=(1:L[1]),(0:L[1]); 1871 } 1872 int i;int k; int indi=0; 1873 int j; 1874 int s=size(v); 1875 poly e; 1876 intvec fore; 1877 for (i=2; i<=size(v); i=i+2) 1878 { 1879 if (v[i]!=0) 1880 { 1881 j=i+1; 1882 while (v[j]!=0) 1883 { 1884 j=j+1; 1885 } 1886 v[i]=0; 1887 v[j]=1; 1888 fore=0; 1889 indi=0; 1890 for (k=1; k<=size(v); k++) 1891 { 1892 if (k!=i and k!=j) 1893 { 1894 if (indi==0) 1895 { 1896 indi=1; 1897 fore[1]=v[k]; 1898 } 1899 else 1900 { 1901 fore[size(fore)+1]=v[k]; 1902 } 1903 } 1904 } 1905 e=e-(j-i)*permutevar(list(fore),n); 1906 } 1907 } 1908 e=e+s(n)^(size(v) div 2); 1909 return (e); 1910 } 1911 1912 /////////////////////////////////////////////////////////////////////////////// 1913 static proc max(int i,int j) 1914 { 1915 if(i>j){return(i);} 1916 return(j); 1917 } 1918 1919 //////////////////////////////////////////////////////////////////////////////////// 1920 version="$Id$"; 2 version="version derham.lib 4.0.0.0 Jun_2013 "; 1921 3 category="Noncommutative"; 1922 4 info=" … … 1954 36 LIB "qhmoduli.lib"; 1955 37 LIB "general.lib"; 1956 LIB "dmod.lib";38 //LIB "dmod.lib"; 1957 39 LIB "bfun.lib"; 1958 40 LIB "dmodapp.lib"; … … 2643 725 annihilators we computed for L[i] ->this will ensure we can compute 2644 726 the maps of the MV complex*/ 2645 minroot= minIntRoot(findminintroot);727 minroot=Derham::minIntRoot(findminintroot); 2646 728 for (j=1; j<=i; j++) 2647 729 { … … 2729 811 s-parametric annihilators we computed for L[i] ->this will 2730 812 ensure we can compute the maps of the MV complex*/ 2731 minroot= minIntRoot(findminintroot);813 minroot=Derham::minIntRoot(findminintroot); 2732 814 for (j=1; j<=i; j++) 2733 815 { … … 2969 1051 } 2970 1052 def B=basering; 2971 int n=nvars(B) div 2 + 1;//+1 m üsste stimmen! bitte kontrollieren!1053 int n=nvars(B) div 2 + 1;//+1 muesste stimmen! bitte kontrollieren! 2972 1054 int d=nvars(B) div 2; 2973 1055 int r=size(L) div 2; … … 3266 1348 matrix zerom; 3267 1349 list rofA; 3268 for (i=3; i<=n+3; i++)////////////////////////////////////////////////////////////////////////////n und si m üssen noch definiert werden1350 for (i=3; i<=n+3; i++)////////////////////////////////////////////////////////////////////////////n und si muessen noch definiert werden 3269 1351 { 3270 1352 if (size(RES)>=i) … … 5796 3878 component of submat(L[j][1],i,(1..ncols(L[j][1])))*/ 5797 3879 i1=(1..ncols(L[j][1])); 5798 out=VdDegTilde(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1);//hier k önnte es evtl noch einen Fehler geben!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!3880 out=VdDegTilde(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1);//hier koennte es evtl noch einen Fehler geben!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 5799 3881 f=L[j][1][i,out[2]]; 5800 3882 if (out[2]==1) … … 6148 4230 for (i=1; i<=size(G4); i++) 6149 4231 { 6150 minmax[i]= minIntRoot(G4[i],1);4232 minmax[i]=Derham::minIntRoot(G4[i],1); 6151 4233 if (size(minmax[i])!=0) 6152 4234 { … … 6196 4278 else 6197 4279 { 6198 out[i div 2]=0;////achtung ge ändert??????????????????????????????????????????????????!!!!!!!!!!!!!!!!!!!!!!!!!war mal out[i-1]4280 out[i div 2]=0;////achtung geaendert??????????????????????????????????????????????????!!!!!!!!!!!!!!!!!!!!!!!!!war mal out[i-1] 6199 4281 im=L[i-1]; 6200 4282 } … … 6534 4616 and submat(Fnew,l,intvec(1..ncols(Fnew)))!=zeromat) 6535 4617 { 6536 //print("Reduzierung des V_d-Grades n ötig");4618 //print("Reduzierung des V_d-Grades noetig"); 6537 4619 /*We need to reduce the V_d-degree. First we homogenize the 6538 4620 lth row of Fnew*/
Note: See TracChangeset
for help on using the changeset viewer.