Great stuff. The incredible identity at the end checks out: import sympy def div...

Great stuff. The incredible identity at the end checks out:

    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

(Yes I see sympy has a `divisor_sigma` already, but I wanted to implement it myself.)

(BTW, one of my earliest well-received comments on HN, back from 2015, was about 1/999999999999999999999998999999999999999999999999, discussed in the post: https://news.ycombinator.com/item?id=9816375)