Great Leap Forward

So, we made the change. Note that we are in fact generating the lists for the combinations!

{-
  Here's our new approach: rather than generate the upper half separately and
  then have to reverse it a digit at a time, we generate both halves at the
  same time, when we have information we need to do so.

  So, out with backwards and {odd,even}DigitsPal, and twoTwos goes back to its
  simple self.
-}

twoTwos :: Int -> [Integer]

twoTwos n
    | n == 1    = []
    | even n    = [twoTwosNoOnes]
    | otherwise = [twoTwosNoOnes, twoTwosNoOnes + 10 ^ (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)) 

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)

{-
  justOnes -- here's where the real action happens. Here we generate palindromes
  with a specified number of digits all of whose non-zero digits are 1s. 1 must
  thus be the most/least significant digit; in addition, between 0 and a
  specified number of 1s are scattered through the rest of each half. If the
  number of digits requested is odd, we will leave a 0 in the middle digit,
  which the caller may or may not replace with another digit.
-}

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]]
        pair i = tenToThe i + tenToThe (n - (i + 1))
    in [shell + sum (map pair c) | c <- innards]

This gives a definite improvement. time output again varies a fair amount; the best we've seen so far is

real    0m0.509s
user    0m0.488s
sys     0m0.016s

and the worst is

real    0m0.581s
user    0m0.564s
sys     0m0.008s

I wish I knew the source of the variation; perhaps other processes running at the time? If we can believe the best result, we're down to 1.27 times the run time of the C++ solution we tested... and we haven't touched the real time sink. Check this out:

        total alloc = 345,155,656 bytes  (excludes profiling overheads)

COST CENTRE            MODULE  %time %alloc

numWithin              Main     44.7   11.3
solve                  Main     24.7   55.5
digitsIn.digitsIn'     Main      6.1    7.8
nDigitYs               Main      5.6    7.8
tenToThe               Main      5.3    0.0
justOnes               Main      3.0    6.4
justOnes.pair          Main      2.2    3.1
floorSqrt.floorSqrt'.y Main      2.1    1.9
noTwos                 Main      1.0    1.0
choices                Main      1.0    2.9

So we've cut down total allocation by over 200MB (not to be confused with the memory usage, which takes a final spike of about a quarter megabyte above the previous version). We could take a look at avoiding actually generating the list of chosen exponents for the half-Ys, but we have to face facts. numWithin is now eating nearly half the time, and if the above is to be believed, solve, of all things, is responsible for over half the allocation! (Maybe I don't understand just what the profiling output is trying to tell me.) We're within spitting distance of the C++ version; a move from O(n) to O(log n) for the range check ought to bring us ahead.

P.S. Note that I've deviated from the define-before-use ordering that I've stuck to in the past, and the results compiled without complaint. I infer from this that Haskell must be waiting until it's read the whole file to do at least some of the type checking, so I should be able to get away with putting definitions in top-down order.

A little finesse...

It took a little bit to get this right:

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)

The gotcha is the n == 0 case; if you have it return [] the list comprehension gives the wrong result (i.e. an empty list!) with n == 1. (I didn't really need to make the m == n case go in descending order, but it looked nicer.)

Shouldn't take long to recast the palindrome program now; the idea, as alluded to earlier, is to represent half-Ys (aside from the easy-to-generate twoTwos) not by an Integer but by a list of the positions where we place 1s. Then backwards doesn't have to take the half-palindrome apart decimal digit by decimal digit, but can just reflect the positions and directly add up the corresponding powers of ten... but OTOH, do I really want to churn out a bunch of lists? For 50-digit Ys, that's up to 25 choose 3, or 2,300. If I go to representing them as a bit mask, I'm just chugging through bits rather than decimal digits--which would be faster, so perhaps that's the thing to do... or we generate both halves at the same time, while we have the information.

Going backwards

I figured I might try the same "loop unrolling" with backwards' as worked so well with digitsIn'. Admittedly I only did it once, but the results didn't go the right way:

COST CENTRE                MODULE  %time %alloc

numWithin                  Main     29.5    7.2
solve                      Main     21.0   35.3
backwards.backwards'       Main     16.6   17.5
backwards.backwards'.(...) Main      9.5   12.6
choices                    Main      5.2    8.9
digitsIn.digitsIn'         Main      5.2    4.9
nDigitYs                   Main      4.6    5.8
tenToThe                   Main      1.8    0.0
floorSqrt.floorSqrt'.y     Main      1.1    1.2
oddDigitsPals              Main      0.9    1.8
noTwos                     Main      0.7    1.3

As you can see by comparison with a previous profile, backwards is now taking up more time and allocation than before.

Guess we'll have to memoize...

backwardsMemo = [backwards' n 0 | n <- [0..]]
    where backwards' n accum
              | n < 10    = accum + n
              | otherwise = let (d, m) = n `divMod` 10
                            in backwards' d (10 * (m + accum))

backwards :: Integer -> Integer

backwards n = backwardsMemo `genericIndex` n

Um, no, we won't, not that way. In theory, at least, we should just be evaluating backwards a bit over 40,000 times, but then there are the calls to backwards'. After about two minutes I killed the process. (Does indexing cause all the intermediate items to be generated?)

We know that the tempting

backwards = read . reverse . show

loses big time in comparison with what we have already as well.

So, let's think about the particular values we're handing backwards:

• For twoTwos, 2 * tenToThe n. Here, backwards will grind down a zero at a time and give us back 2, every time.
• For oneTwo and noTwos,  it will always be a sum of between 1 and 2 (for oneTwo) or 1 and 3 (for noTwos) values of the form tenToThe i.

So... we blew it when we wrote our even and odd number of digit palindrome generators the way we did, and we don't want to generate them the way we are now. We don't want to represent our half-Ys as Integers, but instead as short lists of the powers of ten that you add to make them up. No Integers until we actually churn out the palindromes themselves. Then we can avoid a numeric backwards altogether.

"How low can you go?"

(Hey, with that title I can just keep my mouth shut and young people will think I'm quoting Ludacris, while old people will think I'm quoting Chubby Checker! Oops...)

So I decided to create a dummy version of the program that instead of calculating the number of "fair and square" numbers in a given range, just adds the low and high end of the range and prints it out, using the format the Code Jam problem requires:

numXs :: Integer -> Integer -> Integer

numXs a b = a + b

-- jagd's main, which I believe is set up to read stdin

solve i num = do
    [a, b] <- getLine >>= return . map read . words
    putStrLn $ "Case #" ++ show i ++ ": " ++ (show $ numXs a b)
    if i < num
       then solve (i+1) num
       else return ()

main = do
   num <- getLine >>= return.read
   solve 1 num

Compiled it up with -O2 and ran...

time ./dummy <C-large-practice-2.in >woof.dummy

real    0m0.106s
user    0m0.092s
sys     0m0.012s

Here the variance is a hefty factor; one run took 0.171 seconds!

The surprising part of the profiling was the memory allocation:

total alloc = 196,851,848 bytes  (excludes profiling overheads)

Compare that with the real program:

total alloc = 546,422,264 bytes  (excludes profiling overheads)

The full program is allocating lists of powers of ten, and potentially over 50,00040,000 Y values! There's something very inefficient going on in solve.

UPDATE: There must be a difference between "total alloc" and memory usage. Running dummy with the options to track memory usage shows memory usage peaking at maybe 60K, compared with the palindrome program, which peaks quickly at around 4M. Both programs must be allocating and freeing quite a lot, the difference being that we save up those lists of Ys and powers of ten.

Loop unrolling, Haskell style

We'd written digitsIn' to allow tail call optimization, so that the recursive call gets compiled as a loop:

digitsIn base n = digitsIn' n 1
    where digitsIn' n count
            | n < base  = count
            | otherwise = digitsIn' (n `div` base) (count + 1)

This does remind me of a blog post that Hacker News pointed at, in which another operation was being done in Haskell a digit at a time. The solution there, as here, is something kind of like loop unrolling:

digitsIn base n = digitsIn' n 1
    where base2 = base * base
          base4 = base2 * base2
          digitsIn' n count
            | n < base  = count
            | n > base4 = digitsIn' (n `div` base4) (count + 4)
            | n > base2 = digitsIn' (n `div` base2) (count + 2)
            | otherwise = digitsIn' (n `div` base)  (count + 1)

And the results:

real    0m0.679s
user    0m0.656s
sys     0m0.020s

From the profiler:

COST CENTRE                MODULE  %time %alloc

numWithin                  Main     31.6    7.0
solve                      Main     21.0   34.6
backwards.backwards'       Main     14.7   17.1
backwards.backwards'.(...) Main      9.0   12.9
choices                    Main      5.5   10.8
nDigitYs                   Main      4.8    5.7
digitsIn.digitsIn'         Main      4.5    4.8
oddDigitsPals              Main      1.3    1.8
noTwos                     Main      1.1    1.3
tenToThe                   Main      1.1    0.0
digitsIn                   Main      1.1    0.0
floorSqrt.floorSqrt'.y     Main      0.7    1.2

digitsIn' is down from almost 15% of the time to 4.5%, and from a sixth of the allocation to under a twentieth. We've gotten closer to the runtime of the C++ solution I grabbed and compiled, going from taking not quite twice as long to not quite 1.7 times as long. There ought to be a way to similarly speed up backwards'. (And what's up with solve?)

I'm really procrastinating on that data structure change, huh?

Guess we'll have to take the plunge

Just to check, I ran a version of the program with a revision of ysInRange, the function that actually does the range checking to do some of the alleged improvements mentioned earlier:

numWithin :: Integer -> Integer -> [Integer] -> Int

numWithin m n xs = length (takeWhile (<= n) $ dropWhile (< m) xs)

simpleYsInRange :: Integer -> Integer -> Int -> Int

simpleYsInRange m n digits = numWithin m n (ysByDigits !! (digits - 1))

-- Now for the real thing. As with the palindrome generation, the caller has
-- the number of digits of interest handy, so we take it as a parameter

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

ysInRange m n digits
    | digits < 8 = simpleYsInRange m n digits
    | mMSDs > 20 = 0
    | mMSDs > 11 = numWithin m n (twoTwos digits)
    | otherwise  = simpleYsInRange m n digits
    where mMSDs = m `div` tenToThe (digits - 2)

The results? Inconclusive; there's no clear difference. It's going to take a better data structure.

UPDATE:  here's the "cost center" listing from profiling output.

COST CENTRE                MODULE  %time %alloc

numWithin                  Main     28.2    6.1
solve                      Main     18.1   30.3
digitsIn.digitsIn'         Main     14.8   16.6
backwards.backwards'       Main     12.1   15.0
backwards.backwards'.(...) Main      9.8   11.4
nDigitYs                   Main      4.4    5.0
choices                    Main      4.3    9.4
tenToThe                   Main      1.5    0.0
floorSqrt.floorSqrt'.y     Main      1.3    1.0
digitsIn                   Main      1.3    0.0
noTwos                     Main      1.0    1.1
oddDigitsPals              Main      0.6    1.6

So it is indeed walking the lists of generated palindromes that's the time sink. (I do wonder about digitsIn and backwards', though; I guess arbitrary precision arithmetic will mean heap allocation, but that much?)

A tree? Maybe not...

So, how about that tree? It's tempting to do something like van Laarhoven describes, though I'd be tempted to put maxima on the left links, minima on the right links, to assist with the range check. But wait a minute. Now that we have added the purely combinatorial part, we never actually generate and selectively count more than one power of ten's worth of palindromes at a time, and at least some of those are of the trailing edge of an interval that covers a power of ten. Since the first d-digit Y is 1, for d == 1, or 10^(d-1) + 1, for d > 1, for them it's not like we have to skip a lot of Ys to find those of interest. That leaves the leading part, either [m..n] or [m..tenToThe mdigits - 1]. Just how many Ys are there for a given number of digits, anyway?

*Main> map numNDigitYs [1..50] [3,2,5,3,8,5,13,9,22,16,37,27,60,43,93,65,138,94,197,131,272,177,365,233,478,300,613,379,772,471,957,577,1170,698,1413,835,1688,989,1997,1161,2342,1352,2725,1563,3148,1795,3613,2049,4122,2326] 
*Main>

OK, up to 4122 is a bit long to grind through.

We can do some things:  for small d, might as well just scan as we are now. For bigger d, though, we have the one or two two-2 palindromes at the high end, and at the low end, 10.....01, optionally 10..020..01, and then no-2s and optional one-2s interleaved up to 11110...01111 (d even) or 11110..010...01111 (d odd) for d > 9. If d is 3 or 5, then the largest non-two-2 Y is 121 or 11211, else for d <=9 it's 111...111.

So, for sufficiently large d, if the first two digits of m are > 11, then there are at most two (the two-2 flavor) Ys in [m..n], depending on the parity of d (whether it's odd or even) and the value of n.

This will help when it lets us rule out no-2 and one-2 Ys, but when we have to dig into those, which account for the vast majority of the Ys with a given number of digits, we're back to scanning that list. Doing it right means coming up with a better data structure, one that lets us quickly determine how many values are in a given range.

Hmmm... come to think of it, can we get away with just generating the digits that suffice to determine each Y? Maybe those are easier to work with--it would certainly let us dump all the calls to backwards, which still takes a fair amount of time. The only issue will be the range checking--which we'd have to recast in terms of the half-palindromes. (UPDATE: Sigh, I guess not. You'd be doing it repeatedly rather than just once.)

Speaking of aging...

It's been some time since I've noticed my eyes don't accommodate as they used to--but if you take a look at this blog's layout and hit ctrl-0, that text is pretty darned small. Sorry about that; I've already changed the text color from what I think was the default for this style, #262626, to pure black, which helps, and bumped the text size a little bit. I'll tweak it some more; needing to hit  ctrl-+ a few times to read my own blog is pretty sad.

"I think that I shall never see..."

I'm rummaging around Stack Overflow and Hoogle. Data.Set and Data.Map are implemented underneath as balanced trees, but for my purposes I need the tree nature exposed so I can do something like

countInRange Empty m n = 0
countInRange (Node value (Tree lhs) (Tree rhs)) m n
    | value < m = countInRange rhs m n
    | value > n = countInRange lhs m n
    | otherwise = 1 + (countInRange lhs m n)
                    + (countInRange rhs m n)

The set fold doesn't guarantee ordering, so I don't see a way to let it quit early or otherwise skip portions that can't contribute to the count.

Maybe I will have to roll my own.

UPDATE:  Twar van Laarhoven's "Search trees without sorting" looks promising for my purpose. Dank U wel.

Holy [expletive]!

jejones@eeyore:~/src/haskell_play$ time ./ultimatepalindrome <C-large-practice-2.in >ultimate2.out

real    0m0.866s
user    0m0.852s
sys     0m0.012s

I want to try a thing or two more before posting the source, but that beats the heck out of 1.7 seconds, and is sneaking up on that 0.4 second time from one of the C/C++ solutions.

Things to take away from this:

• Even with Google's explanations of problems, it's worth looking further. I would like to think that the Code Jam folks did this on purpose; it's a good pedagogical technique. (Pamela McCorduck's Machines Who Think quotes Minsky lamenting that in the book he and Papert wrote on perceptrons, they pretty well did all the interesting stuff, leaving nothing for others to follow up on.)
• In conversation about Mr. Reeder's blog post that started all this hoohah, a friend pointed out that whatever one does in some scripting or high-level language could be rewritten in C/assembly language/VLIW microcode and still be faster. That's true, but... to quote Weird Al Yankovic, "You could even cut a tin can with it--but you wouldn't want to!" Nowadays most of the time we accept the overhead of using C or [shudder] C++ rather than assembly language. We're now down to within a factor of a bit overunder two, and that's by a duffer fairly new to Haskell, fiddling with code off and on during spare time over roughly three weeks. Mr. Reeder was right; he just needed to push his analysis a bit further, as did I.

UPDATE: Thinking about it, breaking the problem down as this version does has the advantage that when we are generating and then walking the list of Ys, we're confining ourselves to a single power of ten range rather than walking the list all the way from 1 to wherever.

Then, too, we're not generating the Xs, but operating purely in the realm of Ys, which I'm sure saves some heap.

Profiling output claims that the code is spending 27.4% of its time and 34.7% of its allocations in backwards', and that it is the biggest "cost center" of the program. We made a point of writing it to permit tail call optimization, but perhaps there's another way to improve it.

UPDATE: divMod (or quotRem, depending on the situation) to the rescue!

These return a tuple with the quotient and the remainder (I suspect either the one that C programmers expect or the one that mathematicians expect, respectively), since division frequently hands you the remainder/modulus for free, so a simple tweak of backwards into

backwards n =
    let backwards' n accum
          | n < 10    = accum + n
          | otherwise = let (d, m) = n `divMod` 10
                        in backwards' d (10 * (m + accum))
    in backwards' n 0

took the time output to

real    0m0.779s
user    0m0.756s
sys     0m0.020s

and backwards' down to 13.8% of the time and 15.3% of the allocations. The new major "cost center", as one would expect, is now ysInRange, the function that determines how many of the (generated) Ys are within a specified interval. Let's see whether it's as easy to improve as backwards was. [pause] Probably not, at least not by perusing Hoogle (the Haskell equivalent of Google, i.e. a place to look for information about Haskell's standard modules). Time to create a different data structure, I bet a tree.

UPDATE: Shucks. I remembered a blog post Hacker News linked to about some Haskell code in which the bottleneck was digit-at-a-time conversion, so I took a whack at a fancier digitsIn that went reduced by larger powers of the base first, then took things a digit at a time. It didn't make much difference that I could see.

Looks promising...

A first cut at code that takes the aforementioned approach looks promising. Merging the purely combinatoric code with the code that generates palindromes gave rise to a name clash or two, I left out a couple of needed parentheses, Haskell didn't like my referring to some type by the name of Integber, I wrote [1..3] rather than [1,3] for the list of one digit Ys with no 2s, and I forgot that if m is a power of 10, digitsIn 10 m is one more than the common log of m. That corrected, a few lines from the big input give the same result as the previous version, and the new version was noticeably faster on some of them, which is what we were after. Now to compile and try the whole megillah. There are some things I should probably memoize, but get it right first, right?

By the way, ghci (Glasgow Haskell Compiler Interpreter?!), the interpreter, has a useful, um... OK. "Spell checker" isn't quite right, but if you make a typo, it will notice that it's close to something else that is defined, and say "Maybe you meant..." So far, it has always suggested what I meant.

I am highly impressed with how little Haskell gets in the way, and how easy it is to say what I mean. Of course, I haven't seriously entered the world of do and monads, so I should make a point of pushing myself there, but I expect it won't be the bête noire it's made out to be. I promise that once I do, I will not write yet another "monads are like containers/tortillasburritos/etc." tutorial. (In that context, don't miss Daniel Spiewak's article on the subject.)

So after all that...

What would be our ultimate solution for the problem? Given the previous post, we have two cases:

• m' and n' both have d digits: in this case, if m' is greater than the largest d-digit Y, the result is 0; otherwise, we have to actually generate the d-digit Y values and return the number that are in range.
• m' has d digits and n' has d' digits, d < d': the result then is the sum of the number of Ys between m' and 10^d-1, the sum of the number of k-digit Ys for k from d + 1 to d' - 1, and the number of Ys between 10^d'+1 and n'.

We need smarter memoization of the function that maps d to the list of d-digit Ys and the function that maps d to the count of d-digit Ys, that doesn't grind out all the results for all values from 1 to d. (If they were recursive in d, I wouldn't mind.) Time to do some research.

UPDATE: or maybe I have a misconception about lazy evaluation. If the list elements don't depend on one another save for their position, maybe you don't have to evaluate an element just to skip over it with !!.

(Of course, we can easily generate the largest d-digit Y. If d == 1, it's 3; otherwise it's the largest two-2 d-digit Y.)

Combinatorics

You know, for the Code Jam problem,  you just have to know how many "fair and square" numbers there are in a given interval. That's not necessarily the same thing as having to generate all of them.

Having done (or read) the analysis, we know that the number of "fair and square" numbers between m and n is the number of palindromic numbers between ⌈sqrt(m)⌉ and  ⌊sqrt(n)⌋ (we've shown how to get those values in previous posts)--let's call those values m' and n' respectively--for which the sum of the squares of the digits is less than 10.

Let's just consider d-digit base 10 numbers for d > 1, so we don't have to worry about 0 or 3. If [m'..n'] includes all d-digit base 10 numbers, or, given our theorem, at least the d-digit base 10 numbers from 10...01 to either 20...02 (if d is even) or 20..010...02 (if d is odd), then it has all the d-digit palindromes of the sort we want. (Writing them that way is serious handwaving about how big d is, sorry.) How many is that?

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

numTwoTwos d
    | d == 1    = 0
    | even d    = 1
    | otherwise = 2

numOneTwos d
    | d == 1    = 1
    | even d    = 0
    | otherwise = d `div` 2

numNoTwos d
    | d == 1    = 2
    | otherwise = if even d then n else 2 * n
                  where d' = d `div` 2 - 1
                        n = sum [d' `choose` i | i <- [0..min d' 3]]

numYs :: Int -> Int

numYs d = numTwoTwos d + numOneTwos d + numNoTwos d

(OK, so I went ahead and handled the d = 1 case.)

Then you need only actually generate and filter the Ys from m' to the next higher power of ten, and from the largest power of ten <= n' to n'. The others you can count. Sorry, Little John; if you wanted the numbers themselves you should have taken it up with the Code Jam people.

UPDATE: Fixed my errors, and went ahead and dealt with the one-digit case in the various num*Twos functions.

Now, the only question is what to do about the portion of the interval from ⌈sqrt(m)⌉ to ⌊sqrt(n)⌋ that doesn't cover a whole range of values for a given power of ten? More precisely, how to handle it efficiently, without generating Y values any more than we have to?

UPDATE: revised numNoTwos to omit 0 for the one-digit case, since the statement of the problem says it's not dealing with 0.

Compare and contrast

...as my English teachers used to say.

"Finally, although the subject is not a pleasant one, I must mention PL/1, a programming language for which the defining documentation is of a frightening size and complexity. Using PL/1 must be like flying a plane with 7000 buttons, switches and handles to manipulate in the cockpit. I absolutely fail to see how we can keep our growing programs firmly within our intellectual grip when by its sheer baroqueness the programming language —our basic tool, mind you!— already escapes our intellectual control. And if I have to describe the influence PL/1 can have on its users, the closest metaphor that comes to my mind is that of a drug. I remember from a symposium on higher level programming language a lecture given in defense of PL/1 by a man who described himself as one of its devoted users. But within a one-hour lecture in praise of PL/1. he managed to ask for the addition of about fifty new “features”, little supposing that the main source of his problems could very well be that it contained already far too many “features”. The speaker displayed all the depressing symptoms of addiction, reduced as he was to the state of mental stagnation in which he could only ask for more, more, more... When FORTRAN has been called an infantile disorder, full PL/1, with its growth characteristics of a dangerous tumor, could turn out to be a fatal disease." --Edsger W. Dijkstra, from "The Humble Programmer" (ACM Turing Award Lecture 1972)

"This is the first Part of n, or lets [sic] say many entries in this blog. In total I hope to be able to cover most papers in 3-4 blog posts, giving the reader an overview over the suggestions and changes for C++ at the coming C++ Committee Meeting in April. In total there are 98 Papers, so I will skip some, but try to get as much covered as possible. I'll skip papers with Meeting Minutes for sure, and try to focus on those focusing on C++11 or C++14 features. As the papers are ordered by there [sic] Number (N3522 being the first), I'll go top down, so every blog post will contain different papers from different fields of C++ Standardization. As N352-24 are Reports about Active Issues, Defects and Closed Issues I'll skip them for now, also I will not read through Meeting Minutes and such." --Jens Weller, from "A look at C++14: Papers Part I"

Plus ça change...

Good company

Just saw something on Hacker News pointing to a tweet from John Carmack:

"I want to do a moderate sized Haskell project, and not fall to bad habits. Is there a good forum for getting code review/criticism/coaching?"

From a later tweet from Mr. Carmack: "The Haskell code I started on is a port of the original Wolf 3D...."

I'm sure I'm not at his level, but it's reassuring to be in good company... as I feel I am in my opinion of C++.

UPDATE: here's an HN link. This comment is particularly interesting; I'll have to keep an eye open for new ghc releases.

More attempts at optimization

I set up for profiling, and I made two changes. First, I put back the memoization of powers of ten. Second, it occurred to me that most of the time I invoke oddDigitPal, I'm invoking it twice with just the middle digit changing. So, I went for

oddDigitsPals :: Integer -> Int -> [Integer] -> [Integer]

oddDigitsPals topHalf nDigits middles =
    let noMiddle = topHalf * tenToThe (nDigits + 1) + backwards topHalf
    in [noMiddle + m * tenToThe nDigits | m <- middles]

I also changed over to Int (32-bit integers) for numbers of digits. I'm not worried at this point about going for fair and square numbers of over two billion digits. (Little John, you're on your own after that, OK?)

It did make some difference, based on profiling output. Most notably, it took about 14% off bytes allocated. Time was less impressive; with the profiling pulled, it took execution time down from not quite 1.9 seconds to not quite 1.7 seconds.

A C solution that I downloaded, compiled, and ran took more like 0.4 seconds. I'm bummed by that factor of slightly over four, but OTOH:

• I'm no expert at Haskell optimization.
• I'd much rather try to understand the Haskell I've written than the C or C++ solutions I've seen so far.

Just goes to show that Knuth was right

One other thing occurred to me as a possible optimization: I'm raising 10 to some integer power a lot. Why not try

powersOfTen = [10 ^ i | i <- [0..]]

tenToThe :: Int -> Integer

tenToThe n = powersOfTen !! n

which would memoize those pesky exponentiations? (!! lets you retrieve elements from lists sort of as if they were arrays, with "subscripts" starting at zero.)

It was easy enough to try out, but the results were disappointing. Even on my Eee 900A, with a 32-bit processor that you'd think would get the most benefit, the variations in time output from one run to the next were large enough that I can't say with certainty that it made any difference at all. Time output for the first large data set:

real    0m1.093s
user    0m1.068s
sys     0m0.016s

For the second large data set:

real    0m5.531s
user    0m5.472s
sys     0m0.044s

These are with the program compiled--I still haven't done the file opening code.

One thing that doesn't carry immediately over

One thing that you can do with a compiled Haskell program that doesn't lend itself to ghci is I/O redirection. That's why I have yet to run the large data sets against the code running under the interpreter--I'll have to modify main to take a file name and use it as input source.

I will do it, though. I'm highly motivated, because it's related to what caused me to start all this in the first place.

Since I don't have Code Jam time limits...

...I should improve the style.

Lisp is a wonderful language. It's the first (largely) applicative language. Having a simple form that all language constructs follow means you can write seriously powerful tools to manipulate programs without wasting your time on convoluted parsing.

Compare that with the horrors of parsing C++, which strictly speaking is impossible! Thanks to templates, you have to solve the halting problem to successfully parse C++. Even ignoring that, there are ambiguous constructs, with a rule that says which way to decide in the presence of ambiguity. I suspect it all goes back to Stroustrup's eschewing formal language-based tools when writing that first ad hoc preprocessor for what was then "C with Classes". (BTW, Perl has the same problem. Perhaps there's something about kitchen sink languages.)

All that said, Haskell style isn't Lisp style. Any serious Haskell programmer would look at my palindrome program and say I lean far too heavily on if ... then ... else and I don't take advantage of Haskell's pattern matching and guards when defining functions.

Guards!

Take the backwards function. At first it was

backwards n = backwards' n 0
    where backwards' n accum = if n < 10 then n + accum
                                         else backwards' (n `div` 10) (10 * (n `mod` 10 + accum))

and that's perfectly valid Haskell, but we can do better. How about this?

backwards n =
    let backwards' m accum
          | m < 10    = accum + m
          | otherwise = backwards' (m `div` 10) (10 * (m `mod` 10 + accum))
    in backwards' n 0

Here we are using guards, and while, thanks to the current layout, there is still wraparound, the code now honors the golden 80-column punched card that some things, most notably standards for source code layout, still bow down to in the computer world. (UPDATE: I waved a dead chicken at the layout, and it now fits.) Things to note:

• there's no = immediately after the list of formal parameters of backwards' the way there was in the first version. The = are after the guards.
• it may just be that I'm not familiar enough with Haskell syntax, but I had to turn things around and use a let rather than a where for the definition of the helper function backwards'. 
• speaking of style, did you notice I changed the name of the first parameter of backwards'? I decided I'd better do that, the same  way I'd avoid gratuitously redeclaring a variable in an inner block with the same name as one in an outer block.

Patterns

That's not the only way to avoid if then else, though. Remember choose?

nDigits `choose` nOnes =
    if      nOnes   == 0     then [0]
    else if nDigits == 1     then [1]
    else if nDigits == nOnes then [sum [10 ^ n | n <- [0..nDigits - 1]]]
    else                          concat [[0 + 10 * x | x <- (nDigits - 1) `choose` nOnes],
                                          [1 + 10 * x | x <- (nDigits - 1) `choose` (nOnes - 1)]]

Let's take advantage of pattern matching as well as guards:

_ `choose` 0 = [0]
1 `choose` _ = [1]
m `choose` n
    | m == n    = [sum [10 ^ i | i <- [0..m - 1]]]
    | otherwise = concat [[0 + 10 * x | x <- (m - 1) `choose` n],
                          [1 + 10 * x | x <- (m - 1) `choose` (n - 1)]]

For this, you should know that _ matches anything. An identifier matches anything, too, but you can refer to it on the other side of the =. Use _ when, given the other parameters and/or the previously checked patterns, the value of the function does not depend on the value that _ matches. The patterns are applied in order, by the way, so here 1 `choose` 0 will match  _ `choose` 0 and not 1 `choose` _.

Yes, I must admit that I chose shorter parameter names. While in the case of this test, I am passing in a number of digits, that is distracting from the function, so here I would say that my initial choice was wrong. It's not just a matter of making lines fit in 80 columns.

DRY

DRY here is an acronym, and it applies no matter what your programming language. DRY = Don't Repeat Yourself. I said, Don't Repeat Yourself.

I very much repeat myself in the original program when it comes to the formation of a palindrome from its upper half (or in the case of an odd number of digits, the upper half and the middle digit). That needs its own function, or rather functions.

In theory those are all the parameters needed; you can calculate the number of digits in the upper half. Since we have that value, though, I'm inclined to pass it in.

Also speaking of DRY, I should use these functions to generate the palindromes with two 2s rather than making a special case of them. Say what you mean, as clearly as possible.

[We pause here for some testing and checking of run times for the input file with values up to 10¹⁰⁰.]

It's hard to say whether it makes a difference; just running the same program repeatedly shows variance in the time output. Here's one towards the low end:

real    0m1.822s
user    0m1.804s
sys     0m0.012s

while here's one towards the high end:

real    0m1.905s
user    0m1.872s
sys     0m0.028s

IMHO the result is worth the minimal, if it even exists, increase in run time.

Here's the code after the smoke cleared:

import Prelude
import Data.List

backwards :: Integer -> Integer

backwards n =
    let backwards' n accum
          | n < 10    = accum + n
          | otherwise = backwards' (n `div` 10) (10 * (n `mod` 10 + accum))
    in backwards' n 0

evenDigitsPal :: Integer -> Integer -> Integer

evenDigitsPal topHalf nDigits = topHalf * 10 ^ nDigits + backwards topHalf

oddDigitsPal :: Integer -> Integer -> Integer -> Integer

oddDigitsPal topHalf nDigits middle = (topHalf * 10 + middle) * 10 ^ nDigits + backwards topHalf

choose :: Integer -> Integer-> [Integer]

_ `choose` 0 = [0]
1 `choose` _ = [1]
m `choose` n
    | m == n    = [sum [10 ^ i | i <- [0..m - 1]]]
    | otherwise = concat [[0 + 10 * x | x <- (m - 1) `choose` n],
                          [1 + 10 * x | x <- (m - 1) `choose` (n - 1)]]

chooseRange :: Integer -> [Integer] -> [Integer]

m `chooseRange` ns = concat [m `choose` n |  n <- ns, n <= m]

twoTwos :: Integer -> [Integer]

twoTwos n
    | even n    = [evenDigitsPal msd (n `div` 2)]
    | otherwise = [oddDigitsPal msd (n `div` 2) d | d <- [0,1]]
    where msd = 2 * 10 ^ (n `div` 2 - 1)

oneTwo :: Integer -> [Integer]

oneTwo n
    | even n    = []
    | otherwise = let halfN = n `div` 2
                      msd = 10 ^ (halfN - 1)
                  in [oddDigitsPal (msd + x) halfN 2 |
                        x <- (halfN - 1) `chooseRange` [0..1]]

noTwos :: Integer -> [Integer]

noTwos n
   | even n    = [evenDigitsPal (msd + x) halfN | x <- choices]
   | otherwise = concat [[oddDigitsPal (msd + x) halfN 0,
                          oddDigitsPal (msd + x) halfN 1] | x <- choices]
   where halfN   = n `div` 2
         msd     = 10 ^ (halfN - 1)
         choices = (halfN - 1) `chooseRange` [0..3]

possibles :: Integer -> [Integer]

possibles 1 = [0..3]
possibles 2 = [11, 22]
possibles n = sort (concat [noTwos n, oneTwo n, twoTwos n])

ys = concat [possibles nDigits | nDigits <- [1..]]

xs = [y * y | y <- ys]

palSquares :: Integer -> Integer -> [Integer]

palSquares m n = takeWhile (<= n) (dropWhile (< m) xs)

-- jagd's main, which I believe is set up to read stdin

solve i num = do
        [a, b] <- getLine >>= return . map read . words
        putStrLn $ "Case #" ++ show i ++ ": " ++ (show $ length (palSquares a b))
        if i < num
           then solve (i+1) num
           else return ()

main = do
       num <- getLine >>= return.read
       solve 1 num

UPDATE: Aside from something in a later message, the only things I can think of to do are

• adding concise comments that suffice to describe what the code does and how to the extent that the code itself is unclear
• changing some function names. possibles sort of made sense when we were range checking and/or palindrome checking right there, but the range checking is moved out and given the math, they're not possible, they're certain. Also, choose isn't quite right. We're not counting combinations, as the phrase "m choose n" is meant in standard math nomenclature--we're generating them. I'm not sure what to use in their place, though.

Palindromes again?

Actually, it's more Haskell than palindromes.

I've been looking over other people's solutions to the "Fair and Square" problem, and they are pretty impressive. As a still-inexperienced Haskeller, I am taking inspiration from then rather than feeling discouraged.

In particular, jagd's solution reminded me that I was not thinking Haskell. Haskell is lazy; values aren't calculated until you require them. That lets you define infinite data structures! You only get into trouble if you do something that requires the whole thing.

For example, [1..] is the infinite list of all positive integers. You can safely use it, as long as you only use a finite amount of it. For example, if you try

uptoSqrt n = [x | x <- [1..], x * x <= n]

and then ask for

uptoSqrt 25

you will get back 1, 2, 3, 4, and 5, and then the function call will hang. You and I  realize that \x -> x * x is monotonically increasing for positive numbers, but Haskell doesn't, and obligingly tries to check the infinitely many integers greater than 5 to see whether any have squares no greater than 25. That takes forever, so it never finishes.

To get the result you want, you will have to write something like

uptoSqrt n = [x | x <- (takeWhile (\y -> y * y <= n) [1..])]

What's takeWhile? Here's a declaration for it:

takeWhile :: (a -> Bool) -> [a] -> [a]

Give it a function that takes an a and returns a Bool, and a list of as, and it will hand you back a prefix of that list--all the items up to, but not including, the first one for which the function returns false.

With that definition,

uptoSqrt 25

gives you the

[1,2,3,4,5]

you expect, because we know that [1..] is in ascending order and the square function is monotonically increasing.

(There's nothing magic about takeWhile; if you ask for

takeWhile even [2, 4..]

it will print even numbers forever, because [2, 4..] is the infinite list of positive even integers, and hence takeWhile will never find an odd number to make it stop.)

That means that I can simply define the whole infinite list of "fair and square" numbers. To make it possible to get those in a certain range, I will have to generate them in ascending order, so back comes the sort that we omitted for a while. I'm going to do something that clashes a bit with a Haskell convention for patterns in function definitions, because it matches the nomenclature of the problem discussion, in which the "fair and square" numbers are referred to as X and the palindromes they are squares of as Y.

possibles nDigits =
    if      nDigits == 1 then [0..3]
    else if nDigits == 2 then [11, 22]
                         else sort (concat [noTwos nDigits, oneTwo nDigits, twoTwos nDigits])

ys = concat [possibles n | n <- [1..]]
xs = [y * y | y <- ys]

palSquares m n = takeWhile (<= n) (dropWhile (< m) xs)

dropWhile is kind of like takeWhile, but it gives you a suffix of the list you hand it, namely everything from the first item for which the function returns false on. If you drop the values less than m from the front of a sorted list, and then grab the values <= n from the front of what's left, you have precisely those values between m and n inclusive. (I'm being a bit free with my language here, of course. Haskell doesn't destroy the original list--at least not if anyone's looking.)

Note that we aren't using floorSqrt or ceilSqrt any more, so they can go away, and we thus aren't counting bits any more.  We're defining the whole list, so we aren't counting decimal digits, either, and digitsIn can go away as well.

So do infinite data structures provide any advantage? Well, consider that I have yet to actually solve the Code Jam problem as such. Remember that it requires reading a file with ranges of values m and n for which we have to print a line of a given format with the value of

length (palSquares m n)

Defining the infinite list of Xs is a form of memoization, remembering the values of a function you calculate so you don't have to do the work if you invoke it again with the same arguments. Suppose two of the input lines are

200 1000
40 500

Evaluating

palSquares 200 1000

will evaluate all the Xs between 0 and 1000 (actually one more past that, because you have to find the stopping place), so that when it comes time to process the next line, all the values between 40 and 500 will have already been determined, so all we have to do is select the values in range.

It's not what you'd want for serious memoization; after all, the input lines could all be within a small interval of the upper bound, so we would arguably have wasted the determination of all those lesser "fair and square" numbers. But it does have the potential of saving time. Let's see what happens with the Code Jam inputs (and with jagd's version of main)...

Cool! Of course, I didn't submit these officially--I'm using a portion of someone else's solution, and peeked at the analysis of the problem, and took far longer to write the program than the rules permit, so it would in no way be fair to take credit--but Google let me download the inputs, and confirmed that the output generated is correct. Here's the time output, first for large input 1, 10K lines with ranges a subset of [1, 10¹⁴]:

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

and then for large input 2, 1K lines with a range of [1, 10¹⁰⁰]:

real    0m1.796s
user    0m1.784s
sys     0m0.012s

Not shabby for time--and to refer back to Mr. Reeder's blog post, I'm sure that with those times, an interpreted version would have no trouble running in the eight minute limit Code Jam imposes--but how about memory? Time to instrument again... and sure enough, keeping that list, or perhaps those lists, around eats RAM. With large input 1, it uses 60+ kilobytes of RAM (the graph generated doesn't show the scale all the way up the Y (RAM) axis), and with large input 2, it uses about five megabytes of RAM. Even the five megabytes is small change these days, and now I suspect I wasn't looking closely at which run I was checking memory usage for last time around. I bet I thought a 10¹⁰⁰ run was a 10¹⁴ run, and I shouldn't have blamed the sort at all.

UPDATE: I'm truly amazed at the brevity of jagd's solution; you really should check it out. I will have a good time studying it further.  Just for the heck of it, I compiled it with ghc -O2, and tried it on the large input files.  For large input 1,

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

For large input 2,

real    0m14.410s
user    0m14.200s
sys     0m0.184s

I wonder what it's doing differently to make it run 63.8 times slower on the larger input rather than 7.84 times slower?

UPDATE:  part of it is the way it's keeping the list of Ys both in numeric and string form, with lots of conversions back and forth. That means more memory; on large input 2, jagd's program uses over 55 Mbytes of RAM, with a very interesting shape to the plot of RAM versus time.

I see that jagd also at least thinks he or she is generating some values that might not be valid Ys, given the

filter (ispalindrome . square)

applied to potential Ys. Maybe that's part of the RAM usage fluctuation.

("." is function composition in Haskell, and filter has a signature that should be familiar to you from takeWhile:

filter :: (a -> Bool) -> [a] -> [a]

filter takes a function and a list and hands  you back all the items in the list for which the function returns true.)

Fun with palindromes, and the joys of Haskell

An item in Hacker News caught my eye recently. (Ouch!) It linked to a post, "Ruby is Too Slow for Programming Competitions", by a gentleman named Clif Reeder. A summary: Mr. Reeder chose Ruby to try on a Google Code Jam qualification round problem, and found that his solution ran too slowly to meet the Code Jam time constraints (on the official runs, you have eight minutes from downloading the input to upload your output). A version written in the (compiled) Go language ran much faster.

Commenters on the HN item pointed out

• Mr. Reeder didn't think through the problem fully; there were more possible optimizations. In fact, there's a list of 387 people who used Ruby to successfully handle either the small sample input (which has a four-minute time limit) or the large or really large inputs (with the eight-minute time limit).
• Like C, the size of int in Go may vary, so the Go version might or might not be able to handle the larger inputs correctly.

There was also a long discussion of whether it's really worth writing programs in a higher-level language and dropping to a lower-overhead language for the time-critical parts (something I thought was SOP these days).

Anyway, I thought I'd add my voice to the chorus of "it's the algorithm, not the language" (or more accurately, if the algorithm isn't efficient, compiled vs. interpreted or high versus low level language doesn't matter) by writing in Haskell. The goal: writing code as clear as I can manage while still solving the problem quickly enough.

(I hasten to add that I am no Haskell expert, and I hope to see someone put Haskell to better, more idiomatic use while maintaining or improving clarity.)

Strictly speaking, I'm writing code that solves the main part of the problem:

Given natural numbers M and N, how many numbers between M and N inclusive are both (base ten) palindromes and squares of numbers which are (base ten) palindromes?

More about the full problem later. For now, let's refer to an example of the kind of number we want as X, where X is a palindrome and also equal to the square of a palindrome Y, as the Code Jam analysis of the problem does.

When you look at the problem, some possible optimizations over the brute force method come to mind:

1. Rather than look in [M, N] for Xs, look over values in [⌈sqrt(M)⌉, ⌊sqrt(N)⌋] for Ys. Mr. Reeder took advantage of this in his program; it helps a great deal.
2. Rather than looking at each number in the range for palindromes, better to generate the palindromes in the range. A palindrome is determined by half its digits (unless it has an odd number of digits, in which case there's another digit in the middle that's not constrained to match any other digit), so that gives another reduction on the same order as (1).
3. If Y squared has to be a palindrome, Y is still further constrained. Suppose Y = 4.....4. Then Y squared has the form n....6, where 16 ≤ n < 25, so it can't be a palindrome. A similar argument rules out five through nine, cutting the work down by two thirds.

That, alas, is as far as I took the analysis. With that much, straightforward Haskell code will spit out all the "fair and square" numbers, as the problem statement calls them, between 1 and 10¹⁴ in about a second on my Eee 900A netbook. That range includes all the values in the smaller large input file.

Fortunately, it's possible to do much better than that! Google's full analysis of the problem can be found on the Code Jam site. To summarize it:

1. Given the constraint on Y's leading and trailing digits, there can't be a carry out of the most significant digit when you add up the partial products of Y * Y to get X, so if Y has n digits, X has 2n - 1 digits, an odd number.
2. In fact, there's no carry out of any position in the sum of the partial products of Y and Y; that's both a necessary and sufficient condition for X = Y * Y to be a palindrome if Y is a palindrome.
3. The middle digit of X is thus (since Y is a palindrome) the sum of the squares of the digits of Y, and must be no bigger than 9.

It follows immediately from (3) is that the only Y containing the digit 3 is 3 itself. All other possible Y values can only contain the digits 0, 1, and 2, and can be broken down into three cases:

• No 2s at all. They can have up to nine 1s if they have an odd number of digits, eight if they have an even number of digits.
• One 2. The 2 must be the middle digit, requiring an odd number of digits, and the rest of Y can have no more than four ones (if there were five, Y couldn't be a palindrome). Aside from 2, there must be at least two ones, the most and least significant digits.
• Two 2s. This allows at most one 1, and that only in the middle, which can only happen if there are an odd number of digits.

 Any other digits in Y must be zeroes.

Just constraining the digits to 0, 1, and 2 did speed things up somewhat; the resulting program managed to get up into the 48-digit Xs in about three minutes, but that's nowhere near fast enough to handle the largest test input, which has values up to 10¹⁰⁰, in the eight minutes Code Jam limits you to. To do it right, you have to seriously skip values that don't meet the sum of squares constraint.

Here's the code I ended up with. I split out the one- and two-digit cases so that the noTwos, oneTwo, and twoTwos functions don't have to bother with them.

I should point out that

• I did this without the time constraints of Code Jam, over about eight days after the Hacker News post
• Of course, you're not seeing the mistakes I made over that eight days!
• The integer square root function came from an answer to a Stack Overflow question, which in turn was taken from Crandall and Pomerance, "Prime Numbers: a Computational Perspective"
• The trivial Main is adapted from Chapter 25 of Real World Haskell having to do with profiling (about which more later)

import Prelude
import System.Environment
import Text.Printf

backwards :: Integer -> Integer

backwards n = backwards' n 0
    where backwards' n accum = if n < 10 then n + accum
                                         else backwards' (n `div` 10) (10 * (n `mod` 10 + accum))

numPalindrome :: Integer -> Bool

numPalindrome n = (n == backwards n)

twoTwos :: Integer -> [Integer]

twoTwos nDigits =
    let justTwos = 2 * 10 ^ (nDigits - 1) + 2
    in
        if even nDigits then [justTwos]
                        else [justTwos, justTwos + 10 ^ (nDigits `div` 2)]

choose :: Integer -> Integer-> [Integer]

nDigits `choose` nOnes =
    if      nOnes   == 0     then [0]
    else if nDigits == 1     then [1]
    else if nDigits == nOnes then [sum [10 ^ n | n <- [0..nDigits - 1]]]
    else                          concat [[0 + 10 * x | x <- (nDigits - 1) `choose` nOnes],
                                          [1 + 10 * x | x <- (nDigits - 1) `choose` (nOnes - 1)]]

chooseRange :: Integer -> [Integer] -> [Integer]

chooseRange nDigits onesRange = concat [nDigits `choose` nOnes | nOnes <- onesRange, nOnes <= nDigits]

oneTwo nDigits =
    if even nDigits then []
                    else
                         let msd = 10 ^ (nDigits `div` 2 - 1)
                         in  [((msd + x) * 10 + 2) * 10 ^ (nDigits `div` 2) + backwards (msd + x) |
                                      x <- chooseRange (nDigits `div` 2 - 1) [0..1]]

noTwos :: Integer -> [Integer]

noTwos nDigits =
    let halfN   = nDigits `div` 2
        msd     = 10 ^ (halfN - 1)
        choices = chooseRange (halfN - 1) [0..3]
    in
        if even nDigits then [         (msd + x)           * 10 ^ halfN + backwards (msd + x) | x <- choices]
                        else concat [[((msd + x) * 10 + 1) * 10 ^ halfN + backwards (msd + x),
                                       (msd + x) * 10      * 10 ^ halfN + backwards (msd + x)] | x <- choices]

possibles :: Integer -> [Integer]

possibles nDigits =
    if      nDigits == 1 then [0..3]
    else if nDigits == 2 then [11, 22]
                         else concat [noTwos nDigits, oneTwo nDigits, twoTwos nDigits]

digitsIn :: Integer -> Integer -> Integer

digitsIn base n = digitsIn' n 1
    where digitsIn' n count = if n < base then count else digitsIn' (n `div` base) (count + 1)

candidates :: Integer -> Integer -> [Integer]

candidates m n = [x | x <- concat (map possibles [digitsIn 10 m..digitsIn 10 n]), x >= m, x <= n]

floorSqrt :: Integer -> Integer

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

ceilSqrt :: Integer -> Integer

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

palSquares :: Integer -> Integer -> [Integer]

palSquares m n = [x * x | x <- candidates (ceilSqrt m) (floorSqrt n)]

main = do
    [d] <- map read `fmap` getArgs
    printf "%d\n" (length (palSquares 1 d))

Now, the times mentioned above are for runs without main, under ghci, an interpreter, using a version that sorted the results of possibles to make the output easier to check against the values at "Palindromic Square Numbers". I should also mention that under ghci it was able (again, with the sort) to generate the list of "fair and square" numbers between 1 and 10¹⁰⁰ in about 28 seconds on my 2.8 GHz Athlon II X4 630 processor. Some of that time was I/O.

Profiling meant using the compiler, and what a difference that made! The first version of main didn't bother to read the command line arguments, and instead wired in 10¹⁰⁰ as the upper bound. Compiling with -O2, it ran fast enough that I made the change to read the command line arguments lest it have optimized away all the calculations. That done, time output for a run, again using 10¹⁰⁰ as the upper bound, showed

real  0m0.495s
user  0m0.480s
sys   0m0.012s

With the sort pulled, the time output was

real  0m0.389s
user  0m0.380s
sys   0m0.004s

The sort also turned out to be a major source of memory usage. With sort, it was using nearly four megabytes of RAM. Without, the profile showed it using more like 40K.

Now, you couldn't use the interpreted version in the straightforward way to successfully pass the Code Jam large and larger input tests; you'd want to take advantage of their giving the range of values, generate the list and then convert it to a structure that lends itself to quickly determining how many values are in a given range. (Sorted list, perhaps a search tree...)

The compiled version, on the other hand...  the smaller large input file has no more than 10000 lines, confined to [1, 10¹⁴].  Results of time on that range:

real    0m0.003s
user    0m0.004s
sys     0m0.000s

If every input line covered the whole range, that would be thirty seconds, plus the time for I/O--well within the eight minute range.

The larger large input file has no more than 1000 lines, confined to [1, 10¹⁰⁰]. If every input line covered the whole ranges, that would be 389 seconds plus I/O time. That's cutting it rather close, but it's still possible.

So, have we shown what we set out to show? Clearly compiled Haskell is far faster than interpreted Haskell. (To give a better comparison, since I ran the interpreted cases on my Eee 900A--the compiled version ran in about 2.5 seconds on the larger range on the Eee.) But that's not the point; the optimized algorithm could trivially run the first large data set even interpreted, and with the precalculation could probably do the larger large data set interpreted.

All that said, the revelation for me was this: it was incredibly easy to write the code, straightforward to test and debug in the interpreter, and I had only to add import lines and the simple main to be able to compile it. Total: sixty-three (nonblank) lines of code. It's helped along by Haskell having arbitrary-precision integers and list comprehensions, which I used liberally.

One could argue that this problem is one particularly well-suited to Haskell, but for me this was a sufficiently enjoyable task that I will do more with Haskell, and I urge you to give it a try too.

UPDATE: Looks like Ruby has Bignum, and various people have posted fairly concise-looking constructs to get the effect of list comprehensions.  It might be fun to rewrite the above in Ruby.

UPDATE: Here are the solutions in Haskell from the Code Jam web site. Check out jadg's code; it's very impressively short. I'll be studying that one.

Night of the Living Blog

A good friend nudged me to start putting up stuff I am working/have worked on, so Aging Nerd Notes is shuddering back to activity after what, eight years? I'm amazed Google has kept it around so long.

O would some god the giftie gie us...

Recently I've gotten some email from word-of-mouth-connection.com. People can register that they have experience with a particular person or business, say something about how well they know said entity and how recent their experience is, and then others who want to find out about said entity can go looking through the database and inquire via email; the contacts are done anonymously.

At the time of writing, there are apparently three people who claim to know me well who have registered that fact, and hence, I suppose, their willingness to tell others about me, and one person looking to find out about me.

I guess that it's good that word-of-mouth-connection.com tells me that this is happening...but to be honest, it creeps me out. That, in turn, probably points out my tendency to assume the worst. For all I know, the three people are happily telling anyone who asks that I'm a heck of a guy, but I would have been happy to continue to not think about people talking about me behind my back. If it's anonymous, I suppose I could email the three and ask what they think, and the one and ask what he or she wants to know, but I'm not sure I want to find out.

Ah, well. As the guy said at the end of Asimov's "The Dead Past," happy goldfish bowl, everybody.

OK, now say "shibboleth"...

Computers are far more widespread than they used to be--which perforce means that most of the people who use them are ignorant of how they work, just as most people couldn't tell you how their TVs, cars, or phones work. (Not that there's anything wrong with that.) So, how can you tell the posers from the real deal? One sure mark of someone who doesn't know diddly about computers is this: referring to mass storage space (hard disk, CD-ROM, DVD-ROM) as "memory."

What are some other misstatements that reveal computer ignorance? Someone should compile a list of them--so send your favorites our way, please!

Not the best intro, I guess...

I probably would be more favorably inclined to Joan of Arcadia had I watched it from the beginning, but I started with this evening's episode.

Joan of Arcadia is sort of the mutant offspring of Touched by an Angel and Calvin and Hobbes. Did you ever notice how, if agents of the divine really did all the stuff that the angels on TbaA did, there would be no doubt at all of the existence of a non-denominational Christian God? JoA bypasses that problem by making God an ever-shifting Hobbes to Joan's Calvin. God/Hobbes tells Joan/Calvin what to do, and occasionally gives the stray oblique lecture (apologies to Brian Eno). Joan is supposedly trying to figure out just what God has in mind, but if the producers have any sense, it will never happen; just as the unseen monster is always scarier than the guy in the rubber suit who finally emerges, any supposed divine plan that the writers come up with will necessarily seem lame. (Though come to think of it, if the blatant hint of the title holds, any end is going to be pretty definitive...)

Tonight, there are three plot threads: Joan's family's reactions to her brother's paralysis, Joan's interactions with the in crowd at high school, and the bringing in of a "psychic" in a kidnapping case. Joan's father, the new chief of police in the town they've moved into, is the voice of rationality, and doesn't at all appreciate the psychic. In a tense interaction between the two, the psychic claims to have gained her talents after a near-death experience, and asks him whether it's "the tragedy" that makes him so resistant to her presence.

I am extremely disappointed in the writers of this show; they had and blew a wonderful opportunity for the father to take the "psychic" to task. What he should have said at that point is, "How dare you live off the suffering of others! You've just begun a 'cold reading,' asking vague leading questions to try to pump me for information that you can then claim to have gotten from the spirits or via your supposed powers. Nearly everybody has had some tragedy, and if I happen to be someone who has managed to avoid personal tragedy, you could claim that you were referring to the disappearance. You can sucker a lot of people, even intelligent people, with that scam, but not this time. Now get out!"

Why is it that the media so often pander to pseudoscience like this?
Can't the rational characters get the upper hand once in a while?

Verbal Abuse as Entertainment

When I grew up, my parents always told me that there was a sort of person who needed to tear down others to try to salvage his or her nonexistent self-esteem, and that I should not be that sort of person. Judging by radio and TV today, though, such people seem to be doing awfully well.

Verbal abuse is now a popular entertainment form. Weakest Link, American Idol, Cupid... all “reality shows” that showcase and derive their popularity from gratuitous viciousness. David Letterman has made a career of finding people from “flyover country” and interviewing them in such a way that he and his audience can derive their sick jollies from feeling hip and superior. Such interviews and prank calls are a major feature of, for example, Mancow Muller’s syndicated radio show. We’ve already mentioned “Dr. Laura.” One sees the actions of, say, Simon Cowell rationalized as “honesty” or “tough love”—but one can be honest without being vicious or abusive.

So why are these TV and radio shows doing so well? Have we turned into the people that our parents warned us not to be, or a cowardly variant thereof who just enjoy watching someone else do the dirty work? It doesn’t speak very well for us.

Grumble, Grumble...

I have just one thing to say to the Academy of Motion Picture Arts and Sciences for not giving the Best Animated Feature award to Lilo and Stitch:

Meega na la kweesta!

Sic Transit

The Islamic cultures once were important centers of philosophy, science, mathematics, and literature. Just go through a technical dictionary looking for words that start with al- (Arabic for "the"): algebra, algorithm, alcohol, aldehyde... We have them to thank for the Alhambra, Leyli and Majnun (whose influence stretches through Tristan and Isolde, Romeo and Juliet, and Derek and the Dominoes), and the poetry of Rumi. (Speaking of the Alhambra--the Moors were considerably more tolerant than Los Reyes Católicos, who offered the Jews the choice of conversion or eviction shortly after coming to power.)

So, how on earth could such a culture produce the vile beasts who now rule Iraq, or who make sure schoolgirls dressed "indecently" don't escape burning buildings? Are such people the norm, or are they as reviled in their culture as Nazis and their deluded followers are in the West? I'd like to think the latter.