Annotation of OpenXM/src/math-misc/twistedLogCohomology.m2, Revision 1.1
1.1 ! nisitani 1: -- source code from Dbasic.m2
! 2: -- puts a module or matrix purely in shift degree 0.
! 3:
! 4: zeroize2 = method()
! 5: zeroize2 Module := M -> (
! 6: W := ring M;
! 7: P := presentation M;
! 8: coker map(W^(numgens target P), W^(numgens source P), P)
! 9: )
! 10:
! 11: zeroize2 Matrix := m -> (
! 12: W := ring m;
! 13: map(W^(numgens target m), W^(numgens source m), m)
! 14: )
! 15:
! 16: -- source code from Dresolution.m2
! 17: -- modified Dresolution
! 18: -- return hashTable{Dresolution,Transfermap}
! 19:
! 20: shifts := (m, w, oldshifts) -> (
! 21: tempmat := compress leadTerm m;
! 22: if numgens source tempmat == 0 then newshifts := {}
! 23: else (
! 24: expmat := matrix(apply(toList(0..numgens source tempmat - 1),
! 25: i -> (k := leadComponent tempmat_i;
! 26: append((exponents tempmat_(k,i))#0, oldshifts#k))));
! 27: newshifts = (entries transpose (
! 28: expmat*(transpose matrix{ append(w, 1) })) )#0;
! 29: );
! 30: newshifts)
! 31:
! 32: debug Core
! 33:
! 34: kerGB := m -> (
! 35: -- m should be a matrix which is a GB, and
! 36: -- whose source has the Schreyer order.
! 37: -- The resulting map will have the same form.
! 38: map(ring m, rawKernelOfGB raw m)
! 39: )
! 40:
! 41: Dresolution2 = method( Options => {Strategy => Schreyer, LengthLimit => infinity} )
! 42: Dresolution2 Ideal := options -> I -> (
! 43: Dresolution2((ring I)^1/I, options)
! 44: )
! 45:
! 46: Dresolution2 Module := options -> M -> (
! 47:
! 48: pInfo (1, "ENTERING Dresolution ... ");
! 49:
! 50: W := ring M;
! 51: N := presentation M;
! 52:
! 53: pInfo (1, "Computing usual resolution using Schreyer order ...");
! 54: pInfo (2, "含t Degree " | 0 | "...");
! 55: pInfo (2, "含t含t含t Rank = " | rank target N | "含t time = 0. seconds");
! 56: pInfo (2, "含t Degree " | 1 | "...");
! 57: tInfo := toString first timing (m := schreyerOrder gens gb N);
! 58: pInfo (2, "含t含t含t Rank = " | rank source N | "含t time = " |
! 59: tInfo | " seconds");
! 60:
! 61: M.cache.resolution = new ChainComplex;
! 62: M.cache.resolution.ring = W;
! 63: M.cache.resolution#0 = target m;
! 64: M.cache.resolution#1 = source m;
! 65: M.cache.resolution.dd#0 = map(W^0, target m, 0);
! 66: M.cache.resolution.dd#1 = m;
! 67:
! 68: i := 2;
! 69: while source m != 0 and i <= options.LengthLimit do (
! 70: pInfo (2, "含t Degree " | i | "...");
! 71: tInfo = toString first timing (m = kerGB m);
! 72: M.cache.resolution#i = source m;
! 73: M.cache.resolution.dd#i = m;
! 74: pInfo(2, "含t含t含t Rank = " | rank source m | "含t time = " |
! 75: tInfo | " seconds");
! 76: i = i+1;
! 77: );
! 78: M.cache.resolution.length = i-1;
! 79: M.cache.resolution
! 80: )
! 81:
! 82: Dresolution2 (Ideal, List) := options -> (I, w) -> (
! 83: Dresolution2 ((ring I)^1/I, w, options)
! 84: )
! 85:
! 86: Dresolution2 (Module, List) := options -> (M, w) -> (
! 87:
! 88: pInfo (1, "ENTERING Dresolution ... ");
! 89:
! 90: -- ERROR CHECKING:
! 91: W := ring M;
! 92: k := options.LengthLimit;
! 93:
! 94: -- check that W is a Weyl algebra
! 95: if W.monoid.Options.WeylAlgebra == {}
! 96: then error "expected a Weyl algebra";
! 97: if any(W.monoid.Options.WeylAlgebra, v -> class v =!= Option)
! 98: then error "expected non-homogenized Weyl algebra";
! 99: -- check that w is of the form (-u,u)
! 100: createDpairs W;
! 101: if #w =!= numgens W
! 102: then error ("expected weight vector of length " | numgens W);
! 103: if any( toList(0..#W.dpairInds#0 - 1),
! 104: i -> ( w#(W.dpairInds#0#i) + w#(W.dpairInds#1#i) != 0 ) )
! 105: then error "expected weight vector of the form (-u,u)";
! 106:
! 107: -- PREPROCESSING
! 108: if k == infinity then (
! 109: pInfo (1, "Computing adapted free resolution of length infinity using "
! 110: | toString options.Strategy | " method...");
! 111: if (options.Strategy == Vhomogenize) then
! 112: pInfo(2, "Warning: resolution via Vhomogenize might not terminate");
! 113: )
! 114: else pInfo (1, "Computing adapted free resolution of length " | k |
! 115: " using " | toString options.Strategy | " method...");
! 116:
! 117: homVar := symbol homVar;
! 118: hvw := symbol hvw;
! 119: if options.Strategy == Schreyer then (
! 120: -- Make the homogenizing weight vector in HW
! 121: Hwt := toList(numgens W + 1:1);
! 122: -- Make the V-filtration weight vector in HW
! 123: Vwt := append(w,0);
! 124: -- Make the homogeneous Weyl algebra
! 125: HW := (coefficientRing W)(monoid [(entries vars W)#0, homVar,
! 126: WeylAlgebra => append(W.monoid.Options.WeylAlgebra, homVar),
! 127: MonomialOrder => {Weights=>Hwt, Weights=>Vwt, GRevLex}]);
! 128: homVar = HW_homVar;
! 129: WtoHW := map(HW, W, (vars HW)_{0..numgens W - 1});
! 130: HWtoW := map(W, HW, (vars W)_{0..numgens W - 1} | matrix{{1_W}});
! 131: -- Also make the homogenizing Weyl algebra for shifts
! 132: VW := (coefficientRing W)(monoid [hvw, (entries vars W)#0,
! 133: WeylAlgebra => W.monoid.Options.WeylAlgebra,
! 134: MonomialOrder => Eliminate 1]);
! 135: HWtoVW := map(VW, HW, (vars VW)_{1..numgens W} | matrix{{VW_0}});
! 136: VWtoHW := map(HW, VW, matrix{{homVar}} | (vars HW)_{0..numgens HW - 2});
! 137: hvwVar := VW_0;
! 138: HVWwt := prepend(-1,w);
! 139: VWwt := prepend(0,w);
! 140: )
! 141: else if options.Strategy == Vhomogenize then (
! 142: Hwt = prepend(-1,w);
! 143: Vwt = prepend(0,w);
! 144: -- make the homogenizing Weyl algebra
! 145: HW = (coefficientRing W)(monoid [homVar, (entries vars W)#0,
! 146: WeylAlgebra => W.monoid.Options.WeylAlgebra,
! 147: MonomialOrder => Eliminate 1]);
! 148: homVar = HW_homVar;
! 149: WtoHW = map(HW, W, (vars HW)_{1..numgens W});
! 150: HWtoW = map(W, HW, matrix{{1_W}} | (vars W));
! 151: );
! 152:
! 153: -- CREATE AND INITIALIZE THE CHAIN COMPLEX
! 154: --else
! 155: N := presentation M;
! 156: --if (isSubmodule M) then N := presentation ((ambient M)/M);
! 157: -- get the degree shifts right (need to check this against OT paper)
! 158: if not M.cache.?resolution
! 159: then M.cache.resolution = new MutableHashTable;
! 160: M.cache.resolution#w = new ChainComplex;
! 161: M.cache.resolution#w.ring = W;
! 162: s := rank source N;
! 163: t := rank target N;
! 164: M.cache.resolution#w#0 = target N;
! 165: M.cache.resolution#w.dd#0 = map(W^0, M.cache.resolution#w#0, 0);
! 166:
! 167: -- MAKE THE FIRST STEP OF THE RESOLUTION
! 168: shiftvec := apply(degrees target N, i -> i#0);
! 169: tempMap := map(HW^(-shiftvec), HW^(rank source N), WtoHW N);
! 170: pInfo (2, "含t Degree 0...");
! 171: pInfo (2, "含t含t含t Rank = " | t | "含t time = 0 seconds");
! 172: pInfo (3, "含t Degree 1...");
! 173: tInfo := toString first timing (
! 174: Jgb := gb(homogenize(tempMap, homVar, Hwt), ChangeMatrix => true);
! 175: transfermap := getChangeMatrix Jgb;
! 176: Jgb = homogenize(tempMap, homVar, Hwt)*transfermap;
! 177: if options.Strategy == Schreyer then Jgb = schreyerOrder Jgb;
! 178: if options.Strategy == Schreyer then (
! 179: tempMat := map(VW^(-shiftvec), VW^(numgens source Jgb), HWtoVW(Jgb));
! 180: shiftvec = shifts(homogenize(HWtoVW Jgb, hvwVar, HVWwt),
! 181: VWwt, shiftvec);
! 182: )
! 183: else shiftvec = shifts(Jgb, Vwt, shiftvec);
! 184: transfermap = HWtoW transfermap;
! 185: M.cache.resolution#w#1 = W^(-shiftvec);
! 186: M.cache.resolution#w.dd#1 = map(M.cache.resolution#w#0,
! 187: M.cache.resolution#w#1, HWtoW Jgb);
! 188: );
! 189: pInfo(2, "含t含t含t Rank = " | #shiftvec | "含t time = " |
! 190: tInfo | " seconds");
! 191: startDeg := 2;
! 192:
! 193: -- COMPUTE REST OF THE RESOLUTION
! 194: i := startDeg;
! 195: while i < k+1 and numgens source Jgb != 0 do (
! 196: pInfo (2, "含t Degree " | i | "...");
! 197: tInfo = toString first timing (
! 198: if options.Strategy == Schreyer then Jgb = kerGB Jgb
! 199: else if options.Strategy == Vhomogenize then (
! 200: -- compute the kernel / syzygies
! 201: Jsyz := syz Jgb;
! 202: -- put syzygies in the free module with the correct degree shifts
! 203: Jsyzmap := map(HW^(-shiftvec), HW^(numgens source Jsyz), Jsyz);
! 204: -- compute an adapted (-w,w)-GB of the syzygies module
! 205: Jgb = gens gb homogenize(Jsyzmap, homVar, Hwt);
! 206: );
! 207: if options.Strategy == Schreyer then (
! 208: tempMat = map(VW^(-shiftvec), VW^(numgens source Jgb), HWtoVW(Jgb));
! 209: shiftvec = shifts(homogenize(tempMat, hvwVar, HVWwt),
! 210: VWwt, shiftvec);
! 211: )
! 212: else shiftvec = shifts(Jgb, Vwt, shiftvec);
! 213: M.cache.resolution#w#i = W^(-shiftvec);
! 214: M.cache.resolution#w.dd#i = map(M.cache.resolution#w#(i-1),
! 215: M.cache.resolution#w#i, HWtoW Jgb);
! 216: );
! 217: pInfo(2, "含t含t含t Rank = " | #shiftvec | "含t time = " |
! 218: tInfo | " seconds");
! 219: i = i+1;
! 220: );
! 221: hashTable {VResolution =>M.cache.resolution#w,
! 222: TransferMap => transfermap}
! 223: )
! 224:
! 225:
! 226: -- This routine computes a monomial from a list of variables
! 227: -- and an exponent vector
! 228: List ^ List := (Vars, Exps) -> (
! 229: product (Vars, Exps, (i,j) -> (i^j)) )
! 230:
! 231: -- This routine takes a weight vector w of strictly positive integers
! 232: -- and returns the exponents beta in N^n such that
! 233: -- ( k_0 < w_1 beta_1 + ... + w_n bet_n <= k_1)
! 234: findExps := (w, k0, k1) -> (
! 235: local minimum, local alpha, local tempExps;
! 236: -- base case: weight vector w has dim 1
! 237: if #w == 1 and k0 >= 0 then (
! 238: tempExps = (toList((k0//w#0+1)..k1//w#0) / (i -> {i}) ) )
! 239: else if #w == 1 and k0 < 0 then (
! 240: tempExps = ( toList(0..k1//w#0) / (i -> {i}) ) )
! 241: else ( -- general case
! 242: tempExps = {};
! 243: alpha = 0;
! 244: while alpha <= k1//w#0 do (
! 245: tempExps = join( tempExps,
! 246: apply ( findExps( drop(w,1), k0-alpha*(w#0),
! 247: k1-alpha*(w#0) ), j -> prepend(alpha,j) ) );
! 248: alpha = alpha+1;
! 249: );
! 250: );
! 251: tempExps)
! 252:
! 253:
! 254: -- modified computeRestriction
! 255:
! 256: computeRestriction2 = (M,wt,n0,n1,output,options) -> (
! 257:
! 258: -- ERROR CHECKING
! 259: W := ring M;
! 260: createDpairs W;
! 261: -- check weight vector
! 262: if #wt != #W.dpairInds#0
! 263: then error ("expected weight vector of length " | #W.dpairInds#0);
! 264: if any(wt, i -> (i<0))
! 265: then error "expected non-negative weight vector";
! 266:
! 267: -- PREPROCESSING
! 268: local tInfo;
! 269: local tempvar;
! 270: tempvar = symbol tempvar;
! 271: tInfo = "";
! 272: nW := numgens W;
! 273: -- make the (-w,w) weight vector
! 274: w := new MutableList from join(W.dpairInds#0,W.dpairInds#1);
! 275: i := 0;
! 276: while i < #W.dpairInds#0 do (
! 277: w#(W.dpairInds#0#i) = -wt#i;
! 278: w#(W.dpairInds#1#i) = wt#i;
! 279: i = i+1;
! 280: );
! 281: w = toList w;
! 282: d := #positions(w, i->(i>0));
! 283: -- the variables t_1, ..., t_d
! 284: negVars := (entries vars W)#0_(positions(w, i->(i<0)));
! 285: -- the substitution which sets t_1 = ... = t_d = 0
! 286: resSub := apply( negVars, i -> (i => 0) );
! 287: -- the variables Dt_1, ..., Dt_d, their indices, and their weights
! 288: posVars := (entries vars W)#0_(positions(w, i->(i>0)));
! 289: posInds := positions( w, i->(i>0) );
! 290: posWeights := select( w, i->(i>0) );
! 291: diffSub := apply( posVars, i -> (i => 0) );
! 292: -- the rest of the variables x_1, ..., x_n, D_1, ..., D_n
! 293: otherVars := (entries vars W)#0_(positions(w, i->(i==0)));
! 294:
! 295:
! 296: -- MAKE THE WEYL ALGEBRA "resW" OF THE RESTRICTED SUBSPACE
! 297: -- Case 1: if restriction to pt, then "resW" a field
! 298: if #otherVars == 0 then (
! 299: resW := coefficientRing W;
! 300: WtoresW := map(resW, W, matrix{toList(numgens W: 0_resW)});
! 301: resWtoW := map(W, resW)
! 302: )
! 303: -- Case 2: if restriction to coordinate subspace of dim d, then
! 304: -- resW a Weyl algebra D_d.
! 305: else (i = 0;
! 306: resPairsList := {};
! 307: while i < #otherVars do (
! 308: deriv := select(otherVars, j ->
! 309: (j*otherVars#i - otherVars#i*j == 1));
! 310: if (#deriv == 1) then
! 311: resPairsList = append(resPairsList, otherVars#i=>deriv#0);
! 312: i = i+1;
! 313: );
! 314: resW = (coefficientRing W)(monoid [otherVars, WeylAlgebra=>resPairsList]);
! 315: -- make the inclusion ring map "WtoresW" mapping W --> resW
! 316: counter := 0;
! 317: tempList := {};
! 318: WList := {};
! 319: i = 0;
! 320: while i < numgens W do (
! 321: if w#i == 0 then (
! 322: tempList = append(tempList, resW_counter);
! 323: WList = append(WList, W_i);
! 324: counter = counter+1;)
! 325: else (tempList = append(tempList, 0_resW));
! 326: i = i+1;
! 327: );
! 328: WtoresW = map(resW, W, matrix{tempList});
! 329: resWtoW = map(W, resW, matrix{WList});
! 330: );
! 331:
! 332: -- INITIALIZE THE OUTPUT FORMS
! 333:
! 334: if member(HomologyModules, output) then homologyList := {};
! 335: if member(Cycles, output) then kerList := {};
! 336: if member(Boundaries, output) then imList := {};
! 337: if member(GenCycles, output) then explicitList := {};
! 338: if member(RestrictComplex,output) then (
! 339: restrictComplex := new ChainComplex;
! 340: restrictComplex.ring = resW;
! 341: );
! 342: outputList := {};
! 343: if member(Cycles, output) or member(Boundaries, output) or
! 344: member(GenCycles, output) then explicitFlag := true
! 345: else explicitFlag = false;
! 346:
! 347: -- GET MIN AND MAX ROOTS OF THE B-FUNCTION
! 348: --<< "Computing b-function ";
! 349: if n0 >= d then b := 1_(QQ[tempvar])
! 350: else if rank ambient M == 1 then (
! 351: -- used to use "ideal presentation M" here
! 352: --<< "using ideal bFunction ";
! 353: pInfo(1, "Computing b-function using ideal bFunction... ");
! 354: tInfo = toString first timing(
! 355: b = bFunction(ideal relations M, wt);
! 356: );
! 357: )
! 358: else (
! 359: pInfo(1, "Computing b-function using module bFunction... ");
! 360: tInfo = toString first timing(
! 361: b = bFunctionM(M, wt, toList(rank ambient M: 0));
! 362: );
! 363: );
! 364:
! 365: if b == 0 then (
! 366: error "Module not specializable. Restriction cannot be computed.";
! 367: );
! 368:
! 369: intRoots := getIntRoots b;
! 370:
! 371: -- NO INTEGER ROOTS
! 372: if #intRoots == 0 then (
! 373: k0 := 0;
! 374: k1 := 0;
! 375: pInfo(4, "含t bFunction = " | toString b);
! 376: pInfo(2, "含t No integer roots");
! 377: pInfo(3, "含t time = " | tInfo);
! 378: if member(RestrictComplex, output) then (
! 379: restrictComplex#n0 = resW^0;
! 380: restrictComplex#n1 = resW^0;
! 381: restrictComplex.dd#n0 = map(resW^0,resW^0,0);
! 382: restrictComplex.dd#n1 = map(resW^0,resW^0,0);
! 383: i = n0+1;
! 384: while i < n1 do (
! 385: restrictComplex#i = resW^0;
! 386: restrictComplex.dd#i = map(resW^0,resW^0,0);
! 387: i = i+1;
! 388: );
! 389: );
! 390: if member(HomologyModules, output) then
! 391: homologyList = apply (toList(n0+1..n1-1), i -> (i => resW^0));
! 392: if member(Cycles, output) then
! 393: kerList = apply (toList(n0+1..n1-1), i -> (i => gens W^0));
! 394: if member(Boundaries, output) then
! 395: imList = apply (toList(n0+1..n1-1), i -> (i => gens W^0));
! 396: if member(GenCycles, output) then
! 397: explicitList = apply (toList(n0+1..n1-1), i -> (i => gens W^0));
! 398: )
! 399:
! 400: -- INTEGER ROOTS EXIST
! 401: else (
! 402: k0 = min intRoots - 1;
! 403: k1 = max intRoots;
! 404: pInfo(4, "含t bFunction = " | toString b);
! 405: pInfo(2, "含t min root = " | k0+1 | " , max root = " | k1);
! 406: pInfo(3, "含t time = " | tInfo | " seconds");
! 407: pInfo(2, " ");
! 408:
! 409: ----- SET k0 TO -infinity FOR EXPLICIT COHOMOLOGY CLASSES -----
! 410: if member(Explicit, output) then k0 = -infinity;
! 411:
! 412: -- COMPUTE FREE RESOLUTION ADAPTED TO THE WEIGHT VECTOR "w"
! 413: tInfo = toString first timing (
! 414: Resolution := Dresolution2 (M, w, LengthLimit => n1, Strategy => options.Strategy);
! 415: C := Resolution#VResolution;
! 416: );
! 417: pInfo(2, "含t Finished...");
! 418: pInfo(2, "含t含t含t Total time = " | tInfo | " seconds");
! 419: pInfo(2, " ");
! 420:
! 421:
! 422: -- COMPUTE THE RESTRICTED COMPLEX IN DEGREES "n0" TO "n1"
! 423: tInfo = toString first timing (
! 424: pInfo(1, "Computing induced restriction complex in degrees " |
! 425: n0 | " to " | n1 | "...");
! 426:
! 427: -- INITIALIZE THE ITERATION : making the first differential
! 428: pInfo(2, "含t Degree " | n0+1 | "...");
! 429:
! 430: -- MAKE THE TARGET MODULE AS DIRECT SUM OF D_m MODULES
! 431: -- "targetGens" is a list of s lists, where the i-th list contains
! 432: -- the monomial generators {dx^a} (as left D_m-module) of
! 433: -- F_k1[u_i]((D_n/tD_n) e_i) / F_k0[u_i]((D_n/tD_n) e_i)
! 434: tInfo = toString first timing (
! 435: s := numgens target C.dd#(n0+1);
! 436: targetDeg := degrees target C.dd#(n0+1);
! 437: targetGens := {};
! 438: if explicitFlag then targetMat := map(W^0,W^0,0);
! 439: i = 0;
! 440: while i < s do (
! 441: tempExps := findExps(posWeights, k0-targetDeg#i#0, k1-targetDeg#i#0);
! 442: tempGens := apply(tempExps, j -> posVars^j);
! 443: targetGens = append(targetGens, tempGens);
! 444: if explicitFlag then (
! 445: if tempGens == {} then (
! 446: targetMat = directSum(targetMat, compress matrix{{0_W}}); )
! 447: else (
! 448: targetMat = directSum(targetMat, matrix{tempGens}); );
! 449: );
! 450: i = i+1;
! 451: );
! 452: targetSize := sum(targetGens, i->#i);
! 453:
! 454: -- MAKE THE SOURCE MODULE AS DIRECT SUM OF D_m MODULES
! 455: -- "sourceGens" is a list of r lists, where the i-th list contains
! 456: -- the monomial generators {dx^a} (as left D_m-module) of
! 457: -- F_k1[u_i]((D_n/tD_n) e_i) / F_k0[u_i]((D_n/tD_n) e_i)
! 458: m := C.dd#(n0+1);
! 459: r := numgens C#(n0+1);
! 460: sourceDeg := degrees C#(n0+1);
! 461: sourceGens := {};
! 462: if explicitFlag then sourceMat := map(W^0,W^0,0);
! 463: i = 0;
! 464: while i < r do (
! 465: -- Find generators of the current source
! 466: -- "F_k1(D_n/tD_n)^r/F_k0(D_n/tD_n)^r"
! 467: -- as a left D_m module.
! 468: -- They have the form { 含prod_{i=1}^n D_i^{a_i} }.
! 469: tempExps = findExps(posWeights, k0-sourceDeg#i#0,
! 470: k1-sourceDeg#i#0);
! 471: tempGens = apply(tempExps, j -> posVars^j);
! 472: sourceGens = append(sourceGens, tempGens );
! 473: if explicitFlag then (
! 474: if tempGens == {} then (
! 475: sourceMat = directSum(sourceMat, compress matrix{{0_W}}); )
! 476: else (
! 477: sourceMat = directSum(sourceMat, matrix{tempGens}); );
! 478: );
! 479: i = i+1;
! 480: );
! 481: sourceSize := sum(sourceGens, i -> #i);
! 482:
! 483:
! 484: -- MAKE THE DIFFERENTIAL AS MATRIX OF D_m MODULES
! 485: if sourceSize == 0 and targetSize == 0 then (
! 486: oldDiff := map(resW^0,resW^0,0) )
! 487: else if sourceSize == 0 then ( oldDiff =
! 488: compress matrix toList(targetSize:{0_resW}) )
! 489: else if targetSize == 0 then ( oldDiff =
! 490: transpose compress matrix toList(sourceSize:{0_resW}) )
! 491: else (
! 492: -- For each generator j = 含prod_i D_i^{a_i}, compute its image
! 493: -- j*m_i as an element of the RHS. Get a matrix of image vectors.
! 494: imageMat := matrix join toSequence apply( r, a ->
! 495: apply(sourceGens#a, b -> substitute(b*m_a, resSub) ) );
! 496: -- differentiate with respect to targetGens
! 497: oldDiff = transpose compress matrix toList(sourceSize:{0_W});
! 498: i = 0;
! 499: -- compute the induced image
! 500: while i < s do (
! 501: if targetGens#i =!= {}
! 502: then oldDiff = oldDiff || substitute(
! 503: contract(transpose matrix{targetGens#i}, imageMat^{i}),
! 504: diffSub);
! 505: i = i+1;
! 506: );
! 507: oldDiff = WtoresW oldDiff;
! 508: );
! 509: if member(RestrictComplex, output) then (
! 510: restrictComplex#n0 = resW^(rank target oldDiff);
! 511: restrictComplex#(n0+1) = resW^(rank source oldDiff);
! 512: restrictComplex.dd#(n0+1) = map(restrictComplex#n0,
! 513: restrictComplex#(n0+1), oldDiff);
! 514: );
! 515: if member(Cycles, output) or member(Boundaries, output)
! 516: or member(GenCycles, output) then (
! 517: newKernel := zeroize2 mingens kernel oldDiff;
! 518: -- newKernel := zeroize2 gens kernel oldDiff;
! 519: explicitKernel := compress (sourceMat * resWtoW(newKernel));
! 520: );
! 521: );
! 522: pInfo(2, "含t含t含t Rank = " | sourceSize | "含t time = " | tInfo | " seconds" );
! 523:
! 524: -- DO THE COMPUTATION IN HIGHER COHOMOLOGICAL DEGREES
! 525: s = r;
! 526: targetGens = sourceGens;
! 527: targetSize = sourceSize;
! 528: if explicitFlag then targetMat = sourceMat;
! 529: targetMat = sourceMat;
! 530: currDeg := n0 + 2;
! 531: --newKernel = 0;
! 532: --explicitKernel = 0;
! 533: while currDeg <= n1 and C#?(currDeg) do (
! 534: -- MAKE THE NEXT SOURCE MODULE
! 535: -- "sourceGens" is a list of r lists, where the i-th list contains
! 536: -- the monomial generators {dx^a} (as left D_m-module) of
! 537: -- F_k1[u_i]((D_n/tD_n) e_i) / F_k0[u_i]((D_n/tD_n) e_i)
! 538: pInfo(2, "含t Degree " | currDeg | "...");
! 539: tInfo = toString first timing (
! 540: r = numgens C#currDeg;
! 541: m = C.dd#currDeg;
! 542: sourceDeg = degrees C#(currDeg);
! 543: sourceGens = {};
! 544: if explicitFlag then sourceMat = map(W^0,W^0,0);
! 545: i = 0;
! 546: while i < r do (
! 547: -- Find generators of the current source
! 548: -- "F_k1(D_n/tD_n)^r/F_k0(D_n/tD_n)^r"
! 549: -- as a left D_m module.
! 550: -- They have the form { 含prod_{i=1}^n D_i^{a_i} }.
! 551: tempExps = findExps(posWeights, k0-sourceDeg#i#0,
! 552: k1-sourceDeg#i#0);
! 553: tempGens = apply(tempExps, j -> posVars^j);
! 554: sourceGens = append(sourceGens, tempGens );
! 555: if explicitFlag then (
! 556: if tempGens == {} then (
! 557: sourceMat = directSum(sourceMat, compress matrix{{0_W}}); )
! 558: else (
! 559: sourceMat = directSum(sourceMat, matrix{tempGens}); );
! 560: );
! 561: i = i+1;
! 562: );
! 563: sourceSize = sum(sourceGens, i -> #i);
! 564:
! 565: -- MAKE THE NEXT DIFFERENTIAL OF D_m MODULES
! 566: if sourceSize == 0 and targetSize == 0 then (
! 567: newDiff := map(resW^0,resW^0,0) )
! 568: else if sourceSize == 0 then ( newDiff =
! 569: compress matrix toList(targetSize:{0_resW}) )
! 570: else if targetSize == 0 then ( newDiff =
! 571: transpose compress matrix toList(sourceSize:{0_resW}) )
! 572: else (
! 573: -- For each generator j = 含prod_i D_i^{a_i}, compute its image
! 574: -- j*m_i as an element of the RHS. Get a matrix of image vectors.
! 575: imageMat = matrix join toSequence apply( r, a ->
! 576: apply(sourceGens#a, b -> substitute(b*m_a, resSub) ) );
! 577: -- differentiate with respect to targetGens
! 578: newDiff = transpose compress matrix toList(sourceSize:{0_W});
! 579: i = 0;
! 580: while i < s do (
! 581: if targetGens#i =!= {}
! 582: then newDiff = newDiff || substitute(
! 583: contract(transpose matrix{targetGens#i}, imageMat^{i}),
! 584: diffSub);
! 585: i = i+1;
! 586: );
! 587: newDiff = WtoresW newDiff;
! 588: );
! 589:
! 590: -- UPDATE THE OUTPUT LIST
! 591: if member(RestrictComplex, output) then (
! 592: restrictComplex#currDeg = resW^(rank source newDiff);
! 593: restrictComplex.dd#currDeg = map(restrictComplex#(currDeg-1),
! 594: restrictComplex#currDeg, newDiff);
! 595: );
! 596: if member(HomologyModules, output) then (
! 597: tempHomology := homology(zeroize2 oldDiff, zeroize2 newDiff);
! 598: if tempHomology =!= null then
! 599: tempHomology = cokernel Dprune presentation tempHomology;
! 600: if tempHomology === null then tempHomology = resW^0;
! 601: homologyList = append(homologyList,
! 602: (currDeg-1) => tempHomology);
! 603: );
! 604:
! 605: -- MAKE EXPLICIT COHOMOLOGY CLASSES
! 606: if member(Cycles, output) or member(Boundaries, output) or
! 607: member(GenCycles, output) then (
! 608: oldImage := zeroize2 mingens image newDiff;
! 609: -- oldImage := zeroize2 gens image newDiff;
! 610: if member(GenCycles, output) then (
! 611: if #otherVars == 0 then (
! 612: explicitList = append(explicitList,
! 613: (currDeg-1) => targetMat *
! 614: resWtoW(mingens subquotient(newKernel, oldImage))) )
! 615: else (
! 616: explicitList = append(explicitList,
! 617: (currDeg-1) => explicitKernel);
! 618: );
! 619: );
! 620: if member(Cycles, output) then
! 621: kerList = append(kerList, (currDeg-1) => explicitKernel);
! 622: if member(Boundaries, output) then (
! 623: explicitImage := compress (targetMat * resWtoW(oldImage));
! 624: imList = append(imList, (currDeg-1) => explicitImage);
! 625: );
! 626: newKernel = zeroize2 mingens kernel newDiff;
! 627: -- newKernel = zeroize2 gens kernel newDiff;
! 628: explicitKernel = compress (sourceMat * resWtoW(newKernel));
! 629: );
! 630:
! 631: -- PREPARE FOR NEXT ITERATION
! 632: s = r;
! 633: targetGens = sourceGens;
! 634: targetSize = sourceSize;
! 635: if explicitFlag then targetMat = sourceMat;
! 636: oldDiff = newDiff;
! 637: currDeg = currDeg+1;
! 638: );
! 639: pInfo(2, "含t含t含t Rank = " | targetSize | "含t time = " | tInfo | " seconds");
! 640: );
! 641: );
! 642: pInfo(2, "含t Finished...");
! 643: pInfo(2, "含t含t含t Total time = " | tInfo | " seconds");
! 644: pInfo(2, " ");
! 645: );
! 646:
! 647: -- OUTPUT FORMAT
! 648: if member(HomologyModules, output) then outputList = append(
! 649: outputList, HomologyModules => hashTable homologyList);
! 650: if member(GenCycles, output) then outputList = append(
! 651: outputList, GenCycles => hashTable explicitList);
! 652: if member(Cycles, output) then outputList = append(
! 653: outputList, Cycles => hashTable kerList);
! 654: if member(Boundaries, output) then outputList = append(
! 655: outputList, Boundaries => hashTable imList);
! 656: if member(BFunction, output) then outputList = append(
! 657: outputList, BFunction => factorBFunction b);
! 658: if member(VResolution, output) then outputList = append(
! 659: outputList, VResolution => C);
! 660: if member(RestrictComplex, output) then outputList = append(
! 661: outputList, RestrictComplex => restrictComplex);
! 662: if member(ResToOrigRing, output) then outputList = append(
! 663: outputList, ResToOrigRing => resWtoW);
! 664: if member(TransferMap,output) then
! 665: if Resolution =!= null then
! 666: outputList = append(outputList, TransferMap => Resolution#TransferMap);
! 667:
! 668: hashTable outputList
! 669: )
! 670:
! 671:
! 672: twistedLogCohomology2 = method( Options => {Strategy => Schreyer} );
! 673: twistedLogCohomology2 (List, List) := options -> (f, a) -> (
! 674: k := #f;
! 675: if k != #a
! 676: then error "number of factors is not equal number of parametors";
! 677:
! 678: i := 0;
! 679: prodf := 1;
! 680: while i < k do (
! 681: prodf = prodf * f_i;
! 682: i = i+1;
! 683: );
! 684:
! 685: R = ring prodf;
! 686: if class R =!= PolynomialRing
! 687: then error "expected element of a polynomial ring";
! 688: if R.monoid.Options.WeylAlgebra =!= {}
! 689: then error "expected element of a commutative polynomial ring";
! 690: v := generators(R);
! 691: n := #v;
! 692: if n != 2
! 693: then error "we compute 2-dimention logCohomology";
! 694: d := {diff(v#1,prodf),-diff(v#0,prodf),prodf};
! 695: syzf := kernel matrix{d};
! 696: m := numgens syzf;
! 697: msyz := gens syzf;
! 698: W := makeWeylAlgebra( R, SetVariables=>false);
! 699: phi := map(W,R);
! 700: vw := generators(W);
! 701: i = 0; ell := 0; j := 0, sum := 0;
! 702: op := {};
! 703: while i<m do (
! 704: -- adjoint of differential operators
! 705: j = 0; sum = 0;
! 706: while j < k do (
! 707: sum = sum + a_j * (msyz_(1,i) * diff(v#0,f_j) - msyz_(0,i) * diff(v#1,f_j)) // f_j;
! 708: j = j+1;
! 709: );
! 710: ell = -phi msyz_(1,i) * vw_(2) + phi msyz_(0,i) * vw_(3) - phi msyz_(2,i) + phi sum;
! 711: op = join(op,{ell});
! 712: i = i+1;
! 713: );
! 714: fop := apply(op,Fourier);
! 715:
! 716: outputRequest := {GenCycles, HomologyModules,
! 717: VResolution, BFunction, Explicit, TransferMap};
! 718: M := cokernel matrix{fop};
! 719: outputTable := computeRestriction2(M, {1,1}, -1, n+1,
! 720: outputRequest, options);
! 721: outputList := {BFunction => outputTable#BFunction};
! 722: if outputTable#VResolution =!= null then
! 723: outputList = append(outputList,
! 724: VResolution => FourierInverse outputTable#VResolution)
! 725: else outputList = append(outputList, VResolution => null);
! 726:
! 727: outputList = outputList | {
! 728: PreCycles => hashTable apply(toList(0..n),
! 729: i -> (n-i => FourierInverse outputTable#GenCycles#i)),
! 730: CohomologyGroups => hashTable apply(toList(0..n),
! 731: i -> (n-i => outputTable#HomologyModules#i)) };
! 732:
! 733: intTable := hashTable outputList;
! 734:
! 735: -- make basis of H^0
! 736: -- dim H^0 is -1 or 0
! 737: basis0 := 0;
! 738: if dim intTable#CohomologyGroups#0 == 0 then (
! 739: sumx := 0; sumy := 0; j = 0;
! 740: while j < k do (
! 741: sumx = sumx + a_j * prodf * diff(v#0,f_j) // f_j;
! 742: sumy = sumy + a_j * prodf * diff(v#1,f_j) // f_j;
! 743: j = j+1;
! 744: );
! 745: I := ideal(phi prodf * vw_(2) + phi sumx, phi prodf * vw_(3) + phi sumy );
! 746: basis0 = matrix{PolySols I};
! 747: )
! 748: else basis0 = intTable#PreCycles#0;
! 749:
! 750: -- make basis of H^1
! 751: if outputTable#?TransferMap then(
! 752: transfermap := transpose Dtransposition FourierInverse outputTable#TransferMap;
! 753:
! 754: mat := transpose(intTable#PreCycles#1)*transfermap;
! 755: l := numRows mat;
! 756: mat = mutableMatrix mat;
! 757: i = 0;
! 758: while i < l do(
! 759: j = 0;
! 760: while j < m do(
! 761: mat_(i,j) = substitute(mat_(i,j),{vw_(2) => 0, vw_(3) => 0});
! 762: j = j+1;
! 763: );
! 764: i = i+1;
! 765: );
! 766: mat = matrix mat;
! 767: msyz = phi msyz;
! 768:
! 769: Basis1 := {};
! 770: i = 0;
! 771: while i < l do(
! 772: base = 0; j = 0;
! 773: while j < m do(
! 774: base = base + mat_(i,j) * (msyz_(0,j) * vw_(2) + msyz_(1,j) * vw_(3));
! 775: j = j+1;
! 776: );
! 777: Basis1 = join(Basis1,{base});
! 778: i = i+1;
! 779: );
! 780: Basis1 = matrix{Basis1};
! 781: )
! 782: else Basis1 = intTable#PreCycles#1;
! 783:
! 784: if intTable#VResolution =!= null then (
! 785: createDpairs W;
! 786: Omega := Dres ideal W.dpairVars#1;
! 787: Basis := hashTable {0 => basis0,
! 788: 1 => Basis1,
! 789: 2 => intTable#PreCycles#2 * vw_(2)*vw_(3)};
! 790: outputList = outputList | {LogBasis =>
! 791: Basis, OmegaRes => Omega};
! 792: );
! 793:
! 794: hashTable outputList
! 795: )
! 796:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to twistedLogCohomology.m2 CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / math-misc