## Square Pyramidal Numbers

### March 14, 2017

This is easy; we want the sum of the squares of the first *n* integers:

> (define (sq-pyr n) (sum (map square (range 1 (+ n 1))))) > (sq-pyr 15) 1240

That won’t work when *n* is a million; it takes too much time and space. Fortunately, there is a formula for computing the *n*th square pyramidal number; you might want to look at A000330:

> (define (sq-pyr n) (/ (+ (* 2 n n n) (* 3 n n) n) 6)) > (sq-pyr 15) 1240 > (sq-pyr 1000000) 333333833333500000

You can run the program at http://ideone.com/k1sJUd.

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: