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.

Advertisements

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

    In Haskell…

    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 = [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) * 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
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: