Programming Praxis – Sieve Of Euler

In today’s Programming Praxis exercise, our goal is to implement the Sieve of Euler to generate primes. Let’s get started, shall we?

Some imports:

import Data.List.Ordered
import Data.Time.Clock

Euler’s sieve is pretty simple: just filter out all remaining multiples of each successive prime number. We add the minor optimization of only considering odd numbers.

primes_euler :: [Integer]
primes_euler = 2 : euler [3,5..] where
    euler ~(x:xs) = x : euler (minus xs $ map (* x) (x:xs))

To compare, we use the Sieve of Eratosthenes implementation found here:

primes_eratosthenes :: [Integer]
primes_eratosthenes = 2 : 3 : sieve (tail primes_eratosthenes) 3 [] where
    sieve ~(p:ps) x fs = let q=p*p in
        foldr (flip minus) [x+2,x+4..q-2] [[y+s,y+2*s..q] | (s,y) <- fs]
        ++ sieve ps q ((2*p,q):[(s,q-rem (q-y) s) | (s,y) <- fs])

A quick and dirty benchmarking function:

benchPrimes :: Num a => [a] -> IO ()
benchPrimes f = do start <- getCurrentTime
                   print . sum $ take 10000 f
                   print . flip diffUTCTime start =<< getCurrentTime

Time to run the benchmarks.

main :: IO ()
main = do benchPrimes primes_eratosthenes
          benchPrimes primes_euler

The sieve of Eratosthenes runs in about 0.02 seconds, while Euler’s sieve takes around 6.6 seconds, which means this implementation of Euler’s Sieve is about 300-400 times slower. Clearly you don’t want to use this version in practice.

Programming Praxis – Two Factoring Games

In today’s Programming Praxis exercise, our goal is to implement two algorithms related to prime factors; one to calculate the home prime of a number and one to generate the Euclid-Mullin sequence. Let’s get started, shall we?

Some imports:

import Data.List
import Data.Numbers.Primes

To calculate the home prime we repeatedly concatenate the digits of the prime factors until the number is prime.

homePrime :: Integer -> Integer
homePrime = head . filter isPrime .
            iterate (read . concatMap show . primeFactors)

The Eulic-Mullin sequence can be written as a basic unfold.

euclidMullin :: [Integer]
euclidMullin = unfoldr (\p -> let a = head (primeFactors $ p + 1)
                              in  Just (a, p * a)) 1

Naturally, we need to test if we did everything correctly:

main :: IO ()
main = do print $ homePrime 99 == 71143
          print $ take 10 euclidMullin ==
                  [2,3,7,43,13,53,5,6221671,38709183810571,139]

Looks like everything is working properly.

Programming Praxis – Google Code Jam Qualification Round Africa 2010

In today’s Programming Praxis exercise, our goal is to solve three Google Code Jam assignments. Let’s get started, shall we?

A quick import:

import Data.List

For the store credit exercise, we simply test all unique combinations of two items to see which add up to the correct total and return the first one.

storeCredit :: Num a => a -> [a] -> (Int, Int)
storeCredit c xs = head [ (i, j) | (i, a) <- zip [1..] xs
                        , (j, b) <- drop i $ zip [1..] xs, a + b == c]

Reversing the words in a sentence is trivial: make a list of words, reverse it and assemble them together again.

reverseWords :: String -> String
reverseWords = unwords . reverse . words

The T9 assignment is a bit more complicated. First, we find the correct sequence of digits for each character of the input. Then we add spaces between each consecutive group of identical digits, except between consecutive zeros.

t9 :: String -> String
t9 s = unwords =<< groupBy (\(a:_) (b:_) -> a == b && a > '0') (map f s)
    where f c = maybe "0" id . lookup c . concat $ zipWith zip
                    (words "abc def ghi jkl mno pqrs tuv wxyz")
                    [[replicate n d | n <- [1..]] | d <- ['2'..]]

Some tests to see if we made any mistakes:

main :: IO ()
main = do print $ storeCredit 100 [5,75,25] == (2,3)
          print $ storeCredit 200 [150,24,79,50,88,345,3] == (1,4)
          print $ storeCredit 8 [2,1,9,4,4,56,90,3] == (4,5)
          print $ reverseWords "this is a test" == "test a is this"
          print $ reverseWords "foobar" == "foobar"
          print $ reverseWords "all your base" == "base your all"
          print $ t9 "hi" == "44 444"
          print $ t9 "yes" == "999337777"
          print $ t9 "foo  bar" == "333666 6660022 2777"
          print $ t9 "hello world" == "4433555 555666096667775553"

Nope, looks like everything is working properly.

Programming Praxis – Sums Of Powers

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.

Programming Praxis – Excel Columns

In today’s Programming Praxis exercise our task is to write two functions to convert numbers to and from Excel’s letter-based column names. Let’s get started, shall we?

A quick import:

import Data.Char

The easiest way to covert numbers to Excel notation is recursively. I also wrote a version that uses an unfold to mirror the single fold below, but this version is more concise.

toExcel :: Int -> String
toExcel 0 = ""
toExcel n = toExcel d ++ [chr $ m + 65] where (d,m) = divMod (n - 1) 26

Converting back to numbers is even easier and can be done in a single fold.

fromExcel :: String -> Int
fromExcel = foldl (\a x -> 26 * a + ord x - 64) 0

Some test cases:

main :: IO ()
main = do print $ toExcel   1 == "A"
          print $ toExcel  26 == "Z"
          print $ toExcel  27 == "AA"
          print $ toExcel 256 == "IV"
          print $ fromExcel  "A" ==   1
          print $ fromExcel  "Z" ==  26
          print $ fromExcel "AA" ==  27
          print $ fromExcel "IV" == 256
          print $ all (\n -> n == (fromExcel . toExcel) n) [1..2^16]

Everything seems to be working correctly.