import numberdb.sage as numberdb # initialize Sage before named imports
from math import factorial
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.rational_field import QQ
def B(n):
a = [QQ(0)] * (n + 1)
for m in range(n + 1):
a[m] = QQ(1) / QQ(m + 1)
for j in range(m, 0, -1):
a[j - 1] = QQ(j) * (a[j - 1] - a[j])
return a[0]
def mul(f, g, n):
h = [QQ(0)] * (n + 1)
for i, a in enumerate(f):
for j, b in enumerate(g):
if i + j <= n:
h[i + j] += a * b
return h
def q_coeffs(n):
return [
(QQ(2) ** (1 - 2 * k) - QQ(1)) * B(2 * k) / QQ(factorial(2 * k))
for k in range(n + 1)
]
def log_q(n):
q = q_coeffs(n)
q[0] -= QQ(1)
out = [QQ(0)] * (n + 1)
power = [QQ(1)] + [QQ(0)] * n
for m in range(1, n + 1):
power = mul(power, q, n)
sign = QQ(1) if m % 2 else QQ(-1)
for k in range(1, n + 1):
out[k] += sign * power[k] / QQ(m)
return out
def exps(e):
try:
return tuple(e)
except TypeError:
return (e,)
def weight(e):
return sum((i + 1) * a for i, a in enumerate(exps(e)))
def monomial(R, xs, e):
term = R(1)
for x, a in zip(xs, exps(e)):
term *= x ** a
return term
def terms(f, n, exact):
R = f.parent()
xs = R.gens()
out = R(0)
for e, c in f.dict().items():
w = weight(e)
if (exact and w == n) or (not exact and w <= n):
out += c * monomial(R, xs, e)
return out
def Ahat(n):
R = PolynomialRing(QQ, ["p%s" % i for i in range(1, n + 1)])
p = R.gens()
s = {}
for m in range(1, n + 1):
total = sum(((-1) ** (i + 1) * p[i - 1] * s[m - i]
for i in range(1, m)), R(0))
s[m] = total + (-1) ** (m + 1) * QQ(m) * p[m - 1]
log = log_q(n)
exponent = sum((log[m] * s[m] for m in range(1, n + 1)), R(0))
total = term = R(1)
for k in range(1, n + 1):
term = terms(term * exponent / QQ(k), n, False)
total = terms(total + term, n, False)
return terms(total, n, True)
print(Ahat(6))
Every stored component was checked against the Pontryagin-root definition in Formula (1), evaluated in seven Pontryagin roots. The components $\hat A_1$ to $\hat A_4$ were also checked against the terms printed in [1].