Project Euler #4: Largest Palindrome Product

Spoiler alert! This page is part of a series on solutions to Project Euler problems. If you prefer to solve problems yourself, do not read on!

This post is literate Haskell; you can load the source into GHCi and play along.

First some boilerplate.

module ProjectEuler004 where

import Data.List
import Data.Maybe
import System.Exit

Problem 4 from Project Euler:

“Palindromeness” is very different from most interesting properties of numbers; it is an artifact of one of many possible representations of the number. In this case, the representation in base 10.

Let’s start by nailing down a definition for “palindrome”.

The first thing I’ll do is write helper functions to convert a number to its base 10 digits and back.

-- base 10 digits of n in 'little endian' order
-- (least significant digits first)
toDigits :: Integer -> [Integer]
toDigits n = if n <= 0
  then []
  else (n`rem`10) : toDigits (n`quot`10)

fromDigits :: [Integer] -> Integer
fromDigits = sum . zipWith (*) (map (10^) [0..])

numDigits :: Integer -> Integer
numDigits n = sum $ map (const 1) $ toDigits n

Sanity check:

test_digits :: Integer -> Bool
test_digits n = n == (fromDigits $ toDigits n)

And a test:

$> all test_digits [1..1000000]
True

Note that our lists of digits come out “backward”; that is, least significant first.

$> toDigits 12345
[5,4,3,2,1]

Then we can detect whether a given number is a palindrome in base 10.

is_palindrome_10 :: Integer -> Bool
is_palindrome_10 n = let ds = toDigits n in
  ds == (reverse ds)

Now we’re not just looking for the largest palindrome of a given length; that would be easy – the string of all 9s is the largest palindrome with a given number of digits. Instead, we want the largest palindrome that is the product of two 3-digit numbers. The most obvious solution is to list all the products of two 3-digit numbers, filter for the palindromes, and find the max.

-- the triple (a,b,a*b) which yields the largest
-- palindrome product a*b among the pairs of
-- t-digit numbers a and b
pe4' :: Integer -> (Integer, Integer, Integer)
pe4' t =
  let
    thd (_,_,x) = x
    min = 10^(t-1)
    max = 10^t - 1
  in
    maximumBy (\x y -> compare (thd x) (thd y)) $
    filter (is_palindrome_10 . thd) $ 
    [(a,b,a*b) | a <- [min..max], b <- [a..max]]

This pe4' works well enough:

$> pe4' 1
(3,3,9)
$> pe4' 2
(99,91,9009)
$> pe4' 3
(993,913,906609)
$> pe4' 4
(9901,9999,99000099)

Two observations. First, there appears to be a pattern emerging in these results. Second, unfortunately I can’t feasibly explore that pattern further (yet) because as t gets larger, pe4' t gets way slow. It’s not hard to see why; pe4' t constructs a list of \[\binom{9 \cdot 10^{t-1}}{2} = \frac{(9 \cdot 10^{t-1})(9 \cdot 10^{t-1} - 1)}{2}\] candidates to check by brute force.

To see if we can find a better way, let’s look at the simplest version of this problem: finding the largest palindrome product among products of 1-digit numbers. To this end, consider the multiplication table for 1-digit numbers below.

\[\begin{array}{c|ccccccccc} & 9 & 8 & 7 & 6 & 5 & 4 & 3 & 2 & 1 \\ \hline 9 & 81 & 72 & 63 & 54 & 45 & 36 & 27 & 18 & 9 \\ 8 & & 64 & 56 & 48 & 40 & 32 & 24 & 16 & 8 \\ 7 & & & 49 & 42 & 35 & 28 & 21 & 14 & 7 \\ 6 & & & & 36 & 30 & 24 & 18 & 12 & 6 \\ 5 & & & & & 25 & 20 & 15 & 10 & 5 \\ 4 & & & & & & 16 & 12 & 8 & 4 \\ 3 & & & & & & & 9 & 6 & 3 \\ 2 & & & & & & & & 4 & 2 \\ 1 & & & & & & & & & 1 \\ \end{array}\]

The first thing to jump out at me is that the vast majority of these products are not palindromes. That’s not really surprising I guess – there are \(9 \cdot 10^{2t-1}\) different \(2t\)-digit numbers, but only \(9 \cdot 10^{t-1}\) of them are palindromes, so to a first approximation the probability that an arbitrary \(2t\)-digit number is a palindrome is about \(10^{-t}\).

So maybe instead of searching among the products of \(t\)-digit numbers for palindromes, it would be faster to search among the palindromes for products of \(t\)-digit numbers.

A Digression on Palindromes

The rarity of palindromes highlights a disturbing possibility. The problem of finding the largest palindrome product implicitly assumes that at least one palindrome product must exist, but this is not obvious to me. After some experimenting, I am convinced that palindrome products exist among pairs of \(t\)-digit numbers for any \(t\). I got a little distracted thinking about this probem, so I’ll just dump the results here.

The first result allows us to recognize a family of palindromes constructed from other palindromes.

This result lets us recognize some palindromes.

The next result lets us construct new palindrome products from old ones.

Let’s make \(H_v\) and \(K_v\) executable.

h_ :: Integer -> Integer -> Integer
h_ v a = a * (1 + 10^(2*t+v))
  where t = numDigits a

k_ :: Integer -> Integer -> Integer
k_ v b = b * (1 + 10^(2*t+v))
  where t = numDigits b

And the next result gives us a concrete family of palindrome products with factors of even digit length.

Let’s make \(Q_{t,m}\), \(A_{t,m}\), and \(B_{t,m}\) executable.

q_ :: Integer -> Integer -> Integer
q_ t m = sum $ map (\k -> 10^(m*k)) [0..(t-1)]

a_ :: Integer -> Integer -> Integer
a_ t m = q_ (2*t) m

b_ :: Integer -> Integer -> Integer
b_ t m = 1 + (10^(t*m))*(10^m - 1)*(q_ t m)

c_ :: Integer -> Integer -> Integer
c_ t m = (q_ t m)*(1 + 10^(3*t*m))

Sanity check:

test_abc :: Integer -> Integer -> Bool
test_abc t m =
  let
    a = a_ t m; b = b_ t m; c = c_ t m
  in and
    [ a*b == c
    , (numDigits a) == 2*t*m - m + 1
    , (numDigits b) == 2*t*m
    , (numDigits c) == 4*t*m - m + 1
    , is_palindrome_10 c
    , is_palindrome_10 (2*c)
    ]

And a test:

$> and [test_abc t m | t <- [1..10], m <- [1..10]]
True

This gives an infinite family of palindrome products. Note that if \(m = 1\), then both factors have the same number of digits – \(2t\) – and the product has \(4t\) digits.

Here are the first 5 examples.

Note that row \(k\) in this table has factors \(A\) and \(B\) with \(2k\) digits, and the (palindrome) product has \(4k\) digits.

Taking the first row, \(A = 11\) and \(B = 91\), we can use these in our \(H_v\) and \(K_v\), with \(v = 1\), we get another family of palindrome products.

Note that if \(v = 2k+1\) then this table gives two \((2k+7)\)-digit numbers whose (palindrome) product has \(4k+14\) digits.

Along with the observation that \(993 \cdot 913 = 906609\) and \(11011 \cdot 91091 = 1003003001\), we can say the following.

Right, because the whole point of this digression was to establish this fact about largest palindrome products.

Meanwhile

Our alternate strategy for this problem was to search among the palindromes for \(k\)-digit factorizations; and now we know that it’s enough to look among the \(2k\)-digit palindromes – we are guaranteed to find a factorization there. This simplifies the search space a bit.

So the new strategy is to generate the \(2k\)-digit palindromes in reverse order and look for the first one that factors are a product of two \(k\)-digit numbers.

-- the 2k-digit palindromes
palindromes :: Integer -> [Integer]
palindromes k = map (fromDigits . make) (digits k)
  where
    -- turn a list into a palindrome
    make ds = ds ++ reverse ds

    digits :: Integer -> [[Integer]]
    digits k = foo k []
      where
        foo 0 _ = [[]]
        foo k z = do
          d  <- [9,8,7,6,5,4,3,2,1] ++ z
          ds <- foo (k-1) [0]
          return (d:ds)

Now given a \(2k\)-digit palindrome, we want to know whether it factors as a product of two \(k\)-digit numbers. Note that if \(N = AB\), then without loss of generality \(A \leq \sqrt{N}\) and \(\sqrt{N} \leq B\). So it is enough to search for a factorization of \(N\) where \(10^{t-1} \leq A \leq \lfloor \sqrt{N} \rfloor\).

Here’s a utility function to find \(\lfloor \sqrt{N} \rfloor\) by bisection.

-- find t such that t^2 <= n < (t+1)^2
floor_sqrt :: Integer -> Integer
floor_sqrt n =
  let
    bisect (a,b) = if b-a <= 1
      then a
      else let m = (b+a)`quot`2 in
        if m^2 <= n
          then bisect (m,b)
          else bisect (a,m)
  in
    if n < 0
      then error "floor_sqrt: negative argument"
      else bisect (1,n)

Sanity check:

test_floor_sqrt :: Integer -> Bool
test_floor_sqrt n = let q = floor_sqrt n in
  (q^2 <= n) && (n < (q+1)^2)

And a test:

%> all test_floor_sqrt [1..10000]
True

Also, if we wish to factor a \(2k\)-digit number as a product of \(k\) digit numbers, we should search “down” from \(\lfloor \sqrt{N} \rfloor\) to \(10^{t-1}\) rather than the reverse. This is because smaller factors are likely to give quotients with too many digits.

-- search for a factorization of n into t-digit factors
does_factor :: Integer -> Integer -> Maybe (Integer, Integer)
does_factor t n =
  let
    check m
      | m < 10^(t-1) = Nothing
      | n`rem`m == 0 = if numDigits (n`quot`m) == t
          then Just m
          else check (m-1)
      | otherwise = check (m-1)
  in
    case check (floor_sqrt n) of
      Nothing -> Nothing
      Just m  -> Just (m, n`quot`m)

pe4'' :: Integer -> (Integer, Integer, Integer)
pe4'' k =
  let
    (a,b) = head $ catMaybes $ map (does_factor k) (palindromes k)
  in (a,b,a*b)

This pe4'' is still slow, as we still check an exponential number of cases. But we can squeeze out a couple more values.

$> pe4'' 2
(91,99,9009)
$> pe4'' 3
(913,993,906609)
$> pe4'' 4
(9901,9999,99000099)
$> pe4'' 5
(99681,99979,9966006699)

Interesting! It looks like the form of the largest palindrome product depends on whether the factors have an even or odd number of digits.

Anyway, the final answer is:

pe4 :: Integer
pe4 = let (_,_,x) = pe4'' 3 in x

main :: IO ()
main = do
  let
    success = and
      [ all test_digits [1..1000000]
      , and [ test_abc t m | t <- [1..10], m <- [1..10] ]
      , all test_floor_sqrt [1..10000]
      ]
  if success
    then putStrLn $ show pe4
    else exitFailure