|
(Sage)
Pol.<x> = InfinitePolynomialRing(QQ)
@cached_function
def Z(n):
if n == 0: return Pol.one()
return sum(x[k]*Z(n-k) for k in (1..n))/n
@cached_function
def monprod(M):
p = Pol.one()
V = [m.variables() for m in M]
T = cartesian_product(V)
for t in T:
r = [Pol.varname_key(str(u))[1] for u in t]
j = [Pol(M[u[0]]).degree(u[1]) for u in enumerate(t)]
lcm_r = lcm(r)
p *= x[lcm_r]^(prod(r)/lcm_r*prod(j))
return p
@cached_function
def pol_isotop(n, k):
P = Z(n)
p = Pol.zero()
coeffs = P.coefficients()
mons = P.monomials()
C = cartesian_product(k*[mons])
Csorted = [tuple(sorted(u)) for u in C]
Cset = set(Csorted)
for c in Cset:
p += Csorted.count(c)*prod([coeffs[mons.index(u)] for u in c])*monprod(c)
return p
@cached_function
def rule_sub(r, m):
D = 0
for d in divisors(r):
try: D += d*m.degrees()[-d-1]
except: break
return D
def a(n, k=3):
P = Z(n)
coeffs = P.coefficients()
Q = pol_isotop(n, k)
inds = [Pol.varname_key(str(u))[1] for u in Q.variables()]
p = 0
for mon in enumerate(P.monomials()):
m = Pol(mon[1])
p += coeffs[mon[0]]*Q.subs({x[i]:rule_sub(i, m) for i in inds})
return p
| 