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.
Programming Praxis 
According to oeis.org sequence A000330…
def SPN(n): return n * (n + 1) * (n + n + 1) // 6SPN(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:
-- Falling factorial ff x n = product [x,x-1..x-n+1] -- Stirling numbers of the second kind s2 0 k = fromEnum(k == 0) s2 n k = s2(n-1)(k-1) + k*s2(n-1)k -- (B(n+1,x+1)-B(n+1,0))/(n+1), simplified -- cf. http://oeis.org/A103438 powersum n x = sum (map f [1..n]) where f k = (s2 n k)*(ff(x+1)(k+1) `div` (k+1)) main = (print $ powersum 2 1000000) >> (print $ map (powersum 1) [1..15]) >> (print $ map (powersum 2) [1..15]) >> (print $ map (powersum 3) [1..15])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:
-- Falling factorial ff x n = product [x,x-1..x-n+1] -- Stirling numbers of the second kind s2seq = (1:[0,0..]) : map f s2seq where f s = zipWith (+) (zipWith (*) [0..] s) (0:s) s2 n k = s2seq!!fromIntegral(n)!!fromIntegral(k) -- (B(n+1,x+1)-B(n+1,0))/(n+1), simplified -- cf. http://oeis.org/A103438 powersum :: Integer -> Integer -> Integer powersum n x = sum (map f [0..n]) where f k = (s2 n k)*(ff(x+1)(k+1) `div` (k+1)) main = (mapM print (take 10 (map (take 10) s2seq))) >> (print $ powersum 2 1000000) >> (print $ powersum 2 10000000000) >> (print $ powersum 10 1000) >> (print $ map (powersum 1) [1..15]) >> (print $ map (powersum 2) [1..15]) >> (print $ map (powersum 3) [1..15])Final version:
-- falling factorial ffs x = x : map (x*) (ffs(x-1)) -- Stirling numbers s2s 1 = [1] s2s n = zipWith (+) (zipWith (*) [1..n+1] (s++[0])) (0:s) where s = s2s(n-1) powersum :: Integer -> Integer -> Integer powersum n x = sum (zipWith3 f [1..n] (s2s n) (ffs x)) where f k s2 ff = (x+1) * ff `div` (k+1) * s2In 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: