denominator-series-3.mag

// This program computes a series expansion for the affine Macdonald denominator
// for U_q(sl_3 hat) as defined in the paper:
//
// P. Etingof and A. Kirillov Jr. On the affine analogue of Jack and Macdonald
// polynomials. Duke Math. J., 78(2):229-256, 1995.
//
// It provides computational evidence for the affine Macdonald denominator
// conjecture given in the paper:
//
// E. M. Rains, Y. Sun, and A. Varchenko.  Affine Macdonald conjectures and
// special values of Felder-Varchenko functions.  Preprint, 2016.
//
// The expansion is done to first order in q^{-2\omega}.  The computation is
// done by computing the first order expansion of a singular vector in
// 
// L_{k - 1, 0, -(k - 1), 3(k - 1)} \otimes L_{3(k - 1), 0, 0}(z),
//
// where L_{k - 1, 0, -(k - 1), 3(k - 1)} is an integrable irreducible module for
// U_q(sl_3 hat) and L_{3(k - 1), 0, 0}(z) is an evaluation module for
// U_q(sl_3 hat).

SetEchoInput(true);
ZZ := Integers();
F0 := Rationals();
R := RootDatum("A2");
U := QuantizedUEA(R);
Qq<q> := CoefficientRing(U);
M<f0, f1, f2, p> := FreeAlgebra(Qq, 4);
A<qlam1, qlam2, qomg> := FunctionField(Qq, 3);

qbinom := function(n)
       nn := ZZ!n;
       return (q^nn - 1/q^nn)/(q - 1/q);
end function;

dot := function(tup1, tup2)
    if #tup1 eq 4 then
        return tup1[1] * tup2[1] + tup1[2] * tup2[2] + tup1[3] * tup2[3] + tup1[4] * tup2[4];
    elif #tup1 eq 3 then
        return tup1[1] * tup2[1] + tup1[2] * tup2[2] + tup1[3] * tup2[3];
    end if;
end function;

sub := function(tup1, tup2)
    if #tup1 eq 4 then
        return <tup1[1] - tup2[1], tup1[2] - tup2[2], tup1[3] - tup2[3], tup1[4] - tup2[4]>;
    elif #tup1 eq 3 then
        return <tup1[1] - tup2[1], tup1[2] - tup2[2], tup1[3] - tup2[3]>;
    end if;
end function;

get_e_weight_gen := function(idx)
	       if idx eq 0 then
	       	  return <1, -1, 0, 1>;
	       elif idx eq 1 then
	       	  return <0, 1, -1, 0>;
	       elif idx eq 2 then
	       	  return <0, 0, 1, -1>;
               end if;
end function;

get_weight_vect := function(mono, pweight, plevel)
       ans := <0, 0, 0, 0>;
       for let_idx in [1..Length(mono)] do
       	   if mono[let_idx] eq f0 then
	        ans := ans + <-1, 1, 0, -1>;
	   elif mono[let_idx] eq f1 then
	   	ans := ans + <0, -1, 1, 0>;
	   elif mono[let_idx] eq f2 then
	   	ans := ans + <0, 0, -1, 1>;
	   else
	   end if;
       end for;
       return ans + <plevel, pweight[1], pweight[2], pweight[3]>;
end function;

get_root_weight := function(idx)
	if idx eq 0 then
	   return <1, -1, 0, 1>;
	elif idx eq 1 then
	   return <0, 1, -1, 0>;
	elif idx eq 2 then
	   return <0, 0, 1, -1>;
	else
	end if;
end function;

get_fund_gen_weight := function(idx)
	if idx eq 0 then
	   return <1, 0, 0, 0>;
	elif idx eq 1 then
	   return <0, 1, 0, 0>;
	elif idx eq 2 then
	   return <0, 1, 1, 0>;
	else
	end if;
end function;

get_weight := function(mono, idx, pweight, plevel)
	   return dot(get_root_weight(idx), get_weight_vect(mono, pweight, plevel));
end function;

get_fund_weight := function(mono, idx, pweight, plevel)
	   return dot(get_fund_gen_weight(idx), get_weight_vect(mono, pweight, plevel));
end function;

get_loop_weight := function(mono)
		ans := 0;
		for let_idx in [1..Length(mono)] do
		    if mono[let_idx] eq f0 then
		       ans := ans - 1;
		    end if;
		end for;
		return ans;
end function;

deloop_weight := function(wt)
	      return <wt[2], wt[3], wt[4]>;
end function;

get_ecomm := function(let, idx, weight)
	  if (idx eq 0 and not(let eq f0)) or (idx eq 1 and not(let eq f1)) or (idx eq 2 and not(let eq f2)) then
	     return 0;
	  else
	     return (q^weight - 1/q^weight) / (q - 1/q);
	  end if;
end function;

get_emult_mono := function(mono, idx, pweight, plevel)
	  ans := 0;
	  for let_idx in [1..Length(mono) - 1] do
	      prev := 1; 
	      for idx2 in [1..let_idx - 1] do
	      	  prev := prev * mono[idx2];
	      end for;
	      after := 1;
	      for idx3 in [let_idx + 1..Length(mono)] do
	      	  after := after * mono[idx3];
	      end for;
	      weight := get_weight(after, idx, pweight, plevel);
	      comm := get_ecomm(mono[let_idx], idx, weight);
	      ans := ans + prev * comm * after;
	  end for;
	  return ans;
end function;

get_emult := function(term, idx, pweight, plevel)
	  ans := 0;
	  for mono in Monomials(term) do
	      ans := ans + MonomialCoefficient(term, mono) * get_emult_mono(mono, idx, pweight, plevel);
	  end for;
	  return ans;
end function;

term_to_vector := function(term, monos)
	  ans := [ MonomialCoefficient(term, mono) : mono in monos ];
	  return ans;
end function;

// return true if weight could appear in the second tensor factor
is_valid_sym_weight := function(wt, kk)
		    for i in [1..3] do
		     	 if wt[i] lt - (kk - 1) then
			    return false;
			 end if;
		     end for;
		     return true;
end function;

// returns true if mono could appear as the first tensor factor in a term of the singular vector
is_valid_sing_weight := function(mono, pweight, plevel, kk)
		     needed_wt := sub(pweight, deloop_weight(get_weight_vect(mono, pweight, plevel)));
		     return is_valid_sym_weight(needed_wt, kk);
end function;

// convert weights from basis of fundamental weights to lying in C^3
convert_wt := function(magma_wt)
	   return <magma_wt[1] * 2/3 + magma_wt[2] * 1/3, - magma_wt[1] * 1/3 + magma_wt[2] * 1/3, - magma_wt[1] * 1/3 - magma_wt[2] * 2/3>;
end function;

// convert weights from basis of simple roots to lying in C^3
convert_rt := function(magma_wt)
	   return <magma_wt[1], - magma_wt[1] + magma_wt[2], - magma_wt[2]>;
end function;

// returns the index of the generator of U corresponding to e_idx if idx > 0 or e_theta if idx = 0
get_mag_e_idx := function(idx, U, R)
	  if idx eq 1 then
	     wt := <1, -1, 0>;
	  elif idx eq 2 then
	     wt := <0, 1, -1>;
	  else
	     wt := <1, 0, -1>;
	  end if;

	  prp := PositiveRootsPerm(U);
	  for idx in [1..#(prp)] do 
	      if convert_rt(PositiveRoots(R)[prp[idx]]) eq wt then
	      	 e_idx := idx;
	      end if;
	  end for;
	  if not(idx eq 0) then
                return 5 + e_idx;
          else
		return e_idx;
          end if;
end function;

// returns the index of the generator of U corresponding to f_idx if idx > 0 or f_theta if idx = 0
get_mag_f_idx := function(idx, U, R)
	  if idx eq 1 then
	     wt := <1, -1, 0>;
	  elif idx eq 2 then
	     wt := <0, 1, -1>;
	  else
	     wt := <1, 0, -1>;
	  end if;

	  prp := PositiveRootsPerm(U);
	  for idx in [1..#(prp)] do 
	      if convert_rt(PositiveRoots(R)[prp[idx]]) eq wt then
	      	 f_idx := idx;
	      end if;
	  end for;
	  if not(idx eq 0) then
                return f_idx;
          else
		return 5 + f_idx;
          end if;
end function;

// rep: U-module
// curr_weight: tuple for current weight vector (all weight spaces are dimension 1)
// e_idx: index of e_i
get_e_trans_coeff := function(rep, curr_weight, e_idx)
		  weights, vects := WeightsAndVectors(rep);
		  wt_idx := -1;
		  for idx in [1..#(weights)] do
		      if convert_wt(weights[idx]) eq curr_weight then
		      	 wt_idx := idx;
		      end if;
		  end for;

		  if wt_idx eq -1 then
		     return 0;
		  end if;

		  if not(e_idx eq 0) then
		  	mag_e_idx := get_mag_e_idx(e_idx, U, R);
			u_elt := U.mag_e_idx;
                  else
			u_elt := (U.3 * U.1 - U.1 * U.3 * q) * q^(-dot(<1, 0, 1>, curr_weight));
		  end if;
		  res := ElementToSequence(u_elt^rep.wt_idx);

		  for idx in [1..#(res)] do
		      if not(res[idx] eq 0) then
		      	 return res[idx];
		      end if;
		  end for;
		  return 0;
end function;

// rep: U-module
// curr_weight: tuple for current weight vector (all weight spaces are dimension 1)
// f_idx: index of f_i
get_f_trans_coeff := function(rep, curr_weight, f_idx)
		  weights, vects := WeightsAndVectors(rep);
		  wt_idx := -1;
		  for idx in [1..#(weights)] do
		      if convert_wt(weights[idx]) eq curr_weight then
		      	 wt_idx := idx;
		      end if;
		  end for;

		  if wt_idx eq -1 then
		     return 0;
		  end if;
		  
		  if not(f_idx eq 0) then
		        mag_f_idx := get_mag_f_idx(f_idx, U, R);
			u_elt := U.mag_f_idx;
	          else
			u_elt := (U.6 * U.8 - U.8 * U.6 / q) * q^(dot(<1, 0, 1>, curr_weight));
		  end if;
		  res := ElementToSequence(u_elt^rep.wt_idx);
		  for idx in [1..#(res)] do
		      if not(res[idx] eq 0) then
		      	 return res[idx];
		      end if;
		  end for;
		  return 0;
end function;

// Apply gen to vect, where vect is interpreted as an element of the tensor product by tensoring each monomial by the appropriate basis vector to make the total weight curr_weight
get_tens_mult := function(vect, gen, rep, pweight, plevel, curr_weight, gb)
	      ans := M!0;
	      if gen eq f0 then
	      	 gen_idx := 0;
	      elif gen eq f1 then
	         gen_idx := 1;
	      else
	         gen_idx := 2;
              end if;

	      for mono in Monomials(vect) do
	      	  ans := ans + MonomialCoefficient(vect, mono) * get_f_trans_coeff(rep, deloop_weight(sub(curr_weight, get_weight_vect(mono, pweight, plevel))), gen_idx) * mono;

		  weight_exp := dot(get_root_weight(gen_idx),
		                    sub(curr_weight, get_weight_vect(mono, pweight, plevel)));
		  ans := ans + MonomialCoefficient(vect, mono) * q^(- weight_exp) * NormalForm(gen * mono, gb);
	      end for;
	      return ans;
end function;

// sing_vect: the singular vector, as an algebra element
// mono: the monomial to get the value of the intertwiner on
// rep: module object
// pweight: finite-type h.w.
// plevel: level
// gb: Groebner basis
get_inter_val := function(sing_vect, mono, rep, pweight, plevel, gb)
	      curr_weight := <plevel, pweight[1], pweight[2], pweight[3]>;
	      ans := sing_vect;
	      for idx in [1..Length(mono)] do
	          idx_mod := Length(mono) - idx + 1;
	      	  if not(mono[idx_mod] eq p) then
		     ans := get_tens_mult(ans, mono[idx_mod], rep, pweight, plevel, curr_weight, gb);
		     if mono[idx_mod] eq f0 then
		     	curr_weight := curr_weight + <-1, 1, 0, -1>;
	 	     elif mono[idx_mod] eq f1 then
		        curr_weight := curr_weight + <0, -1, 1, 0>;
	             elif mono[idx_mod] eq f2 then
                        curr_weight := curr_weight + <0, 0, -1, 1>;
		     end if;	
		  end if;
	      end for;
	      return ans;
end function;

get_tnorm := function(kk)
	  ans := qlam1^(kk - 1) * qlam2^(kk - 1);
	  for i in [1..kk - 1] do
	  ans := ans * (1 - q^(2 * i) / qlam1) * (1 - q^(2 * i) / qlam2) * (1 - q^(2 * i) / qlam1 / qlam2);
	  end for;
	  return ans;
end function;

// Return the trace associated to some kk.
get_trace := function(kk)
	  pweight := <(kk - 1), 0, - (kk - 1)>;
	  plevel := 3 * (kk - 1);
	  SV := HighestWeightRepresentation(U, [3 * (kk - 1), 0]);
	  SM := HighestWeightModule(U, [3 * (kk - 1), 0]);

	  // Find Groebner basis for first few degrees of irreducible integrable
	  // module using Serre relations.
	  max_deg := 30 + 10 * kk;
	  urels := [f0^2 * f1 - GaussianBinomial(2, 1, q)*f0*f1*f0 + f1*f0^2,
               	    f0^2 * f2 - GaussianBinomial(2, 1, q)*f0*f2*f0 + f2*f0^2,
		    f1^2 * f0 - GaussianBinomial(2, 1, q)*f1*f0*f1 + f0*f1^2,
		    f1^2 * f2 - GaussianBinomial(2, 1, q)*f1*f2*f1 + f2*f1^2,
		    f2^2 * f0 - GaussianBinomial(2, 1, q)*f2*f0*f2 + f0*f2^2,
		    f2^2 * f1 - GaussianBinomial(2, 1, q)*f2*f1*f2 + f1*f2^2];
	  irels := [f1^(dot(pweight, <1, -1, 0>) + 1),
	            f2^(dot(pweight, <0, 1, -1>) + 1)];
          trels := [f0^2,
          	    f0 * f1 * f0, f0 * f2 * f0,
          	    f0 * f1^2 * f0, f0 * f1 * f2 * f0, f0 * f2 * f1 * f0,
		    f0 * f2^2 * f0, f0 * f1^3 * f0, f0 * f1^2 * f2 * f0,
		    f0 * f2 * f1^2 * f0, f0 * f2^3 * f0, f0 * f1^4 * f0,
		    f0 * f1^3 * f2 * f0, f0 * f1^2 * f2 * f1 * f0,
                    f0 * f1 * f2 * f1^2 * f0, f0 * f1^2 * f2^2 * f0,
		    f0 * f2^2 * f1^2 * f0, f0 * f1 * f2 * f1 * f2 * f0,
		    f0 * f2 * f1 * f2 * f1 * f0, f0 * f1 * f2^3 * f0,
		    f0 * f2 * f1 * f2^2 * f0, f0 * f2^2 * f1 * f2 * f0,
		    f0 * f2^3 * f1 * f0, f0 * f2^4 * f0];
          I_gens := [p^2, p * f0, p * f1, p * f2] cat urels cat [i * p : i in irels] cat trels;
	  I := ideal<M | I_gens>;
	  gb := GroebnerBasis(I_gens, max_deg);
	  lgb := [LeadingMonomial(p) : p in gb];

	  // Find basis of monomials for first few degrees of irreducible 
	  // integrable module.
	  monomials := [[p]];
	  vars := [f2, f1, f0, p];
	  for i in [2..max_deg] do
	      monomials[i] := [r * s : r in vars, s in monomials[i - 1] | NormalForm(r * s, lgb) ne 0];
	  end for;

	  // Compute weights of each monomial.
	  mono_list := [];
	  weights := [];
	  loop_weights := [];
	  mono_idx := 1;
	  for deg_mono in monomials do
	      for mono in deg_mono do
	          mono_list[mono_idx] := mono;
		  weights[mono_idx] := get_weight_vect(mono, pweight, plevel);
		  loop_weights[mono_idx] := get_loop_weight(mono);
		  mono_idx := mono_idx + 1;
	      end for;
	  end for;
	  num_mono := mono_idx - 1;

	  // Compute actions of generators e_0, e_1, and e_2.
	  e_mat := [ ];
	  for i in [0..2] do
	     e_mat[i + 1] := [];
	     for idx in [1..num_mono] do
	       	mult_res := get_emult(mono_list[idx], i, pweight, plevel);
		if not(mult_res eq 0) then
		   norm_mult_res := NormalForm(mult_res, gb);
		   e_mat[i + 1][idx] := [];
		   e_mat[i + 1][idx] := term_to_vector(norm_mult_res, mono_list);
		else
		   e_mat[i + 1][idx] := [];
		   for idx2 in [1..num_mono] do
		       e_mat[i + 1][idx][idx2] := Qq!0;
		   end for;
		end if;
	    end for;
	 end for;

	 // Compute actions of generators e_0, e_1, and e_2 in tensor product.
	 // The coproduct is \Delta(e_i) = e_i \ot 1 + q^{e_i} \ot e_i.
	 e_tensmat := [];
	 for i in [0..2] do
	     e_tensmat[i + 1] := [];
	     for idx in [1..num_mono] do
	       	e_tensmat[i + 1][idx] := [];
	    	for j in [1..num_mono] do 
    		   e_tensmat[i + 1][idx][j] := e_mat[i + 1][idx][j];
		end for;
		shifted_wt := deloop_weight(get_e_weight_gen(i)) + sub(pweight, deloop_weight(get_weight_vect(mono_list[idx], pweight, plevel)));
		if is_valid_sym_weight(shifted_wt, kk) and is_valid_sing_weight(mono_list[idx], pweight, plevel, kk) then
	    	   e_tensmat[i + 1][idx][idx] := q^(get_weight(mono_list[idx], i, pweight, plevel)) * get_e_trans_coeff(SM, sub(pweight, deloop_weight(get_weight_vect(mono_list[idx], pweight, plevel))), i);
		end if;
	    end for;
	end for;

	// Filter rows of the matrix actions of e_0, e_1, e_2 to include only
	// those landing in the degree 0 piece.
	y_monos := [];
	y_seq := [];
	for idx in [1..3] do
	    y_seq[idx] := [];
	    i_count := 1;
	    for i in [1..num_mono] do 
	       	if is_valid_sing_weight(mono_list[i], pweight, plevel, kk) then
		   y_seq[idx][i_count] := [];
		   y_monos[i_count] := mono_list[i];
		   j_count := 1;
		   for j in [1..num_mono] do
		       if is_valid_sing_weight(mono_list[j], pweight, plevel, kk) then
		       	  if (loop_weights[j] eq 0) or not(idx - 1 eq 0) then
			     y_seq[idx][i_count][j_count] := e_tensmat[idx][i][j];
			  else
			     y_seq[idx][i_count][j_count] := Qq!0;
			  end if;
			  j_count := j_count + 1;
		      end if;
		   end for;
		   i_count := i_count + 1;
	        end if;
	    end for;
	 end for;

	 // Compute the singular vector up to degree 1.
	 y0 := Matrix(y_seq[1]);
	 y1 := Matrix(y_seq[2]);
	 y2 := Matrix(y_seq[3]);
	 yy := HorizontalJoin(HorizontalJoin(y0, y1), y2);
	 yy_trig := HorizontalJoin(y1, y2);
	 kernel := Kernel(yy);
	 ker_trig := Kernel(HorizontalJoin(y1, y2));
	 sv_temp := Basis(kernel)[1];
	 sing_vect := M!0;
	 for i in [1.. (# y_monos)] do
	     sing_vect := sing_vect + sv_temp[i] * y_monos[i];
	 end for;

	 // Compute the intertwiner up to first order
	 inter_mat := [];
	 for i in [1..num_mono] do
	     inter_mat[i] := term_to_vector(get_inter_val(sing_vect, mono_list[i], SM, pweight, plevel, gb), mono_list);
	 end for;     

	 // Compute trace up to first order
	 trace := A!0;
	 trace_parts := [];
	 for i in [1..num_mono] do
	     trace_parts[i] := inter_mat[i][i] * qlam1^(get_fund_weight(mono_list[i], 1, pweight, plevel)) * qlam2^(get_fund_weight(mono_list[i], 2, pweight, plevel)) * qomg^(-loop_weights[i]);
	     trace := trace + trace_parts[i];
	 end for;
	 return trace, trace_parts;
end function;

for kk in [2..4] do
    time print kk, get_trace(kk) / get_tnorm(kk); 
end for;