Great stuff. The incredible identity at the end checks out:<p><pre><code> import sympy
def divisor_power_sum(n, k):
return sum(d**k for d in sympy.divisors(n))
n = 0
while True:
n += 1
lhs = divisor_power_sum(n, 7)
rhs = divisor_power_sum(n, 3) + 120 * sum(
divisor_power_sum(m, 3) * divisor_power_sum(n - m, 3)
for m in range(1, n))
print(n, lhs, rhs)
assert lhs == rhs
</code></pre>
(Yes I see sympy has a `divisor_sigma` already, but I wanted to implement it myself.)<p>(BTW, one of my earliest well-received comments on HN, back from 2015, was about 1/999999999999999999999998999999999999999999999999, discussed in the post: <a href="https://news.ycombinator.com/item?id=9816375" rel="nofollow">https://news.ycombinator.com/item?id=9816375</a>)
Also see generatingfunctionology: <a href="https://www2.math.upenn.edu/~wilf/DownldGF.html" rel="nofollow">https://www2.math.upenn.edu/~wilf/DownldGF.html</a>