Changeset 0e8a5a in git
 Timestamp:
 Sep 26, 2013, 2:16:29 PM (10 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
 Children:
 12d61707b29797f7bf3e42e5b9f0a01e937da42b
 Parents:
 89b2d289ea9d528c84c54a5e203ee55bb49e600c
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/derham.lib
r89b2d28 r0e8a5a 1 /////////////////////////////////////////////////////////////////////////////// //2 version=" version derham.lib 4.0.0.0 Jun_2013 "; // $Id$1 /////////////////////////////////////////////////////////////////////////////// 2 version=""; 3 3 category="Noncommutative"; 4 4 info=" 5 5 LIBRARY: derham.lib Computation of deRham cohomology 6 7 6 AUTHORS: Cornelia Rottner, rottner@mathematik.unikl.de 8 9 OVERVIEW: 10 A library for computing the de Rham cohomology of complements of complex affine 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, #[i1] , #[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 Vdstrict 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:(ncMsize(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=l1; 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=flead(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 1xrmatrix, d int, v intvec of size r 355 RETURN:int; the Vddegree 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 0modul 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 Vdstrict 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=s1; 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 22; 725 for (k=1; k<=(size(L)+nvars(basering) div 23); 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[i1][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[i1][k][4]=matrix(out[i1][k][4],nrows(out[i1][k][4]),out[i][k][1]); 923 } 924 for (k=newj+1; k<=oldj; k++) 925 { 926 out[i1][k]=delete(out[i1][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*i2]=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*i2]=0; 959 v=0; 960 } 961 962 for (j=i+1; j<=size(L); j++) 963 { 964 if (size(L[j])>=ji+1) 965 { 966 out[3*i2]=out[3*i2]+L[j][ji+1][1]; 967 if (emp==0) 968 { 969 rem1=L[j][ji+1][2]; 970 emp=1; 971 } 972 else 973 { 974 rem1=rem1,L[j][ji+1][2]; 975 } 976 v[size(v)+1]=L[j][ji+1][1]; 977 } 978 else 979 { 980 v[size(v)+1]=0; 981 } 982 } 983 out[3*i1]=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*i2])!=0) 994 { 995 o1=out[3*i2]; 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][ji+1]!=0) 1014 { 1015 for (k=c+1; k<=c+remsize[i][ji+1]; k++) 1016 { 1017 for (l=d+1; l<=d+remsize[i+1][ji];l++) 1018 { 1019 M[k,l]=L[j][ji+1][3][(kc),(ld)]; 1020 } 1021 } 1022 d=d+remsize[i+1][ji]; 1023 if (remsize[i+1][ji+1]!=0) 1024 { 1025 for (k=c+1; k<=c+remsize[i][ji+1]; k++) 1026 { 1027 for (l=d+1; l<=d+remsize[i+1][ji+1];l++) 1028 { 1029 M[k,l]=L[j][ji+1][4][kc,ld]; 1030 } 1031 } 1032 c=c+remsize[i][ji+1]; 1033 } 1034 } 1035 else 1036 { 1037 d=d+remsize[i+1][ji]; 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*n1]=v; 1077 out[3*n2]=ncols(L[1]); 1078 out[3*n]=L[2]; 1079 out[3*n3]=VdstrictGB(L[1],d,v); 1080 for (i=n1; i>=1; i) 1081 { 1082 out[3*i2]=nrows(out[3*i]); 1083 v=0; 1084 for (j=1; j<=out[3*i2]; 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*i1]=v; 1090 if (i!=1) 1091 { 1092 out[3*i3]=transpose(syz(transpose(out[3*i]))); 1093 if (out[3*i3]!=matrix(0,nrows(out[3*i3]),ncols(out[3*i3]))) 1094 { 1095 out[3*i3]=VdstrictGB(out[3*i3],d,out[3*i1]); 1096 } 1097 else 1098 { 1099 out[3*i3]=matrix(0,1,ncols(out[3*i3])); 1100 out[3*i4]=intvec(0); 1101 out[3*i5]=int(0); 1102 for (j=i2; j>=1; j) 1103 { 1104 out[3*j]=matrix(0,1,1); 1105 out[3*j1]=intvec(0); 1106 out[3*j2]=int(0); 1107 } 1108 break; 1109 } 1110 } 1111 } 1112 outall[1]=out; 1113 outall[2]=list(list(out[3*n3],out[3*n1])); 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)+n1; 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[i1]);j++) 1175 { 1176 for (k=maxnum[i1][j]+1; k<=r; k++) 1177 { 1178 maxnum[i][count]=k; 1179 v=subsets[i1][j],k; 1180 subsets[i][count]=v; 1181 Fi[i][count]=lcm(Fi[i1][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*i1]=transpose(cFourier(sup)); 1247 for (j=2; j<=size(annideals[i]); j++) 1248 { 1249 sup=matrix(annideals[i][j]); 1250 fortoVdstrict[2*i1]=dsum(fortoVdstrict[2*i1],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*i1]=list(); 1266 for (j=1; j<=size(newcomplex[3*i1]); j++) 1267 { 1268 shorten[3*i1][j]=list(minmaxk[1]newcomplex[3*i1][j]+1,minmaxk[2]newcomplex[3*i1][j]+1); 1269 upto[size(upto)+1]=shorten[3*i1][j][2]; 1270 if (shorten[3*i1][j][2]<=0) 1271 { 1272 shorten[3*i1][j]=list(); 1273 } 1274 else 1275 { 1276 if (shorten[3*i1][j][1]<=0) 1277 { 1278 shorten[3*i1][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<=iupto1; 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*i1]=list(); 1309 sizetruncom[2*i1]=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*i2]; j++) 1315 { 1316 if (size(shorten[3*i1][j])!=0) 1317 { 1318 fortrun=sublist(allpolys,shorten[3*i1][j][1],shorten[3*i1][j][2]); 1319 truncatedcomplex[2*i1][size(truncatedcomplex[2*i1])+1]=fortrun[1]; 1320 count=count+fortrun[2]; 1321 sizetruncom[2*i1][size(sizetruncom[2*i1])+1]=list(int(shorten[3*i1][j][1])1,int(shorten[3*i1][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*(i1)]=submat(truncatedcomplex[2*(i1)],1..nrows(truncatedcomplex[2*(i1)]),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*(i1)]=matrix(0,nrows(truncatedcomplex[2*(i1)]),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*i1]);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*i1][k]); j++) 1362 { 1363 for (b=1; b<=size(truncatedcomplex[2*i1][k][j]); b++) 1364 { 1365 form=truncatedcomplex[2*i1][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*i1][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=formlform; 1387 } 1388 } 1389 } 1390 } 1391 } 1392 } 1393 truncatedcomplex[2*i]=M; 1394 truncatedcomplex[2*i1]=sizetruncom[2*i][size(sizetruncom[2*i])]; 1395 } 1396 truncatedcomplex[2*i1]=sizetruncom[2*i][size(sizetruncom[2*i])]; 1397 if (truncatedcomplex[2*i1]!=0) 1398 { 1399 truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i1],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)^(i1)); 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*n2; i++) 1492 { 1493 mord[(3+i)*4*ni]=1; 1494 mord[(2*n+2+i)*4*n2*ni]=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]=1v(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,tgetintvecs[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=frest; 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,(i1))==0) 1689 { 1690 maxk=i1; 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,(i1))==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[i1]==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[i1]; 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[i1]; 1781 N1=transpose(syz(transpose(L[i]))); 1782 N=transpose(syz(transpose(N1))); 1783 lift1=matrixlift(N1,L[i2]); 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=fl; 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(ji)*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$"; 1921 category="Noncommutative"; 1922 info=" 1923 LIBRARY: derham.lib Computation of deRham cohomology 1924 1925 AUTHORS: Cornelia Rottner, rottner@mathematik.unikl.de 1926 1927 OVERVIEW: 1928 A library for computing the de Rham cohomology of complements of complex affine 11 1929 varieties. 12 1930 13 1931 14 REFERENCES: 15 [OT] Oaku, T.; Takayama, N.: Algorithms of Dmodules  restriction, tensor product, 16 localzation, and local cohomology groups , J. Pure Appl. Algebra 156, 2673081932 REFERENCES: 1933 [OT] Oaku, T.; Takayama, N.: Algorithms of Dmodules  restriction, tensor product, 1934 localzation, and local cohomology groups}, J. Pure Appl. Algebra 156, 267308 17 1935 (2001) 18 1936 [R] Rottner, C.: Computing de Rham Cohomology,diploma thesis (2012) 19 [W1] Walther, U.: Algorithmic computation of local cohomology modules and the local 20 cohomological dimension of algebraic varieties , J. Pure Appl. Algebra 139,1937 [W1] Walther, U.: Algorithmic computation of local cohomology modules and the local 1938 cohomological dimension of algebraic varieties}, J. Pure Appl. Algebra 139, 21 1939 303321 (1999) 22 [W2] Walther, U.: Algorithmic computation of de Rham Cohomology of Complements of 23 Complex Affine Varieties, J. Symbolic Computation 29, 796839 (2000) 1940 [W2] Walther, U.: Algorithmic computation of de Rham Cohomology of Complements of 1941 Complex Affine Varieties}, J. Symbolic Computation 29, 796839 (2000) 1942 [W3] Walther, U.: Computing the cup product structure for complements of complex 1943 affine varieties, J. Pure Appl. Algebra 164, 247273 (2001) 24 1944 25 1945 … … 39 1959 LIB "poly.lib"; 40 1960 LIB "schreyer.lib"; 1961 LIB "dmodloc.lib"; 1962 41 1963 42 1964 //////////////////////////////////////////////////////////////////////////////////// 43 1965 44 1966 proc deRhamCohomology(list L,list #) 45 "USAGE: deRhamCohomology(L[,choices]); L a list consisting of polynomials, choices 1967 "USAGE: deRhamCohomology(L[,choices]); L a list consisting of polynomials, choices 46 1968 optional list consisting of one up to three strings @* 47 1969 The optional strings may be one of the strings@* 48 'Vdres': compute the CartanEilenberg resolutions via V_dhomogenization 1970 'noCE': compute quasiisomorphic complexes without using CartanEilenberg 1971 resolutionsq@* 1972 'Vdres': compute quasiisomorphic complexes using CartanEilenberg 1973 resolutions; the CE resolutions are computed via V__dhomogenization 49 1974 and without using Schreyer's method @* 50 'Sres': compute the CartanEilenberg resolutions in the homogenized Weyl51 algebra usingSchreyer's method@*1975 'Sres': compute quasiisomorphic complexes using CartanEilenberg 1976 resolutions in the homogenized Weyl algebra via Schreyer's method@* 52 1977 one of the strings@* 53 'iterativeloc': compute localizations by factorizing the polynomials and 1978 'iterativeloc': compute localizations by factorizing the polynomials and 54 1979 sucessive localization of the factors @* 55 1980 'no iterativeloc': compute localizations by directly localizing the … … 60 1985 'exactroots' computes the minimal and maximal integer root of the global 61 1986 bfunction 62 The default is ' Sres', 'iterativeloc' and 'onlybounds'.1987 The default is 'noCE', 'iterativeloc' and 'onlybounds'. 63 1988 ASSUME: The basering must be a polynomial ring over the field of rational numbers@* 64 RETURN: list, where the ith entry is the (i1)st de Rham cohomology group of the 65 complement of the complex affine variety given by the polynomials in L 1989 RETURN: list, where the ith entry is the (i1)st de Rham cohomology group of the 1990 complement of the complex affine variety given by the polynomials in L 66 1991 EXAMPLE:example deRhamCohomology; shows an example 67 1992 " … … 76 2001 int n=nvars(R); 77 2002 int le=size(L)+n; 78 string Syzstring=" Sres";2003 string Syzstring="noCE"; 79 2004 int onlybounds=1; 2005 int diffforms; 80 2006 for (i=1; i<=size(#); i++) 81 { 82 if (#[i]=="Vdres") 83 { 84 Syzstring="Vdres"; 85 } 86 if (#[i]=="noiterativeloc") 87 { 88 recursiveloc=0; 89 } 90 if (#[i]=="exactroots") 91 { 92 onlybounds=0; 93 } 94 } 2007 { 2008 if (#[i]=="Sres") 2009 { 2010 Syzstring="Sres"; 2011 } 2012 if (#[i]=="Vdres") 2013 { 2014 Syzstring="Vdres"; 2015 } 2016 if (#[i]=="noiterativeloc") 2017 { 2018 recursiveloc=0; 2019 } 2020 if (#[i]=="exactroots") 2021 { 2022 onlybounds=0; 2023 } 2024 if (#[i]=="diffforms") 2025 { 2026 diffforms=1; 2027 } 2028 } 95 2029 for (i=1; i<=size(L); i++) 96 {97 if (L[i]==0)98 {99 L=delete(L,i);100 i=i1;101 }102 }2030 { 2031 if (L[i]==0) 2032 { 2033 L=delete(L,i); 2034 i=i1; 2035 } 2036 } 103 2037 if (size(L)==0) 104 {105 return (list(0));106 }2038 { 2039 return (list(0));//////////////////////////////////////////////////////////////////stimmt das jetzt?!?????????????????????????????????? 2040 } 107 2041 for (i=1; i<= size(L); i++) 108 {109 if (leadcoef(L[i])L[i]==0)110 {111 return(list(1));112 }113 }2042 { 2043 if (leadcoef(L[i])L[i]==0) 2044 { 2045 return(list(1)); ///////////////////////////////////////////////////////////////stimmt das jetzt?!???????????????????????????????????? 2046 } 2047 } 114 2048 if (size(L)==0) 115 {116 /*the complement of the variety given by the input is the whole space*/117 return(list(1));118 }2049 { 2050 /*the complement of the variety given by the input is the whole space*/ 2051 return(list(1)); 2052 } 119 2053 for (i=1; i<=size(L); i++) 120 { 121 if (typeof(L[i])!="poly") 122 { 123 print("The input list must consist of polynomials"); 124 retrun(); 125 } 126 } 2054 { 2055 if (typeof(L[i])!="poly") 2056 { 2057 print("The input list must consist of polynomials"); 2058 return(); 2059 } 2060 } 2061 if (size(L)==1 and Syzstring=="noCE") 2062 { 2063 Syzstring="Sres"; 2064 } 127 2065 /* 1st step: compute the MayerVietoris Complex and its Fourier transform*/ 128 2066 def W=MVComplex(L,recursiveloc);//new ring that contains the MV complex 129 2067 setring W; 130 2068 list fortoVdstrict=MV; 131 ideal IFourier=var(n+1); 132 for (i=2;i<=n;i++) 133 { 134 IFourier=IFourier,var(n+i); 135 } 136 for (i=1; i<=n;i++) 137 { 138 IFourier=IFourier,var(i); 139 } 140 map cFourier=W,IFourier; 141 matrix sup; 142 for (i=1; i<=size(MV); i++) 143 { 144 sup=fortoVdstrict[i]; 145 /*takes the Fourier transform of the MV complex*/ 146 fortoVdstrict[i]=cFourier(sup); 147 } 148 /* 2nd step: Compute a V_dstrict free complex that is quasiisomorphic to the 149 complex fortoVdstrict 150 The 1st entry of the list rem will be the quasiisomorphic complex, the 2nd 151 entry contains the cohomology modules and is needed for the computation of the 152 global bfunction*/ 153 list rem=toVdStrictFreeComplex(fortoVdstrict,Syzstring); 2069 if (diffforms==0) 2070 { 2071 ideal IFourier=var(n+1); 2072 for (i=2;i<=n;i++) 2073 { 2074 IFourier=IFourier,var(n+i); 2075 } 2076 for (i=1; i<=n;i++) 2077 { 2078 IFourier=IFourier,var(i); 2079 } 2080 map cFourier=W,IFourier; 2081 matrix sup; 2082 for (i=1; i<=size(MV); i++) 2083 { 2084 sup=fortoVdstrict[i]; 2085 /*takes the Fourier transform of the MV complex*/ 2086 fortoVdstrict[i]=cFourier(sup); 2087 } 2088 } 2089 /* 2nd step: Compute a V_dstrict free complex that is quasiisomorphic to the 2090 complex fortoVdstrict 2091 The 1st entry of the list rem will be the quasiisomorphic complex, the 2nd 2092 entry contains the cohomology modules and is needed for the computation of the 2093 global bfunction*/ 2094 if (Syzstring=="noCE") 2095 { 2096 list rem=quasiisomorphicVdComplex(fortoVdstrict,diffforms); 2097 list quasiiso=rem[3]; 2098 } 2099 else 2100 { 2101 list rem=toVdStrictFreeComplex(fortoVdstrict,Syzstring,diffforms); 2102 if (diffforms==1) 2103 { 2104 list quasiiso=list(matrix(1,1,1)); 2105 } 2106 } 154 2107 list newcomplex=rem[1]; 155 156 /* 3rd step: Compute the bounds for the minimal and maximal integer root of the 157 global bfunction of newcomplex(i.e. compute the lcm of the bfunctions of its158 cohomology modules)(if onlybouns=1). Else we compute the minimal and maximal159 integer root.160 161 If we compute only the bounds, we omit additional Groebner basis computations.162 However this leads to a higherdimensional truncated complex.163 164 Note that the cohomology modules are already contained in rem[2].165 minmaxk[1] and minmaxk[2] will contain the bounds resp exact roots.*/166 if ( onlybounds==1)167 {168 list minmaxk=globalBFun(rem[2],Syzstring);169 }2108 //////////////////////////////////////////////////////////////////////////////////// 2109 /* 3rd step: Compute the bounds for the minimal and maximal integer root of the 2110 global bfunction of newcomplex(i.e. compute the lcm of the bfunctions of its 2111 cohomology modules)(if onlybouns=1). Else we compute the minimal and maximal 2112 integer root. 2113 2114 If we compute only the bounds, we omit additional Groebner basis computations. 2115 However this leads to a higherdimensional truncated complex. 2116 2117 Note that the cohomology modules are already contained in rem[2]. 2118 minmaxk[1] and minmaxk[2] will contain the bounds resp exact roots.*/ 2119 if (diffforms==1) 2120 { 2121 list minmaxk=exactGlobalBFunIntegration(rem[2]); 2122 } 170 2123 else 171 { 172 list minmaxk=exactGlobalBFun(rem[2],Syzstring); 173 } 2124 { 2125 if (onlybounds==1) 2126 { 2127 list minmaxk=globalBFun(rem[2],Syzstring); 2128 } 2129 else 2130 { 2131 list minmaxk=exactGlobalBFun(rem[2],Syzstring); 2132 } 2133 } 174 2134 if (size(minmaxk)==0) 175 { 176 return (0); 177 } 178 /*4th step: Truncate the complex D_n/(x_1,...,x_n)\otimes C, (where 179 C=(C^i[m^i],d^i) is given by newcomplex, i.e. C^i=D_n^newcomplex[3*i2], 180 m^i=newcomplex[3*i1], d^i=newcomplex[3*i]), using Thm 5.7 in [W1]: 181 The truncated module D_n/(x_1,..,x_n)\otimes C[i] is generated by the set 182 (0,...,P_(i_j),0,...), where P_(i_j) is a monomial in C[D(1),...,D(n)] and 183 if it is placed in component k it holds that 184 minmaxk[1]m^i[k]<=deg(P_(i_j))<=minmaxk[2]m^i[k]*/ 2135 { 2136 return (0); 2137 } 2138 ///////////////////////////////////////////////////////////////////////////Bis hierhin angepasst 2139 /*4th step: Truncate the complex D_n/(x_1,...,x_n)\otimes C, (where 2140 C=(C^i[m^i],d^i) is given by newcomplex, i.e. C^i=D_n^newcomplex[3*i2], 2141 m^i=newcomplex[3*i1], d^i=newcomplex[3*i]), using Thm 5.7 in [W1]: 2142 The truncated module D_n/(x_1,..,x_n)\otimes C[i] is generated by the set 2143 (0,...,P_(i_j),0,...), where P_(i_j) is a monomial in C[D(1),...,D(n)] and 2144 if it is placed in component k it holds that 2145 minmaxk[1]m^i[k]<=deg(P_(i_j))<=minmaxk[2]m^i[k]*/ 185 2146 int k,l; 186 2147 list truncatedcomplex,shorten,upto; 187 2148 for (i=1; i<=size(newcomplex) div 3; i++) 188 { 189 shorten[3*i1]=list(); 190 for (j=1; j<=size(newcomplex[3*i1]); j++) 191 { 192 /*shorten[3*i1][j][k]=minmaxk[k]m^i[j]+1 (for k=1,2) if this value is 193 positive otherwise we will set it to be list(); 194 we added +1, because we will use a list, where we put in position l 195 polys of degree l+1*/ 196 shorten[3*i1][j]=list(minmaxk[1]newcomplex[3*i1][j]+1); 197 shorten[3*i1][j][2]=minmaxk[2]newcomplex[3*i1][j]+1; 198 upto[size(upto)+1]=shorten[3*i1][j][2]; 199 if (shorten[3*i1][j][2]<=0) 200 { 201 shorten[3*i1][j]=list(); 202 } 203 else 204 { 205 if (shorten[3*i1][j][1]<=0) 206 { 207 shorten[3*i1][j][1]=1; 208 } 209 } 210 } 211 } 2149 { 2150 shorten[3*i1]=list(); 2151 for (j=1; j<=size(newcomplex[3*i1]); j++) 2152 { 2153 /*shorten[3*i1][j][k]=minmaxk[k]m^i[j]+1 (for k=1,2) if this value is 2154 positive otherwise we will set it to be list(); 2155 . we added +1, because we will use a list, where we put in position l 2156 polys of degree l+1*/ 2157 shorten[3*i1][j]=list(minmaxk[1]newcomplex[3*i1][j]+1); 2158 if (diffforms==1) 2159 { 2160 shorten[3*i1][j][1]=1; 2161 } 2162 shorten[3*i1][j][2]=minmaxk[2]newcomplex[3*i1][j]+1; 2163 upto[size(upto)+1]=shorten[3*i1][j][2]; 2164 if (shorten[3*i1][j][2]<=0) 2165 { 2166 shorten[3*i1][j]=list(); 2167 } 2168 else 2169 { 2170 if (shorten[3*i1][j][1]<=0) 2171 { 2172 shorten[3*i1][j][1]=1; 2173 } 2174 } 2175 } 2176 } 212 2177 int iupto=Max(upto);//maximal degree +1 of the polynomials we have to consider 213 2178 if (iupto<=0) 214 {215 return(list(0));216 }2179 { 2180 return(list(0)); 2181 } 217 2182 list allpolys; 218 2183 /*allpolys[i] will consist list of all monomials in D(1),...,D(n) of degree i1*/ 219 2184 allpolys[1]=list(1); 220 2185 list minvar; 2186 list keepv; 221 2187 minvar[1]=list(1); 222 2188 for (i=1; i<=iupto1; i++) 223 { 224 allpolys[i+1]=list(); 225 minvar[i+1]=list(); 226 for (k=1; k<=size(allpolys[i]); k++) 227 { 228 for (j=minvar[i][k]; j<=nvars(W) div 2; j++) 229 { 230 allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j); 231 minvar[i+1][size(minvar[i+1])+1]=j; 232 } 233 } 234 } 2189 { 2190 allpolys[i+1]=list(); 2191 minvar[i+1]=list(); 2192 for (k=1; k<=size(allpolys[i]); k++) 2193 { 2194 for (j=minvar[i][k]; j<=nvars(W) div 2; j++) 2195 { 2196 if (diffforms==0) 2197 { 2198 allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j); 2199 } 2200 else 2201 { 2202 allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*x(j); 2203 } 2204 minvar[i+1][size(minvar[i+1])+1]=j; 2205 } 2206 } 2207 } 235 2208 list keepformatrix,sizetruncom,fortrun,fst; 236 2209 int count,stc; 237 2210 intvec v,forin; 238 2211 matrix subm; 2212 list keepcount; 2213 list passendespoly; 239 2214 /*now we compute the truncation*/ 240 2215 for (i=1; i<=size(newcomplex) div 3; i++) 241 { 242 /*truncatedcomplex[2*i1] will contain all the generators for the truncation 243 of D_n/(x(1),..,x(n))\otimes C[i]*/ 244 truncatedcomplex[2*i1]=list(); 245 sizetruncom[2*i1]=list(); 246 sizetruncom[2*i]=list(); 247 /*truncatedcomplex[2*i] will be the map trunc(D_n/(x(1),..,x(n))\otimes C[i]) 248 >trunc(D_n/(x(1),..,x(n))\otimes C[i+1])*/ 249 truncatedcomplex[2*i]=newcomplex[3*i]; 250 v=0;count=0; 251 sizetruncom[2*i][1]=0; 252 for (j=1; j<=newcomplex[3*i2]; j++) 253 { 254 if (size(shorten[3*i1][j])!=0) 255 { 256 fortrun=sublist(allpolys,shorten[3*i1][j][1],shorten[3*i1][j][2]); 257 truncatedcomplex[2*i1][size(truncatedcomplex[2*i1])+1]=fortrun[1]; 258 count=count+fortrun[2]; 259 fst=list(int(shorten[3*i1][j][1])1,int(shorten[3*i1][j][2])1); 260 sizetruncom[2*i1][size(sizetruncom[2*i1])+1]=fst; 261 sizetruncom[2*i][size(sizetruncom[2*i])+1]=count; 262 if (v!=0) 263 { 264 v[size(v)+1]=j; 265 } 266 else 267 { 268 v[1]=j; 269 } 270 } 271 } 272 if (v!=0) 273 { 274 subm=submat(truncatedcomplex[2*i],v,1..ncols(truncatedcomplex[2*i])); 275 truncatedcomplex[2*i]=subm; 276 if (i!=1) 277 { 278 i1=1..nrows(truncatedcomplex[2*(i1)]); 279 subm=submat(truncatedcomplex[2*(i1)],i1,v); 280 truncatedcomplex[2*(i1)]=subm; 281 } 282 } 283 else 284 { 285 truncatedcomplex[2*i]=matrix(0,1,ncols(truncatedcomplex[2*i])); 286 if (i!=1) 287 { 288 nr=nrows(truncatedcomplex[2*(i1)]); 289 truncatedcomplex[2*(i1)]=matrix(0,nr,1); 290 } 291 } 292 } 2216 { 2217 /*truncatedcomplex[2*i1] will contain all the generators for the truncation 2218 of D_n/(x(1),..,x(n))\otimes C[i]*/ 2219 truncatedcomplex[2*i1]=list(); 2220 sizetruncom[2*i1]=list(); 2221 sizetruncom[2*i]=list(); 2222 passendespoly[i]=list(); 2223 /*truncatedcomplex[2*i] will be the map trunc(D_n/(x(1),..,x(n))\otimes C[i]) 2224 >trunc(D_n/(x(1),..,x(n))\otimes C[i+1])*/ 2225 truncatedcomplex[2*i]=newcomplex[3*i]; 2226 v=0;count=0; 2227 sizetruncom[2*i][1]=0; 2228 for (j=1; j<=newcomplex[3*i2]; j++) 2229 { 2230 if (size(shorten[3*i1][j])!=0) 2231 { 2232 fortrun=sublist(allpolys,shorten[3*i1][j][1],shorten[3*i1][j][2]); 2233 truncatedcomplex[2*i1][size(truncatedcomplex[2*i1])+1]=fortrun[1]; 2234 for (k=1; k<=size(fortrun[1]); k++) 2235 { 2236 for (l=1; l<=size(fortrun[1][k]); l++) 2237 { 2238 passendespoly[i][size(passendespoly[i])+1]=list(fortrun[1][k][l][1],j); 2239 } 2240 } 2241 count=count+fortrun[2]; 2242 fst=list(int(shorten[3*i1][j][1])1,int(shorten[3*i1][j][2])1); 2243 sizetruncom[2*i1][size(sizetruncom[2*i1])+1]=fst; 2244 sizetruncom[2*i][size(sizetruncom[2*i])+1]=count; 2245 if (v!=0) 2246 { 2247 v[size(v)+1]=j; 2248 } 2249 else 2250 { 2251 v[1]=j; 2252 } 2253 } 2254 } 2255 if (v!=0) 2256 { 2257 keepv[i]=v; 2258 subm=submat(truncatedcomplex[2*i],v,1..ncols(truncatedcomplex[2*i])); 2259 truncatedcomplex[2*i]=subm; 2260 if (i!=1) 2261 { 2262 i1=1..nrows(truncatedcomplex[2*(i1)]); 2263 subm=submat(truncatedcomplex[2*(i1)],i1,v); 2264 truncatedcomplex[2*(i1)]=subm; 2265 } 2266 } 2267 else 2268 { 2269 keepv[i]=list(); 2270 truncatedcomplex[2*i]=matrix(0,1,ncols(truncatedcomplex[2*i])); 2271 if (i!=1) 2272 { 2273 nr=nrows(truncatedcomplex[2*(i1)]); 2274 truncatedcomplex[2*(i1)]=matrix(0,nr,1); 2275 } 2276 } 2277 } 2278 list keeptruncatedcomplex=truncatedcomplex; 293 2279 matrix M; 294 2280 int st,pi,pj; 295 2281 poly ptc; 296 2282 int b,d,ideg,kplus,lplus; 2283 int z; 297 2284 poly form,lform,nform; 298 2285 /*computation of the maps*/ 2286 if (diffforms==1) 2287 { 2288 def ConvWeyl=makeConverseWeyl(nvars(basering) div 2); 2289 setring ConvWeyl; 2290 poly form,lform,nform; 2291 poly ptc; 2292 list truncatedcomplex; 2293 matrix M; 2294 ideal I=x(1); 2295 for (i=2; i<=nvars(basering) div 2; i++) 2296 { 2297 I=I,var(nvars(basering) div 2 + i); 2298 } 2299 for (i=1; i<=nvars(basering) div 2; i++) 2300 { 2301 I=I,var(i); 2302 } 2303 map transtc=W,I; 2304 truncatedcomplex=transtc(truncatedcomplex); 2305 } 299 2306 for (i=1; i<size(truncatedcomplex) div 2; i++) 300 { 301 nr=max(1,sizetruncom[2*i][size(sizetruncom[2*i])]); 302 nc=max(1,sizetruncom[2*i+2][size(sizetruncom[2*i+2])]); 303 M=matrix(0,nr,nc); 304 for (k=1; k<=size(truncatedcomplex[2*i1]);k++) 305 { 306 for (l=1; l<=size(truncatedcomplex[2*(i+1)1]); l++) 307 { 308 if (size(sizetruncom[2*i])!=1) 309 { 310 for (j=1; j<=size(truncatedcomplex[2*i1][k]); j++) 311 { 312 for (b=1; b<=size(truncatedcomplex[2*i1][k][j]); b++) 313 { 314 form=truncatedcomplex[2*i1][k][j][b][1]; 315 form=form*truncatedcomplex[2*i][k,l]; 316 while (form!=0) 317 { 318 lform=lead(form); 319 v=leadexp(lform); 320 v=v[1..n]; 321 if (v==(0:n)) 322 { 323 ideg=deg(lform)sizetruncom[2*(i+1)1][l][1]; 324 if (ideg>=0) 325 { 326 nr=ideg+1; 327 st=size(truncatedcomplex[2*(i+1)1][l][nr]); 328 for (d=1; d<=st;d++) 2307 { 2308 nr=max(1,sizetruncom[2*i][size(sizetruncom[2*i])]); 2309 nc=max(1,sizetruncom[2*i+2][size(sizetruncom[2*i+2])]); 2310 M=matrix(0,nr,nc); 2311 for (k=1; k<=size(truncatedcomplex[2*i1]);k++) 2312 { 2313 for (l=1; l<=size(truncatedcomplex[2*(i+1)1]); l++) 2314 { 2315 if (size(sizetruncom[2*i])!=1) 2316 { 2317 for (j=1; j<=size(truncatedcomplex[2*i1][k]); j++) 329 2318 { 330 nc=2*(i+1)1; 331 ptc=truncatedcomplex[nc][l][ideg+1][d][1]; 332 if (leadmonom(lform)==ptc) 333 { 334 nr=2*i1; 335 pi=truncatedcomplex[nr][k][j][b][2]; 336 pi=pi+sizetruncom[2*i][k]; 337 nc=2*(i+1)1; 338 nr=ideg+1; 339 pj=truncatedcomplex[nc][l][nr][d][2]; 340 pj=pj+sizetruncom[2*(i+1)][l]; 341 M[pi,pj]=leadcoef(lform); 342 break; 343 } 2319 for (b=1; b<=size(truncatedcomplex[2*i1][k][j]); b++) 2320 { 2321 form=truncatedcomplex[2*i1][k][j][b][1]; 2322 form=form*truncatedcomplex[2*i][k,l]; 2323 2324 2325 for (z=1; z<=nvars(basering) div 2; z++)//neu 2326 {// 2327 form=subst(form,var(z),0);// 2328 }// 2329 2330 while (form!=0) 2331 { 2332 lform=lead(form); 2333 v=leadexp(lform); 2334 v=v[1..n]; 2335 // if (v==(0:n)) 2336 //{ 2337 ideg=deg(lform)sizetruncom[2*(i+1)1][l][1]; 2338 if (ideg>=0) 2339 { 2340 nr=ideg+1; 2341 st=size(truncatedcomplex[2*(i+1)1][l][nr]); 2342 for (d=1; d<=st;d++) 2343 { 2344 nc=2*(i+1)1; 2345 ptc=truncatedcomplex[nc][l][ideg+1][d][1]; 2346 if (leadmonom(lform)==ptc) 2347 { 2348 nr=2*i1; 2349 pi=truncatedcomplex[nr][k][j][b][2]; 2350 pi=pi+sizetruncom[2*i][k]; 2351 nc=2*(i+1)1; 2352 nr=ideg+1; 2353 pj=truncatedcomplex[nc][l][nr][d][2]; 2354 pj=pj+sizetruncom[2*(i+1)][l]; 2355 M[pi,pj]=leadcoef(lform); 2356 break; 2357 } 2358 } 2359 } 2360 // } 2361 2362 form=formlform; 2363 } 2364 } 344 2365 } 345 } 346 } 347 form=formlform; 348 } 349 } 350 } 351 } 352 } 353 } 354 truncatedcomplex[2*i]=M; 355 truncatedcomplex[2*i1]=sizetruncom[2*i][size(sizetruncom[2*i])]; 356 } 2366 } 2367 } 2368 } 2369 truncatedcomplex[2*i]=M; 2370 truncatedcomplex[2*i1]=sizetruncom[2*i][size(sizetruncom[2*i])]; 2371 } 357 2372 truncatedcomplex[2*i1]=sizetruncom[2*i][size(sizetruncom[2*i])]; 358 2373 if (truncatedcomplex[2*i1]!=0) 359 { 360 truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i1],1); 2374 { 2375 truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i1],1); 2376 } 2377 if (diffforms==1) 2378 { 2379 setring W; 2380 truncatedcomplex=imap(ConvWeyl,truncatedcomplex); 361 2381 } 362 2382 setring R; 363 364 2383 list truncatedcomplex=imap(W,truncatedcomplex); 2384 /*computes the cohomology of the complex (D^i,d^i) given by truncatedcomplex, 365 2385 i.e. D^i=C^truncatedcomplex[2*i1] and d^i=truncatedcomplex[2*i]*/ 366 list derhamhom=findCohomology(truncatedcomplex,le); 367 option(set,saveoptions); 368 return (derhamhom); 2386 if (diffforms==0) 2387 { 2388 list derhamhom=findCohomology(truncatedcomplex,le); 2389 option(set,saveoptions); 2390 return (derhamhom); 2391 } 2392 list outall=findCohomologyDiffForms(truncatedcomplex,le); 2393 setring W; 2394 list dimanddiff=imap(R,outall); 2395 list alldiffforms=dimanddiff[2]; 2396 while(size(alldiffforms)<size(passendespoly)) 2397 { 2398 passendespoly=delete(passendespoly,1); 2399 } 2400 list newdiffforms; 2401 matrix Diff; 2402 for (i=1; i<=size(alldiffforms); i++) 2403 { 2404 newdiffforms[i]=list(); 2405 for (j=1; j<=size(alldiffforms[i]); j++) 2406 { 2407 Diff=matrix(0,1,newcomplex[3*(i+size(newcomplex) div 3  size(alldiffforms))2]); 2408 for (k=1; k<=ncols(alldiffforms[i][j]); k++) 2409 { 2410 if (alldiffforms[i][j][1,k]!=0) 2411 { 2412 Diff[1,passendespoly[i][k][2]]=Diff[1,passendespoly[i][k][2]]+alldiffforms[i][j][1,k]*passendespoly[i][k][1]; 2413 } 2414 } 2415 newdiffforms[i][j]=Diff; 2416 } 2417 } 2418 list omegacomplex=makeOmega(nvars(W) div 2); 2419 list newcomplexmod; 2420 for (i=1; i<=size(newcomplex) div 3; i++) 2421 { 2422 newcomplexmod[2*i1]=newcomplex[3*i2]; 2423 newcomplexmod[2*i]=newcomplex[3*i]; 2424 } 2425 while (size(dimanddiff[1])<size(newcomplexmod) div 2) 2426 { 2427 newcomplexmod=delete(newcomplexmod,1); 2428 newcomplexmod=delete(newcomplexmod,1); 2429 } 2430 while (size(dimanddiff[1])<size(quasiiso)) 2431 { 2432 quasiiso=delete(quasiiso,1); 2433 } 2434 while (size(dimanddiff[1])>size(generators)) 2435 { 2436 generators=insert(generators,list()); 2437 } 2438 while (size(dimanddiff[1])>size(quasiiso)) 2439 { 2440 quasiiso=insert(quasiiso,list()); 2441 } 2442 int keepsign; 2443 list derhamdiff; 2444 list doublecom=makeDoubleComplex(newcomplexmod,omegacomplex,quasiiso,generators); 2445 matrix diffform; 2446 int stopping; 2447 int p; 2448 matrix convert; 2449 list interim; 2450 list correspondingposition; 2451 list allforms=list(); 2452 for (i=1; i<=size(newdiffforms); i++) 2453 { 2454 derhamdiff[i]=list(); 2455 allforms[i]=list(); 2456 for (j=1; j<=size(newdiffforms[i]); j++) 2457 { 2458 allforms[i][j]=list(); 2459 keepsign=1; 2460 derhamdiff[i][j]=0; 2461 diffform=newdiffforms[i][j];//Zeilenform 2462 correspondingposition=doublecom[i][1];//needed fpr transformation process 2463 interim=transferDiffforms(diffform,correspondingposition); 2464 if (size(interim)!=0) 2465 { 2466 allforms[i][j][size(allforms[i][j])+1]=interim; 2467 } 2468 stopping=0; 2469 p=1; 2470 for (k=i; k<=size(newdiffforms); k++) 2471 { 2472 keepsign=(1)*keepsign; 2473 if (stopping==0) 2474 { 2475 if (size(doublecom[k][p][2])==0) 2476 { 2477 stopping=1; 2478 } 2479 else 2480 { 2481 if (size(doublecom[k+1][p][3])!=0) 2482 { 2483 diffform=diffform*doublecom[k][p][2];//Spaltenform 2484 if (diffform!=matrix(0,nrows(diffform),ncols(diffform))) 2485 { 2486 diffform=findPreimage(doublecom[k+1][p][3],transpose(diffform));//Zeilenform 2487 correspondingposition=doublecom[k+1][p+1];//needed for transformation process 2488 interim=transferDiffforms(keepsign*diffform,correspondingposition); 2489 if (size(interim)!=0) 2490 { 2491 allforms[i][j][size(allforms[i][j])+1]=interim; 2492 } 2493 p=p+1; 2494 } 2495 else 2496 { 2497 stopping=1; 2498 } 2499 } 2500 else 2501 { 2502 stopping=1; 2503 } 2504 } 2505 } 2506 } 2507 } 2508 } 2509 setring R; 2510 list allforms=fetch(W,allforms); 2511 option(set,saveoptions); 2512 return (allforms); 369 2513 } 370 2514 … … 385 2529 L is a list of nonconstant polynomials 386 2530 RETURN: ring W: the nth Weyl algebra @* 387 W contains a list MV, which represents the MayerVietrois complex (C^i,d^i) of the 2531 W contains a list MV, which represents the MayerVietrois complex (C^i,d^i) of the 388 2532 polynomials contained in L as follows:@* 389 2533 the C^i are given by D_n^ncols(C[2*i1])/im(C[2*i1]) and the differentials … … 393 2537 { 394 2538 /* We follow algorithm 3.2.5 in [R],if #!=0 we use also Remark 3.2.6 in [R] for 395 an additional iterative localization*/2539 an additional iterative localization*/ 396 2540 def R=basering; 397 2541 int i; 398 2542 int iterative=1; 399 2543 if (size(#)!=0) 400 {401 iterative=#[1];402 }2544 { 2545 iterative=#[1]; 2546 } 403 2547 for (i=1; i<=size(L); i++) 404 {405 if (L[i]==0)406 {407 print("localization with respect to 0 not possible");408 return();409 }410 if (leadcoef(L[i])L[i]==0)411 {412 print("polynomials must be nonconstant");413 return();414 }415 }2548 { 2549 if (L[i]==0) 2550 { 2551 print("localization with respect to 0 not possible"); 2552 return(); 2553 } 2554 if (leadcoef(L[i])L[i]==0) 2555 { 2556 print("polynomials must be nonconstant"); 2557 return(); 2558 } 2559 } 416 2560 if (iterative==1) 417 {418 /*compute the localizations by factorizing the polynomials and iterative419 localization of the factors */420 for (i=1; i<=size(L); i++)421 {422 L[i]=factorize(L[i],1);423 }424 }2561 { 2562 /*compute the localizations by factorizing the polynomials and iterative 2563 localization of the factors */ 2564 for (i=1; i<=size(L); i++) 2565 { 2566 L[i]=factorize(L[i],1); 2567 } 2568 } 425 2569 int r=size(L); 426 2570 int n=nvars(basering); … … 428 2572 /*construct the ring Ws*/ 429 2573 def W=makeWeyl(n); 430 setring W; 2574 setring W; 431 2575 list man=ringlist(W); 432 2576 if (n==1) 433 {434 man[2][1]="x(1)";435 man[2][2]="D(1)";436 def Wi=ring(man);437 setring Wi;438 kill W;439 def W=Wi;440 setring W;441 list man=ringlist(W);442 }2577 { 2578 man[2][1]="x(1)"; 2579 man[2][2]="D(1)"; 2580 def Wi=ring(man); 2581 setring Wi; 2582 kill W; 2583 def W=Wi; 2584 setring W; 2585 list man=ringlist(W); 2586 } 443 2587 man[2][size(man[2])+1]="s";; 444 2588 man[3][3]=man[3][2]; … … 448 2592 matrix M[1][1]; 449 2593 man[6]=transpose(concat(transpose(concat(man[6],M)),M)); 450 def Ws=ring(man); 2594 def Ws=ring(man); 451 2595 setring Ws; 452 int j,k,l,c; 453 list L=fetch(R,L); 454 list Cech; 2596 int j,k,l,c; 2597 list L=fetch(R,L); 2598 list Cech; 455 2599 ideal J=var(1+n); 456 2600 for (i=2; i<=n; i++) 457 {458 J=J,var(i+n);459 }2601 { 2602 J=J,var(i+n); 2603 } 460 2604 Cech[1]=list(J); 461 2605 list Theta, remminroots; … … 468 2612 int rmr; 469 2613 if (iterative==0) 470 {/*computation of the modules of the MV complex*/ 471 for (i=1; i<=r; i++) 472 { 473 findminintroot=list(); 474 Cech[i+1]=list(); 475 Theta[i+1]=list(); 476 k1=1; 477 for (j=1; j<=i; j++) 478 { 479 k1[size(k1)+1]=size(Theta[j+1]); 480 for (k=1; k<=k1[j]; k++) 481 { 482 Theta[j+1][size(Theta[j+1])+1]=list(Theta[j][k][1]+list(i)); 483 Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i]; 484 /*We compute the sparametric annihilator J(s) and the bfunction 485 of the polynomial L[i] and Cech[i][k] to localize the module 486 D_n/(D(1),...,D(n))[L[i]^(1)]\otimes D_n^c/im(Cech[i][k]), 487 where c=ncols(Cech[i][k]) and the im(Cech[i][k]) is generated by 488 the rows of the matrix. 489 If we plug the minimal integer root r(or a smaller integer 490 value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic to 491 the above localization*/ 492 rem=SannfsIBM(L[i],Cech[j][k]); 493 Cech[j+1][size(Cech[j+1])+1]=rem[1]; 494 findminintroot[size(findminintroot)+1]=rem[2]; 495 } 496 } 497 /* we compute the minimal root of all bfunctions of L[i] computed above, 498 because we want to plug in the same root r in all sparametric 499 annihilators we computed for L[i] >this will ensure we can compute 500 the maps of the MV complex*/ 501 minroot=minIntRoot0(findminintroot); 502 for (j=1; j<=i; j++) 503 { 504 for (k=1; k<=k1[j]; k++) 505 { 506 sk=size(Cech[j+1])+1k; 507 Cech[j+1][size(Cech[j+1])+1k]=subst(Cech[j+1][sk],s,minroot); 508 } 509 } 510 remminroots[i]=minroot; 511 } 512 Cech=delete(Cech,1); 513 Theta=delete(Theta,1); 514 list zw; 515 poly reme; 516 /*computation of the maps of the MV complex*/ 517 for (i=1; i<r; i++) 518 { 519 diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1])); 520 for (j=1; j<=size(Cech[i]); j++) 521 { 522 for (k=1; k<=size(Cech[i+1]); k++) 523 { 524 zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]); 525 if (zw[2]!=0) 526 { 527 rmr=remminroots[zw[1]]; 528 reme=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr); 529 zw[2]=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr); 530 diffmaps[i][j,k]=zw[2]; 531 } 532 } 533 } 534 } 535 diffmaps[r]=matrix(0,1,1); 536 } 2614 {/*computation of the modules of the MV complex*/ 2615 for (i=1; i<=r; i++) 2616 { 2617 findminintroot=list(); 2618 Cech[i+1]=list(); 2619 Theta[i+1]=list(); 2620 k1=1; 2621 for (j=1; j<=i; j++) 2622 { 2623 k1[size(k1)+1]=size(Theta[j+1]); 2624 for (k=1; k<=k1[j]; k++) 2625 { 2626 Theta[j+1][size(Theta[j+1])+1]=list(Theta[j][k][1]+list(i)); 2627 Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i]; 2628 /*We compute the sparametric annihilator J(s) and the bfunction 2629 of the polynomial L[i] and Cech[i][k] to localize the module 2630 D_n/(D(1),...,D(n))[L[i]^(1)]\otimes D_n^c/im(Cech[i][k]), 2631 where c=ncols(Cech[i][k]) and the im(Cech[i][k]) is generated by 2632 the rows of the matrix. 2633 If we plug the minimal integer root r(or a smaller integer 2634 value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic to 2635 the above localization*/ 2636 rem=SannfsIBM(L[i],Cech[j][k]); 2637 Cech[j+1][size(Cech[j+1])+1]=rem[1]; 2638 findminintroot[size(findminintroot)+1]=rem[2]; 2639 } 2640 } 2641 /* we compute the minimal root of all bfunctions of L[i] computed above, 2642 because we want to plug in the same root r in all sparametric 2643 annihilators we computed for L[i] >this will ensure we can compute 2644 the maps of the MV complex*/ 2645 minroot=minIntRoot(findminintroot); 2646 for (j=1; j<=i; j++) 2647 { 2648 for (k=1; k<=k1[j]; k++) 2649 { 2650 sk=size(Cech[j+1])+1k; 2651 Cech[j+1][size(Cech[j+1])+1k]=subst(Cech[j+1][sk],s,minroot); 2652 } 2653 } 2654 remminroots[i]=minroot; 2655 } 2656 Cech=delete(Cech,1); 2657 Theta=delete(Theta,1); 2658 list zw; 2659 poly reme; 2660 /*computation of the maps of the MV complex*/ 2661 for (i=1; i<r; i++) 2662 { 2663 diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1])); 2664 for (j=1; j<=size(Cech[i]); j++) 2665 { 2666 for (k=1; k<=size(Cech[i+1]); k++) 2667 { 2668 zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]); 2669 if (zw[2]!=0) 2670 { 2671 rmr=remminroots[zw[1]]; 2672 reme=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr); 2673 zw[2]=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr); 2674 diffmaps[i][j,k]=zw[2]; 2675 } 2676 } 2677 } 2678 } 2679 diffmaps[r]=matrix(0,1,1); 2680 } 2681 list generators; 537 2682 if (iterative==1) 538 { 539 for (i=1; i<=r;i++) 540 { 541 Cech[i+1]=list(); 542 Theta[i+1]=list(); 543 k1=1; 544 for (c=1; c<=size(L[i]); c++) 545 { 546 findminintroot=list(); 547 for (j=1; j<=i; j++) 548 { 549 if (c==1) 550 { 551 k1[size(k1)+1]=size(Theta[j+1]); 552 } 553 for (k=1; k<=k1[j]; k++) 554 { 555 /*We compute the sparametric annihilator J(s) und the b 556 function of the polynomial L[i][c] and Cech[i][k] to 557 localize the module D_n/(D(1),...,D(n))[L[i][c]^(1)]\otimes 558 D_n^c/im(Cech[i][k]), where c=ncols(Cech[i][k]). 559 If we plug the minimal integer root r(or a smaller integer 560 value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic 561 to the above localization*/ 562 if (c==1) 563 { 564 rmr=size(Theta[j+1])+1; 565 Theta[j+1][rmr]=list(Theta[j][k][1]+list(i)); 566 Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i][c]; 567 rem=SannfsIBM(L[i][c],Cech[j][k]); 568 Cech[j+1][size(Cech[j+1])+1]=rem[1]; 569 findminintroot[size(findminintroot)+1]=rem[2]; 570 } 571 else 572 { 573 st=size(Theta[j+1])k1[j]+k; 574 Theta[j+1][st][2]=Theta[j+1][st][2]*L[i][c]; 575 rem=SannfsIBM(L[i][c],Cech[j+1][size(Cech[j+1])k1[j]+k]); 576 Cech[j+1][size(Cech[j+1])k1[j]+k]=rem[1]; 577 findminintroot[size(findminintroot)+1]=rem[2]; 578 } 579 } 580 } 581 /* we compute the minimal root of all bfunctions of L[i][c] 582 computed above,because we want to plug in the same root r in all 583 sparametric annihilators we computed for L[i] >this will 584 ensure we can compute the maps of the MV complex*/ 585 minroot=minIntRoot0(findminintroot); 586 for (j=1; j<=i; j++) 587 { 588 for (k=1; k<=k1[j]; k++) 589 { 590 st=size(Cech[j+1])+1k; 591 Cech[j+1][st]=subst(Cech[j+1][st],s,minroot); 592 } 593 } 594 if (c==1) 595 { 596 remminroots[i]=list(); 597 } 598 remminroots[i][c]=minroot; 599 } 600 } 601 Cech=delete(Cech,1); 602 Theta=delete(Theta,1); 603 list zw; 604 poly reme; 605 /*maps of the MV Complex*/ 606 for (i=1; i<r; i++) 607 { 608 diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1])); 609 for (j=1; j<=size(Cech[i]); j++) 610 { 611 for (k=1; k<=size(Cech[i+1]); k++) 612 { 613 zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]); 614 if (zw[2]!=0) 615 { 616 reme=1; 617 for (c=1; c<=size(L[zw[1]]);c++) 618 { 619 reme=reme*L[zw[1]][c]^(remminroots[zw[1]][c]); 620 } 621 diffmaps[i][j,k]=zw[2]*reme; 622 } 623 } 624 } 625 } 626 diffmaps[r]=matrix(0,1,1); 627 } 2683 { 2684 for (i=1; i<=r;i++) 2685 { 2686 generators[i]=list();//////////////////////////////////////////////////////////////////// 2687 Cech[i+1]=list(); 2688 Theta[i+1]=list(); 2689 k1=1; 2690 for (c=1; c<=size(L[i]); c++) 2691 { 2692 findminintroot=list(); 2693 for (j=1; j<=i; j++) 2694 { 2695 if (c==1) 2696 { 2697 k1[size(k1)+1]=size(Theta[j+1]); 2698 } 2699 for (k=1; k<=k1[j]; k++) 2700 { 2701 /*We compute the sparametric annihilator J(s) und the b 2702 function of the polynomial L[i][c] and Cech[i][k] to 2703 localize the module D_n/(D(1),...,D(n))[L[i][c]^(1)]\otimes 2704 D_n^c/im(Cech[i][k]), where c=ncols(Cech[i][k]). 2705 If we plug the minimal integer root r(or a smaller integer 2706 value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic 2707 to the above localization*/ 2708 if (c==1) 2709 { 2710 rmr=size(Theta[j+1])+1; 2711 Theta[j+1][rmr]=list(Theta[j][k][1]+list(i)); 2712 Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i][c]; 2713 rem=SannfsIBM(L[i][c],Cech[j][k]); 2714 Cech[j+1][size(Cech[j+1])+1]=rem[1]; 2715 findminintroot[size(findminintroot)+1]=rem[2]; 2716 } 2717 else 2718 { 2719 st=size(Theta[j+1])k1[j]+k; 2720 Theta[j+1][st][2]=Theta[j+1][st][2]*L[i][c]; 2721 rem=SannfsIBM(L[i][c],Cech[j+1][size(Cech[j+1])k1[j]+k]); 2722 Cech[j+1][size(Cech[j+1])k1[j]+k]=rem[1]; 2723 findminintroot[size(findminintroot)+1]=rem[2]; 2724 } 2725 } 2726 } 2727 /* we compute the minimal root of all bfunctions of L[i][c] 2728 computed above,because we want to plug in the same root r in all 2729 sparametric annihilators we computed for L[i] >this will 2730 ensure we can compute the maps of the MV complex*/ 2731 minroot=minIntRoot(findminintroot); 2732 for (j=1; j<=i; j++) 2733 { 2734 for (k=1; k<=k1[j]; k++) 2735 { 2736 st=size(Cech[j+1])+1k; 2737 Cech[j+1][st]=subst(Cech[j+1][st],s,minroot); 2738 } 2739 } 2740 if (c==1) 2741 { 2742 remminroots[i]=list(); 2743 } 2744 remminroots[i][c]=minroot; 2745 } 2746 } 2747 Cech=delete(Cech,1); 2748 Theta=delete(Theta,1); 2749 list zw; 2750 poly reme; 2751 /*maps of the MV Complex*/ 2752 for (i=1; i<r; i++) 2753 { 2754 diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1])); 2755 for (j=1; j<=size(Cech[i]); j++) 2756 { 2757 for (k=1; k<=size(Cech[i+1]); k++) 2758 { 2759 zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]); 2760 if (zw[2]!=0) 2761 { 2762 reme=1; 2763 for (c=1; c<=size(L[zw[1]]);c++) 2764 { 2765 reme=reme*L[zw[1]][c]^(remminroots[zw[1]][c]); 2766 } 2767 diffmaps[i][j,k]=zw[2]*reme; 2768 } 2769 } 2770 } 2771 } 2772 diffmaps[r]=matrix(0,1,1); 2773 for (i=1; i<=r; i++) 2774 { 2775 for (j=1; j<=size(Theta[i]); j++) 2776 { 2777 generators[i][j]=1; 2778 for (c=1; c<=size(Theta[i][j][1]); c++) 2779 { 2780 for (k=1; k<=size(L[Theta[i][j][1][c]]); k++) 2781 { 2782 generators[i][j]=generators[i][j]*L[Theta[i][j][1][c]][k]^((1)*remminroots[Theta[i][j][1][c]][k]); 2783 } 2784 } 2785 } 2786 } 2787 } 628 2788 setring W; 629 2789 /*map the modules and maps to the Weyl algebra*/ 630 2790 list diffmaps=imap(Ws,diffmaps); 631 2791 list Cechmodules=imap(Ws,Cech); 2792 if (iterative==1) 2793 { 2794 list Theta=imap(Ws,Theta); 2795 list generators=imap(Ws,generators); 2796 } 632 2797 list Cech; 633 2798 matrix sup; 634 2799 for (i=1; i<=r; i++) 635 {636 sup=transpose(matrix(Cechmodules[i][1]));637 Cech[2*i1]=sup;638 for (j=2; j<=size(Cechmodules[i]); j++)639 {640 sup=transpose(matrix(Cechmodules[i][j]));641 Cech[2*i1]=dsum(Cech[2*i1],sup);642 }643 sup=matrix(diffmaps[i]);644 Cech[2*i]=sup;645 }2800 { 2801 sup=transpose(matrix(Cechmodules[i][1])); 2802 Cech[2*i1]=sup; 2803 for (j=2; j<=size(Cechmodules[i]); j++) 2804 { 2805 sup=transpose(matrix(Cechmodules[i][j])); 2806 Cech[2*i1]=dsum(Cech[2*i1],sup); 2807 } 2808 sup=matrix(diffmaps[i]); 2809 Cech[2*i]=sup; 2810 } 646 2811 list MV=Cech; 2812 if (iterative==1) 2813 { 2814 export Theta; 2815 export generators; 2816 } 647 2817 export MV; 2818 648 2819 return (W); 649 2820 } 650 2821 651 2822 example 652 { "EXAMPLE:"; 2823 { "EXAMPLE:"; 653 2824 ring r = 0,(x,y,z),dp; 654 2825 list L=xy,xz; … … 662 2833 static proc SannfsIBM(poly F,ideal myJ) 663 2834 "USAGE: SannfsIBM(f,J), F poly, J ideal 664 ASSUME: basering is D_n[s], where D_n is the Weyl algebra and s and extra 2835 ASSUME: basering is D_n[s], where D_n is the Weyl algebra and s and extra 665 2836 commutative variable@* 666 f is a polynomial in the variables x(1),...,x(n) with rational coefficients 2837 f is a polynomial in the variables x(1),...,x(n) with rational coefficients 667 2838 @* 668 J is holonomic and fsaturated 2839 J is holonomic and fsaturated 669 2840 RETURN AlList of the form (K,g), where K is an ideal and g a univariant polynomial 670 2841 in the variable s. K is the sparametric annihilator of F and J and g is … … 673 2844 { 674 2845 /*modified version of the procedure SannfsBM from the library dmod.lib: SannfsBM 675 computes the sparametric annihilator for J=(x_1,...,x_n)*/676 /* We use Algorithm 3.1.12 in[R] to compute the sparametric 677 678 2846 computes the sparametric annihilator for J=(x_1,...,x_n)*/ 2847 /* We use Algorithm 3.1.12 in[R] to compute the sparametric 2848 annihilator. Then we use the sparametric annihilator to compute the bfunction 2849 via Algorithm 4.7 in [W1].*/ 679 2850 /* We assume that the basering the the nth Weyl algebra D_n. We create the ring 680 2851 D_n[s,t], where t*s=s*tt*/ 681 2852 def save = basering; 682 2853 int N = nvars(basering)1; … … 688 2859 list tmp; 689 2860 intvec iv; 690 L[1] = RL[1]; 691 L[4] = RL[4]; 2861 L[1] = RL[1]; 2862 L[4] = RL[4]; 692 2863 list Name = RL[2]; 693 2864 Name=delete(Name,size(Name)); … … 696 2867 RName[2] = "s"; 697 2868 list DName; 698 2869 for(i=1;i<=N div 2;i++) 699 2870 { 700 DName[i] = var(N div 2+i); 2871 DName[i] = var(N div 2+i); 701 2872 Name=delete(Name,N div 2+1); 702 2873 } … … 706 2877 L[2] = NName; 707 2878 kill NName; 708 tmp[1] = "lp"; 2879 tmp[1] = "lp"; 709 2880 iv = 1,1; 710 tmp[2] = iv; 2881 tmp[2] = iv; 711 2882 Lord[1] = tmp; 712 tmp[1] = "dp"; 2883 tmp[1] = "dp"; 713 2884 s = "iv="; 714 2885 for(i=1;i<=Nnew;i++) … … 742 2913 ideal myJ=imap(save,myJ); 743 2914 for (i=1; i<=N div 2; i++) 744 {745 myJ=subst(myJ,D(i),D(i)+diff(F,x(i))*t);746 }2915 { 2916 myJ=subst(myJ,D(i),D(i)+diff(F,x(i))*t); 2917 } 747 2918 ideal I = t*F+s; 748 2919 I=I,myJ;//the sparametric annihilator in D_n[s,t] … … 763 2934 764 2935 //////////////////////////////////////////////////////////////////////////////////// 765 //COMPUTATION OF QAUASIISOMORPHIC V_DSTRICT FREE COMPLEX2936 //COMPUTATION OF A QUASIISOMORPHIC V_DSTRICT FREE COMPLEX 766 2937 //////////////////////////////////////////////////////////////////////////////////// 767 2938 2939 static proc quasiisomorphicVdComplex(list L,list #) 2940 "USAGE: quasiisomorphicVdComplex(L[,df]); L a list of the form (M_1,f_1,...,M_s,f_s), 2941 where the M_i and f_i are matrices 2942 ASSUME: Basering is the Weyl algebra D_n @* 2943 (M_1,f_1,...,M_s,f_s) represents a complex 0>D_n^(r_1)/im(M_1)> 2944 D_n^(r_2)/im(M_2)>...>D_n^(r_s)>0 with differentials f_i, where im(M_i) 2945 is generated by the rows of M_i. In particular it hold:@* 2946  The M_i are m_i x r_imatrices and the f_iare r_i x r_(i+1)matrices @* 2947 the image of M_1*f_i is contained in the image of M_(i+1) @* 2948 d is an integer between 1 and n. If no value for d is given, it is assumed 2949 to be n @* 2950 df is an optional int, if df equals 1 a \tilde(V_d)strict complex 2951 will be computed (instead of a V_dstrict one) (for a definition see [W3]) 2952 RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @* 2953 L_1 is of the form (i_(n1),g_(n1),m_(n1),...,i_s,g_s,m_s) such that:@* 2954 the i_j are integers, the g_j are i_j x i_(j+1)matrices, the m_j intvecs 2955 of size i_j@* 2956 D_n^(i_(n1))[m_(n1)]>...>D_n^(i_s)[m_s]>0 is a V_dstrict complex 2957 with differentials m_i that is quasiisomorphic to the complex given by L@* 2958 L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and 2959 the n_i are shift vectors such that:@* 2960 coker(H_i) is the ith cohomology group of the complex given by L_1@* 2961 the n_i are the shift vectors of the coker(H_i) 2962 THEORY: We follow Proposition 3.2 and Corollary 3.3 in [W3] 2963 " 2964 { 2965 int tilde; 2966 if (size(#)!=0) 2967 { 2968 tilde=#[1]; 2969 } 2970 def B=basering; 2971 int n=nvars(B) div 2 + 1;//+1 müsste stimmen! bitte kontrollieren! 2972 int d=nvars(B) div 2; 2973 int r=size(L) div 2; 2974 int lonc=n+r; 2975 int Kiold=0; 2976 matrix kerold; 2977 // matrix kernew=out[r][2][2]; 2978 matrix kernew=diag(1,ncols(L[size(L)1])); 2979 module mL; 2980 int i; 2981 int k; 2982 matrix testm; 2983 int Kinew=nrows(kernew); 2984 int Jiold=0; 2985 int Jinew=0; 2986 matrix Niold; 2987 matrix Ninew; 2988 list newcomplex; 2989 int Aiold=Kinew; 2990 matrix savediv; 2991 newcomplex[3*lonc2]=Kinew; 2992 newcomplex[3*lonc1]=intvec(0:Kinew); 2993 newcomplex[3*lonc]=matrix(0,Kinew,1); 2994 list quasiiso; 2995 quasiiso[lonc]=diag(1,Kinew); 2996 matrix invimage; 2997 matrix keralpha; 2998 intvec v; 2999 int j; 3000 matrix sc; 3001 matrix fnc; 3002 int indk; 3003 int indj; 3004 int Aiold; 3005 list saveres; 3006 matrix Liplus; 3007 for (i=r1; i>=0; i) 3008 { 3009 indk=0; 3010 indj=0; 3011 Kiold=Kinew; 3012 kerold=kernew; 3013 if (i!=0) 3014 { 3015 // kernew=divdr(L[2*i],L[2*i+1],1); 3016 kernew=divdr(L[2*i],L[2*i+1]); 3017 mL=slimgb(transpose(L[2*i1])); 3018 for (k=1; k<=nrows(kernew); k++) 3019 { 3020 testm=reduce(transpose(submat(kernew,k,intvec(1..ncols(kernew)))),mL); 3021 if (testm==matrix(0,nrows(testm),ncols(testm))) 3022 { 3023 kernew=transpose(deletecol(transpose(kernew),k)); 3024 k=k1; 3025 } 3026 } 3027 Kinew=nrows(kernew); 3028 if (kernew==matrix(0,nrows(kernew),ncols(kernew))) 3029 { 3030 Kinew=0; 3031 indk=1; 3032 } 3033 } 3034 else 3035 { 3036 Kinew=0; 3037 indk=1; 3038 } 3039 Jiold=Jinew; 3040 Niold=Ninew; 3041 keralpha=transpose(syz(transpose(newcomplex[3*(i+n)+3]))); 3042 if (i!=0) 3043 { 3044 invimage=divdr(quasiiso[n+i+1],transpose(concat(transpose(L[2*i]),transpose(L[2*i+1])))); 3045 Ninew=vdStrictIntersect(keralpha,invimage,newcomplex[3*(n+i+1)1],tilde);////////////// 3046 } 3047 else 3048 { 3049 invimage=divdr(quasiiso[n+i+1],L[2*i+1]); 3050 saveres=vdStrictIntersectPlus(keralpha,invimage,newcomplex[3*(n+i+1)1],tilde);//////////////////////// 3051 3052 ///////////////////BIS HIERHIN VERALLGEMEINERT//////////////////////////////////////////////////////////////////// 3053 3054 3055 Ninew=saveres[1]; 3056 } 3057 Jinew=nrows(Ninew); 3058 if (Ninew==matrix(0,nrows(Ninew),ncols(Ninew))) 3059 { 3060 Jinew=0; 3061 indk=1; 3062 } 3063 newcomplex[3*(n+i)2]=Kinew+Jinew; 3064 v=0; 3065 if (indk==0) 3066 { 3067 v=(0:Kinew); 3068 if (indj==0) 3069 { 3070 fnc=transpose(concat(transpose(matrix(0,Kinew,Kiold+Jiold)),transpose(Ninew))); 3071 } 3072 else 3073 { 3074 fnc=matrix(0,Kinew,Kiold+Jiold); 3075 } 3076 } 3077 else 3078 { 3079 if (indj==0) 3080 { 3081 fnc=Ninew; 3082 } 3083 else 3084 { 3085 fnc=matrix(0,1,Kiold+Jiold); 3086 newcomplex[3*(n+i)2]=1; 3087 } 3088 } 3089 Aiold=Jinew+Kinew; 3090 if (Aiold==0) 3091 { 3092 Aiold=1; 3093 } 3094 newcomplex[3*(n+i)]=fnc; 3095 for (j=1; j<=Jinew; j++) 3096 { 3097 if (tilde==0) 3098 { 3099 v[Kinew+j]=VdDeg(submat(Ninew,j,(1..ncols(Ninew))),nvars(B) div 2,newcomplex[3*(n+i)+2]); 3100 } 3101 else 3102 { 3103 v[Kinew+j]=VdDegTilde(submat(Ninew,j,(1..ncols(Ninew))),nvars(B) div 2,newcomplex[3*(n+i)+2]); 3104 } 3105 } 3106 newcomplex[3*(n+i)1]=v; 3107 if (i==0) 3108 { 3109 quasiiso[n+i]=matrix(0,Jinew,1); 3110 } 3111 else 3112 { 3113 if (indj==0) 3114 { 3115 sc=submat(fnc,intvec(Kinew+1..nrows(fnc)),intvec(1..ncols(fnc)))*quasiiso[n+i+1]; 3116 Liplus=transpose(concat(transpose(L[2*i]),transpose(L[2*i+1]))); 3117 sc=matrixLift(Liplus,sc);//stimmt das jetzt 3118 sc=submat(sc,intvec(1..nrows(sc)),intvec(1..nrows(L[2*i]))); 3119 if (indk==0) 3120 { 3121 //pi=kernew 3122 quasiiso[n+i]=transpose(concat(transpose(kernew),transpose(sc))); 3123 } 3124 else 3125 { 3126 quasiiso[n+i]=sc; 3127 } 3128 } 3129 else 3130 { 3131 if (indk==0) 3132 { 3133 quasiiso[n+i]=kernew; 3134 } 3135 else 3136 { 3137 quasiiso[n+i]=matrix(0,1,ncols(kernew)); 3138 } 3139 } 3140 } 3141 } 3142 for (i=1; i<=n1; i++) 3143 { 3144 quasiiso[ni]=list(); 3145 if (size(saveres[2][i])!=0) 3146 { 3147 newcomplex[3*(ni)]=saveres[2][i]; 3148 newcomplex[3*(ni)2]=nrows(saveres[2][i]); 3149 v=0; 3150 for (j=1; j<=newcomplex[3*(ni)2]; j++) 3151 { 3152 if (tilde==0) 3153 { 3154 v[j]=VdDeg(submat(saveres[2][i],j,(1..ncols(saveres[2][i]))),nvars(B) div 2, newcomplex[3*(ni)+2]); 3155 } 3156 else 3157 { 3158 v[j]=VdDegTilde(submat(saveres[2][i],j,(1..ncols(saveres[2][i]))),nvars(B) div 2, newcomplex[3*(ni)+2]); 3159 } 3160 } 3161 newcomplex[3*(ni)1]=v; 3162 } 3163 else 3164 { 3165 newcomplex[3*(ni)]=matrix(0,1,1); 3166 if (newcomplex[3*(ni)+1]!=0) 3167 { 3168 newcomplex[3*(ni)]=matrix(0,1,newcomplex[3*(ni)+1]); 3169 } 3170 newcomplex[3*(ni)2]=int(0); 3171 newcomplex[3*(ni)1]=intvec(0); 3172 } 3173 } 3174 list result; 3175 result[1]=newcomplex; 3176 result[2]=list(); 3177 list forsep; 3178 for (i=1; i<=size(L) div 2+1; i++) 3179 { 3180 forsep[2*i]=newcomplex[3*(n+i1)]; 3181 forsep[2*i1]=matrix(0,1,nrows(forsep[2*i])); 3182 } 3183 forsep=shortExactPieces(forsep); 3184 list listofHis; 3185 matrix forVd; 3186 for (i=1; i<=size(L) div 2; i++) 3187 { 3188 v=0; 3189 listofHis[i]=list(forsep[i+1][1][5]); 3190 forVd=forsep[i+1][2][2]; 3191 for (j=1; j<=nrows(forVd); j++) 3192 { 3193 if (tilde==0) 3194 { 3195 v[j]=VdDeg(submat(forVd,j,intvec(1..ncols(forVd))),nvars(B) div 2, newcomplex[3*(n+i)1]); 3196 } 3197 else 3198 { 3199 v[j]=VdDegTilde(submat(forVd,j,intvec(1..ncols(forVd))),nvars(B) div 2, newcomplex[3*(n+i)1]); 3200 } 3201 } 3202 listofHis[i][2]=v; 3203 } 3204 result[2]=listofHis; 3205 result[3]=quasiiso; 3206 return(result); 3207 } 3208 3209 //////////////////////////////////////////////////////////////////////////////////// 3210 3211 static proc vdStrictIntersect(matrix M, matrix N, intvec v, int tilde) 3212 { 3213 def B=basering; 3214 option(returnSB);// alternative:erst intersect und dann SBBerechung mit slimgb 3215 if (tilde==0) 3216 { 3217 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2,v); 3218 } 3219 else 3220 { 3221 def HomWeyl=makeHomogenizedWeylTilde(nvars(B) div 2,v); 3222 } 3223 setring HomWeyl; 3224 matrix M=fetch(B,M); 3225 matrix N=fetch(B,N); 3226 M=nHomogenize(M); 3227 N=nHomogenize(N); 3228 matrix vdintersection=transpose(intersect(transpose(M),transpose(N))); 3229 vdintersection=subst(vdintersection,h,1); 3230 setring B; 3231 matrix vdintersection=fetch(HomWeyl,vdintersection); 3232 option(noreturnSB); 3233 return(vdintersection); 3234 } 3235 3236 //////////////////////////////////////////////////////////////////////////////////// 3237 3238 static proc vdStrictIntersectPlus(matrix M, matrix N, intvec v, int tilde) 3239 { 3240 def B=basering; 3241 int n=nvars(B) div 2; 3242 matrix vdint=transpose(intersect(transpose(M),transpose(N))); 3243 if (tilde==0) 3244 { 3245 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2,v); 3246 } 3247 else 3248 { 3249 def HomWeyl=makeHomogenizedWeylTilde(nvars(B) div 2,v); 3250 } 3251 setring HomWeyl; 3252 matrix vdint=fetch(B,vdint); 3253 matrix N=fetch(B,N); 3254 vdint=nHomogenize(vdint); 3255 intvec i1; 3256 intvec i2; 3257 int i; 3258 int nr; 3259 int nc; 3260 def ringofSyz=Sres(transpose(vdint),n);//////////////////////////////////////////////////////////////// 3261 setring ringofSyz; 3262 matrix vdint=transpose(matrix(RES[2])); 3263 vdint=subst(vdint,h,1); 3264 int logens=ncols(vdint)+1; 3265 int omitemptylist; 3266 matrix zerom; 3267 list rofA; 3268 for (i=3; i<=n+3; i++)////////////////////////////////////////////////////////////////////////////n und si müssen noch definiert werden 3269 { 3270 if (size(RES)>=i) 3271 { 3272 zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))); 3273 if (RES[i]!=zerom) 3274 { 3275 rofA[i2]=(matrix(RES[i])); 3276 if (i==3) 3277 { 3278 if (nrows(rofA[i2])logens+1!=nrows(vdint)) 3279 { 3280 //build the resolution 3281 nr=nrows(vdint)+logens1; 3282 nc=ncols(rofA[i2]); 3283 rofA[i2]=matrix(rofA[i2],nr,nc); 3284 } 3285 3286 } 3287 if (i!=3) 3288 { 3289 if (nrows(rofA[i2])logens+1!=nrows(rofA[i3])) 3290 { 3291 nr=nrows(rofA[i3])+logens1; 3292 nc=ncols(rofA[i2]); 3293 rofA[i2]=matrix(rofA[i2],nr,nc); 3294 } 3295 } 3296 i1=intvec(logens..nrows(rofA[i2])); 3297 i2=intvec(1..ncols(rofA[i2])); 3298 rofA[i2]=transpose(submat(rofA[i2],i1,i2)); 3299 logens=logens+ncols(rofA[i2]); 3300 rofA[i2]=subst(rofA[i2],h,1); 3301 } 3302 else 3303 { 3304 rofA[i2]=list(); 3305 } 3306 } 3307 else 3308 { 3309 rofA[i2]=list(); 3310 } 3311 } 3312 if(size(rofA[1])==0) 3313 { 3314 omitemptylist=1; 3315 } 3316 setring B; 3317 vdint=fetch(ringofSyz,vdint); 3318 if (omitemptylist!=1) 3319 { 3320 list rofA=fetch(ringofSyz,rofA); 3321 } 3322 kill HomWeyl; 3323 kill ringofSyz; 3324 return(list(vdint,rofA)); 3325 } 3326 3327 //////////////////////////////////////////////////////////////////////////////////// 3328 768 3329 static proc toVdStrictFreeComplex(list L,string Syzstring,list #) 769 "USAGE: toVdStrictFreeComplex(L, Syzstring [,d]); L a list of the form 770 (M_1,f_1,...,M_s,f_s), where the M_i and f_i are matrices, Syzstring a 3330 "USAGE: toVdStrictFreeComplex(L, Syzstring [,d]); L a list of the form 3331 (M_1,f_1,...,M_s,f_s), where the M_i and f_i are matrices, Syzstring a 771 3332 string, d an optional integer 772 3333 ASSUME: Basering is the Weyl algebra D_n @* 773 3334 (M_1,f_1,...,M_s,f_s) represents a complex 0>D_n^(r_1)/im(M_1)> 774 D_n^(r_2)/im(M_2)>...>D_n^(r_s)>0 with differentials f_i, where im(M_i) 3335 D_n^(r_2)/im(M_2)>...>D_n^(r_s)>0 with differentials f_i, where im(M_i) 775 3336 is generated by the rows of M_i. In particular it hold:@* 776 3337  The M_i are m_i x r_imatrices and the f_iare r_i x r_(i+1)matrices @* 777 3338 the image of M_1*f_i is contained in the image of M_(i+1) @* 778 d is an integer between 1 and n. If no value for d is given, it is assumed779 to be n@*3339 d is an optional integer which indices in the case size(L)=2, whether a 3340 V_dstrict or \tilde(V_d)strict will be computed@* 780 3341 Syzstring is either: @* 781 'Sres' (computes the resolutions and Groebner bases in the homogenized 782 Weyl algebra using Schreyer's method)@* 3342 'Sres' (computes the resolutions and Groebner bases in the homogenized 3343 Weyl algebra using Schreyer's method)@* 783 3344 or @* 784 'Vdres' (computes the resolutions via V_dhomogenization and without 3345 'Vdres' (computes the resolutions via V_dhomogenization and without 785 3346 Schreyer's method)@* 786 3347 RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @* 787 3348 L_1 is of the form (i_(n1),g_(n1),m_(n1),...,i_s,g_s,m_s) such that:@* 788 the i_j are integers, the g_j are i_j x i_(j+1)matrices, the m_j intvecs 3349 the i_j are integers, the g_j are i_j x i_(j+1)matrices, the m_j intvecs 789 3350 of size i_j@* 790 3351 D_n^(i_(n1))[m_(n1)]>...>D_n^(i_s)[m_s]>0 is a V_dstrict complex 791 3352 with differentials m_i that is quasiisomorphic to the complex given by L@* 792 L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and 3353 L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and 793 3354 the n_i are shift vectors such that:@* 794 coker(H_i) is the ith cohomology group of the complex given by L_1@* 3355 coker(H_i) is the ith cohomology group of the complex given by L_1@* 795 3356 the n_i are the shift vectors of the coker(H_i) 796 3357 THEORY: We follow Algorithm 3.8 in [W2] 797 3358 " 798 3359 { 799 def B=basering; 3360 def B=basering; 800 3361 int n=nvars(B) div 2+2; 801 3362 int d=nvars(B) div 2; … … 805 3366 matrix mem; 806 3367 intvec i1,i2; 3368 int tilde; 807 3369 if (size(#)!=0) 808 { 809 for (i=1; i<=size(#); i++) 810 { 811 if (typeof(#[i])=="int") 812 { 813 if (#[1]>=1 and #[1]<=n) 814 { 815 d=#[i]; 816 } 817 } 818 } 819 } 3370 { 3371 for (i=1; i<=size(#); i++) 3372 { 3373 if (typeof(#[i])=="int") 3374 { 3375 tilde=#[i]; 3376 } 3377 } 3378 } 820 3379 /* If size(L)=2, our complex consists for only one nontrivial module. 821 Therefore, we just have to compute a V_dstrict resolution of this module.*/ 822 if (size(L)==2) 823 { 824 v=(0:ncols(L[1])); 825 out[3*n1]=v; 826 out[3*n2]=ncols(L[1]); 827 out[3*n]=L[2]; 828 if (Syzstring=="Vdres") 829 { 830 /*if Syzstring="Vdres", we compute a V_dstrict Groebner basis of L[1] 831 using Fhomogenization (Prop. 3.9 in [OT]); then we compute the syzygies 832 and make them V_dstrict using Prop 3.9[OT] and so on*/ 833 out[3*n3]=VdStrictGB(L[1],d,v); 834 for (i=n1; i>=1; i) 835 { 836 out[3*i2]=nrows(out[3*i]); 837 v=0; 838 for (j=1; j<=out[3*i2]; j++) 839 { 840 mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); 841 v[j]=VdDeg(mem,d, out[3*i+2]);//next shift vector 842 } 843 out[3*i1]=v; 844 if (i!=1) 845 { 846 /*next step in the resolution*/ 847 out[3*i3]=transpose(syz(transpose(out[3*i]))); 848 if (out[3*i3]!=matrix(0,nrows(out[3*i3]),ncols(out[3*i3]))) 849 { 850 /*makes the resolution V_dstrict*/ 851 out[3*i3]=VdStrictGB(out[3*i3],d,out[3*i1]); 852 } 853 else 854 { 855 /*resolution is already computed*/ 856 out[3*i3]=matrix(0,1,ncols(out[3*i3])); 857 out[3*i4]=intvec(0); 858 out[3*i5]=int(0); 859 for (j=i2; j>=1; j) 860 { 861 out[3*j]=matrix(0,1,1); 862 out[3*j1]=intvec(0); 863 out[3*j2]=int(0); 864 } 865 break; 866 } 867 } 868 } 869 } 870 else 871 { 872 /*in the case Syzstring!="Vdres" we compute the resolution in the 873 homogenized Weyl algebra using Thm 9.10 in[OT]*/ 874 def HomWeyl=makeHomogenizedWeyl(d); 875 setring HomWeyl; 876 list L=fetch(B,L); 877 L[1]=nHomogenize(L[1]); 878 list out=fetch(B,out); 879 out[3*n3]=L[1]; 880 /*computes a ring with a list RES; RES is a V_dstrict resolution of 881 L[1]*/ 882 def ringofSyz=Sres(transpose(L[1]),d); 883 setring ringofSyz; 884 int logens=2; 885 matrix mem; 886 list out=fetch(HomWeyl,out); 887 out[3*n3]=transpose(matrix(RES[2])); 888 out[3*n3]=subst(out[3*n3],h,1); 889 for (i=n1; i>=1; i) 890 { 891 out[3*i2]=nrows(out[3*i]); 892 v=0; 893 for (j=1; j<=out[3*i2]; j++) 894 { 895 mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); 896 v[j]=VdDeg(mem,d, out[3*i+2]); 897 } 898 out[3*i1]=v;//shift vector such that the resolution RES is V_dstrict 899 if (i!=1) 900 { 901 indi=0; 902 if (size(RES)>=ni+2) 903 { 904 nr=nrows(matrix(RES[ni+2])); 905 mem=matrix(0,nr,ncols(matrix(RES[ni+2]))); 906 if (matrix(RES[ni+2])!=mem) 907 { 908 indi=1; 909 out[3*i3]=(matrix(RES[ni+2])); 910 if (nrows(out[3*i3])logens+1!=nrows(out[3*i])) 911 { 912 mem=out[3*i3]; 913 out[3*i3]=matrix(mem,nrows(mem)+logens1,ncols(mem)); 914 } 915 mem=out[3*i3]; 916 i1=intvec(logens..nrows(mem)); 917 mem=submat(mem,i1,intvec(1..ncols(mem))); 918 out[3*i3]=transpose(mem); 919 out[3*i3]=subst(out[3*i3],h,1); 920 logens=logens+ncols(out[3*i3]); 921 } 922 } 923 if(indi==0) 924 { 925 out[3*i3]=matrix(0,1,nrows(out[3*i])); 926 out[3*i4]=intvec(0); 927 out[3*i5]=int(0); 928 for (j=i2; j>=1; j) 929 { 930 out[3*j]=matrix(0,1,1); 931 out[3*j1]=intvec(0); 932 out[3*j2]=int(0); 933 } 934 break; 935 } 936 } 937 } 938 setring B; 939 out=fetch(ringofSyz,out);//contains the V_dstrict resolution 940 kill ringofSyz; 941 } 942 outall[1]=out; 943 outall[2]=list(list(out[3*n3],out[3*n1])); 944 return(outall); 945 } 946 /*case size(L)>2: We compute a quasiisomorphic free complex following Alg 3.8 in 947 [W2]*/ 948 /* We denote the complex given by L as (C^i,d^i). 949 We start by computing in the proc shortExaxtPieces representations for the 950 short exact sequences B^i>Z^i>H^i and Z^i>C^i>B^(i+1), where the B^i, Z^i 951 and H^i are coboundaries, cocycles and cohomology groups, respectively.*/ 952 out=shortExactPieces(L); 3380 Therefore, we just have to compute a V_dstrict resolution of this module.*/ 3381 if (size(L)==2) 3382 { 3383 v=(0:ncols(L[1])); 3384 out[3*n1]=v; 3385 out[3*n2]=ncols(L[1]); 3386 out[3*n]=L[2]; 3387 if (Syzstring=="Vdres") 3388 { 3389 /*if Syzstring="Vdres", we compute a V_dstrict Groebner basis of L[1] 3390 using Fhomogenization (Prop. 3.9 in [OT]); then we compute the syzygies 3391 and make them V_dstrict using Prop 3.9[OT] and so on*/ 3392 out[3*n3]=VdStrictGB(L[1],d,v); 3393 for (i=n1; i>=1; i) 3394 { 3395 out[3*i2]=nrows(out[3*i]); 3396 v=0; 3397 for (j=1; j<=out[3*i2]; j++) 3398 { 3399 mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); 3400 v[j]=VdDeg(mem,d, out[3*i+2]);//next shift vector 3401 } 3402 out[3*i1]=v; 3403 if (i!=1) 3404 { 3405 /*next step in the resolution*/ 3406 out[3*i3]=transpose(syz(transpose(out[3*i]))); 3407 if (out[3*i3]!=matrix(0,nrows(out[3*i3]),ncols(out[3*i3]))) 3408 { 3409 /*makes the resolution V_dstrict*/ 3410 out[3*i3]=VdStrictGB(out[3*i3],d,out[3*i1]); 3411 } 3412 else 3413 { 3414 /*resolution is already computed*/ 3415 out[3*i3]=matrix(0,1,ncols(out[3*i3])); 3416 out[3*i4]=intvec(0); 3417 out[3*i5]=int(0); 3418 for (j=i2; j>=1; j) 3419 { 3420 out[3*j]=matrix(0,1,1); 3421 out[3*j1]=intvec(0); 3422 out[3*j2]=int(0); 3423 } 3424 break; 3425 } 3426 } 3427 } 3428 } 3429 else 3430 { 3431 /*in the case Syzstring!="Vdres" we compute the resolution in the 3432 homogenized Weyl algebra using Thm 9.10 in[OT]*/ 3433 if (tilde==0) 3434 { 3435 def HomWeyl=makeHomogenizedWeyl(d); 3436 } 3437 else 3438 { 3439 def HomWeyl=makeHomogenizedWeylTilde(d); 3440 } 3441 setring HomWeyl; 3442 list L=fetch(B,L); 3443 L[1]=nHomogenize(L[1]); 3444 list out=fetch(B,out); 3445 out[3*n3]=L[1]; 3446 /*computes a ring with a list RES; RES is a V_dstrict resolution of 3447 L[1]*/ 3448 def ringofSyz=Sres(transpose(L[1]),d); 3449 setring ringofSyz; 3450 int logens=2; 3451 matrix mem; 3452 list out=fetch(HomWeyl,out); 3453 out[3*n3]=transpose(matrix(RES[2])); 3454 out[3*n3]=subst(out[3*n3],h,1); 3455 for (i=n1; i>=1; i) 3456 { 3457 out[3*i2]=nrows(out[3*i]); 3458 v=0; 3459 for (j=1; j<=out[3*i2]; j++) 3460 { 3461 mem=submat(out[3*i],j,intvec(1..ncols(out[3*i]))); 3462 if (tilde==0) 3463 { 3464 v[j]=VdDeg(mem,d, out[3*i+2]); 3465 } 3466 else 3467 { 3468 v[j]=VdDegTilde(mem,d, out[3*i+2]); 3469 } 3470 } 3471 out[3*i1]=v;//shift vector such that the resolution RES is V_dstrict 3472 if (i!=1) 3473 { 3474 indi=0; 3475 if (size(RES)>=ni+2) 3476 { 3477 nr=nrows(matrix(RES[ni+2])); 3478 mem=matrix(0,nr,ncols(matrix(RES[ni+2]))); 3479 if (matrix(RES[ni+2])!=mem) 3480 { 3481 indi=1; 3482 out[3*i3]=(matrix(RES[ni+2])); 3483 if (nrows(out[3*i3])logens+1!=nrows(out[3*i])) 3484 { 3485 mem=out[3*i3]; 3486 out[3*i3]=matrix(mem,nrows(mem)+logens1,ncols(mem)); 3487 } 3488 mem=out[3*i3]; 3489 i1=intvec(logens..nrows(mem)); 3490 mem=submat(mem,i1,intvec(1..ncols(mem))); 3491 out[3*i3]=transpose(mem); 3492 out[3*i3]=subst(out[3*i3],h,1); 3493 logens=logens+ncols(out[3*i3]); 3494 } 3495 } 3496 if(indi==0) 3497 { 3498 out[3*i3]=matrix(0,1,nrows(out[3*i])); 3499 out[3*i4]=intvec(0); 3500 out[3*i5]=int(0); 3501 for (j=i2; j>=1; j) 3502 { 3503 out[3*j]=matrix(0,1,1); 3504 out[3*j1]=intvec(0); 3505 out[3*j2]=int(0); 3506 } 3507 break; 3508 } 3509 } 3510 } 3511 setring B; 3512 out=fetch(ringofSyz,out);//contains the V_dstrict resolution 3513 kill ringofSyz; 3514 } 3515 outall[1]=out; 3516 outall[2]=list(list(out[3*n3],out[3*n1])); 3517 return(outall); 3518 } 3519 /*case size(L)>2: We compute a quasiisomorphic free complex following Alg 3.8 in 3520 [W2]*/ 3521 /* We denote the complex given by L as (C^i,d^i). 3522 We start by computing in the proc shortExaxtPieces representations for the 3523 short exact sequences B^i>Z^i>H^i and Z^i>C^i>B^(i+1), where the B^i, Z^i 3524 and H^i are coboundaries, cocycles and cohomology groups, respectively.*/ 3525 out=shortExactPieces(L); 953 3526 list rem; 954 /* shortExactpiecesToVdStrict makes the sequences B^i>Z^i>H^i and 955 Z^i>C^i>B^(i+1) V_dstrict*/3527 /* shortExactpiecesToVdStrict makes the sequences B^i>Z^i>H^i and 3528 Z^i>C^i>B^(i+1) V_dstrict*/ 956 3529 rem=shortExactPiecesToVdStrict(out,d,Syzstring); 957 /*VdStrictDoubleComplexes computes V_dstrict resolutions over the seqeunces from 958 proc shortExactPiecesToVdstrict*/3530 /*VdStrictDoubleComplexes computes V_dstrict resolutions over the seqeunces from 3531 proc shortExactPiecesToVdstrict*/ 959 3532 out=VdStrictDoubleComplexes(rem[1],d,Syzstring); 960 3533 for (i=1;i<=size(out); i++) 961 {962 rem[2][i][1]=out[i][1][5][1];963 rem[2][i][2]=out[i][1][8][1];964 }965 /* AssemblingDoubleComplexes puts the resolution of the C^i (from the sequences 966 Z^i>C^i>B^(i+1)) together to obtain a CartanEilenberg resolution of967 (C^i,d^i)*/3534 { 3535 rem[2][i][1]=out[i][1][5][1]; 3536 rem[2][i][2]=out[i][1][8][1]; 3537 } 3538 /* AssemblingDoubleComplexes puts the resolution of the C^i (from the sequences 3539 Z^i>C^i>B^(i+1)) together to obtain a CartanEilenberg resolution of 3540 (C^i,d^i)*/ 968 3541 out=assemblingDoubleComplexes(out); 969 3542 /*the proc totalComplex takes the total complex of the double complex from the 970 proc assemblingDoubleComplexes*/3543 proc assemblingDoubleComplexes*/ 971 3544 out=totalComplex(out); 972 3545 outall[1]=out; 973 outall[2]=rem[2];//contains the cohomology groups and their shift vectors 3546 outall[2]=rem[2];//contains the cohomology groups and their shift vectors 974 3547 return (outall); 975 3548 } … … 984 3557 int count; 985 3558 for (i=m; i<=n; i++) 986 {987 out[size(out)+1]=list();988 for (j=1; j<=size(L[i]); j++)989 {990 count=count+1;991 out[size(out)][j]=list(L[i][j],count);992 }993 }3559 { 3560 out[size(out)+1]=list(); 3561 for (j=1; j<=size(L[i]); j++) 3562 { 3563 count=count+1; 3564 out[size(out)][j]=list(L[i][j],count); 3565 } 3566 } 994 3567 list o=list(out,count); 995 3568 return(o); … … 998 3571 //////////////////////////////////////////////////////////////////////////////////// 999 3572 1000 static proc LMSubset(list L,list M )3573 static proc LMSubset(list L,list M, list #) 1001 3574 { 1002 3575 int i; 1003 3576 int j=1; 1004 list position=(M[size(M)],(1)^(size(L))); 3577 if (size(#)==0) 3578 { 3579 list position=(M[size(M)],(1)^(size(L))); 3580 } 3581 else 3582 { 3583 list position=(M[size(M)],1); 3584 } 1005 3585 for (i=1; i<=size(L); i++) 1006 { 1007 if (L[i]!=M[j]) 1008 { 1009 if (L[i]!=M[i+1] or j!=i) 1010 { 1011 return (L[i],0); 1012 } 1013 else 1014 { 1015 position=(M[i],(1)^(i1)); 1016 j=j+1; 1017 } 1018 } 1019 j=j+1; 1020 1021 } 3586 { 3587 if (L[i]!=M[j]) 3588 { 3589 if (L[i]!=M[i+1] or j!=i) 3590 { 3591 return (L[i],0); 3592 } 3593 else 3594 { 3595 if (size(#)==0) 3596 { 3597 position=(M[i],(1)^(i1)); 3598 } 3599 else 3600 { 3601 position=(M[i],(1)^(size(L)+1i)); 3602 } 3603 j=j+1; 3604 } 3605 } 3606 j=j+1; 3607 3608 } 1022 3609 return (position); 1023 3610 } … … 1028 3615 { 1029 3616 /*we follow Section 3.3 in [W2]*/ 1030 /* we assume that L=(M_1,f_1,...,M_s,f_s) defines the complex C=(C^i,d^i) 1031 as in the procedure toVdstrictcomplex*/3617 /* we assume that L=(M_1,f_1,...,M_s,f_s) defines the complex C=(C^i,d^i) 3618 as in the procedure toVdstrictcomplex*/ 1032 3619 matrix Bnew= divdr(L[2],L[3]); 1033 3620 matrix Bold=Bnew; … … 1038 3625 list sep; 1039 3626 /* the list sep will be of size s such that 1040 sep[i]=(sep[i][1],sep[i][2]) is a list of two lists1041 sep[i][1]=(B^i,f^(BZi),Z^i,f_^(ZHi),H^i) such that coker(B^i)>coker(Z^i)1042 >coker(H^i) represents the short exact seqeuence B^i(C)>Z^i(C)>H^i(C)1043 sep[i][2]=(Z^i,f^(ZCi),C^i,f^(CBi),B^(i+1)) such that coker(Z^i)>coker(C^i)>1044 coker(B^(i+1)) represents the short exact seqeuence Z^i(C)>C^i>B^(i+1)(C)*/3627 sep[i]=(sep[i][1],sep[i][2]) is a list of two lists 3628 sep[i][1]=(B^i,f^(BZi),Z^i,f_^(ZHi),H^i) such that coker(B^i)>coker(Z^i) 3629 >coker(H^i) represents the short exact seqeuence B^i(C)>Z^i(C)>H^i(C) 3630 sep[i][2]=(Z^i,f^(ZCi),C^i,f^(CBi),B^(i+1)) such that coker(Z^i)>coker(C^i)> 3631 coker(B^(i+1)) represents the short exact seqeuence Z^i(C)>C^i>B^(i+1)(C)*/ 1045 3632 sep[1]=list(bzh,zcb); 1046 3633 int i; 1047 3634 list out; 1048 3635 for (i=3; i<=size(L)2; i=i+2) 1049 {1050 /*the proc bzhzcb computes representations for the short exact seqeunces */1051 out=bzhzcb(Bold, L[i1] , L[i], L[i+1], L[i+2]);1052 sep[size(sep)+1]=out[1];1053 Bold=out[2];1054 }3636 { 3637 /*the proc bzhzcb computes representations for the short exact seqeunces */ 3638 out=bzhzcb(Bold, L[i1] , L[i], L[i+1], L[i+2]); 3639 sep[size(sep)+1]=out[1]; 3640 Bold=out[2]; 3641 } 1055 3642 bzh=(divdr(L[size(L)2], L[size(L)1]),L[size(L)2], L[size(L)1]); 1056 3643 bzh[4]=unitmat(ncols(L[size(L)1])); … … 1079 3666 static proc shortExactPiecesToVdStrict(list C,int d,list #) 1080 3667 {/* We transform the short exact pieces from procedure shortExactPieces to V_d 1081 strict short exact sequences. For this, we use Algorithm 3.11 and Lemma 4.2 in 1082 [W2].*/ 3668 strict short exact sequences. For this, we use Algorithm 3.11 and Lemma 4.2 in 3669 [W2].*/ 1083 3670 /* If we compute our Groebner bases in the homogenized Weyl algebra, we already 1084 compute some resolutions it omit additional Groebner basis computations later1085 on.*/3671 compute some resolutions it omit additional Groebner basis computations later 3672 on.*/ 1086 3673 int s =size(C);int i; int j; 1087 3674 string Syzstring="Sres"; 1088 intvec v=0:ncols(C[s][1][5]); 3675 intvec v=0:ncols(C[s][1][5]); 1089 3676 if (size(#)!=0) 1090 {1091 for (i=1; i<=size(#); i++)1092 {1093 if (typeof(#[i])=="string")1094 {1095 Syzstring=#[i];1096 }1097 if (typeof(#[i])=="intvec")1098 {1099 v=#[i];1100 }1101 }1102 }3677 { 3678 for (i=1; i<=size(#); i++) 3679 { 3680 if (typeof(#[i])=="string") 3681 { 3682 Syzstring=#[i]; 3683 } 3684 if (typeof(#[i])=="intvec") 3685 { 3686 v=#[i]; 3687 } 3688 } 3689 } 1103 3690 list out; 1104 3691 list forout; 1105 3692 if (Syzstring=="Vdres") 1106 {1107 out[s]=list(toVdStrictSequence(C[s][1],d,v, Syzstring,s));1108 }3693 { 3694 out[s]=list(toVdStrictSequence(C[s][1],d,v, Syzstring,s)); 3695 } 1109 3696 else 1110 {1111 forout=toVdStrictSequence(C[s][1],d,v, Syzstring,s);1112 list resolutionofA=forout[9];1113 list resolutionofC=forout[10];1114 forout=delete(forout,10);1115 forout=delete(forout,9);1116 out[s]=list(forout);1117 for (i=1; i<=size(resolutionofC); i++)1118 {1119 out[s][1][5][i+1]=resolutionofC[i];//save the resolutions1120 out[s][1][1][i+1]=resolutionofA[i];1121 }1122 }3697 { 3698 forout=toVdStrictSequence(C[s][1],d,v, Syzstring,s); 3699 list resolutionofA=forout[9]; 3700 list resolutionofC=forout[10]; 3701 forout=delete(forout,10); 3702 forout=delete(forout,9); 3703 out[s]=list(forout); 3704 for (i=1; i<=size(resolutionofC); i++) 3705 { 3706 out[s][1][5][i+1]=resolutionofC[i];//save the resolutions 3707 out[s][1][1][i+1]=resolutionofA[i]; 3708 } 3709 } 1123 3710 out[s][2]=list(list(out[s][1][3][1])); 1124 3711 out[s][2][2]=list(unitmat(ncols(out[s][1][3][1]))); … … 1132 3719 list resolutionofF; 1133 3720 for (i=s1; i>=2; i) 1134 {1135 C[i][2][5]=out[i+1][1][1][1];1136 forout=toVdStrictSequences(C[i],d,out[i+1][1][6][1],Syzstring,s);1137 if (Syzstring=="Sres")1138 {1139 resolutionofD=forout[3];//save the resolutions1140 resolutionofF=forout[4];1141 forout=delete(forout,4);1142 forout=delete(forout,3);1143 }1144 out[i]=forout;1145 if(Syzstring=="Sres")1146 {1147 for (j=2; j<=size(out[i+1][1][1]); j++)1148 {1149 out[i][2][5][j]=out[i+1][1][1][j];1150 }1151 for (j=1; j<=size(resolutionofD);j++)1152 {1153 out[i][1][1][j+1]=resolutionofD[j];1154 out[i][1][5][j+1]=resolutionofF[j];1155 }1156 }1157 }3721 { 3722 C[i][2][5]=out[i+1][1][1][1]; 3723 forout=toVdStrictSequences(C[i],d,out[i+1][1][6][1],Syzstring,s); 3724 if (Syzstring=="Sres") 3725 { 3726 resolutionofD=forout[3];//save the resolutions 3727 resolutionofF=forout[4]; 3728 forout=delete(forout,4); 3729 forout=delete(forout,3); 3730 } 3731 out[i]=forout; 3732 if(Syzstring=="Sres") 3733 { 3734 for (j=2; j<=size(out[i+1][1][1]); j++) 3735 { 3736 out[i][2][5][j]=out[i+1][1][1][j]; 3737 } 3738 for (j=1; j<=size(resolutionofD);j++) 3739 { 3740 out[i][1][1][j+1]=resolutionofD[j]; 3741 out[i][1][5][j+1]=resolutionofF[j]; 3742 } 3743 } 3744 } 1158 3745 out[1]=list(list());//initalize our list 1159 3746 C[1][2][5]=out[2][1][1][1]; 1160 3747 /*Compute the last V_dstrict seqeunce*/ 1161 3748 if (Syzstring=="Vdres") 1162 {1163 out[1][2]=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");1164 }3749 { 3750 out[1][2]=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv"); 3751 } 1165 3752 else 1166 {1167 forout=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");1168 out[1][2]=delete(forout,9);1169 list resolutionofA2=forout[9];1170 for (i=1; i<=size(out[2][1][1]); i++)1171 {1172 /*put the modules for the resolutions in the right spot*/1173 out[1][2][5][i]=out[2][1][1][i];1174 }1175 for (i=1; i<=size(resolutionofA2); i++)1176 {1177 out[1][2][1][i+1]=resolutionofA2[i];1178 }1179 }3753 { 3754 forout=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv"); 3755 out[1][2]=delete(forout,9); 3756 list resolutionofA2=forout[9]; 3757 for (i=1; i<=size(out[2][1][1]); i++) 3758 { 3759 /*put the modules for the resolutions in the right spot*/ 3760 out[1][2][5][i]=out[2][1][1][i]; 3761 } 3762 for (i=1; i<=size(resolutionofA2); i++) 3763 { 3764 out[1][2][1][i+1]=resolutionofA2[i]; 3765 } 3766 } 1180 3767 out[1][1][3]=list(out[1][2][1][1]); 1181 3768 out[1][1][5]=list(out[1][2][1][1]); … … 1187 3774 out[1][1][6]=list(list()); 1188 3775 if (Syzstring=="Sres") 1189 {1190 for (i=1; i<=size(out[1][2][1]); i++)1191 {1192 out[1][1][3][i]=out[1][2][1][i];1193 out[1][1][5][i]=out[1][2][1][i];1194 }1195 }3776 { 3777 for (i=1; i<=size(out[1][2][1]); i++) 3778 { 3779 out[1][1][3][i]=out[1][2][1][i]; 3780 out[1][1][5][i]=out[1][2][1][i]; 3781 } 3782 } 1196 3783 list Hi; 1197 3784 for (i=1; i<=size(out); i++) 1198 {1199 Hi[i]=list(out[i][1][5][1],out[i][1][8][1]);1200 }3785 { 3786 Hi[i]=list(out[i][1][5][1],out[i][1][8][1]); 3787 } 1201 3788 list outall; 1202 3789 outall[1]=out; … … 1207 3794 //////////////////////////////////////////////////////////////////////////////////// 1208 3795 1209 static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #) 3796 static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #) 1210 3797 { 1211 3798 /*this is the Algorithm 3.11 in [W2]*/ … … 1215 3802 def B=basering; 1216 3803 matrix bi=slimgb(transpose(C[5])); 1217 /* Computation of a V_dstrict Groebner basis of C[5]: 1218 if Syzstring=="Vdres" this is done using the method of weighted homogenization1219 (Prop. 3.9 [OT])1220 else we use the homogenized Weyl algebra for Groebner basis computations1221 (Prop 9.9 [OT]),1222 in this case we already compute someresolutions (Thm. 9.10 [OT]) to omit1223 extra Groebner basis computations later on*/3804 /* Computation of a V_dstrict Groebner basis of C[5]: 3805 if Syzstring=="Vdres" this is done using the method of weighted homogenization 3806 (Prop. 3.9 [OT]) 3807 else we use the homogenized Weyl algebra for Groebner basis computations 3808 (Prop 9.9 [OT]), 3809 in this case we already compute someresolutions (Thm. 9.10 [OT]) to omit 3810 extra Groebner basis computations later on*/ 1224 3811 int nr,nc; 1225 3812 intvec i1,i2; 1226 3813 if (Syzstring=="Vdres") 1227 {1228 if(size(#)==0)1229 {1230 matrix J_C=VdStrictGB(C[5],n,list(v));1231 }1232 else1233 {1234 matrix J_C=C[5];//C[5] is already a V_dstrict Groebner basis1235 }1236 }3814 { 3815 if(size(#)==0) 3816 { 3817 matrix J_C=VdStrictGB(C[5],n,list(v)); 3818 } 3819 else 3820 { 3821 matrix J_C=C[5];//C[5] is already a V_dstrict Groebner basis 3822 } 3823 } 1237 3824 else 1238 {1239 if (size(#)==0)1240 {1241 matrix MC=C[5];1242 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, v);1243 setring HomWeyl;1244 matrix J_C=fetch(B,MC);1245 J_C=nHomogenize(J_C);1246 /*computation of V_dstrict resolution of C[5]>needed for proc1247 VdstrictDoubleComplexes*/1248 def ringofSyz=Sres(groebner(transpose(J_C)),lengthofres);1249 setring ringofSyz;1250 matrix J_C=transpose(matrix(RES[2]));1251 J_C=subst(J_C,h,1);1252 logens=ncols(J_C)+1;1253 matrix zerom;1254 list rofC;//will contain resolution of C1255 for (i=3; i<=n+si+1; i++)1256 {1257 if (size(RES)>=i)1258 {1259 zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));1260 if (RES[i]!=zerom)1261 {1262 rofC[i2]=(matrix(RES[i]));1263 1264 if (i==3)1265 {1266 if (nrows(rofC[i2])logens+1!=nrows(J_C))1267 {1268 //build the resolution1269 nr=nrows(J_C)+logens1;1270 nc=ncols(rofC[i2]);1271 rofC[i2]=matrix(rofC[i2],nr,nc);1272 }1273 1274 }1275 if (i!=3)1276 {1277 if (nrows(rofC[i2])logens+1!=nrows(rofC[i3]))1278 {1279 nr=nrows(rofC[i3])+logens1;1280 nc=ncols(rofC[i2]);1281 rofC[i2]=matrix(rofC[i2],nr,nc);1282 }1283 }1284 i1=intvec(logens..nrows(rofC[i2]));1285 i2=intvec(1..ncols(rofC[i2]));1286 rofC[i2]=transpose(submat(rofC[i2],i1,i2));1287 logens=logens+ncols(rofC[i2]);1288 rofC[i2]=subst(rofC[i2],h,1);1289 }1290 else1291 {1292 rofC[i2]=list();1293 }1294 }1295 else1296 {1297 rofC[i2]=list();1298 }1299 }1300 if(size(rofC[1])==0)1301 {1302 omitemptylist=1;1303 }1304 setring B;1305 matrix J_C=fetch(ringofSyz,J_C);1306 if (omitemptylist!=1)1307 {1308 list rofC=fetch(ringofSyz,rofC);1309 }1310 omitemptylist=0;1311 kill HomWeyl;1312 kill ringofSyz;1313 }1314 else1315 {1316 matrix J_C=C[5];//C[5] is already a V_dstrict Groebner basis1317 }1318 }3825 { 3826 if (size(#)==0) 3827 { 3828 matrix MC=C[5]; 3829 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, v); 3830 setring HomWeyl; 3831 matrix J_C=fetch(B,MC); 3832 J_C=nHomogenize(J_C); 3833 /*computation of V_dstrict resolution of C[5]>needed for proc 3834 VdstrictDoubleComplexes*/ 3835 def ringofSyz=Sres(transpose(J_C),lengthofres); 3836 setring ringofSyz; 3837 matrix J_C=transpose(matrix(RES[2])); 3838 J_C=subst(J_C,h,1); 3839 logens=ncols(J_C)+1; 3840 matrix zerom; 3841 list rofC;//will contain resolution of C 3842 for (i=3; i<=n+si+1; i++) 3843 { 3844 if (size(RES)>=i) 3845 { 3846 zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))); 3847 if (RES[i]!=zerom) 3848 { 3849 rofC[i2]=(matrix(RES[i])); 3850 3851 if (i==3) 3852 { 3853 if (nrows(rofC[i2])logens+1!=nrows(J_C)) 3854 { 3855 //build the resolution 3856 nr=nrows(J_C)+logens1; 3857 nc=ncols(rofC[i2]); 3858 rofC[i2]=matrix(rofC[i2],nr,nc); 3859 } 3860 3861 } 3862 if (i!=3) 3863 { 3864 if (nrows(rofC[i2])logens+1!=nrows(rofC[i3])) 3865 { 3866 nr=nrows(rofC[i3])+logens1; 3867 nc=ncols(rofC[i2]); 3868 rofC[i2]=matrix(rofC[i2],nr,nc); 3869 } 3870 } 3871 i1=intvec(logens..nrows(rofC[i2])); 3872 i2=intvec(1..ncols(rofC[i2])); 3873 rofC[i2]=transpose(submat(rofC[i2],i1,i2)); 3874 logens=logens+ncols(rofC[i2]); 3875 rofC[i2]=subst(rofC[i2],h,1); 3876 } 3877 else 3878 { 3879 rofC[i2]=list(); 3880 } 3881 } 3882 else 3883 { 3884 rofC[i2]=list(); 3885 } 3886 } 3887 if(size(rofC[1])==0) 3888 { 3889 omitemptylist=1; 3890 } 3891 setring B; 3892 matrix J_C=fetch(ringofSyz,J_C); 3893 if (omitemptylist!=1) 3894 { 3895 list rofC=fetch(ringofSyz,rofC); 3896 } 3897 omitemptylist=0; 3898 kill HomWeyl; 3899 kill ringofSyz; 3900 } 3901 else 3902 { 3903 matrix J_C=C[5];//C[5] is already a V_dstrict Groebner basis 3904 } 3905 } 1319 3906 /* we compute a V_dstrict Groebner basis for C[3]*/ 1320 3907 matrix J_A=C[1]; … … 1326 3913 matrix Pi[1][ncols(J_AC)]; 1327 3914 for (i=1; i<=nrows(J_C); i++) 1328 {1329 for (j=1; j<=nrows(J_AC);j++)1330 {1331 Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC)));1332 }1333 storePi[i]=Pi;1334 Pi=0;1335 }3915 { 3916 for (j=1; j<=nrows(J_AC);j++) 3917 { 3918 Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC))); 3919 } 3920 storePi[i]=Pi; 3921 Pi=0; 3922 } 1336 3923 /*we compute the shift vector for C[1]*/ 1337 3924 intvec m_a; 1338 3925 list findMin; 1339 int comMin; 3926 int comMin; 1340 3927 for (i=1; i<=ncols(J_A); i++) 1341 {1342 for (j=1; j<=size(storePi);j++)1343 {1344 if (storePi[j][1,i]!=0)1345 {1346 comMin=VdDeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v);1347 comMin=comMinVdDeg(storePi[j][1,i],n,intvec(0));1348 findMin[size(findMin)+1]=comMin;1349 }1350 }1351 if (size(findMin)!=0)1352 {1353 m_a[i]=Min(findMin);1354 findMin=list();1355 }1356 else1357 {1358 m_a[i]=0;1359 }1360 }3928 { 3929 for (j=1; j<=size(storePi);j++) 3930 { 3931 if (storePi[j][1,i]!=0) 3932 { 3933 comMin=VdDeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v); 3934 comMin=comMinVdDeg(storePi[j][1,i],n,intvec(0)); 3935 findMin[size(findMin)+1]=comMin; 3936 } 3937 } 3938 if (size(findMin)!=0) 3939 { 3940 m_a[i]=Min(findMin); 3941 findMin=list(); 3942 } 3943 else 3944 { 3945 m_a[i]=0; 3946 } 3947 } 1361 3948 matrix zero[ncols(J_A)][ncols(J_C)]; 1362 3949 matrix g_AB=concat(unitmat(ncols(J_A)),zero); 1363 3950 matrix g_BC= transpose(concat(transpose(zero),transpose(unitmat(ncols(J_C))))); 1364 3951 intvec m_b=m_a,v; 1365 /* computation of a V_dstrict Groebner basis of C[1] (and resolution if 1366 Syzstring=='Vdres') */3952 /* computation of a V_dstrict Groebner basis of C[1] (and resolution if 3953 Syzstring=='Vdres') */ 1367 3954 if (Syzstring=="Vdres") 1368 {1369 J_A=VdStrictGB(J_A,n,m_a);1370 }3955 { 3956 J_A=VdStrictGB(J_A,n,m_a); 3957 } 1371 3958 else 1372 {1373 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_a);1374 setring HomWeyl;1375 matrix J_A=fetch(B,J_A);1376 J_A=nHomogenize(J_A);1377 def ringofSyz=Sres(groebner(transpose(J_A)),lengthofres);1378 setring ringofSyz;1379 matrix J_A=transpose(matrix(RES[2]));1380 matrix zerom;1381 J_A=subst(J_A,h,1);1382 logens=ncols(J_A)+1;1383 list rofA;1384 for (i=3; i<=n+si+1; i++)1385 {1386 if (size(RES)>=i)1387 {1388 zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));1389 if (RES[i]!=zerom)1390 {1391 rofA[i2]=matrix(RES[i]);// resolution for C[1]1392 if (i==3)1393 {1394 if (nrows(rofA[i2])logens+1!=nrows(J_A))1395 {1396 nr=nrows(J_A)+logens1;1397 nc=ncols(rofA[i2]);1398 rofA[i2]=matrix(rofA[i2],nr,nc);1399 }1400 }1401 if (i!=3)1402 {1403 if (nrows(rofA[i2])logens+1!=nrows(rofA[i3]))1404 {1405 nr=nrows(rofA[i3])+logens1;1406 nc=ncols(rofA[i2]);1407 rofA[i2]=matrix(rofA[i2],nr,nc);1408 }1409 }1410 i1=intvec(logens..nrows(rofA[i2]));1411 i2=intvec(1..ncols(rofA[i2]));1412 rofA[i2]=transpose(submat(rofA[i2],i1,i2));1413 logens=logens+ncols(rofA[i2]);1414 rofA[i2]=subst(rofA[i2],h,1);1415 }1416 else1417 {1418 rofA[i2]=list();1419 }1420 }1421 else1422 {1423 rofA[i2]=list();1424 }1425 }1426 if(size(rofA[1])==0)1427 {1428 omitemptylist=1;1429 }1430 setring B;1431 J_A=fetch(ringofSyz,J_A);1432 if (omitemptylist!=1)1433 {1434 list rofA=fetch(ringofSyz,rofA);1435 }1436 omitemptylist=0;1437 kill HomWeyl;1438 kill ringofSyz;1439 }3959 { 3960 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_a); 3961 setring HomWeyl; 3962 matrix J_A=fetch(B,J_A); 3963 J_A=nHomogenize(J_A); 3964 def ringofSyz=Sres(transpose(J_A),lengthofres); 3965 setring ringofSyz; 3966 matrix J_A=transpose(matrix(RES[2])); 3967 matrix zerom; 3968 J_A=subst(J_A,h,1); 3969 logens=ncols(J_A)+1; 3970 list rofA; 3971 for (i=3; i<=n+si+1; i++) 3972 { 3973 if (size(RES)>=i) 3974 { 3975 zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))); 3976 if (RES[i]!=zerom) 3977 { 3978 rofA[i2]=matrix(RES[i]);// resolution for C[1] 3979 if (i==3) 3980 { 3981 if (nrows(rofA[i2])logens+1!=nrows(J_A)) 3982 { 3983 nr=nrows(J_A)+logens1; 3984 nc=ncols(rofA[i2]); 3985 rofA[i2]=matrix(rofA[i2],nr,nc); 3986 } 3987 } 3988 if (i!=3) 3989 { 3990 if (nrows(rofA[i2])logens+1!=nrows(rofA[i3])) 3991 { 3992 nr=nrows(rofA[i3])+logens1; 3993 nc=ncols(rofA[i2]); 3994 rofA[i2]=matrix(rofA[i2],nr,nc); 3995 } 3996 } 3997 i1=intvec(logens..nrows(rofA[i2])); 3998 i2=intvec(1..ncols(rofA[i2])); 3999 rofA[i2]=transpose(submat(rofA[i2],i1,i2)); 4000 logens=logens+ncols(rofA[i2]); 4001 rofA[i2]=subst(rofA[i2],h,1); 4002 } 4003 else 4004 { 4005 rofA[i2]=list(); 4006 } 4007 } 4008 else 4009 { 4010 rofA[i2]=list(); 4011 } 4012 } 4013 if(size(rofA[1])==0) 4014 { 4015 omitemptylist=1; 4016 } 4017 setring B; 4018 J_A=fetch(ringofSyz,J_A); 4019 if (omitemptylist!=1) 4020 { 4021 list rofA=fetch(ringofSyz,rofA); 4022 } 4023 omitemptylist=0; 4024 kill HomWeyl; 4025 kill ringofSyz; 4026 } 1440 4027 J_AC=transpose(storePi[1]); 1441 4028 for (i=2; i<= size(storePi); i++) 1442 {1443 J_AC=concat(J_AC, transpose(storePi[i]));1444 }4029 { 4030 J_AC=concat(J_AC, transpose(storePi[i])); 4031 } 1445 4032 J_AC=transpose(concat(transpose(matrix(J_A,nrows(J_A),nrows(J_AC))),J_AC)); 1446 4033 list Vdstrict=(list(J_A),list(g_AB),list(J_AC),list(g_BC),list(J_C),list(m_a)); 1447 4034 Vdstrict[7]=list(m_b); 1448 Vdstrict[8]=list(v); 4035 Vdstrict[8]=list(v); 1449 4036 if(Syzstring=="Sres") 1450 {1451 Vdstrict[9]=rofA;1452 if(size(#)==0)1453 {1454 Vdstrict[10]=rofC;1455 }1456 }4037 { 4038 Vdstrict[9]=rofA; 4039 if(size(#)==0) 4040 { 4041 Vdstrict[10]=rofC; 4042 } 4043 } 1457 4044 return (Vdstrict); 1458 4045 } … … 1460 4047 //////////////////////////////////////////////////////////////////////////////////// 1461 4048 1462 static proc toVdStrictSequences (list L,int d,intvec v,string Syzstring,int sizeL) 4049 static proc toVdStrictSequences (list L,int d,intvec v,string Syzstring,int sizeL) 1463 4050 { 1464 /* this is Argorithm 3.11 combined with Lemma 4.2 in [W2] for two short exact 1465 pieces.1466 We asume that we are given two sequences of the form coker(L[i][1])>1467 coker(L[i][3])>coker(L[i][5]) with differentials L[i][2] and L[i][4] such1468 that L[1][3]=L[2][1].We are going to transform them to V_dstrict sequences1469 J_D>J_A>J_F and J_A>J_B>J_C*/4051 /* this is Argorithm 3.11 combined with Lemma 4.2 in [W2] for two short exact 4052 pieces. 4053 We asume that we are given two sequences of the form coker(L[i][1])> 4054 coker(L[i][3])>coker(L[i][5]) with differentials L[i][2] and L[i][4] such 4055 that L[1][3]=L[2][1].We are going to transform them to V_dstrict sequences 4056 J_D>J_A>J_F and J_A>J_B>J_C*/ 1470 4057 int omitemptylist; 1471 4058 int lengthofres=sizeL+d1; … … 1475 4062 matrix J_D=L[1][1]; 1476 4063 matrix f_FA=L[1][4]; 1477 /*We find new presentations coker(J_DF) and coker(J_DFC) for L[1][4]=L[2][1] 1478 and L[2][4],resp. such that ncols(L[i][1])+ncols(L[i][5])=ncols(L[i][3]) */4064 /*We find new presentations coker(J_DF) and coker(J_DFC) for L[1][4]=L[2][1] 4065 and L[2][4],resp. such that ncols(L[i][1])+ncols(L[i][5])=ncols(L[i][3]) */ 1479 4066 matrix f_DFA=transpose(concat(transpose(L[1][2]),transpose(f_FA))); 1480 4067 matrix J_DF=divdr(f_DFA,L[1][3]);//coker(J_DF) is isomorphic to coker(L[2][1]); … … 1483 4070 matrix f_DFCB=transpose(concat(transpose(f_DFA*L[2][2]),transpose(f_CB))); 1484 4071 matrix J_DFC=divdr(f_DFCB,L[2][3]);//coker(J_DFC) are coker(L[2][3)]) isomorphic 1485 /* find a shift vector on the range of J_F such that the first sequence is 1486 exact*/4072 /* find a shift vector on the range of J_F such that the first sequence is 4073 exact*/ 1487 4074 matrix P=matrixLift(J_DFC*prodr(ncols(J_DF),ncols(L[2][5])),J_C); 1488 4075 list storePi; … … 1490 4077 int i; int j; 1491 4078 for (i=1; i<=nrows(J_C); i++) 1492 {1493 for (j=1; j<=nrows(J_DFC);j++)1494 {1495 Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC)));1496 }1497 storePi[i]=Pi;1498 Pi=0;1499 }4079 { 4080 for (j=1; j<=nrows(J_DFC);j++) 4081 { 4082 Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC))); 4083 } 4084 storePi[i]=Pi; 4085 Pi=0; 4086 } 1500 4087 intvec m_a; 1501 4088 list findMin; … … 1505 4092 intvec i1,i2; 1506 4093 for (i=1; i<=ncols(J_DF); i++) 1507 {1508 for (j=1; j<=size(storePi);j++)1509 {1510 if (storePi[j][1,i]!=0)1511 {1512 comMin=VdDeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v);1513 comMin=comMinVdDeg(storePi[j][1,i],d,intvec(0));1514 findMin[size(findMin)+1]=comMin;1515 }1516 }1517 if (size(findMin)!=0)1518 {1519 m_a[i]=Min(findMin);// shift vector for L[2][1]1520 findMin=list();1521 noMin[i]=0;1522 }1523 else1524 {1525 noMin[i]=1;1526 }1527 }4094 { 4095 for (j=1; j<=size(storePi);j++) 4096 { 4097 if (storePi[j][1,i]!=0) 4098 { 4099 comMin=VdDeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v); 4100 comMin=comMinVdDeg(storePi[j][1,i],d,intvec(0)); 4101 findMin[size(findMin)+1]=comMin; 4102 } 4103 } 4104 if (size(findMin)!=0) 4105 { 4106 m_a[i]=Min(findMin);// shift vector for L[2][1] 4107 findMin=list(); 4108 noMin[i]=0; 4109 } 4110 else 4111 { 4112 noMin[i]=1; 4113 } 4114 } 1528 4115 if (size(m_a) < ncols(J_DF)) 1529 {1530 m_a[ncols(J_DF)]=0;1531 }4116 { 4117 m_a[ncols(J_DF)]=0; 4118 } 1532 4119 intvec m_f=m_a[ncols(J_D)+1..size(m_a)]; 1533 /* Computation of a V_dstrict Groebner basis of J_F=L[1][5]: 1534 if Syzstring=="Vdres" this is done using the method of weighted homogenization1535 (Prop. 3.9 [OT])1536 else we use the homogenized Weyl algerba for Groebner basis computations1537 (Prop 9.9 [OT]), in this case we already compute resolutions1538 (Thm. 9.10 in [OT]) to omit extra Groebner basis computations later on*/4120 /* Computation of a V_dstrict Groebner basis of J_F=L[1][5]: 4121 if Syzstring=="Vdres" this is done using the method of weighted homogenization 4122 (Prop. 3.9 [OT]) 4123 else we use the homogenized Weyl algerba for Groebner basis computations 4124 (Prop 9.9 [OT]), in this case we already compute resolutions 4125 (Thm. 9.10 in [OT]) to omit extra Groebner basis computations later on*/ 1539 4126 if (Syzstring=="Vdres") 1540 {1541 J_F=VdStrictGB(J_F,d,m_f);1542 }4127 { 4128 J_F=VdStrictGB(J_F,d,m_f); 4129 } 1543 4130 else 1544 {1545 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_f);1546 setring HomWeyl;1547 matrix J_F=fetch(B,J_F);1548 J_F=nHomogenize(J_F);1549 def ringofSyz=Sres(groebner(transpose(J_F)),lengthofres);1550 setring ringofSyz;1551 matrix J_F=transpose(matrix(RES[2]));1552 J_F=subst(J_F,h,1);1553 logens=ncols(J_F)+1;1554 list rofF;1555 for (i=3; i<=d+sizeL+1; i++)1556 {1557 if (size(RES)>=i)1558 {1559 if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))1560 {1561 rofF[i2]=(matrix(RES[i]));// resolution for J_F1562 if (i==3)1563 {1564 if (nrows(rofF[i2])logens+1!=nrows(J_F))1565 {1566 nr=nrows(J_F)+logens1;1567 nc=ncols(rofF[i2]);1568 rofF[i2]=matrix(rofF[i2],nr,nc);1569 }1570 }1571 if (i!=3)1572 {1573 if (nrows(rofF[i2])logens+1!=nrows(rofF[i3]))1574 {1575 nr=nrows(rofF[i3])+logens1;1576 rofF[i2]=matrix(rofF[i2],nr,ncols(rofF[i2]));1577 }1578 }1579 i1=intvec(logens..nrows(rofF[i2]));1580 i2=intvec(1..ncols(rofF[i2]));1581 rofF[i2]=transpose(submat(rofF[i2],i1,i2));1582 logens=logens+ncols(rofF[i2]);1583 rofF[i2]=subst(rofF[i2],h,1);1584 }1585 else1586 {1587 rofF[i2]=list();1588 }1589 }1590 else1591 {1592 rofF[i2]=list();1593 }1594 }1595 if(size(rofF[1])==0)1596 {1597 omitemptylist=1;1598 }1599 setring B;1600 J_F=fetch(ringofSyz,J_F);1601 if (omitemptylist!=1)1602 {1603 list rofF=fetch(ringofSyz,rofF);1604 }1605 omitemptylist=0;1606 kill HomWeyl;1607 kill ringofSyz;1608 }4131 { 4132 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_f); 4133 setring HomWeyl; 4134 matrix J_F=fetch(B,J_F); 4135 J_F=nHomogenize(J_F); 4136 def ringofSyz=Sres(transpose(J_F),lengthofres); 4137 setring ringofSyz; 4138 matrix J_F=transpose(matrix(RES[2])); 4139 J_F=subst(J_F,h,1); 4140 logens=ncols(J_F)+1; 4141 list rofF; 4142 for (i=3; i<=d+sizeL+1; i++) 4143 { 4144 if (size(RES)>=i) 4145 { 4146 if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])))) 4147 { 4148 rofF[i2]=(matrix(RES[i]));// resolution for J_F 4149 if (i==3) 4150 { 4151 if (nrows(rofF[i2])logens+1!=nrows(J_F)) 4152 { 4153 nr=nrows(J_F)+logens1; 4154 nc=ncols(rofF[i2]); 4155 rofF[i2]=matrix(rofF[i2],nr,nc); 4156 } 4157 } 4158 if (i!=3) 4159 { 4160 if (nrows(rofF[i2])logens+1!=nrows(rofF[i3])) 4161 { 4162 nr=nrows(rofF[i3])+logens1; 4163 rofF[i2]=matrix(rofF[i2],nr,ncols(rofF[i2])); 4164 } 4165 } 4166 i1=intvec(logens..nrows(rofF[i2])); 4167 i2=intvec(1..ncols(rofF[i2])); 4168 rofF[i2]=transpose(submat(rofF[i2],i1,i2)); 4169 logens=logens+ncols(rofF[i2]); 4170 rofF[i2]=subst(rofF[i2],h,1); 4171 } 4172 else 4173 { 4174 rofF[i2]=list(); 4175 } 4176 } 4177 else 4178 { 4179 rofF[i2]=list(); 4180 } 4181 } 4182 if(size(rofF[1])==0) 4183 { 4184 omitemptylist=1; 4185 } 4186 setring B; 4187 J_F=fetch(ringofSyz,J_F); 4188 if (omitemptylist!=1) 4189 { 4190 list rofF=fetch(ringofSyz,rofF); 4191 } 4192 omitemptylist=0; 4193 kill HomWeyl; 4194 kill ringofSyz; 4195 } 1609 4196 /*find shift vectors on the range of J_D*/ 1610 4197 P=matrixLift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F); … … 1612 4199 matrix Pidf[1][ncols(J_DF)]; 1613 4200 for (i=1; i<=nrows(J_F); i++) 1614 {1615 for (j=1; j<=nrows(J_DF);j++)1616 {1617 Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF)));1618 }1619 storePinew[i]=Pidf;1620 Pidf=0;1621 }4201 { 4202 for (j=1; j<=nrows(J_DF);j++) 4203 { 4204 Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF))); 4205 } 4206 storePinew[i]=Pidf; 4207 Pidf=0; 4208 } 1622 4209 intvec m_d; 1623 4210 for (i=1; i<=ncols(J_D); i++) 1624 { 1625 for (j=1; j<=size(storePinew);j++) 1626 { 1627 if (storePinew[j][1,i]!=0) 1628 { 1629 comMin=VdDeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f); 1630 comMin=comMinVdDeg(storePinew[j][1,i],d,intvec(0)); 1631 findMin[size(findMin)+1]=comMin; 1632 } 1633 } 1634 if (size(findMin)!=0) 1635 { 1636 if (noMin[i]==0) 1637 { 1638 m_d[i]=Min(insert(findMin,m_a[i])); 1639 m_a[i]=m_d[i]; 1640 } 4211 { 4212 for (j=1; j<=size(storePinew);j++) 4213 { 4214 if (storePinew[j][1,i]!=0) 4215 { 4216 comMin=VdDeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f); 4217 comMin=comMinVdDeg(storePinew[j][1,i],d,intvec(0)); 4218 findMin[size(findMin)+1]=comMin; 4219 } 4220 } 4221 if (size(findMin)!=0) 4222 { 4223 if (noMin[i]==0) 4224 { 4225 m_d[i]=Min(insert(findMin,m_a[i])); 4226 m_a[i]=m_d[i]; 4227 } 4228 else 4229 { 4230 m_d[i]=Min(findMin); 4231 m_a[i]=m_d[i]; 4232 } 4233 } 1641 4234 else 1642 { 1643 m_d[i]=Min(findMin); 1644 m_a[i]=m_d[i]; 1645 } 1646 } 1647 else 1648 { 1649 m_d[i]=m_a[i]; 1650 } 1651 findMin=list(); 1652 } 1653 /* compute a V_dstrict Groebner basis (and resolution of J_D if 1654 Syzstring!='Vdres') for J_D*/ 4235 { 4236 m_d[i]=m_a[i]; 4237 } 4238 findMin=list(); 4239 } 4240 /* compute a V_dstrict Groebner basis (and resolution of J_D if 4241 Syzstring!='Vdres') for J_D*/ 1655 4242 if (Syzstring=="Vdres") 1656 {1657 J_D=VdStrictGB(J_D,d,m_d);1658 }4243 { 4244 J_D=VdStrictGB(J_D,d,m_d); 4245 } 1659 4246 else 1660 {1661 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_d);1662 setring HomWeyl;1663 matrix J_D=fetch(B,J_D);1664 J_D=nHomogenize(J_D);1665 def ringofSyz=Sres(groebner(transpose(J_D)),lengthofres);1666 setring ringofSyz;1667 matrix J_D=transpose(matrix(RES[2]));1668 J_D=subst(J_D,h,1);1669 logens=ncols(J_D)+1;1670 list rofD;1671 for (i=3; i<=d+sizeL+1; i++)1672 {1673 if (size(RES)>=i)1674 {1675 if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))1676 {1677 rofD[i2]=(matrix(RES[i]));// resolution for J_D1678 if (i==3)1679 {1680 if (nrows(rofD[i2])logens+1!=nrows(J_D))1681 {1682 nr=nrows(J_D)+logens1;1683 rofD[i2]=matrix(rofD[i2],nr,ncols(rofD[i2]));1684 }1685 }1686 if (i!=3)1687 {1688 if (nrows(rofD[i2])logens+1!=nrows(rofD[i3]))1689 {1690 nr=nrows(rofD[i3])+logens1;1691 rofD[i2]=matrix(rofD[i2],nr,ncols(rofD[i2]));1692 }1693 }1694 i1=intvec(logens..nrows(rofD[i2]));1695 i2=intvec(1..ncols(rofD[i2]));1696 rofD[i2]=transpose(submat(rofD[i2],i1,i2));1697 logens=logens+ncols(rofD[i2]);1698 rofD[i2]=subst(rofD[i2],h,1);1699 }1700 else1701 {1702 rofD[i2]=list();1703 }1704 }1705 else1706 {1707 rofD[i2]=list();1708 }1709 }1710 if(size(rofD[1])==0)1711 {1712 omitemptylist=1;1713 }1714 setring B;1715 J_D=fetch(ringofSyz,J_D);1716 if (omitemptylist!=1)1717 {1718 list rofD=fetch(ringofSyz,rofD);1719 }1720 omitemptylist=0;1721 kill HomWeyl;1722 kill ringofSyz;1723 }1724 /* compute new matrices for J_A and J_B such that their rows form a V_dstrict 1725 Groebner basis and nrows(J_A)=nrows(J_D)+nrows(J_F) and1726 nrows(J_B)=nrows(J_A)+nrows(J_C)*/4247 { 4248 def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_d); 4249 setring HomWeyl; 4250 matrix J_D=fetch(B,J_D); 4251 J_D=nHomogenize(J_D); 4252 def ringofSyz=Sres(transpose(J_D),lengthofres); 4253 setring ringofSyz; 4254 matrix J_D=transpose(matrix(RES[2])); 4255 J_D=subst(J_D,h,1); 4256 logens=ncols(J_D)+1; 4257 list rofD; 4258 for (i=3; i<=d+sizeL+1; i++) 4259 { 4260 if (size(RES)>=i) 4261 { 4262 if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])))) 4263 { 4264 rofD[i2]=(matrix(RES[i]));// resolution for J_D 4265 if (i==3) 4266 { 4267 if (nrows(rofD[i2])logens+1!=nrows(J_D)) 4268 { 4269 nr=nrows(J_D)+logens1; 4270 rofD[i2]=matrix(rofD[i2],nr,ncols(rofD[i2])); 4271 } 4272 } 4273 if (i!=3) 4274 { 4275 if (nrows(rofD[i2])logens+1!=nrows(rofD[i3])) 4276 { 4277 nr=nrows(rofD[i3])+logens1; 4278 rofD[i2]=matrix(rofD[i2],nr,ncols(rofD[i2])); 4279 } 4280 } 4281 i1=intvec(logens..nrows(rofD[i2])); 4282 i2=intvec(1..ncols(rofD[i2])); 4283 rofD[i2]=transpose(submat(rofD[i2],i1,i2)); 4284 logens=logens+ncols(rofD[i2]); 4285 rofD[i2]=subst(rofD[i2],h,1); 4286 } 4287 else 4288 { 4289 rofD[i2]=list(); 4290 } 4291 } 4292 else 4293 { 4294 rofD[i2]=list(); 4295 } 4296 } 4297 if(size(rofD[1])==0) 4298 { 4299 omitemptylist=1; 4300 } 4301 setring B; 4302 J_D=fetch(ringofSyz,J_D); 4303 if (omitemptylist!=1) 4304 { 4305 list rofD=fetch(ringofSyz,rofD); 4306 } 4307 omitemptylist=0; 4308 kill HomWeyl; 4309 kill ringofSyz; 4310 } 4311 /* compute new matrices for J_A and J_B such that their rows form a V_dstrict 4312 Groebner basis and nrows(J_A)=nrows(J_D)+nrows(J_F) and 4313 nrows(J_B)=nrows(J_A)+nrows(J_C)*/ 1727 4314 J_DF=transpose(storePinew[1]); 1728 4315 for (i=2; i<=nrows(J_F); i++) 1729 {1730 J_DF=concat(J_DF,transpose(storePinew[i]));1731 }4316 { 4317 J_DF=concat(J_DF,transpose(storePinew[i])); 4318 } 1732 4319 J_DF=transpose(concat(transpose(matrix(J_D,nrows(J_D),nrows(J_DF))),J_DF)); 1733 4320 J_DFC=transpose(storePi[1]); 1734 4321 for (i=2; i<=nrows(J_C); i++) 1735 {1736 J_DFC=concat(J_DFC,transpose(storePi[i]));1737 }4322 { 4323 J_DFC=concat(J_DFC,transpose(storePi[i])); 4324 } 1738 4325 J_DFC=transpose(concat(transpose(matrix(J_DF,nrows(J_DF),nrows(J_DFC))),J_DFC)); 1739 4326 intvec m_b=m_a,v; … … 1744 4331 matrix g_AB=concat(unitmat(ncols(J_DF)),zero1); 1745 4332 matrix g_BC=transpose(concat(transpose(zero1),unitmat(ncols(J_C)))); 1746 list out; 4333 list out; 1747 4334 out[1]=list(list(J_D),list(g_DA),list(J_DF),list(g_AF),list(J_F)); 1748 4335 out[1]=out[1]+list(list(m_d),list(m_a),list(m_f)); … … 1750 4337 out[2]=out[2]+list(list(m_a),list(m_b),list(v)); 1751 4338 if (Syzstring=="Sres") 1752 {1753 out[3]=rofD;1754 out[4]=rofF;1755 }4339 { 4340 out[3]=rofD; 4341 out[4]=rofF; 4342 } 1756 4343 return(out); 1757 4344 } … … 1761 4348 static proc VdStrictDoubleComplexes(list L,int d,string Syzstring) 1762 4349 { 1763 /* We compute V_dstrict resolutions over the V_dstrict short exact pieces from 1764 the procedure shortExactPiecesToVdStrict.1765 We use Algorithms 3.14 and 3.15 in [W2]*/4350 /* We compute V_dstrict resolutions over the V_dstrict short exact pieces from 4351 the procedure shortExactPiecesToVdStrict. 4352 We use Algorithms 3.14 and 3.15 in [W2]*/ 1766 4353 int i,k,c,j,l,totaldeg,comparedegs,SBcom,verk; 1767 intvec fordegs; 4354 intvec fordegs; 1768 4355 intvec n_b,i1,i2; 1769 4356 matrix rem,forML,subm,zerom,unitm,subm2; … … 1777 4364 list forhW; 1778 4365 if (Syzstring=="Sres") 1779 { 1780 /*we already computed some of the resolutions in the procedure 1781 shortExactPiecesToVdStrict*/ 1782 matrix Pold,Pnew,Picombined; intvec containsndeg; matrix Pinew; 1783 for (k=1; k<=(size(L)+d1); k++) 1784 { 4366 { 4367 /*we already computed some of the resolutions in the procedure 4368 shortExactPiecesToVdStrict*/ 4369 matrix Pold,Pnew,Picombined; intvec containsndeg; matrix Pinew; 4370 for (k=1; k<=(size(L)+d1); k++) 4371 { 4372 L[1][1][1][k+1]=list(); 4373 L[1][1][2][k+1]=list(); 4374 L[1][1][6][k+1]=list(); 4375 } 4376 L[1][1][6][size(L)+d+1]=list(); 4377 matrix mem; 4378 for (i=2; i<=d+size(L)+1; i++) 4379 {; 4380 v=0; 4381 if(size(L[1][1][3][i1])!=0) 4382 { 4383 if(i!=d+size(L)+1) 4384 { 4385 /*horizontal differential*/ 4386 L[1][1][4][i1]=unitmat(nrows(L[1][1][3][i1])); 4387 } 4388 for (j=1; j<=nrows(L[1][1][3][i1]); j++) 4389 { 4390 mem=submat(L[1][1][3][i1],j,intvec(1..ncols(L[1][1][3][i1]))); 4391 v[j]=VdDeg(mem,d,L[1][1][7][i1]); 4392 } 4393 L[1][1][7][i]=v;//new shift vector 4394 L[1][1][8][i]=v; 4395 L[1][2][6][i]=v; 4396 } 4397 else 4398 { 4399 if (i!=d+size(L)+1) 4400 { 4401 L[1][1][4][i1]=list(); 4402 } 4403 L[1][1][7][i]=list(); 4404 L[1][1][8][i]=list(); 4405 L[1][2][6][i]=list(); 4406 } 4407 } 4408 if (size(L[1][1][3][d+size(L)])!=0) 4409 { 4410 /*horizontal differential*/ 4411 L[1][1][4][d+size(L)]=unitmat(nrows(L[1][1][3][d+size(L)])); 4412 } 4413 else 4414 { 4415 L[1][1][4][d+size(L)]=list(); 4416 } 4417 for (k=1; k<size(L); k++) 4418 { 4419 /* We build a V_dstrict resolution for coker(L[k][2][1][1])> 4420 coker(L[k][2][3][1])>coker(L[k][2][5][1]) using the resolution 4421 obtained for coker(L[k][1][3][1]). 4422 L[k][2][i][j] will be the jth module in the resolution of L[k][2][i][1] 4423 for i=1,3,5. 4424 L[k][2][i+5][j] will be the jth shift vector in the resolution of 4425 L[k][2][i][1](this holds also for the case Syzstring=="Vdres")*/ 4426 for (i=2; i<=d+size(L); i++) 4427 { 4428 v=0; 4429 if (size(L[k][2][5][i1])!=0) 4430 { 4431 for (j=1; j<=nrows(L[k][2][5][i1]); j++) 4432 { 4433 i1=intvec(1..ncols(L[k][2][5][i1])); 4434 mem=submat(L[k][2][5][i1],j,i1); 4435 v[j]=VdDeg(mem,d,L[k][2][8][i1]); 4436 } 4437 /*next shift vector in th resolution of coker(L[k][2][5][1])*/ 4438 L[k][2][8][i]=v; 4439 } 4440 else 4441 { 4442 L[k][2][8][i]=list(); 4443 } 4444 /* we build step by step a resolution for coker(L[k][2][5][1]) using 4445 the resolutions of coker(L[k][2][1][1]) and coker(L[k][2][5][1])*/ 4446 if (size(L[k][2][5][i])!=0) 4447 { 4448 if (size(L[k][2][1][i])!=0 or size(L[k][2][1][i1])!=0) 4449 { 4450 L[k][2][3][i]=transpose(syz(transpose(L[k][2][3][i1]))); 4451 nr= nrows(L[k][2][1][i1]); 4452 nc=ncols(L[k][2][5][i]); 4453 Pold=matrixLift(L[k][2][3][i]*prodr(nr,nc), L[k][2][5][i]); 4454 matrix Pi[1][ncols(L[k][2][3][i])]; 4455 for (l=1; l<=nrows(L[k][2][5][i]); l++) 4456 { 4457 for (j=1; j<=nrows(L[k][2][3][i]); j++) 4458 { 4459 i2=intvec(1..ncols(L[k][2][3][i])); 4460 Pi=Pi+Pold[l,j]*submat(L[k][2][3][i],j,i2); 4461 } 4462 if (l==1) 4463 { 4464 Picombined=transpose(Pi); 4465 } 4466 else 4467 { 4468 Picombined=concat(Picombined,transpose(Pi)); 4469 } 4470 Pi=0; 4471 } 4472 kill Pi; 4473 Picombined=transpose(Picombined); 4474 if (size(L[k][2][1][i])!=0) 4475 { 4476 if (i==2) 4477 { 4478 containsndeg=(0:ncols(L[k][2][1][1])); 4479 } 4480 containsndeg=nDeg(L[k][2][1][i1],containsndeg); 4481 forhW=list(L[k][2][6][i],containsndeg); 4482 def HomWeyl=makeHomogenizedWeyl(n,forhW); 4483 setring HomWeyl; 4484 list L=fetch(B,L); 4485 matrix M=L[k][2][1][i]; 4486 module Mmod; 4487 list forM=nHomogenize(M,containsndeg,1); 4488 M=forM[1]; 4489 totaldeg=forM[2]; 4490 kill forM; 4491 matrix Maorig=fetch(B,Picombined); 4492 matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M))); 4493 matrix mem,subm,zerom; 4494 matrix Pinew; 4495 M=transpose(M); 4496 SBcom=0; 4497 for (l=1; l<=nrows(Ma); l++) 4498 { 4499 zerom=matrix(0,1,(ncols(Maorig)ncols(Ma))); 4500 i1=(ncols(Ma)+1..ncols(Maorig)); 4501 if (submat(Maorig,l,i1)==zerom) 4502 { 4503 for (cc=1; cc<=ncols(Ma); cc++) 4504 { 4505 Maorig[l,cc]=0; 4506 } 4507 } 4508 i2=(ncols(Ma)+1..ncols(Maorig)); 4509 i1=(1..ncols(Ma)); 4510 if (VdDeg(submat(Maorig,l,i1),d,L[k][2][6][i])> 4511 VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]) and 4512 submat(Maorig,l,i1)!=matrix(0,1,ncols(Ma))) 4513 { 4514 /*V_dGrad is to big> we make it smaller using 4515 Vdnormal form computations*/ 4516 if (SBcom==0) 4517 { 4518 Mmod=slimgb(M); 4519 M=Mmod; 4520 SBcom=1; 4521 } 4522 //print("Reduzierung des V_dGrades(Stelle1)"); 4523 i2=(ncols(Ma)+1..ncols(Maorig)); 4524 vd1=VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]); 4525 mem=submat(Ma,l,(1..ncols(Ma))); 4526 mem=nHomogenize(mem,containsndeg); 4527 mem=h^totaldeg*mem; 4528 mem=transpose(mem); 4529 mem=reduce(mem,Mod);////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 4530 matrix jt=transpose(subst(mem,h,1)); 4531 setring B; 4532 matrix jt=fetch(HomWeyl,jt); 4533 matrix need=fetch(HomWeyl,Maorig); 4534 need=submat(need,l,(1..ncols(need))); 4535 i1=L[k][2][6][i]; 4536 i2=L[k][2][8][i]; 4537 jt=VdNormalForm(need,L[k][2][1][i],d,i1,i2); 4538 setring HomWeyl; 4539 mem=fetch(B,jt); 4540 mem=transpose(mem); 4541 if (l==1) 4542 { 4543 Pinew=mem; 4544 } 4545 else 4546 { 4547 Pinew=concat(Pinew,mem); 4548 } 4549 vd2=VdDeg(transpose(mem),d,L[k][2][6][i]); 4550 if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem))) 4551 {//should not happen!! 4552 //print("Reduzierung fehlgeschlagen!!(Stelle1)"); 4553 } 4554 } 4555 else 4556 { 4557 if (l==1) 4558 { 4559 Pinew=transpose(submat(Ma,l,(1..ncols(Ma)))); 4560 } 4561 else 4562 { 4563 subm=transpose(submat(Ma,l,(1..ncols(Ma)))); 4564 Pinew=concat(Pinew,subm); 4565 } 4566 } 4567 } 4568 Pinew=subst(Pinew,h,1); 4569 Pinew=transpose(Pinew); 4570 setring B; 4571 Pinew=fetch(HomWeyl,Pinew); 4572 kill HomWeyl; 4573 L[k][2][3][i]=concat(Pinew,L[k][2][5][i]); 4574 subm=transpose(L[k][2][3][i]); 4575 subm=concat(transpose(L[k][2][1][i]),subm); 4576 L[k][2][3][i]=transpose(subm); 4577 } 4578 else 4579 { 4580 L[k][2][3][i]=Picombined; 4581 } 4582 L[k+1][1][1][i]=L[k][2][5][i]; 4583 nr=nrows(L[k][2][1][i1]); 4584 nc=ncols(L[k][2][5][i]); 4585 L[k][2][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); 4586 L[k][2][4][i]=prodr(nrows(L[k][2][1][i1]),nc); 4587 v=L[k][2][6][i],L[k][2][8][i]; 4588 L[k][2][7][i]=v; 4589 L[k+1][1][6][i]=L[k][2][8][i]; 4590 } 4591 else 4592 { 4593 L[k][2][3][i]=L[k][2][5][i]; 4594 L[k][2][2][i]=list(); 4595 L[k][2][7][i]=L[k][2][8][i]; 4596 L[k][2][4][i]=unitmat(nrows(L[k][2][5][i1])); 4597 L[k+1][1][6][i]=L[k][2][8][i]; 4598 L[k+1][1][1][i]=L[k][2][5][i]; 4599 } 4600 } 4601 else 4602 { 4603 if (size(L[k][2][1][i])!=0) 4604 { 4605 if (size(L[k][2][5][i1])!=0) 4606 { 4607 nr=nrows(L[k][2][5][i1]); 4608 L[k][2][3][i]=concat(L[k][2][1][i],matrix(0,1,nr)); 4609 v=L[k][2][6][i],L[k][2][8][i]; 4610 L[k][2][7][i]=v; 4611 nc=nrows(L[k][2][1][i1]); 4612 L[k][2][2][i]=concat(unitmat(nc),matrix(0,nc,nr)); 4613 L[k][2][4][i]=prodr(nrows(L[k][2][1][i1]),nr); 4614 } 4615 else 4616 { 4617 L[k][2][3][i]=L[k][2][1][i]; 4618 L[k][2][7][i]=L[k][2][6][i]; 4619 L[k][2][2][i]=unitmat(nrows(L[k][2][1][i1])); 4620 L[k][2][4][i]=list(); 4621 } 4622 L[k+1][1][1][i]=L[k][2][5][i]; 4623 L[k+1][1][6][i]=L[k][2][8][i]; 4624 } 4625 else 4626 { 4627 L[k][2][3][i]=list(); 4628 if (size(L[k][2][6][i])!=0) 4629 { 4630 if (size(L[k][2][8][i])!=0) 4631 { 4632 v=L[k][2][6][i],L[k][2][8][i]; 4633 L[k][2][7][i]=v; 4634 nr=nrows(L[k][2][1][i1]); 4635 nc=nrows(L[k][2][5][i1]); 4636 L[k][2][2][i]=concat(unitmat(nc),matrix(0,nr,nc)); 4637 L[k][2][4][i]=prodr(nr,nrows(L[k][2][5][i1])); 4638 } 4639 else 4640 { 4641 L[k][2][7][i]=L[k][2][6][i]; 4642 L[k][2][2][i]=unitmat(nrows(L[k][2][1][i1])); 4643 L[k][2][4][i]=list(); 4644 } 4645 } 4646 else 4647 { 4648 if (size(L[k][2][8][i])!=0) 4649 { 4650 L[k][2][7][i]=L[k][2][8][i]; 4651 L[k][2][2][i]=list(); 4652 L[k][2][4][i]=unitmat(nrows(L[k][2][5][i1])); 4653 } 4654 else 4655 { 4656 L[k][2][7][i]=list(); 4657 L[k][2][2][i]=list(); 4658 L[k][2][4][i]=list(); 4659 } 4660 } 4661 L[k+1][1][1][i]=L[k][2][5][i]; 4662 L[k+1][1][6][i]=L[k][2][8][i]; 4663 } 4664 } 4665 } 4666 i=d+size(L)+1; 4667 v=0; 4668 if (size(L[k][2][5][i1])!=0) 4669 { 4670 for (j=1; j<=nrows(L[k][2][5][i1]); j++) 4671 { 4672 mem=submat(L[k][2][5][i1],j,intvec(1..ncols(L[k][2][5][i1]))); 4673 v[j]=VdDeg(mem,d,L[k][2][8][i1]); 4674 } 4675 L[k][2][8][i]=v; 4676 if (size(L[k][2][6][i])!=0) 4677 { 4678 v=L[k][2][6][i],L[k][2][8][i]; 4679 L[k][2][7][i]=v; 4680 } 4681 else 4682 { 4683 L[k][2][7][i]=L[k][2][8][i]; 4684 } 4685 } 4686 else 4687 { 4688 L[k][2][8][i]=list(); 4689 L[k][2][7][i]=L[k][2][6][i]; 4690 } 4691 L[k+1][1][6][i]=L[k][2][8][i]; 4692 /* now we build V_dstrict resolutions for the sequences 4693 coker(L[k+1][1][1][1])>coker(L[k+1][1][3][1])>coker(L[k+1][1][5][i]) 4694 using the resolutions for coker(L[k][2][5][1]) we just obtained 4695 (works exactly the same as above)*/ 4696 for (i=2; i<=d+size(L); i++) 4697 { 4698 v=0; 4699 if (size(L[k+1][1][5][i1])!=0) 4700 { 4701 for (j=1; j<=nrows(L[k+1][1][5][i1]); j++) 4702 { 4703 i1=intvec(1..ncols(L[k+1][1][5][i1])); 4704 mem=submat(L[k+1][1][5][i1],j,i1); 4705 v[j]=VdDeg(mem,d,L[k+1][1][8][i1]); 4706 } 4707 L[k+1][1][8][i]=v; 4708 } 4709 else 4710 { 4711 L[k+1][1][8][i]=list(); 4712 } 4713 if (size(L[k+1][1][5][i])!=0) 4714 { 4715 if (size(L[k+1][1][1][i])!=0 or size(L[k+1][1][1][i1])!=0) 4716 { 4717 L[k+1][1][3][i]=transpose(syz(transpose(L[k+1][1][3][i1]))); 4718 nr=nrows(L[k+1][1][1][i1]); 4719 nc=ncols(L[k+1][1][5][i]); 4720 Pold=matrixLift(L[k+1][1][3][i]*prodr(nr,nc),L[k+1][1][5][i]); 4721 matrix Pi[1][ncols(L[k+1][1][3][i])]; 4722 for (l=1; l<=nrows(L[k+1][1][5][i]); l++) 4723 { 4724 for (j=1; j<=nrows(L[k+1][1][3][i]); j++) 4725 { 4726 i2=intvec(1..ncols(L[k+1][1][3][i])); 4727 Pi=Pi+Pold[l,j]*submat(L[k+1][1][3][i],j,i2); 4728 } 4729 if (l==1) 4730 { 4731 Picombined=transpose(Pi); 4732 } 4733 else 4734 { 4735 Picombined=concat(Picombined,transpose(Pi)); 4736 } 4737 Pi=0; 4738 } 4739 kill Pi; 4740 Picombined=transpose(Picombined); 4741 if(size(L[k+1][1][1][i])!=0) 4742 { 4743 if (i==2) 4744 { 4745 containsndeg=(0:ncols(L[k+1][1][1][i1])); 4746 } 4747 containsndeg=nDeg(L[k+1][1][1][i1],containsndeg); 4748 forhW=list(L[k+1][1][6][i], containsndeg); 4749 def HomWeyl=makeHomogenizedWeyl(n,forhW); 4750 setring HomWeyl; 4751 list L=fetch(B,L); 4752 matrix M=L[k+1][1][1][i]; 4753 module Mmod; 4754 list forM=nHomogenize(M,containsndeg,1); 4755 M=forM[1]; 4756 totaldeg=forM[2]; 4757 kill forM; 4758 matrix Maorig=fetch(B,Picombined); 4759 matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M))); 4760 Ma=nHomogenize(Ma,containsndeg); 4761 matrix mem,subm,zerom,subm2; 4762 matrix Pinew; 4763 M=transpose(M); 4764 SBcom=0; 4765 for (l=1; l<=nrows(Ma); l++) 4766 { 4767 i2=(ncols(Ma)+1..ncols(Maorig)); 4768 nc=ncols(Maorig)ncols(Ma); 4769 if (submat(Maorig,l,i2)==matrix(0,1,nc)) 4770 { 4771 for (cc=1; cc<=ncols(Ma); cc++) 4772 { 4773 Maorig[l,cc]=0; 4774 } 4775 } 4776 i1=(1..ncols(Ma)); 4777 i2=L[k+1][1][8][i]; 4778 subm=submat(Maorig,l,i1); 4779 subm2=submat(Maorig,l,(ncols(Ma)+1..ncols(Maorig))); 4780 if (VdDeg(subm,d,L[k+1][1][6][i])>VdDeg(subm2,d,i2) 4781 and subm!=matrix(0,1,ncols(Ma))) 4782 { 4783 //print("Reduzierung des VdGrades (Stelle2)"); 4784 if (SBcom==0) 4785 { 4786 Mmod=slimgb(M); 4787 M=Mmod; 4788 SBcom=1; 4789 } 4790 vd1=VdDeg(subm2,d,L[k+1][1][8][i]); 4791 mem=submat(Ma,l,(1..ncols(Ma))); 4792 mem=nHomogenize(mem,containsndeg); 4793 mem=h^totaldeg*mem; 4794 mem=transpose(mem); 4795 mem=reduce(mem,Mmod); 4796 if (l==1) 4797 { 4798 Pinew=mem; 4799 } 4800 else 4801 { 4802 Pinew=concat(Pinew,mem); 4803 } 4804 vd2=VdDeg(transpose(mem),d,L[k+1][1][6][i]); 4805 if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem))) 4806 {//should not happen 4807 //print("Reduzierung fehlgeschlagen!!!!(Stelle2)"); 4808 } 4809 } 4810 else 4811 { 4812 if (l==1) 4813 { 4814 Pinew=transpose(submat(Ma,l,(1..ncols(Ma)))); 4815 } 4816 else 4817 { 4818 subm=transpose(submat(Ma,l,(1..ncols(Ma)))); 4819 Pinew=concat(Pinew,subm); 4820 } 4821 } 4822 } 4823 Pinew=subst(Pinew,h,1); 4824 Pinew=transpose(Pinew); 4825 setring B; 4826 Pinew=fetch(HomWeyl,Pinew); 4827 kill HomWeyl; 4828 L[k+1][1][3][i]=concat(Pinew,L[k+1][1][5][i]); 4829 subm=transpose(L[k+1][1][1][i]); 4830 subm2=transpose(L[k+1][1][3][i]); 4831 L[k+1][1][3][i]=transpose(concat(subm,subm2)); 4832 } 4833 else 4834 { 4835 L[k+1][1][3][i]=Picombined; 4836 } 4837 L[k+1][2][1][i]=L[k+1][1][3][i]; 4838 nr=nrows(L[k+1][1][1][i1]); 4839 nc=ncols(L[k+1][1][5][i]); 4840 L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); 4841 L[k+1][1][4][i]=prodr(nr,nc); 4842 v=L[k+1][1][6][i],L[k+1][1][8][i]; 4843 L[k+1][1][7][i]=v; 4844 L[k+1][2][6][i]=L[k+1][1][7][i]; 4845 } 4846 else 4847 { 4848 L[k+1][1][3][i]=L[k+1][1][5][i]; 4849 L[k+1][1][2][i]=list(); 4850 L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i1])); 4851 L[k+1][1][7][i]=L[k+1][1][8][i]; 4852 L[k+1][2][6][i]=L[k+1][1][7][i]; 4853 L[k+1][2][1][i]=L[k+1][1][3][i]; 4854 } 4855 } 4856 else 4857 { 4858 if (size(L[k+1][1][1][i])!=0) 4859 { 4860 if (size(L[k+1][1][5][i1])!=0) 4861 { 4862 zerom=matrix(0,1,nrows(L[k+1][1][5][i1])); 4863 L[k+1][1][3][i]=concat(L[k+1][1][1][i],zerom); 4864 v=L[k+1][1][6][i],L[k+1][1][8][i]; 4865 L[k+1][1][7][i]=v; 4866 nr=nrows(L[k+1][1][1][i1]); 4867 nc=nrows(L[k+1][1][5][i1]); 4868 L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); 4869 L[k+1][1][4][i]=prodr(nr,nc); 4870 } 4871 else 4872 { 4873 L[k+1][1][3][i]=L[k+1][1][1][i]; 4874 &