[7fe9f8b] | 1 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 2 | version="$Id$"; |
---|
| 3 | category="Noncommutative"; |
---|
| 4 | info=" |
---|
| 5 | LIBRARY: derham.lib Computation of deRham cohomology |
---|
| 6 | |
---|
| 7 | AUTHORS: Cornelia Rottner, rottner@mathematik.uni-kl.de |
---|
| 8 | |
---|
| 9 | OVERVIEW: |
---|
| 10 | A library for computing the de Rham cohomology of complements of complex affine |
---|
| 11 | varieties. |
---|
| 12 | |
---|
| 13 | |
---|
| 14 | REFERENCES: |
---|
| 15 | [OT] Oaku, T.; Takayama, N.: Algorithms of D-modules - restriction, tensor product, |
---|
| 16 | localzation, and local cohomology groups, J. Pure Appl. Algebra 156, 267-308 |
---|
| 17 | (2001) |
---|
| 18 | [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, |
---|
| 21 | 303-321 (1999) |
---|
| 22 | [W2] Walther, U.: Algorithmic computation of de Rham Cohomology of Complements of |
---|
| 23 | Complex Affine Varieties, J. Symbolic Computation 29, 796-839 (2000) |
---|
| 24 | |
---|
| 25 | |
---|
| 26 | PROCEDURES: |
---|
| 27 | |
---|
| 28 | deRhamCohomology(list[,opt]); computes the de Rham cohomology |
---|
| 29 | MVComplex(list); computes the Mayer-Vietoris complex |
---|
| 30 | "; |
---|
| 31 | |
---|
| 32 | LIB "nctools.lib"; |
---|
| 33 | LIB "matrix.lib"; |
---|
| 34 | LIB "qhmoduli.lib"; |
---|
| 35 | LIB "general.lib"; |
---|
| 36 | LIB "dmod.lib"; |
---|
| 37 | LIB "bfun.lib"; |
---|
| 38 | LIB "dmodapp.lib"; |
---|
| 39 | LIB "poly.lib"; |
---|
| 40 | LIB "schreyer.lib"; |
---|
| 41 | |
---|
| 42 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 43 | |
---|
| 44 | proc deRhamCohomology(list L,list #) |
---|
| 45 | "USAGE: deRhamCohomology(L[,choices]); L a list consisting of polynomials, choices |
---|
| 46 | optional list consisting of one up to three strings @* |
---|
| 47 | The optional strings may be one of the strings@* |
---|
| 48 | -'Vdres': compute the Cartan-Eilenberg resolutions via V_d-homogenization |
---|
| 49 | and without using Schreyer's method @* |
---|
| 50 | -'Sres': compute the Cartan-Eilenberg resolutions in the homogenized Weyl |
---|
| 51 | algebra using Schreyer's method@* |
---|
| 52 | one of the strings@* |
---|
| 53 | -'iterativeloc': compute localizations by factorizing the polynomials and |
---|
| 54 | sucessive localization of the factors @* |
---|
| 55 | -'no iterativeloc': compute localizations by directly localizing the |
---|
| 56 | product@* |
---|
| 57 | and one of the strings |
---|
| 58 | -'onlybounds': computes bounds for the minimal and maximal interger roots |
---|
| 59 | of the global b-function |
---|
| 60 | -'exactroots' computes the minimal and maximal integer root of the global |
---|
| 61 | b-function |
---|
| 62 | The default is 'Sres', 'iterativeloc' and 'onlybounds'. |
---|
| 63 | ASSUME: -The basering must be a polynomial ring over the field of rational numbers@* |
---|
| 64 | RETURN: list, where the ith entry is the (i-1)st de Rham cohomology group of the |
---|
| 65 | complement of the complex affine variety given by the polynomials in L |
---|
| 66 | EXAMPLE:example deRhamCohomology; shows an example |
---|
| 67 | " |
---|
| 68 | { |
---|
| 69 | intvec saveoptions=option(get); |
---|
| 70 | intvec i1,i2; |
---|
| 71 | option(none); |
---|
| 72 | int recursiveloc=1; |
---|
| 73 | int i,j,nr,nc; |
---|
| 74 | def R=basering; |
---|
| 75 | poly islcm, forlcm; |
---|
| 76 | int n=nvars(R); |
---|
| 77 | int le=size(L)+n; |
---|
| 78 | string Syzstring="Sres"; |
---|
| 79 | int onlybounds=1; |
---|
| 80 | 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 | } |
---|
| 95 | for (i=1; i<=size(L); i++) |
---|
| 96 | { |
---|
| 97 | if (L[i]==0) |
---|
| 98 | { |
---|
| 99 | L=delete(L,i); |
---|
| 100 | i=i-1; |
---|
| 101 | } |
---|
| 102 | } |
---|
| 103 | if (size(L)==0) |
---|
| 104 | { |
---|
| 105 | return (list(0)); |
---|
| 106 | } |
---|
| 107 | for (i=1; i<= size(L); i++) |
---|
| 108 | { |
---|
| 109 | if (leadcoef(L[i])-L[i]==0) |
---|
| 110 | { |
---|
| 111 | return(list(1)); |
---|
| 112 | } |
---|
| 113 | } |
---|
| 114 | if (size(L)==0) |
---|
| 115 | { |
---|
| 116 | /*the complement of the variety given by the input is the whole space*/ |
---|
| 117 | return(list(1)); |
---|
| 118 | } |
---|
| 119 | 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 | } |
---|
| 127 | /* 1st step: compute the Mayer-Vietoris Complex and its Fourier transform*/ |
---|
| 128 | def W=MVComplex(L,recursiveloc);//new ring that contains the MV complex |
---|
| 129 | setring W; |
---|
| 130 | 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_d-strict free complex that is quasi-isomorphic to the |
---|
| 149 | complex fortoVdstrict |
---|
| 150 | The 1st entry of the list rem will be the quasi-isomorphic complex, the 2nd |
---|
| 151 | entry contains the cohomology modules and is needed for the computation of the |
---|
| 152 | global b-function*/ |
---|
| 153 | list rem=toVdStrictFreeComplex(fortoVdstrict,Syzstring); |
---|
| 154 | list newcomplex=rem[1]; |
---|
| 155 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 156 | /* 3rd step: Compute the bounds for the minimal and maximal integer root of the |
---|
| 157 | global b-function of newcomplex(i.e. compute the lcm of the b-functions of its |
---|
| 158 | cohomology modules)(if onlybouns=1). Else we compute the minimal and maximal |
---|
| 159 | integer root. |
---|
| 160 | |
---|
| 161 | If we compute only the bounds, we omit additional Groebner basis computations. |
---|
| 162 | However this leads to a higher-dimensional 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 | } |
---|
| 170 | else |
---|
| 171 | { |
---|
| 172 | list minmaxk=exactGlobalBFun(rem[2],Syzstring); |
---|
| 173 | } |
---|
| 174 | 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*i-2], |
---|
| 180 | m^i=newcomplex[3*i-1], 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]*/ |
---|
| 185 | int k,l; |
---|
| 186 | list truncatedcomplex,shorten,upto; |
---|
| 187 | for (i=1; i<=size(newcomplex) div 3; i++) |
---|
| 188 | { |
---|
| 189 | shorten[3*i-1]=list(); |
---|
| 190 | for (j=1; j<=size(newcomplex[3*i-1]); j++) |
---|
| 191 | { |
---|
| 192 | /*shorten[3*i-1][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*i-1][j]=list(minmaxk[1]-newcomplex[3*i-1][j]+1); |
---|
| 197 | shorten[3*i-1][j][2]=minmaxk[2]-newcomplex[3*i-1][j]+1; |
---|
| 198 | upto[size(upto)+1]=shorten[3*i-1][j][2]; |
---|
| 199 | if (shorten[3*i-1][j][2]<=0) |
---|
| 200 | { |
---|
| 201 | shorten[3*i-1][j]=list(); |
---|
| 202 | } |
---|
| 203 | else |
---|
| 204 | { |
---|
| 205 | if (shorten[3*i-1][j][1]<=0) |
---|
| 206 | { |
---|
| 207 | shorten[3*i-1][j][1]=1; |
---|
| 208 | } |
---|
| 209 | } |
---|
| 210 | } |
---|
| 211 | } |
---|
| 212 | int iupto=Max(upto);//maximal degree +1 of the polynomials we have to consider |
---|
| 213 | if (iupto<=0) |
---|
| 214 | { |
---|
| 215 | return(list(0)); |
---|
| 216 | } |
---|
| 217 | list allpolys; |
---|
| 218 | /*allpolys[i] will consist list of all monomials in D(1),...,D(n) of degree i-1*/ |
---|
| 219 | allpolys[1]=list(1); |
---|
| 220 | list minvar; |
---|
| 221 | minvar[1]=list(1); |
---|
| 222 | for (i=1; i<=iupto-1; 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 | } |
---|
| 235 | list keepformatrix,sizetruncom,fortrun,fst; |
---|
| 236 | int count,stc; |
---|
| 237 | intvec v,forin; |
---|
| 238 | matrix subm; |
---|
| 239 | /*now we compute the truncation*/ |
---|
| 240 | for (i=1; i<=size(newcomplex) div 3; i++) |
---|
| 241 | { |
---|
| 242 | /*truncatedcomplex[2*i-1] will contain all the generators for the truncation |
---|
| 243 | of D_n/(x(1),..,x(n))\otimes C[i]*/ |
---|
| 244 | truncatedcomplex[2*i-1]=list(); |
---|
| 245 | sizetruncom[2*i-1]=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*i-2]; j++) |
---|
| 253 | { |
---|
| 254 | if (size(shorten[3*i-1][j])!=0) |
---|
| 255 | { |
---|
| 256 | fortrun=sublist(allpolys,shorten[3*i-1][j][1],shorten[3*i-1][j][2]); |
---|
| 257 | truncatedcomplex[2*i-1][size(truncatedcomplex[2*i-1])+1]=fortrun[1]; |
---|
| 258 | count=count+fortrun[2]; |
---|
| 259 | fst=list(int(shorten[3*i-1][j][1])-1,int(shorten[3*i-1][j][2])-1); |
---|
| 260 | sizetruncom[2*i-1][size(sizetruncom[2*i-1])+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*(i-1)]); |
---|
| 279 | subm=submat(truncatedcomplex[2*(i-1)],i1,v); |
---|
| 280 | truncatedcomplex[2*(i-1)]=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*(i-1)]); |
---|
| 289 | truncatedcomplex[2*(i-1)]=matrix(0,nr,1); |
---|
| 290 | } |
---|
| 291 | } |
---|
| 292 | } |
---|
| 293 | matrix M; |
---|
| 294 | int st,pi,pj; |
---|
| 295 | poly ptc; |
---|
| 296 | int b,d,ideg,kplus,lplus; |
---|
| 297 | poly form,lform,nform; |
---|
| 298 | /*computation of the maps*/ |
---|
| 299 | 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*i-1]);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*i-1][k]); j++) |
---|
| 311 | { |
---|
| 312 | for (b=1; b<=size(truncatedcomplex[2*i-1][k][j]); b++) |
---|
| 313 | { |
---|
| 314 | form=truncatedcomplex[2*i-1][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++) |
---|
| 329 | { |
---|
| 330 | nc=2*(i+1)-1; |
---|
| 331 | ptc=truncatedcomplex[nc][l][ideg+1][d][1]; |
---|
| 332 | if (leadmonom(lform)==ptc) |
---|
| 333 | { |
---|
| 334 | nr=2*i-1; |
---|
| 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 | } |
---|
| 344 | } |
---|
| 345 | } |
---|
| 346 | } |
---|
| 347 | form=form-lform; |
---|
| 348 | } |
---|
| 349 | } |
---|
| 350 | } |
---|
| 351 | } |
---|
| 352 | } |
---|
| 353 | } |
---|
| 354 | truncatedcomplex[2*i]=M; |
---|
| 355 | truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])]; |
---|
| 356 | } |
---|
| 357 | truncatedcomplex[2*i-1]=sizetruncom[2*i][size(sizetruncom[2*i])]; |
---|
| 358 | if (truncatedcomplex[2*i-1]!=0) |
---|
| 359 | { |
---|
| 360 | truncatedcomplex[2*i]=matrix(0,truncatedcomplex[2*i-1],1); |
---|
| 361 | } |
---|
| 362 | setring R; |
---|
| 363 | list truncatedcomplex=imap(W,truncatedcomplex); |
---|
| 364 | /*computes the cohomology of the complex (D^i,d^i) given by truncatedcomplex, |
---|
| 365 | i.e. D^i=C^truncatedcomplex[2*i-1] and d^i=truncatedcomplex[2*i]*/ |
---|
| 366 | list derhamhom=findCohomology(truncatedcomplex,le); |
---|
| 367 | option(set,saveoptions); |
---|
| 368 | return (derhamhom); |
---|
| 369 | } |
---|
| 370 | |
---|
| 371 | example |
---|
| 372 | { "EXAMPLE:"; |
---|
| 373 | ring r = 0,(x,y,z),dp; |
---|
| 374 | list L=(xy,xz); |
---|
| 375 | deRhamCohomology(L); |
---|
| 376 | } |
---|
| 377 | |
---|
| 378 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 379 | //COMPUTATION OF THE MAYER-VIETORIS COMPLEX |
---|
| 380 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 381 | |
---|
| 382 | proc MVComplex(list L,list #) |
---|
| 383 | "USAGE:MVComplex(L); L a list of polynomials |
---|
| 384 | ASSUME: -Basering is a polynomial ring with n vwariables and rational coefficients |
---|
| 385 | -L is a list of non-constant polynomials |
---|
| 386 | RETURN: ring W: the nth Weyl algebra @* |
---|
| 387 | W contains a list MV, which represents the Mayer-Vietrois complex (C^i,d^i) of the |
---|
| 388 | polynomials contained in L as follows:@* |
---|
| 389 | the C^i are given by D_n^ncols(C[2*i-1])/im(C[2*i-1]) and the differentials |
---|
| 390 | d^i are given by C[2*i] |
---|
| 391 | EXAMPLE:example MVComplex; shows an example |
---|
| 392 | " |
---|
| 393 | { |
---|
| 394 | /* 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*/ |
---|
| 396 | def R=basering; |
---|
| 397 | int i; |
---|
| 398 | int iterative=1; |
---|
| 399 | if (size(#)!=0) |
---|
| 400 | { |
---|
| 401 | iterative=#[1]; |
---|
| 402 | } |
---|
| 403 | 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 non-constant"); |
---|
| 413 | return(); |
---|
| 414 | } |
---|
| 415 | } |
---|
| 416 | if (iterative==1) |
---|
| 417 | { |
---|
| 418 | /*compute the localizations by factorizing the polynomials and iterative |
---|
| 419 | localization of the factors */ |
---|
| 420 | for (i=1; i<=size(L); i++) |
---|
| 421 | { |
---|
| 422 | L[i]=factorize(L[i],1); |
---|
| 423 | } |
---|
| 424 | } |
---|
| 425 | int r=size(L); |
---|
| 426 | int n=nvars(basering); |
---|
| 427 | int le=size(L)+n; |
---|
| 428 | /*construct the ring Ws*/ |
---|
| 429 | def W=makeWeyl(n); |
---|
| 430 | setring W; |
---|
| 431 | list man=ringlist(W); |
---|
| 432 | 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 | } |
---|
| 443 | man[2][size(man[2])+1]="s";; |
---|
| 444 | man[3][3]=man[3][2]; |
---|
| 445 | man[3][2]=list("dp",intvec(1)); |
---|
| 446 | matrix N=UpOneMatrix(size(man[2])); |
---|
| 447 | man[5]=N; |
---|
| 448 | matrix M[1][1]; |
---|
| 449 | man[6]=transpose(concat(transpose(concat(man[6],M)),M)); |
---|
| 450 | def Ws=ring(man); |
---|
| 451 | setring Ws; |
---|
| 452 | int j,k,l,c; |
---|
| 453 | list L=fetch(R,L); |
---|
| 454 | list Cech; |
---|
| 455 | ideal J=var(1+n); |
---|
| 456 | for (i=2; i<=n; i++) |
---|
| 457 | { |
---|
| 458 | J=J,var(i+n); |
---|
| 459 | } |
---|
| 460 | Cech[1]=list(J); |
---|
| 461 | list Theta, remminroots; |
---|
| 462 | Theta[1]=list(list(list(),1,1)); |
---|
| 463 | list rem,findminintroot,diffmaps; |
---|
| 464 | int minroot,st,sk; |
---|
| 465 | intvec k1; |
---|
| 466 | poly fred,forfetch; |
---|
| 467 | matrix subm; |
---|
| 468 | int rmr; |
---|
| 469 | 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 s-parametric annihilator J(s) and the b-function |
---|
| 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 b-functions of L[i] computed above, |
---|
| 498 | because we want to plug in the same root r in all s-parametric |
---|
| 499 | annihilators we computed for L[i] ->this will ensure we can compute |
---|
| 500 | the maps of the MV complex*/ |
---|
| 501 | minroot=minIntRoot(findminintroot); |
---|
| 502 | for (j=1; j<=i; j++) |
---|
| 503 | { |
---|
| 504 | for (k=1; k<=k1[j]; k++) |
---|
| 505 | { |
---|
| 506 | sk=size(Cech[j+1])+1-k; |
---|
| 507 | Cech[j+1][size(Cech[j+1])+1-k]=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 | } |
---|
| 537 | 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 s-parametric 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 b-functions of L[i][c] |
---|
| 582 | computed above,because we want to plug in the same root r in all |
---|
| 583 | s-parametric annihilators we computed for L[i] ->this will |
---|
| 584 | ensure we can compute the maps of the MV complex*/ |
---|
| 585 | minroot=minIntRoot(findminintroot); |
---|
| 586 | for (j=1; j<=i; j++) |
---|
| 587 | { |
---|
| 588 | for (k=1; k<=k1[j]; k++) |
---|
| 589 | { |
---|
| 590 | st=size(Cech[j+1])+1-k; |
---|
| 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 | } |
---|
| 628 | setring W; |
---|
| 629 | /*map the modules and maps to the Weyl algebra*/ |
---|
| 630 | list diffmaps=imap(Ws,diffmaps); |
---|
| 631 | list Cechmodules=imap(Ws,Cech); |
---|
| 632 | list Cech; |
---|
| 633 | matrix sup; |
---|
| 634 | for (i=1; i<=r; i++) |
---|
| 635 | { |
---|
| 636 | sup=transpose(matrix(Cechmodules[i][1])); |
---|
| 637 | Cech[2*i-1]=sup; |
---|
| 638 | for (j=2; j<=size(Cechmodules[i]); j++) |
---|
| 639 | { |
---|
| 640 | sup=transpose(matrix(Cechmodules[i][j])); |
---|
| 641 | Cech[2*i-1]=dsum(Cech[2*i-1],sup); |
---|
| 642 | } |
---|
| 643 | sup=matrix(diffmaps[i]); |
---|
| 644 | Cech[2*i]=sup; |
---|
| 645 | } |
---|
| 646 | list MV=Cech; |
---|
| 647 | export MV; |
---|
| 648 | return (W); |
---|
| 649 | } |
---|
| 650 | |
---|
| 651 | example |
---|
| 652 | { "EXAMPLE:"; |
---|
| 653 | ring r = 0,(x,y,z),dp; |
---|
| 654 | list L=xy,xz; |
---|
| 655 | def C=MVComplex(L); |
---|
| 656 | setring C; |
---|
| 657 | MV; |
---|
| 658 | } |
---|
| 659 | |
---|
| 660 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 661 | |
---|
| 662 | static proc SannfsIBM(poly F,ideal myJ) |
---|
| 663 | "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 |
---|
| 665 | commutative variable@* |
---|
| 666 | f is a polynomial in the variables x(1),...,x(n) with rational coefficients |
---|
| 667 | @* |
---|
| 668 | J is holonomic and f-saturated |
---|
| 669 | RETURN AlList of the form (K,g), where K is an ideal and g a univariant polynomial |
---|
| 670 | in the variable s. K is the s-parametric annihilator of F and J and g is |
---|
| 671 | the b-function of F and J. |
---|
| 672 | " |
---|
| 673 | { |
---|
| 674 | /*modified version of the procedure SannfsBM from the library dmod.lib: SannfsBM |
---|
| 675 | computes the s-parametric annihilator for J=(x_1,...,x_n)*/ |
---|
| 676 | /* We use Algorithm 3.1.12 in[R] to compute the s-parametric |
---|
| 677 | annihilator. Then we use the s-parametric annihilator to compute the b-function |
---|
| 678 | via Algorithm 4.7 in [W1].*/ |
---|
| 679 | /* We assume that the basering the the nth Weyl algebra D_n. We create the ring |
---|
| 680 | D_n[s,t], where t*s=s*t-t*/ |
---|
| 681 | def save = basering; |
---|
| 682 | int N = nvars(basering)-1; |
---|
| 683 | int Nnew = N+2; |
---|
| 684 | int i,j; |
---|
| 685 | string s; |
---|
| 686 | list RL = ringlist(basering); |
---|
| 687 | list L, Lord; |
---|
| 688 | list tmp; |
---|
| 689 | intvec iv; |
---|
| 690 | L[1] = RL[1]; |
---|
| 691 | L[4] = RL[4]; |
---|
| 692 | list Name = RL[2]; |
---|
| 693 | Name=delete(Name,size(Name)); |
---|
| 694 | list RName; |
---|
| 695 | RName[1] = "t"; |
---|
| 696 | RName[2] = "s"; |
---|
| 697 | list DName; |
---|
| 698 | for(i=1;i<=N div 2;i++) |
---|
| 699 | { |
---|
| 700 | DName[i] = var(N div 2+i); |
---|
| 701 | Name=delete(Name,N div 2+1); |
---|
| 702 | } |
---|
| 703 | tmp[1] = "t"; |
---|
| 704 | tmp[2] = "s"; |
---|
| 705 | list NName = tmp +Name+DName; |
---|
| 706 | L[2] = NName; |
---|
| 707 | kill NName; |
---|
| 708 | tmp[1] = "lp"; |
---|
| 709 | iv = 1,1; |
---|
| 710 | tmp[2] = iv; |
---|
| 711 | Lord[1] = tmp; |
---|
| 712 | tmp[1] = "dp"; |
---|
| 713 | s = "iv="; |
---|
| 714 | for(i=1;i<=Nnew;i++) |
---|
| 715 | { |
---|
| 716 | s = s+"1,"; |
---|
| 717 | } |
---|
| 718 | s[size(s)]= ";"; |
---|
| 719 | execute(s); |
---|
| 720 | kill s; |
---|
| 721 | tmp[2] = iv; |
---|
| 722 | Lord[2] = tmp; |
---|
| 723 | tmp[1] = "C"; |
---|
| 724 | iv = 0; |
---|
| 725 | tmp[2] = iv; |
---|
| 726 | Lord[3] = tmp; |
---|
| 727 | tmp = 0; |
---|
| 728 | L[3] = Lord; |
---|
| 729 | def @R@ = ring(L); |
---|
| 730 | setring @R@; |
---|
| 731 | matrix @D[Nnew][Nnew]; |
---|
| 732 | @D[1,2]=t; |
---|
| 733 | for(i=1; i<=N div 2; i++) |
---|
| 734 | { |
---|
| 735 | @D[2+i, N div 2+2+i]=1; |
---|
| 736 | } |
---|
| 737 | def @R = nc_algebra(1,@D); |
---|
| 738 | setring @R; |
---|
| 739 | kill @R@; |
---|
| 740 | /*we start with the computation of the s-parametric annihilator*/ |
---|
| 741 | poly F = imap(save,F); |
---|
| 742 | ideal myJ=imap(save,myJ); |
---|
| 743 | for (i=1; i<=N div 2; i++) |
---|
| 744 | { |
---|
| 745 | myJ=subst(myJ,D(i),D(i)+diff(F,x(i))*t); |
---|
| 746 | } |
---|
| 747 | ideal I = t*F+s; |
---|
| 748 | I=I,myJ;//the s-parametric annihilator in D_n[s,t] |
---|
| 749 | /*we compute the intersection of I and D_n[s]*/ |
---|
| 750 | ideal J = slimgb(I); |
---|
| 751 | ideal K = nselect(J,1); |
---|
| 752 | K = slimgb(K);//the s-parametric annihilator |
---|
| 753 | /*we use K to compute the b-function*/ |
---|
| 754 | ideal B=K,F; |
---|
| 755 | B=slimgb(B); |
---|
| 756 | vector p=pIntersect(s,B); |
---|
| 757 | poly f=vec2poly(p,2); |
---|
| 758 | setring save; |
---|
| 759 | poly f=imap(@R,f); |
---|
| 760 | ideal K=imap(@R,K); |
---|
| 761 | return (list(K,f)); |
---|
| 762 | } |
---|
| 763 | |
---|
| 764 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 765 | //COMPUTATION OF QA UASI-ISOMORPHIC V_D-STRICT FREE COMPLEX |
---|
| 766 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 767 | |
---|
| 768 | 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 |
---|
| 771 | string, d an optional integer |
---|
| 772 | ASSUME: Basering is the Weyl algebra D_n @* |
---|
| 773 | (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) |
---|
| 775 | is generated by the rows of M_i. In particular it hold:@* |
---|
| 776 | - The M_i are m_i x r_i-matrices and the f_iare r_i x r_(i+1)-matrices @* |
---|
| 777 | -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 assumed |
---|
| 779 | to be n @* |
---|
| 780 | Syzstring is either: @* |
---|
| 781 | -'Sres' (computes the resolutions and Groebner bases in the homogenized |
---|
| 782 | Weyl algebra using Schreyer's method)@* |
---|
| 783 | or @* |
---|
| 784 | -'Vdres' (computes the resolutions via V_d-homogenization and without |
---|
| 785 | Schreyer's method)@* |
---|
| 786 | RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @* |
---|
| 787 | L_1 is of the form (i_(-n-1),g_(-n-1),m_(-n-1),...,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 |
---|
| 789 | of size i_j@* |
---|
| 790 | -D_n^(i_(-n-1))[m_(-n-1)]->...->D_n^(i_s)[m_s]->0 is a V_d-strict complex |
---|
| 791 | with differentials m_i that is quasi-isomorphic 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 |
---|
| 793 | the n_i are shift vectors such that:@* |
---|
| 794 | -coker(H_i) is the ith cohomology group of the complex given by L_1@* |
---|
| 795 | -the n_i are the shift vectors of the coker(H_i) |
---|
| 796 | THEORY: We follow Algorithm 3.8 in [W2] |
---|
| 797 | " |
---|
| 798 | { |
---|
| 799 | def B=basering; |
---|
| 800 | int n=nvars(B) div 2+2; |
---|
| 801 | int d=nvars(B) div 2; |
---|
| 802 | intvec v; |
---|
| 803 | list out, outall; |
---|
| 804 | int i,j,k,indi,nc,nr; |
---|
| 805 | matrix mem; |
---|
| 806 | intvec i1,i2; |
---|
| 807 | 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 | } |
---|
| 820 | /* If size(L)=2, our complex consists for only one non-trivial module. |
---|
| 821 | Therefore, we just have to compute a V_d-strict resolution of this module.*/ |
---|
| 822 | if (size(L)==2) |
---|
| 823 | { |
---|
| 824 | v=(0:ncols(L[1])); |
---|
| 825 | out[3*n-1]=v; |
---|
| 826 | out[3*n-2]=ncols(L[1]); |
---|
| 827 | out[3*n]=L[2]; |
---|
| 828 | if (Syzstring=="Vdres") |
---|
| 829 | { |
---|
| 830 | /*if Syzstring="Vdres", we compute a V_d-strict Groebner basis of L[1] |
---|
| 831 | using F-homogenization (Prop. 3.9 in [OT]); then we compute the syzygies |
---|
| 832 | and make them V_d-strict using Prop 3.9[OT] and so on*/ |
---|
| 833 | out[3*n-3]=VdStrictGB(L[1],d,v); |
---|
| 834 | for (i=n-1; i>=1; i--) |
---|
| 835 | { |
---|
| 836 | out[3*i-2]=nrows(out[3*i]); |
---|
| 837 | v=0; |
---|
| 838 | for (j=1; j<=out[3*i-2]; 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*i-1]=v; |
---|
| 844 | if (i!=1) |
---|
| 845 | { |
---|
| 846 | /*next step in the resolution*/ |
---|
| 847 | out[3*i-3]=transpose(syz(transpose(out[3*i]))); |
---|
| 848 | if (out[3*i-3]!=matrix(0,nrows(out[3*i-3]),ncols(out[3*i-3]))) |
---|
| 849 | { |
---|
| 850 | /*makes the resolution V_d-strict*/ |
---|
| 851 | out[3*i-3]=VdStrictGB(out[3*i-3],d,out[3*i-1]); |
---|
| 852 | } |
---|
| 853 | else |
---|
| 854 | { |
---|
| 855 | /*resolution is already computed*/ |
---|
| 856 | out[3*i-3]=matrix(0,1,ncols(out[3*i-3])); |
---|
| 857 | out[3*i-4]=intvec(0); |
---|
| 858 | out[3*i-5]=int(0); |
---|
| 859 | for (j=i-2; j>=1; j--) |
---|
| 860 | { |
---|
| 861 | out[3*j]=matrix(0,1,1); |
---|
| 862 | out[3*j-1]=intvec(0); |
---|
| 863 | out[3*j-2]=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*n-3]=L[1]; |
---|
| 880 | /*computes a ring with a list RES; RES is a V_d-strict 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*n-3]=transpose(matrix(RES[2])); |
---|
| 888 | out[3*n-3]=subst(out[3*n-3],h,1); |
---|
| 889 | for (i=n-1; i>=1; i--) |
---|
| 890 | { |
---|
| 891 | out[3*i-2]=nrows(out[3*i]); |
---|
| 892 | v=0; |
---|
| 893 | for (j=1; j<=out[3*i-2]; 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*i-1]=v;//shift vector such that the resolution RES is V_d-strict |
---|
| 899 | if (i!=1) |
---|
| 900 | { |
---|
| 901 | indi=0; |
---|
| 902 | if (size(RES)>=n-i+2) |
---|
| 903 | { |
---|
| 904 | nr=nrows(matrix(RES[n-i+2])); |
---|
| 905 | mem=matrix(0,nr,ncols(matrix(RES[n-i+2]))); |
---|
| 906 | if (matrix(RES[n-i+2])!=mem) |
---|
| 907 | { |
---|
| 908 | indi=1; |
---|
| 909 | out[3*i-3]=(matrix(RES[n-i+2])); |
---|
| 910 | if (nrows(out[3*i-3])-logens+1!=nrows(out[3*i])) |
---|
| 911 | { |
---|
| 912 | mem=out[3*i-3]; |
---|
| 913 | out[3*i-3]=matrix(mem,nrows(mem)+logens-1,ncols(mem)); |
---|
| 914 | } |
---|
| 915 | mem=out[3*i-3]; |
---|
| 916 | i1=intvec(logens..nrows(mem)); |
---|
| 917 | mem=submat(mem,i1,intvec(1..ncols(mem))); |
---|
| 918 | out[3*i-3]=transpose(mem); |
---|
| 919 | out[3*i-3]=subst(out[3*i-3],h,1); |
---|
| 920 | logens=logens+ncols(out[3*i-3]); |
---|
| 921 | } |
---|
| 922 | } |
---|
| 923 | if(indi==0) |
---|
| 924 | { |
---|
| 925 | out[3*i-3]=matrix(0,1,nrows(out[3*i])); |
---|
| 926 | out[3*i-4]=intvec(0); |
---|
| 927 | out[3*i-5]=int(0); |
---|
| 928 | for (j=i-2; j>=1; j--) |
---|
| 929 | { |
---|
| 930 | out[3*j]=matrix(0,1,1); |
---|
| 931 | out[3*j-1]=intvec(0); |
---|
| 932 | out[3*j-2]=int(0); |
---|
| 933 | } |
---|
| 934 | break; |
---|
| 935 | } |
---|
| 936 | } |
---|
| 937 | } |
---|
| 938 | setring B; |
---|
| 939 | out=fetch(ringofSyz,out);//contains the V_d-strict resolution |
---|
| 940 | kill ringofSyz; |
---|
| 941 | } |
---|
| 942 | outall[1]=out; |
---|
| 943 | outall[2]=list(list(out[3*n-3],out[3*n-1])); |
---|
| 944 | return(outall); |
---|
| 945 | } |
---|
| 946 | /*case size(L)>2: We compute a quasi-isomorphic 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); |
---|
| 953 | list rem; |
---|
| 954 | /* shortExactpiecesToVdStrict makes the sequences B^i->Z^i->H^i and |
---|
| 955 | Z^i->C^i->B^(i+1) V_d-strict*/ |
---|
| 956 | rem=shortExactPiecesToVdStrict(out,d,Syzstring); |
---|
| 957 | /*VdStrictDoubleComplexes computes V_d-strict resolutions over the seqeunces from |
---|
| 958 | proc shortExactPiecesToVdstrict*/ |
---|
| 959 | out=VdStrictDoubleComplexes(rem[1],d,Syzstring); |
---|
| 960 | 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 Cartan-Eilenberg resolution of |
---|
| 967 | (C^i,d^i)*/ |
---|
| 968 | out=assemblingDoubleComplexes(out); |
---|
| 969 | /*the proc totalComplex takes the total complex of the double complex from the |
---|
| 970 | proc assemblingDoubleComplexes*/ |
---|
| 971 | out=totalComplex(out); |
---|
| 972 | outall[1]=out; |
---|
| 973 | outall[2]=rem[2];//contains the cohomology groups and their shift vectors |
---|
| 974 | return (outall); |
---|
| 975 | } |
---|
| 976 | |
---|
| 977 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 978 | |
---|
| 979 | |
---|
| 980 | static proc sublist(list L,int m,int n) |
---|
| 981 | { |
---|
| 982 | list out; |
---|
| 983 | int i; int j; |
---|
| 984 | int count; |
---|
| 985 | 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 | } |
---|
| 994 | list o=list(out,count); |
---|
| 995 | return(o); |
---|
| 996 | } |
---|
| 997 | |
---|
| 998 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 999 | |
---|
| 1000 | static proc LMSubset(list L,list M) |
---|
| 1001 | { |
---|
| 1002 | int i; |
---|
| 1003 | int j=1; |
---|
| 1004 | list position=(M[size(M)],(-1)^(size(L))); |
---|
| 1005 | 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)^(i-1)); |
---|
| 1016 | j=j+1; |
---|
| 1017 | } |
---|
| 1018 | } |
---|
| 1019 | j=j+1; |
---|
| 1020 | |
---|
| 1021 | } |
---|
| 1022 | return (position); |
---|
| 1023 | } |
---|
| 1024 | |
---|
| 1025 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 1026 | |
---|
| 1027 | static proc shortExactPieces(list L) |
---|
| 1028 | { |
---|
| 1029 | /*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*/ |
---|
| 1032 | matrix Bnew= divdr(L[2],L[3]); |
---|
| 1033 | matrix Bold=Bnew; |
---|
| 1034 | matrix Z=divdr(Bnew,L[1]); |
---|
| 1035 | list bzh,zcb; |
---|
| 1036 | bzh=list(list(),list(),Z,unitmat(ncols(Z)),Z); |
---|
| 1037 | zcb=(Z, Bnew, L[1], unitmat(ncols(L[1])), Bnew); |
---|
| 1038 | list sep; |
---|
| 1039 | /* the list sep will be of size s such that |
---|
| 1040 | -sep[i]=(sep[i][1],sep[i][2]) is a list of two lists |
---|
| 1041 | -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)*/ |
---|
| 1045 | sep[1]=list(bzh,zcb); |
---|
| 1046 | int i; |
---|
| 1047 | list out; |
---|
| 1048 | 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[i-1] , L[i], L[i+1], L[i+2]); |
---|
| 1052 | sep[size(sep)+1]=out[1]; |
---|
| 1053 | Bold=out[2]; |
---|
| 1054 | } |
---|
| 1055 | bzh=(divdr(L[size(L)-2], L[size(L)-1]),L[size(L)-2], L[size(L)-1]); |
---|
| 1056 | bzh[4]=unitmat(ncols(L[size(L)-1])); |
---|
| 1057 | bzh[5]=transpose(concat(transpose(L[size(L)-2]),transpose(L[size(L)-1]))); |
---|
| 1058 | zcb=(L[size(L)-1], unitmat(ncols(L[size(L)-1])), L[size(L)-1],list(),list()); |
---|
| 1059 | sep[size(sep)+1]=list(bzh,zcb); |
---|
| 1060 | return(sep); |
---|
| 1061 | } |
---|
| 1062 | |
---|
| 1063 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 1064 | |
---|
| 1065 | static proc bzhzcb (matrix Bold,matrix f0,matrix C1,matrix f1,matrix C2) |
---|
| 1066 | { |
---|
| 1067 | matrix Bnew=divdr(f1,C2); |
---|
| 1068 | matrix Z= divdr(Bnew,C1); |
---|
| 1069 | matrix lift1= matrixLift(Bnew,f0); |
---|
| 1070 | matrix H=transpose(concat(transpose(lift1),transpose(Z))); |
---|
| 1071 | list bzh=(Bold, lift1, Z, unitmat(ncols(Z)),H); |
---|
| 1072 | list zcb=(Z, Bnew, C1, unitmat(ncols(C1)),Bnew); |
---|
| 1073 | list out=(list(bzh, zcb), Bnew); |
---|
| 1074 | return(out); |
---|
| 1075 | } |
---|
| 1076 | |
---|
| 1077 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 1078 | |
---|
| 1079 | static proc shortExactPiecesToVdStrict(list C,int d,list #) |
---|
| 1080 | {/* 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].*/ |
---|
| 1083 | /* If we compute our Groebner bases in the homogenized Weyl algebra, we already |
---|
| 1084 | compute some resolutions it omit additional Groebner basis computations later |
---|
| 1085 | on.*/ |
---|
| 1086 | int s =size(C);int i; int j; |
---|
| 1087 | string Syzstring="Sres"; |
---|
| 1088 | intvec v=0:ncols(C[s][1][5]); |
---|
| 1089 | 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 | } |
---|
| 1103 | list out; |
---|
| 1104 | list forout; |
---|
| 1105 | if (Syzstring=="Vdres") |
---|
| 1106 | { |
---|
| 1107 | out[s]=list(toVdStrictSequence(C[s][1],d,v, Syzstring,s)); |
---|
| 1108 | } |
---|
| 1109 | 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 resolutions |
---|
| 1120 | out[s][1][1][i+1]=resolutionofA[i]; |
---|
| 1121 | } |
---|
| 1122 | } |
---|
| 1123 | out[s][2]=list(list(out[s][1][3][1])); |
---|
| 1124 | out[s][2][2]=list(unitmat(ncols(out[s][1][3][1]))); |
---|
| 1125 | out[s][2][3]=list(out[s][1][3][1]); |
---|
| 1126 | out[s][2][4]=list(list()); |
---|
| 1127 | out[s][2][5]=list(list()); |
---|
| 1128 | out[s][2][6]=list(out[s][1][7][1]); |
---|
| 1129 | out[s][2][7]=list(out[s][2][6][1]); |
---|
| 1130 | out[s][2][8]=list(list()); |
---|
| 1131 | list resolutionofD; |
---|
| 1132 | list resolutionofF; |
---|
| 1133 | for (i=s-1; 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 resolutions |
---|
| 1140 | 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 | } |
---|
| 1158 | out[1]=list(list());//initalize our list |
---|
| 1159 | C[1][2][5]=out[2][1][1][1]; |
---|
| 1160 | /*Compute the last V_d-strict seqeunce*/ |
---|
| 1161 | if (Syzstring=="Vdres") |
---|
| 1162 | { |
---|
| 1163 | out[1][2]=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv"); |
---|
| 1164 | } |
---|
| 1165 | 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 | } |
---|
| 1180 | out[1][1][3]=list(out[1][2][1][1]); |
---|
| 1181 | out[1][1][5]=list(out[1][2][1][1]); |
---|
| 1182 | out[1][1][4]=list(unitmat(ncols(out[1][1][3][1]))); |
---|
| 1183 | out[1][1][7]=list(out[1][2][6][1]); |
---|
| 1184 | out[1][1][8]=list(out[1][2][6][1]); |
---|
| 1185 | out[1][1][1]=list(list()); |
---|
| 1186 | out[1][1][2]=list(list()); |
---|
| 1187 | out[1][1][6]=list(list()); |
---|
| 1188 | 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 | } |
---|
| 1196 | list Hi; |
---|
| 1197 | for (i=1; i<=size(out); i++) |
---|
| 1198 | { |
---|
| 1199 | Hi[i]=list(out[i][1][5][1],out[i][1][8][1]); |
---|
| 1200 | } |
---|
| 1201 | list outall; |
---|
| 1202 | outall[1]=out; |
---|
| 1203 | outall[2]=Hi; |
---|
| 1204 | return(outall); |
---|
| 1205 | } |
---|
| 1206 | |
---|
| 1207 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 1208 | |
---|
| 1209 | static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #) |
---|
| 1210 | { |
---|
| 1211 | /*this is the Algorithm 3.11 in [W2]*/ |
---|
| 1212 | int omitemptylist; |
---|
| 1213 | int lengthofres=si+n-1; |
---|
| 1214 | int i,j,logens; |
---|
| 1215 | def B=basering; |
---|
| 1216 | matrix bi=slimgb(transpose(C[5])); |
---|
| 1217 | /* Computation of a V_d-strict Groebner basis of C[5]: |
---|
| 1218 | -if Syzstring=="Vdres" this is done using the method of weighted homogenization |
---|
| 1219 | (Prop. 3.9 [OT]) |
---|
| 1220 | -else we use the homogenized Weyl algebra for Groebner basis computations |
---|
| 1221 | (Prop 9.9 [OT]), |
---|
| 1222 | in this case we already compute someresolutions (Thm. 9.10 [OT]) to omit |
---|
| 1223 | extra Groebner basis computations later on*/ |
---|
| 1224 | int nr,nc; |
---|
| 1225 | intvec i1,i2; |
---|
| 1226 | if (Syzstring=="Vdres") |
---|
| 1227 | { |
---|
| 1228 | if(size(#)==0) |
---|
| 1229 | { |
---|
| 1230 | matrix J_C=VdStrictGB(C[5],n,list(v)); |
---|
| 1231 | } |
---|
| 1232 | else |
---|
| 1233 | { |
---|
| 1234 | matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis |
---|
| 1235 | } |
---|
| 1236 | } |
---|
| 1237 | 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_d-strict resolution of C[5]->needed for proc |
---|
| 1247 | VdstrictDoubleComplexes*/ |
---|
| 1248 | def ringofSyz=Sres(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 C |
---|
| 1255 | 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[i-2]=(matrix(RES[i])); |
---|
| 1263 | |
---|
| 1264 | if (i==3) |
---|
| 1265 | { |
---|
| 1266 | if (nrows(rofC[i-2])-logens+1!=nrows(J_C)) |
---|
| 1267 | { |
---|
| 1268 | //build the resolution |
---|
| 1269 | nr=nrows(J_C)+logens-1; |
---|
| 1270 | nc=ncols(rofC[i-2]); |
---|
| 1271 | rofC[i-2]=matrix(rofC[i-2],nr,nc); |
---|
| 1272 | } |
---|
| 1273 | |
---|
| 1274 | } |
---|
| 1275 | if (i!=3) |
---|
| 1276 | { |
---|
| 1277 | if (nrows(rofC[i-2])-logens+1!=nrows(rofC[i-3])) |
---|
| 1278 | { |
---|
| 1279 | nr=nrows(rofC[i-3])+logens-1; |
---|
| 1280 | nc=ncols(rofC[i-2]); |
---|
| 1281 | rofC[i-2]=matrix(rofC[i-2],nr,nc); |
---|
| 1282 | } |
---|
| 1283 | } |
---|
| 1284 | i1=intvec(logens..nrows(rofC[i-2])); |
---|
| 1285 | i2=intvec(1..ncols(rofC[i-2])); |
---|
| 1286 | rofC[i-2]=transpose(submat(rofC[i-2],i1,i2)); |
---|
| 1287 | logens=logens+ncols(rofC[i-2]); |
---|
| 1288 | rofC[i-2]=subst(rofC[i-2],h,1); |
---|
| 1289 | } |
---|
| 1290 | else |
---|
| 1291 | { |
---|
| 1292 | rofC[i-2]=list(); |
---|
| 1293 | } |
---|
| 1294 | } |
---|
| 1295 | else |
---|
| 1296 | { |
---|
| 1297 | rofC[i-2]=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 | else |
---|
| 1315 | { |
---|
| 1316 | matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis |
---|
| 1317 | } |
---|
| 1318 | } |
---|
| 1319 | /* we compute a V_d-strict Groebner basis for C[3]*/ |
---|
| 1320 | matrix J_A=C[1]; |
---|
| 1321 | matrix f_CB=C[4]; |
---|
| 1322 | matrix f_ACB=transpose(concat(transpose(C[2]),transpose(f_CB))); |
---|
| 1323 | matrix J_AC=divdr(f_ACB,C[3]); |
---|
| 1324 | matrix P=matrixLift(J_AC * prodr(ncols(C[1]),ncols(C[5])) ,J_C); |
---|
| 1325 | list storePi; |
---|
| 1326 | matrix Pi[1][ncols(J_AC)]; |
---|
| 1327 | 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 | } |
---|
| 1336 | /*we compute the shift vector for C[1]*/ |
---|
| 1337 | intvec m_a; |
---|
| 1338 | list findMin; |
---|
| 1339 | int comMin; |
---|
| 1340 | 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=comMin-VdDeg(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 | else |
---|
| 1357 | { |
---|
| 1358 | m_a[i]=0; |
---|
| 1359 | } |
---|
| 1360 | } |
---|
| 1361 | matrix zero[ncols(J_A)][ncols(J_C)]; |
---|
| 1362 | matrix g_AB=concat(unitmat(ncols(J_A)),zero); |
---|
| 1363 | matrix g_BC= transpose(concat(transpose(zero),transpose(unitmat(ncols(J_C))))); |
---|
| 1364 | intvec m_b=m_a,v; |
---|
| 1365 | /* computation of a V_d-strict Groebner basis of C[1] (and resolution if |
---|
| 1366 | Syzstring=='Vdres') */ |
---|
| 1367 | if (Syzstring=="Vdres") |
---|
| 1368 | { |
---|
| 1369 | J_A=VdStrictGB(J_A,n,m_a); |
---|
| 1370 | } |
---|
| 1371 | 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(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[i-2]=matrix(RES[i]);// resolution for C[1] |
---|
| 1392 | if (i==3) |
---|
| 1393 | { |
---|
| 1394 | if (nrows(rofA[i-2])-logens+1!=nrows(J_A)) |
---|
| 1395 | { |
---|
| 1396 | nr=nrows(J_A)+logens-1; |
---|
| 1397 | nc=ncols(rofA[i-2]); |
---|
| 1398 | rofA[i-2]=matrix(rofA[i-2],nr,nc); |
---|
| 1399 | } |
---|
| 1400 | } |
---|
| 1401 | if (i!=3) |
---|
| 1402 | { |
---|
| 1403 | if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3])) |
---|
| 1404 | { |
---|
| 1405 | nr=nrows(rofA[i-3])+logens-1; |
---|
| 1406 | nc=ncols(rofA[i-2]); |
---|
| 1407 | rofA[i-2]=matrix(rofA[i-2],nr,nc); |
---|
| 1408 | } |
---|
| 1409 | } |
---|
| 1410 | i1=intvec(logens..nrows(rofA[i-2])); |
---|
| 1411 | i2=intvec(1..ncols(rofA[i-2])); |
---|
| 1412 | rofA[i-2]=transpose(submat(rofA[i-2],i1,i2)); |
---|
| 1413 | logens=logens+ncols(rofA[i-2]); |
---|
| 1414 | rofA[i-2]=subst(rofA[i-2],h,1); |
---|
| 1415 | } |
---|
| 1416 | else |
---|
| 1417 | { |
---|
| 1418 | rofA[i-2]=list(); |
---|
| 1419 | } |
---|
| 1420 | } |
---|
| 1421 | else |
---|
| 1422 | { |
---|
| 1423 | rofA[i-2]=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 | } |
---|
| 1440 | J_AC=transpose(storePi[1]); |
---|
| 1441 | for (i=2; i<= size(storePi); i++) |
---|
| 1442 | { |
---|
| 1443 | J_AC=concat(J_AC, transpose(storePi[i])); |
---|
| 1444 | } |
---|
| 1445 | J_AC=transpose(concat(transpose(matrix(J_A,nrows(J_A),nrows(J_AC))),J_AC)); |
---|
| 1446 | list Vdstrict=(list(J_A),list(g_AB),list(J_AC),list(g_BC),list(J_C),list(m_a)); |
---|
| 1447 | Vdstrict[7]=list(m_b); |
---|
| 1448 | Vdstrict[8]=list(v); |
---|
| 1449 | if(Syzstring=="Sres") |
---|
| 1450 | { |
---|
| 1451 | Vdstrict[9]=rofA; |
---|
| 1452 | if(size(#)==0) |
---|
| 1453 | { |
---|
| 1454 | Vdstrict[10]=rofC; |
---|
| 1455 | } |
---|
| 1456 | } |
---|
| 1457 | return (Vdstrict); |
---|
| 1458 | } |
---|
| 1459 | |
---|
| 1460 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 1461 | |
---|
| 1462 | static proc toVdStrictSequences (list L,int d,intvec v,string Syzstring,int sizeL) |
---|
| 1463 | { |
---|
| 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] such |
---|
| 1468 | that L[1][3]=L[2][1].We are going to transform them to V_d-strict sequences |
---|
| 1469 | J_D->J_A->J_F and J_A->J_B->J_C*/ |
---|
| 1470 | int omitemptylist; |
---|
| 1471 | int lengthofres=sizeL+d-1; |
---|
| 1472 | int logens; |
---|
| 1473 | def B=basering; |
---|
| 1474 | matrix J_F=L[1][5]; |
---|
| 1475 | matrix J_D=L[1][1]; |
---|
| 1476 | 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]) */ |
---|
| 1479 | matrix f_DFA=transpose(concat(transpose(L[1][2]),transpose(f_FA))); |
---|
| 1480 | matrix J_DF=divdr(f_DFA,L[1][3]);//coker(J_DF) is isomorphic to coker(L[2][1]); |
---|
| 1481 | matrix J_C=L[2][5]; |
---|
| 1482 | matrix f_CB=L[2][4]; |
---|
| 1483 | matrix f_DFCB=transpose(concat(transpose(f_DFA*L[2][2]),transpose(f_CB))); |
---|
| 1484 | 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*/ |
---|
| 1487 | matrix P=matrixLift(J_DFC*prodr(ncols(J_DF),ncols(L[2][5])),J_C); |
---|
| 1488 | list storePi; |
---|
| 1489 | matrix Pi[1][ncols(J_DFC)]; |
---|
| 1490 | int i; int j; |
---|
| 1491 | 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 | } |
---|
| 1500 | intvec m_a; |
---|
| 1501 | list findMin; |
---|
| 1502 | list noMin; |
---|
| 1503 | int comMin; |
---|
| 1504 | int nr,nc; |
---|
| 1505 | intvec i1,i2; |
---|
| 1506 | 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=comMin-VdDeg(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 | else |
---|
| 1524 | { |
---|
| 1525 | noMin[i]=1; |
---|
| 1526 | } |
---|
| 1527 | } |
---|
| 1528 | if (size(m_a) < ncols(J_DF)) |
---|
| 1529 | { |
---|
| 1530 | m_a[ncols(J_DF)]=0; |
---|
| 1531 | } |
---|
| 1532 | intvec m_f=m_a[ncols(J_D)+1..size(m_a)]; |
---|
| 1533 | /* Computation of a V_d-strict Groebner basis of J_F=L[1][5]: |
---|
| 1534 | if Syzstring=="Vdres" this is done using the method of weighted homogenization |
---|
| 1535 | (Prop. 3.9 [OT]) |
---|
| 1536 | else we use the homogenized Weyl algerba for Groebner basis computations |
---|
| 1537 | (Prop 9.9 [OT]), in this case we already compute resolutions |
---|
| 1538 | (Thm. 9.10 in [OT]) to omit extra Groebner basis computations later on*/ |
---|
| 1539 | if (Syzstring=="Vdres") |
---|
| 1540 | { |
---|
| 1541 | J_F=VdStrictGB(J_F,d,m_f); |
---|
| 1542 | } |
---|
| 1543 | 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(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[i-2]=(matrix(RES[i]));// resolution for J_F |
---|
| 1562 | if (i==3) |
---|
| 1563 | { |
---|
| 1564 | if (nrows(rofF[i-2])-logens+1!=nrows(J_F)) |
---|
| 1565 | { |
---|
| 1566 | nr=nrows(J_F)+logens-1; |
---|
| 1567 | nc=ncols(rofF[i-2]); |
---|
| 1568 | rofF[i-2]=matrix(rofF[i-2],nr,nc); |
---|
| 1569 | } |
---|
| 1570 | } |
---|
| 1571 | if (i!=3) |
---|
| 1572 | { |
---|
| 1573 | if (nrows(rofF[i-2])-logens+1!=nrows(rofF[i-3])) |
---|
| 1574 | { |
---|
| 1575 | nr=nrows(rofF[i-3])+logens-1; |
---|
| 1576 | rofF[i-2]=matrix(rofF[i-2],nr,ncols(rofF[i-2])); |
---|
| 1577 | } |
---|
| 1578 | } |
---|
| 1579 | i1=intvec(logens..nrows(rofF[i-2])); |
---|
| 1580 | i2=intvec(1..ncols(rofF[i-2])); |
---|
| 1581 | rofF[i-2]=transpose(submat(rofF[i-2],i1,i2)); |
---|
| 1582 | logens=logens+ncols(rofF[i-2]); |
---|
| 1583 | rofF[i-2]=subst(rofF[i-2],h,1); |
---|
| 1584 | } |
---|
| 1585 | else |
---|
| 1586 | { |
---|
| 1587 | rofF[i-2]=list(); |
---|
| 1588 | } |
---|
| 1589 | } |
---|
| 1590 | else |
---|
| 1591 | { |
---|
| 1592 | rofF[i-2]=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 | } |
---|
| 1609 | /*find shift vectors on the range of J_D*/ |
---|
| 1610 | P=matrixLift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F); |
---|
| 1611 | list storePinew; |
---|
| 1612 | matrix Pidf[1][ncols(J_DF)]; |
---|
| 1613 | 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 | } |
---|
| 1622 | intvec m_d; |
---|
| 1623 | 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=comMin-VdDeg(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 | } |
---|
| 1641 | 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_d-strict Groebner basis (and resolution of J_D if |
---|
| 1654 | Syzstring!='Vdres') for J_D*/ |
---|
| 1655 | if (Syzstring=="Vdres") |
---|
| 1656 | { |
---|
| 1657 | J_D=VdStrictGB(J_D,d,m_d); |
---|
| 1658 | } |
---|
| 1659 | 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(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[i-2]=(matrix(RES[i]));// resolution for J_D |
---|
| 1678 | if (i==3) |
---|
| 1679 | { |
---|
| 1680 | if (nrows(rofD[i-2])-logens+1!=nrows(J_D)) |
---|
| 1681 | { |
---|
| 1682 | nr=nrows(J_D)+logens-1; |
---|
| 1683 | rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2])); |
---|
| 1684 | } |
---|
| 1685 | } |
---|
| 1686 | if (i!=3) |
---|
| 1687 | { |
---|
| 1688 | if (nrows(rofD[i-2])-logens+1!=nrows(rofD[i-3])) |
---|
| 1689 | { |
---|
| 1690 | nr=nrows(rofD[i-3])+logens-1; |
---|
| 1691 | rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2])); |
---|
| 1692 | } |
---|
| 1693 | } |
---|
| 1694 | i1=intvec(logens..nrows(rofD[i-2])); |
---|
| 1695 | i2=intvec(1..ncols(rofD[i-2])); |
---|
| 1696 | rofD[i-2]=transpose(submat(rofD[i-2],i1,i2)); |
---|
| 1697 | logens=logens+ncols(rofD[i-2]); |
---|
| 1698 | rofD[i-2]=subst(rofD[i-2],h,1); |
---|
| 1699 | } |
---|
| 1700 | else |
---|
| 1701 | { |
---|
| 1702 | rofD[i-2]=list(); |
---|
| 1703 | } |
---|
| 1704 | } |
---|
| 1705 | else |
---|
| 1706 | { |
---|
| 1707 | rofD[i-2]=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_d-strict |
---|
| 1725 | Groebner basis and nrows(J_A)=nrows(J_D)+nrows(J_F) and |
---|
| 1726 | nrows(J_B)=nrows(J_A)+nrows(J_C)*/ |
---|
| 1727 | J_DF=transpose(storePinew[1]); |
---|
| 1728 | for (i=2; i<=nrows(J_F); i++) |
---|
| 1729 | { |
---|
| 1730 | J_DF=concat(J_DF,transpose(storePinew[i])); |
---|
| 1731 | } |
---|
| 1732 | J_DF=transpose(concat(transpose(matrix(J_D,nrows(J_D),nrows(J_DF))),J_DF)); |
---|
| 1733 | J_DFC=transpose(storePi[1]); |
---|
| 1734 | for (i=2; i<=nrows(J_C); i++) |
---|
| 1735 | { |
---|
| 1736 | J_DFC=concat(J_DFC,transpose(storePi[i])); |
---|
| 1737 | } |
---|
| 1738 | J_DFC=transpose(concat(transpose(matrix(J_DF,nrows(J_DF),nrows(J_DFC))),J_DFC)); |
---|
| 1739 | intvec m_b=m_a,v; |
---|
| 1740 | matrix zero[ncols(J_D)][ncols(J_F)]; |
---|
| 1741 | matrix g_DA=concat(unitmat(ncols(J_D)),zero); |
---|
| 1742 | matrix g_AF=transpose(concat(transpose(zero),unitmat(ncols(J_F)))); |
---|
| 1743 | matrix zero1[ncols(J_DF)][ncols(J_C)]; |
---|
| 1744 | matrix g_AB=concat(unitmat(ncols(J_DF)),zero1); |
---|
| 1745 | matrix g_BC=transpose(concat(transpose(zero1),unitmat(ncols(J_C)))); |
---|
| 1746 | list out; |
---|
| 1747 | out[1]=list(list(J_D),list(g_DA),list(J_DF),list(g_AF),list(J_F)); |
---|
| 1748 | out[1]=out[1]+list(list(m_d),list(m_a),list(m_f)); |
---|
| 1749 | out[2]=list(list(J_DF),list(g_AB),list(J_DFC),list(g_BC),list(J_C)); |
---|
| 1750 | out[2]=out[2]+list(list(m_a),list(m_b),list(v)); |
---|
| 1751 | if (Syzstring=="Sres") |
---|
| 1752 | { |
---|
| 1753 | out[3]=rofD; |
---|
| 1754 | out[4]=rofF; |
---|
| 1755 | } |
---|
| 1756 | return(out); |
---|
| 1757 | } |
---|
| 1758 | |
---|
| 1759 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 1760 | |
---|
| 1761 | static proc VdStrictDoubleComplexes(list L,int d,string Syzstring) |
---|
| 1762 | { |
---|
| 1763 | /* We compute V_d-strict resolutions over the V_d-strict short exact pieces from |
---|
| 1764 | the procedure shortExactPiecesToVdStrict. |
---|
| 1765 | We use Algorithms 3.14 and 3.15 in [W2]*/ |
---|
| 1766 | int i,k,c,j,l,totaldeg,comparedegs,SBcom,verk; |
---|
| 1767 | intvec fordegs; |
---|
| 1768 | intvec n_b,i1,i2; |
---|
| 1769 | matrix rem,forML,subm,zerom,unitm,subm2; |
---|
| 1770 | matrix J_B; |
---|
| 1771 | list store; |
---|
| 1772 | int t=size(L)+d; |
---|
| 1773 | int vd1,vd2,nr,nc; |
---|
| 1774 | def B=basering; |
---|
| 1775 | int n=nvars(B) div 2; |
---|
| 1776 | intvec v; |
---|
| 1777 | list forhW; |
---|
| 1778 | 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)+d-1); k++) |
---|
| 1784 | { |
---|
| 1785 | L[1][1][1][k+1]=list(); |
---|
| 1786 | L[1][1][2][k+1]=list(); |
---|
| 1787 | L[1][1][6][k+1]=list(); |
---|
| 1788 | } |
---|
| 1789 | L[1][1][6][size(L)+d+1]=list(); |
---|
| 1790 | matrix mem; |
---|
| 1791 | for (i=2; i<=d+size(L)+1; i++) |
---|
| 1792 | {; |
---|
| 1793 | v=0; |
---|
| 1794 | if(size(L[1][1][3][i-1])!=0) |
---|
| 1795 | { |
---|
| 1796 | if(i!=d+size(L)+1) |
---|
| 1797 | { |
---|
| 1798 | /*horizontal differential*/ |
---|
| 1799 | L[1][1][4][i-1]=unitmat(nrows(L[1][1][3][i-1])); |
---|
| 1800 | } |
---|
| 1801 | for (j=1; j<=nrows(L[1][1][3][i-1]); j++) |
---|
| 1802 | { |
---|
| 1803 | mem=submat(L[1][1][3][i-1],j,intvec(1..ncols(L[1][1][3][i-1]))); |
---|
| 1804 | v[j]=VdDeg(mem,d,L[1][1][7][i-1]); |
---|
| 1805 | } |
---|
| 1806 | L[1][1][7][i]=v;//new shift vector |
---|
| 1807 | L[1][1][8][i]=v; |
---|
| 1808 | L[1][2][6][i]=v; |
---|
| 1809 | } |
---|
| 1810 | else |
---|
| 1811 | { |
---|
| 1812 | if (i!=d+size(L)+1) |
---|
| 1813 | { |
---|
| 1814 | L[1][1][4][i-1]=list(); |
---|
| 1815 | } |
---|
| 1816 | L[1][1][7][i]=list(); |
---|
| 1817 | L[1][1][8][i]=list(); |
---|
| 1818 | L[1][2][6][i]=list(); |
---|
| 1819 | } |
---|
| 1820 | } |
---|
| 1821 | if (size(L[1][1][3][d+size(L)])!=0) |
---|
| 1822 | { |
---|
| 1823 | /*horizontal differential*/ |
---|
| 1824 | L[1][1][4][d+size(L)]=unitmat(nrows(L[1][1][3][d+size(L)])); |
---|
| 1825 | } |
---|
| 1826 | else |
---|
| 1827 | { |
---|
| 1828 | L[1][1][4][d+size(L)]=list(); |
---|
| 1829 | } |
---|
| 1830 | for (k=1; k<size(L); k++) |
---|
| 1831 | { |
---|
| 1832 | /* We build a V_d-strict resolution for coker(L[k][2][1][1])-> |
---|
| 1833 | coker(L[k][2][3][1])->coker(L[k][2][5][1]) using the resolution |
---|
| 1834 | obtained for coker(L[k][1][3][1]). |
---|
| 1835 | L[k][2][i][j] will be the jth module in the resolution of L[k][2][i][1] |
---|
| 1836 | for i=1,3,5. |
---|
| 1837 | L[k][2][i+5][j] will be the jth shift vector in the resolution of |
---|
| 1838 | L[k][2][i][1](this holds also for the case Syzstring=="Vdres")*/ |
---|
| 1839 | for (i=2; i<=d+size(L); i++) |
---|
| 1840 | { |
---|
| 1841 | v=0; |
---|
| 1842 | if (size(L[k][2][5][i-1])!=0) |
---|
| 1843 | { |
---|
| 1844 | for (j=1; j<=nrows(L[k][2][5][i-1]); j++) |
---|
| 1845 | { |
---|
| 1846 | i1=intvec(1..ncols(L[k][2][5][i-1])); |
---|
| 1847 | mem=submat(L[k][2][5][i-1],j,i1); |
---|
| 1848 | v[j]=VdDeg(mem,d,L[k][2][8][i-1]); |
---|
| 1849 | } |
---|
| 1850 | /*next shift vector in th resolution of coker(L[k][2][5][1])*/ |
---|
| 1851 | L[k][2][8][i]=v; |
---|
| 1852 | } |
---|
| 1853 | else |
---|
| 1854 | { |
---|
| 1855 | L[k][2][8][i]=list(); |
---|
| 1856 | } |
---|
| 1857 | /* we build step by step a resolution for coker(L[k][2][5][1]) using |
---|
| 1858 | the resolutions of coker(L[k][2][1][1]) and coker(L[k][2][5][1])*/ |
---|
| 1859 | if (size(L[k][2][5][i])!=0) |
---|
| 1860 | { |
---|
| 1861 | if (size(L[k][2][1][i])!=0 or size(L[k][2][1][i-1])!=0) |
---|
| 1862 | { |
---|
| 1863 | L[k][2][3][i]=transpose(syz(transpose(L[k][2][3][i-1]))); |
---|
| 1864 | nr= nrows(L[k][2][1][i-1]); |
---|
| 1865 | nc=ncols(L[k][2][5][i]); |
---|
| 1866 | Pold=matrixLift(L[k][2][3][i]*prodr(nr,nc), L[k][2][5][i]); |
---|
| 1867 | matrix Pi[1][ncols(L[k][2][3][i])]; |
---|
| 1868 | for (l=1; l<=nrows(L[k][2][5][i]); l++) |
---|
| 1869 | { |
---|
| 1870 | for (j=1; j<=nrows(L[k][2][3][i]); j++) |
---|
| 1871 | { |
---|
| 1872 | i2=intvec(1..ncols(L[k][2][3][i])); |
---|
| 1873 | Pi=Pi+Pold[l,j]*submat(L[k][2][3][i],j,i2); |
---|
| 1874 | } |
---|
| 1875 | if (l==1) |
---|
| 1876 | { |
---|
| 1877 | Picombined=transpose(Pi); |
---|
| 1878 | } |
---|
| 1879 | else |
---|
| 1880 | { |
---|
| 1881 | Picombined=concat(Picombined,transpose(Pi)); |
---|
| 1882 | } |
---|
| 1883 | Pi=0; |
---|
| 1884 | } |
---|
| 1885 | kill Pi; |
---|
| 1886 | Picombined=transpose(Picombined); |
---|
| 1887 | if (size(L[k][2][1][i])!=0) |
---|
| 1888 | { |
---|
| 1889 | if (i==2) |
---|
| 1890 | { |
---|
| 1891 | containsndeg=(0:ncols(L[k][2][1][1])); |
---|
| 1892 | } |
---|
| 1893 | containsndeg=nDeg(L[k][2][1][i-1],containsndeg); |
---|
| 1894 | forhW=list(L[k][2][6][i],containsndeg); |
---|
| 1895 | def HomWeyl=makeHomogenizedWeyl(n,forhW); |
---|
| 1896 | setring HomWeyl; |
---|
| 1897 | list L=fetch(B,L); |
---|
| 1898 | matrix M=L[k][2][1][i]; |
---|
| 1899 | list forM=nHomogenize(M,containsndeg,1); |
---|
| 1900 | M=forM[1]; |
---|
| 1901 | totaldeg=forM[2]; |
---|
| 1902 | kill forM; |
---|
| 1903 | matrix Maorig=fetch(B,Picombined); |
---|
| 1904 | matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M))); |
---|
| 1905 | matrix mem,subm,zerom; |
---|
| 1906 | matrix Pinew; |
---|
| 1907 | M=transpose(M); |
---|
| 1908 | SBcom=0; |
---|
| 1909 | for (l=1; l<=nrows(Ma); l++) |
---|
| 1910 | { |
---|
| 1911 | zerom=matrix(0,1,(ncols(Maorig)-ncols(Ma))); |
---|
| 1912 | i1=(ncols(Ma)+1..ncols(Maorig)); |
---|
| 1913 | if (submat(Maorig,l,i1)==zerom) |
---|
| 1914 | { |
---|
| 1915 | for (cc=1; cc<=ncols(Ma); cc++) |
---|
| 1916 | { |
---|
| 1917 | Maorig[l,cc]=0; |
---|
| 1918 | } |
---|
| 1919 | } |
---|
| 1920 | i2=(ncols(Ma)+1..ncols(Maorig)); |
---|
| 1921 | i1=(1..ncols(Ma)); |
---|
| 1922 | if (VdDeg(submat(Maorig,l,i1),d,L[k][2][6][i])> |
---|
| 1923 | VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]) and |
---|
| 1924 | submat(Maorig,l,i1)!=matrix(0,1,ncols(Ma))) |
---|
| 1925 | { |
---|
| 1926 | /*V_d-Grad is to big--> we make it smaller using |
---|
| 1927 | Vdnormal form computations*/ |
---|
| 1928 | if (SBcom==0) |
---|
| 1929 | { |
---|
| 1930 | M=slimgb(M); |
---|
| 1931 | SBcom=1; |
---|
| 1932 | } |
---|
| 1933 | //print("Reduzierung des V_d-Grades(Stelle1)"); |
---|
| 1934 | i2=(ncols(Ma)+1..ncols(Maorig)); |
---|
| 1935 | vd1=VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]); |
---|
| 1936 | mem=submat(Ma,l,(1..ncols(Ma))); |
---|
| 1937 | mem=nHomogenize(mem,containsndeg); |
---|
| 1938 | mem=h^totaldeg*mem; |
---|
| 1939 | mem=transpose(mem); |
---|
| 1940 | mem=reduce(mem,M); |
---|
| 1941 | matrix jt=transpose(subst(mem,h,1)); |
---|
| 1942 | setring B; |
---|
| 1943 | matrix jt=fetch(HomWeyl,jt); |
---|
| 1944 | matrix need=fetch(HomWeyl,Maorig); |
---|
| 1945 | need=submat(need,l,(1..ncols(need))); |
---|
| 1946 | i1=L[k][2][6][i]; |
---|
| 1947 | i2=L[k][2][8][i]; |
---|
| 1948 | jt=VdNormalForm(need,L[k][2][1][i],d,i1,i2); |
---|
| 1949 | setring HomWeyl; |
---|
| 1950 | mem=fetch(B,jt); |
---|
| 1951 | mem=transpose(mem); |
---|
| 1952 | if (l==1) |
---|
| 1953 | { |
---|
| 1954 | Pinew=mem; |
---|
| 1955 | } |
---|
| 1956 | else |
---|
| 1957 | { |
---|
| 1958 | Pinew=concat(Pinew,mem); |
---|
| 1959 | } |
---|
| 1960 | vd2=VdDeg(transpose(mem),d,L[k][2][6][i]); |
---|
| 1961 | if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem))) |
---|
| 1962 | {//should not happen!! |
---|
| 1963 | //print("Reduzierung fehlgeschlagen!!(Stelle1)"); |
---|
| 1964 | } |
---|
| 1965 | } |
---|
| 1966 | else |
---|
| 1967 | { |
---|
| 1968 | if (l==1) |
---|
| 1969 | { |
---|
| 1970 | Pinew=transpose(submat(Ma,l,(1..ncols(Ma)))); |
---|
| 1971 | } |
---|
| 1972 | else |
---|
| 1973 | { |
---|
| 1974 | subm=transpose(submat(Ma,l,(1..ncols(Ma)))); |
---|
| 1975 | Pinew=concat(Pinew,subm); |
---|
| 1976 | } |
---|
| 1977 | } |
---|
| 1978 | } |
---|
| 1979 | Pinew=subst(Pinew,h,1); |
---|
| 1980 | Pinew=transpose(Pinew); |
---|
| 1981 | setring B; |
---|
| 1982 | Pinew=fetch(HomWeyl,Pinew); |
---|
| 1983 | kill HomWeyl; |
---|
| 1984 | L[k][2][3][i]=concat(Pinew,L[k][2][5][i]); |
---|
| 1985 | subm=transpose(L[k][2][3][i]); |
---|
| 1986 | subm=concat(transpose(L[k][2][1][i]),subm); |
---|
| 1987 | L[k][2][3][i]=transpose(subm); |
---|
| 1988 | } |
---|
| 1989 | else |
---|
| 1990 | { |
---|
| 1991 | L[k][2][3][i]=Picombined; |
---|
| 1992 | } |
---|
| 1993 | L[k+1][1][1][i]=L[k][2][5][i]; |
---|
| 1994 | nr=nrows(L[k][2][1][i-1]); |
---|
| 1995 | nc=ncols(L[k][2][5][i]); |
---|
| 1996 | L[k][2][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); |
---|
| 1997 | L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nc); |
---|
| 1998 | v=L[k][2][6][i],L[k][2][8][i]; |
---|
| 1999 | L[k][2][7][i]=v; |
---|
| 2000 | L[k+1][1][6][i]=L[k][2][8][i]; |
---|
| 2001 | } |
---|
| 2002 | else |
---|
| 2003 | { |
---|
| 2004 | L[k][2][3][i]=L[k][2][5][i]; |
---|
| 2005 | L[k][2][2][i]=list(); |
---|
| 2006 | L[k][2][7][i]=L[k][2][8][i]; |
---|
| 2007 | L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1])); |
---|
| 2008 | L[k+1][1][6][i]=L[k][2][8][i]; |
---|
| 2009 | L[k+1][1][1][i]=L[k][2][5][i]; |
---|
| 2010 | } |
---|
| 2011 | } |
---|
| 2012 | else |
---|
| 2013 | { |
---|
| 2014 | if (size(L[k][2][1][i])!=0) |
---|
| 2015 | { |
---|
| 2016 | if (size(L[k][2][5][i-1])!=0) |
---|
| 2017 | { |
---|
| 2018 | nr=nrows(L[k][2][5][i-1]); |
---|
| 2019 | L[k][2][3][i]=concat(L[k][2][1][i],matrix(0,1,nr)); |
---|
| 2020 | v=L[k][2][6][i],L[k][2][8][i]; |
---|
| 2021 | L[k][2][7][i]=v; |
---|
| 2022 | nc=nrows(L[k][2][1][i-1]); |
---|
| 2023 | L[k][2][2][i]=concat(unitmat(nc),matrix(0,nc,nr)); |
---|
| 2024 | L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nr); |
---|
| 2025 | } |
---|
| 2026 | else |
---|
| 2027 | { |
---|
| 2028 | L[k][2][3][i]=L[k][2][1][i]; |
---|
| 2029 | L[k][2][7][i]=L[k][2][6][i]; |
---|
| 2030 | L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1])); |
---|
| 2031 | L[k][2][4][i]=list(); |
---|
| 2032 | } |
---|
| 2033 | L[k+1][1][1][i]=L[k][2][5][i]; |
---|
| 2034 | L[k+1][1][6][i]=L[k][2][8][i]; |
---|
| 2035 | } |
---|
| 2036 | else |
---|
| 2037 | { |
---|
| 2038 | L[k][2][3][i]=list(); |
---|
| 2039 | if (size(L[k][2][6][i])!=0) |
---|
| 2040 | { |
---|
| 2041 | if (size(L[k][2][8][i])!=0) |
---|
| 2042 | { |
---|
| 2043 | v=L[k][2][6][i],L[k][2][8][i]; |
---|
| 2044 | L[k][2][7][i]=v; |
---|
| 2045 | nr=nrows(L[k][2][1][i-1]); |
---|
| 2046 | nc=nrows(L[k][2][5][i-1]); |
---|
| 2047 | L[k][2][2][i]=concat(unitmat(nc),matrix(0,nr,nc)); |
---|
| 2048 | L[k][2][4][i]=prodr(nr,nrows(L[k][2][5][i-1])); |
---|
| 2049 | } |
---|
| 2050 | else |
---|
| 2051 | { |
---|
| 2052 | L[k][2][7][i]=L[k][2][6][i]; |
---|
| 2053 | L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1])); |
---|
| 2054 | L[k][2][4][i]=list(); |
---|
| 2055 | } |
---|
| 2056 | } |
---|
| 2057 | else |
---|
| 2058 | { |
---|
| 2059 | if (size(L[k][2][8][i])!=0) |
---|
| 2060 | { |
---|
| 2061 | L[k][2][7][i]=L[k][2][8][i]; |
---|
| 2062 | L[k][2][2][i]=list(); |
---|
| 2063 | L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1])); |
---|
| 2064 | } |
---|
| 2065 | else |
---|
| 2066 | { |
---|
| 2067 | L[k][2][7][i]=list(); |
---|
| 2068 | L[k][2][2][i]=list(); |
---|
| 2069 | L[k][2][4][i]=list(); |
---|
| 2070 | } |
---|
| 2071 | } |
---|
| 2072 | L[k+1][1][1][i]=L[k][2][5][i]; |
---|
| 2073 | L[k+1][1][6][i]=L[k][2][8][i]; |
---|
| 2074 | } |
---|
| 2075 | } |
---|
| 2076 | } |
---|
| 2077 | i=d+size(L)+1; |
---|
| 2078 | v=0; |
---|
| 2079 | if (size(L[k][2][5][i-1])!=0) |
---|
| 2080 | { |
---|
| 2081 | for (j=1; j<=nrows(L[k][2][5][i-1]); j++) |
---|
| 2082 | { |
---|
| 2083 | mem=submat(L[k][2][5][i-1],j,intvec(1..ncols(L[k][2][5][i-1]))); |
---|
| 2084 | v[j]=VdDeg(mem,d,L[k][2][8][i-1]); |
---|
| 2085 | } |
---|
| 2086 | L[k][2][8][i]=v; |
---|
| 2087 | if (size(L[k][2][6][i])!=0) |
---|
| 2088 | { |
---|
| 2089 | v=L[k][2][6][i],L[k][2][8][i]; |
---|
| 2090 | L[k][2][7][i]=v; |
---|
| 2091 | } |
---|
| 2092 | else |
---|
| 2093 | { |
---|
| 2094 | L[k][2][7][i]=L[k][2][8][i]; |
---|
| 2095 | } |
---|
| 2096 | } |
---|
| 2097 | else |
---|
| 2098 | { |
---|
| 2099 | L[k][2][8][i]=list(); |
---|
| 2100 | L[k][2][7][i]=L[k][2][6][i]; |
---|
| 2101 | } |
---|
| 2102 | L[k+1][1][6][i]=L[k][2][8][i]; |
---|
| 2103 | /* now we build V_d-strict resolutions for the sequences |
---|
| 2104 | coker(L[k+1][1][1][1])->coker(L[k+1][1][3][1])->coker(L[k+1][1][5][i]) |
---|
| 2105 | using the resolutions for coker(L[k][2][5][1]) we just obtained |
---|
| 2106 | (works exactly the same as above)*/ |
---|
| 2107 | for (i=2; i<=d+size(L); i++) |
---|
| 2108 | { |
---|
| 2109 | v=0; |
---|
| 2110 | if (size(L[k+1][1][5][i-1])!=0) |
---|
| 2111 | { |
---|
| 2112 | for (j=1; j<=nrows(L[k+1][1][5][i-1]); j++) |
---|
| 2113 | { |
---|
| 2114 | i1=intvec(1..ncols(L[k+1][1][5][i-1])); |
---|
| 2115 | mem=submat(L[k+1][1][5][i-1],j,i1); |
---|
| 2116 | v[j]=VdDeg(mem,d,L[k+1][1][8][i-1]); |
---|
| 2117 | } |
---|
| 2118 | L[k+1][1][8][i]=v; |
---|
| 2119 | } |
---|
| 2120 | else |
---|
| 2121 | { |
---|
| 2122 | L[k+1][1][8][i]=list(); |
---|
| 2123 | } |
---|
| 2124 | if (size(L[k+1][1][5][i])!=0) |
---|
| 2125 | { |
---|
| 2126 | if (size(L[k+1][1][1][i])!=0 or size(L[k+1][1][1][i-1])!=0) |
---|
| 2127 | { |
---|
| 2128 | L[k+1][1][3][i]=transpose(syz(transpose(L[k+1][1][3][i-1]))); |
---|
| 2129 | nr=nrows(L[k+1][1][1][i-1]); |
---|
| 2130 | nc=ncols(L[k+1][1][5][i]); |
---|
| 2131 | Pold=matrixLift(L[k+1][1][3][i]*prodr(nr,nc),L[k+1][1][5][i]); |
---|
| 2132 | matrix Pi[1][ncols(L[k+1][1][3][i])]; |
---|
| 2133 | for (l=1; l<=nrows(L[k+1][1][5][i]); l++) |
---|
| 2134 | { |
---|
| 2135 | for (j=1; j<=nrows(L[k+1][1][3][i]); j++) |
---|
| 2136 | { |
---|
| 2137 | i2=intvec(1..ncols(L[k+1][1][3][i])); |
---|
| 2138 | Pi=Pi+Pold[l,j]*submat(L[k+1][1][3][i],j,i2); |
---|
| 2139 | } |
---|
| 2140 | if (l==1) |
---|
| 2141 | { |
---|
| 2142 | Picombined=transpose(Pi); |
---|
| 2143 | } |
---|
| 2144 | else |
---|
| 2145 | { |
---|
| 2146 | Picombined=concat(Picombined,transpose(Pi)); |
---|
| 2147 | } |
---|
| 2148 | Pi=0; |
---|
| 2149 | } |
---|
| 2150 | kill Pi; |
---|
| 2151 | Picombined=transpose(Picombined); |
---|
| 2152 | if(size(L[k+1][1][1][i])!=0) |
---|
| 2153 | { |
---|
| 2154 | if (i==2) |
---|
| 2155 | { |
---|
| 2156 | containsndeg=(0:ncols(L[k+1][1][1][i-1])); |
---|
| 2157 | } |
---|
| 2158 | containsndeg=nDeg(L[k+1][1][1][i-1],containsndeg); |
---|
| 2159 | forhW=list(L[k+1][1][6][i], containsndeg); |
---|
| 2160 | def HomWeyl=makeHomogenizedWeyl(n,forhW); |
---|
| 2161 | setring HomWeyl; |
---|
| 2162 | list L=fetch(B,L); |
---|
| 2163 | matrix M=L[k+1][1][1][i]; |
---|
| 2164 | list forM=nHomogenize(M,containsndeg,1); |
---|
| 2165 | M=forM[1]; |
---|
| 2166 | totaldeg=forM[2]; |
---|
| 2167 | kill forM; |
---|
| 2168 | matrix Maorig=fetch(B,Picombined); |
---|
| 2169 | matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M))); |
---|
| 2170 | Ma=nHomogenize(Ma,containsndeg); |
---|
| 2171 | matrix mem,subm,zerom,subm2; |
---|
| 2172 | matrix Pinew; |
---|
| 2173 | M=transpose(M); |
---|
| 2174 | SBcom=0; |
---|
| 2175 | for (l=1; l<=nrows(Ma); l++) |
---|
| 2176 | { |
---|
| 2177 | i2=(ncols(Ma)+1..ncols(Maorig)); |
---|
| 2178 | nc=ncols(Maorig)-ncols(Ma); |
---|
| 2179 | if (submat(Maorig,l,i2)==matrix(0,1,nc)) |
---|
| 2180 | { |
---|
| 2181 | for (cc=1; cc<=ncols(Ma); cc++) |
---|
| 2182 | { |
---|
| 2183 | Maorig[l,cc]=0; |
---|
| 2184 | } |
---|
| 2185 | } |
---|
| 2186 | i1=(1..ncols(Ma)); |
---|
| 2187 | i2=L[k+1][1][8][i]; |
---|
| 2188 | subm=submat(Maorig,l,i1); |
---|
| 2189 | subm2=submat(Maorig,l,(ncols(Ma)+1..ncols(Maorig))); |
---|
| 2190 | if (VdDeg(subm,d,L[k+1][1][6][i])>VdDeg(subm2,d,i2) |
---|
| 2191 | and subm!=matrix(0,1,ncols(Ma))) |
---|
| 2192 | { |
---|
| 2193 | //print("Reduzierung des Vd-Grades (Stelle2)"); |
---|
| 2194 | if (SBcom==0) |
---|
| 2195 | { |
---|
| 2196 | M=slimgb(M); |
---|
| 2197 | SBcom=1; |
---|
| 2198 | } |
---|
| 2199 | vd1=VdDeg(subm2,d,L[k+1][1][8][i]); |
---|
| 2200 | mem=submat(Ma,l,(1..ncols(Ma))); |
---|
| 2201 | mem=nHomogenize(mem,containsndeg); |
---|
| 2202 | mem=h^totaldeg*mem; |
---|
| 2203 | mem=transpose(mem); |
---|
| 2204 | mem=reduce(mem,M); |
---|
| 2205 | if (l==1) |
---|
| 2206 | { |
---|
| 2207 | Pinew=mem; |
---|
| 2208 | } |
---|
| 2209 | else |
---|
| 2210 | { |
---|
| 2211 | Pinew=concat(Pinew,mem); |
---|
| 2212 | } |
---|
| 2213 | vd2=VdDeg(transpose(mem),d,L[k+1][1][6][i]); |
---|
| 2214 | if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem))) |
---|
| 2215 | {//should not happen |
---|
| 2216 | //print("Reduzierung fehlgeschlagen!!!!(Stelle2)"); |
---|
| 2217 | } |
---|
| 2218 | } |
---|
| 2219 | else |
---|
| 2220 | { |
---|
| 2221 | if (l==1) |
---|
| 2222 | { |
---|
| 2223 | Pinew=transpose(submat(Ma,l,(1..ncols(Ma)))); |
---|
| 2224 | } |
---|
| 2225 | else |
---|
| 2226 | { |
---|
| 2227 | subm=transpose(submat(Ma,l,(1..ncols(Ma)))); |
---|
| 2228 | Pinew=concat(Pinew,subm); |
---|
| 2229 | } |
---|
| 2230 | } |
---|
| 2231 | } |
---|
| 2232 | Pinew=subst(Pinew,h,1); |
---|
| 2233 | Pinew=transpose(Pinew); |
---|
| 2234 | setring B; |
---|
| 2235 | Pinew=fetch(HomWeyl,Pinew); |
---|
| 2236 | kill HomWeyl; |
---|
| 2237 | L[k+1][1][3][i]=concat(Pinew,L[k+1][1][5][i]); |
---|
| 2238 | subm=transpose(L[k+1][1][1][i]); |
---|
| 2239 | subm2=transpose(L[k+1][1][3][i]); |
---|
| 2240 | L[k+1][1][3][i]=transpose(concat(subm,subm2)); |
---|
| 2241 | } |
---|
| 2242 | else |
---|
| 2243 | { |
---|
| 2244 | L[k+1][1][3][i]=Picombined; |
---|
| 2245 | } |
---|
| 2246 | L[k+1][2][1][i]=L[k+1][1][3][i]; |
---|
| 2247 | nr=nrows(L[k+1][1][1][i-1]); |
---|
| 2248 | nc=ncols(L[k+1][1][5][i]); |
---|
| 2249 | L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); |
---|
| 2250 | L[k+1][1][4][i]=prodr(nr,nc); |
---|
| 2251 | v=L[k+1][1][6][i],L[k+1][1][8][i]; |
---|
| 2252 | L[k+1][1][7][i]=v; |
---|
| 2253 | L[k+1][2][6][i]=L[k+1][1][7][i]; |
---|
| 2254 | } |
---|
| 2255 | else |
---|
| 2256 | { |
---|
| 2257 | L[k+1][1][3][i]=L[k+1][1][5][i]; |
---|
| 2258 | L[k+1][1][2][i]=list(); |
---|
| 2259 | L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i-1])); |
---|
| 2260 | L[k+1][1][7][i]=L[k+1][1][8][i]; |
---|
| 2261 | L[k+1][2][6][i]=L[k+1][1][7][i]; |
---|
| 2262 | L[k+1][2][1][i]=L[k+1][1][3][i]; |
---|
| 2263 | } |
---|
| 2264 | } |
---|
| 2265 | else |
---|
| 2266 | { |
---|
| 2267 | if (size(L[k+1][1][1][i])!=0) |
---|
| 2268 | { |
---|
| 2269 | if (size(L[k+1][1][5][i-1])!=0) |
---|
| 2270 | { |
---|
| 2271 | zerom=matrix(0,1,nrows(L[k+1][1][5][i-1])); |
---|
| 2272 | L[k+1][1][3][i]=concat(L[k+1][1][1][i],zerom); |
---|
| 2273 | v=L[k+1][1][6][i],L[k+1][1][8][i]; |
---|
| 2274 | L[k+1][1][7][i]=v; |
---|
| 2275 | nr=nrows(L[k+1][1][1][i-1]); |
---|
| 2276 | nc=nrows(L[k+1][1][5][i-1]); |
---|
| 2277 | L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); |
---|
| 2278 | L[k+1][1][4][i]=prodr(nr,nc); |
---|
| 2279 | } |
---|
| 2280 | else |
---|
| 2281 | { |
---|
| 2282 | L[k+1][1][3][i]=L[k+1][1][1][i]; |
---|
| 2283 | L[k+1][1][7][i]=L[k+1][1][6][i]; |
---|
| 2284 | L[k+1][1][2][i]=unitmat(nrows(L[k+1][1][1][i-1])); |
---|
| 2285 | L[k+1][1][4][i]=list(); |
---|
| 2286 | } |
---|
| 2287 | L[k+1][2][1][i]=L[k+1][1][3][i]; |
---|
| 2288 | L[k+1][2][6][i]=L[k+1][1][7][i]; |
---|
| 2289 | } |
---|
| 2290 | else |
---|
| 2291 | { |
---|
| 2292 | L[k+1][1][3][i]=list(); |
---|
| 2293 | if (size(L[k+1][1][6][i])!=0) |
---|
| 2294 | { |
---|
| 2295 | if (size(L[k+1][1][8][i])!=0) |
---|
| 2296 | { |
---|
| 2297 | v=L[k+1][1][6][i],L[k+1][1][8][i]; |
---|
| 2298 | L[k+1][1][7][i]=v; |
---|
| 2299 | nr=nrows(L[k+1][1][1][i-1]); |
---|
| 2300 | nc=nrows(L[k+1][1][5][i-1]); |
---|
| 2301 | L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc)); |
---|
| 2302 | L[k+1][1][4][i]=prodr(nr,nrows(L[k+1][1][5][i-1])); |
---|
| 2303 | } |
---|
| 2304 | else |
---|
| 2305 | { |
---|
| 2306 | L[k+1][1][7][i]=L[k+1][1][6][i]; |
---|
| 2307 | L[k+1][1][2][i]=unitmat(nrows(L[k+1][1][1][i-1])); |
---|
| 2308 | L[k+1][1][4][i]=list(); |
---|
| 2309 | } |
---|
| 2310 | } |
---|
| 2311 | else |
---|
| 2312 | { |
---|
| 2313 | if (size(L[k+1][1][8][i])!=0) |
---|
| 2314 | { |
---|
| 2315 | L[k+1][1][7][i]=L[k+1][1][8][i]; |
---|
| 2316 | L[k+1][1][2][i]=list(); |
---|
| 2317 | L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i-1])); |
---|
| 2318 | } |
---|
| 2319 | else |
---|
| 2320 | { |
---|
| 2321 | L[k+1][1][7][i]=list(); |
---|
| 2322 | L[k+1][1][2][i]=list(); |
---|
| 2323 | L[k+1][1][4][i]=list(); |
---|
| 2324 | } |
---|
| 2325 | } |
---|
| 2326 | |
---|
| 2327 | L[k+1][2][1][i]=L[k+1][1][3][i]; |
---|
| 2328 | L[k+1][2][6][i]=L[k+1][1][7][i]; |
---|
| 2329 | } |
---|
| 2330 | } |
---|
| 2331 | } |
---|
| 2332 | i=size(L)+d+1; |
---|
| 2333 | v=0; |
---|
| 2334 | if (size(L[k+1][1][5][i-1])!=0) |
---|
| 2335 | { |
---|
| 2336 | for (j=1; j<=nrows(L[k+1][1][5][i-1]); j++) |
---|
| 2337 | { |
---|
| 2338 | i1=intvec(1..ncols(L[k+1][1][5][i-1])); |
---|
| 2339 | mem=submat(L[k+1][1][5][i-1],j,i1); |
---|
| 2340 | v[j]=VdDeg(mem,d,L[k+1][1][8][i-1]); |
---|
| 2341 | } |
---|
| 2342 | L[k+1][1][8][i]=v; |
---|
| 2343 | if (size(L[k+1][1][6][i])!=0) |
---|
| 2344 | { |
---|
| 2345 | v=L[k+1][1][6][i],L[k+1][1][8][i]; |
---|
| 2346 | L[k+1][1][7][i]=v; |
---|
| 2347 | } |
---|
| 2348 | else |
---|
| 2349 | { |
---|
| 2350 | L[k+1][1][7][i]=L[k+1][1][8][i]; |
---|
| 2351 | } |
---|
| 2352 | } |
---|
| 2353 | else |
---|
| 2354 | { |
---|
| 2355 | L[k+1][1][8][i]=list(); |
---|
| 2356 | L[k+1][1][7][i]=L[k+1][1][8][i]; |
---|
| 2357 | } |
---|
| 2358 | L[k+1][2][6][i]=L[k+1][1][7][i]; |
---|
| 2359 | } |
---|
| 2360 | for (k=1; k<=(size(L)+d); k++) |
---|
| 2361 | { |
---|
| 2362 | L[size(L)][2][5][k]=list(); |
---|
| 2363 | L[size(L)][2][4][k]=list(); |
---|
| 2364 | L[size(L)][2][8][k]=list(); |
---|
| 2365 | L[size(L)][2][3][k]=L[size(L)][2][1][k]; |
---|
| 2366 | L[size(L)][2][7][k]=L[size(L)][2][6][k]; |
---|
| 2367 | } |
---|
| 2368 | L[size(L)][2][7][size(L)+d+1]=L[size(L)][2][6][size(L)+d+1]; |
---|
| 2369 | L[size(L)][2][8][size(L)+d+1]=list(); |
---|
| 2370 | /* building the resolution of the last short exact piece*/ |
---|
| 2371 | for (i=2; i<=d+size(L); i++) |
---|
| 2372 | { |
---|
| 2373 | v=0; |
---|
| 2374 | if(size(L[size(L)][2][1][i-1])!=0) |
---|
| 2375 | { |
---|
| 2376 | L[size(L)][2][2][i]=unitmat(nrows(L[size(L)][2][1][i-1])); |
---|
| 2377 | } |
---|
| 2378 | else |
---|
| 2379 | { |
---|
| 2380 | L[size(L)][2][2][i-1]=list(); |
---|
| 2381 | } |
---|
| 2382 | } |
---|
| 2383 | return(L); |
---|
| 2384 | } |
---|
| 2385 | /*case Syzstring=="Vdres"*/ |
---|
| 2386 | list forVd; |
---|
| 2387 | for (k=1; k<=(size(L)+d); k++)//????? |
---|
| 2388 | { |
---|
| 2389 | /* we compute a V_d-strict resolution for the first short exact piece*/ |
---|
| 2390 | L[1][1][1][k+1]=list(); |
---|
| 2391 | L[1][1][2][k+1]=list(); |
---|
| 2392 | L[1][1][6][k+1]=list(); |
---|
| 2393 | if (size(L[1][1][3][k])!=0) |
---|
| 2394 | { |
---|
| 2395 | for (i=1; i<=nrows(L[1][1][3][k]); i++) |
---|
| 2396 | { |
---|
| 2397 | rem=submat(L[1][1][3][k],i,(1..ncols(L[1][1][3][k]))); |
---|
| 2398 | n_b[i]=VdDeg(rem,d,L[1][1][7][k]); |
---|
| 2399 | } |
---|
| 2400 | J_B=transpose(syz(transpose(L[1][1][3][k]))); |
---|
| 2401 | L[1][1][7][k+1]=n_b; |
---|
| 2402 | L[1][1][8][k+1]=n_b; |
---|
| 2403 | L[1][1][4][k+1]=unitmat(nrows(L[1][1][3][k])); |
---|
| 2404 | if (J_B!=matrix(0,nrows(J_B),ncols(J_B))) |
---|
| 2405 | { |
---|
| 2406 | J_B=VdStrictGB(J_B,d,n_b); |
---|
| 2407 | L[1][1][3][k+1]=J_B; |
---|
| 2408 | L[1][1][5][k+1]=J_B; |
---|
| 2409 | } |
---|
| 2410 | else |
---|
| 2411 | { |
---|
| 2412 | L[1][1][3][k+1]=list(); |
---|
| 2413 | L[1][1][5][k+1]=list(); |
---|
| 2414 | } |
---|
| 2415 | n_b=0; |
---|
| 2416 | } |
---|
| 2417 | else |
---|
| 2418 | { |
---|
| 2419 | L[1][1][3][k+1]=list(); |
---|
| 2420 | L[1][1][5][k+1]=list(); |
---|
| 2421 | L[1][1][7][k+1]=list(); |
---|
| 2422 | L[1][1][8][k+1]=list(); |
---|
| 2423 | L[1][1][4][k+1]=list(); |
---|
| 2424 | } |
---|
| 2425 | /* we compute step by step V_d-strict resolutions over |
---|
| 2426 | coker(L[i][2][1][1])->coker(L[i][2][3][1])->coker(L[i][2][1][5]) |
---|
| 2427 | and coker(L[i+1][1][1][1])->coker(L[i+1][1][3][1])->coker(L[i+1][1][1][5]) |
---|
| 2428 | using the already computed resolutions for coker(L[i][2][1][1])= |
---|
| 2429 | coker(L[i][1][3][1]) and coker(L[i+1][1][1][1])=coker(L[i][2][5][1])*/ |
---|
| 2430 | for (i=1; i<size(L); i++) |
---|
| 2431 | { |
---|
| 2432 | forVd[1]=L[i][2][1][k]; |
---|
| 2433 | forVd[2]=L[i][2][2][k]; |
---|
| 2434 | forVd[3]=L[i][2][3][k]; |
---|
| 2435 | forVd[4]=L[i][2][4][k]; |
---|
| 2436 | forVd[5]=L[i][2][5][k]; |
---|
| 2437 | forVd[6]=L[i][2][6][k]; |
---|
| 2438 | forVd[7]=L[i][2][7][k]; |
---|
| 2439 | forVd[8]=L[i][2][8][k]; |
---|
| 2440 | store=toVdStrict2x3Complex(forVd,d,L[i][1][3][k+1],L[i][1][7][k+1]); |
---|
| 2441 | for (j=1; j<=8; j++) |
---|
| 2442 | { |
---|
| 2443 | L[i][2][j][k+1]=store[j]; |
---|
| 2444 | } |
---|
| 2445 | forVd[1]=L[i+1][1][1][k]; |
---|
| 2446 | forVd[2]=L[i+1][1][2][k]; |
---|
| 2447 | forVd[3]=L[i+1][1][3][k]; |
---|
| 2448 | forVd[4]=L[i+1][1][4][k]; |
---|
| 2449 | forVd[5]=L[i+1][1][5][k]; |
---|
| 2450 | forVd[6]=L[i+1][1][6][k]; |
---|
| 2451 | forVd[7]=L[i+1][1][7][k]; |
---|
| 2452 | forVd[8]=L[i+1][1][8][k]; |
---|
| 2453 | store=toVdStrict2x3Complex(forVd,d,L[i][2][5][k+1],L[i][2][8][k+1]); |
---|
| 2454 | for (j=1; j<=8; j++) |
---|
| 2455 | { |
---|
| 2456 | L[i+1][1][j][k+1]=store[j]; |
---|
| 2457 | } |
---|
| 2458 | } |
---|
| 2459 | if (size(L[size(L)][1][7][k+1])!=0) |
---|
| 2460 | { |
---|
| 2461 | L[size(L)][2][4][k+1]=list(); |
---|
| 2462 | L[size(L)][2][5][k+1]=list(); |
---|
| 2463 | L[size(L)][2][6][k+1]=L[size(L)][1][7][k+1]; |
---|
| 2464 | L[size(L)][2][7][k+1]=L[size(L)][1][7][k+1]; |
---|
| 2465 | L[size(L)][2][8][k+1]=list(); |
---|
| 2466 | L[size(L)][2][2][k+1]=unitmat(size(L[size(L)][1][7][k+1])); |
---|
| 2467 | if (size(L[size(L)][1][3][k+1])!=0) |
---|
| 2468 | { |
---|
| 2469 | L[size(L)][2][1][k+1]=L[size(L)][1][3][k+1]; |
---|
| 2470 | L[size(L)][2][3][k+1]=L[size(L)][1][3][k+1]; |
---|
| 2471 | } |
---|
| 2472 | else |
---|
| 2473 | { |
---|
| 2474 | L[size(L)][2][1][k+1]=list(); |
---|
| 2475 | L[size(L)][2][3][k+1]=list(); |
---|
| 2476 | } |
---|
| 2477 | } |
---|
| 2478 | else |
---|
| 2479 | { |
---|
| 2480 | for (j=1; j<=8; j++) |
---|
| 2481 | { |
---|
| 2482 | L[size(L)][2][j][k+1]=list(); |
---|
| 2483 | } |
---|
| 2484 | } |
---|
| 2485 | } |
---|
| 2486 | k=t; |
---|
| 2487 | intvec n_c; |
---|
| 2488 | intvec vn_b; |
---|
| 2489 | list N_b; |
---|
| 2490 | int n; |
---|
| 2491 | /*computation of the shift vectors*/ |
---|
| 2492 | for (i=1; i<=size(L); i++) |
---|
| 2493 | { |
---|
| 2494 | for (n=1; n<=2; n++) |
---|
| 2495 | { |
---|
| 2496 | if (i==1 and n==1) |
---|
| 2497 | { |
---|
| 2498 | L[i][n][6][k+1]=list(); |
---|
| 2499 | } |
---|
| 2500 | else |
---|
| 2501 | { |
---|
| 2502 | if (n==1) |
---|
| 2503 | { |
---|
| 2504 | L[i][1][6][k+1]=L[i-1][2][8][k+1]; |
---|
| 2505 | } |
---|
| 2506 | else |
---|
| 2507 | { |
---|
| 2508 | L[i][2][6][k+1]=L[i][1][7][k+1]; |
---|
| 2509 | } |
---|
| 2510 | } |
---|
| 2511 | N_b[1]=L[i][n][6][k+1]; |
---|
| 2512 | if (size(L[i][n][5][k])!=0) |
---|
| 2513 | { |
---|
| 2514 | for (j=1; j<=nrows(L[i][n][5][k]); j++) |
---|
| 2515 | { |
---|
| 2516 | rem=submat(L[i][n][5][k],j,(1..ncols(L[i][n][5][k]))); |
---|
| 2517 | n_c[j]=VdDeg(rem,d,L[i][n][8][k]); |
---|
| 2518 | } |
---|
| 2519 | L[i][n][8][k+1]=n_c; |
---|
| 2520 | } |
---|
| 2521 | else |
---|
| 2522 | { |
---|
| 2523 | L[i][n][8][k+1]=list(); |
---|
| 2524 | } |
---|
| 2525 | N_b[2]=L[i][n][8][k+1]; |
---|
| 2526 | n_c=0; |
---|
| 2527 | if (size(N_b[1])!=0) |
---|
| 2528 | { |
---|
| 2529 | vn_b=N_b[1]; |
---|
| 2530 | if (size(N_b[2])!=0) |
---|
| 2531 | { |
---|
| 2532 | vn_b=vn_b,N_b[2]; |
---|
| 2533 | } |
---|
| 2534 | L[i][n][7][k+1]=vn_b; |
---|
| 2535 | } |
---|
| 2536 | else |
---|
| 2537 | { |
---|
| 2538 | if (size(N_b[2])!=0) |
---|
| 2539 | { |
---|
| 2540 | L[i][n][7][k+1]=N_b[2]; |
---|
| 2541 | } |
---|
| 2542 | else |
---|
| 2543 | { |
---|
| 2544 | L[i][n][7][k+1]=list(); |
---|
| 2545 | } |
---|
| 2546 | } |
---|
| 2547 | } |
---|
| 2548 | } |
---|
| 2549 | return(L); |
---|
| 2550 | } |
---|
| 2551 | |
---|
| 2552 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 2553 | |
---|
| 2554 | static proc toVdStrict2x3Complex(list L,int d,list #) |
---|
| 2555 | { |
---|
| 2556 | /* We build a one-step free resolution over a V_d-strict short exact piece |
---|
| 2557 | (Algorithm 3.14 in [W2]). |
---|
| 2558 | This procedure is called from the procedure VdStrictDoubleComplexes |
---|
| 2559 | if Syzstring=='Vdres'*/ |
---|
| 2560 | matrix rem; |
---|
| 2561 | int i,j,cc; |
---|
| 2562 | list J_A=list(list()); |
---|
| 2563 | list J_B=list(list()); |
---|
| 2564 | list J_C=list(list()); |
---|
| 2565 | list g_AB=list(list()); |
---|
| 2566 | list g_BC=list(list()); |
---|
| 2567 | list n_a=list(list()); |
---|
| 2568 | list n_b=list(list()); |
---|
| 2569 | list n_c=list(list()); |
---|
| 2570 | intvec n_b1; |
---|
| 2571 | matrix fromnf; |
---|
| 2572 | intvec i1,i2; |
---|
| 2573 | /* compute a one step V_d-strict resolution for L[5]*/ |
---|
| 2574 | if (size(L[5])!=0) |
---|
| 2575 | { |
---|
| 2576 | intvec n_c1; |
---|
| 2577 | for (i=1; i<=nrows(L[5]); i++) |
---|
| 2578 | { |
---|
| 2579 | rem=submat(L[5],i,intvec(1..ncols(L[5]))); |
---|
| 2580 | n_c1[i]=VdDeg(rem,d, L[8]);//new shift vector |
---|
| 2581 | } |
---|
| 2582 | n_c[1]=n_c1; |
---|
| 2583 | J_C[1]=transpose(syz(transpose(L[5]))); |
---|
| 2584 | if (J_C[1]!=matrix(0,nrows(J_C[1]),ncols(J_C[1]))) |
---|
| 2585 | { |
---|
| 2586 | J_C[1]=VdStrictGB(J_C[1],d,n_c1); |
---|
| 2587 | if (size(#[2])!=0)// new shift vector for the resolution of L[1] |
---|
| 2588 | { |
---|
| 2589 | n_a[1]=#[2]; |
---|
| 2590 | n_b1=n_a[1],n_c[1]; |
---|
| 2591 | n_b[1]=n_b1; |
---|
| 2592 | matrix zero[nrows(L[1])][nrows(L[5])]; |
---|
| 2593 | g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5]))); |
---|
| 2594 | if (size(#[1])!=0) |
---|
| 2595 | { |
---|
| 2596 | J_A=#[1];// one step V_d-strict resolution for L[1] |
---|
| 2597 | /* use resolutions of L[1] and L[5] to build a resolution for |
---|
| 2598 | L[3]*/ |
---|
| 2599 | J_B[1]=transpose(matrix(syz(transpose(L[3])))); |
---|
| 2600 | matrix P=matrixLift(J_B[1]*prodr(nrows(L[1]),nrows(L[5])),J_C[1]); |
---|
| 2601 | matrix Pi[1][ncols(J_B[1])]; |
---|
| 2602 | matrix Picombined; |
---|
| 2603 | for (i=1; i<=nrows(J_C[1]); i++) |
---|
| 2604 | { |
---|
| 2605 | for (j=1; j<=nrows(J_B[1]);j++) |
---|
| 2606 | { |
---|
| 2607 | Pi=Pi+P[i,j]*submat(J_B[1],j,intvec(1..ncols(J_B[1]))); |
---|
| 2608 | } |
---|
| 2609 | if (i==1) |
---|
| 2610 | { |
---|
| 2611 | Picombined=transpose(Pi); |
---|
| 2612 | } |
---|
| 2613 | else |
---|
| 2614 | { |
---|
| 2615 | Picombined=concat(Picombined,transpose(Pi)); |
---|
| 2616 | } |
---|
| 2617 | Pi=0; |
---|
| 2618 | } |
---|
| 2619 | Picombined=transpose(Picombined); |
---|
| 2620 | fromnf=VdNormalForm(Picombined,J_A[1],d,n_a[1],n_c[1]); |
---|
| 2621 | i1=intvec(1..nrows(Picombined)); |
---|
| 2622 | i2=intvec((ncols(J_A[1])+1)..ncols(Picombined)); |
---|
| 2623 | Picombined=concat(fromnf,submat(Picombined,i1,i2)); |
---|
| 2624 | J_B[1]=transpose(matrix(J_A[1],nrows(J_A[1]),ncols(J_B[1]))); |
---|
| 2625 | J_B[1]=transpose(concat(J_B[1],transpose(Picombined))); |
---|
| 2626 | g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5])))); |
---|
| 2627 | } |
---|
| 2628 | else//L[1] is already a resolution |
---|
| 2629 | { |
---|
| 2630 | //compute a resolution for L[3] |
---|
| 2631 | J_B=transpose(matrix(syz(transpose(L[3])))); |
---|
| 2632 | matrix P=matrixLift(J_B[1]*prodr(nrows(L[1]),nrows(L[5])),J_C[1]); |
---|
| 2633 | matrix Pi[1][ncols(J_B[1])]; |
---|
| 2634 | matrix Picombined; |
---|
| 2635 | for (i=1; i<=nrows(J_C[1]); i++) |
---|
| 2636 | { |
---|
| 2637 | for (j=1; j<=nrows(J_B[1]);j++) |
---|
| 2638 | { |
---|
| 2639 | Pi=Pi+P[i,j]*submat(J_B[1],j,intvec(1..ncols(J_B[1]))); |
---|
| 2640 | } |
---|
| 2641 | if (i==1) |
---|
| 2642 | { |
---|
| 2643 | Picombined=transpose(Pi); |
---|
| 2644 | } |
---|
| 2645 | else |
---|
| 2646 | { |
---|
| 2647 | Picombined=concat(Picombined,transpose(Pi)); |
---|
| 2648 | } |
---|
| 2649 | Pi=0; |
---|
| 2650 | } |
---|
| 2651 | Picombined=transpose(Picombined); |
---|
| 2652 | J_B[1]=Picombined; |
---|
| 2653 | g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5])))); |
---|
| 2654 | } |
---|
| 2655 | } |
---|
| 2656 | else |
---|
| 2657 | { |
---|
| 2658 | n_b=n_c[1]; |
---|
| 2659 | J_B[1]=J_C[1]; |
---|
| 2660 | g_BC=unitmat(ncols(J_C[1])); |
---|
| 2661 | } |
---|
| 2662 | } |
---|
| 2663 | else |
---|
| 2664 | { |
---|
| 2665 | J_C=list(list());// L[5] is already a resolution |
---|
| 2666 | if (size(#[2])!=0) |
---|
| 2667 | { |
---|
| 2668 | matrix zero[nrows(L[1])][nrows(L[5])]; |
---|
| 2669 | g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5])))); |
---|
| 2670 | n_a[1]=#[2]; |
---|
| 2671 | n_b1=n_a[1],n_c[1]; |
---|
| 2672 | n_b[1]=n_b1; |
---|
| 2673 | g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5]))); |
---|
| 2674 | if (size(#[1])!=0) |
---|
| 2675 | { |
---|
| 2676 | J_A=#[1]; |
---|
| 2677 | /*resolution of L[3]*/ |
---|
| 2678 | nr=nrows(J_A[1]); |
---|
| 2679 | J_B=concat(J_A[1],matrix(0,nr,nrows(L[3])-nrows(L[1]))); |
---|
| 2680 | } |
---|
| 2681 | } |
---|
| 2682 | else |
---|
| 2683 | { |
---|
| 2684 | n_b=n_c[1]; |
---|
| 2685 | g_BC=unitmat(ncols(L[5])); |
---|
| 2686 | } |
---|
| 2687 | } |
---|
| 2688 | } |
---|
| 2689 | else// L[5]=list(); |
---|
| 2690 | { |
---|
| 2691 | if (size(#[2])!=0) |
---|
| 2692 | { |
---|
| 2693 | n_a[1]=#[2]; |
---|
| 2694 | n_b=n_a[1]; |
---|
| 2695 | g_AB=unitmat(size(n_b[1])); |
---|
| 2696 | if (size(#[1])!=0) |
---|
| 2697 | { |
---|
| 2698 | J_A=#[1]; |
---|
| 2699 | J_B[1]=J_A[1];// resolution of L[3] equals that of L[1] |
---|
| 2700 | } |
---|
| 2701 | } |
---|
| 2702 | } |
---|
| 2703 | 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]); |
---|
| 2704 | return (out); |
---|
| 2705 | } |
---|
| 2706 | |
---|
| 2707 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 2708 | |
---|
| 2709 | static proc assemblingDoubleComplexes(list L) |
---|
| 2710 | { |
---|
| 2711 | /* The input is the output of VdStrictDoubleComplexes, we assemble the |
---|
| 2712 | resolutions of the L[i][2][3][1] to obtain a V_d-strict free Cartan-Eilenberg |
---|
| 2713 | resolution with modules P^i_j (1<=i<=size(L), j>=0) for the seqeunce |
---|
| 2714 | coker(L[1][2][3][1])->...->coker(L[size(L)][2][3][1])*/ |
---|
| 2715 | list out; |
---|
| 2716 | int i,j,k,l,oldj,newj,nr,nc; |
---|
| 2717 | for (i=1; i<=size(L); i++) |
---|
| 2718 | { |
---|
| 2719 | out[i]=list(list()); |
---|
| 2720 | out[i][1][1]=ncols(L[i][2][3][1]);//rank of module P^i_0 |
---|
| 2721 | if (size(L[i][2][5][1])!=0) |
---|
| 2722 | { |
---|
| 2723 | /*horizontal differential P^i_0->P^(i+1)_0*/ |
---|
| 2724 | nc=ncols(L[i][2][5][1]); |
---|
| 2725 | out[i][1][4]=prodr(ncols(L[i][2][3][1])-ncols(L[i][2][5][1]),nc); |
---|
| 2726 | } |
---|
| 2727 | else |
---|
| 2728 | { |
---|
| 2729 | /*horizontal differential P^i_0->0*/ |
---|
| 2730 | out[i][1][4]=matrix(0,ncols(L[i][2][3][1]),1); |
---|
| 2731 | } |
---|
| 2732 | oldj=newj; |
---|
| 2733 | for (j=1; j<=size(L[i][2][3]);j++) |
---|
| 2734 | { |
---|
| 2735 | out[i][j][2]=L[i][2][7][j];//shift vector of P^i_{j-1} |
---|
| 2736 | if (size(L[i][2][3][j])==0) |
---|
| 2737 | { |
---|
| 2738 | newj =j; |
---|
| 2739 | break; |
---|
| 2740 | } |
---|
| 2741 | out[i][j+1]=list(); |
---|
| 2742 | out[i][j+1][1]=nrows(L[i][2][3][j]);//rank of the module P^i_j |
---|
| 2743 | out[i][j+1][3]=L[i][2][3][j];//vertical differential P^i_j->P^(i+1)_j |
---|
| 2744 | if (size(L[i][2][5][j])!=0) |
---|
| 2745 | { |
---|
| 2746 | //horizonal differential P^i_j->P^(i-1)_j |
---|
| 2747 | nr=nrows(L[i][2][3][j])-nrows(L[i][2][5][j]); |
---|
| 2748 | out[i][j+1][4]=(-1)^j*prodr(nr,nrows(L[i][2][5][j])); |
---|
| 2749 | } |
---|
| 2750 | else |
---|
| 2751 | { |
---|
| 2752 | /*horizontal differential P^i_j->P^(i-1)_j*/ |
---|
| 2753 | out[i][j+1][4]=matrix(0,nrows(L[i][2][3][j]),1); |
---|
| 2754 | } |
---|
| 2755 | if(j==size(L[i][2][3])) |
---|
| 2756 | { |
---|
| 2757 | out[i][j+1][2]=L[i][2][7][j+1];//shift vector of P^i_j |
---|
| 2758 | newj=j+1; |
---|
| 2759 | } |
---|
| 2760 | } |
---|
| 2761 | if (i>1) |
---|
| 2762 | { |
---|
| 2763 | |
---|
| 2764 | for (k=1; k<=Min(list(oldj,newj)); k++) |
---|
| 2765 | { |
---|
| 2766 | /*horizonal differential P^(i-1)_(k-1)->P^i_(k-1)*/ |
---|
| 2767 | nr=nrows(out[i-1][k][4]); |
---|
| 2768 | out[i-1][k][4]=matrix(out[i-1][k][4],nr,out[i][k][1]); |
---|
| 2769 | } |
---|
| 2770 | for (k=newj+1; k<=oldj; k++) |
---|
| 2771 | { |
---|
| 2772 | /*no differential needed*/ |
---|
| 2773 | out[i-1][k]=delete(out[i-1][k],4); |
---|
| 2774 | } |
---|
| 2775 | } |
---|
| 2776 | } |
---|
| 2777 | return (out); |
---|
| 2778 | } |
---|
| 2779 | |
---|
| 2780 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 2781 | |
---|
| 2782 | static proc totalComplex(list L); |
---|
| 2783 | { |
---|
| 2784 | /* Input is the output of assemblingDoubleComplexes. |
---|
| 2785 | We obtain a complex C^1[m^1]->...->C^(r)[m^r] with differentials d^i and |
---|
| 2786 | shift vectors m^i (where C^r is placed in degree size(L)-1). |
---|
| 2787 | This complex is dercribed in the list out as follows: |
---|
| 2788 | rank(C^i)=out[3*i-2]; m_i=out[3*i-1] and d^i=out[3*i]*/ |
---|
| 2789 | list out;intvec rem1;intvec v; list remsize; int emp; |
---|
| 2790 | int i; int j; int c; int d; matrix M; int k; int l; |
---|
| 2791 | int n=nvars(basering) div 2; |
---|
| 2792 | list K; |
---|
| 2793 | for (i=1; i<=n+1; i++) |
---|
| 2794 | { |
---|
| 2795 | K[i]=list(); |
---|
| 2796 | } |
---|
| 2797 | L=K+L; |
---|
| 2798 | for (i=1; i<=size(L); i++) |
---|
| 2799 | { |
---|
| 2800 | emp=0; |
---|
| 2801 | if (size(L[i])!=0) |
---|
| 2802 | { |
---|
| 2803 | out[3*i-2]=L[i][1][1]; |
---|
| 2804 | v=L[i][1][1]; |
---|
| 2805 | rem1=L[i][1][2]; |
---|
| 2806 | emp=1; |
---|
| 2807 | } |
---|
| 2808 | else |
---|
| 2809 | { |
---|
| 2810 | out[3*i-2]=0; |
---|
| 2811 | v=0; |
---|
| 2812 | } |
---|
| 2813 | for (j=i+1; j<=size(L); j++) |
---|
| 2814 | { |
---|
| 2815 | if (size(L[j])>=j-i+1) |
---|
| 2816 | { |
---|
| 2817 | out[3*i-2]=out[3*i-2]+L[j][j-i+1][1]; |
---|
| 2818 | if (emp==0) |
---|
| 2819 | { |
---|
| 2820 | rem1=L[j][j-i+1][2]; |
---|
| 2821 | emp=1; |
---|
| 2822 | } |
---|
| 2823 | else |
---|
| 2824 | { |
---|
| 2825 | rem1=rem1,L[j][j-i+1][2]; |
---|
| 2826 | } |
---|
| 2827 | v[size(v)+1]=L[j][j-i+1][1]; |
---|
| 2828 | } |
---|
| 2829 | else |
---|
| 2830 | { |
---|
| 2831 | v[size(v)+1]=0; |
---|
| 2832 | } |
---|
| 2833 | } |
---|
| 2834 | out[3*i-1]=rem1; |
---|
| 2835 | v[size(v)+1]=0; |
---|
| 2836 | remsize[i]=v; |
---|
| 2837 | } |
---|
| 2838 | int o1; |
---|
| 2839 | int o2; |
---|
| 2840 | for (i=1; i<=size(L)-1; i++) |
---|
| 2841 | { |
---|
| 2842 | o1=1; |
---|
| 2843 | o2=1; |
---|
| 2844 | if (size(out[3*i-2])!=0) |
---|
| 2845 | { |
---|
| 2846 | o1=out[3*i-2]; |
---|
| 2847 | } |
---|
| 2848 | if (size(out[3*i+1])!=0) |
---|
| 2849 | { |
---|
| 2850 | o2=out[3*i+1]; |
---|
| 2851 | } |
---|
| 2852 | M=matrix(0,o1,o2); |
---|
| 2853 | if (size(L[i])!=0) |
---|
| 2854 | { |
---|
| 2855 | if (size(L[i][1][4])!=0) |
---|
| 2856 | { |
---|
| 2857 | M=matrix(L[i][1][4],o1,o2); |
---|
| 2858 | } |
---|
| 2859 | } |
---|
| 2860 | c=remsize[i][1]; |
---|
| 2861 | for (j=i+1; j<=size(L); j++) |
---|
| 2862 | { |
---|
| 2863 | if (remsize[i][j-i+1]!=0) |
---|
| 2864 | { |
---|
| 2865 | for (k=c+1; k<=c+remsize[i][j-i+1]; k++) |
---|
| 2866 | { |
---|
| 2867 | for (l=d+1; l<=d+remsize[i+1][j-i];l++) |
---|
| 2868 | { |
---|
| 2869 | M[k,l]=L[j][j-i+1][3][(k-c),(l-d)]; |
---|
| 2870 | } |
---|
| 2871 | } |
---|
| 2872 | d=d+remsize[i+1][j-i]; |
---|
| 2873 | if (remsize[i+1][j-i+1]!=0) |
---|
| 2874 | { |
---|
| 2875 | for (k=c+1; k<=c+remsize[i][j-i+1]; k++) |
---|
| 2876 | { |
---|
| 2877 | for (l=d+1; l<=d+remsize[i+1][j-i+1];l++) |
---|
| 2878 | { |
---|
| 2879 | M[k,l]=L[j][j-i+1][4][k-c,l-d]; |
---|
| 2880 | } |
---|
| 2881 | } |
---|
| 2882 | c=c+remsize[i][j-i+1]; |
---|
| 2883 | } |
---|
| 2884 | } |
---|
| 2885 | else |
---|
| 2886 | { |
---|
| 2887 | d=d+remsize[i+1][j-i]; |
---|
| 2888 | } |
---|
| 2889 | } |
---|
| 2890 | out[3*i]=M; |
---|
| 2891 | d=0; c=0; |
---|
| 2892 | } |
---|
| 2893 | out[3*size(L)]=matrix(0,out[3*size(L)-2],1); |
---|
| 2894 | return (out); |
---|
| 2895 | |
---|
| 2896 | } |
---|
| 2897 | |
---|
| 2898 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 2899 | //COMPUTATION OF THE BLOBAL B-FUNCTION |
---|
| 2900 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 2901 | |
---|
| 2902 | static proc globalBFun(list L,list #) |
---|
| 2903 | { |
---|
| 2904 | /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]), |
---|
| 2905 | where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an |
---|
| 2906 | intvec of size n_i. |
---|
| 2907 | We compute bounds for the minimal and maximal integer roots of the b-functions |
---|
| 2908 | of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def. |
---|
| 2909 | 6.1.1 in [R]) by combining Algorithm 6.1.6 in [R] and the method of principal |
---|
| 2910 | intersection (cf. Remark 6.1.7 in [R] 2012). |
---|
| 2911 | This works ONLY IF ALL B-FUNCTIONS ARE NON-ZERO, but this is the case since this |
---|
| 2912 | proc is only called from the procedure deRhamCohomology and the input comes |
---|
| 2913 | originally from the procedure toVdstrictFreeComplex*/ |
---|
| 2914 | if (size(#)==0)//# may contain the Syzstring |
---|
| 2915 | { |
---|
| 2916 | string Syzstring="Sres"; |
---|
| 2917 | } |
---|
| 2918 | else |
---|
| 2919 | { |
---|
| 2920 | string Syzstring=#[1]; |
---|
| 2921 | } |
---|
| 2922 | int i,j; |
---|
| 2923 | def W=basering; |
---|
| 2924 | int n=nvars(W) div 2; |
---|
| 2925 | list G0; |
---|
| 2926 | ideal I; |
---|
| 2927 | for (j=1; j<=size(L); j++) |
---|
| 2928 | { |
---|
| 2929 | G0[j]=list(); |
---|
| 2930 | for (i=1; i<=ncols(L[j][1]); i++) |
---|
| 2931 | { |
---|
| 2932 | G0[j][i]=I; |
---|
| 2933 | } |
---|
| 2934 | } |
---|
| 2935 | list out; |
---|
| 2936 | ideal I; poly f; |
---|
| 2937 | intvec i1; |
---|
| 2938 | for (j=1; j<=size(L); j++) |
---|
| 2939 | { |
---|
| 2940 | /*if the shift vector L[j][2] is non-zero we have to compute a V_d-strict |
---|
| 2941 | Groebner basis of L[j][1] with respect to the zero shift; otherwise L[i][1] |
---|
| 2942 | is already a V_d-strict Groebner basis, because it was obtained by the |
---|
| 2943 | procedure toVdStrictFreeComplex*/ |
---|
| 2944 | if (L[j][2]!=intvec(0:size(L[j][2]))) |
---|
| 2945 | { |
---|
| 2946 | if (Syzstring=="Vdres") |
---|
| 2947 | { |
---|
| 2948 | L[j][1]=VdStrictGB(L[j][1],n); |
---|
| 2949 | } |
---|
| 2950 | else |
---|
| 2951 | { |
---|
| 2952 | def HomWeyl=makeHomogenizedWeyl(n); |
---|
| 2953 | setring HomWeyl; |
---|
| 2954 | list L=fetch(W,L); |
---|
| 2955 | L[j][1]=nHomogenize(L[j][1]); |
---|
| 2956 | L[j][1]=transpose(matrix(slimgb(transpose(L[j][1])))); |
---|
| 2957 | L[j][1]=subst(L[j][1],h,1); |
---|
| 2958 | setring W; |
---|
| 2959 | L=fetch(HomWeyl,L); |
---|
| 2960 | kill HomWeyl; |
---|
| 2961 | } |
---|
| 2962 | } |
---|
| 2963 | for (i=1; i<=ncols(L[j][1]); i++) |
---|
| 2964 | { |
---|
| 2965 | G0[j][i]=I; |
---|
| 2966 | } |
---|
| 2967 | for (i=1; i<=nrows(L[j][1]); i++) |
---|
| 2968 | { |
---|
| 2969 | /*computes the terms of maximal V_d-degree of the biggest non-zero |
---|
| 2970 | component of submat(L[j][1],i,(1..ncols(L[j][1])))*/ |
---|
| 2971 | i1=(1..ncols(L[j][1])); |
---|
| 2972 | out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1); |
---|
| 2973 | f=L[j][1][i,out[2]]; |
---|
| 2974 | G0[j][out[2]]=G0[j][out[2]],f; |
---|
| 2975 | G0[j][out[2]]=compress(G0[j][out[2]]); |
---|
| 2976 | } |
---|
| 2977 | } |
---|
| 2978 | list save; |
---|
| 2979 | int l; |
---|
| 2980 | list weights; |
---|
| 2981 | /*bFctIdealModified computes the intersection of G0[j][i] and |
---|
| 2982 | x(1)D(1)+...+x(n)D(n) using the method of principal intersection*/ |
---|
| 2983 | for (j=1; j<=size(G0); j++) |
---|
| 2984 | { |
---|
| 2985 | for (i=1; i<=size(G0[j]); i++) |
---|
| 2986 | { |
---|
| 2987 | G0[j][i]=bFctIdealModified(G0[j][i]); |
---|
| 2988 | } |
---|
| 2989 | for (i=1; i<=size(G0[j]); i++) |
---|
| 2990 | { |
---|
| 2991 | weights=list(); |
---|
| 2992 | if (size(G0[j][i])!=0) |
---|
| 2993 | { |
---|
| 2994 | for (l=i; l<=size(G0[j]); l++) |
---|
| 2995 | { |
---|
| 2996 | if (size(G0[j][l])!=0) |
---|
| 2997 | { |
---|
| 2998 | weights[size(weights)+1]=L[j][2][l]; |
---|
| 2999 | } |
---|
| 3000 | } |
---|
| 3001 | G0[j][i]=list(G0[j][i][1]+Min(weights),G0[j][i][2]+Max(weights)); |
---|
| 3002 | } |
---|
| 3003 | } |
---|
| 3004 | } |
---|
| 3005 | list allmin; |
---|
| 3006 | list allmax; |
---|
| 3007 | for (j=1; j<=size(G0); j++) |
---|
| 3008 | { |
---|
| 3009 | for (i=1; i<=size(G0[j]); i++) |
---|
| 3010 | { |
---|
| 3011 | if (size(G0[j][i])!=0) |
---|
| 3012 | { |
---|
| 3013 | allmin[size(allmin)+1]=G0[j][i][1]; |
---|
| 3014 | allmax[size(allmax)+1]=G0[j][i][2]; |
---|
| 3015 | } |
---|
| 3016 | } |
---|
| 3017 | } |
---|
| 3018 | list minmax=list(Min(allmin),Max(allmax)); |
---|
| 3019 | return(minmax); |
---|
| 3020 | } |
---|
| 3021 | |
---|
| 3022 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3023 | |
---|
| 3024 | static proc exactGlobalBFun(list L,list #) |
---|
| 3025 | { |
---|
| 3026 | /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]), |
---|
| 3027 | where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an |
---|
| 3028 | intvec of size n_i. |
---|
| 3029 | We compute bounds for the minimal and maximal integer roots of the b-functions |
---|
| 3030 | of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def. |
---|
| 3031 | 6.1.1 in [R]) by combining Algorithm 6.1.6 in [R] and the method of principal |
---|
| 3032 | intersection (cf. Remark 6.1.7 in [R] 2012). |
---|
| 3033 | This works ONLY IF ALL B-FUNCTIONS ARE NON-ZERO, but this is the case since this |
---|
| 3034 | proc is only called from the procedure deRhamCohomology and the input comes |
---|
| 3035 | originally from the procedure toVdstrictFreeComplex*/ |
---|
| 3036 | if (size(#)==0)//# may contain the Syzstring |
---|
| 3037 | { |
---|
| 3038 | string Syzstring="Sres"; |
---|
| 3039 | } |
---|
| 3040 | else |
---|
| 3041 | { |
---|
| 3042 | string Syzstring=#[1]; |
---|
| 3043 | } |
---|
| 3044 | int i,j,k; |
---|
| 3045 | def W=basering; |
---|
| 3046 | int n=nvars(W) div 2; |
---|
| 3047 | list G0; |
---|
| 3048 | ideal I; |
---|
| 3049 | for (j=1; j<=size(L); j++) |
---|
| 3050 | { |
---|
| 3051 | G0[j]=list(); |
---|
| 3052 | for (i=1; i<=ncols(L[j][1]); i++) |
---|
| 3053 | { |
---|
| 3054 | G0[j][i]=I; |
---|
| 3055 | } |
---|
| 3056 | } |
---|
| 3057 | list out; |
---|
| 3058 | matrix M; |
---|
| 3059 | ideal I; poly f; |
---|
| 3060 | intvec i1; |
---|
| 3061 | for (j=1; j<=size(L); j++) |
---|
| 3062 | { |
---|
| 3063 | M=L[j][1]; |
---|
| 3064 | /*if the shift vector L[j][2] is non-zero we have to compute a V_d-strict |
---|
| 3065 | Groebner basis of L[j][1] with respect to the zero shift; otherwise L[i][1] |
---|
| 3066 | is already a V_d-strict Groebner basis, because it was obtained by the |
---|
| 3067 | procedure toVdStrictFreeComplex*/ |
---|
| 3068 | for (k=1; k<=ncols(L[j][1]); k++) |
---|
| 3069 | { |
---|
| 3070 | L[j][1]=permcol(M,1,k); |
---|
| 3071 | if (Syzstring=="Vdres") |
---|
| 3072 | { |
---|
| 3073 | L[j][1]=VdStrictGB(L[j][1],n); |
---|
| 3074 | } |
---|
| 3075 | else |
---|
| 3076 | { |
---|
| 3077 | def HomWeyl=makeHomogenizedWeyl(n); |
---|
| 3078 | setring HomWeyl; |
---|
| 3079 | list L=fetch(W,L); |
---|
| 3080 | L[j][1]=nHomogenize(L[j][1]); |
---|
| 3081 | L[j][1]=transpose(matrix(slimgb(transpose(L[j][1])))); |
---|
| 3082 | L[j][1]=subst(L[j][1],h,1); |
---|
| 3083 | setring W; |
---|
| 3084 | L=fetch(HomWeyl,L); |
---|
| 3085 | kill HomWeyl; |
---|
| 3086 | } |
---|
| 3087 | for (i=1; i<=nrows(L[j][1]); i++) |
---|
| 3088 | { |
---|
| 3089 | /*computes the terms of maximal V_d-degree of the biggest non-zero |
---|
| 3090 | component of submat(L[j][1],i,(1..ncols(L[j][1])))*/ |
---|
| 3091 | i1=(1..ncols(L[j][1])); |
---|
| 3092 | out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1); |
---|
| 3093 | f=L[j][1][i,out[2]]; |
---|
| 3094 | if (out[2]==1) |
---|
| 3095 | { |
---|
| 3096 | G0[j][k]=G0[j][k],f; |
---|
| 3097 | G0[j][k]=compress(G0[j][k]); |
---|
| 3098 | } |
---|
| 3099 | } |
---|
| 3100 | } |
---|
| 3101 | } |
---|
| 3102 | list save; |
---|
| 3103 | int l; |
---|
| 3104 | list weights; |
---|
| 3105 | /*bFctIdealModified computes the intersection of G0[j][i] and |
---|
| 3106 | x(1)D(1)+...+x(n)D(n) using the method of principal intersection*/ |
---|
| 3107 | for (j=1; j<=size(G0); j++) |
---|
| 3108 | { |
---|
| 3109 | for (i=1; i<=size(G0[j]); i++) |
---|
| 3110 | { |
---|
| 3111 | G0[j][i]=bFctIdealModified(G0[j][i]); |
---|
| 3112 | } |
---|
| 3113 | for (i=1; i<=size(G0[j]); i++) |
---|
| 3114 | { |
---|
| 3115 | if (size(G0[j][i])!=0) |
---|
| 3116 | { |
---|
| 3117 | G0[j][i]=list(G0[j][i][1]+L[j][2][i],G0[j][i][2]+L[j][2][i]); |
---|
| 3118 | } |
---|
| 3119 | } |
---|
| 3120 | } |
---|
| 3121 | list allmin; |
---|
| 3122 | list allmax; |
---|
| 3123 | for (j=1; j<=size(G0); j++) |
---|
| 3124 | { |
---|
| 3125 | for (i=1; i<=size(G0[j]); i++) |
---|
| 3126 | { |
---|
| 3127 | if (size(G0[j][i])!=0) |
---|
| 3128 | { |
---|
| 3129 | allmin[size(allmin)+1]=G0[j][i][1]; |
---|
| 3130 | allmax[size(allmax)+1]=G0[j][i][2]; |
---|
| 3131 | } |
---|
| 3132 | } |
---|
| 3133 | } |
---|
| 3134 | list minmax=list(Min(allmin),Max(allmax)); |
---|
| 3135 | return(minmax); |
---|
| 3136 | } |
---|
| 3137 | |
---|
| 3138 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3139 | |
---|
| 3140 | static proc bFctIdealModified (ideal I) |
---|
| 3141 | {/*modified version of the procedure bfunIdeal from bfun.lib*/ |
---|
| 3142 | def B= basering; |
---|
| 3143 | int n = nvars(B) div 2; |
---|
| 3144 | intvec w=(1:n); |
---|
| 3145 | I= initialIdealW(I,-w,w); |
---|
| 3146 | poly s; int i; |
---|
| 3147 | for (i=1; i<=n; i++) |
---|
| 3148 | { |
---|
| 3149 | s=s+x(i)*D(i); |
---|
| 3150 | } |
---|
| 3151 | /*pIntersect computes the intersection on s and I*/ |
---|
| 3152 | vector b = pIntersect(s,I); |
---|
| 3153 | list RL = ringlist(B); RL = RL[1..4]; |
---|
| 3154 | RL[2] = list(safeVarName("s")); |
---|
| 3155 | RL[3] = list(list("dp",intvec(1)),list("C",intvec(0))); |
---|
| 3156 | def @S = ring(RL); setring @S; |
---|
| 3157 | vector b = imap(B,b); |
---|
| 3158 | poly bs = vec2poly(b); |
---|
| 3159 | ring r=0,s,dp; |
---|
| 3160 | poly bs=imap(@S,bs); |
---|
| 3161 | /*find minimal and maximal integer root*/ |
---|
| 3162 | ideal allfac=factorize(bs,1); |
---|
| 3163 | list allfacs; |
---|
| 3164 | for (i=1; i<=ncols(allfac); i++) |
---|
| 3165 | { |
---|
| 3166 | allfacs[i]=allfac[i]; |
---|
| 3167 | } |
---|
| 3168 | number testzero; |
---|
| 3169 | list zeros; |
---|
| 3170 | for (i=1; i<=size(allfacs); i++) |
---|
| 3171 | { |
---|
| 3172 | if (deg(allfacs[i])==1) |
---|
| 3173 | { |
---|
| 3174 | testzero=number(subst(allfacs[i],s,0))/leadcoef(allfacs[i]); |
---|
| 3175 | if (testzero-int(testzero)==0) |
---|
| 3176 | { |
---|
| 3177 | zeros[size(zeros)+1]=int(-1)*int(testzero); |
---|
| 3178 | } |
---|
| 3179 | } |
---|
| 3180 | } |
---|
| 3181 | if (size(zeros)!=0) |
---|
| 3182 | { |
---|
| 3183 | list minmax=(Min(zeros),Max(zeros)); |
---|
| 3184 | } |
---|
| 3185 | else |
---|
| 3186 | { |
---|
| 3187 | list minmax=list(); |
---|
| 3188 | } |
---|
| 3189 | setring B; |
---|
| 3190 | return(minmax); |
---|
| 3191 | } |
---|
| 3192 | |
---|
| 3193 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3194 | |
---|
| 3195 | static proc safeVarName (string s) |
---|
| 3196 | {/* from the library "bfun.lib"*/ |
---|
| 3197 | string S = "," + charstr(basering) + "," + varstr(basering) + ","; |
---|
| 3198 | s = "," + s + ","; |
---|
| 3199 | while (find(S,s) <> 0) |
---|
| 3200 | { |
---|
| 3201 | s[1] = "@"; |
---|
| 3202 | s = "," + s; |
---|
| 3203 | } |
---|
| 3204 | s = s[2..size(s)-1]; |
---|
| 3205 | return(s) |
---|
| 3206 | } |
---|
| 3207 | |
---|
| 3208 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3209 | |
---|
| 3210 | static proc globalBFunOT(list L,list #) |
---|
| 3211 | { |
---|
| 3212 | /*this proc is currently not used since globalBFun computes the same output and is |
---|
| 3213 | faster, however globalBFun works only for non-zero b-functions!*/ |
---|
| 3214 | /*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]), |
---|
| 3215 | where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an |
---|
| 3216 | intvec of size n_i. |
---|
| 3217 | We compute bounds for the minimal and maximal integer roots of the b-functions |
---|
| 3218 | of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def. |
---|
| 3219 | 6.1.1 in [R]) using Algorithm 6.1.6 in [R].*/ |
---|
| 3220 | if (size(#)==0) |
---|
| 3221 | { |
---|
| 3222 | string Syzstring="Sres"; |
---|
| 3223 | } |
---|
| 3224 | else |
---|
| 3225 | { |
---|
| 3226 | string Syzstring=#[1]; |
---|
| 3227 | } |
---|
| 3228 | int i; int j; |
---|
| 3229 | def W=basering; |
---|
| 3230 | int n=nvars(W) div 2; |
---|
| 3231 | list G0; |
---|
| 3232 | ideal I; |
---|
| 3233 | intvec i1; |
---|
| 3234 | for (j=1; j<=size(L); j++) |
---|
| 3235 | { |
---|
| 3236 | G0[j]=list(); |
---|
| 3237 | for (i=1; i<=ncols(L[j][1]); i++) |
---|
| 3238 | { |
---|
| 3239 | G0[j][i]=I; |
---|
| 3240 | } |
---|
| 3241 | } |
---|
| 3242 | list out; |
---|
| 3243 | for (j=1; j<=size(L); j++) |
---|
| 3244 | { |
---|
| 3245 | if (L[j][2]!=intvec(0:size(L[j][2]))) |
---|
| 3246 | { |
---|
| 3247 | if (Syzstring=="Vdres") |
---|
| 3248 | { |
---|
| 3249 | L[j][1]=VdStrictGB(L[j][1],n); |
---|
| 3250 | } |
---|
| 3251 | else |
---|
| 3252 | { |
---|
| 3253 | def HomWeyl=makeHomogenizedWeyl(n); |
---|
| 3254 | setring HomWeyl; |
---|
| 3255 | list L=fetch(W,L); |
---|
| 3256 | L[j][1]=nHomogenize(L[j][1]); |
---|
| 3257 | L[j][1]=transpose(matrix(slimgb(transpose(L[j][1])))); |
---|
| 3258 | L[j][1]=subst(L[j][1],h,1); |
---|
| 3259 | setring W; |
---|
| 3260 | L=fetch(HomWeyl,L); |
---|
| 3261 | kill HomWeyl; |
---|
| 3262 | } |
---|
| 3263 | } |
---|
| 3264 | for (i=1; i<=nrows(L[j][1]); i++) |
---|
| 3265 | { |
---|
| 3266 | i1=(1..ncols(L[j][1])); |
---|
| 3267 | out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1); |
---|
| 3268 | G0[j][out[2]][size(G0[j][out[2]])+1]=(out[1]); |
---|
| 3269 | } |
---|
| 3270 | } |
---|
| 3271 | list Data=ringlist(W); |
---|
| 3272 | for (i=1; i<=n; i++) |
---|
| 3273 | { |
---|
| 3274 | Data[2][2*n+i]=Data[2][i]; |
---|
| 3275 | Data[2][3*n+i]=Data[2][n+i]; |
---|
| 3276 | Data[2][i]="v("+string(i)+")"; |
---|
| 3277 | Data[2][n+i]="w("+string(i)+")"; |
---|
| 3278 | } |
---|
| 3279 | Data[3][1][1]="M"; |
---|
| 3280 | intvec mord=(0:16*n^2); |
---|
| 3281 | mord[1..2*n]=(1:2*n); |
---|
| 3282 | mord[6*n+1..8*n]=(1:2*n); |
---|
| 3283 | for (i=0; i<=2*n-2; i++) |
---|
| 3284 | { |
---|
| 3285 | mord[(3+i)*4*n-i]=-1; |
---|
| 3286 | mord[(2*n+2+i)*4*n-2*n-i]=-1; |
---|
| 3287 | } |
---|
| 3288 | Data[3][1][2]=mord; |
---|
| 3289 | matrix Ones=UpOneMatrix(4*n); |
---|
| 3290 | Data[5]=Ones; |
---|
| 3291 | matrix con[2*n][2*n]; |
---|
| 3292 | Data[6]=transpose(concat(con,transpose(concat(con,Data[6])))); |
---|
| 3293 | def Wuv=ring(Data); |
---|
| 3294 | setring Wuv; |
---|
| 3295 | list G0=imap(W,G0); list G3; poly lterm;intvec lexp; |
---|
| 3296 | list G1,G2,LL; |
---|
| 3297 | intvec e,f; |
---|
| 3298 | int kapp,k,l; |
---|
| 3299 | poly h; |
---|
| 3300 | ideal I; |
---|
| 3301 | for (l=1; l<=size(G0); l++) |
---|
| 3302 | { |
---|
| 3303 | G1[l]=list(); G2[l]=list(); G3[l]=list(); |
---|
| 3304 | for (i=1; i<=size(G0[l]); i++) |
---|
| 3305 | { |
---|
| 3306 | for (j=1; j<=ncols(G0[l][i]);j++) |
---|
| 3307 | { |
---|
| 3308 | G0[l][i][j]=mHom(G0[l][i][j]); |
---|
| 3309 | } |
---|
| 3310 | for (j=1; j<=nvars(Wuv) div 4; j++) |
---|
| 3311 | { |
---|
| 3312 | G0[l][i][size(G0[l][i])+1]=1-v(j)*w(j); |
---|
| 3313 | } |
---|
| 3314 | G1[l][i]=slimgb(G0[l][i]); |
---|
| 3315 | G2[l][i]=I; |
---|
| 3316 | G3[l][i]=list(); |
---|
| 3317 | for (j=1; j<=ncols(G1[l][i]); j++) |
---|
| 3318 | { |
---|
| 3319 | e=leadexp(G1[l][i][j]); |
---|
| 3320 | f=e[1..2*n]; |
---|
| 3321 | if (f==intvec(0:(2*n))) |
---|
| 3322 | { |
---|
| 3323 | for (k=1; k<=n; k++) |
---|
| 3324 | { |
---|
| 3325 | kapp=-e[2*n+k]+e[3*n+k]; |
---|
| 3326 | if (kapp>0) |
---|
| 3327 | { |
---|
| 3328 | G1[l][i][j]=(x(k)^kapp)*G1[l][i][j]; |
---|
| 3329 | } |
---|
| 3330 | if (kapp<0) |
---|
| 3331 | { |
---|
| 3332 | G1[l][i][j]=(D(k)^(-kapp))*G1[l][i][j]; |
---|
| 3333 | } |
---|
| 3334 | } |
---|
| 3335 | G2[l][i][size(G2[l][i])+1]=G1[l][i][j]; |
---|
| 3336 | G3[l][i][size(G3[l][i])+1]=list(); |
---|
| 3337 | while (G1[l][i][j]!=0) |
---|
| 3338 | { |
---|
| 3339 | lterm=lead(G1[l][i][j]); |
---|
| 3340 | G1[l][i][j]=G1[l][i][j]-lterm; |
---|
| 3341 | lexp=leadexp(lterm); |
---|
| 3342 | lexp=lexp[2*n+1..3*n]; |
---|
| 3343 | LL=list(lexp,leadcoef(lterm)); |
---|
| 3344 | G3[l][i][size(G3[l][i])][size(G3[l][i][size(G3[l][i])])+1]=LL; |
---|
| 3345 | } |
---|
| 3346 | } |
---|
| 3347 | } |
---|
| 3348 | } |
---|
| 3349 | } |
---|
| 3350 | ring r=0,(s(1..n)),dp; |
---|
| 3351 | ideal I; |
---|
| 3352 | map G3forr=Wuv,I; |
---|
| 3353 | list G3=G3forr(G3); |
---|
| 3354 | poly fs,gs; |
---|
| 3355 | int a; |
---|
| 3356 | list G4; |
---|
| 3357 | for (l=1; l<=size(G3); l++) |
---|
| 3358 | { |
---|
| 3359 | G4[l]=list(); |
---|
| 3360 | for (i=1; i<=size(G3[l]);i++) |
---|
| 3361 | { |
---|
| 3362 | G4[l][i]=I; |
---|
| 3363 | |
---|
| 3364 | for (j=1; j<=size(G3[l][i]); j++) |
---|
| 3365 | { |
---|
| 3366 | fs=0; |
---|
| 3367 | for (k=1; k<=size(G3[l][i][j]); k++) |
---|
| 3368 | { |
---|
| 3369 | gs=1; |
---|
| 3370 | for (a=1; a<=n; a++) |
---|
| 3371 | { |
---|
| 3372 | if (G3[l][i][j][k][1][a]!=0) |
---|
| 3373 | { |
---|
| 3374 | gs=gs*permuteVar(list(G3[l][i][j][k][1][a]),a); |
---|
| 3375 | } |
---|
| 3376 | } |
---|
| 3377 | gs=gs*G3[l][i][j][k][2]; |
---|
| 3378 | fs=fs+gs; |
---|
| 3379 | } |
---|
| 3380 | G4[l][i]=G4[l][i],fs; |
---|
| 3381 | } |
---|
| 3382 | } |
---|
| 3383 | } |
---|
| 3384 | if (n==1) |
---|
| 3385 | { |
---|
| 3386 | ring rnew=0,t,dp; |
---|
| 3387 | } |
---|
| 3388 | else |
---|
| 3389 | { |
---|
| 3390 | ring rnew=0,(t,s(2..n)),dp; |
---|
| 3391 | } |
---|
| 3392 | ideal Iformap; |
---|
| 3393 | Iformap[1]=t; |
---|
| 3394 | poly forel=1; |
---|
| 3395 | for (i=2; i<=n; i++) |
---|
| 3396 | { |
---|
| 3397 | Iformap[1]=Iformap[1]-s(i); |
---|
| 3398 | Iformap[i]=s(i); |
---|
| 3399 | forel=forel*s(i); |
---|
| 3400 | } |
---|
| 3401 | map rtornew=r,Iformap; |
---|
| 3402 | list G4=rtornew(G4); |
---|
| 3403 | list getintvecs=fetch(W,L); |
---|
| 3404 | ideal J; |
---|
| 3405 | option(redSB); |
---|
| 3406 | for (l=1; l<=size(G4); l++) |
---|
| 3407 | { |
---|
| 3408 | J=1; |
---|
| 3409 | for (i=1; i<=size(G4[l]); i++) |
---|
| 3410 | { |
---|
| 3411 | G4[l][i]=eliminate(G4[l][i],forel); |
---|
| 3412 | J=intersect(J,G4[l][i]); |
---|
| 3413 | } |
---|
| 3414 | G4[l]=poly(std(J)[1]); |
---|
| 3415 | } |
---|
| 3416 | list minmax; |
---|
| 3417 | list mini=list(); |
---|
| 3418 | list maxi=list(); |
---|
| 3419 | list L=fetch(W,L); |
---|
| 3420 | for (i=1; i<=size(G4); i++) |
---|
| 3421 | { |
---|
| 3422 | minmax[i]=minIntRoot(G4[i],1); |
---|
| 3423 | if (size(minmax[i])!=0) |
---|
| 3424 | { |
---|
| 3425 | mini=insert(mini,minmax[i][1]+Min(L[i][2])); |
---|
| 3426 | maxi=insert(maxi,minmax[i][2]+Max(L[i][2])); |
---|
| 3427 | } |
---|
| 3428 | } |
---|
| 3429 | mini=Min(mini); |
---|
| 3430 | maxi=Max(maxi); |
---|
| 3431 | minmax=list(mini[1],maxi[1]); |
---|
| 3432 | option(none); |
---|
| 3433 | return(minmax); |
---|
| 3434 | } |
---|
| 3435 | |
---|
| 3436 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3437 | //COMPUTATION OF THE COHOMOLOGY |
---|
| 3438 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3439 | |
---|
| 3440 | static proc findCohomology(list L,int le) |
---|
| 3441 | { |
---|
| 3442 | /*computes the cohomology of the complex (D^i,d^i) given by D^i=C^L[2*i-1] and |
---|
| 3443 | d^i=L[2*i]*/ |
---|
| 3444 | def R=basering; |
---|
| 3445 | ring r=0,(x),dp; |
---|
| 3446 | list L=imap(R,L); |
---|
| 3447 | list out; |
---|
| 3448 | int i, ker, im; |
---|
| 3449 | matrix S; |
---|
| 3450 | option(returnSB); |
---|
| 3451 | option(redSB); |
---|
| 3452 | for (i=2; i<=size(L); i=i+2) |
---|
| 3453 | { |
---|
| 3454 | if (L[i-1]==0) |
---|
| 3455 | { |
---|
| 3456 | out[i div 2]=0; |
---|
| 3457 | im=0; |
---|
| 3458 | } |
---|
| 3459 | else |
---|
| 3460 | { |
---|
| 3461 | S=matrix(syz(transpose(L[i]))); |
---|
| 3462 | if (S!=matrix(0,nrows(S),ncols(S))) |
---|
| 3463 | { |
---|
| 3464 | ker=ncols(S); |
---|
| 3465 | out[i div 2]=ker-im; |
---|
| 3466 | im=L[i-1]-ker; |
---|
| 3467 | } |
---|
| 3468 | else |
---|
| 3469 | { |
---|
| 3470 | out[i-1]=0; |
---|
| 3471 | im=L[i-1]; |
---|
| 3472 | } |
---|
| 3473 | } |
---|
| 3474 | } |
---|
| 3475 | option(none); |
---|
| 3476 | while (size(out)>le) |
---|
| 3477 | { |
---|
| 3478 | out=delete(out,1); |
---|
| 3479 | } |
---|
| 3480 | setring R; |
---|
| 3481 | return(out); |
---|
| 3482 | } |
---|
| 3483 | |
---|
| 3484 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3485 | //AUXILIARY PROCEDURES |
---|
| 3486 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3487 | |
---|
| 3488 | static proc divdr(matrix m,matrix n) |
---|
| 3489 | { |
---|
| 3490 | m=transpose(m); |
---|
| 3491 | n=transpose(n); |
---|
| 3492 | matrix con=concat(m,n); |
---|
| 3493 | matrix s=syz(con); |
---|
| 3494 | s=submat(s,1..ncols(m),1..ncols(s)); |
---|
| 3495 | s=transpose(compress(s)); |
---|
| 3496 | return(s); |
---|
| 3497 | } |
---|
| 3498 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3499 | |
---|
| 3500 | static proc matrixLift(matrix M,matrix N) |
---|
| 3501 | { |
---|
| 3502 | intvec v=option(get); |
---|
| 3503 | option(none); |
---|
| 3504 | matrix l=transpose(lift(transpose(M),transpose(N))); |
---|
| 3505 | option(set,v); |
---|
| 3506 | return(l); |
---|
| 3507 | } |
---|
| 3508 | |
---|
| 3509 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3510 | |
---|
| 3511 | static proc VdStrictGB (matrix M,int d,list #) |
---|
| 3512 | "USAGE:VdStrictGB(M,d[,v]); M a matrix, d an integer, v an optional intvec |
---|
| 3513 | ASSUME:-basering is the nth Weyl algebra D_n @* |
---|
| 3514 | -1<=d<=n @* |
---|
| 3515 | -v (if given) is the shift vector on the range of M (in particular, |
---|
| 3516 | size(v)=ncols(M)); otherwise v is assumed to be the zero shift vector |
---|
| 3517 | RETURN:matrix N; the rows of N form a V_d-strict Groebner basis with respect to v |
---|
| 3518 | for the module generated by the rows of M |
---|
| 3519 | " |
---|
| 3520 | { |
---|
| 3521 | if (M==matrix(0,nrows(M),ncols(M))) |
---|
| 3522 | { |
---|
| 3523 | return (matrix(0,1,ncols(M))); |
---|
| 3524 | } |
---|
| 3525 | intvec op=option(get); |
---|
| 3526 | def W =basering; |
---|
| 3527 | int ncM=ncols(M); |
---|
| 3528 | list Data=ringlist(W); |
---|
| 3529 | Data[2]=list("nhv")+Data[2]; |
---|
| 3530 | Data[3][3]=Data[3][1]; |
---|
| 3531 | Data[3][1]=list("dp",intvec(1)); |
---|
| 3532 | matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2])); |
---|
| 3533 | Data[5]=re; |
---|
| 3534 | int k,l; |
---|
| 3535 | Data[6]=transpose(concat(matrix(0,1,1),transpose(concat(matrix(0,1,1),Data[6])))); |
---|
| 3536 | def Whom=ring(Data);// D_n[nhv] with the new commuative variable nhv |
---|
| 3537 | setring Whom; |
---|
| 3538 | matrix Mnew=imap(W,M); |
---|
| 3539 | intvec v; |
---|
| 3540 | if (size(#)!=0) |
---|
| 3541 | { |
---|
| 3542 | v=#[1]; |
---|
| 3543 | } |
---|
| 3544 | if (size(v) < ncM) |
---|
| 3545 | { |
---|
| 3546 | v=v,0:(ncM-size(v)); |
---|
| 3547 | } |
---|
| 3548 | Mnew=homogenize(Mnew, d, v);//homogenization of M with respect to the new variable |
---|
| 3549 | Mnew=transpose(Mnew); |
---|
| 3550 | Mnew=slimgb(Mnew);// computes a Groebner basis of the homogenzition of M |
---|
| 3551 | Mnew=subst(Mnew,nhv,1);// substitution of 1 gives V_d-strict Groebner basis of M |
---|
| 3552 | Mnew=compress(Mnew); |
---|
| 3553 | Mnew=transpose(Mnew); |
---|
| 3554 | setring W; |
---|
| 3555 | M=imap(Whom,Mnew); |
---|
| 3556 | option(set,op); |
---|
| 3557 | return(M); |
---|
| 3558 | } |
---|
| 3559 | |
---|
| 3560 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3561 | |
---|
| 3562 | static proc VdNormalForm(matrix F,matrix M,int d,intvec v,list #) |
---|
| 3563 | "USAGE:VdNormalForm(F,M,d,v[,w]); F and M matrices, d int, v intvec, w an optional |
---|
| 3564 | intvec |
---|
| 3565 | ASSUME:-basering is the nth Weyl algebra D_n @* |
---|
| 3566 | -F a n_1 x n_2-matrix and M a m_1 x m_2-matrix with m_2<=n_2 @* |
---|
| 3567 | -d is an integer between 1 and n @* |
---|
| 3568 | -v is a shift vector for D_n^(m_2) and hence size(v)=m_2 @* |
---|
| 3569 | -w is a shift vector for D_n^(m_1-m_2) and hence size(v)=m_1-m_2 @* |
---|
| 3570 | RETURN:a n_1 x n_2-matrix N such that:@* |
---|
| 3571 | -If no optional intvec w is given:(N[i,1],..,N[i,m_2]) is a V_d-strict normal |
---|
| 3572 | form of (F[i,1],...,F[i,m_2]) with respect to a V_d-strict Groebner basis of |
---|
| 3573 | the rows of M and the shift vector v |
---|
| 3574 | -If w is given:(N[i,1],..,N[i,m_2]) is chosen such that |
---|
| 3575 | Vddeg((N[i,1],...,N[i,m_2])[v])<=Vddeg((F[i,m_2+1],...,F[i,m_1])[v]); |
---|
| 3576 | -N[i,j]=F[i,j] for j>m_2 |
---|
| 3577 | " |
---|
| 3578 | { |
---|
| 3579 | int SBcom; |
---|
| 3580 | def W =basering; |
---|
| 3581 | int c=ncols(M); |
---|
| 3582 | matrix keepF=F; |
---|
| 3583 | if (size(#)!=0) |
---|
| 3584 | { |
---|
| 3585 | intvec w=#[1]; |
---|
| 3586 | } |
---|
| 3587 | F=submat(F,intvec(1..nrows(F)),intvec(1..c)); |
---|
| 3588 | list Data=ringlist(W); |
---|
| 3589 | Data[2]=list("nhv")+Data[2]; |
---|
| 3590 | Data[3][3]=Data[3][1]; |
---|
| 3591 | Data[3][1]=list("dp",intvec(1)); |
---|
| 3592 | matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2])); |
---|
| 3593 | Data[5]=re; |
---|
| 3594 | int k,l,nr,nc; |
---|
| 3595 | matrix rep[size(Data[2])][size(Data[2])]; |
---|
| 3596 | for (l=size(Data[2])-1;l>=1; l--) |
---|
| 3597 | { |
---|
| 3598 | for (k=l-1; k>=1;k--) |
---|
| 3599 | { |
---|
| 3600 | rep[k+1,l+1]=Data[6][k,l]; |
---|
| 3601 | } |
---|
| 3602 | } |
---|
| 3603 | Data[6]=rep; |
---|
| 3604 | def Whom=ring(Data);//new ring D_n[nvh] this new commuative variable nhv |
---|
| 3605 | setring Whom; |
---|
| 3606 | matrix Mnew=imap(W,M); |
---|
| 3607 | list forMnew=homogenize(Mnew,d,v,1);//commputes homogenization of M; |
---|
| 3608 | Mnew=forMnew[1]; |
---|
| 3609 | int rightexp=forMnew[2]; |
---|
| 3610 | matrix Fnew=imap(W,F); |
---|
| 3611 | matrix keepF=imap(W,keepF); |
---|
| 3612 | matrix Fb; |
---|
| 3613 | int cc; |
---|
| 3614 | intvec i1,i2; |
---|
| 3615 | matrix zeromat,subm1,subm2,zeromat2; |
---|
| 3616 | for (l=1; l<=nrows(Fnew); l++) |
---|
| 3617 | { |
---|
| 3618 | if (size(#)!=0) |
---|
| 3619 | { |
---|
| 3620 | subm2=submat(keepF,l,((ncols(Fnew)+1)..ncols(keepF))); |
---|
| 3621 | zeromat2=matrix(0,1,ncols(subm2)); |
---|
| 3622 | if (submat(keepF,l,((ncols(Fnew)+1)..ncols(keepF)))==zeromat2) |
---|
| 3623 | { |
---|
| 3624 | for (cc=1; cc<=ncols(Fnew); c++) |
---|
| 3625 | { |
---|
| 3626 | Fnew[l,cc]=0; |
---|
| 3627 | } |
---|
| 3628 | } |
---|
| 3629 | i1=intvec(1..ncols(Fnew)); |
---|
| 3630 | subm1=submat(Fnew,l,i1); |
---|
| 3631 | subm2=submat(keepF,l,(ncols(Fnew)+1)..ncols(keepF)); |
---|
| 3632 | zeromat=matrix(0,1,ncols(Fnew)); |
---|
| 3633 | if (VdDegnhv(subm1,d,v)>VdDegnhv(subm2,d,w) |
---|
| 3634 | and submat(Fnew,l,intvec(1..ncols(Fnew)))!=zeromat) |
---|
| 3635 | { |
---|
| 3636 | //print("Reduzierung des V_d-Grades nötig"); |
---|
| 3637 | /*We need to reduce the V_d-degree. First we homogenize the |
---|
| 3638 | lth row of Fnew*/ |
---|
| 3639 | Fb=homogenize(subm1,d,v)*(nhv^rightexp); |
---|
| 3640 | if (SBcom==0) |
---|
| 3641 | { |
---|
| 3642 | /*computes a V_d-strict standard basis*/ |
---|
| 3643 | Mnew=slimgb(transpose(Mnew));// |
---|
| 3644 | SBcom=1; |
---|
| 3645 | } |
---|
| 3646 | /*computes a V_d-strict normal form for FB*/ |
---|
| 3647 | Fb=transpose(reduce(transpose(Fb),Mnew)); |
---|
| 3648 | if (VdDegnhv(Fb,d,v)> VdDegnhv(subm2,d,w) |
---|
| 3649 | and Fb!=matrix(0,nrows(Fb),ncols(Fb)))//should not happen |
---|
| 3650 | { |
---|
| 3651 | //print("Reduzierung fehlgeschlagen!!!!!!!!!!!!!!!!"); |
---|
| 3652 | } |
---|
| 3653 | } |
---|
| 3654 | else |
---|
| 3655 | { |
---|
| 3656 | /*condition on V_ddeg already satisfied -> no normal form |
---|
| 3657 | computation is needed*/ |
---|
| 3658 | Fb=submat(Fnew,l,intvec(1..ncols(Fnew))); |
---|
| 3659 | } |
---|
| 3660 | } |
---|
| 3661 | else |
---|
| 3662 | { |
---|
| 3663 | Fb=homogenize(submat(Fnew,l,intvec(1..ncols(Fnew))),d,v); |
---|
| 3664 | if (SBcom==0) |
---|
| 3665 | { |
---|
| 3666 | Mnew=slimgb(transpose(Mnew));// computes a V_d-strict Groebner basis |
---|
| 3667 | SBcom=1; |
---|
| 3668 | } |
---|
| 3669 | Fb=transpose(reduce(transpose(Fb),Mnew));//normal form |
---|
| 3670 | } |
---|
| 3671 | for (k=1; k<=ncols(Fnew);k++) |
---|
| 3672 | { |
---|
| 3673 | Fnew[l,k]=Fb[1,k]; |
---|
| 3674 | } |
---|
| 3675 | } |
---|
| 3676 | Fnew=subst(Fnew,nhv,1);//obtain normal form in D_n |
---|
| 3677 | setring W; |
---|
| 3678 | F=imap(Whom,Fnew); |
---|
| 3679 | return(F); |
---|
| 3680 | } |
---|
| 3681 | |
---|
| 3682 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3683 | |
---|
| 3684 | static proc homogenize (matrix M,int d,intvec v,list #) |
---|
| 3685 | { |
---|
| 3686 | /* we compute the F[v]-homogenization of each row of M (cf. Def. 3.4 in [OT])*/ |
---|
| 3687 | if (M==matrix(0,nrows(M),ncols(M))) |
---|
| 3688 | { |
---|
| 3689 | return(M); |
---|
| 3690 | } |
---|
| 3691 | int i,l,s, kmin, nhvexp; |
---|
| 3692 | poly f; |
---|
| 3693 | intvec vnm; |
---|
| 3694 | list findmin,maxnhv,rempoly,remk,rem1,rem2; |
---|
| 3695 | int n=(nvars(basering)-1) div 2; |
---|
| 3696 | for (int k=1; k<=nrows(M); k++) |
---|
| 3697 | { |
---|
| 3698 | for (l=1; l<=ncols (M); l++) |
---|
| 3699 | { |
---|
| 3700 | f=M[k,l]; |
---|
| 3701 | s=size(f); |
---|
| 3702 | for (i=1; i<=s; i++) |
---|
| 3703 | { |
---|
| 3704 | vnm=leadexp(f); |
---|
| 3705 | vnm=vnm[n+2..n+d+1]-vnm[2..d+1]; |
---|
| 3706 | kmin=sum(vnm)+v[l]; |
---|
| 3707 | rem1[size(rem1)+1]=lead(f); |
---|
| 3708 | rem2[size(rem2)+1]=kmin; |
---|
| 3709 | findmin=insert(findmin,kmin); |
---|
| 3710 | f=f-lead(f); |
---|
| 3711 | } |
---|
| 3712 | rempoly[l]=rem1; |
---|
| 3713 | remk[l]=rem2; |
---|
| 3714 | rem1=list(); |
---|
| 3715 | rem2=list(); |
---|
| 3716 | } |
---|
| 3717 | if (size(findmin)!=0) |
---|
| 3718 | { |
---|
| 3719 | kmin=Min(findmin); |
---|
| 3720 | } |
---|
| 3721 | for (l=1; l<=ncols(M); l++) |
---|
| 3722 | { |
---|
| 3723 | if (M[k,l]!=0) |
---|
| 3724 | { |
---|
| 3725 | M[k,l]=0; |
---|
| 3726 | for (i=1; i<=size(rempoly[l]);i++) |
---|
| 3727 | { |
---|
| 3728 | nhvexp=remk[l][i]-kmin; |
---|
| 3729 | M[k,l]=M[k,l]+nhv^(nhvexp)*rempoly[l][i]; |
---|
| 3730 | maxnhv[size(maxnhv)+1]=nhvexp; |
---|
| 3731 | } |
---|
| 3732 | } |
---|
| 3733 | } |
---|
| 3734 | rempoly=list(); |
---|
| 3735 | remk=list(); |
---|
| 3736 | findmin=list(); |
---|
| 3737 | } |
---|
| 3738 | maxnhv=Max(maxnhv); |
---|
| 3739 | nhvexp=maxnhv[1]; |
---|
| 3740 | if (size(#)!=0) |
---|
| 3741 | { |
---|
| 3742 | return(list(M,nhvexp));//only needed for normal form computations |
---|
| 3743 | } |
---|
| 3744 | return(M); |
---|
| 3745 | } |
---|
| 3746 | |
---|
| 3747 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3748 | |
---|
| 3749 | static proc soldr (matrix M,matrix N) |
---|
| 3750 | { |
---|
| 3751 | /* We compute a ncols(M) x nrows(M)-matrix C such that |
---|
| 3752 | C[i,1]M_1+...+C[i,nrows(M)]M_(nrows(M))= e_i mod im(N), |
---|
| 3753 | where e_i is the ith basis element on the range of M, M_j denotes the jth row |
---|
| 3754 | of M and im(N) is generated by the rows of N */ |
---|
| 3755 | int n=nrows(M); |
---|
| 3756 | int q=ncols(M); |
---|
| 3757 | matrix S=concat(transpose(M),transpose(N)); |
---|
| 3758 | def W=basering; |
---|
| 3759 | list Data=ringlist(W); |
---|
| 3760 | list Save=Data[3]; |
---|
| 3761 | Data[3]=list(list("c",0),list("dp",intvec(1..nvars(W)))); |
---|
| 3762 | def Wmod=ring(Data); |
---|
| 3763 | setring Wmod; |
---|
| 3764 | matrix Smod=imap(W,S); |
---|
| 3765 | matrix E[q][1]; |
---|
| 3766 | matrix Smod2,Smodnew; |
---|
| 3767 | option(returnSB); |
---|
| 3768 | int i,j; |
---|
| 3769 | for (i=1;i<=q;i++) |
---|
| 3770 | { |
---|
| 3771 | E[i,1]=1; |
---|
| 3772 | Smod2=concat(E,Smod); |
---|
| 3773 | Smod2=syz(Smod2); |
---|
| 3774 | E[i,1]=0; |
---|
| 3775 | for (j=1;j<=ncols(Smod2);j++) |
---|
| 3776 | { |
---|
| 3777 | if (Smod2[1,j]==1) |
---|
| 3778 | { |
---|
| 3779 | Smodnew=concat(Smodnew,(-1)*(submat(Smod2,intvec(2..n+1),j))); |
---|
| 3780 | break; |
---|
| 3781 | } |
---|
| 3782 | } |
---|
| 3783 | } |
---|
| 3784 | Smodnew=transpose(submat(Smodnew,intvec(1..n),intvec(2..q+1))); |
---|
| 3785 | setring W; |
---|
| 3786 | matrix Snew=imap(Wmod,Smodnew); |
---|
| 3787 | option(none); |
---|
| 3788 | return (Snew); |
---|
| 3789 | } |
---|
| 3790 | |
---|
| 3791 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3792 | |
---|
| 3793 | static proc prodr (int k,int l) |
---|
| 3794 | { |
---|
| 3795 | if (k==0) |
---|
| 3796 | { |
---|
| 3797 | matrix P=unitmat(l); |
---|
| 3798 | return (P); |
---|
| 3799 | } |
---|
| 3800 | matrix O[l][k]; |
---|
| 3801 | matrix P=transpose(concat(O,unitmat(l))); |
---|
| 3802 | return (P); |
---|
| 3803 | } |
---|
| 3804 | |
---|
| 3805 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3806 | |
---|
| 3807 | static proc VdDeg(matrix M,int d,intvec v,list #) |
---|
| 3808 | { |
---|
| 3809 | /* We assume that the basering it the nth Weyl algebra and that M is a 1 x r- |
---|
| 3810 | matrix. |
---|
| 3811 | We compute the V_d-deg of M with respect to the shift vector v, |
---|
| 3812 | i.e V_ddeg(M)=max (V_ddeg(M_i)+v[i]), where k=V_ddeg(M_i) if k is the minimal |
---|
| 3813 | integer, such that M_i can be expressed as a sum of operators |
---|
| 3814 | x(1)^(a_1)*...*x(n)^(a_n)*D(1)^(b_1)*...*D(n)^(b_n) with |
---|
| 3815 | a_1+..+a_d+k>=b_1+..+b_d*/ |
---|
| 3816 | int i, j, etoint; |
---|
| 3817 | int n=nvars(basering) div 2; |
---|
| 3818 | intvec e; |
---|
| 3819 | list findmax; |
---|
| 3820 | int c=ncols(M); |
---|
| 3821 | poly l; |
---|
| 3822 | list positionpoly,positionVd; |
---|
| 3823 | for (i=1; i<=c; i++) |
---|
| 3824 | { |
---|
| 3825 | positionpoly[i]=list(); |
---|
| 3826 | positionVd[i]=list(); |
---|
| 3827 | while (M[1,i]!=0) |
---|
| 3828 | { |
---|
| 3829 | l=lead(M[1,i]); |
---|
| 3830 | positionpoly[i][size(positionpoly[i])+1]=l; |
---|
| 3831 | e=leadexp(l); |
---|
| 3832 | e=-e[1..d]+e[n+1..n+d]; |
---|
| 3833 | e=sum(e)+v[i]; |
---|
| 3834 | etoint=e[1]; |
---|
| 3835 | positionVd[i][size(positionVd[i])+1]=etoint; |
---|
| 3836 | findmax[size(findmax)+1]=etoint; |
---|
| 3837 | M[1,i]=M[1,i]-l; |
---|
| 3838 | } |
---|
| 3839 | } |
---|
| 3840 | if (size(findmax)!=0) |
---|
| 3841 | { |
---|
| 3842 | int maxVd=Max(findmax); |
---|
| 3843 | if (size(#)==0) |
---|
| 3844 | { |
---|
| 3845 | return (maxVd); |
---|
| 3846 | } |
---|
| 3847 | } |
---|
| 3848 | else // M is 0-modul |
---|
| 3849 | { |
---|
| 3850 | return(int(0)); |
---|
| 3851 | } |
---|
| 3852 | l=0; |
---|
| 3853 | for (i=c; i>=1; i--) |
---|
| 3854 | { |
---|
| 3855 | for (j=1; j<=size(positionVd[i]); j++) |
---|
| 3856 | { |
---|
| 3857 | if (positionVd[i][j]==maxVd) |
---|
| 3858 | { |
---|
| 3859 | l=l+positionpoly[i][j]; |
---|
| 3860 | } |
---|
| 3861 | } |
---|
| 3862 | if (l!=0) |
---|
| 3863 | { |
---|
| 3864 | /*returns the largest component that has maximal V_d-degree |
---|
| 3865 | and its terms of maximal V_d-deg (needed for globalBFun)*/ |
---|
| 3866 | return (list(l,i)); |
---|
| 3867 | } |
---|
| 3868 | } |
---|
| 3869 | } |
---|
| 3870 | |
---|
| 3871 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3872 | |
---|
| 3873 | static proc VdDegnhv(matrix M,int d,intvec v,list #) |
---|
| 3874 | { |
---|
| 3875 | /* As the procedure VdDeg, but the basering is the nth Weyl algebra |
---|
| 3876 | with a commutative variable nhv*/ |
---|
| 3877 | int i,j,etoint; |
---|
| 3878 | int n=nvars(basering) div 2; |
---|
| 3879 | intvec e; |
---|
| 3880 | int etoint; |
---|
| 3881 | list findmax; |
---|
| 3882 | int c=ncols(M); |
---|
| 3883 | poly l; |
---|
| 3884 | list positionpoly; |
---|
| 3885 | list positionVd; |
---|
| 3886 | for (i=1; i<=c; i++) |
---|
| 3887 | { |
---|
| 3888 | positionpoly[i]=list(); |
---|
| 3889 | positionVd[i]=list(); |
---|
| 3890 | while (M[1,i]!=0) |
---|
| 3891 | { |
---|
| 3892 | l=lead(M[1,i]); |
---|
| 3893 | positionpoly[i][size(positionpoly[i])+1]=l; |
---|
| 3894 | e=leadexp(l); |
---|
| 3895 | e=-e[2..d+1]+e[n+2..n+d+1]; |
---|
| 3896 | e=sum(e)+v[i]; |
---|
| 3897 | etoint=e[1]; |
---|
| 3898 | positionVd[i][size(positionVd[i])+1]=etoint; |
---|
| 3899 | findmax[size(findmax)+1]=etoint; |
---|
| 3900 | M[1,i]=M[1,i]-l; |
---|
| 3901 | } |
---|
| 3902 | } |
---|
| 3903 | if (size(findmax)!=0) |
---|
| 3904 | { |
---|
| 3905 | int maxVd=Max(findmax); |
---|
| 3906 | if (size(#)==0) |
---|
| 3907 | { |
---|
| 3908 | return (maxVd); |
---|
| 3909 | } |
---|
| 3910 | } |
---|
| 3911 | else // M is 0-modul |
---|
| 3912 | { |
---|
| 3913 | return(int(0)); |
---|
| 3914 | } |
---|
| 3915 | } |
---|
| 3916 | |
---|
| 3917 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3918 | |
---|
| 3919 | static proc deletecol(matrix M,int l) |
---|
| 3920 | { |
---|
| 3921 | int s=ncols(M); |
---|
| 3922 | if (l==1) |
---|
| 3923 | { |
---|
| 3924 | M=submat(M,(1..nrows(M)),(2..ncols(M))); |
---|
| 3925 | return(M); |
---|
| 3926 | } |
---|
| 3927 | if (l==s) |
---|
| 3928 | { |
---|
| 3929 | M=submat(M,(1..nrows(M)),(1..(ncols(M)-1))); |
---|
| 3930 | return(M); |
---|
| 3931 | } |
---|
| 3932 | intvec v=(1..(l-1)),((l+1)..s); |
---|
| 3933 | M=submat(M,(1..nrows(M)),v); |
---|
| 3934 | return(M); |
---|
| 3935 | } |
---|
| 3936 | |
---|
| 3937 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3938 | |
---|
| 3939 | static proc mHom(poly f) |
---|
| 3940 | {/*for globalBFunOT*/ |
---|
| 3941 | poly g; |
---|
| 3942 | poly l; |
---|
| 3943 | poly add; |
---|
| 3944 | intvec e; |
---|
| 3945 | list minint; |
---|
| 3946 | list remf; |
---|
| 3947 | int i; |
---|
| 3948 | int j; |
---|
| 3949 | int n=nvars(basering) div 4; |
---|
| 3950 | if (f==0) |
---|
| 3951 | { |
---|
| 3952 | return(f); |
---|
| 3953 | } |
---|
| 3954 | while (f!=0) |
---|
| 3955 | { |
---|
| 3956 | l=lead(f); |
---|
| 3957 | e=leadexp(l); |
---|
| 3958 | remf[size(remf)+1]=list(); |
---|
| 3959 | remf[size(remf)][1]=l; |
---|
| 3960 | for (i=1; i<=n; i++) |
---|
| 3961 | { |
---|
| 3962 | remf[size(remf)][i+1]=-e[2*n+i]+e[3*n+i]; |
---|
| 3963 | if (size(minint)<i) |
---|
| 3964 | { |
---|
| 3965 | minint[i]=list(); |
---|
| 3966 | } |
---|
| 3967 | minint[i][size(minint[i])+1]=-e[2*n+i]+e[3*n+i]; |
---|
| 3968 | } |
---|
| 3969 | f=f-l; |
---|
| 3970 | } |
---|
| 3971 | for (i=1; i<=n; i++) |
---|
| 3972 | { |
---|
| 3973 | minint[i]=Min(minint[i]); |
---|
| 3974 | } |
---|
| 3975 | for (i=1; i<=size(remf); i++) |
---|
| 3976 | { |
---|
| 3977 | add=remf[i][1]; |
---|
| 3978 | for (j=1; j<=n; j++) |
---|
| 3979 | { |
---|
| 3980 | add=v(j)^(remf[i][j+1]-minint[j])*add; |
---|
| 3981 | } |
---|
| 3982 | g=g+add; |
---|
| 3983 | } |
---|
| 3984 | return (g); |
---|
| 3985 | } |
---|
| 3986 | |
---|
| 3987 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 3988 | |
---|
| 3989 | static proc permuteVar(list L,int n) |
---|
| 3990 | {/*for globalBFunOT*/ |
---|
| 3991 | if (typeof(L[1])=="intvec") |
---|
| 3992 | { |
---|
| 3993 | intvec v=L[1]; |
---|
| 3994 | } |
---|
| 3995 | else |
---|
| 3996 | { |
---|
| 3997 | intvec v=(1:L[1]),(0:L[1]); |
---|
| 3998 | } |
---|
| 3999 | int i;int k; int indi=0; |
---|
| 4000 | int j; |
---|
| 4001 | int s=size(v); |
---|
| 4002 | poly e; |
---|
| 4003 | intvec fore; |
---|
| 4004 | for (i=2; i<=size(v); i=i+2) |
---|
| 4005 | { |
---|
| 4006 | |
---|
| 4007 | if (v[i]!=0) |
---|
| 4008 | { |
---|
| 4009 | j=i+1; |
---|
| 4010 | while (v[j]!=0) |
---|
| 4011 | { |
---|
| 4012 | j=j+1; |
---|
| 4013 | } |
---|
| 4014 | v[i]=0; |
---|
| 4015 | v[j]=1; |
---|
| 4016 | fore=0; |
---|
| 4017 | indi=0; |
---|
| 4018 | for (k=1; k<=size(v); k++) |
---|
| 4019 | { |
---|
| 4020 | if (k!=i and k!=j) |
---|
| 4021 | { |
---|
| 4022 | if (indi==0) |
---|
| 4023 | { |
---|
| 4024 | indi=1; |
---|
| 4025 | fore[1]=v[k]; |
---|
| 4026 | } |
---|
| 4027 | else |
---|
| 4028 | { |
---|
| 4029 | fore[size(fore)+1]=v[k]; |
---|
| 4030 | } |
---|
| 4031 | } |
---|
| 4032 | } |
---|
| 4033 | e=e-(j-i)*permutevar(list(fore),n); |
---|
| 4034 | } |
---|
| 4035 | } |
---|
| 4036 | e=e+s(n)^(size(v) div 2); |
---|
| 4037 | return (e); |
---|
| 4038 | } |
---|
| 4039 | |
---|
| 4040 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4041 | |
---|
| 4042 | static proc makeHomogenizedWeyl(int n,list #) |
---|
| 4043 | { |
---|
| 4044 | /*modified version of the procedure makeWeyl() from the library nctools.lib*/ |
---|
| 4045 | /*Creates the nth homogenized Weyl algebra with variables x(1),..,x(n),D(1),.., |
---|
| 4046 | D(n) and homogenization variable h, i.e. it holds x(i)*D(i)=D(i)*x(1)+h^2. |
---|
| 4047 | If # contains on intvec v, we assign weight v[i] to the ith module component.*/ |
---|
| 4048 | if (n<1) |
---|
| 4049 | { |
---|
| 4050 | print("Incorrect input"); |
---|
| 4051 | return(); |
---|
| 4052 | } |
---|
| 4053 | if (n ==1) |
---|
| 4054 | { |
---|
| 4055 | ring @rr = 0,(x(1),D(1),h),dp; |
---|
| 4056 | } |
---|
| 4057 | else |
---|
| 4058 | { |
---|
| 4059 | ring @rr = 0,(x(1..n),D(1..n),h),dp; |
---|
| 4060 | } |
---|
| 4061 | setring @rr; |
---|
| 4062 | if (size(#)==0) |
---|
| 4063 | { |
---|
| 4064 | def @rrr = homogenizedWeyl(); |
---|
| 4065 | } |
---|
| 4066 | else |
---|
| 4067 | { |
---|
| 4068 | def @rrr=homogenizedWeyl(#); |
---|
| 4069 | } |
---|
| 4070 | return(@rrr); |
---|
| 4071 | } |
---|
| 4072 | |
---|
| 4073 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4074 | |
---|
| 4075 | static proc homogenizedWeyl (list #) |
---|
| 4076 | { |
---|
| 4077 | /*modified version of the procedure Weyl() from the library nctools.lib*/ |
---|
| 4078 | /*Creates a homogenized Weyl algebra structure on the basering. We assume |
---|
| 4079 | n=nvars(basering) is odd. The first (n-1)/2 variables will be treated as the |
---|
| 4080 | x(i), the next (n-1)/2 as the corresponding differentials D(i) and the last as |
---|
| 4081 | the homogenization variable h, i.e. it holds x(i)*D(i)=D(i)*x(1)+h^2. |
---|
| 4082 | If # contains on intvec v, we assign weight v[i] to the ith module component.*/ |
---|
| 4083 | string rname=nameof(basering); |
---|
| 4084 | if ( rname == "basering") // i.e. no ring has been set yet |
---|
| 4085 | { |
---|
| 4086 | "You have to call the procedure from the ring"; |
---|
| 4087 | return(); |
---|
| 4088 | } |
---|
| 4089 | int nv = nvars(basering); |
---|
| 4090 | int N = (nv-1) div 2; |
---|
| 4091 | if (((nv-1) % 2) != 0) |
---|
| 4092 | { |
---|
| 4093 | "Cannot create homogenized Weyl structure for an even number of generators"; |
---|
| 4094 | return(); |
---|
| 4095 | } |
---|
| 4096 | matrix @D[nv][nv]; |
---|
| 4097 | int i; |
---|
| 4098 | for ( i=1; i<=N; i++ ) |
---|
| 4099 | { |
---|
| 4100 | @D[i,N+i]=h^2; |
---|
| 4101 | } |
---|
| 4102 | def @R = nc_algebra(1,@D); |
---|
| 4103 | setring @R; |
---|
| 4104 | list RL=ringlist(@R); |
---|
| 4105 | intvec v; |
---|
| 4106 | /*we need this ordering for Groebner basis computations*/ |
---|
| 4107 | for (i=1; i<=N; i++) |
---|
| 4108 | { |
---|
| 4109 | v[i]=-1; |
---|
| 4110 | v[N+i]=1; |
---|
| 4111 | } |
---|
| 4112 | v[nv]=0; |
---|
| 4113 | /* we assign weights to module components*/ |
---|
| 4114 | if (size(#)!=0) |
---|
| 4115 | { |
---|
| 4116 | if (typeof(#[1])=="intvec") |
---|
| 4117 | { |
---|
| 4118 | intvec m=#[1]; |
---|
| 4119 | for (i=1; i<=size(m); i++) |
---|
| 4120 | { |
---|
| 4121 | v[size(v)+1]=m[i];//assigns weight m[i] to the ith module component |
---|
| 4122 | } |
---|
| 4123 | RL[3]=insert(RL[3],list("am",v)); |
---|
| 4124 | } |
---|
| 4125 | else |
---|
| 4126 | { |
---|
| 4127 | RL[3]=insert(RL[3],list("a",v)); |
---|
| 4128 | } |
---|
| 4129 | } |
---|
| 4130 | else |
---|
| 4131 | { |
---|
| 4132 | RL[3]=insert(RL[3],list("a",v)); |
---|
| 4133 | } |
---|
| 4134 | intvec w=(1:nv); |
---|
| 4135 | if (size(#)>=2) |
---|
| 4136 | { |
---|
| 4137 | if (typeof(#[2])=="intvec") |
---|
| 4138 | { |
---|
| 4139 | intvec n=#[2]; |
---|
| 4140 | for (i=1; i<=size(n); i++) |
---|
| 4141 | { |
---|
| 4142 | w[size(w)+1]=n[i]; |
---|
| 4143 | } |
---|
| 4144 | RL[3]=insert(RL[3],list("am",w)); |
---|
| 4145 | } |
---|
| 4146 | else |
---|
| 4147 | { |
---|
| 4148 | RL[3]=insert(RL[3],list("a",w)); |
---|
| 4149 | } |
---|
| 4150 | } |
---|
| 4151 | else |
---|
| 4152 | { |
---|
| 4153 | RL[3]=insert(RL[3],list("a",w)); |
---|
| 4154 | } |
---|
| 4155 | /*this ordering is needed for globalBFun and globalBFunOT*/ |
---|
| 4156 | list saveord=RL[3][3]; |
---|
| 4157 | RL[3][3]=RL[3][4]; |
---|
| 4158 | RL[3][4]=saveord; |
---|
| 4159 | intvec notforh=(1:(size(RL[3][4][2])-1)); |
---|
| 4160 | RL[3][4][2]=notforh; |
---|
| 4161 | RL[3][5]=list("dp",1); |
---|
| 4162 | def @@R=ring(RL); |
---|
| 4163 | return(@@R); |
---|
| 4164 | } |
---|
| 4165 | |
---|
| 4166 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4167 | |
---|
| 4168 | static proc nHomogenize (matrix M,list #) |
---|
| 4169 | { |
---|
| 4170 | /* # may contain an intvec v, if no intvec is given, we assume that v=(0:ncols(M)) |
---|
| 4171 | We compute the h[v]-homogenization of the rows of M as in Definition 9.2 [OT]*/ |
---|
| 4172 | int l; poly f; int s; int i; intvec vnm;int kmin; list findmax; |
---|
| 4173 | int n=(nvars(basering)-1) div 2; |
---|
| 4174 | list rempoly; |
---|
| 4175 | list remk; |
---|
| 4176 | list rem1; |
---|
| 4177 | list rem2; |
---|
| 4178 | list maxhexp; |
---|
| 4179 | int hexp; |
---|
| 4180 | intvec v=(0:ncols(M)); |
---|
| 4181 | if (size(#)!=0) |
---|
| 4182 | { |
---|
| 4183 | if (typeof(#[1])=="intvec") |
---|
| 4184 | { |
---|
| 4185 | v=#[1]; |
---|
| 4186 | } |
---|
| 4187 | } |
---|
| 4188 | if (size(v)<ncols(M)) |
---|
| 4189 | { |
---|
| 4190 | for (i=size(v)+1; i<=ncols(M); i++) |
---|
| 4191 | { |
---|
| 4192 | v[i]=0; |
---|
| 4193 | } |
---|
| 4194 | } |
---|
| 4195 | for (int k=1; k<=nrows(M); k++) |
---|
| 4196 | { |
---|
| 4197 | for (l=1; l<=ncols (M); l++) |
---|
| 4198 | { |
---|
| 4199 | f=M[k,l]; |
---|
| 4200 | s=size(f); |
---|
| 4201 | for (i=1; i<=s; i++) |
---|
| 4202 | { |
---|
| 4203 | vnm=leadexp(f); |
---|
| 4204 | kmin=sum(vnm)+v[l]; |
---|
| 4205 | rem1[size(rem1)+1]=lead(f); |
---|
| 4206 | rem2[size(rem2)+1]=kmin; |
---|
| 4207 | findmax=insert(findmax,kmin); |
---|
| 4208 | f=f-lead(f); |
---|
| 4209 | } |
---|
| 4210 | rempoly[l]=rem1; |
---|
| 4211 | remk[l]=rem2; |
---|
| 4212 | rem1=list(); |
---|
| 4213 | rem2=list(); |
---|
| 4214 | } |
---|
| 4215 | if (size(findmax)!=0) |
---|
| 4216 | { |
---|
| 4217 | kmin=Max(findmax); |
---|
| 4218 | } |
---|
| 4219 | else |
---|
| 4220 | { |
---|
| 4221 | kmin=0; |
---|
| 4222 | } |
---|
| 4223 | for (l=1; l<=ncols(M); l++) |
---|
| 4224 | { |
---|
| 4225 | if (M[k,l]!=0) |
---|
| 4226 | { |
---|
| 4227 | M[k,l]=0; |
---|
| 4228 | for (i=1; i<=size(rempoly[l]);i++) |
---|
| 4229 | { |
---|
| 4230 | hexp=kmin-remk[l][i]; |
---|
| 4231 | maxhexp[size(maxhexp)+1]=hexp; |
---|
| 4232 | M[k,l]=M[k,l]+h^hexp*rempoly[l][i]; |
---|
| 4233 | } |
---|
| 4234 | } |
---|
| 4235 | } |
---|
| 4236 | rempoly=list(); |
---|
| 4237 | remk=list(); |
---|
| 4238 | findmax=list(); |
---|
| 4239 | } |
---|
| 4240 | if (size(maxhexp)!=0) |
---|
| 4241 | { |
---|
| 4242 | maxhexp=Max(maxhexp); |
---|
| 4243 | hexp=maxhexp[1]; |
---|
| 4244 | } |
---|
| 4245 | else |
---|
| 4246 | { |
---|
| 4247 | hexp=0; |
---|
| 4248 | } |
---|
| 4249 | if (size(#)>1) |
---|
| 4250 | { |
---|
| 4251 | list forreturn=M,hexp; |
---|
| 4252 | |
---|
| 4253 | return(forreturn); |
---|
| 4254 | } |
---|
| 4255 | return(M); |
---|
| 4256 | } |
---|
| 4257 | |
---|
| 4258 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4259 | |
---|
| 4260 | static proc max(int i,int j) |
---|
| 4261 | { |
---|
| 4262 | if(i>j){return(i);} |
---|
| 4263 | return(j); |
---|
| 4264 | } |
---|
| 4265 | |
---|
| 4266 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4267 | |
---|
| 4268 | static proc nDeg (matrix M,intvec m) |
---|
| 4269 | {/*we compute an intvec n such that n[i]=max(deg(M[i,j])+m[j]|M[i,j]!=0) (where deg |
---|
| 4270 | stands for the total degree) if (M[i,j]!=0 for some j) and n[i]=0 else*/ |
---|
| 4271 | int i; int j; |
---|
| 4272 | intvec n; |
---|
| 4273 | list L; |
---|
| 4274 | for (i=1; i<=nrows(M); i++) |
---|
| 4275 | { |
---|
| 4276 | L=list(); |
---|
| 4277 | for (j=1; j<=ncols(M); j++) |
---|
| 4278 | { |
---|
| 4279 | if (M[i,j]!=0) |
---|
| 4280 | { |
---|
| 4281 | L=insert(L,deg(M[i,j])+m[j]); |
---|
| 4282 | } |
---|
| 4283 | } |
---|
| 4284 | if (size(L)==0) |
---|
| 4285 | { |
---|
| 4286 | n[i]=0; |
---|
| 4287 | } |
---|
| 4288 | else |
---|
| 4289 | { |
---|
| 4290 | n[i]=Max(L); |
---|
| 4291 | } |
---|
| 4292 | } |
---|
| 4293 | return(n); |
---|
| 4294 | } |
---|
| 4295 | |
---|
| 4296 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4297 | |
---|
| 4298 | static proc minIntRoot(list L,list #) |
---|
| 4299 | "USAGE:minIntRoot(L [,M]); L list, M optinonal list |
---|
| 4300 | ASSUME:L a list of univariate polynomials with rational coefficients @* |
---|
| 4301 | the variable of the polynomial is s if size(#)==0 (needed for proc |
---|
| 4302 | MVComplex) and t else (needed for globalBFun) |
---|
| 4303 | RETURN:-if size(#)==0: int i, where i is an integer root of one of the polynomials |
---|
| 4304 | and it is minimal with respect to that property@* |
---|
| 4305 | -if size(#)!=0: list L=(i,j), where i is as above and j is an integer root |
---|
| 4306 | of one of the polynomials and is maximal with respect to that property (if |
---|
| 4307 | an integer root exists) or L=list() else |
---|
| 4308 | " |
---|
| 4309 | { |
---|
| 4310 | def B=basering; |
---|
| 4311 | if (size(#)==0) |
---|
| 4312 | { |
---|
| 4313 | ring rnew=0,s,dp; |
---|
| 4314 | } |
---|
| 4315 | else |
---|
| 4316 | { |
---|
| 4317 | ring rnew=0,t,dp; |
---|
| 4318 | } |
---|
| 4319 | list L=imap(B,L); |
---|
| 4320 | |
---|
| 4321 | int i; |
---|
| 4322 | int j; |
---|
| 4323 | number isint; |
---|
| 4324 | list possmin; |
---|
| 4325 | ideal allfac; |
---|
| 4326 | list allfacs; |
---|
| 4327 | for (i=1; i<=size(L); i++) |
---|
| 4328 | { |
---|
| 4329 | allfac=factorize(L[i],1); |
---|
| 4330 | for (j=1; j<=ncols(allfac); j++) |
---|
| 4331 | { |
---|
| 4332 | allfacs[j]=allfac[j]; |
---|
| 4333 | } |
---|
| 4334 | for (j=1; j<=size(allfacs); j++) |
---|
| 4335 | { |
---|
| 4336 | if (deg(allfacs[j])==1) |
---|
| 4337 | { |
---|
| 4338 | isint=number(subst(allfacs[j],var(1),0)/leadcoef(allfacs[j])); |
---|
| 4339 | if (isint-int(isint)==0) |
---|
| 4340 | { |
---|
| 4341 | possmin[size(possmin)+1]=int(isint); |
---|
| 4342 | } |
---|
| 4343 | } |
---|
| 4344 | } |
---|
| 4345 | allfacs=list(); |
---|
| 4346 | } |
---|
| 4347 | int zerolist; |
---|
| 4348 | if (size(possmin)!=0) |
---|
| 4349 | { |
---|
| 4350 | int miniroot=(-1)*Max(possmin); |
---|
| 4351 | int maxiroot=(-1)*Min(possmin); |
---|
| 4352 | } |
---|
| 4353 | else |
---|
| 4354 | { |
---|
| 4355 | zerolist=1; |
---|
| 4356 | } |
---|
| 4357 | setring B; |
---|
| 4358 | if (size(#)==0) |
---|
| 4359 | { |
---|
| 4360 | return(miniroot); |
---|
| 4361 | } |
---|
| 4362 | else |
---|
| 4363 | { |
---|
| 4364 | if (zerolist==0) |
---|
| 4365 | { |
---|
| 4366 | return(list(miniroot,maxiroot)); |
---|
| 4367 | } |
---|
| 4368 | else |
---|
| 4369 | { |
---|
| 4370 | return(list()); |
---|
| 4371 | } |
---|
| 4372 | } |
---|
| 4373 | } |
---|
| 4374 | |
---|
| 4375 | |
---|
| 4376 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4377 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4378 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4379 | /* |
---|
| 4380 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4381 | FURTHER EXAMPLES FOR TESTING THE PROCEDURES |
---|
| 4382 | //////////////////////////////////////////////////////////////////////////////////// |
---|
| 4383 | LIB "derham.lib"; |
---|
| 4384 | |
---|
| 4385 | //---------------------------------------- |
---|
| 4386 | //EXAMPLE 1 |
---|
| 4387 | //---------------------------------------- |
---|
| 4388 | ring r=0,(x,y),dp; |
---|
| 4389 | poly f=y2-x3-2x+3; |
---|
| 4390 | list L=deRhamCohomology(f); |
---|
| 4391 | L; |
---|
| 4392 | kill r; |
---|
| 4393 | |
---|
| 4394 | //---------------------------------------- |
---|
| 4395 | //EXAMPLE 2 |
---|
| 4396 | //---------------------------------------- |
---|
| 4397 | ring r=0,(x,y),dp; |
---|
| 4398 | poly f=y2-x3-x; |
---|
| 4399 | list L=deRhamCohomology(f); |
---|
| 4400 | L; |
---|
| 4401 | kill r; |
---|
| 4402 | |
---|
| 4403 | //---------------------------------------- |
---|
| 4404 | //EXAMPLE 3 |
---|
| 4405 | //---------------------------------------- |
---|
| 4406 | ring r=0,(x,y),dp; |
---|
| 4407 | list C=list(x2-1,(x+1)*y,y*(y2+2x+1)); |
---|
| 4408 | list L=deRhamCohomology(C); |
---|
| 4409 | L; |
---|
| 4410 | kill r; |
---|
| 4411 | |
---|
| 4412 | //---------------------------------------- |
---|
| 4413 | //EXAMPLE 4 |
---|
| 4414 | //---------------------------------------- |
---|
| 4415 | ring r=0,(x,y,z),dp; |
---|
| 4416 | list C=list(x*(x-1),y,z*(z-1),z*(x-1)); |
---|
| 4417 | list L=deRhamCohomology(C); |
---|
| 4418 | L; |
---|
| 4419 | kill r; |
---|
| 4420 | |
---|
| 4421 | //---------------------------------------- |
---|
| 4422 | //EXAMPLE 5 |
---|
| 4423 | //---------------------------------------- |
---|
| 4424 | ring r=0,(x,y,z),dp; |
---|
| 4425 | list C=list(x*y,y*z); |
---|
| 4426 | list L=deRhamCohomology(C,"Vdres"); |
---|
| 4427 | L; |
---|
| 4428 | kill r; |
---|
| 4429 | |
---|
| 4430 | //---------------------------------------- |
---|
| 4431 | //EXAMPLE 6 |
---|
| 4432 | //---------------------------------------- |
---|
| 4433 | ring r=0,(x,y,z,u),dp; |
---|
| 4434 | list C=list(x,y,z,u); |
---|
| 4435 | list L=deRhamCohomology(C); |
---|
| 4436 | L; |
---|
| 4437 | kill r; |
---|
| 4438 | |
---|
| 4439 | //---------------------------------------- |
---|
| 4440 | //EXAMPLE 7 |
---|
| 4441 | //---------------------------------------- |
---|
| 4442 | ring r=0,(x,y,z),dp; |
---|
| 4443 | poly f=x3+y3+z3; |
---|
| 4444 | list L=deRhamCohomology(f); |
---|
| 4445 | L; |
---|
| 4446 | kill r; |
---|
| 4447 | |
---|
| 4448 | //---------------------------------------- |
---|
| 4449 | //EXAMPLE 8 |
---|
| 4450 | //---------------------------------------- |
---|
| 4451 | ring r=0,(x,y,z),dp; |
---|
| 4452 | poly f=x2+y2+z2; |
---|
| 4453 | list L=deRhamCohomology(f,"Vdres"); |
---|
| 4454 | L; |
---|
| 4455 | kill r; |
---|
| 4456 | |
---|
| 4457 | //---------------------------------------- |
---|
| 4458 | //EXAMPLE 9 |
---|
| 4459 | //---------------------------------------- |
---|
| 4460 | ring r=0,(x,y,z,u),dp; |
---|
| 4461 | list C=list(x2+y2+z2,u); |
---|
| 4462 | list L=deRhamCohomology(C); |
---|
| 4463 | L; |
---|
| 4464 | kill r; |
---|
| 4465 | |
---|
| 4466 | |
---|
| 4467 | //---------------------------------------- |
---|
| 4468 | //EXAMPLE 10 |
---|
| 4469 | //---------------------------------------- |
---|
| 4470 | ring r=0,(x,y,z),dp; |
---|
| 4471 | list C=list((x*(y-1),y2-1)); |
---|
| 4472 | list L=deRhamCohomology(C); |
---|
| 4473 | L; |
---|
| 4474 | kill r; |
---|
| 4475 | |
---|
| 4476 | |
---|
| 4477 | */ |
---|