attachment:q_characters.sage of days116

   1 from sage.structure.sage_object import SageObject
   2 from sage.combinat.root_system.cartan_type import CartanType
   3 from sage.data_structures.blas_dict import add as blas_dict_add
   4 from sage.data_structures.blas_dict import axpy
   5 
   6 class Monomial(SageObject):
   7     """
   8     A monomial of a `q`-character.
   9     """
  10     def __init__(self, d):
  11         self._data = dict(d)
  12 
  13     def __hash__(self):
  14         return hash(tuple(sorted(self._data.items())))
  15 
  16     def __eq__(self, other):
  17         return type(self) == type(other) and self._data == other._data
  18 
  19     def __ne__(self, other):
  20         return not (self == other)
  21 
  22     def __mul__(self, other):
  23         if len(self._data) < len(other._data):
  24             self, other = other, self
  25         ret = Monomial(self._data)
  26         for k in other._data:
  27             if k in ret._data:
  28                 ret._data[k] += other._data[k]
  29                 if ret._data[k] == 0:
  30                     del ret._data[k]
  31             else:
  32                 ret._data[k] = other._data[k]
  33         return ret
  34 
  35     def _repr_(self):
  36         if not self._data:
  37             return '1'
  38         def exp_repr(e):
  39             if e == 1:
  40                 return ''
  41             return '^%s' % e
  42         return '*'.join('Y{}[{}]'.format(i, qpow) + exp_repr(self._data[i,qpow])
  43                         for i,qpow in sorted(self._data))
  44 
  45     def _latex_(self):
  46         if not self._data:
  47             return '1'
  48         def exp_repr(e):
  49             if e == 1:
  50                 return ''
  51             return '^{%s}' % e
  52         def q_exp(e):
  53             if e == 0:
  54                 return '1'
  55             if e == 1:
  56                 return 'q'
  57             return 'q^{{{}}}'.format(e)
  58         return ' '.join('Y_{{{},{}}}'.format(i, q_exp(qpow)) + exp_repr(self._data[i,qpow])
  59                         for i,qpow in sorted(self._data))
  60 
  61     def Y_iq_powers(self, i):
  62         """
  63         Return a ``dict`` whose keys are the powers of `q` for the variables
  64         `Y_{i,q^k}` occuring in ``self`` and the values are the exponents.
  65         """
  66         return {k[1]: self._data[k] for k in self._data if k[0] == i}
  67 
  68     def mult_Y_iq(self, i, k, e):
  69         """
  70         Multiply (and mutate) ``self`` by `Y_{i,q^k}^e`.
  71         """
  72         if (i, k) in self._data:
  73             if self._data[i,k] == -e:
  74                 del self._data[i,k]
  75             else:
  76                 self._data[i,k] += e
  77         else:
  78             self._data[i,k] = e
  79 
  80     def is_dominant(self):
  81         return all(k >= 0 for k in self._data.values())
  82 
  83 class qCharacter(SageObject):
  84     """
  85     A `q`-character.
  86     """
  87     def __init__(self, d):
  88         self._poly_dict = dict(d)
  89 
  90     def _repr_(self):
  91         def coeff(c):
  92             if c == 1:
  93                 return ''
  94             return repr(c) + '*'
  95         return " + ".join(coeff(self._poly_dict[m]) + repr(m) if m._data
  96                           else repr(self._poly_dict[m])
  97                           for m in self._poly_dict)
  98 
  99     def _latex_(self):
 100         def coeff(c):
 101             if c == 1:
 102                 return ''
 103             return latex(c)
 104         return " + ".join(coeff(self._poly_dict[m]) + latex(m) if m._data
 105                           else latex(self._poly_dict[m])
 106                           for m in self._poly_dict)
 107 
 108     def __eq__(self, other):
 109         return type(self) == type(other) and self._poly_dict == other._poly_dict
 110 
 111     def __ne__(self, other):
 112         return not (self == other)
 113 
 114     def __len__(self):
 115         """
 116         Return the number of monomials of ``self``.
 117         """
 118         return sum(self._poly_dict.values())
 119 
 120     def __getitem__(self, m):
 121         if m not in self._poly_dict:
 122             return 0
 123         else:
 124             return self._poly_dict[m]
 125 
 126     def __add__(self, other):
 127         return qCharacter(blas_dict_add(self._poly_dict, other._poly_dict))
 128 
 129     def __sub__(self, other):
 130         return qCharacter(axpy(-1, other._poly_dict, self._poly_dict))
 131 
 132     def __mul__(self, other):
 133         ret = qCharacter({})
 134         for m in self._poly_dict:
 135             for mp in other._poly_dict:
 136                 mpp = m * mp
 137                 if mpp in ret._poly_dict:
 138                     ret._poly_dict[mpp] += self._poly_dict[m] * other._poly_dict[mp]
 139                 else:
 140                     ret._poly_dict[mpp] = self._poly_dict[m] * other._poly_dict[mp]
 141         return ret
 142 
 143     @cached_method
 144     def dominant_monomials(self):
 145         return tuple([m for m in self._poly_dict if m.is_dominant()])
 146 
 147 def FM_algorithm(ct, initial):
 148     """
 149     Return the `q`-character of ``m`` of the Cartan type ``ct``
 150     using the FM algorithm.
 151 
 152     INPUT:
 153 
 154     - ``ct`` -- a Cartan type
 155     - ``initial`` -- the initial monomial data as a ``dict`` whose keys
 156       are pairs `(i, k)` corresponding to `Y_{i,q^k}` and the value is
 157       the corresponding expontent
 158 
 159     EXAMPLES:
 160 
 161     We create the `q`-character whose dominant monomial
 162     is `Y_{1,q^2} Y_{2,q^{-1}}` in type `A_2`::
 163 
 164         sage: FM_algorithm(['A',2], {(1,2):1, (2,-1):1})
 165         Y1[0]*Y2[1]^-1*Y2[5]^-1 + Y2[-1]*Y2[5]^-1 + Y1[4]^-1*Y2[-1]*Y2[3]
 166          + Y1[2]*Y2[-1] + Y1[2]^-1*Y1[4]^-1*Y2[3] + Y1[0]*Y1[4]^-1*Y2[1]^-1*Y2[3]
 167          + Y1[0]*Y1[2]*Y2[1]^-1 + Y1[2]^-1*Y2[5]^-1
 168 
 169     Next, we compute the `q`-character of `Y_{2,q^{-1}}` in type `C_2`::
 170 
 171         sage: FM_algorithm(['C',2], {(2,-1):1})
 172         Y1[0]*Y1[4]^-1 + Y2[-1] + Y2[5]^-1 + Y1[0]*Y1[2]*Y2[3]^-1
 173          + Y1[2]^-1*Y1[4]^-1*Y2[1]
 174 
 175     Now `Y_{2,q^{-1}} Y_{2,q^1}`::
 176 
 177         sage: FM_algorithm(['C',2], {(2,-1):1, (2,1):1})
 178         Y1[2]*Y1[6]^-1*Y2[-1] + Y1[0]*Y1[4]^-1*Y2[1] + Y1[0]*Y1[2]*Y2[3]^-1*Y2[7]^-1
 179          + Y2[-1]*Y2[1] + Y1[2]*Y1[6]^-1*Y2[5]^-1 + Y1[0]*Y1[2]^2*Y1[6]^-1*Y2[3]^-1
 180          + Y1[2]*Y1[4]*Y2[-1]*Y2[5]^-1 + 2*Y1[0]*Y1[2]*Y1[4]^-1*Y1[6]^-1
 181          + Y1[2]^-1*Y1[4]^-2*Y1[6]^-1*Y2[1]*Y2[3] + Y1[4]^-1*Y1[6]^-1*Y2[3]*Y2[5]^-1
 182          + Y1[2]^-1*Y1[4]^-1*Y2[1]*Y2[7]^-1 + Y1[0]*Y1[2]^2*Y1[4]*Y2[3]^-1*Y2[5]^-1
 183          + Y1[0]*Y1[4]^-1*Y2[7]^-1 + Y1[0]*Y1[2]*Y2[5]^-1 + Y1[2]^-1*Y1[4]^-1*Y2[1]^2
 184          + Y1[4]^-1*Y1[6]^-1*Y2[1] + Y1[0]*Y1[2]*Y2[1]*Y2[3]^-1
 185          + Y1[4]^-1*Y1[6]^-1*Y2[-1]*Y2[3] + Y2[-1]*Y2[7]^-1 + Y2[5]^-1*Y2[7]^-1
 186          + 2*Y2[1]*Y2[5]^-1 + Y1[2]*Y1[4]*Y2[5]^-2 + Y1[0]*Y1[4]^-2*Y1[6]^-1*Y2[3]
 187 
 188     Finally, `Y_{2,q^{-1}} Y_{2,q^3}`, which gives the KR module `W^{2,2}`::
 189 
 190         sage: FM_algorithm(['C',2], {(2,-1):1, (2,3):1})
 191         Y1[2]^-1*Y1[4]^-1*Y2[1]*Y2[9]^-1 + Y1[4]*Y1[8]^-1*Y2[-1]
 192          + Y1[2]^-1*Y1[4]^-1*Y1[6]^-1*Y1[8]^-1*Y2[1]*Y2[5]
 193          + Y1[0]*Y1[2]*Y1[4]*Y1[8]^-1*Y2[3]^-1 + Y2[5]^-1*Y2[9]^-1
 194          + Y1[0]*Y1[4]^-1*Y1[6]^-1*Y1[8]^-1*Y2[5]
 195          + Y1[0]*Y1[2]*Y1[4]*Y1[6]*Y2[3]^-1*Y2[7]^-1 + Y1[0]*Y1[2]*Y2[3]^-1*Y2[9]^-1
 196          + Y2[-1]*Y2[9]^-1 + Y2[-1]*Y2[3] + Y1[0]*Y1[4]^-1*Y2[9]^-1
 197          + Y1[6]^-1*Y1[8]^-1*Y2[-1]*Y2[5] + Y1[4]*Y1[6]*Y2[-1]*Y2[7]^-1
 198          + Y1[0]*Y1[2]*Y1[6]^-1*Y1[8]^-1*Y2[3]^-1*Y2[5]
 199     """
 200     from itertools import product
 201     ct = CartanType(ct)
 202     I = ct.index_set()
 203     d = ct.symmetrizer()
 204     m = Monomial(initial)
 205     ret = {m: 1}
 206     coloring = {m: {i: 0 for i in I}}
 207     CM = ct.cartan_matrix()
 208     adjacent = {i: {j: CM[ij,ii] for ij,j in enumerate(I) if i != j and CM[ij,ii]}
 209                 for ii,i in enumerate(I)}
 210     # We go through weights (and hence monomials) by depth
 211     next = [[m]]
 212     num = -1
 213     while next:
 214         cur = next.pop(0)
 215         num += 1
 216         for m in cur:
 217             for i in I:
 218                 # There is nothing to check/do for i and m such that
 219                 #   the i-color is equal to the coefficient.
 220                 if coloring[m][i] == ret[m]:
 221                     #print("coloring equals coefficient {} and {}".format(m, i))
 222                     continue
 223 
 224                 # Get the powers k of q for the variables Y_{i,q^k}
 225                 powers = m.Y_iq_powers(i)
 226 
 227                 # Check for failure of FM algorithm by verifying m is i-dominant
 228                 assert all(k >= 0 for k in powers.values()), "FM algorithm failed at %s" % m
 229 
 230                 ##### Perform the i-expansion #####
 231                 #print("Performing {}-expansion of {}".format(i,m))
 232 
 233                 # Compute the q-strings in general position to determine
 234                 #   which sl_2 representations m corresponds to.
 235                 di = d[i]
 236                 adj = adjacent[i]
 237                 top_q = []
 238                 q_str_len = []
 239                 while powers:
 240                     cur_pow = min(powers)
 241                     cur_len = 0
 242                     while cur_pow in powers:
 243                         if powers[cur_pow] == 1:
 244                             del powers[cur_pow]
 245                         else:
 246                             powers[cur_pow] -= 1
 247                         cur_len += 1
 248                         cur_pow += 2 * di
 249                     top_q.append(cur_pow - 2*di)
 250                     q_str_len.append(cur_len)
 251                 #print("q-strings: {}, {}".format(top_q, q_str_len))
 252 
 253                 # Build the i-string and add in the corresponding variables
 254                 t = (ret[m] - coloring[m][i])
 255                 for exponent in product(*[range(ell+1) for ell in q_str_len]):
 256                     depth = sum(exponent)
 257                     if depth == 0: # Nothing to do for this monomial
 258                         continue
 259 
 260                     while depth > len(next):
 261                         next.append([])
 262 
 263                     # Compute the new monomial
 264                     mp = Monomial(m._data)
 265                     for ind, e in enumerate(exponent):
 266                         # Multiply by prod_{k=0}^e A_{i, a q_i^{1-2k}}
 267                         for k in range(e):
 268                             a = top_q[ind] + di*(1 - 2*k)
 269                             mp.mult_Y_iq(i, a-di, -1)
 270                             mp.mult_Y_iq(i, a+di, -1)
 271                             for j in adj:
 272                                 if adj[j] == -1:
 273                                     mp.mult_Y_iq(j, a, 1)
 274                                 elif adj[j] == -2:
 275                                     mp.mult_Y_iq(j, a-1, 1)
 276                                     mp.mult_Y_iq(j, a+1, 1)
 277                                 elif adj[j] == -3:
 278                                     mp.mult_Y_iq(j, a-2, 1)
 279                                     mp.mult_Y_iq(j, a, 1)
 280                                     mp.mult_Y_iq(j, a+2, 1)
 281 
 282                     # Now add it to the colored polynomial
 283                     if mp not in ret:
 284                         coloring[mp] = {i: 0 for i in I}
 285                         coloring[mp][i] = t
 286                         ret[mp] = t
 287                         next[depth-1].append(mp)
 288                     else:
 289                         coloring[mp][i] += t
 290                         if coloring[mp][i] > ret[mp]:
 291                             ret[mp] = coloring[mp][i]
 292 
 293                 #print(qCharacter(ret))
 294 
 295     return qCharacter(ret)
 296 
 297 def product_dominant_monomials(V1, V2):
 298     ret = []
 299     for m in V1._poly_dict:
 300         for mp in V2._poly_dict:
 301             mpp = m * mp
 302             if mpp.is_dominant():
 303                 ret.extend([mpp]*(V1._poly_dict[m]*V2._poly_dict[mp]))
 304     return ret
 305 
 306 def is_product_simple(V1, V2):
 307     count = 0
 308     for m in V1._poly_dict:
 309         for mp in V2._poly_dict:
 310             mpp = m * mp
 311             if mpp.is_dominant():
 312                 count += 1
 313                 if count > 1:
 314                     return False
 315     return True