Add experimental support for wreath Macdonald polynomials

docs/oscar_references.bib

@@ -232,6 +232,17 @@ @Article{BES23
   doi           = {10.4171/jems/1350}
 }
 
+@Article{BF14,
+  author        = {Bezrukavnikov, Roman and Finkelberg, Michael},
+  title         = {Wreath Macdonald polynomials and the categorical McKay correspondence},
+  journal       = {Camb. J. Math.},
+  volume        = {2},
+  number        = {2},
+  pages         = {163--190},
+  year          = {2014},
+  doi           = {10.4310/CJM.2014.v2.n2.a1}
+}
+
 @Misc{BG24,
   author        = {Simon Brandhorst and Víctor González-Alonso},
   title         = {527 elliptic fibrations on Enriques surfaces},
@@ -1560,6 +1571,17 @@ @Article{HT17
   primaryclass  = {hep-th}
 }
 
+@Article{Hai02,
+  author        = {Haiman, Mark},
+  title         = {Combinatorics, symmetric functions and Hilbert schemes},
+  journal       = {Current developments in mathematics},
+  volume        = {2002},
+  number        = {1},
+  publisher     = {International Press of Boston},
+  pages         = {39--111},
+  year          = {2002}
+}
+
 @Book{Har77,
   author        = {Hartshorne, Robin},
   title         = {Algebraic Geometry},
@@ -2291,6 +2313,13 @@ @Article{OMdCS00
   doi           = {10.1006/jsco.1999.0413}
 }
 
+@Article{OS23,
+  author        = {Orr, Daniel and Shimozono, Mark},
+  title         = {Wreath Macdonald polynomials, a survey},
+  journal       = {arXiv preprint arXiv:2308.12166},
+  year          = {2023}
+}
+
 @Book{Oxl11,
   author        = {Oxley, James},
   title         = {Matroid theory},

experimental/WreathMacdonaldpols/README.md

@@ -0,0 +1,3 @@
+# Wreath Macdonald polynomials
+
+The main purpose of this code is to compute the wreath Macdonald polynomials in the Schur basis.

experimental/WreathMacdonaldpols/docs/doc.main

@@ -0,0 +1,5 @@
+[
+   "Wreath Macdonald polynomials" => [
+      "introduction.md",
+   ]
+]

experimental/WreathMacdonaldpols/docs/src/introduction.md

@@ -0,0 +1,42 @@
+```@meta
+CurrentModule = Oscar
+CollapsedDocStrings = true
+DocTestSetup = Oscar.doctestsetup()
+```
+
+# Wreath Macdonald polynomials
+
+The existence, integrality and positivity of wreath Macdonald polynomials
+ has been conjectured by Haiman [Hai02](@cite) and proved by Bezrukavnikov
+ and Finkelberg [BF14](@cite). When ``r=1``, wreath Macdonald polynomials are
+ equal to the Haiman-Macdonald polynomials, used to prove the Macdonald positivity conjecture.
+
+Here we have implemented an algorithm computing the wreath Macdonald
+ polynomials as defined in the survey by Orr and Shimozono on this topic [OS23](@cite).
+
+Wreath Macdonald polynomials depend on two parameters. The first parameter is
+ an ``r``-multipartition of ``n``. The second parameter is an element of the affine Weyl group
+ of type ``A^{(1)}_{r-1}`` which is isomorphic to the semi-direct product of the finite Weyl group
+ of type ``A_{r-1}`` (the symmetric group on ``r`` letters) and of the coroot lattice of type ``A_{r-1}``.
+ The element of the coroot lattice is given in the canonical basis. It is then the sublattice
+of ``\mathbb{Z}^r`` of elements summing up to zero.
+
+```@docs
+wreath_macdonald_polynomial
+wreath_macdonald_polynomials
+```
+
+Compare the following computation with Example 3.15 in [OS23](@cite).
+
+```jldoctest
+julia> collect(multipartitions(1,3))
+3-element Vector{Multipartition{Int64}}:
+ Partition{Int64}[[], [], [1]]
+ Partition{Int64}[[], [1], []]
+ Partition{Int64}[[1], [], []]
+
+julia> wreath_macdonald_polynomials(1,3,cperm(1:3),[0,1,-1])[[3, 2, 1],[3, 2, 1]]
+[1   q^2     q]
+[1     t     q]
+[1     t   t^2]
+```

experimental/WreathMacdonaldpols/src/WreathMacdonaldpols.jl

@@ -0,0 +1,255 @@
+# I would like to credit Dario Mathiä who produced an initial version of the following code
+
+export wreath_macdonald_polynomials, wreath_macdonald_polynomial
+
+# Tools
+
+# Computes the b-invariant
+function b1(lambda::Vector{Int})
+  return sum((i-1)*l_i for (i, l_i) in enumerate(lambda); init=0)
+end
+
+# Computes the b-invariant of a multipartition
+function b_inv(lbb::Multipartition)
+  return b1(map(sum,lbb)) + length(lbb)*sum(lambda -> b1(data(lambda)),lbb)
+end
+
+function beta_number_to_partition(beta::Vector{Int})
+  lb=Int[]
+  for j in 1:length(beta)
+    nbholes=beta[j]-beta[1]-(j-1)
+    if nbholes >= 1
+      pushfirst!(lb, nbholes)
+    end
+  end
+  return partition(lb)
+end
+
+# Computes the core associated to the coroot
+#The charge of the core is equal to the coroot
+#in the canonical basis (e_i) of Z^I.
+function core(coroot::Vector{Int},r::Int)
+  beta=Int[]
+  m=minimum(coroot)
+  M=maximum(coroot)
+  beta=[k*r+(j-1) for k in m:M for j in 1:r if k <= coroot[j]]
+  return beta_number_to_partition(beta)
+end
+
+#Inspired by Sagemath's algorithm to obtain a partition from core and quotient
+#Be careful: it is not the usual quotient but the one defined in [Gordon 2008]
+function tau_om(lbb::Multipartition, wperm::PermGroupElem, coroot::Vector{Int})
+   r=length(lbb)
+   w0=perm([r-i for i in 0:(r-1)])  #cf Remark 3.5 Orr and Shimozono
+   lbb_perm=multipartition(permuted(lbb.mp,wperm*w0))
+   gamma=core(coroot,r)
+   lg=length(gamma)
+   k=r*maximum(length(lbb_perm[i]) for i in 1:r) + lg
+   v=[[gamma[i]-(i-1) for i in 1:lg]; [-i for i in lg:(k-1)]]
+   w=[[x for x in v if mod((x-i),r) == 0] for i in 1:r]
+   new_w=Int[]
+   for i in 1:r
+     lw=length(w[i])
+     lq=length(lbb_perm[i])
+     append!(new_w, w[i][1:lq] + r*lbb_perm[i])
+     append!(new_w, w[i][lq+1:lw])
+   end
+   sort!(new_w,rev=true)
+   new_w=[new_w[i]+(i-1) for i in 1:length(new_w)]
+   filter!(x-> x!=0,new_w)
+   return partition(new_w)
+end
+
+#Reoders the CharTable
+function reorder(charTirr::Matrix, Modules::Vector{Multipartition{Int}},  mps::Vector{Multipartition{Int}})
+  new_ord=indexin(mps,Modules)
+  return charTirr[new_ord,:]
+end
+
+#Tools from representation theory
+
+# Computes the fake degree of a multipartition cf. Ste89 Thm 5.3 and Prop. 3.3.2 Haiman cdm
+function fake_deg(lbb::Multipartition, Q::FracField{T}, var::T) where {T <: MPolyRingElem}
+  r=length(lbb)
+  n=sum(lbb)
+  res=Q(1)
+  for p in 1:n
+    res=res*(1-var^(r*p))
+  end
+  for lambda in lbb
+    for i in 1:length(lambda)
+      for j in 1:lambda[i]
+        hookfactor=(Q(1)-var^(r*(1+lambda[i]+count(>=(j),lambda)-j-i)))
+        res=res//hookfactor
+      end
+    end
+  end
+  res=res*var^b_inv(lbb)
+  return res
+end
+
+# computes C_Delta defined in PhD relation (1.45)
+function C_Delta(r::Int, n::Int, Q::FracField{T}, var::T, mps::Vector{Multipartition{Int}}) where {T <: MPolyRingElem}
+  charTable=character_table_complex_reflection_group(r,1,n)
+  charTirr=[charTable[i,j] for i in 1:nrows(charTable), j in 1:ncols(charTable)]
+  Modules=[multipartition([lbb...]) for lbb in class_parameters(charTable)]
+  l=length(Modules)
+  charTirr=reorder(charTirr,Modules,mps)
+  charTirrT=solve_init(matrix(Q,transpose(charTirr)))
+  rows=zero_matrix(Q,l,0)
+  for i in 1:l
+    col=zero_matrix(Q,l,1)
+    for j in 1:l
+      v=matrix(Q,l,1,map(k->charTirr[i,k]*charTirr[j,k],1:l))
+      x=solve(charTirrT,v,side=:right)
+      col=col+fake_deg(mps[j],Q,var)*x
+    end
+    rows=hcat(rows,col)
+  end
+  return rows
+end
+
+# computes the character of the simple representations of the Cherednik algebra.
+function C_L(r::Int, n::Int, Q::FracField{T}, var::T, mps::Vector{Multipartition{Int}}) where {T <: MPolyRingElem}
+  C_D=C_Delta(r,n,Q,var,mps)
+  d = [var^b_inv(mp)//fake_deg(mp,Q,var) for mp in mps]
+  diag=diagonal_matrix(Q,d)
+  return diag * C_D
+end
+
+function bigger_ord(lbb::Multipartition, wperm::PermGroupElem, coroot::Vector{Int}, mps::Vector{Multipartition{Int}})
+  r=length(lbb)
+  n=sum(lbb)
+  lbbquot=tau_om(lbb,wperm,coroot)
+  res=map(x->x,Iterators.filter(x-> dominates(tau_om(x,wperm,coroot),lbbquot),mps))
+  return res
+end
+
+function smaller_ord(lbb::Multipartition, wperm::PermGroupElem, coroot::Vector{Int}, mps::Vector{Multipartition{Int}})
+  r=length(lbb)
+  n=sum(lbb)
+  lbbquot=tau_om(lbb,wperm,coroot)
+  res=map(x->x,Iterators.filter(x-> dominates(lbbquot,tau_om(x,wperm,coroot)),mps))
+  return res
+end
+
+@doc raw"""
+    wreath_macdonald_polynomials(n::Int,
+                                 r::Int,
+                                 wperm::PermGroupElem,
+                                 coroot::Vector{Int};
+                                 parent::MPolyRing{<:QQAbFieldElem}=
+                                 polynomial_ring(abelian_closure(QQ)[1], [:q,:t];cached=true)[1])
+
+Given two integers `n` and `r` and an element of the affine Weyl group of type ``A^{(1)}_{r-1}``
+(seen as the semi-direct product of the symmetric group on `r` letters with the coroot lattice
+of the finite type ``A_{r-1}``), this function returns the square matrix of coefficients of the wreath
+ Macdonald polynomials associated with all multipartitions of size `n` and length `r` in the standard
+Schur basis of multisymmetric functions. Each row of this matrix is a wreath Macdonald polynomial.
+Here is an example of how to use it:
+
+```jldoctest
+julia> wreath_macdonald_polynomials(1,3,cperm(1:3),[0,1,-1])
+[t^2     t   1]
+[  q     t   1]
+[  q   q^2   1]
+```
+"""
+function wreath_macdonald_polynomials(n::Int,
+                                      r::Int,
+                                      wperm::PermGroupElem,
+                                      coroot::Vector{Int};
+                                      parent::MPolyRing{<:QQAbFieldElem}=
+                                      polynomial_ring(abelian_closure(QQ)[1], [:q,:t];cached=true)[1])
+
+  Q = fraction_field(parent)
+  q,t = gens(parent)
+
+  mps=collect(multipartitions(n,r))
+  l=length(mps)
+
+  c_L=C_L(r,n,Q,t,mps)
+  c_L_q=map_entries(f->numerator(f)(0,q)//denominator(f)(0,q),c_L)
+  c_L_tinv=map_entries(f->numerator(f)(0,1//t)//denominator(f)(0,1//t),c_L)
+
+  rows=[]
+  for lbb in mps
+    smallers=smaller_ord(lbb, wperm, coroot,mps) #tinv
+    biggers=bigger_ord(lbb, wperm, coroot,mps) #q
+    smaller_indices=sort!(indexin(smallers,mps))
+    bigger_indices=sort!(indexin(biggers,mps))
+    sub_smaller_tinv=vcat(map(i->c_L_tinv[i:i,:],smaller_indices)...)
+    sub_bigger_q=vcat(map(i->c_L_q[i:i,:],bigger_indices)...)
+    M=vcat(sub_smaller_tinv,sub_bigger_q)
+    B=kernel(M, side=:left)
+    rows=vcat(rows,B[1,1:length(smaller_indices)]*sub_smaller_tinv)
+  end
+  c_L_qt = matrix(Q,l,l,rows)
+
+  triv=[partition([n])]
+  for i in 2:r
+    push!(triv,partition([]))
+  end
+  triv=multipartition(triv)
+  index_triv=findfirst(x-> x==triv,mps)
+  d=[1//c_L_qt[i,index_triv] for i in 1:l]
+  diag=diagonal_matrix(Q,d)
+  c_L_qt_H=diag*c_L_qt
+  return c_L_qt_H
+end
+
+@doc raw"""
+    wreath_macdonald_polynomial(lbb::Multipartition,
+                                wperm::PermGroupElem,
+                                coroot::Vector{Int};
+                                parent::MPolyRing{<:QQAbFieldElem}=
+                                polynomial_ring(abelian_closure(QQ)[1], [:q,:t];cached=true)[1])
+
+Given a multipartition `lbb` of size ``n`` and length ``r`` and an element of the affine Weyl group of type ``A^{(1)}_{r-1}``
+(seen as the semi-direct product of the symmetric group on ``r`` letters with the coroot lattice
+of the finite type ``A_{r-1}``), this function returns the coefficients of the wreath Macdonald polynomial
+associated with `lbb` and the affine Weyl group element in the standard Schur basis of multisymmetric functions. Here is an example of how to use it:
+
+```jldoctest
+julia> wreath_macdonald_polynomial(multipartition([[1],[],[]]),cperm(1:3),[0,1,-1])
+[q   q^2   1]
+```
+"""
+function wreath_macdonald_polynomial(lbb::Multipartition,
+                                     wperm::PermGroupElem,
+                                     coroot::Vector{Int};
+                                     parent::MPolyRing{<:QQAbFieldElem}=
+                                     polynomial_ring(abelian_closure(QQ)[1], [:q,:t];cached=true)[1])
+
+  Q = fraction_field(parent)
+  q,t = gens(parent)
+
+  r=length(lbb)
+  n=sum(lbb)
+  mps=collect(multipartitions(n,r))
+  l=length(mps)
+
+  c_L=C_L(r,n,Q,t,mps)
+  c_L_q=map_entries(f->numerator(f)(0,q)//denominator(f)(0,q),c_L)
+  c_L_tinv=map_entries(f->numerator(f)(0,1//t)//denominator(f)(0,1//t),c_L)
+
+  smallers=smaller_ord(lbb, wperm, coroot,mps) #tinv
+  biggers=bigger_ord(lbb, wperm, coroot,mps) #q
+  smaller_indices=sort!(indexin(smallers,mps))
+  bigger_indices=sort!(indexin(biggers,mps))
+  sub_smaller_tinv=vcat(map(i->c_L_tinv[i:i,:],smaller_indices)...)
+  sub_bigger_q=vcat(map(i->c_L_q[i:i,:],bigger_indices)...)
+  M=vcat(sub_smaller_tinv,sub_bigger_q)
+  B=kernel(M, side=:left)
+  row=B[1,1:length(smaller_indices)]*sub_smaller_tinv
+  c_L_qt = matrix(Q,1,l,row)
+
+  triv=[partition([n])]
+  for i in 2:r
+    push!(triv,partition([]))
+  end
+  triv=multipartition(triv)
+  index_triv=findfirst(x-> x==triv,mps)
+  c_L_qt_H=1//c_L_qt[index_triv]*c_L_qt
+  return c_L_qt_H
+end

experimental/WreathMacdonaldpols/test/runtests.jl

@@ -0,0 +1,23 @@
+@testset "wreath macdonald polynomials" begin
+  K,_ = abelian_closure(QQ)
+  parent, = polynomial_ring(K, [:q,:t], cached=false)
+  Q = fraction_field(parent)
+  (q,t) = gens(parent)
+
+  result_p1 = matrix(Q,[t^2 t 1; q t 1; q q^2 1])
+  @test result_p1 == wreath_macdonald_polynomials(1,3,(@perm 3 (1,2,3)),[0,1,-1];parent)
+
+  result_p2 =matrix(Q,[[1//q^4, 1//q, (q^3 + 1)//q^3, 1//q^2, q, (q^3 + 1)//q^2, (q^3 + 1)//q, 1, q^3],
+              [t//q^3, t^3//q^4, (q*t + t^3)//q^3, t//q,  t^3//q^2, (q^2*t + t^2)//q^3, (q^2*t + t^2)//q^2, 1, t^2//q],
+              [1//(q^2*t), 1//(q^2*t), (q + t^2)//(q^2*t), t//q, t//q, (q^2 + t)//(q^2*t), (q*t + 1)//q, 1, 1],
+              [1//t^2, q//t^4, (q + t^2)//t^2, t^2, q, (q + t^2)//t^3, (q + t^2)//t, 1, q//t^2],
+              [1//t^2, t, t^3 + 1, t^2, t^5, (t^3 + 1)//t, t^4 + t, 1, t^3],
+              [1//q, t, q*t + 1, q, q^2*t, (q^2 + t)//q, q^2 + t, 1, q*t],
+              [1//t^2, 1//t^2, (q + t^2)//t^2, q, q, 2//t, (q + t^2)//t, 1, 1],
+              [q//t, q^4//t, (q^5 + q^2)//t, q^3//t, q^6//t, (q^3 + q*t)//t, (q^4 + q^2*t)//t, 1, q^3],
+              [1//t^2, 1//q, (q + t^2)//t^2, q^2//t^2, q, (q + t^2)//(q*t), (q + t^2)//t, 1, t^2//q]])
+  @test result_p2 == wreath_macdonald_polynomials(2,3,(@perm 3 ()),[1,-1,0];parent)
+
+  result_p3 = matrix(Q,[q^2   t   q^3 + q*t   q^4   q^2*t   (q^2 + t)//q   q^2 + t   1   t//q^2])
+  @test result_p3 == wreath_macdonald_polynomial(multipartition([[1,1],[],[]]),(@perm 3 (1,2,3)), [0,1,-1];parent)
+end

src/exports.jl

@@ -1273,7 +1273,7 @@ export number_of_standard_tableaux
 export number_of_transitive_groups, has_number_of_transitive_groups
 export number_of_weak_compositions
 export numerator
-export numerical_lattice 
+export numerical_lattice
 export numerical_lattice_of_K3_cover
 export objective_function
 export omega_group