## Square Pyramidal Numbers

### March 14, 2017

Cannonballs are traditionally stacked in a pyramid with a square base. A stack with a square base of 15 cannonballs has 15 × 15 = 225 on the bottom level, 14 × 14 = 196 on the level above that, and so on, a total of 1240 cannonballs.

Your task is to write a program to compute the number of cannonballs in a stack; use it to compute the number of cannonballs in a stack with a base of 1,000,000 cannonballs on a side. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

According to oeis.org sequence A000330…

SPN(1000000) = 333333833333500000

This is probably the most minimalistic challenge I’ve seen here so far. I guess it’s more of a challenge in programming efficiency than anything else, using the right types and the a memory-efficient process overall.

Anyway, here is my take on it using Julia: https://app.box.com/s/1ipc05ylpj5f2029r2vratnehdbyttuv

Since the 1mil input was a bit trivial, I decided to go with something higher to test it, since even the most inefficient algorithm could work out the spn of 1 million:

julia> spn(10000000000)

5614543182559713920072066560

This result was obtain in about 6.9 seconds, using 150 memory allocations. The spn(1 million) gives the same result as the formula that Milbrae used previously (333333833333500000), in about 2 milliseconds and 5 memory allocations.

In Haskell…

Let’s generalize a bit and find the sum of any power of the first n integers. https://en.wikipedia.org/wiki/Bernoulli_polynomials#Sums_of_pth_powers tells us we can do this with Bernoulli polynomials, which in turn can be calculated with Stirling numbers and falling factorials. Here’s some Haskell:

Better make that:

else the type inference makes everything an Int rather than an arbitrary precision Integer, so Globules’ big example comes out as 5575454574161208832 rather than 333333333383333333335000000000. I’m surprised Haskell doesn’t give a runtime error on overflow here (maybe there is a runtime option).

@Zack: My blog has readers at all levels. I will always include a mix of simple and not-so-simple exercises.

@Matthew: You might enjoy https://programmingpraxis.com/2011/02/11/sums-of-powers.

@Praxis: Very interesting. Using the Bernoulli polynomials & Stirling numbers rather than Bernoulli numbers (essentially the same calculation, but rearranged a bit) is nice because it just requires integer operations – the divide by k+1 is always exact as the falling factorial is a multiple of k+1 consecutive numbers. Here’s a revised Haskell version with the the Stirling number calculation memoized and some more test cases:

Final version:

In Java:

return LongStream.rangeClosed(1, 1_000_000).map(y->y*y).sum();

I don’t see Forth here much:

The last takes about 16 seconds on my computer, or much faster with the algebraic formula: