Flashback, Part Two

Last time we found that our program to count ways to make change, while much better than the engineering class's brute force nested DO loops, isn't up to snuff.

There is a better way, it turns out. If you look at the Rosetta Code entry for this problem, which they call "Count the Coins", you'll find a number of approaches. For Haskell, there are two. The first "naive" program is essentially what we wrote (save that there's no explicit case for a one-element list; it's a little less efficient, but the recalculation may well swamp that difference).

The other Haskell program takes a different approach:

count = foldr addCoin (1:repeat 0)
        where addCoin c oldlist = newlist
              where newlist = (take c oldlist) ++ zipWith (+) newlist (drop c oldlist)

and if you look at the main for the program,

main = do
    print (count [25,10,5,1] !! 100)
    print (count [100,50,25,10,5,1] !! 100000)

you see that count returns a list that one indexes by the amount one wants to know how many ways to make change for (remember that the head of a list has index 0), reminiscent of the elegant Haskell Fibonacci sequence definition

fibs = 0 : 1 : (zipWith (+) fibs (tail fibs))

If you try it, you'll find it runs a lot faster than our program.

Consider the base case, no denominations, and then how to add a denomination. With no denominations, the only amount you can make change for is zero, so if count is correct it should hand back an infinite list whose head is 1 and all the other elements zero. That's exactly the

1 : repeat 0

that is the starting value for the fold.

So, suppose we have the list that corresponds to the first n denominations. How do we use that to generate the list for the first n + 1?

Well... if there are w ways to make change for c cents (just to have a name, we'll call the smallest denomination "cents") using the first n denominations, and the (n+1)st denomination is d cents, then there are w more ways to make change for c + i * d cents, for i in [1..] when you add the new denomination. So the new list is the current list plus the current list shifted right d places plus the current list shifted right 2 * d places plus... etc. That does look like what addCoin does. (Apologies for picking names that clash with addCoin; what they call "c", we're using "d" for.) Looking like, though, isn't proof. Back when I have proof.... and why foldr  rather than foldl?

Flashback

Back in the Cretaceous era I worked at the University of Oklahoma as a student assistant at Remote 1. OU was a [shudder] IBM big iron shop at the time, and beginning students typically got a job class that limited CPU time to thirty seconds, lest a beginner submit a program with an endless loop and chew up serious resources. (And when the computer everyone's sharing is a 370 of mid-1970s vintage, it doesn't take much to be considered serious resources. According to the Hercules FAQ, a Celeron 300 can emulate a 370 at speeds greater than a 3033, which was at least twice as fast as what OU had at the time.)

One day one of the instructors assigned his charges the following problem: count the number of ways to make change for a dollar assuming you have a generous supply of half dollars, quarters, dimes, nickels, and pennies. (The language is a bit loose; there wasn't any buying going on, just figuring out how to get coinage adding up to a dollar.) They all headed off to punch their FORTRAN programs to run using WATFIV, successor to WATFOR, the University of Waterloo FORTRAN compiler designed to give good error messages rather than worrying about optimization.

The student assistants were getting swamped with students coming in and wanting to know how to get more computer time; their change-making programs were running thirty seconds and getting killed. Every one of them had written something like this, modulo the length of that logical IF; I'm not going to bother to churn out the appropriate "continuation card"):

       IWAYS = 0
       DO 100 ICENTS = 0,100
         DO 200 INICKELS = 0,20
           DO 300 IDIMES = 0,10
             DO 400 IQUARTERS = 0,4
               DO 500 IHALFS = 0,2
                  IF (ICENTS + 5 * INICKELS + 10 * IDIMES + 25 * IQUARTERS + 50 * IHALFS .EQ. 100) IWAYS = IWAYS + 1
 500           CONTINUE
 400         CONTINUE
 300       CONTINUE
 200     CONTINUE
 100   CONTINUE
       WRITE(6,600) IWAYS
 600   FORMAT(I8,' WAYS')
       STOP
       END
       
so that no attention was paid to how much you'd accumulated in outer loops--no matter whether ICENTS was 0 or 99, it would doggedly try counts of other coins sure to push the total far over a dollar.

I wrote a recursive program in Algol W, my then language of choice, with a function that took an arbitrary array of coin denominations and amount and a main program passing the particular values for the problem the instructor set. It ran in 0.01 seconds. A friend and coworker, Raymond Schlecht (alevasholem), wrote a FORTRAN version that took less than 0.01 seconds, so from the output, which only gave run time to two decimal places, it looked like it took no time at all. We handed the listings and output to our boss, who had a chat with the instructor. I like to think the discussion was educational.

Nowadays, of course, just about everyone has computers so much faster than that 370/158 that the above brute force code would run to completion in no time at all. The old fogey in me hopes that the up and coming programmers serve time with low-end hardware, and indeed systems like Arduino mean that does happen at least some of the time. OTOH, another way to do it is to crank the input size so you still can't get away with an inefficient method, vide the small and large inputs for Code Jam problems... and we'll see that shortly as I get my comeuppance.

So, how to do it in Haskell?

waysToMakeChange :: [Int] -> Int -> Int

Obviously there's exactly one way to make change for 0.

waysToMakeChange _ 0  = 1

For a nonzero amount, you need at least one denomination!

waysToMakeChange [] _ = 0

With just one denomination, either you can do or do not, as Yoda would say.

waysToMakeChange [denom] amount
    | amount `mod` denom == 0 = 1
    | otherwise               = 0

With more than one, consider the first one. The number should be the sum of the ways to make change for the amount minus 0, 1, 2, ... of the coins of that denomination using the remaining denominations as long as the difference is non-negative.

waysToMakeChange (denom:denoms) amount
    = sum $ map (waysToMakeChange denoms) [amount, amount - denom .. 0]

That should do the trick... shouldn't it?

*Main> waysToMakeChange [50, 25, 10, 5, 1] 100
292

Do we believe that? Turns out that Rosetta Code has the same problem, only their example doesn't bother with half-dollars. The result they get for that variation is 242, and sure enough, that's what we get as well.

For completeness's sake, the Rosetta Code page shows an example with considerably larger amounts, for us it would be

waysToMakeChange [100, 50, 25, 10, 5, 1] 100000

Let's try it. (Given the output shown on Rosetta Code, we have to switch over to Integer for the result.) Now we're in the same boat as the engineering students! Even compiled, ten minutes isn't enough given the above implementation.

We need a better method. Why? This version of the code is doing a lot of recalculation. For example, if we have a half-dollar and a quarter, or three quarters, or two quarters, two dimes, and a nickel, or... we will keep evaluating, among other things, the number of ways to make change for 25 cents just using pennies. So while we have a nice recursive function that is pretty clearly correct, it's inefficient. We tried one dollar, ten dollars, and one hundred dollars, running time on the compiled program, and looking at the user line of the time output, we see

one dollar:                 0.002 s
ten dollars:                 0.009 s
one hundred dollars:  6.794 s

Definitely the larger problem isn't practical this way... but this is already long enough and has been sitting for months waiting to be published, so let's make this Part One. Stay tuned.

2010 Code Jam problem: "Store Credit"

A fellow on the haskell-beginners mailing list asked for suggestions on how to improve code he'd written to solve a 2010 Code Jam qualifying problem, "Store Credit". Here's the problem statement:

"You receive a credit C at a local store and would like to buy two items. You first walk through the store and create a list L of all available items. From this list you would like to buy two items that add up to the entire value of the credit. The solution you provide will consist of the two integers indicating the positions of the items in your list (smaller number first)."

Then the format of the input is described, and among the things we learn are:

• The prices are all integers.
• "Each test case will have exactly one solution."

But then one of the examples for which a solution is shown has a credit amount of 8 and this list of prices:

2 1 9 4 4 56 90 3

and the solution "4 5". At first one would--well, I did--think that violates the promise of only one solution, since two of item 4 or two of item 5 would also add up to the credit amount of 8. I would say that the problem statement is at least a little misleading, and should have said "...would like to buy two different items that add up to the entire value of the credit." It's arguably implied by the parenthetical "(smaller number first)", and perhaps the intent is to show the need to read specs carefully.

Clearly you want to keep around two items, or for a functional version, you want to generate two items as you trudge down the list of prices:

• an array or vector of integers with C/2 elements
• a mapping that, given a price, returns the position (or positions, given the one example) of the items with that price

The array/vector starts out as all zeroes, and for each item, you increment the element whose index is the minimum of the price and C - the price. (Ignore anything with price > C.) You've found a solution when the result of the increment is 2.

It's pretty clear how to implement this in an imperative language. We'll see what I come up with for a Haskell version.

Look and say sequence

If I don't already have the Haskell subreddit link over on the right, I'll add it ASAP.

This evening a Haskell beginner posted about some trouble he was having writing code to generate a particular sequence. I didn't catch on to the sequence he was going for, but I should have from a comment in his code:

enunBlock :: [Int] -> [Int] -- [2,2,2] -> [3,2] | [3] -> [1,3]  

Someone did catch on, though, and asked "Are you trying to make a look and say sequence?" The poster said yes... and off to Wikipedia I went.

Said sequence starts 1, 11, 21, 1211, 111221, ... and the way you get the next term is to take the current one and sort of run-length encode it. The first term would be "one one", i.e. a run of ones of length one, so the second term is 11. That in turn would be described as "two ones", hence the third term is 21, or "one two, one one", giving 1211, and so on.

If each digit were on its own line, you could get the next term of the sequence by piping it through uniq -c and playing some sed games to pull out some whitespace. (Hey, wait; will we ever have a run length needing two digits? Actually, you can bound it even more tightly than that, unless you pick a starting value that forces the issue. The result, and a reference, is in the Wikipedia article.)

We'll go along with the poster--as you can see from that declaration, he's using a [Int] for a term of the sequence, saving some hassle. If we can just write a function that takes a term and generates the next, then, if you've read much here or done much Haskell at all, then you know where we're headed for the sequence, namely iterate.

lookAndSay :: [Int] -> [[Int]]
lookAndSay xs = iterate nextTerm xs
     where nextTerm xs = ...

So, how to define nextTerm (that's what we're calling what the original poster called enunBlock)? Clearly our base case is

nextTerm [] = []

the more interesting one is

nextTerm (x:xs) = ....

and for it we're going to take advantage of span from Data.List:

nextTerm (x:xs) = length run  + 1 : x : (nextTerm leftovers)
      where (run, leftovers) = span (== x)  xs

put it all together after a pass through hlint, and you have

import Data.List

lookAndSay :: [Int] -> [[Int]]
lookAndSay = iterate nextTerm
    where nextTerm []     = []
          nextTerm (x:xs) = length run + 1 : x : nextTerm leftovers
                            where (run, leftovers) = span (== x) xs

UPDATE: Aha! Should've looked closer at what the poster was trying to do. He or she was looking to use group which takes a list and chops it up into a list of lists made up of the runs of equal elements; for example, if you hand it "balloon", you get back ["b", "a", "ll", "oo", "n"]. That does come closer to what we want, but we need the counts and to get rid of the duplicates where they exist, something like

nextTerm xs = concat [[length ys, head ys] | ys <- group xs]

or, if you're feeling pointless--er, pointfree,

nextTerm = concat . (map (\xs -> [length xs, head xs]) . group

and sure enough, that does the trick.

UPDATE: Darn it! You'd think that after having found and used sequence, I'd remember it. Make that

nextTerm = concat . map (sequence [length, head]) . group

and we'll be that much more concise. Ah, even better is

nextTerm = concatMap (sequence [length, head]) . group

and once you have it boiled down this far, why bother giving it a name? Once the smoke clears, the whole thing is just

import Data.List

lookAndSay :: [Int] -> [[Int]]
lookAndSay = iterate $ concatMap (sequence [length, head]) . group

Clearly I need to spend some time just rummaging through the available functions... and now I'll have Cream's "I Feel Free" stuck in my head until I take a stab at a filk, working title "I'm Pointfree".

Fun with pattern matching

I just noticed that I did something almost without thinking.

I'm still being lazy about dealing with inputs. I haven't gotten around to using one of the many elegant parser packages that Haskell makes possible, and the Code Jam problems I've tried so far haven't required much more than skipping the leading line--kids, ask your grandparents about the card decks they would feed to their FORTRAN programs that started with a card that said how many sets of input values followed!--and handing back a list of values per line of input afterwards.

A particular Cookie Clicker Alpha input has three values, and so I passed cookieTime a list I knew would always have three elements, so that the pattern

cookieTime [farmCost, farmRate, goal] = ...

matches it and also splits out the elements of the list and binds them to names--just the way that introductory Haskell texts always mention that you could make yourself comfortable if you're used to languages that want parameter lists parenthesized by passing tuples, but you shouldn't. I will try to wean myself away from it.

Well, drat... *spoiler alert*

I didn't do well in the qualifying round of Code Jam this year, alas. Part of it I can thank the IRS for, in a way, because we paid our annual last-minute visit to H&R Block today. But really it's my fault; I didn't arrange to have the interval during which one could enter solutions clear to devote to it.

I didn't get my improved version of Cookie Clicker Alpha done in time to count, but it did check out, so I have a moral victory at least. I probably should have devoted the time needed to do the simplest of the qualifying problems, Magic Trick, rather than just sketching out the solution.

A quick summary of the problem: in the problem's version of the game "Cookie Clicker" (inspired by a real game of that name) is parametrized by three values:

• X, a goal number of cookies to accumulate
• C, the cost in cookies of a "cookie farm"
• F, the rate at which one can "harvest" cookies from a cookie farm once bought

OK, actually, there's a fourth value: if you click on a giant cookie, you can gain two cookies/second. (You have to have a way to get started accumulating cookies, or you couldn't buy a cookie farm.)

One thing to note: in real life, cookies are discrete entities; a whole cookie is either there or not. In this problem, though, you should think of cookies as being produced continuously; for example, with no farms, you can lay claim to two thirds of a cookie in a third of a second.

So, the problem as posed: given X, C, and F, figure out the soonest you can accumulate the goal of X cookies. Of course, cookies you spend on farms don't count towards X, though of course it will usually pay off to buy farms. (Not necessarily, though! If X < C, best to just click away.)

As usual for Code Jam, the input is a line with the number of cases of the problem to solve, followed by lines with the data for the cases. For this one, it's one line per case, with C, F, and X (which are "real" numbers) on the line in that order, separated by spaces. The output format? You should know by now:

Case #<number>: <time>

where <time> is the minimum time to get X cookies with the given cost and rate per cookie farm.

Magic Trick is simple enough that you don't have to worry about optimization the way you do for problems like Fair and Square last year... but Cookie Clicker has two inputs to try, and my initial solution was far too slow to handle the larger input within the time limit.

Spoilers behind the break...

Letting the compiler do it for you

A while back I posted about how looking at the type can give you hints or even tell you what a function has to do, giving the simple example of function composition and pointing at a video by Matthew Brecknell showing how you can enlist ghci to help you in that respect.

Here's another example, provoked by a comment on a reddit question, namely map. If you only knew

map :: (a -> b) -> [a] -> [b]

could you write map? The signature says "give me a function that maps as to bs and a list of as, and I'll give you back a list of bs". Well, you could be a wise guy and write

map _ _ = []

and if anyone complained, say "Hey, the empty list is a list of bs", but in your heart you'd know you were cheating. The only way you're given to get a b given an a is to use the function, so the obvious thing to write is

map f xs = [f x | x <- xs]

(Yes, you could be a wise guy at the next level up and write

map f xs = reverse [f x | x <- xs]

but again, your conscience should bother you.)

Well, it turns out that at least sometimes, a compiler can do some of that kind of reasoning for you. Not in Haskell (though there's been a proposal), but in languages that support what are called "dependent types". Dependent types depend on a value; the canonical example is parametrizing lists by their length. That lets the language's type checker catch more errors.

"Hey, wait!" you say, "Does that mean I can't write my cool

fibs = 1 : 2 : (zipWith (+) fibs (tail fibs))

to neatly specify the Fibonacci sequence?" Alas, I suspect it does--after all, it doesn't have a finite length. There may be a way around that in this case, but the parametrization means you have to do calculations when type checking: "is this list empty?" "is this list at least as long as that list?" and in the general case, it can be undecidable.

All that said, it's still pretty darned impressive what can be done, and there's an excellent video from David Raymond Christiansen that shows it being done in the language Idris (which recently had a new version released). Check it out.

It's in the papers...

I don't necessarily agree with the New York Times's editorial position, but I applaud this fellow.
"Haskell in the Newsroom" is worth your time.

IMVU and Haskell

Check out "What it's Like to Use Haskell" over on the IMVU Engineering Blog. The author is enthusiastic, with enthusiasm based on, but also tempered by, the results of real-world use. Thank you, Mr. Friesen.

Another thing someone's probably done before

In part it comes from still being fairly new to the language, but I think it also has to do with the seriously cool high level of abstraction Haskell lends itself to. When I come up with a useful function, it's already been implemented, and far more generally. The first time around I wrote a special case of sequence. Here's my next function, and I am sure I just have to come up with the right monadic way to think of it to get the right string to hand to Hoogle to find the official version.

maybeIterate :: (a -> Maybe a) -> a -> [a]
maybeIterate f x = x : maybeIterate' (f x)
  where maybeIterate' Nothing   = []
        maybeIterate' (Just fx) = fx : maybeIterate' (f fx)

So in the background I will keep an eye out for blog posts like the one that clued me in to sequence, and I will continue the study of monads.

Numerical methods and C++

I sometimes have to use C++. What I'm curious about is numerical methods code for C++. Functions for finding zeroes or fixed points of functions, or numerical integration or differentiation, typically take as a parameter (a pointer to) the function being integrated, differentiated, or whatever. Here comes the fun part: what if that function is a member function?

You can't pass a pointer to a member function to something wanting a plain old pointer to function. The C++ FAQ section 33.2 answers the question "How do I pass a pointer-to-member-function to a signal handler, X event callback, system call that starts a thread/task, etc?" with the one word "Don't", followed by a hack to give the effect involving a function that returns the result of calling the member function via a global pointer to an object of the appropriate class that you have to assign to. Note the use of the "G word" there--so you can kiss the hack goodbye if multiple threads try it.

Googling "numerical methods in C++" turns up a fair number of book titles. I hesitate to buy a C++ book for the same reason I scrupulously avoid buying Microsoft products, but I may have to make an exception.

NixOS

If you are a Linux user, you probably have heard at least a bit about the package manager/format hoohah. .rpm, .deb, etc.... It's been blamed for dissuading companies from porting their apps to Linux, and one must admit that it's a pain when a package isn't available using the package manager your distribution uses (*cough*Java*cough*).

It's only recently that I've found out about the Nix package manager. I'm somewhat embarrassed, because it's based on a lazy, pure functional language. You can have multiple versions of a package installed at the same time, eliminating what in the Windows world is called "DLL hell". It supports multiple users, so non-privileged users can install packages. (And if two users happen to install the same package, there will only be one instance of it.)

The Nix package manager is the basis for the Linux distribution NixOS. It may not be for everybody; there are programs out there that still have wired-in assumptions about where things live that NixOS doesn't necessarily follow. OTOH, it may well be for me, judging by a blog post titled "How I Develop with Nix". (Do also check out the comments from where a link was posted on Reddit.) In particular, check out the section titled "Profiling", which talks about getting around the very problem that made me swear off cabal, i.e. not having profiling versions of packages.

Think I'll start out with it on a virtual machine, and once I get some more experience, maybe it will become my distro of choice. (Admittedly I have yet to migrate to Sabayon as I've been making noises about...)

Getting fed up...

Sigh. First streams, now matrices. The Ubuntu repositories don't have all the Haskell packages one would like. I guess it's going to come to removing the Ubuntu packages and building ghc and associated code myself or via.... [cue ominous music]... cabal.

I may combine that with a switch to Sabayon Linux, which will let me get back into Gentoo with minimal hassle, and also return to Enlightenment.

Once I can get the packages I want, it should be pretty easy to do Project Euler #11; the main trick for me is reading the input, and I think something like this should do the trick:

import qualified Data.Matrix as M
import qualified Data.ByteString.Char8 as B

parseLine :: B.ByteString -> [Int]
parseLine s = case B.parseInt s of
                Nothing      -> []
                Just (n, s') -> n : parseLine s'

main = do
         s <- getContents
         let m = M.fromLists $ map (parseLine . B.pack) (lines s)

and then it's just how to generate the contiguous portions of the matrix so I can find which one has the maximum product.

Someone has done this before...

Happened across it while looking for Haskell-related blogs to link to: from Derek Wyatt's blog, a post titled "Haskell sequence over functions -- explained" gives this example:

sequence [(> 4), (< 10), odd] 7

gives the result we would get from our wonkyMap with the same arguments:

[True, True, True]

What is this sequence?

:t sequence
sequence :: Monad m => [m a] -> m [a]

If m is a monad, then sequence will map a [m a] to a m [a].

So, our functions must have the type m a for some m and a. Let's inquire with ghci once again.

:t sequence [(> 4), (< 10), odd]
sequence [(> 4), (< 10), odd] :: Integral a => a -> [Bool]

Aha. That reminds me of the exercise in Typeclassopedia: show that (->) r is a functor.  (WTF? That's the type of functions that take an argument of type r.) Turns out that not only is it a functor, it's a monad; Typeclassopedia refers to it as the "reader monad".

Unfortunately, we have dueling definitions here. Mr. Wyatt's blog post uses the Learn You a Haskell for Great Good! definition, as opposed to what we see from ghci, though clearly the end result is the same. Let's carry on with the version that ghci uses, and you neither want nor need to see my long, rambling, not-yet-successful attempt to convince myself that (->) r really is a monad. We'll just suppose it is for now, and see what the consequences are.

Haskell's type inference takes the most general possible type. The functions we have here are all Integral a => a -> Bool.  In fact, I suspect Integral b => b is r for this example, using b because the sequence declaration has already used a, and I don't need help getting more confused. For us, the a mentioned in sequence's declaration is Bool, and in our case, we have

sequence :: Integral b => [b -> Bool] -> b -> [Bool]

which is exactly what we wanted for our solution to the problem in that online Haskell course.

So, bye-bye wonkyMap; we can rewrite strong as

strong password = all $ sequence constraints password
    where constraints = [(>= 15) . length,
                         any isUpper,
                         any isLower,
                         any isNumber]

...and I definitely need to get my head around monads, if for no other reason than to come up with alternative Hoogle search strings. :)

Only a few years late...

So I was looking at a list of requirements, trying to make sense of them and make sure I would come up with something that would satisfy them, when it occurred to me that you should bypass the middleman. Requirements should be written in the form of tests that pass iff the software satisfies them.

Turns out I'm a few years late. Dan North has already come up with the notion of "behavior-driven development". There are now parsers for "ubiquitous languages" in which one writes descriptions of how a piece of software should behave if you do certain things, and from those descriptions one generates tests. One example: Cucumber. (Why Cucumber? I think it's from the convention of having output from tests that fail come out in red while tests that pass generate output in green--like a cucumber.)

Haskell has some very neat facilities for testing, e.g. QuickCheck and HUnit. Is there a behavior-driven development package for Haskell? Why yes, there is: HSpec. Sweet.

You'd think someone's done this before

Saw a link to an Introduction to Haskell course set up by a student at the University of Virginia (which inspired others), and saw one of the homework assignments: write a function

strong :: String -> Bool

which returns True iff the string passed to it is a "strong" password, which for purposes of the assignment means

• it's at least 15 characters long
• it includes upper case letters
• it includes lower case letters
• it includes digits

Easy enough to write, especially if you import Data.Char, but the lesson is about higher-order functions, and you're urged to use the things taught in the lesson (among which are the character classification functions in Data.Char). So it occurred to me that it would be good to write the function in a way that makes it easy to add other constraints--perhaps it shouldn't be in some list of bad password choices, say.

So, we'd like to write it to take a list of String -> Bool, apply them to the String, and confirm that they all return True. The lesson teaches about all and about  map, but we don't want map here.

map :: (a -> b) -> [a] -> [b]

but we want

wonkyMap :: [a -> b] -> a -> [b]

An obvious implementation is

wonkyMap [] _      = []

wonkyMap (f:fs) x = (f x) : (wonkyMap fs x)

though I'm sure it can be written more elegantly as a fold.

(UPDATE: Duh... list comprehensions are your friend.

wonkyMap fs x = [f x | f <- fs]

just as you could define

map f xs = [f x | x <- xs]

There's a pleasing symmetry there.)

Then you have

strong password = all $ wonkyMap constraints password

    where constraints = [(> 15) . length,

                          any isUpper,

                          any isLower,

                          any isNumber]

OTOH, someone has to have come up with wonkyMap before--and given it a better name. Hoogling the signature, though, didn't turn up anything. Does it look familiar to anyone?

UPDATE: Oops... make that (>= 15) . length

Project Euler #14 revisited

Recently I came a cross a link from... aargh, was it Reddit or Hacker News? I should have saved the link. In any case, it was from someone who had decided to try to optimize a solution to Project Euler #14 in Haskell, and it induced me to revisit the program I'd written for it. Euler #14 asks the musical question, which number less than a million takes the longest for its Collatz sequence to get to 1?

{-# LANGUAGE BangPatterns #-}

import Data.List
import Data.Ord

collatzChainLen :: Int -> Int
collatzChainLen n = collatzChainLen' n 1
    where collatzChainLen' n !l
            | n == 1    = l
            | otherwise = collatzChainLen' (collatz n) (l + 1)
          collatz n = if even n then n `div` 2 else 3 * n + 1

pairMap :: (a -> b) -> [a] -> [(a, b)]
pairMap f xs = [(x, f x) | x <- xs]

main :: IO ()
main = print $ fst (maximumBy (comparing snd) (pairMap collatzChainLen [1..999999]))

Compiled with ghc -O2, about 3.5 seconds on my 2.8 GHz Athlon 64.
Of course, the thing to do is memoize collatzChainLen; if you have a sequence {x₁, x₂, ..., x_n = 1, ...} you know the values for all the x_i, not just x₁.

But the first thing I did was say to myself, "Self, that collatzChainLen function is yet another example of grinding out yet another boilerplate recursive function rather than using the Prelude and getting in the groove of using what the language gives you, combining functisns as Mr. Hughes and Mr. Backus recommend. Surely we can be more Haskelly than that." (I swear I hadn't been watching Looney Tunes.)

The result:

collatzChainLen :: Int -> Int
collatzChainLen n = 1 + (length . takeWhile (/= 1) . collatzChain) n
    where collatzChain = iterate collatz
          collatz n = if even n then n `div` 2 else 3 * n + 1

We've made the acquaintance of iterate before; it generates as much of the endless chain as we wish to grab, which we wish to do until we find a 1 with takewhile (/= 1). length counts how many that is, and while it doesn't affect the comparison, we add 1 rather than fib about what we're calculating. One add shouldn't affect much (or should it? If we're not memoizing, we're doing it rather a lot).

Does that affect the time? It sure does, but not as we'd like. With that version of collatzChainLen, the program takes about 9.7 seconds to run, almost three times as long! What's the difference? "It's profilin' time", as Ben Grimm would say if he learned Haskell. The evidence is pretty clear:

first version:   total alloc = 192,049,472 bytes
second version: total alloc = 21,277,532,952 bytes

Over a hundred times as much RAM, and who knows how much GC? I could believe that's the issue.

The lower-level version of collatzChainLen generates and counts values one at a time... but then, Haskell is lazy, so the higher-level one should, too--but perhaps something about creating a list of those values, or rather a couple, since takeWhile can't destructively chop off the list, but has to copy, is increasing the overhead. I recall that there's supposed to be something called "stream fusion" that gets rid of that overhead; let's look up how that's done.

Oops. haskell.org seems to be down, but apt-get install is my friend. Looks like there's a Data.Stream library, and you create streams instead of lists. Here goes:

import qualified Data.Stream as S

collatzChainLen :: Int -> Int
collatzChainLen n = 1 + (length . S.takeWhile (/= 1) . collatzChain) n
    where collatzChain = S.iterate collatz
          collatz n = if even n then n `div` 2 else 3 * n + 1

Time? Drat; still around 9.7 seconds. Doesn't seem very fused to me with that time, and a blog post about this very problem had example code that had

import qualified Data.List.Stream as S

and used S.length instead of just length. OK, rummage through synaptic for the libghc packages (sorry, but I've been burned by cabal), and...??? Darned if I can find it, and the blog post dates from August 2013, so it's not all that old.

While we're trying to figure out how to get to Data.List.Stream in Ubuntu, we can do one other thing:

collatz n = case n `divMod` 2 of
                (n', 0) -> n'
                _       -> 3 * n + 1

turns out to make a very respectable difference, cutting runtime down by 1.2 seconds! More news as it happens.

UPDATE: clearly

colllatzChainLen (2 * n) == 1 + collatzChainLen n

so, given a number < 500000, you could double it until it was in [500000, 999999] and have a number with a longer chain, so we can start with 500000.

Results: 4.6 seconds for the "Haskelly" version, 1.2 seconds for the explicitly written out collatzChainLen.

UPDATE: -fllvm is your friend; it took the runtimes down to 3.5 seconds for the higher-order function flavor... and 0.29 seconds for the version with the written-out collatzChainLen.

Project Euler #78

If you have a pile of n coins, there's some number of distinct ways to divide them up into piles; let's refer to the function that maps a number of coins to the number of ways to divide them up into piles as p. The example they give is for five coins; as you can see, p 5 == 7.

OOOOO
OOOO O
OOO OO
OOO O O
OO OO O
OO O O O
O O O O O

Things to note:

1. the coins are all alike; it's not like you have different denominations that could distinguish one pile from another with the same number of coins.
2. order doesn't matter; OOO OO and OO OOO don't count as different.

One's first inclination is to say "If you create a pile of m coins, that leaves n - m, and there are p (n - m) ways to divide them up, right?" Not exactly; if you don't confine yourself to piles of m or fewer coins, you'll violate (2). So we need what we'll call p', where p' n m is the number of ways to divide up n coins into piles of size m or less. Then you have

p :: Int -> Int
p n = sum $ map (p' n) [1..n]

That leaves p'. There are some obvious base cases:

p :: Int -> Int -> Int
p' n m
    | m >= n - 1 = 1
    | m == 1     = 1

If not those, then it turns out to be

    | otherwise = sum [p' (n - m) i | i <- [1..m']]
                      where m' = min (n - m) m

That matches up with examples I did by hand, and in fact corrected some I managed to get wrong, so I am fairly sure it's correct.

Unfortunately, the full Euler Project problem asks for

head [n | n <- [1..], (p n) `mod` 1000000 == 0]

and this way of calculating p' and thence p is up there with the canonical recursive Fibonacci series function for slow. We need a faster way to evaluate  p', or even better a fast way to calculate p without having to bother with p'.

So we generated some lines of output, the i^th line being map (p' i) [1..i],  Here it is, left justified:

[1]
[1,1]
[1,1,1]
[1,2,1,1]
[1,2,2,1,1]
[1,3,3,2,1,1]
[1,3,4,3,2,1,1]
[1,4,5,5,3,2,1,1]
[1,4,7,6,5,3,2,1,1]
[1,5,8,9,7,5,3,2,1,1]
[1,5,10,11,10,7,5,3,2,1,1]
[1,6,12,15,13,11,7,5,3,2,1,1]
[1,6,14,18,18,14,11,7,5,3,2,1,1]
[1,7,16,23,23,20,15,11,7,5,3,2,1,1]
[1,7,19,27,30,26,21,15,11,7,5,3,2,1,1]
[1,8,21,34,37,35,28,22,15,11,7,5,3,2,1,1]
[1,8,24,39,47,44,38,29,22,15,11,7,5,3,2,1,1]
[1,9,27,47,57,58,49,40,30,22,15,11,7,5,3,2,1,1]
[1,9,30,54,70,71,65,52,41,30,22,15,11,7,5,3,2,1,1]
[1,10,33,64,84,90,82,70,54,42,30,22,15,11,7,5,3,2,1,1]
[1,10,37,72,101,110,105,89,73,55,42,30,22,15,11,7,5,3,2,1,1]
[1,11,40,84,119,136,131,116,94,75,56,42,30,22,15,11,7,5,3,2,1,1]

And here it is right justified:

[1]
[1,1]
[1,1,1]
[1,2,1,1]
[1,2,2,1,1]
[1,3,3,2,1,1]
[1,3,4,3,2,1,1]
[1,4,5,5,3,2,1,1]
[1,4,7,6,5,3,2,1,1]
[1,5,8,9,7,5,3,2,1,1]
[1,5,10,11,10,7,5,3,2,1,1]
[1,6,12,15,13,11,7,5,3,2,1,1]
[1,6,14,18,18,14,11,7,5,3,2,1,1]
[1,7,16,23,23,20,15,11,7,5,3,2,1,1]
[1,7,19,27,30,26,21,15,11,7,5,3,2,1,1]
[1,8,21,34,37,35,28,22,15,11,7,5,3,2,1,1]
[1,8,24,39,47,44,38,29,22,15,11,7,5,3,2,1,1]
[1,9,27,47,57,58,49,40,30,22,15,11,7,5,3,2,1,1]
[1,9,30,54,70,71,65,52,41,30,22,15,11,7,5,3,2,1,1]
[1,10,33,64,84,90,82,70,54,42,30,22,15,11,7,5,3,2,1,1]
[1,10,37,72,101,110,105,89,73,55,42,30,22,15,11,7,5,3,2,1,1]
[1,11,40,84,119,136,131,116,94,75,56,42,30,22,15,11,7,5,3,2,1,1]

We see some patterns.

• p' n 2 == n `div` 2 for n > 1. There are n `div` 2 possible piles of size two, and each is either there or split into two piles of size one, so that makes sense.
• Check out that trailing portion. There's a steadily growing hunk of the end that never changes. It's because the maximum possible pile of what's left is no longer constrained, so past a certain value of n, p' n (n - m) = p m.

Let's check that last one by looking at the initial hunk of map p [1..]:

[1,2,3,5,7,11,15,22,30,42,56,77,101,135,176,231,297,385,490,627]

So, it would seem that p n == p' (2 * n) (2 * n - n) == p' (2 * n) n. That gets rid of the sum at the p end, but adds to the work on the p' side; it's not any faster, and maybe even a little slower. Worse yet, this isn't like factorial, where the number of trailing zeroes never shrinks--we're liable to end up using Integer instead of Int for the values the problem wants, which will make things slower still.

What we're really doing is counting "integer partitions". The integer being partitioned is the total number of coins, and each pile is represented by the number of coins in it (all the coins are alike). I will have to study this further.

UPDATE: maybe we can generalize the argument we did for p' n 2. What can we say about p' n k? There are n `div` k possible piles of size k, all right, but the leftovers have more than one possible state, and it's not as easy to enumerate them so that the pile sizes stay monotonically non-increasing.

Data Structures

The conventional wisdom on data structures in functional languages is that immutability costs you. It was a big advance when Okasaki came up with algorithms that have an amortized complexity for lazy functional languages the same as the non-amortized complexity for the data structures with destructive update.

Work is going on in that area; check out Edward Kmett's post "Part I: Deamortized ST" about how to come up with a way to make the amortization unnecessary. (It's a bit reminiscent of Asimov's essay "Behind the Teacher's Back".) I'm hoping that will ultimately mean Haskell will be usable in still more situations.

(Actually, I should just say "read Edward Kmett's blog" period. It's all very good stuff.)

Let ghc help you write your code--must {read, watch} from Matthew Brecknell

Suppose you had no idea what function composition should do, but you were given the type:

(.) :: (b -> c) -> (a -> b) -> a -> c

and had to define it:

(f . g) x = ???

We want a c. The only way we know to get a c is with f, which will give us a c if we give it a b. The only way we're given to get a b is with g, which will give us a b if we give it an a. Hey, x is an a!

(f . g) x = f (g x)

is thus the only way to honor the declaration.

Yes, this is a simple example... but not only can the method be applied more generally, you can get ghc to help you with error messages of the form "I expected a here but {you gave me something else, I can't infer that}". For details and a video working through examples, check out Matthew Brecknell's blog post "Hole-Driven Haskell".