back to table · edit · history · where entries came from · files · download
5460 bytes, as of the version from 2026-10-28 10:24 (current). Recorded here, not run.
"""Macdonald polynomials -- numberdb.org/T578.
Run with SageMath:
sage -pip install numberdb # once
sage -python generate.py # verify every stored entry (key for draft)
sage -python generate.py --fill # fill the existing draft, never publish it
The integral J normalization uses c_lambda=prod(1-q^arm*t^(leg+1)).
Exact Gram-Schmidt in the power-sum basis computes P, then multiplies by c.
The monomial-to-power-sum conversion uses set-partition inclusion-exclusion.
Only at the end is the alphabet specialized and coerced to ZZ[x1,...,xn,q,t].
There are no numerical arguments q,t and no approximation or precision setting.
See argument-selection.md for the extent and checks.py for independent checks.
"""
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from collections import Counter
from functools import lru_cache
from math import factorial
import os
import sys
QT = PolynomialRing(QQ, ['q', 't'])
K = QT.fraction_field()
q, t = K.gens()
# Complete degree slices selected after measuring expanded entries.
DEGREE_BOUNDS = ((2, 4), (3, 3), (4, 2))
def partitions(d, largest=None):
"""All partitions, descending lexicographically, including () at d=0."""
if d == 0:
yield ()
return
for first in range(min(d, d if largest is None else largest), 0, -1):
for tail in partitions(d-first, first):
yield (first,) + tail
def set_partitions(items):
if not items:
yield ()
return
first, *rest = items
for blocks in set_partitions(rest):
yield ((first,),) + blocks
for k in range(len(blocks)):
yield blocks[:k] + ((first,) + blocks[k],) + blocks[k+1:]
def monomial_in_power(lam):
result = {}
denominator = ZZ(1)
for count in Counter(lam).values():
denominator *= factorial(count)
for blocks in set_partitions(list(range(len(lam)))):
mu = tuple(sorted((sum(lam[i] for i in b) for b in blocks), reverse=True))
coefficient = QQ(1) / denominator
for block in blocks:
coefficient *= (-1)**(len(block)-1)*factorial(len(block)-1)
result[mu] = result.get(mu, K.zero()) + coefficient
return {mu: c for mu, c in result.items() if c}
def inner(left, right):
value = K.zero()
for mu, a in left.items():
if mu not in right:
continue
norm = K.one()
for r, count in Counter(mu).items():
norm *= ZZ(r)**count*factorial(count)
for r in mu:
norm *= (1-q**r)/(1-t**r)
value += a*right[mu]*norm
return value
@lru_cache(None)
def monic_power_basis(d):
basis = {}
for lam in reversed(list(partitions(d))):
v = monomial_in_power(lam)
for previous in basis.values():
factor = inner(v, previous)/inner(previous, previous)
if factor:
for mu, coefficient in previous.items():
v[mu] = v.get(mu, K.zero()) - factor*coefficient
basis[lam] = {mu: c for mu, c in v.items() if c}
return basis
def hook(lam):
c = K.one()
for i, row in enumerate(lam):
for j in range(row):
arm = row-j-1
leg = sum(other > j for other in lam[i+1:])
c *= 1-q**arm*t**(leg+1)
return c
@lru_cache(None)
def polynomial(n, lam):
if len(lam) > n:
raise ValueError('Partition length exceeds alphabet size')
names = ['x%d' % (i+1) for i in range(n)] + ['q','t']
RQ = PolynomialRing(QQ, names)
RZ = PolynomialRing(ZZ, names)
x = RQ.gens()[:n]
d = sum(lam)
power = {r: sum((g**r for g in x), RQ.zero()) for r in range(1,d+1)}
result = RQ.zero()
c = hook(lam)
for mu, coefficient in monic_power_basis(d)[lam].items():
coefficient = QT(c*coefficient) # no q,t denominator is allowed
term = RQ(str(coefficient))
for r in mu:
term *= power[r]
result += term
return RZ(result) # catches any residual rational coefficient
class MacdonaldPolynomials(numberdb.Generator):
table = 'T578'
parameters = ('n', 'lambda')
type = 'Z[]'
rigour = 'exact'
files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md',
'skill-manifest.json', 'T122.json', 'T261.json')
def enumerate(self):
for n, maximum in DEGREE_BOUNDS:
for d in range(2, maximum+1):
for lam in partitions(d):
if len(lam) <= n:
yield {'n': n, 'lambda': ','.join(map(str, lam))}
def value(self, params, digits=None):
lam = tuple(map(int, str(params['lambda']).split(',')))
return polynomial(int(params['n']), lam)
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key = sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].strip().isidentifier():
key = key.split('=',1)[1].strip().strip('\"\'')
numberdb.configure(api_key=key)
generator = MacdonaldPolynomials()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Fill exact integral Macdonald polynomials'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)