In the great tradition of nibbling at the edges...

...and because I am still puzzling over why pairsum is grabbing so much RAM, let's contemplate floorSqrt. We're just calling it 2,000 times, and yet it's collectively taking up over six percent of the CPU time and five percent of allocation?!

Let's remind ourselves of floorSqrt (tweaked because I forgot that I'd memoized powers of two)

floorSqrt :: Integer -> Integer

floorSqrt 0 = 0
floorSqrt n = floorSqrt' (twoToThe ((1 + bitsIn n) `div` 2))
    where floorSqrt' x =
              let y = (x + n `div` x) `div` 2
              in if y >= x then x else floorSqrt' y

This is a special case of the Newton-Raphson method that was known to the ancient Greeks. We treat zero specially to avoid one of the gotchas of Newton-Raphson, but aside from that, this particular example is well-behaved, and has what they call quadratic convergence: once you're sufficiently close, each successive guess is good to twice as many digits as the one before. Our first guess is good to one bit, so the next is good to two bits, then four, then eight, etc. We're dealing with Xs in ranges up to a googol, or around 2³³⁰, which would me no more than nine iterations per value. Checking the profiling output, we see a total of 7,175 + 7,472 = 14,647 calls to floorSqrt', which is within our constraint, but darned close to it, so the inputs are tending to the high end of the range.

I kind of hate to cheat, but let's see what happens with

floorSqrt n = floorSqrt' $ (floor. sqrt . fromIntegral) n

[pause to compile, run, check profile output, and check... uh-oh.]

What happens is that it doesn't give the correct results; the output is not correct, which is weird. You'd think that would give you 52 good bits, and thus just need three iterations, but then, look at that termination condition. There has to be something about it that needs that first trial value to make it work, so you can't just plug in another, even better, guess. Taking the first value from the large input, the floorSqrt we've been using returns

10000010000100100100100001000000

while using the supposedly better first guess gives us

10000010000100100017503961350144

We could change the function to take ceil log2 (# bits in n / 52) iterations, but I think I should get back to dealing with memory usage and pairsum. (Even with that change you probably would end up flipping a coin for whether you got the floor or the ceiling of the square root of n.)

Now that we have better data...

...let's see what's happening.

First, the single sample point faked us out, making us think we were using less RAM than we really were. We are chewing up pretty nearly ten megabytes--still less than that C program, but a bit disappointing.

Or maybe more disappointing than we thought. Here's the output describing heap usage from a run with the options +RTS -sstderr:

     202,392,800 bytes allocated in the heap
      84,601,736 bytes copied during GC
      12,550,176 bytes maximum residency (8 sample(s))
         220,824 bytes maximum slop
              30 MB total memory in use (0 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0       387 colls,     0 par    0.10s    0.10s     0.0003s    0.0009s
  Gen  1         8 colls,     0 par    0.06s    0.06s     0.0080s    0.0230s

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    0.19s  (  0.20s elapsed)
  GC      time    0.16s  (  0.16s elapsed)
  RP      time    0.00s  (  0.00s elapsed)
  PROF    time    0.00s  (  0.00s elapsed)
  EXIT    time    0.00s  (  0.00s elapsed)
  Total   time    0.36s  (  0.36s elapsed)

  %GC     time      45.2%  (45.3% elapsed)

  Alloc rate    1,028,021,469 bytes per MUT second

  Productivity  54.8% of total user, 55.0% of total elapsed

Whoa. Nearly half our time spent garbage collecting? There's got to be a way to improve that. (And on the graph, nearly two megabytes taken by pairsum?)

Learning to avoid cabal--or not

Well... last night I thought I'd install EclipseFP, a package for Eclipse to support Haskell development. When you fire it up, it goes looking for packages it wants, and apparently uses cabal to install them. It did; I watched it do so for some time.

It turned out to be a waste of time; when I fired up Eclipse (which I'm rather new to) and clicked on the little lambda over to the left, a window opened up that looked half-drawn and very broken. I suspect that was issues with Eclipse--perhaps I should wipe the latest version that I grabbed and installed, and settle for the ancient version that, for some reason, is what Ubuntu has in its repositories.

OK, so I'll pass on an IDE for Haskell for now, or start up with leksah.

This morning, I had one of those sudden realizations that you get that make you laugh at yourself. Why do those memory usage graphs look like pyramids? Because the default sample interval is 0.1 seconds, and I have the run time down around 0.2 seconds, sort of like taking a sample every couple of minutes and expecting to get an accurate playback of a song.

A peek at Real World Haskell and I see the option to change the interval, so I recompile and run... and it claims that I should link with an option to turn on the RTS capabilities.

What?!

Another compile or two to make sure I did indeed specify the right options... and then some Googling, because the same thing still happened.

From the Google results, I suspect that all those cabal installs pulled in versions of libraries that don't support profiling, and that ghc is pulling them in. Following people's advice, I rm -rf ~/.ghc. No luck, still can't profile.

I have learned one thing: I am not going to let cabal touch my computer ever again. I would have sworn that I specified that packages should just be installed for me, so the rm should have done the trick. I guess I can look for Haskell libraries dated yesterday and delete every single one of them.

UPDATE: BZZT! Turns out the issue comes from trying to use that -i option to override the heap measurement interval. That's what's giving me the problem. I was wrong... and maybe I'll consider using cabal... sometime when I know Haskell much better than I do now.

UPDATE: Found it... since ghc 6.x, you have to compile with -rtsopts=all to be able to use some (OK, most) RTS options, lest the logging they permit be used to breach security. I will have to try to figure out how -i would give you the opportunity to do so over and above the output that -h*, which doesn't require -rtsopts=all, allows. (OK, maybe it would be possible to allocate a lot of RAM or not every [interval], a lot of RAM means a 1, not much means 0.)

hlint

Back in the early days of Unix, when the PDP--11/70's main advantage over the later 8/16-bit 6809 was having separate I/D (instruction and data) space, so that you could have 64K of code with access to 64K of data, the virtue of the Unix Way of small programs that did one thing and did it well was a necessity. One of the things it gave rise to was a separate program, "lint", to check C source code for constructs that might be evidence of a coding error, so that the compiler could concentrate on simply generating code.

Nowadays, C compilers often do some of the checking that was once delegated to lint (though separate lint programs still exist, and are very good--e.g. splint, or Gimpel Software's excellent products).

What about Haskell? Well, for Haskell there's hlint. It will give you advice on the Haskell source that you feed it.

jejones@eeyore:~/src/haskell_play$ hlint ultimatepalindrome14.hs
ultimatepalindrome14.hs:236:13: Warning: Use zipWith
Found:
  map showsResult $ zip [1 ..] (map process (tail $ B.lines s))
Why not:
  zipWith (curry showsResult) [1 ..] (map process (tail $ B.lines s))

ultimatepalindrome14.hs:237:5: Error: Use .
Found:
  mapM putStr $ map ($ "\n") r
Why not:
  mapM (putStr . ($ "\n")) r

2 suggestions
jejones@eeyore:~/src/haskell_play$

So there are (at least!) a couple of places I could have written arguably better, easier to read, more idiomatic Haskell. Pretty cool.

What's with that "."? In that context, "." is the function composition operator; given two functions, one of which returns values of the type the other takes as input, composing them gives you a function that applies first one function, then the other. Huh? Better to write it in Haskell:

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

f . g  x = f (g x)

or, equivalently,

f . g = \x -> f (g x)

\ isn't a character escape; it's the closest the creators of Haskell could come to  the Greek letter lambda, as in Alonzo Church's "lambda calculus". A lambda expression is an "anonymous function"; the way to read the right hand side of that last line is "the function that, given an argument x, returns f (g x)".

(Ironically, the lambda was in turn a sort of best approximation; it was inspired by the use of the circumflex in the notation of Russell and Whitehead's Principia Mathematica. Details here.)

Anyway... there's also a Haskell "style scanner" that one can either just get suggestions from or use in the fashion of indent. I will have to check it out.

It's not just an idiom...

In this program I've taken a list of values and fed it to

zip [1..]

You'll recall that zip takes two lists and returns a new list as long as the shortest of the lists handed to it. Each element of the new list is a pair of values at corresponding positions in the lists, so that, for example,

zip [1..] "hiya" == [(1, 'h'), (2, 'i'), (3, 'y'), (4, 'a')]

Remind you of anything you learned to do as a kid? Yes, it's constructing an bijection from a set (here represented as a list) to the first however many counting numbers, aka counting. Little did I know back then that I was preparing for Haskell.

Further restructuring

It occurred to me that I could push the code even further, making the connection between counting and generation more apparent. Recall the counting code:

numNDigitYs 1 = 3
numNDigitYs n = numTwoTwos n + numOneTwos n + numNoTwos n
    where numTwoTwos n = if even n then 1 else 2
          numOneTwos n = if even n then 0 else n `div` 2
          numNoTwos  n = if even n then s else 2 * s
                         where h = n `div` 2 - 1
                               s = sum [h `choose` i | i <- [0..min h 3]]

Here's the new generation (which sounds like some cheesy 60s thing--sorry!); we've gotten rid of justOnes.

nDigitYs 1 = [1,2,3]
nDigitYs n = sort (noTwos n ++ oneTwos n ++ twoTwos n)
    where halfN  = n `div` 2
          pair i = tenToThe i + tenToThe (n - (i + 1))
          twoTwos n
              | even n    = [twoTwosNoOnes]
              | otherwise = [twoTwosNoOnes, twoTwosNoOnes + tenToThe halfN]
              where twoTwosNoOnes = 2 * tenToThe (n - 1) + 2
          oneTwos n
              | even n    = []
              | otherwise =  map (+ common)
                                 (0 : [pair i | i <- [1..halfN - 1]])
              where common = pair 0 + 2 * tenToThe halfN
          noTwos n
              | even n    = base
              | otherwise = concat [[p, p + tenToThe halfN] | p <- base]
              where pairsum xs = foldl' (+) (pair 0) (map pair xs)
                    base       = map pairsum (noTwosChoices !! (halfN - 1))

noTwosChoices = [concat [n `choices` k | k <- [0..min 3 n]] | n <- [0..]]

We explicitly generate the oneTwos values in a way that makes clear that there are n `div` 2 of them, and it's similarly clear that numNoTwos is correct. Life is beautiful, right?

Well, not quite. We're now eating up ten megabytes of memory instead of six and change; total allocation is up from 127 MB to 140 MB, and total time, according to the profiling output, is up by .03 seconds. Saving lists of lists of Ints will definitely add to the allocation, so we would expect that. We're really beating on pair and pairsum, so that's probably where to look to get resource usage pared (no word play intended) back down.

UPDATE: We didn't memoize pair; doing that took us back down to 0.2 seconds, cutting the difference back down to .01 seconds. Total allocation is back down to 130 MB, but memory is up somewhere between 10 and 11 MB.

UPDATE: I wonder whether I forgot something, perhaps -fllvm. Recompiling and rerunning shows 0.19 seconds run time in profiler output, and just about 8 MB of memory usage (as opposed to total allocation, still around 130 MB).

I hate to do something as crass as packing three values in an Int or maybe an Int32, but it's grating to keep a list of lists of Ints hanging around, too. I'll give it a try. (The results were slower and used more RAM. So much for that.)

Maybe we're saving the wrong thing

Maybe we should keep the lists of combinations around instead of the base values; that would get us out of the divMod business, though OTOH it means hanging on to those lists; as it stands, I bet that they're being garbage collected.

The problem is that the bitmap moves us back to O(n), where n is the number of digits, rather than a fixed amount of work (because there's a fixed upper bound on the number of nonzero digits). We could roll our own bit fields, and stuff up to three bitsIn (n `div` 2) values into an Int or an Integer (actually, Haskell has Int8, Int16, Int32, and Int64 types to choose from; 21 bits a shot would allow for up to two million digits, times two for both halves of the palindrome--that's a lot bigger than the second official data set needs), but that seems kind of cheesy.

I'll have to think about that some more. In the meantime, I restructured the code a little bit to make the code that generates Ys more clearly parallel to the code that just counts them. The counting code:

numNDigitYs :: Int -> Int

numNDigitYs 1 = 3
numNDigitYs n = numTwoTwos n + numOneTwos n + numNoTwos n
    where numTwoTwos n = if even n then 1 else 2
          numOneTwos n = if even n then 0 else n `div` 2
          numNoTwos  n = if even n then s else 2 * s
                         where h = n `div` 2 - 1
                               s = sum [h `choose` i | i <- [0..min h 3]]

The generating code:

nDigitYs :: Int -> [Integer]

nDigitYs 1 = [1,2,3]
nDigitYs n = sort (noTwos n ++ oneTwos n ++ twoTwos n)
    where twoTwos n
              | even n    = [twoTwosNoOnes]
              | otherwise = [twoTwosNoOnes,
                             twoTwosNoOnes + tenToThe (n `div` 2)]
              where twoTwosNoOnes = 2 * tenToThe (n - 1) + 2
          oneTwos n
              | even n    = []
              | otherwise =  map (+ 2 * tenToThe halfN)
                                 (justOnes n (min 1 (halfN - 1)))
              where halfN = n `div` 2
          noTwos n
              | even n    = base
              | otherwise = concat [[p, p + tenToThe halfN] | p <- base]
              where halfN = n `div` 2
                    base = justOnes n (min 3 (halfN - 1))

Yeah, I guess that shows I'm seriously at the point of diminishing returns.

Nomenclature, and data structures

First, I should correct my terminology. All computers are binary these days; Setun was a one-shot, and the days of BCD-based computers are long gone.

This is the case of noTwos, and we've already dealt with the one case where a 3 can appear, so all the digits are either zero or one. So there's an obvious correspondence between these decimal palindromes and binary palindromes, pairing palindromes whose printable representations in their respective bases are equal.

Of course, n-digit binary numbers map exactly to the elements of the power set of a set of n elements, so we can represent our generated combinations (vide the function choices) as binary numbers. This is pretty tempting as an alternative to a [Int], but I don't want to iterate over all those bits to see which are set. After all, if we're just doing half of the palindrome, and we know the most significant digit has to be 1, that leaves at most three bits on.

Maybe we can use a trick: consider a number n of a type in the Bits class. What can we say about

n .&. (n - 1)?

(The Haskell bitwise "and" and "or" are a little funky-looking.)

Well, if n == 0, the result is zero; otherwise, n has some bits set, and one of them is the least-significant bit. n - 1 will have a zero in that position, and the more significant bits will be unchanged. The "and" gets rid of the trailing ones in case the least significant bit of n isn't in the ones position.

So, you can write an lsb function:

import Data.Bits

lsb :: (Bits a) => a -> a

lsb n = n `xor` (n .&. (n - 1))

and then peel the bits off one at a time. To do that, we shouldn't just give back the LSB, because we're generating the next value in the sequence while we're at it, so

bitStrip :: (Bits a) => a -> (a, a)

bitStrip n = (lsb, leftovers)
    where  leftovers = n .&. (n - 1)
           lsb       = n `xor` leftovers

and we can pull them off until we end up with zero...but we still have to map over to base 10; that lsb is 2^i for some i. I'm not sure whether the LSB trick will help any.

Of course there's a pattern

I can't believe I didn't see it before. Here's the deal:

• You can get from the base values for 2 * k digits to those for 2 * k + 1 digits as we described earlier.
• You can get from the base values for 2 * k digits to some of those for 2 * k + 2 digits.

If you take all the base values for 2 * k digits and make a two-digit gap between their halves, you have all the base values for 2 * k + 2 digits that don't have ones in the gap. The ones that do have ones in the gap have, outside the gap, the halves of 2 * k digit palindromes that have one less one per half than the base values for 2 * k digits, because only that way do you get ones left over to fill the gap.

The 65536-dollar question, then, is whether it's worthwhile to generate them this way. (And of course, you can generate the ones with one less one per half for 2 * k digits from those with two less ones per half for 2 * k - 2 digits. That doesn't go on forever, to be sure; since the added ones per half is no bigger than three, you're down to the base case of 100...001 pretty quickly. Again, though, is it worth it?)

One thing to consider that might help out: doing justOnes and noTwos sorts of things purely in binary and moving back to decimal only at the very end. If we take this route, it would also be worth dealing purely with the upper halves of the no-one Ys until it's time to convert. Consider: dealing with the Ys takes you down from at most 100 decimal digits to at most 50 decimal digits. Dealing with just the upper half takes you down to at most 25 decimal digits. Doing the stuff where the decimal digits are restricted to ones and zeros means you can operate on values with at most 25 bits until the very end. Save the arbitrary precision for where you need it, right? Even if, like me, you kind of hate setting arbitrary limits, in practice it will mean that your Integers will actually, for purposes of this problem, be small enough that operations on them will be, if they're implemented by a list or array of fixed-size hunks, the fastest kind, because there will only be one hunk to mess with.

(OK, I admit I've limited some things in my code to Int rather than Integer. OTOH, I think the (29-bit) Int values suffice, since there are only something like 51,000 of these palindromes up to a googol, for heaven's sake, I suspect, though I haven't calculated, that there aren't 2^29, or around 512 million of the palindromes, until a pretty darned big upper bound. OK, maybe it's 256 million--signed, right?--but that's still on up there.)

I'll have to think about this some more.

UPDATE: I pulled the counting routines, tweaked them to take and return Integer, and then did

last $ takeWhile (\(x, y) -> y <= 2^29) (zip [1..] nDigitYsPartialSums)

Turns out you have to get past 10^515 before the number of Ys risks going past 2^29 (and hence needing an Integer to count them), and thus the corresponding Xs would be up around 10^1030.

Not sure whether it helped

With that change, we now have the following for noTwos and associated functions. (Oh, yeah... we made the one-digit case special as we did for counting the d-digit Ys, for the same reason, i.e. putting the special case in one place.)

evenNoTwos = [justOnes n (min 3 (n `div` 2 - 1)) | n <- [2,4..]]

noTwos :: Int -> [Integer]

noTwos n
    | even n     = base
    | otherwise  = concat [[p, p + tenToThe halfN] | p <- map spread base]
    where halfN    = n `div` 2
          base     = evenNoTwos !! (halfN - 1)
          spread x = let (a, b) = x `divMod` (tenToThe halfN)
                     in a * tenToThe (halfN + 1) + b

Did it make a difference? Yes, but... it was a bit slower. By profiling output, 193 ticks instead of 187; time output looked worse, with the best apparently

real    0m0.232s
user    0m0.200s
sys     0m0.028s

Recall that for the previous version, the best was

real    0m0.224s
user    0m0.176s
sys     0m0.044s

On the other hand, looking deeper into the profiler output, before, noTwos and things it called was eating 19.3% of the CPU time, while now that's 7.8%, but to be fair we should include evenNoTwos, which snarfs 7.3%--but still that's a total of 15.1%, which should be a win. Despite total allocation being down, actual memory usage is up, from maybe 6.25 MB to close to 7 MB; we are, after all, keeping more stuff around.

Meanwhile, I'm still thinking...

I took a quick look at choices, and I don't see an obvious way to generate them in what would correspond to ascending order. (Though that doesn't mean there isn't one.)

One thing we can do, though... take a look at noTwos, the source of by far the most of the d-digit Ys for all but the smallest d values.

noTwos n
    | n == 1     = [1,3]
    | even n     = base
    | otherwise  = concat [[p, p + tenToThe halfN] | p <- base]
    where halfN = n `div` 2
          base = justOnes n (min 3 (halfN - 1))

In particular, look at how base is calculated. Suppose n == 2 * k for some k. Then for n and n + 1, halfN - 1 has the same value and thus min 3 (halfN - 1) has the same value. Sure, n is different, but that just means for n + 1 there's a one digit gap left to either put a one in or leave at zero. Check it out.

*Main> justOnes 2 0
[11]
*Main> justOnes 3 0
[101]
*Main> justOnes 4 1
[1001,1111]
*Main> justOnes 5 1
[10001,11011]
*Main> justOnes 6 2
[100001,101101,110011,111111]
*Main> justOnes 7 2
[1000001,1010101,1100011,1110111]

So, corresponding to a base value x for noTwos (2 * k), the corresponding value for noTwos (2 * k + 1) is

let (a, b) = x `divMod` (tenToThe k)
in a * tenToThe (k + 1) + b

Less work than going through justOnes, I'd say... and it applies whether we figure out a clever way to generate the values in order or not. Note also that this mapping preserves numerical order, so if we do find a way to cleverly generate these values in order, taking advantage of this mapping won't mess that up.

We can save the base values for even numbers, though it would be nice if we could somehow notice when we do their successors and render them garbage collectable; otherwise that's 15,275 Integers hanging around in lists. Haskell has some memoization types and functions around that go beyond the simple lists as arrays we've been doing, but I'm not sure whether they support something like that.

"You can see a lot by just looking" --Yogi Berra

There's a reason Code Jam urges one to look at examples--in this problem, it points you towards the constraint on the Ys that I didn't see because I only considered the most/least significant digit. And now, we'll see whether it can point us at how we might generate noTwos in ascending numerical order.

noTwos, like oneTwos, is based on the result of justOnes. For now, at least, let's ignore the case of an odd number of digits; for them, justOnes simply returns the same values as for the next lower number of digits with a zero in the middle. Now, to take a look at  

putStr $ unlines $ map show (sort (justOnes (2 * n) (min (n - 1) 3))

for some values of n.  Hey, wait; these are palindromes, so let's just look at the top half. After all, the rest is redundant, and it's the top half that determines the ordering.

4:

10
11

6:

100
101
110
111

8:

1000
1001
1010
1011
1100
1101
1110
1111
10:

10000
10001
10010
10011
10100
10101
10110
10111
11000
11001
11010
11011
11100
11101
11110

Oops. 10 breaks the pattern, because you can't have five ones. 12 shows a skip in the middle where there would otherwise be five ones, and as n increases, we'll presumably see more of those.

12:

100000
100001
100010
100011
100100
100101
100110
100111
101000
101001
101010
101011
101100
101101
101110
110000 <-
110001
110010
110011
110100
110101
110110
111000
111001
111010
111100

So, if we pretend these are binary numbers, up to eight digits (total), they go from 10...0 to 11..1. Past eight digits, they go from 10...0 to 11110..0, skipping any value that has more than four ones.

So there is a pattern. How to generate it efficiently?

UPDATE: I'm not sure there is a way; it's tempting to take the interval including all the numbers and then filter (\x -> popCount x <= o + 1) but that's a lot of numbers to generate and then mostly throw away; up at the high end for numbers up to a googol, that'd be around 28 million values that you'd grind out and then throw away all but a few thousand... and then you get to create the palindrome from the top half and switch back to base 10 on top of that. I don't think that's going to cut it... so, back to getting the list of lists of ones to come out in what corresponds to ascending order for the values justOnes generates.

Down with Gray Text!

It wasn't that long ago that having white text on a black background was le dernier cri. (Maybe it still is in some circles; a web page on the preannounced new Mac Pro was done up that way. After five seconds trying to read it on my phone, I gave it up as a lost cause.)

It seems to have been replaced lately by gray--in some cases pretty darned light gray--text on a white background.

Do web designers not care about whether anyone can read the pages they create?

Must-read about three monads

...and nary a burrito in sight. "Three Useful Monads"

Overkill (not really...)

So I thought that

ysInDigitRange d1 d2 = sum [numNDigitYs d | d <- [d1..d2]]

did some needless recalculation, and tried

memoNumNDigitYs = map numNDigitYs [1..]
nDigitYsPartialSums = scanl (+) 0 memoNumNDigitYs

ysInDigitRange d1 d2 = nDigitYsPartialSums !! d2 - nDigitYsPartialSums !! (d1 - 1)

scanl? Ah, that's like foldl, but instead of giving you the final result of the fold, it gives you a list: at the head is the "initial" value, followed by the results of the operation at each stage. nDigitYsPartialSums !! i is thus equal to

sum map numNDigitYs [1..i]

so that the two formulations of ysInDigitRange give the same result. Trying it, though, showed the scanl version taking a little more time, a little less total allocation--but actually allocating about nine megabytes! I didn't expect that. Just goes to show that not everything you think will be an optimization actually is.

UPDATE: OK, now I'm puzzled. The profiler output shows ysInDigitRange being invoked 325 times. Ditto for numNDigitYs, but... but... the only calls to numNDigitYs should be where values in memoNumNDigitYs are calculated, and we shouldn't ask for more than 50 of those, so what's with the other 275 invocations?

UPDATE: Comparing profiler output showed what I was overlooking. Before this change, there were a lot more calls to numNDigitYs. We're talking 2,732 + 325 instead of just 50 + 325, down by a factor of more than eight. (OK, that's where the 50 went, but where do the 325 come from?)

Moreover,  the time percentage of the top cost centers falls off a lot quicker with the change. In the past, that's tended to happen along with a decrease in execution time, but here it's gone up a bit... and why on earth does keeping around what should be just an extra 101 Integers cause such a jump in memory usage?! Time to generate some of those other memory graphs to see whether we can figure it out.

D'Oh! I forgot... I'd switched to -fllvm. LLVM is the abstract machine used as an intermediate form by clang, and if you have it installed, you can ask ghc to use it. That, or rather not doing that, made the difference in memory usage and time. We're back to a bit over 6 MB, and with profile output claiming .19 seconds. Best time output:

real    0m0.227s
user    0m0.184s
sys     0m0.040s

Worst time output:

real    0m0.233s
user    0m0.220s
sys     0m0.008s

I'm amazed that llvm makes that much difference. (Compilation takes a bit longer, but darn, it's worth it.)

"This time, fer shur!" --Bullwinkle Moose

OK. I really will go on to doing something else in Haskell; I'm thinking of something that will actually get me into the M-word.

(Though in the back of my mind I will still wonder whether there's a way to make justOnes generate values in ascending order other than just sorting the darn things...)

I want to give a little more data.

I went back and compiled and ran the C solution to the "Fair and Square" problem I'd grabbed, to get time output so I could remember something other than just the vague "about .4 seconds". The variance seems to be less than it is for "my" Haskell program, (I have to use the scare quotes given my the integer square root and semilattice search tree code) so I'll just give one output, as usual from the data set with values up to a googol, from the middle:

real    0m0.455s
user    0m0.432s
sys     0m0.020s

In comparison, the best value from a Haskell run:

real    0m0.224s
user    0m0.176s
sys     0m0.044s

If I had to guess, I'd say maybe the difference in sys time is due to memory allocation; the C program visits the heap for buffers for I/O for input and output files, but that's about it. The rest it does via fixed-size declared arrays.

How much memory, then? What dominates for the C program is a single array, a two-dimensional 50000 x 200 array of ints. That comes to 24 MB. The graph generated from profiling the current version of the Haskell program shows it using somewhere between 6.25 and 6.5 MB (I really wish they'd show the scale all the way up on that). So half the time, and given the other arrays in the C code, I think it's safe to say about a quarter of the RAM.

The C code definitely wins on size of compiled program, though; with both stripped, it's 860K for Haskell, a hair over 10K for the C program. To be sure, that's taking advantage of shared libraries. I wonder whether ghc pays much attention to only including what's necessary in a compiled program? (Or maybe it's my doing, and I should be more specific about what I import.)

Skating around the edges

Still thinking a bit about how to churn out those palindromes in order, or at least take advantage of what order there is. If you take a look at justOnes, it may be that each of the k values gives a hunk in ascending order, so you may be able to foldl' merge [] them and get them in order for cheaper... (UPDATE: no such luck) though that may not actually be cheaper than a single O(n log n) sort.

But now we take another peek at the output. The slowness of piecing together Strings to generate output has been noticed... oh, yes, it has been noticed... and people have come up with a way to cut down on it: using shows or something like it.

shows returns a function that takes a String as parameter and returns a String--the returned String is a printable representation of the thing (which must be in the ShowS type class) with the String parameter appended at the end. Instead of generating the String, you create a function that will generate it when called.

So.. remember

format :: (Int, Int) -> String

format (lineNum, result) = "Case #" ++ show lineNum ++ ": " ++ show result

? Well, now we have

showsResult :: (Int, Int) -> String -> String

showsResult (c, n) = ("Case #" ++) . shows c . (": " ++) . shows n

The use of "sections" makes the correspondence clear. Since we're now passing along functions instead of strings, they're used a little differently:

main = do
    s <- B.getContents
    let r = map format $ zip [1..] (map process (tail $ B.lines s))
    putStr (unlines r)

becomes

main = do
    s <- B.getContents
    let r = map showsResult $ zip [1..] (map process (tail $ B.lines s))
    mapM putStr $ map ($ "\n") r

Now, putStr, if you look it up, has the type String -> IO (), but it has a reassuring name. mapM, though, hints at the dreaded "M" word... monad. [insert dramatic chipmunk sound clip here]. All it is, though, is a monad-flavored map. Before, r was a [String]; now it's a [String -> String]. To get the [String], we use map with a section. (In Haskell, $ is the function application operator. We've used it before; its low priority makes it attractive as a way to avoid parentheses, but here, it's more than a syntactic convenience.) The "\n" means we don't have to use unlines, mapM putStr just runs putStr on each of the Strings on the list.

Does it help? Well... where time output is concerned, it's a wash. The profiler output shows it taking .01 second longer. Ah, but the memory graph! The -hc graph shows the format version somewhere around 7 MB, but the showsResult version at maybe 6.25 MB. That's got to count for something.

"Blessed rage for order, pale Ramon..."

So, the big cost center is nDigitYs. All it does is concatenate our three flavors of Y (the ones with no 2s, the ones with one 2, and the ones with two 2s) and sort them. Of the d-digit Ys, there are at most two with two 2s, and there are d `div` 2 with one 2. By far the majority are those with no 2s; for d == 49, of the 4,122 Ys, 4,122 - 2 - 24 = 4,096 have no 2s.

We know that the two-2 Ys are greater than all the others, so we could go with

nDigitYs n = sort (noTwos n ++ oneTwos n) ++ twoTwos n

...for what little good that does.

If you take a look at the result of oneTwos, you'll find that they're always in ascending order; we grind them out so that the optional 1 is first not there, and then in ascending order in the most signfiicant half, which suffices to order them. noTwos, on the other hand, isn't, for d > 7.  If there's some way to generate them in ascending order, we could go with

nDigitYs n = merge (noTwos n) (oneTwos n) ++ twoTwos n

(be sure to import Data.List.Utils) and O(n) beats O(n log n) any day. Unfortunately, I don't know of a good way to do that.

OK, I have to take one more look...

What can I say? I can sympathize with Lot's wife.

I have yet to show you the spiffy graphic you can get showing memory usage.

What does this tell us?

First, we can say that it looks like we're down about a megabyte from earlier. (Alas, they don't give the scale all the way up the vertical axis; if only they'd move the word "bytes" out of the way!) The other thing is, it looks like it's accumulating thunks until about halfway through when something finally demands evaluation and the RAM soufflé collapses.

(What's a thunk? The name goes back to the days when the world was new and people were figuring out how to generate code for Algol 60 programs. When the defining committee set up the rules for "call by name" they perhaps didn't realize the full consequences of their actions (vide Jensen's Device). Implementing it required generating a little function for each call by name parameter that returned a reference to it, and those little functions were called "thunks". In the context of Haskell and other lazy programming languages, a thunk represents a value that hasn't yet been evaluated.)

Real World Haskell has a chapter on optimization that we've leaned heavily on, and it has a section on memory usage,. thunks, and strictness issues. I will look more closely at that section and take one last shot at the palindrome program--after all, these issues will come up in the future as well. (No, really, this will be the last time...)

After the smoke clears...

OK. I think I've taken this problem as far as I can go... though someday I will set it up to read and write named files so I can run it in ghci. After all, this all started because I wanted to demonstrate that it's the algorithm, not the language or whether it's interpreted or compiled after reading Clif Reeder's "Ruby is Too Slow for Programming Competitions". On the other hand, I now have a Haskell program that's within spitting distance of running twice as fast as a C++ solution for the problem, taking a bit over 0.2 seconds to run the hairier official input file. Eight minutes, the time Code Jam gives you to run and send back results when it counts, is 480 seconds; I tend to doubt that interpreter overhead would come near a factor of 2,400, in which case the point is made.

This started with the less-than-noble motivation of "Well, I'll show him...", but quickly changed, for a couple of reasons.

• The unalloyed pleasure of problem solving always wins.
• Haskell is a sheer delight to use.

It's been quite a while since I've had this much fun programming. Haskell made it easy to clearly and simply say what I meant, and when it came time to speed it up, the tools were there, and even the optimized versions were still intelligible. Not to mention that while there were bugs, there weren't many and they were easy to find. I could pull pieces into ghci and test if need be, and need wasn't all that often. All this (aside from Twan van Laarhoven's semilattice search trees) done by a Haskell newbie and duffer pushing 60 (with 40+ years of indoctrination in imperative programming aside from a bit of Lisp here and there), in spare time on evenings and (mostly) weekends.

Here's the result. You'll need the semilattice code, which is available at the link above. I compiled with -O2 in the interests of speed.

{-
  OK, here's what I hope is the ultimate "Fair and Square" program.
  (See https://code.google.com/codejam/contest/2270488/dashboard#s=p2;
  the core of the problem is, given two positive integers, M and N, determine
  the number of base 10 palindromes between M and N inclusive that are squares
  of base 10 palindromes.)

  Best time output, reading the L2 input file, on a 2.8 GHz AMD Athlon II
  X4 630:

  real    0m0.227s
  user    0m0.208s
  sys     0m0.016s

  [Ultimate? Silly me; there's always something else to try. However, I think
  I'm done with it, for now at least.]

  The main insights to speed up the process follow.
  (4)-(6) are directly from the Google Code Jam analysis of the problem at
  https://code.google.com/codejam/contest/2270488/dashboard#s=a&a=2
  (9) is from Twan van Laarhoven's "Search trees without sorting", on his blog
  at http://twanvl.nl/blog/haskell/SemilatticeSearchTree,

  1. Look for palindromes in [ceil(sqrt(M))..floor(sqrt(N))] whose squares are
     palindromes. (We'll refer to one of the former set as X, and one of the
     latter set as Y, as is done in Google's analysis.)
  2. Don't look for palindromes, generate them. They're determined (aside from
     the middle digit if there are an odd number) completely by half their
     digits, so in a range of 10^(2*n) integers, there are on the order of
     10^n palindromes, a reduction on the order of (1).
  3. If Y^2 is to be a palindome, its most/least significant digit can't be
     greater than 3.
  4. In particular, the Ys are precisely those that, when you square them via
     long multiplication, have no carries in the addition of partial products.
  5. (3) implies that the Xs have an odd number of digits.
  6. (3) also implies that the middle digit of X, which is the sum of the
     squares of the digits of the palindrome being squared, must be no larger
     than 9. Therefore, *none* of the digits of Y can be greater than 3, and
     3 only appears in Y = 3.
  7. (6) implies that Ys come in three flavors: those with two 2s, those with
     one 2, and those with no 2s.
  8. Save for very small n, it's easier to count d-digit Ys than to generate
     them--and remember, the problem only asks for the count.
  9. By treating a binary tree as a sublattice, labeling branches with the
     upper bounds of the corresponding subtrees, one can quickly search for
     the first value in the tree >= a specified value.
 10. By labelling the sorted Ys for a given number of digits with their ordinal
     positions, one can use those tree searches to quickly determine how many
     Ys are in a given interval.

  NOTE: 0 is a one digit palindrome, and its square is, too, but since we are
  only dealing with positive integers, we'll ignore zero.
-}

import Prelude
import Data.List
import SemilatticeSearchTree
import qualified Data.ByteString.Char8 as B

-- It's easy to calculate the number of n-digit Ys.

numNDigitYs :: Int -> Int

numNDigitYs n = if n == 1 then 3 else numTwoTwos n + numOneTwos n + numNoTwos n
    where numTwoTwos d = if even d then 1 else 2
          numOneTwos d = if even d then 0 else d `div` 2
          numNoTwos  d = let d' = d `div` 2 - 1
                             n' = sum [d' `choose` i | i <- [0..min d' 3]]
                         in if even d then n' else 2 * n'
                                   
choose :: Int -> Int -> Int

m `choose` 0 = 1
m `choose` n = product [m - n + 1..m] `div` product [1..n]                            

-- We may wish to memoize numNDigitYs...

ysInDigitRange d1 d2 = sum [numNDigitYs d | d <- [d1..d2]]

{-
  Now the slow part: generating and counting Ys in a range other than "all
  d-digit numbers".
-}

powersOfTen = [10 ^ n | n <- [0..]] -- or, if you wish, map (10 ^) [0..]

tenToThe :: Int -> Integer

tenToThe n = powersOfTen !! n

twoTwos :: Int -> [Integer]

twoTwos n
    | n == 1    = []
    | even n    = [twoTwosNoOnes]
    | otherwise = [twoTwosNoOnes, twoTwosNoOnes + tenToThe (n `div` 2)]               
    where twoTwosNoOnes = 2 * tenToThe (n - 1) + 2

oneTwos :: Int -> [Integer]

oneTwos n
    | n == 1     = [2]
    | even n     = []
    | otherwise  = [p + 2 * tenToThe halfN
                      | p <- justOnes n (min 1 (halfN - 1))]
    where halfN = n `div` 2

noTwos :: Int -> [Integer]

noTwos n
    | n == 1     = [1,3]
    | even n     = base
    | otherwise  = concat [[p, p + tenToThe halfN] | p <- base]
    where halfN = n `div` 2
          base = justOnes n (min 3 (halfN - 1))

{-
  justOnes -- generate the palindromes with a specified number of digits all of
  whose non-zero digits are 1s. 1 must be the most/least significant digit; up
  to a specified number of 1s are scattered through the rest of each half.
  If the number of digits requested is odd, we leave a 0 in the middle digit,
  for the caller to either replace with another digit or not.
-}

justOnes :: Int -> Int -> [Integer]

justOnes n o =
    let halfN = n `div` 2
        shell = tenToThe (n - 1) + 1
        innards = concat [(halfN - 1) `choices` k | k <- [0..o]]
        -- originally pair i = tenToThe i + tenToThe (n - (i + 1))
        ascending  = take (halfN - 1) (drop 1 powersOfTen)
        descending = take (halfN - 1) (reverse (take (n - 1) powersOfTen))
        memoPair = zipWith (+) ascending descending
        pair i = memoPair !! (i - 1)
    in [foldl' (+) shell (map pair c) | c <- innards]

choices :: Int -> Int-> [[Int]]

m `choices` n
    | n == 0           = [[]]
    | m == n           = [[m, m-1..1]]
    | otherwise        = [m : c | c <- (m - 1) `choices` (n - 1)]
                         ++ ((m - 1) `choices` n)

nDigitYs :: Int -> [Integer]

nDigitYs n = sort (noTwos n ++ oneTwos n ++ twoTwos n)

-- We turn the sorted lists from nDigitYs into search trees.
-- Here we use van Laarhoven's semilattice search tree routines/types.
-- Sigh... we're doing a bit of cargo cult programming here, because even
-- though we'll never search on the ordinal index, we'll follow the blog
-- post example.

yTree :: Int -> SearchTree (Max Integer, Max Int)

yTree n = fromList [(Max x, Max y) | (x, y) <- zip (nDigitYs n) [0..]]

yTrees = map yTree [1..]

dDigitYTree :: Int -> SearchTree (Max Integer, Max Int)

dDigitYTree d = yTrees !! (d - 1)

-- ...leading up to this function.

ysInRange :: Integer -> Integer -> Int -> Int

ysInRange m n d
      | nVal == n = nPos - mPos + 1
      | otherwise = nPos - mPos
      where (mVal, 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)

-- counting decimal digits and bits
-- the numbers here are big enough to nudge us past the obvious

powersOfTwo = map (2^) [0..]

bigTwos = map (2^) powersOfTwo
bigTens = map (10^) powersOfTwo

bitsIn n   = bDigits n bigTwos
digitsIn n = bDigits n bigTens

bDigits :: Integer -> [Integer] -> Int

bDigits n xs = bDigits' n 1
    where bDigits' n s = let upToN = takeWhile (<= n) xs
                         in case upToN of
                            [] -> s
                            _  -> bDigits' (n `div` (last upToN))
                                           (s + powersOfTwo !! (length upToN - 1))

-- The following is derived from a response to a Stack Overflow question
-- http://stackoverflow.com/questions/1623375/writing-your-own-square-root-function
-- which in turn is taken from Crandall & Pomerance, "Prime Numbers: A
-- Computational Perspective".

floorSqrt :: Integer -> Integer

floorSqrt 0 = 0
floorSqrt n = floorSqrt' (2 ^ ((1 + bitsIn n) `div` 2))
    where floorSqrt' x =
               let y = (x + n `div` x) `div` 2
               in if y >= x then x else floorSqrt' y

-- then the obvious...

ceilSqrt :: Integer -> Integer

ceilSqrt n =
    let y = floorSqrt n
    in if y * y == n then y else y + 1

-- Finally, the top level, which counts where it can, and generates, filters,
-- and counts where it must. Note that we break the interval down into ranges
-- that are subsets of d-digit numbers for appropriate values of d.

numYs :: Integer -> Integer -> Int

numYs m n
    | mDigits == nDigits   = ysInRange m n mDigits
    | mPower10 && nPower10 = ysInDigitRange mDigits (nDigits - 1)
    | mPower10             = ysInDigitRange mDigits (nDigits - 1)
                             + ysInRange (tenToThe (nDigits - 1) + 1) n nDigits
    | otherwise            = ysInRange m (tenToThe mDigits - 1) mDigits
                             + ysInDigitRange (mDigits + 1) (nDigits - 1)
                             + ysInRange (tenToThe (nDigits - 1) + 1) n nDigits
    where mDigits    = digitsIn m
          nDigits    = digitsIn n
          mPower10 = m == tenToThe (mDigits - 1)
          nPower10 = n == tenToThe (nDigits - 1)

numXs :: Integer -> Integer -> Int

numXs m n = numYs (ceilSqrt m) (floorSqrt n)

-- A new main, heavily inspired by "Haskell I/O for Imperative Programmers"
-- http://www.haskell.org/haskellwiki/Haskell_IO_for_Imperative_Programmers
-- (What's with that tail, you ask? We're skipping that first line that says
-- how many more lines there are. They could have been persnickety and padded
-- the file with additional lines to make sure you have to pay attention to
-- that first line, but they weren't.

main = do
    s <- B.getContents
    let r = map format $ zip [1..] (map process (tail $ B.lines s))
    putStr (unlines r)

process :: B.ByteString -> Int

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

format :: (Int, Int) -> String

format (lineNum, result) = "Case #" ++ show lineNum ++ ": " ++ show result