[BACK]Return to twistedLogCohomology2.m2 CVS log [TXT][DIR] Up to [local] / OpenXM / src / math-misc

Annotation of OpenXM/src/math-misc/twistedLogCohomology2.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: )

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>