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".

A new problem

OK, let's try a new problem this time: Project Euler #92. I happened across a YouTube video of a "live coding" session by Daniel Silverstone. It's worth your while. I admire the courage of someone willing to show such an unedited session, typos and all, and besides, the video shows just how expressive Haskell is.

So, having seen Mr. Silverstone's video, this is another case in which I can't pretend to having had all the insights myself. I will write the code, though; it's not copy and paste.

The problem has to do with a function that maps n to the sum of the squares of the digits of n expressed in base ten. I can't quite bring myself to use show, so...

sumSqDigits :: Int -> Int
sumSqDigits n = sumSqDigits' n 0
    where sumSqDigits' n !s
           | n < 10    = n * n + s
           | otherwise = let (q, r) = divMod n 10 in sumSqDigits' q (r * r + s)

Now, it can be shown that for non-zero n, iterate sumSqDigits n will "end" with one of

• repeat 1
• concat $ repeat [89, 145, 42, 20, 4, 16, 37, 58]

So, Project Euler problem #92 is the musical question "how many of the positive integers less than ten million are of the second kind?"

We'll borrow Mr. Silverstone's function name...

terminator :: Int -> Int
terminator n
    | n == 1 || n == 89 = n
    | otherwise         = terminator $ sumSqDigits n

Well, then just type

length . filter (== 89) $ map terminator [1...9999999]

and you're done, right? Well, you are if you're willing to wait.  It is, after all, evaluating sumSqDigits rather a lot. Mr. Silverstone uses the State monad to memoize terminator and speed things up quite a bit.

How best to do this? There are only two values, so one could argue that we're wasting space saving terminator values; we can recast the function as

terminatesWith89 :: Int -> Bool
terminatesWith89 n = n == 89 || terminatesWith89 (sumSqDigits n)

which in turn suggests a bitmap, since we have a bound on the values we're looking at. Ten million bits isn't all that much--at least not these days. But we can do better than that. After all, sumSqDigits takes an n-digit number to at most 81 * n. We're just interested in numbers of up to seven digits, so we need only memoize up to 7 * 81 = 567. Seventy-one bytes beats 1.2 megabytes any day.

For that matter, we can decide on how much space we're willing to devote to it, and just apply sumSqDigits until it's in range for our bitmap, which won't take too long. Say we start with n = a googol - 1. One hundred digits, so sumSqDigits n is at most 8100. That's four digits, so sumSqDgits of that is at most 243, and we're good to go.

With that it didn't take long to write up a 31-line (including "declarations" and blank lines) Haskell program that, when compiled with ghci -O2, solved the problem in a little over three seconds. We wrote a little helper function that followed the problem terminology:

chain :: Int -> [Int]
chain n = iterate sumSqDigits n

It proved nice to use along with dropWhile... but could we make it better? Well, we might have been able to speed up the generation of our bitmap, because everything on the same chain has the same terminator value, so we could have filled in a whole chain at once.

(Come to think of it, a chain that ends with 1, 1, 1, ... has to get there by way of a power of ten. Other values can get there, e.g. 86 (64 + 36 = 100);  is there some way they can be characterized? Maybe there's a closed form solution.)

You call that a cleanup?

OK, over on Reddit they make a good point, echoing Simon Thompson in the intro to the second edition of Haskell: The Craft of Functional Programming. He mentions there emphasizing the higher-order functions of the Standard Prelude because (I'm paraphrasing here, and working from memory) otherwise students tend to stay at the lower level, churning out the same recursive schema over and over again... and I'm certainly People's Exhibit #1 of that with bdigits and lastLE. I mean, good grief.

So, let's review the main idea: given two Integer values n and b, where b > 1, we want to know how many digits it takes to write n in base b. The thing is that n may be big, so rather than the stock counting digits one at a time, comparing against b each time, we compare with the values of iterate square b, i.e. [b ^ (2 ^ i) | i <- [0..]].  Only a finite number will be less than or equal to n. If none are, then obviously one base b digit will do. Otherwise, divide by the last one, which will have the form b ^ (2 ^ j) for some j, add 2 ^ j  to a running total of digits required, and do it again.

What drove us down the road we took is that one's first thought,

takeWhile (<= n)

has to copy a prefix of the list of powers of b, and we really only care about the last element of that prefix.

So perhaps what we need is a foldWhile. That would let us prime the pump with what we want for the base case, and walks the list rather than chopping off a prefix.

More when I have time

lastLE cleanup

Remember lastLE? Back when we were trying to determine how many digits/bits it takes to represent a non-negative Integer, we used it to avoid copying a prefix of a list just to grab the value at its end.

lastLE :: Integer -> [Integer] -> Maybe (Integer, Int)

lastLE n xs =
    let lastLE' xs prevVal prevIndex
           | head xs <= n = lastLE' (tail xs) (head xs) (prevIndex + 1)
           | otherwise    = if prevIndex < 0 then Nothing
                                             else Just (prevVal, prevIndex)
    in lastLE' xs (-1) (-1)

It still bugs me. First, the mixed guards and if/else are less than idiomatic Haskell:

lastLE n xs =
    let lastLE' xs prevVal prevIndex
           | head xs <= n  = lastLE' (tail xs) (head xs) (prevIndex + 1)
           | prevIndex < 0 = Nothing
           | otherwise     = Just (prevVal, prevIndex)
    in lastLE' xs (-1) (-1)

Second, the head and tail are clumsy:

lastLE n xs =
    let lastLE' (x:xs) prevVal prevIndex
           | x <= n        = lastLE' xs x (prevIndex + 1)
           | prevIndex < 0 = Nothing
           | otherwise     = Just (prevVal, prevIndex)
    in lastLE' xs (-1) (-1)

Third, that accumulator should be strict, so pretend there's a ! in front of prevIndex in that let.

I think I'd rather see Nothing as the otherwise case.

lastLE n xs =
    let lastLE' (x:xs) prevVal !prevIndex
           | x <= n         = lastLE' xs x (prevIndex + 1)
           | prevIndex >= 0 = Just (prevVal, prevIndex)
           | otherwise      = Nothing
    in lastLE' xs (-1) (-1)

Come to think of it, we can be consistent in our use of pairs:

lastLE n xs =
    let lastLE' (x:xs) !prev@(val, index)
           | x <= n         = lastLE' xs (x, index + 1)
           | prevIndex >= 0 = Just prev
           | otherwise      = Nothing
    in lastLE' xs (-1, -1)

Now, one more thing. Can I get rid of the magic numbers there, the -1s? The -1 as initial index needs to be there to get the correct position in the list of powers of the base. The initial val will never be used. If the whole pair is strict, though, we couldn't pass in ⊥ ("bottom", the undefined value, which is undefined in Haskell). Perhaps

lastLE n xs =
    let lastLE' (x:xs) prev@(val, !index)
           | x <= n     = lastLE' xs (x, index + 1)
           | index >= 0 = Just prev
           | otherwise  = Nothing
    in lastLE' xs (undefined, -1)

would be the way to go. Hey! We never even use that value here--it's used in the caller--so we can write

lastLE n xs =
    let lastLE' (x:xs) prev@(_, !index)
           | x <= n     = lastLE' xs (x, index + 1)
           | index >= 0 = Just prev
           | otherwise  = Nothing
    in lastLE' xs (undefined, -1)

and it should work (sure enough, it does), as well as expressing explicitly that we don't use it. Well, that's one magic number down, at least.

Small things...

Good grief... I've let this go for a month. (UPDATE: No, I haven't. It helps to read the dates carefully.)

An option in pattern matching that I wasn't fully aware of: you can use _ when you don't care about what is there, and have no need to refer to it. I guess I did use it on the ByteString readInteger that returns Maybe(Integer, ByteString):

process line = case map B.readInteger (B.words line) of
               [Just (m, _), Just (n, _)] -> numXs m n
               _                          -> -1

But I didn't realize that I could use it in

ysInRange m n d
    | nVal == n = nPos - mPos + 1
    | otherwise = nPos - mPos
    where (_,    mPos) = findFirst' m
          (nVal, nPos) = findFirst' n
          findFirst' x = case findFirst (Ge x, Any) (dDigitYTree d) of
                             Just (Max i, Max j) -> (i, j)
                             Nothing             -> (tenToThe d, numNDigitYs d)

because we need only check for an exact match at the high end.

Look for elegance

Check out nomeata's blog post "On taking the last n elements of a list". There's the obvious

takeLast n = (reverse . take n . reverse)

but all that reversing is a lot of work. The temptation is to think imperative and haul out the heavy artillery... but then he realized how it could be done efficiently and idiomatically. I won't put a spoiler here; read the original.

"Evelyn was a modified dog..."

So, how can we modify the tree to do what we'd like?

As the code stands, it's set up to say that if a is in the type class Semilattice, you can create a SearchTree of as that you can efficiently search. What we want to say is that if you have some type a that has a function key :: a -> b where b is in the type class Semilattice, you can create a SearchTree of as that you can efficiently search based on the result of key. Then, if we have a list of (palindrome, position) pairs, our key function is just fst (the function that gives you the first item in a pair). Now, how to express that in Haskell?

UPDATE: come to think of it, wouldn't that subsume the existing cases? For them,

let key = id
 
UPDATE: OK, I should have remembered. You can create relationships between type classes. (Take a look at this diagram of relationships among types and type classes in the Standard Prelude.) So, it looks like I can say something like this:

{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}

class Semilattice a where
    meet :: a -> a -> a

class Semilattice k => Keyed r k | r -> k where
   key :: r -> k

The first type class definition is straight from Mr. van Laarhoven's code; a type is a Semilattice type if it has a meet function which satisfies the semilattice requirements for meet.

Keyed is a "multiparameter type class"; it sets up a relationship between a type r (intended to suggest "record") and a type k (intended to suggest "key"). "r -> k" is a "functional dependency". Those two features aren't in the Haskell 98 standard, but GHC supports them (if they're enabled in a pragma) and I suspect they're in Haskell'.

(Historical aside: you might not realize that the preceding paragraphs use terms that originate with Algol 68: "standard prelude" and "pragma". Algol 68 is a language that suffered from a lot of undeserved bad press; you really should read C.H. Lindsey's excellent paper on its history, and check out Algol 68 Genie. For those of you interested in functional programming, I should note that partial parametrization was up for addition to the standard and is implemented in Algol 68 Genie.)

So, I will proceed on this basis. It should be a simple rewrite of the Semilattice tree code, though I wish there were a way to have one version of it to do everything; code reuse is a Good Thing. The above lines make it past ghci, so I hope it's a good start. (Then, of course, the question is whether it really does help the Fair and Square program! OTOH, in a way I don't care, because the real goal was to learn more Haskell, not to mention that I was just lucky that the original blog post gave an example in which the Semilattice type was a pair.)

What next?

Yesterday it occurred to me that the profiling output might be more accurate if I just asked for -p, rather than always asking for the heap tracking. The first part of the profile output:

    Mon Jul  1 22:23 2013 Time and Allocation Profiling Report  (Final)

       ultimatepalindrome20 +RTS -p -RTS

    total time  =        0.17 secs   (167 ticks @ 1000 us, 1 processor)
    total alloc = 116,700,576 bytes  (excludes profiling overheads)

COST CENTRE            MODULE                %time %alloc

fromList.(...)         SemilatticeSearchTree  19.2   26.8
ysPair.noTwos          Main                   11.4    8.9
floorSqrt.floorSqrt'.y Main                   10.2    5.6
choices                Main                    7.8   10.7
process                Main                    7.2   11.0
fromList               SemilatticeSearchTree   6.6    0.6
ysPair.spread.(...)    Main                    4.8    2.4
ysPair.noTwos'         Main                    4.2    4.2
meet                   SemilatticeSearchTree   3.6    0.0
satisfy                SemilatticeSearchTree   3.0    0.0
main                   Main                    1.8    1.0
numYs                  Main                    1.8    0.1
ysPair.pairSum         Main                    1.8    3.5
makeYTree              Main                    1.8    4.1
bDigits                Main                    1.8    0.0
meet                   SemilatticeSearchTree   1.8    3.5
mkBranch               SemilatticeSearchTree   1.8    3.2
ysPair.noTwosChoices   Main                    1.2    0.7
bDigits.bDigits'       Main                    1.2    1.9
bound                  SemilatticeSearchTree   1.2    0.0
ysPair.noTwos'.base    Main                    0.6    1.2
ysPair.spread          Main                    0.6    2.5
ysPair.twoNPlus1List   Main                    0.0    2.9

The cost of ysPair is spread out over several lines; add them up and you get 32.4% of the time (choices count; there's only one caller) and 37% of the allocation... but then, the rest of the profile does the adding for us; let's check it out. OK, I was pretty close. fromList with its descendants takes up about another third of the time and allocation.

floorSqrt and ceilSqrt should be low-hanging fruit. As we've said before, it converges quadratically, doubling the number of valid bits each time around... but for numbers of up to 330 bits, that's as many as nine iterations starting from just one good bit as floorSqrt does. So, the temptation is to use the double precision square root directly where possible, and use it as a first approximation otherwise. With 52 good bits as a starting point, three iterations will do the job for values up to a googol.

On the other hand, coming up with some way to use the semilattice search trees to just search on one value (and only store the maximum value searched for as the meet!), with the rest of the node along for the ride, would save time and memory and be more generally applicable, and would help me with Haskell in general. That's the way to go.

Dueling with the input

Still thinking about how to keep those combinations (or palindromes, now that we see we can avoid some sorting if we go ahead and generate them) around only as long as necessary. One thing that's kind of tempting is to take advantage of what we know about when we can reuse these values--we can generate them for one of  2 * k and 2  * k + 1 and then reuse them (with some modification) for the other. So, why not generate the trees two at a time?

"Well, you don't know that you'll need both of them," you object. "The Code Jam people could make you waste your time by generating input that only asks about ranges corresponding to Ys with even numbers of digits (or odd, take your pick)." And you're right, they could. Our attempts to speed things up can be subverted... but I think this one is worth a try.

At least initially, I'd be inclined to generate a list of pairs of trees, with the first being for 2 * k digits, the second being for 2 * k + 1. (Yes, one-digit Ys are a special case again.) More news as it happens.

UPDATE: new best time output:

real    0m0.192s
user    0m0.160s
sys     0m0.028s

Total allocation is down to not quite 117 MB, compared with 124 MB before, and "maximum residency" is around 9.6 MB. GC time is down to 41.2% of total execution time. I guess that's better than almost half, but that still seems high. One thing that is gratifying: the first time I used -sstderr to save GC info, it listed almost 160 MB copied during GC and over 16 MB maximum residency. Now the maximum copying during GC is just 63 MB, less than half what it was before.

And that's using spread, which is expensive because it does Integer divides. We'll do that next, and then I really should try to ditch the gratuitous meets of values in the trees that we never search for. That should be pure gravy.

UPDATE: spread is worth it to avoid having to sort again. (Before we weren't taking advantage of that, so spread was overhead.)

Improvements

Chaddaï Fouché kindly responded to a query I put out on the haskell-beginners mailing list, suggesting:

• using iterate
• using Data.Vector or MemoTrie instead of lists where possible

What's iterate, you ask?

iterate :: (a -> a) -> a -> [a]

Give it a function f and a starting value s and it will hand you back

[s, f s, f . f s, f . f . f s, ...]

So, for example, rather than

powersOfTen = map (10 ^) [0..]

we can write

powersOfTen = iterate (10 *) 1

and instead of

bigTwos = map (2 ^) powersOfTwo

we can write

square n = n * n

bigTwos = iterate square 2

Not shabby, eh?

About those Vectors: they're a data structure that makes for O(1) (i.e. constant time) indexing, as opposed to the O(n) time for lists. There are two flavors: Data.Vector and Data.Vector.Unboxed. The unboxed version has lower overhead, but can't be used on all types.

So I added

import qualified Data.Vector as V

and changed a little bit of code:

          memoPair   = V.generate halfN (\i -> tenToThe i + tenToThe (n - 1 - i))
          pair i     = memoPair V.! i
          pairSum v = V.foldl' (+) shell (V.map pair v)

choices :: Int -> Int-> [V.Vector Int]

m `choices` n
    | n == 0           = [V.empty]
    | m == n           = [V.enumFromStepN m (-1) m]
    | otherwise        = [m `V.cons` c | c <- (m - 1) `choices` (n - 1)]
                         ++ ((m - 1) `choices` n)

and indeed it helped. I hope it will help more soon; I can't have an unboxed vector of Integers (or is that a vector of unboxed Integers?), so if I want the choices Ints unboxed, I can't use the unboxed vector map; I'll have to roll my own function for that.

A bigger payoff, for now, had to do with sorting, or the minimizing thereof. It's easy to make oneTwos come out in order; noTwos is the gotcha. We did the folliowing:

nDigitYs n = (merge (oneTwos n) (noTwos n)) ++ twoTwos n
    where twoTwos n
              | even n    = [twoShell]
              | otherwise = [twoShell, twoShell + tenToThe halfN]
              where twoShell = 2 * shell
          oneTwos n
              | even n    = []
              | otherwise = map (+ (shell + 2 * tenToThe halfN))
                                (0 : map pair [halfN - 1,halfN - 2..1])
          noTwos n
              | even n    = base
              | otherwise = merge base [p + tenToThe halfN | p <- base]
              where base  = sort $ map pairSum (noTwosChoices !! (halfN - 1))

and poof! The profiler output shows total time as 174 ticks, 0.17 seconds.  The best time output is now

real    0m0.214s
user    0m0.176s
sys     0m0.028s

No reduction in the percentage of time used for garbage collection, though.

You know... now it may be worth keeping palindromes rather than combinations, because we would only have to sort them once instead of twice--and that might cut the memory usage as well.