How could you guess the formula for the sum of the first n fifth powers?

The following three formulae are reasonably well known:

1 + 2 + 3 + ... + n = n(n+1)/2

1² + 2² + 3² + ... + n ² = n(n+1)(2n+1)/6

1³ + 2³ + 3³ + ... + n ³ = (n(n+1)/2)²

(The sums of first and second powers arise pretty often naturally; the sum of cubes is rare, but it's easy to remember because the sum of the first n cubes is the square of the sum of the first n natural numbers.)

The first member of this series I can't remember is the following:

1⁴ + 2⁴ + 3⁴ + ... + n ⁴ = n(n+1)(2n+1)(3n²+3n-1)/30

and generally, the sum of the first n kth powers is a polynomial of degree k+1.

I ran into these formulas, which I'd seen plenty of times before while perusing the book Gamma: Exploring Euler's Constant by Julian Havil, which I had heard of a few years ago, forgotten about, and then found while browsing the library shelves. (Also in German, under the title GAMMA: Eulers Konstante, Primzahlstrände und die Riemannsche Vermutung.) This is a most interesting book, at least if you're someone like me who pretends to know number theory but really doesn't.

Anyway, back to the main story. Say I wanted to know the sum of the first n fifth powers. Well, there's a general method for finding the formula of the first k powers; it involves the Bernoulli numbers. But let's say I didn't know that. Let's say somebody hands me the sequence

1, 33, 276, 1300, 4425, 12201, 29008, 61776, 120825, 220825

in which the nth term, s_n, is the sum 1⁵ + 2⁵ + ... + n⁵ -- but doesn't tell me that's where the sequence comes from -- and challenges me to guess a formula for it in "closed form". (Smart-asses who will say that there are infinitely many formulas are hereby asked to leave.) How would I guess it?

Well, it can't hurt to find factorizations for these numbers. And if you do that you get

1, 3¹ 11¹, 2² 3¹ 23¹, 2² 5² 13¹, 3¹ 5² 59¹, 3¹ 7² 83¹, 2⁴ 7² 37¹, 2⁴ 3³ 11¹ 13¹, 3³ 5² 179¹, 5² 11² 73¹

and this seems interesting; these numbers seem to have lots of small factors. Furthermore, a fair number of them seem to have one largeish prime factor, which I've bolded. (Yes, I realize, 11 times 13 isn't prime, but I actually did think of it as a large factor.) What are the large factors that I observe in these numbers? They are

?, 11, 23, ?, 59, 83, ?, 143, 179, ?

and the nth of these is easily seen to be 2n² + 2n - 1. (Some terms don't show up from inspection of the factorizations because they get "lost in the noise", as it were.)

From there the rest is pretty easy. We see that often (but not always), s_n is divisible by 2n² + 2n - 1. You can check that for n = 1, 2, ..., 10, the term s_n is always divisible by (2n² + 2n - 1)/3. So now consider the sequence t_n = s_n / ((2n² + 2n - 1)/3). The numbers t₁ through t₁₀ are

1, 9, 36, 100, 225, 441, 784, 1296, 2025, 3025

and I recognized that all of these are squares; in particular t_n = (n(n+1)/2)².

Putting everything together, I get the conjecture that the sum of the first n fifth powers

s_n = n²(n+1)²(2n²+2n-1)/12

which could be proven by induction, but actually writing out the proof is best left to undergrads.

The method here is reminiscent of Enumeration of Matchings: Problems and Progress by James Propp. In that article, Propp lists various unsolved problems in the enumeration of tilings, and conjectures that some of them might have answers which are given by simple product formulas, because actually counting the tilings in question gave numbers with nice prime factorizations.

Edit, 9:13 pm: of course this is not the only method, or even the best method; it's just the method I played around with this morning. See the comments for other methods.