## Archive for May, 2009

### Programming Praxis – Double Transposition Cipher

May 30, 2009

Yesterday’s Programming Praxis problem is about the double transposition cypher.  Our target is 34 lines (the length of the provided solution minus prelude code, blank lines and testing), so let’s dive in.

Our import:

`import GHC.Exts`

First up we have zipSort, which I also used in the Rail-Fence Cipher assignment:

```zipSort :: Ord a => [a] -> [b] -> [b]
zipSort ks = map snd . sortWith fst . zip ks```

Unsort is the opposite of zipSort. It takes a list sorted with zipSort and returns it to its original state:

```unsort :: Ord a => [a] -> [b] -> [b]
unsort ks xs = zipSort (zipSort ks [1..length xs]) xs```

With that out of the way, encrypting and decrypting becomes trivial. We zip the key with [0..] to avoid mix ups if there are identical letters in the key.

```encrypt :: Ord a => [a] -> [b] -> [b]
encrypt key = zipSort (cycle \$ zip key [0..])

decrypt :: Ord a => [a] -> [b] -> [b]
decrypt key = unsort (cycle \$ zip key [0..])```

And naturally we have to test if everything works ok:

```main :: IO ()
main = do
print . encrypt "STRIPE" \$ encrypt "COACH" "PROGRAMMINGPRAXIS"
print . decrypt "COACH" \$ decrypt "STRIPE" "GNPAPARSRIMOIXMGR"```

1 import and 4 lines of code. That will do nicely.

### Programming Praxis – Word Search Solver

May 26, 2009

Today’s Programming Praxis problem is about word search solvers. The provided solution is 77 lines, so let’s see if we can improve on that.

Our imports:

```import Data.List
import Data.Map (Map, fromList, member, keys, (!))
import Text.Printf
```

First let’s define the 8 directions that we can search in. The puzzle is going to be represented as a Map with a tuple of Ints as the key, so the directions are functions for transforming these keys.

```dirs :: [(String, (Int, Int) -> (Int, Int))]
dirs = zip ["down right", "up right", "right", "down left",
"up left", "left", "down", "up"] \$
[\(x,y) -> (x+h, y+v) | h <- [1,-1,0], v <- [1,-1,0]]
```

We’re going to enter the puzzle as a list of strings, but since that would make access an O(n2) operation we’re going to turn it into a Map instead, since that gives us O(log n2) access.

```toGrid :: [[a]] -> Map (Int, Int) a
toGrid = fromList . concat .
zipWith (\y -> zipWith (\x c -> ((x,y), c)) [1..]) [1..]
```

Next we need a function to check whether the search word appears at the given position in the given direction.

```found :: (Eq a, Ord k) => k -> (k -> k) -> Map k a -> [a] -> Bool
found pos dir g w = isPrefixOf w . map (g !) .
takeWhile (flip member g) \$ iterate dir pos
```

Finding the location and direction of a search word is then simply a matter of checking every direction for every position:

```findWord :: Map (Int, Int) Char -> String -> String
findWord g w = head [printf "%s row %d column %d %s" w y x ds |
(x,y) <- keys g, (ds, dir) <- dirs,
found (x,y) dir g w]
```

That’s all we need, so let’s test if it works.

```puzzle :: [String]
puzzle = ["FYYHNRD",
"RLJCINU",
"AAWAAHR",
"NTKLPNE",
"CILFSAP",
"EOGOTPN",
"HPOLAND"]

main :: IO ()
main = mapM_ (putStrLn . findWord (toGrid puzzle)) \$ words
"ITALY HOLLAND POLAND SPAIN FRANCE JAPAN TOGO PERU"
```

And indeed it does, using less than half the lines of the provided solution (almost a third if you ignore the lines required to define the puzzle). That will do nicely.

### Programming Praxis – The Next Palindrome

May 22, 2009

Today’s Programming Praxis problem is about palindromic numbers, i.e. numbers that read the same backwards and forwards, such as 1001 or 35753. As mentioned in the comments, the provided solution is not particularly elegant. Let’s see if we can do better.

First our import:

`import Data.List.Split`

This is all we need to get the next palindrome. The reason I’m returning a String instead of an Integer is because converting a one million-digit string is not a very fast process. The read function from the prelude takes forever, and a custom function took about 4 seconds. Since nextPalindrome (10^1000000) takes about 2 seconds, that would triple the execution time.

```nextPalindrome :: Integer -> String
nextPalindrome n = palindrome r where
l = length . show \$ n + 1
[s, m] = splitPlaces [div (l - 1) 2, 2 - mod l 2] \$ show n
palindrome x = x ++ drop (mod l 2) (reverse x)
r = if head m > last m then s ++ [head m]
else take (div (l + 1) 2) . show \$ div n (10 ^ div l 2) + 1```

Does it work? Let’s test. I’m only checking the first 10^5 instead of 10^6 numbers here so codepad doesn’t time out, but naturally it works correctly for all numbers.

```main :: IO ()
main = do
print \$ map nextPalindrome [0, 88, 808, 1999, 2133, 9999, 99999]
print \$ (takeWhile (<= 10^5) \$ iterate (read . nextPalindrome) 0)
== filter (\x -> show x == reverse (show x)) [0..10^5]```

Seven times shorter than the provided solution and about half the length of the other Haskell solution. Good enough for me.

### Programming Praxis – Fermat’s Method

May 19, 2009

Today’s Programming Praxis problem brings us yet another factorization problem, this time using Fermat’s method. Originally I tried using the algorithm from the assignment, but that resulted in pretty much a direct translation of his Scheme solution, since I couldn’t find a way to express it that was both more elegant and not slower. Since that would be boring, we’re going to use a slightly different approach to Fermat’s method that is slightly shorter and faster.

First we need a way to take the square root of integers:

```isqrt :: Integer -> Integer
isqrt = ceiling . sqrt . fromIntegral```

What this algorithm does is check whether x2 – n is square for any xs between sqrt(n) and n. If so, we return the factors of n (which are further factorized if needed).

```factor :: Integer -> [Integer] -> [Integer]
factor _ []     = []
factor n (x:xs) | y*y /= x*x - n = factor n xs
| x - y == 1     = [x + y]
| otherwise      = concatMap factors [x - y, x + y]
where y = isqrt (x*x - n)```

Like the Scheme solution, we divide by 2 until we have an odd number to work with, since Fermat’s method only works on odd numbers.

```factors :: Integer -> [Integer]
factors n | even n    = 2 : factors (div n 2)
| otherwise = factor n [isqrt n..n]```

And, as usual, we test our code:

```main :: IO ()
main = print \$ factors 600851475143```

Nine lines of code. Not bad.

### Programming Praxis – Cellular Automata

May 15, 2009

Today’s Programming Praxis problem is about cellular automata. Let’s dive in.

As usual, our imports:

```import Data.Bits
import Data.List```

Programming Praxis’ author uses ones and zeros to represent on and off cells. That works, but you run the risk of accidentally having another number somewhere and crashing your program. We Haskell programmers like our type safety, so we’re going to use booleans. First we need a function to get the successor of a group of cells:

```successor :: Int -> [Bool] -> Bool
successor r bs = testBit r . sum \$
zipWith (shiftL . fromEnum) (reverse bs) [0..]```

Using that, we can generate the next row as follows:

```nextRow :: Int -> [Bool] -> [Bool]
nextRow r xs = map (successor r) . take (length xs) . transpose .
take 3 . iterate (drop 1) \$ [head xs] ++ xs ++ [last xs]```

Since lists of booleans are not that easy to read, let’s use spaces and Xs.

```displayRow :: [Bool] -> String
displayRow = intersperse ' ' . map (\b -> if b then 'X' else ' ')```

Each run starts with one active cell and simply consists of repeating nextRow as often as needed.

```cells :: Int -> Int -> [String]
cells r h = map displayRow . take (h + 1) . iterate (nextRow r) \$
replicate h False ++ [True] ++ replicate h False```

And finally we test out program. Let’s make a lovely Sierpinski triangle:

```main :: IO ()
main = mapM_ putStrLn \$ cells 82 15```

Done. 10 lines of code. Not bad.

### Programming Praxis – Loan Amortization

May 12, 2009

Today’s Programming Praxis problem is about loan amortization. For those of you who, like me, had no idea how to calculate that, here‘s a good step by step guide. Let’s begin.

First our import:

`import Text.Printf`

With a bit of simplification we can reduce the step-by-step guide to this:

```amortize :: Fractional a => a -> a -> Int -> [(Int, a, a, a)]
amortize _ _ 0 = []
amortize b r l = (l - 1, prn, int, b') : amortize b' r (l - 1)
where int = b * r
prn = int / ((1 + r) ^ l - 1)
b' = b - prn```

Creating a table is then simply a matter of printing the header, starting month and formatting the output:

```amortization :: Double -> Double -> Int -> IO ()
amortization b r l = do
putStrLn "Left Principal Interest   Balance"
mapM_ table \$ (l, 0, 0, b) : amortize b (r / 12 / 100) l
where table (m, p, i, t) = printf "%4d %9.2f %8.2f %9.2f\n" m p i t```

And to test our function:

```main :: IO ()
main = amortization 10000 7 36```

And that’s all there is to it.

### Programming Praxis – Wheel Factorization

May 8, 2009

In today’s Programming Praxis problem we’re supposed to factor numbers both the naive way and using Wheel factorization. Let’s go.

First, let’s define a function that finds factors using trial division:

```tdFactors :: Integer -> [Integer] -> [Integer]
tdFactors _ [] = []
tdFactors r (x:xs) | x * x > r    = [r]
| mod r x == 0 = x : tdFactors (div r x) (x:xs)
| otherwise    = tdFactors r xs```

Using that, the function for naive factorization is trivial:

```factorNaive :: Integer -> [Integer]
factorNaive n = tdFactors n [2..]```

Now for the wheel factorization. First we define a function that tells us which columns of our wheel to use given the circumference.

```spokes :: Integer -> [Bool]
spokes c = map ((== 1) . gcd c) [1..c]```

To construct our wheel, we keep repeating this pattern and use only the numbers from the correct columns (and correct the start of the wheel).

```wheel :: [Integer] -> [Integer]
wheel ps = (ps ++) . drop 1 . map snd . filter fst \$
zip (cycle . spokes \$ product ps) [1..]```

The function to do wheel factorization is practically identical to the naive one:

```factorWheel :: Integer -> [Integer]
factorWheel n = tdFactors n \$ wheel [2,3,5,7]```

And of course we need to test if everything works:

```main :: IO ()
main = do print \$ factorNaive 600851475143
print \$ factorWheel 600851475143```

And we’re done. Pretty simple.

### Programming Praxis – Priority Queues

May 5, 2009

Today’s Programming Praxis problem is about priority queues. Specifically, we have to implement one using a Leftist Heap.

We define a priority queue as follows. It’s basically a binary tree, but with an extra field in which we store the rank.

`data PQueue a = Node Int a (PQueue a) (PQueue a) | Empty`

Empty nodes have rank 0.

```rank :: PQueue a -> Int
rank Empty          = 0
rank (Node r _ _ _) = r```

A convenience function for node creation that calculates the rank automatically:

```node :: a -> PQueue a -> PQueue a -> PQueue a
node i l r = if rank l > rank r then node i r l else Node (1 + rank r) i l r```

Two priority queues can be merged as follows:

```merge :: (a -> a -> Bool) -> PQueue a -> PQueue a -> PQueue a
merge _ Empty q = q
merge _ q Empty = q
merge p l@(Node _ il _ _) r@(Node _ ir lc rc) =
if p ir il then node ir lc (merge p l rc) else merge p r l```

To insert an item into a priority queue we make a new queue out of it and merge it into our original queue.

```insert :: (a -> a -> Bool) -> a -> PQueue a -> PQueue a
insert p i = merge p (node i Empty Empty)```

We convert a list to a priority queue by inserting all the items into an empty queue.

```fromList :: (a -> a -> Bool) -> [a] -> PQueue a
fromList p = foldr (insert p) Empty```

And to do the opposite we keep taking the root of the queue and merging its branches.

```toList :: (a -> a -> Bool) -> PQueue a -> [a]
toList _ Empty          = []
toList p (Node _ i l r) = i : toList p (merge p l r)```

With these functions, we can easily sort a list on priority by converting it to a priority queue and back.

```pqSort :: (a -> a -> Bool) -> [a] -> [a]
pqSort p = toList p . fromList p```

And finally we test if everything works ok.

```main :: IO ()
main = print \$ pqSort (<) [3, 7, 8, 1, 2, 9, 6, 4, 5]```

30 lines counting white space. Not bad.

### Programming Praxis – Primality Checking

May 1, 2009

Today’s Programming Praxis problem is about checking whether or not a number is prime. We’re supposed to implement the Miller-Rabin primality test, which is probably the most well-known algorithm to do this. Let’s go.

First the imports:

```import Control.Arrow
import Data.Bits
import Data.List
import System.Random```

The Miller-Rabin algorithm calls for random numbers. Since I didn’t want to make isPrime an IO function, we have to pass it a random number generator as well.

```isPrime :: Integer -> StdGen -> Bool
isPrime n g =
let (s, d) = (length *** head) . span even \$ iterate (flip div 2) (n - 1)
xs = map (expm n d) . take 50 \$ randomRs (2, n - 2) g
in all (\x -> elem x [1, n - 1] ||
any (== n - 1) (take s \$ iterate (expm n 2) x)) xs```

You may have noticed the expm in there. In theory all you have to do is mod (random_number ^ d) n. However, in our test case of 2^89 – 1, both the random number and d can easily be 10 digits. And 10-digit-number ^ 10-digit-number is a very big number, which takes a while to calculate. As a result, your algorithm is going to be terribly slow, as I discovered in my initial attempt. The solution lies in using modular exponentiation. This is a much faster way of calculating mod (a ^ b) m. In Scheme this function is apparently in the standard Prelude, but I couldn’t find a Haskell library that has it defined. No matter, we’ll just make our own:

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

And there’s our primality checking code. To test it, use:

```main :: IO ()
main = print . isPrime (2 ^ 89 - 1) =<< getStdGen```