Posts Tagged ‘bernouilli’

Programming Praxis – Sums Of Powers

February 11, 2011

In today’s Programming Praxis exercise, our goal is to implement an algorithm that calculates the bernoulli numbers and one that uses them to quickly calculate the sum of the mth power of numbers 1 through n. Let’s get started, shall we?

A quick import:

import Data.Ratio

To calculate the Bernouilli numbers I initially used the naive version, which simply uses the given mathematical formula. This is quick enough for the test case of 1000 numbers, but too slow for the test case that has a million, so we have to do some memoization. A closer look at the formula reveals that any row in the table depends only on the previous row. Since for the end result we are only interested in the last row, we can use iterate to produce the rows of the table. The value of a given column depends only on the number directly above it and the one to the upper right, so we can use a simple zip to calculate the new row.

a :: (Integral a, Integral b) => a -> a -> Ratio b
a i j = iterate (\xs -> zipWith (*) [1..] $ zipWith (-) xs (tail xs))
                (map (1 %) [1..]) !! fromIntegral i !! fromIntegral j

With this function calculating the Bernoulli numbers is trivial.

bernoullis :: (Integral a, Integral b) => a -> [Ratio b]
bernoullis upto = map (flip a 0) [0..upto]

For the algorithm we also need to calculate binomial coefficients, i.e. the amount of different ways you can choose k objects from a group of size n.

choose :: Integral a => a -> a -> Ratio a
choose n k = product [1..n] % (product [1..n-k] * product [1..k])

And some more executable math for the function that calculates the sum of powers.

s :: Integral a => a -> a -> Ratio a
s m n = 1 % (m+1) * sum [choose (m+1) k * a k 0 * (n%1)^(m+1-k) | k <- [0..m]]

We have one test case to test if the algorithm works correctly and one to judge the speed.

main :: IO ()
main = do print $ bernoullis 6 == [1, 1%2, 1%6, 0, -1%30, 0, 1%42]
          print $ s 10 1000 == 91409924241424243424241924242500
          print $ s 100 1000000

The program runs in about 150-170 ms, so we get the same speed as the Scheme version. Good enough for me.