## 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.

Pages: 1 2

### 10 Responses to “Square Pyramidal Numbers”

1. Milbrae said

According to oeis.org sequence A000330…

```def SPN(n):
return n * (n + 1) * (n + n + 1) // 6
```

SPN(1000000) = 333333833333500000

2. 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.

3. Globules said

```spn :: Integer -> Integer
spn n = (n * (n+1) * (2*n+1)) `quot` 6

main :: IO ()
main = do
print \$ spn 1000000
print \$ spn 10000000000
```
```\$ time ./pyrim
333333833333500000
333333333383333333335000000000

real    0m0.011s
user    0m0.002s
sys     0m0.007s
```
4. matthew said

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])
```
```333333833333500000
[1,3,6,10,15,21,28,36,45,55,66,78,91,105,120]
[1,5,14,30,55,91,140,204,285,385,506,650,819,1015,1240]
[1,9,36,100,225,441,784,1296,2025,3025,4356,6084,8281,11025,14400]
```
5. matthew said

Better make that:

```s2 0 k = toInteger(fromEnum(k == 0))
```

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).

6. programmingpraxis said

@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.

7. matthew said

@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])

```
```[1,0,0,0,0,0,0,0,0,0]
[0,1,0,0,0,0,0,0,0,0]
[0,1,1,0,0,0,0,0,0,0]
[0,1,3,1,0,0,0,0,0,0]
[0,1,7,6,1,0,0,0,0,0]
[0,1,15,25,10,1,0,0,0,0]
[0,1,31,90,65,15,1,0,0,0]
[0,1,63,301,350,140,21,1,0,0]
[0,1,127,966,1701,1050,266,28,1,0]
[0,1,255,3025,7770,6951,2646,462,36,1]
333333833333500000
333333333383333333335000000000
91409924241424243424241924242500
[1,3,6,10,15,21,28,36,45,55,66,78,91,105,120]
[1,5,14,30,55,91,140,204,285,385,506,650,819,1015,1240]
[1,9,36,100,225,441,784,1296,2025,3025,4356,6084,8281,11025,14400]
```
8. matthew said

Final version:

```-- falling factorial
ffs x = x : map (x*) (ffs(x-1))
-- Stirling numbers
s2s 1 = 
s2s n = zipWith (+) (zipWith (*) [1..n+1] (s++)) (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) * s2
```
9. Aurelio said

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

10. bavier said

I don’t see Forth here much:

```: spn   0 0 rot  1+ 1 do i i * 0 d+ loop ;  ok
15 spn d. 1240  ok
1000000 spn d. 333333833333500000  ok
1000000000 spn d. 333333333833333333500000000  ok
```

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

```: spn'   dup 1 + ( n n+1) over 2 um* 1 m+ ( n n+1 2n+1) rot 1 m*/ rot 6 m*/ ;
15 spn' d. 1240  ok
1000000000 spn' d. 333333333833333333500000000  ok
10000000000 spn' d. 333333333383333333335000000000  ok
```