Jacob Bernoulli’s numbers and sums of powers

Computes the Bernoulli numbers exactly, uses them to evaluate sums of p-th powers in closed form, and repeats Bernoulli’s own boast: the sum of the tenth powers of the first 1000 integers.

from mathematicskit.combinatorics import bernoulli_numbers, sum_of_powers

The first Bernoulli numbers

for m, b in enumerate(bernoulli_numbers(12)):
    print(f"B_{m:<2d} = {b}")

B_0  = 1
B_1  = 1/2
B_2  = 1/6
B_3  = 0
B_4  = -1/30
B_5  = 0
B_6  = 1/42
B_7  = 0
B_8  = -1/30
B_9  = 0
B_10 = 5/66
B_11 = 0
B_12 = -691/2730

Closed-form power sums

for p in range(5):
    formula = sum_of_powers(100, p)
    direct = sum(k**p for k in range(1, 101))
    print(f"sum of k^{p} for k = 1..100: formula {formula}, direct {direct}")

print(f"\n1^10 + ... + 1000^10 = {sum_of_powers(1000, 10):,}")

sum of k^0 for k = 1..100: formula 100, direct 100
sum of k^1 for k = 1..100: formula 5050, direct 5050
sum of k^2 for k = 1..100: formula 338350, direct 338350
sum of k^3 for k = 1..100: formula 25502500, direct 25502500
sum of k^4 for k = 1..100: formula 2050333330, direct 2050333330

1^10 + ... + 1000^10 = 91,409,924,241,424,243,424,241,924,242,500