Archive for March, 2010

Programming Praxis – Passover

March 30, 2010

In today’s Programming Praxis exercise we have to write functions to calculate the dates of the Jewish holidays Rosh Hashanah and Passover. Let’s get started.

As usual, some imports:

import Data.Time.Calendar
import Data.Time.Calendar.WeekDate
import Data.Fixed

We start with the function to calculate on which day of September (or October) Rosh Hashanah falls. Unfortunately, because in Haskell we can’t just add integers and floats together, there is a bit of inelegance due to the two fromIntegrals. Nothing too serious, though.

roshDay :: Integer -> Integer
roshDay year | elem d [3,5,7] = n + 1
             | d == 1 && f >= 23269/25920 && r > 11 = n + 1
             | d == 2 && f >= 1367/2160   && r > 6 = n + 2
             | otherwise = n
    where
    y = fromIntegral year
    g = mod' y 19 + 1
    r = mod' (12 * g) 19
    a // b = fromIntegral $ floor (a / b)
    (n, f) = properFraction (y // 100 - y // 400 - 2 +
                             765433/492480 * r + mod' y 4 / 4 -
                             (313 * y + 89081) / 98496)
    (_,_,d) = toWeekDate . addDays n $ fromGregorian year 8 31

Once we have that, calculating the actual dates of Passover and Rosh Hashanah is a simple matter of adding the found number of days to the correct starting dates.

roshHashanah :: Integer -> Day
roshHashanah year = addDays (roshDay year) (fromGregorian year 8 31)

passover :: Integer -> Day
passover year = addDays (roshDay year) (fromGregorian year 3 21)

A quick test demonstrates we’re getting the correct dates.

main :: IO ()
main = do print $ roshHashanah 2010
          print $ passover 2010

Programming Praxis – Texas Hold ‘Em

March 23, 2010

In today’s Programming Praxis exercise we have to write a program to rank poker hands. The provided Scheme solution clocks in at 62 lines. Let’s see if we can’t bring that number down a little.

As usual, some imports:

import Data.Char
import Data.List
import qualified Data.List.Key as K

We’ll need to convert the two-character strings that are used to define the card to something that’s easier to work with.

toCard :: String -> (Int, Char)
toCard ~[v,s] = maybe undefined (flip (,) $ toUpper s) .
                lookup v $ zip "23456789TJQKA" [2..]

If a hand is a flush, grouping on the suit should produce a list with only one element.

flush :: [(Int, Char)] -> Bool
flush = (== 1) . length . K.group snd

If a hand is a straight, it should consist of ascending values or be ace through five.

straight :: [Int] -> Bool
straight xs = (l:r) == [2,3,4,5,14] || isPrefixOf r [l + 1..]
              where (l:r) = reverse $ map fromEnum xs

As with flushes, we use grouping to find how many cards we have with the same value. This function assumes that the list of values is already sorted.

same :: [Int] -> [Int]
same = map length . group

Calculating the rank of a hand then becomes a fairly simple matter of using the functions above to find the highest possible rank. We also add the sorted card values to the result so that multiple hands with the same rank will be compared on highest card value.

UPDATE: The original version had a bug where n-of-a-kinds were compared on the highest value of the other cards. I altered the definition of s to remedy this.

rank :: [String] -> (Int, [Int])
rank xs | straight s && flush cs     = (9, s)
        | elem 4 $ same s            = (8, s)
        | same s == [3,2]            = (7, s)
        | flush cs                   = (6, s)
        | straight s                 = (5, s)
        | elem 3 $ same s            = (4, s)
        | elem 2 . delete 2 $ same s = (3, s)
        | elem 2 $ same s            = (2, s)
        | otherwise                  = (1, s)
        where s = concat . reverse . K.sort (\x -> (length x, head x)) .
                  group . sort $ map fst cs
              cs = map toCard xs

To determine the best hand, we need a way to look at every possible hand that can be chosen from the 7 available cards.

choose :: Int -> [a] -> [[a]]
choose 0 _      = [[]]
choose _ []     = []
choose n (x:xs) = map (x: ) (choose (n-1) xs) ++ choose n xs

The best hand is simply the one with the highest rank.

bestHand :: [String] -> [String]
bestHand = K.maximum rank . choose 5

Nothing left but to test if everything works correctly:

main :: IO ()
main = do
    mapM_ print [fst (rank ["AH", "KH", "QH", "JH", "TH"]) == 9
                ,fst (rank ["7H", "7C", "3H", "7S", "7D"]) == 8
                ,fst (rank ["TH", "JC", "TS", "JD", "TC"]) == 7
                ,fst (rank ["4H", "7H", "AH", "KH", "9H"]) == 6
                ,fst (rank ["AH", "2C", "3S", "4D", "5H"]) == 5
                ,fst (rank ["9C", "4S", "KD", "9D", "9H"]) == 4
                ,fst (rank ["6D", "6C", "8H", "TD", "8D"]) == 3
                ,fst (rank ["9C", "3S", "4D", "7C", "3D"]) == 2
                ,fst (rank ["4C", "KD", "8S", "6D", "2D"]) == 1]
    print $ bestHand ["AC", "JD", "8S", "8C", "AS", "3D", "4D"]
    --Added after Phil's bug report
    print $ rank ["AC", "7D", "7S", "7H", "7C"] <
            rank ["4C", "8D", "8S", "8H", "8C"]

Yup. And at around a third of the line count of the Scheme version, that will do just fine.

Programming Praxis – Extending Pollard’s P-1 Factorization Algorithm

March 19, 2010

In today’s Programming Praxis exercise we need to write an improved version of a factorization algorithm. I was on vacation when the original exercise was posted, so let’s see what we can do with it.

As usual, some imports:

import Data.Bits
import Data.List

We need the same expm function we have used in several previous exercises. Alternatively we could use the expmod function from Codec.Encryption.RSA.NumberTheory, but it’s a lot slower than this version.

expm :: Integer -> Integer -> Integer -> Integer
expm b e m = foldl' (\r (b', _) -> mod (r * b') m) 1 .
             filter (flip testBit 0 . snd) .
             zip (iterate (flip mod m . (^ 2)) b) .
             takeWhile (> 0) $ iterate (`shiftR` 1) e

The scheme solution has some duplication in it: the start of the pollard1 and pollard2 functions is nearly identical. Since programmers hate repeating themselves, let’s factor that out into a separate function.

pollard :: (Integer -> t) -> (Integer -> t) -> Integer -> Integer -> t
pollard found notFound n b1 = f 2 2 where
    f a i | i < b1         = f (expm a i n) (i + 1)
          | 1 < d && d < n = found d
          | otherwise      = notFound a
          where d = gcd (a - 1) n

pollard1 then becomes very simple: if we don’t find anything we stop, otherwise we return the result.

pollard1 :: Integer -> Integer -> Maybe Integer
pollard1 = pollard Just (const Nothing)

pollard2 is a bit more involved, because we now have an extra step if we don’t find anything. The structure of this is very similar to the pollard function, but there are enough differences that it’s not worth the bother of abstracting it.

pollard2 :: Integer -> Integer -> Integer -> Maybe (String, Integer)
pollard2 n b1 b2 = pollard (Just . (,) "stage1") (f b1) n b1 where
    f j a | j == b2        = Nothing
          | 1 < d && d < n = Just ("stage2", d)
          | otherwise      = f (j + 1) a
          where d = gcd (expm a j n - 1) n

And of course the test to see if everything’s working correctly:

main :: IO ()
main = do print $ pollard1 15770708441 150
          print $ pollard1 15770708441 180
          print $ pollard2 15770708441 150 180

Programming Praxis – Traveling Salesman: Nearest Neighbor

March 16, 2010

In today’s Programming Praxis exercise we have to implement a significantly faster algorithm for the traveling salesman problem than in the previous exercise. Let’s get started, shall we?

As usual, some imports:

import Control.Monad
import Data.List
import qualified Data.List.Key as K
import System.Random

The functions for calculating distance and total distance are largely the same as in the previous exercise, though because I’ve switched to using lists of integers we now need an extra fromIntegral, and repeating the first point in order to complete the loop has been moved to the cost function.

dist :: (Int, Int) -> (Int, Int) -> Float
dist (x1, y1) (x2, y2) = sqrt (f (x1 - x2) ** 2 + f (y1 - y2) ** 2)
    where f = fromIntegral

cost :: [(Int, Int)] -> Float
cost xs = sum $ zipWith dist xs (tail xs ++ xs)

Generating a series of random points is a bit longer than it could be because we have to make sure all the points are unique.

randomPoints :: Int -> IO [(Int, Int)]
randomPoints n = f [] where
    f ps = if length ps == n then return ps else
           do p <- liftM2 (,) rnd rnd
              if elem p ps then f ps else f (p:ps)
    rnd = randomRIO (0, 10 * n)

Determining the tour to take using the nearest neighbor algorithm is not that difficult. Again, we index the points for similarity to the Programming Praxis solution, not because it is needed.

nearTour :: [(Int, Int)] -> [(Integer, (Int, Int))]
nearTour = f . zip [0..] where
    f [] = []
    f [x] = [x]
    f ((i,p):ps) = (i,p) : f (nxt : delete nxt ps) where
        nxt = K.minimum (dist p . snd) ps

To test, we check both a random set of points, as well as the set from the Programming Praxis solution.

main :: IO ()
main = do rndTour <- fmap nearTour $ randomPoints 25
          print (cost $ map snd rndTour, rndTour)
          let praxisTour = nearTour
                [(139, 31),( 41,126),(108, 49),(112,193),(179,188),
                 (212, 24),(245, 50),(167,187),(159,236),(185, 78),
                 ( 27, 63),(101,188),(195,167),( 30, 10),(238,110),
                 (221, 60),( 27,231),(146, 67),(249,172),( 36, 71),
                 ( 37,203),(118, 38),(241,226),(197, 29),(220,186)]
          print (cost $ map snd praxisTour, praxisTour)

We get the same tour as the Programming Praxis solution (Ok, the reverse to be exact. Again, this doesn’t matter and I think starting with the first point is more logical), and at a third of the line count, so I think we can call this one done.

Programming Praxis – Traveling Salesman: Brute Force

March 12, 2010

In today’s Programming Praxis exercise we have to implement a brute-force algorithm for solving the well-known traveling salesman algorithm. The provided Scheme solution clocks in at 28 lines. Let’s see if we can come up with something a tad more compact.

A quick import or two:

import Data.List
import qualified Data.List.Key as K

In order to calculate the total distance traveled, we need to calculate the distance between two points.

dist :: Floating a => (a, a) -> (a, a) -> a
dist (x1, y1) (x2, y2) = sqrt ((x1 - x2) ** 2 + (y1 - y2) ** 2)

Calculating all possible tours is just a matter of generating all the permutations of the points, making sure to duplicate the first element of a tour at the end in order to get back home. We also add an index to the points to make the result a bit easier to read (and mainly because the Scheme solution does it too).

tours :: [b] -> [[(Int, b)]]
tours = map (\(x:xs) -> x:xs ++ [x]) . permutations . zip [0..]

Calculating the total cost of a tour is a simple matter of summing the distances between every consecutive pair of points. This function resembles the typical definition of the Fibonacci sequence, since that works in much the same way.

cost :: Floating a => [(b, (a, a))] -> a
cost xs = sum $ zipWith dist xs (tail xs)

Finally, showing the shortest path is a simple matter of generating all the tours, taking the one with the lowest cost, showing the indices of the points and getting rid of the duplicated starting point.

shortestPath :: [(Double, Double)] -> [Int]
shortestPath = init . map fst . K.minimum cost . tours

As usual, a test to see if everything is working correctly:

main :: IO ()
main = print $ shortestPath [(5, 2), (19, 13), (4, 8), (6, 32),
                             (23, 7), (57, 54), (55, 8), (70, 59)]

You’ll notice that the result is different from the Scheme solution. Specifically, it’s reversed and cycled a few positions. Since we’re walking a loop, this means that we’re walking clockwise instead of counterclockwise and starting in a different town. A trivial difference, since the route is the same and can be reversed and cycled at will. The reason is the order in which the permutations are generated.

That brings us to 4 lines total, an 85% reduction compared to the Scheme solution. That’s not half bad in my book.

Programming Praxis – Lexicographic Permutations

March 9, 2010

In today’s Programming Praxis exercise we have to generate all the lexicographic permutations of a list. Let’s get started, shall we?

Some imports:

import Data.List
import qualified Data.List.Key as K

We’re going to be taking a different approach than the Scheme solution, which relies on explicit element swapping. We annotate each element of the list with how many elements are bigger than it. Next, we generate all the permutations of the sequence (not in lexicographic order). We then sort the permutations based on the annotations, reversing them because the most significant number is on the right. Because the biggest element will have an annotation of 0 and the rest are all higher, this results in a lexicographic ordering.

perms :: (a -> a -> Bool) -> [a] -> [[a]]
perms cmp xs = map fst . K.sort (reverse . snd) . map unzip $
               permutations [(x, length $ filter (cmp x) xs) | x <- xs]

And of course a test to show that our alternative approach is working correctly:

main :: IO ()
main = do print $ perms (<) [1,1]
          print $ perms (<) [1,2,3,4]

Programming Praxis – Binary Search Tree

March 5, 2010

In today’s Programming Praxis exercise we have to implement a Binary Search Tree. Let’s get started, shall we?

We need two imports:

import Control.Monad
import System.Random

The data structure is your run-of-the-mill binary tree.

data BTree k v = Node k v (BTree k v) (BTree k v) | Empty

Finding an element is pretty straightforward. Just keep taking the correct branch until we exhaust the tree or find what we want.

find :: (k -> k -> Ordering) -> k -> BTree k v -> Maybe v
find _   _ Empty          = Nothing
find cmp k (Node k' v' l r) = case cmp k k' of EQ -> Just v'
                                               LT -> find cmp k l
                                               GT -> find cmp k r

Inserting works the same way as find: move to the correct position and insert or replace the new value.

insert :: (k -> k -> Ordering) -> k -> v -> BTree k v -> BTree k v
insert _   k v Empty            = Node k v Empty Empty
insert cmp k v (Node k' v' l r) = case cmp k k' of
    EQ -> Node k v l r
    LT -> Node k' v' (insert cmp k v l) r
    GT -> Node k' v' l (insert cmp k v r)

Since the deletion algorithm calls for a random number, delete is an IO action. You can consider using unsafePerformIO to hide this (I did in my first draft), but I decided to stick with the honest, safer (though less convenient) version. Alternatively you could accept the occasional imbalance and just always start on the left.

delete :: (k -> k -> Ordering) -> k -> BTree k v -> IO (BTree k v)
delete _   _ Empty              = return Empty
delete cmp k t@(Node k' v' l r) = case cmp k k' of
    EQ -> fmap (flip deroot t . (== 0)) $ randomRIO (0,1 :: Int)
    LT -> fmap (flip (Node k' v') r) $ delete cmp k l
    GT -> fmap (      Node k' v'  l) $ delete cmp k r

For the deroot function we use a slightly different approach than the Scheme version. I’m not sure how that version deals with the case of one of the two branches being empty, but here they are explicitly included in the patterns. The rot-left and rot-right functions are rewritten as patterns.

deroot :: Bool -> BTree k v -> BTree k v
deroot _    Empty              = Empty
deroot _    (Node _ _ l Empty) = l
deroot _    (Node _ _ Empty r) = r
deroot True (Node k v l (Node rk rv rl rr)) =
    Node rk rv (deroot False $ Node k v l rl) rr
deroot _    (Node k v (Node lk lv ll lr) r) =
    Node lk lv ll (deroot True $ Node k v lr r)

Converting the search tree to a list is trivial.

toList :: BTree k v -> [(k, v)]
toList Empty          = []
toList (Node k v l r) = toList l ++ (k, v) : toList r

And, as always, a test to see if everything is working correctly:

main :: IO ()
main = do let t = foldr (uncurry $ insert compare) Empty $
                  [(n, n) | n <- [4,1,3,5,2]]
          print $ toList t
          print $ find compare 3 t
          print $ find compare 9 t
          print . toList =<< foldM (flip $ delete compare) t [4,2,3,5,1]

Programming Praxis – Goldbach’s Conjecture

March 2, 2010

In today’s Programming Praxis problem we have to test Golbach’s conjecture that every even number greater than 2 can be written as the sum of two primes. Let’s get started.

A quick import (note that both the Numbers and primes packages define this module, but only Numbers contains the isPrime function we’re using):

import Data.Numbers.Primes

The function for testing Goldbach’s conjecture is pretty simple: just look at all the primes p less than n and check if n – p is prime.

goldbach :: Integer -> (Integer, Integer)
goldbach n = head [(p, n - p) | p <- takeWhile (< n) primes, isPrime (n - p)]

A quick test reveals everything’s working correctly.

main :: IO ()
main = do print $ goldbach 28
          print $ goldbach 986332

The Scheme solution implements two versions, citing speed as a reason. I also implemented his second version, but the speed was nearly identical to the algorithm above (about 11 seconds to determine that 532 is indeed the highest first prime in the first one million numbers), and given the choice I’d much rather have one line than thirteen.