That works, but its exponential running time renders it unsuitable for use. Here’s a dynamic programming version that builds the array in such a way that all intermediate calculations are available as they are needed, one row at a time:

Bernoulli’s formula is strikingly better than the brute-force calculation when n is large. For instance, on my computer (sum-powers 100 1000000) takes about 12 seconds, but (bernoulli-formula 100 1000000) makes the same calculation in 150 milliseconds, about two orders of magnitude better. That’s how Jakob Bernoulli was able to calculate S10(1000) in half of a quarter of an hour.

I believe the formula given for s in the exercise is wrong: it should be sum of k = 1 to m + 1 (or 0 to m), not 1 to n. At least that’s the version I found on wikipedia that worked for me after the 1 to n one didn’t.