Skip to content

Commit

Permalink
Merge pull request #2940 from pplt/localCohom
Browse files Browse the repository at this point in the history 
Improvements in local cohomology computations
  • Loading branch information
@DanGrayson
DanGrayson authored Oct 10, 2023
2 parents 5ccd34f + 1716a0f commit 00c0bb9
Show file tree
Hide file tree
Showing 9 changed files with 104 additions and 100 deletions.
26 changes: 13 additions & 13 deletions M2/Macaulay2/packages/BernsteinSato.m2
View file
Original file line number Diff line number Diff line change
Expand Up @@ -271,26 +271,26 @@ export {
-- Tests
--------------------------------------------------------------------------------

load "BernsteinSato/TST/tests.m2"
load "./BernsteinSato/TST/tests.m2"

--------------------------------------------------------------------------------
-- Documentation
--------------------------------------------------------------------------------

beginDocumentation()

load "BernsteinSato/DOC/main.m2"
load "BernsteinSato/DOC/DHom.m2"
load "BernsteinSato/DOC/Dlocalize.m2"
load "BernsteinSato/DOC/Drestriction.m2"
load "BernsteinSato/DOC/WeylClosure.m2"
load "BernsteinSato/DOC/annFs.m2"
load "BernsteinSato/DOC/bFunctions.m2"
load "BernsteinSato/DOC/localCohom.m2"
load "BernsteinSato/DOC/intersectionCohom.m2"
load "BernsteinSato/DOC/multiplierIdeals.m2"
load "BernsteinSato/DOC/paco-anton-paper.m2"
load "BernsteinSato/DOC/other.m2"
load "./BernsteinSato/DOC/main.m2"
load "./BernsteinSato/DOC/DHom.m2"
load "./BernsteinSato/DOC/Dlocalize.m2"
load "./BernsteinSato/DOC/Drestriction.m2"
load "./BernsteinSato/DOC/WeylClosure.m2"
load "./BernsteinSato/DOC/annFs.m2"
load "./BernsteinSato/DOC/bFunctions.m2"
load "./BernsteinSato/DOC/localCohom.m2"
load "./BernsteinSato/DOC/intersectionCohom.m2"
load "./BernsteinSato/DOC/multiplierIdeals.m2"
load "./BernsteinSato/DOC/paco-anton-paper.m2"
load "./BernsteinSato/DOC/other.m2"

--------------------------------------------------------------------------------

Expand Down
18 changes: 3 additions & 15 deletions M2/Macaulay2/packages/BernsteinSato/DOC/annFs.m2
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,38 +2,26 @@ doc ///
Key
AnnFs
(AnnFs, RingElement)
(AnnFs, List)
Headline
differential annihilator of a polynomial in a Weyl algebra
Usage
AnnFs(f)
AnnFs(L)
Inputs
f:RingElement
polynomial in a Weyl algebra $D$
L:List
of polynomials in the Weyl algebra $D$
Outputs
:Ideal
the differential annihilator of $f$ in $D[s]$, for a new variable $s$, or the differential
annihilator of $L$ in $D[t_0,..,t_0,dt_0,..,dt_k]$, where $L$ is a list of $k+1$ elements
and the $t_i$ are new variables.
the differential annihilator of $f$ in $D[s]$, for a new variable $s$.
Description
Text
This routine computes the ideal of the differential annihilator of a polynomial or list of
polynomials in a Weyl algebra $D$. This ideal is a left ideal of the ring $D[s]$
or $D[t_0,..,t_k,dt_0,..,dt_k]$. More
This routine computes the ideal of the differential annihilator of a polynomial. This ideal is a left ideal of the ring $D[s]$. More
details can be found in
[@HREF("https://mathscinet.ams.org/mathscinet/pdf/1734566.pdf","SST")@, Chapter 5].
The computation in the case of the element $f$ is via Algorithm 5.3.6, and the computation
in the case of the list $L$ is via the Algorithm 5.3.15.
The computation in the case of the element $f$ is via Algorithm 5.3.6.
Example
makeWA(QQ[x,y])
f = x^2+y
AnnFs f
Example
makeWA(QQ[x,y,z])
L = {x^3,y+5*z}
Caveat
Must be over a ring of characteristic $0$.
///
Expand Down
38 changes: 22 additions & 16 deletions M2/Macaulay2/packages/BernsteinSato/DOC/localCohom.m2
View file
Original file line number Diff line number Diff line change
@@ -1,36 +1,42 @@
document {
Key => [localCohom,Strategy],
Headline => "specify strategy for local cohomology",
"This option together with ", TO "LocStrategy", " determines a strategy for ",
"There are two main strategies, Walther and OaTa. If the user selects Walther, which is the default, then ", TO "LocStrategy", " determines the localization strategy for ",
TT "localCohom(...Ideal...)", " and ", TT "localCohom(...Ideal, Module...)", ".",
UL {
{BOLD "Walther", " -- the algorithm of U. Walther that uses Cech complex."},
{BOLD "LocStrategy => null",
" -- used only for ", TT "localCohom(...Ideal...)",
", localizations are done by straitforward computation of
annihilators and b-polynomials as described in [1]."},
{BOLD "LocStrategy => OaTaWa",
" -- localizations are done following Oaku-Takayama-Walther method."},
{BOLD "LocStrategy => Oaku",
" -- localizations are done following Oaku's algorithm."},
{BOLD "OaTa", " -- restriction algorithm is used,
which is due to T. Oaku and N. Takayama [2]"}
UL {
{BOLD "LocStrategy => null",
" -- used only for ", TT "localCohom(...Ideal...)",
", localizations are done by straigthforward computation of
annihilators and b-polynomials as described in [1]."},
{BOLD "LocStrategy => OaTaWa",
" -- localizations are done following Oaku-Takayama-Walther method [2]."},
{BOLD "LocStrategy => Oaku",
" -- localizations are done following Oaku's algorithm."},
},
{BOLD "OaTa", " -- restriction from the graph embedding is used,
which is due to T. Oaku and N. Takayama [3]. See ", TO "Drestriction", "."}
},
--Caveat => {"When WaltherOTW strategy is used the error 'Bad luck!'
Caveat => {"localCohom(...Ideal, Module...) with the default strategy computes presentations for all the terms in the Cech complex regardless of the requested homological degrees. All strategies use the given generators of the ideal; the user is advised to call ", TO "mingens", " before calling localCohom."},
--Caveat => {"When OaTaWa strategy is used the error 'Bad luck!'
--may appear. This means your are not a lucky individual...
--The glitch is due to the fact that the localizations are iterated
--for this particular strategy; it was resolved for WaltherOaku,
--a strategy that considers everyone lucky."
--},
"For detailed description of the algorithms see",
UL {
{BOLD "[1]", "U. Walther, ",
{BOLD "[1] ", "Walther, ",
EM "Algorithmic computation of local cohomology
modules and the local cohomological dimension of algebraic
varieties (JPAA (139), 1999.)"
},
{BOLD "[2]", "Oaku, Takayama",
EM "Algorithms for D-modules..."
{BOLD "[2] ", "Oaku, Takayama, Walther, ",
EM "A Localization Algorithm for D-modules (J. Symbolic Computation (29), 2000.)"
},
{BOLD "[3] ", "Oaku, Takayama, ",
EM "Algorithms for D-modules -- restriction, tensor product, localization, and local cohomology groups (JPAA (156), 2001.)"
}
}
}
Expand All @@ -40,7 +46,7 @@ document {
document {
Key => [localCohom,LocStrategy],
Headline => "specify localization strategy for local cohomology",
"See ", TO [localCohom,Strategy]
"These strategies determine how presentations of localization in the Cech complex are calculated when selecting Walther's strategy. See ", TO [localCohom,Strategy]
}
document {
Key => Walther,
Expand Down
9 changes: 8 additions & 1 deletion M2/Macaulay2/packages/BernsteinSato/annFs.m2
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
-- Copyright 1999-2008 by Anton Leykin and Harrison Tsai

-- AnnFs uses the algorithm of Oaku-Takayama which involves adding 4 extra variables to the Weyl algebra and performing two eliminations. TODO: Implement the Briancon-Maisonobe algorithm. This would require support for working with Ore algebras (for instance [dt, s] = s) in the Macaulay2 core.

---------------------------------------------------------------------------------
AnnFs = method()
AnnFs RingElement := Ideal => f -> (
Expand All @@ -24,7 +26,11 @@ AnnFs RingElement := Ideal => f -> (

AnnIFs(ideal dpV#1, f)
);
AnnFs List := Ideal => F -> (

-- This function used to be called AnnFs but was incomplete, it only computes the Malgrange ideal.

MalgrangeIdeal = method()
MalgrangeIdeal List := Ideal => F -> (
-- Input: F = {f_1,...,f_r}, a list of polynomials in n variables
-- (f_i has to be an element of A_n, the Weyl algebra).
-- Output: Ann f_1^{s_1}...f_r^{s_r}, an ideal in A_n<t_1,..., t_r,dt_1,...,dt_r>.
Expand Down Expand Up @@ -57,6 +63,7 @@ AnnFs List := Ideal => F -> (
sub(DXj,WT) + sum(r, i->sub(DXj*F#i-F#i*DXj,WT)*WT_(i+r+2*n) )
))
);

------------------------------------------------------------------------------
--This needs documentation.
AnnIFs = method()
Expand Down
11 changes: 9 additions & 2 deletions M2/Macaulay2/packages/BernsteinSato/globalBFunction.m2
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,9 +71,16 @@ globalBFunctionIdeal RingElement := RingElement => o -> f -> (
-- REDUCED global b-function: b(s)/(s+1)
------------------------------------------------------------------
globalRB := method()
globalRBAnnFs = method()

globalRB (RingElement,Boolean) := RingElement => (f,isRed) -> (
W := ring f;
AnnI := AnnFs f;
globalRBAnnFs(f, AnnI, isRed);
)

globalRBAnnFs (RingElement,Ideal,Boolean) := RingElement => (f,AnnI,isRed) -> (
W := ring f;
Ws := ring AnnI;
ns := numgens Ws;
createDpairs W;
Expand Down Expand Up @@ -193,7 +200,7 @@ generalB (List, RingElement) := RingElement => o->(F,g) -> (
g = sub(g,D);
);
r := #F;
AnnI := AnnFs F;
AnnI := MalgrangeIdeal F;
DY := ring AnnI;
K := coefficientRing DY;
n := numgens DY // 2 - r; -- DY = k[x_1,...,x_n,t_1,...,t_r,dx_1,...,dx_n,dt_1,...,dt_r]
Expand Down Expand Up @@ -258,7 +265,7 @@ generalBideal (List, RingElement) := RingElement => o->(F,g) -> (
-- g, a polynomial
-- Output: the generalized B-S ideal, an ideal of QQ[s_1..s_r]
r := #F;
AnnI := AnnFs F;
AnnI := MalgrangeIdeal F;
DY := ring AnnI;
K := coefficientRing DY;
n := numgens DY // 2 - r; -- DY = k[x_1,...,x_n,t_1,...,t_r,dx_1,...,dx_n,dt_1,...,dt_r]
Expand Down
2 changes: 1 addition & 1 deletion M2/Macaulay2/packages/BernsteinSato/localBFunction.m2
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ computeJf RingElement := Ideal => f -> (
K := coefficientRing D;

w := toList(n:0) | {1};
inIf := inw(AnnFs {f}, -w|w);
inIf := inw(MalgrangeIdeal {f}, -w|w);

Dt := ring inIf;
x := take(gens Dt,{0,n-1});
Expand Down
90 changes: 43 additions & 47 deletions M2/Macaulay2/packages/BernsteinSato/localCohom.m2
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ localCohom(List, Ideal) := HashTable => o -> (l,I) -> (
-- Promote I to the Weyl algebra if it is not already there
if #(ring I).monoid.Options.WeylAlgebra == 0
then I = sub(I, makeWA ring I); -- TODO: what degrees are best for the differentials?

if (o.Strategy == Walther and o.LocStrategy === null)
then localCohomUli (l,I)
else (
Expand All @@ -36,66 +37,61 @@ localCohom(List, Ideal) := HashTable => o -> (l,I) -> (
localCohomUli = (l, I) -> (
-- error checking to be added
-- I is assumed to be an ideal in WA generated by polynomials
f := first entries gens I;

f := I_*;
r := #f;
W := ring I;
subISets := select(subsets toList (0..r-1), s -> s =!= {});


n := (dim ring I) // 2; -- dimension ambient ring
L := sort unique flatten apply(l, x -> {x - 1, x, x + 1});
Lmaps := sort unique flatten apply(l, x -> {x - 1, x});

L = select(L, x -> x >= 0 and x <= min(r, n)); -- cohomological degrees
Lmaps = select(Lmaps, x -> x >= 0 and x < min(r, n)); -- source degrees


subISets := new HashTable from apply(L, x -> x => subsets(toList (0..r-1), x));

-- Step1.
-- Calculate J^/delta( (F_/theta)^s ) and b^/delta_(F_/theta)(s) for all /theta
pInfo(1, "localCohom: Computing b-functions and annihilators...");
J := new MutableHashTable;
bF := new MutableHashTable;
scan(subISets, theta->(
Ftheta := product(theta, i->f#i);
J#theta = AnnFs(Ftheta);
bF#theta = globalBFunction(Ftheta);
));

prodS := new HashTable from apply(flatten values subISets, theta -> theta => (product(theta, i -> f#i))_W);
J := applyValues(prodS, f -> (f, AnnFs f, true));
bF := applyValues(J, globalRBAnnFs);
J = applyValues(J, x -> x#1); -- Now J is a hash table of annihilators

-- Step 2.
-- a = min integer root of all bF-s
a := min flatten (subISets / (theta -> getIntRoots(bF#theta)));
a := min append(flatten ((values bF) / getIntRoots), -1);
if a == infinity then a = 0;
pInfo(666, {"BEST POWER = " , a});
-- Substitute s = a in J-s
scan(subISets, theta -> (
AforS := map(W, ring J#theta, vars W | matrix {{a_W}});
J#theta = AforS J#theta;
));


Specialize := X -> (Phi := map(W, ring X, vars W | matrix {{a_W}}); Phi X);
J = applyValues(J, Specialize);

-- Step 3.
-- Compute the Cech complex
C := new MutableList from toList ((r+1):());
M := new MutableList from toList (r:());
pInfo(1, "Constructing Cech complex...");
C#0 = directSum { {} => W^1 / ideal W.dpairVars#1 };
apply(toList(1..r), k->(
C#k = directSum (select(subISets, u -> #u == k) / (theta -> (
theta => W^1 / J#theta
)));

M#(k-1) = map (C#k, C#(k-1), (i,j) -> (
i0 := (indices C#(k-1))_j;
j0 := (indices C#k)_i;
if isSubset(i0, j0)
then (
l := first toList (set j0 - set i0);
(-1)^(position(j0, u -> u == l)) * (f_l)^(-a)
)
else 0_W
)
);
));
C := applyValues(subISets, S -> directSum apply(S, theta -> (theta => W^1 / J#theta)));
M := new HashTable from apply(Lmaps, k -> k => map (C#(k+1), C#k, (i,j) -> (
i0 := (indices C#k)_j;
j0 := (indices C#(k+1))_i;
if isSubset(i0, j0) then
(
m := first toList (set j0 - set i0);
(-1)^(position(j0, u -> u == m)) * (f_m)^(-a)
)
else 0_W)
));

-- Step 4.
-- Compute the homology of the complex
pInfo(1, "Computing cohomology...");
ret := new HashTable from apply(l, k -> k=>
if k<0 or k>r then W^0
else if k==0 then kernel M#0
else if k==r then cokernel M#(r-1)
else homology(M#k, M#(k-1))
);
ret
new HashTable from apply(l, k -> k =>
if k < 0 or k > min(r, n) then W^0
else if k == 0 then kernel M#0
else if k == min(r, n) then cokernel M#(min(r, n)-1)
else homology(M#k, M#(k-1))
)
);

----------------------------------------------------------------------------------------
Expand Down Expand Up @@ -391,7 +387,7 @@ localCohomOT(Ideal, Ideal) := (I, J) -> (
if not J.?quotient then J.quotient = (ring J)^1/J;
localCohomOT(I, J.quotient)
)
localCohomOT(Ideal, Module) := (I, M) -> computeLocalCohomOT(I, M, 0, numgens I)
localCohomOT(Ideal, Module) := (I, M) -> computeLocalCohomOT(I, M, 0, min(numgens I, (dim ring I)//2))

localCohomOT(List, Ideal, Module) := (l, I, M) -> (
locOut := computeLocalCohomOT(I, M, min l, max l);
Expand Down
4 changes: 2 additions & 2 deletions M2/Macaulay2/packages/BernsteinSato/multiplierIdeals.m2
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ multiplierIdeal (Ideal, List) := o -> (a,cs) -> (
F = apply(F, f->sub(f,D));
);
r := #F;
I0 := AnnFs F;
I0 := MalgrangeIdeal F;
DY := ring I0;
K := coefficientRing DY;
n := numgens DY // 2 - r; -- DY = k[x_1,...,x_n,t_1,...,t_r,dx_1,...,dx_n,dt_1,...,dt_r]
Expand Down Expand Up @@ -163,7 +163,7 @@ lct Ideal := RingElement => o -> I -> (
F := (sub(I,W))_*;
if o.Strategy === ViaBFunction then (
w := toList (numgens ring I:0) | toList(#F:1);
b := bFunction(AnnFs F, w);
b := bFunction(MalgrangeIdeal F, w);
S := ring b;
r := numgens I;
-- lct(I) = min root of b(s-r)
Expand Down
6 changes: 3 additions & 3 deletions M2/Macaulay2/packages/QuillenSuslin.m2
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ newPackage(
Date => "May 10, 2013",
Authors => {
{Name => "Brett Barwick", Email => "[email protected]", HomePage => "http://faculty.uscupstate.edu/bbarwick/"},
{Name => "Branden Stone", Email => "bstone@adelphi.edu", HomePage => "http://math.adelphi.edu/~bstone/"}
{Name => "Branden Stone", Email => "[email protected].edu", HomePage => "http://bstone.github.io/"}
},
Headline => "the Quillen-Suslin algorithm for bases of projective modules",
Keywords => {"Commutative Algebra"},
Expand Down Expand Up @@ -1457,7 +1457,7 @@ monicPolySubs(RingElement,List) := opts -> (f,varList) -> (
print "The element had degree zero in the last variable.";
degZeroSub = mutableMatrix vars R;
degZeroSub = columnSwap(degZeroSub,last usedVarPosition,lastVarPosition);
f = sub(f,degZeroSub); -- Interchange variables so that last varList is involved in f. Now f has positive degree in last varList.
f = sub(f,matrix degZeroSub); -- Interchange variables so that last varList is involved in f. Now f has positive degree in last varList.
);

-- Now we enter the general algorithm.
Expand Down Expand Up @@ -1495,7 +1495,7 @@ monicPolySubs(RingElement,List) := opts -> (f,varList) -> (
if degZeroSub =!= null then (
print("degZeroSub: "|toString(degZeroSub));
tempSub = columnSwap(tempSub,last usedVarPosition,lastVarPosition);
tempInvSub = sub(tempInvSub,degZeroSub);
tempInvSub = sub(matrix tempInvSub,matrix degZeroSub);
);

return (matrix tempSub,matrix tempInvSub);
Expand Down

0 comments on commit 00c0bb9

Please sign in to comment.