Project Euler Problem 4

A palindromic number reads the same both ways. The largest palindrome made from the product of two $2$-digit numbers is $9009 = 91 \times 99$.

Find the largest palindrome made from the product of two $3$-digit numbers.

The naive approach

First let's write a function that can quickly test whether an integer $n$ is a palindrome. This can be done pretty easily and elegantly:

function is_palindrome(n)
    n = abs(n)  # handle negative numbers
    return string(n) == reverse(string(n))
end

The most naive solution checks every product $ij$ of two $n$-digit numbers this way:

integer_type(num_digits) = num_digits <= 18 ? Int64 : num_digits <= 38 ? Int128 : BigInt

function largest_palindrome_product_naive(n, ::Type{T}=integer_type(2n)) where {T}
    lower_limit, upper_limit = T(10)^(n - 1), T(10)^n - 1
    max_palindrome = zero(T)
    best_i, best_j = zero(T), zero(T)

    for i in lower_limit:upper_limit
        for j in lower_limit:upper_limit
            product = i * j
            s = string(product)
            if s == reverse(s) && product > max_palindrome
                max_palindrome = product
                best_i, best_j = i, j
            end
        end
    end

    return (palindrome=max_palindrome, factors=(best_i, best_j))
end

The products have up to $2n$ digits, so the integer type T defaults to the smallest one they fit in: Int64 up to $9$ digits, and Int128 after that.

It finds the answer in 34.875 ms but allocates 74.16 MiB along the way, since it builds two new strings for each of the $900^2 = 810,000$ products. The 6-digit case has a million times more products to check, so it would take hours.

Searching the largest products first

Testing for palindromes with strings allocates memory for two strings every time. We can instead take a numerical approach where we reverse the number by extracting digits from right to left and rebuilding the number:

@inline function is_palindrome(n; base=10)
    n = abs(n)  # Handle negative numbers
    base = oftype(n, base)
    original = n
    reversed = zero(typeof(n))

    while n > 0
        n, digit = divrem(n, base)
        reversed = reversed * base + digit
    end

    return reversed == original
end

We can then use this function to search for the largest palindrome made from the product of two $n$-digit numbers:

function largest_palindrome_product_pruned(n, ::Type{T}=integer_type(2n); max_product=nothing) where {T}
    lower_limit, upper_limit = T(10)^(n - 1), T(10)^n - 1
    max_palindrome = zero(T)
    best_i, best_j = zero(T), zero(T)

    for i in upper_limit:-1:lower_limit
        # break early if we can't find a larger palindrome
        if i * upper_limit < max_palindrome
            break
        end

        # Stop at j = i to avoid duplicate combinations as i*j == j*i
        for j in upper_limit:-1:i
            product = i * j

            # Skip if product exceeds max_product constraint
            if !isnothing(max_product) && product >= max_product
                continue
            end

            if product < max_palindrome
                break
            end

            if is_palindrome(product) && product > max_palindrome
                max_palindrome = product
                best_i, best_j = i, j
            end
        end
    end

    return (palindrome=max_palindrome, factors=(best_i, best_j))
end

So we search through all products $ij$ in descending order to find the largest palindrome. The search is sped up in a few ways. First, we iterate from largest to smallest values since we're searching for a maximum. This lets us prune the search, terminating early once $ij$ can no longer exceed the current maximum palindrome found so far. We also added a max_product option that only considers products below a given value, which is what the HackerRank version asks for. With up to $100$ queries per run though, the submission instead precomputes every palindrome product of two $3$-digit numbers once and binary searches that sorted list for each query.

Benchmarking the 3-digit case we find the solution largest_palindrome_product_pruned(3) in 15.924 μs, which is 2,190× faster than the naive solution.

For the 6-digit case we call largest_palindrome_product_pruned(6) to find a maximum palindrome of $999,000,000,999 = 999,001 \times 999,999$ in 4.568 ms.

The 8-digit case takes 641.611 ms. We can still do the 9-digit case, where largest_palindrome_product_pruned(9) finds $999,900,665,566,009,999 = 999,920,317 \times 999,980,347$, but that takes almost a minute.

Beyond that, the product of two 10-digit numbers no longer fits in a 64-bit integer, so T switches to Int128. The 12-digit case will probably also take a lot longer than 1 minute! Besides Int128 operations taking more CPU cycles than Int64 operations, the solution time scales superlinearly.

There are $R = 9 \times 10^{n-1}$ numbers with $n$ digits. In the worst case our solution takes $\mathcal{O}(R^2 n)$ time as it checks all $(i, j)$ pairs of which there are $R^2$ pairs and calls to is_palindrome take $\mathcal{O}(n)$ time. The early termination optimizations will reduce the number of pairs to check by a lot but in the worst case you're still searching a large space.

Fermat's factorization method

I asked Claude how we can do better and it suggested Fermat's factorization method, an idea it most likely picked up from solutions to LeetCode 479.

Both solutions so far loop over pairs of factors and check whether their product is a palindrome. We can flip the search around and instead check each palindrome starting with the largest one, and for each one ask whether it's the product of two $n$-digit numbers.

There are far fewer palindromes than products, and the answer is always close to the top. So hopefully we don't have to search for long. Getting the factors means trial division, which is slow, but because both factors are close to $10^n$, we can check each palindrome with just an integer square root or two!

Let $2n$ be the number of digits in the palindrome. Then let $B = 10^n$ and write both factors as distances below $B$, so that $x = B - u$ is one factor and $y = B - v$ is the other. Multiplying them out we get

Multiplying by $B = 10^n$ just shifts a number $n$ digits to the left. So as long as $uv < B$, the upper $n$ digits of $xy$ are $B - (u + v)$ and the lower $n$ digits are $uv$.

A $2n$-digit palindrome is determined by its upper half $U$, since the lower half is just $U$ with its digits reversed. So let's write the upper half as $U = B - m$ and the lower half as $L = \mathrm{reverse}(U)$.

Then $m = 1, 2, 3, \dots$ lists the palindromes from the largest down: $999999$, $998899$, $997799$, and so on. The palindrome with top half $B - m$ is a product $xy$ exactly when both halves match, i.e.

Knowing the sum and product of two numbers is enough to work out the numbers themselves! Since $(t - u)(t - v) = t^2 - (u + v)t + uv$, the numbers $u$ and $v$ are the two roots of the quadratic $t^2 - mt + L = 0$, so

These are whole numbers only when the discriminant $D$ is a perfect square, which a single isqrt can tell us. In fact $D = (u + v)^2 - 4uv = (u - v)^2$, so $\sqrt{D}$ is just the gap between $u$ and $v$.

We assumed $uv < B$, but what if $uv \geq B$? a larger $uv$ spills into the upper half, just like carrying when adding by hand. If the carry $c$ is how many whole $B$'s fit into $uv$, the upper half is $B - (u + v) + c$ and the lower half is $uv - cB$. So the palindrome with upper half $B - m$ is $xy$ when

We don't know the carry ahead of time, so we try $c = 0, 1, 2, \dots$, each with its own discriminant

$D$ only shrinks as $c$ grows so we can stop at the first negative $D$ as we cannot get a perfect square by increasing $c$ anymore.

Solutions to the LeetCode version of this problem, which only goes up to 8 digits, use the same idea without the carry. As far as I can tell it goes back to a 2017 reply by nizametdinov on LeetCode's discussion board, which CerebrumMaize and Minghan turned into popular write-ups. As hvsavage pointed out in the comments on Minghan's LeetCode post, dropping the carry breaks at 9 digits. I guess Claude's contribution is the carry.

Putting it all together:

function reverse_digits(h)
    ten = oftype(h, 10)
    reversed = zero(h)

    while h > 0
        h, digit = divrem(h, ten)
        reversed = reversed * ten + digit
    end

    return reversed
end

function largest_palindrome_product_fermat(n, ::Type{T}=integer_type(n)) where {T}
    n >= 2 || throw(ArgumentError("expected n ≥ 2 digits, but got n = $n"))

    B = T(10)^n
    uv_max = B - B ÷ 10  # largest u or v that keeps x and y at n digits

    for m in one(T):uv_max
        H = B - m              # top half of the palindrome, from 99...9 downwards
        L = reverse_digits(H)  # bottom half

        c = zero(T) # carry
        while true
            s = m + c                  # u + v
            D = s^2 - 4 * (c * B + L)  # (u - v)^2
            D < 0 && break             # D shrinks as c grows, so no larger carry can work either

            r = isqrt(D)
            if r^2 == D
                u, v = (s - r) ÷ 2, (s + r) ÷ 2
                if u >= 1 && v <= uv_max  # x and y must both have n digits
                    return (palindrome=widemul(H, B) + L, factors=(B - v, B - u))
                end
            end

            c += 1
        end
    end

    # No 2n-digit palindrome is a product of two n-digit numbers
    return (palindrome=zero(widen(T)), factors=(zero(T), zero(T)))
end

uv_max keeps both factors at $n$ digits, and the inner loop tries carries until $D$ goes negative. This search can't miss a factorization because every way of writing a palindrome as a product of two $n$-digit numbers gives some $u$, $v$, and $c \ge 0$ that satisfy both equations, and for that carry $D = (u - v)^2 \ge 0$, so it gets tested before the loop stops. And since the palindromes get smaller as $m$ grows, the first one that passes is the largest.

This is Fermat's factorization method in disguise. Fermat factors a number $N$ by looking for an $a$ that makes $a^2 - N$ a perfect square $b^2$, because then $N = (a - b)(a + b)$. It's fast when the two factors are close together. Here $x + y = 2B - (m + c)$, so our discriminant is $D = (x + y)^2 - 4xy = (x - y)^2$, which is the same test with $a = (x + y)/2$. The digits of the palindrome pin down $x + y$ to one candidate per carry, and both factors are within about $10^{n/2}$ of $10^n$.

Benchmarking largest_palindrome_product_fermat(8) we find the answer in 63.725 μs, which is 10,100× faster than the pruned search. For the original 3-digit problem it takes 613.451 ns, which is 26.0× faster. The 12-digit case, which would take many hours by searching products runs in 9.242 ms, and the 15-digit case runs in 301.083 ms.

Skipping palindromes that can't be products

The Fermat search reverses the upper half of every palindrome and then computes a discriminant and an integer square root for every carry it tries. But most palindromes can be ruled out by looking at just a few of their digits, before doing any of that work, and the ones that are left can be checked more cheaply.

Let $s = u + v = m + c$. Near the top, the upper half starts with $999$, so the palindrome ends in $999$. The lower half is $uv - cB$ and $B = 10^n$ is a multiple of $1000$, so $uv$ ends in $999$ as well. That tells us a lot about $s$.

First, $uv$ ends in $9$, so $u$ and $v$ end in $1$ and $9$, $3$ and $3$, or $7$ and $7$. So $s$ ends in $0$, $6$ or $4$. Second, $uv \equiv 999 \equiv 7 \pmod 8$, so $u$ and $v$ are both odd. The residue of a number modulo $k$ is its remainder after dividing by $k$. The only odd residues modulo $8$ that multiply to $7$ are $1 \times 7$ and $3 \times 5$. Both add up to $8$, so $s \equiv 0 \pmod 8$.

These rules are about $s$ rather than $m$, so we'll loop over $s$ instead and set $m = s - c$ for each carry.

Both rules only depend on $s \bmod 8$ and $s \bmod 10$, which repeat every $40$ sums since $40$ is the smallest number that $8$ and $10$ both divide. Of $0, 1, \dots, 39$, only $0$, $16$ and $24$ pass, so of any $40$ consecutive sums only the ones congruent to those modulo $40$ can work. So rather than testing every $s$, we can step through the sums in blocks of $40$ and only try those three, skipping the other $37$ without doing any work.

In blocks of $120$ the same rules allow these residues:

const SUM_RESIDUES = (0, 16, 24, 40, 56, 64, 80, 96, 104)

Blocks of $40$ would keep the same $3$ in $40$ sums, but since $40 \equiv 1 \pmod 3$, the same residue would give a different $s \bmod 3$ in each of three consecutive blocks. In blocks of $120 = 8 \times 3 \times 5$ each residue also fixes $s \bmod 3$, so the next rule can drop whole residues instead of checking every sum.

A number is congruent to its digit sum modulo $3$, and the upper and lower halves have the same digits, so $L \equiv U \pmod 3$. Since $B = 10^n \equiv 1 \pmod 3$ and $U = B - m = B - s + c$, the discriminant becomes

A perfect square is always $0$ or $1$ modulo $3$, so we can skip any $s$ and $c$ that make this $2$. Here r is $s \bmod 120$:

passes_mod_3_rule(r, c) = mod(r^2 + r + c - 1, 3) != 2

Without a carry that only leaves $s \equiv 1 \pmod 3$, i.e. $s \equiv 16$, $40$ or $64 \pmod{120}$. And there can't be a carry while $s^2 < 4B$, since $cB \le uv \le s^2/4$. So near the top only $1$ in $40$ values of $s$ needs checking.

The middle of the palindrome narrows it down even further. The last digit of $U$ is $d = (c - s) \bmod 10$, and it's also the first digit of $L$, so $L \ge d \times 10^{n-1}$. For $D \ge 0$ we need $s^2 \ge 4(cB + L)$, so we can check $s^2 \ge 4(cB + d \times 10^{n-1})$ before working out $L$ at all:

function passes_middle_digit_rule(s_squared, r, c, search)
    d = mod(c - r, 10)
    return 4 * (c * search.B + d * search.leading_place) <= s_squared
end

The palindromes that get through can also be checked more cheaply:

• Near the top the upper half is a run of $9$s followed by a short tail, so the lower half is the reversed tail followed by $9$s, and only the tail needs reversing. Reversing takes a division per digit, and once the numbers need Int128 those divisions are slow, but the tail always fits in an Int64.
• Only $12$ of the $64$ residues modulo $64$ are perfect squares, so checking $D \bmod 64$ against a 64-bit mask first rules out about $81%$ of the non-squares before taking the square root.
• Julia turns s^2 into s * s for integers up to 64 bits, but for Int128 it calls the generic power_by_squaring, so we write s * s ourselves.

function lower_half(m, search)
    reversed_tail = reverse_tail(search.ten_to_tail - Int64(m), search.tail_digits)
    return reversed_tail * search.nines_shift + (search.nines_shift - 1)
end

function reverse_tail(t::Int64, k)
    reversed = 0
    for _ in 1:k
        t, digit = divrem(t, 10)
        reversed = 10reversed + digit
    end
    return reversed
end

const SQUARES_MOD_64 = reduce(|, UInt64(1) << (i^2 % 64) for i in 0:63)

function perfect_square_root(x)
    (x >= 0 && isodd(SQUARES_MOD_64 >> ((x % UInt8) & 0x3f))) || return nothing
    root = isqrt(x)
    return root * root == x ? root : nothing
end

Looping over $s$ has one catch: we no longer visit the palindromes in order, since a larger $s$ with a larger carry can still give a smaller $m$. So instead of stopping at the first hit, we keep the smallest $m$ found so far. The carry is at most $s^2/4B$, so every $s$ from here on gives $m \ge s - \lfloor s^2/4B \rfloor$, and once that's larger than the best $m$ we can stop:

no_smaller_m_from(s, m_best, search) = s - fld(s * s, 4 * search.B) > m_best

Putting it all together:

struct FilteredSearch{W}
    B::W                # 10^n
    leading_place::W    # 10^(n-1), the place value of the lower half's first digit
    tail_digits::Int    # k: while m ≤ 10^k the upper half is n - k nines followed by a k-digit tail
    ten_to_tail::Int64  # 10^k
    nines_shift::W      # 10^(n-k), which shifts the reversed tail past the n - k nines that end the lower half
    max_sum::W          # the largest s to check, since the residue rules need m ≤ 10^(n-3) and the tail needs m ≤ 10^k
end

function FilteredSearch(n, ::Type{W}) where W
    k = min(n, cld(n, 2) + 3, 18)  # deeper than any answer, and short enough for an Int64
    max_sum = n >= 3 ? min(W(10)^(n - 3), W(10)^k) : zero(W)
    return FilteredSearch(W(10)^n, W(10)^(n - 1), k, 10^k, W(10)^(n - k), max_sum)
end

function largest_palindrome_product_fermat_filtered(n, ::Type{T}=integer_type(n)) where T
    n >= 2 || throw(ArgumentError("expected n ≥ 2 digits, but got n = $n"))

    search = FilteredSearch(n, promote_type(T, Int128))
    best = find_smallest_m(search)

    # Fall back to the plain search unless no sum past max_sum can beat the best m, which only happens for n ≤ 5
    if isnothing(best) || !no_smaller_m_from(search.max_sum + 1, best.m, search)
        return largest_palindrome_product_fermat(n, T)
    end

    B = T(10)^n
    U = B - T(best.m)
    return (palindrome=widemul(U, B) + reverse_digits(U), factors=(B - T(best.v), B - T(best.u)))
end

function find_smallest_m(search::FilteredSearch{W}) where {W}
    best = nothing
    for block_start in zero(W):W(120):search.max_sum
        !isnothing(best) && no_smaller_m_from(block_start, best.m, search) && break
        block_end = block_start + 119
        no_carry = block_end * block_end < 4 * search.B  # cB ≤ uv ≤ s²/4

        for r in SUM_RESIDUES
            no_carry && r % 3 != 1 && continue  # with c = 0 the mod 3 rule needs s ≡ 1 (mod 3)
            s = block_start + r
            s <= search.max_sum || break

            hit = smallest_m_for_sum(s, r, search)
            if !isnothing(hit) && (isnothing(best) || hit.m < best.m)
                best = hit
            end
        end
    end
    return best
end

function smallest_m_for_sum(s::W, r, search::FilteredSearch{W}) where {W}
    s_squared = s * s  # not s^2, which isn't a plain multiplication for Int128
    best = nothing
    c = 0
    while 4 * c * search.B <= s_squared  # cB ≤ uv ≤ s²/4
        if passes_mod_3_rule(r, c) && passes_middle_digit_rule(s_squared, r, c, search)
            m = s - c
            root = perfect_square_root(s_squared - 4 * (c * search.B + lower_half(m, search)))
            if !isnothing(root)
                best = (m=m, u=(s - root) ÷ 2, v=(s + root) ÷ 2)  # a larger carry gives a smaller m
            end
        end
        c += 1
    end
    return best
end

The rules only hold while the upper half starts with $999$, i.e. for $s \le 10^{n-3}$, and the tail has $k = \min(n, \lceil n/2 \rceil + 3, 18)$ digits, which is deeper than any answer has been. If the search can't prove its best palindrome is the largest within those limits, it falls back to the plain Fermat search. That only happens for 5 digits or fewer, where the plain search is fast anyway.

Benchmarking largest_palindrome_product_fermat_filtered(12) we find the same answer in 460.289 μs, which is 20.1× faster than the plain Fermat search, and the 15-digit case takes 11.801 ms, which is 25.5× faster.

Running the search on a GPU

The filtered search still checks one sum at a time. The answer is usually around $10^{n/2}$ palindromes from the top, so each extra digit means about $\sqrt{10} \approx 3$ times as many sums to check. But checking one sum doesn't depend on any other, and all we need from them is the smallest $m$, so we can check many sums at the same time. It's an embarassingly parallel problem. That's what GPUs are good at!

Each GPU thread can take groups of $120$ consecutive sums, like the blocks in find_smallest_m, and checks the ones whose residues pass the rules. When a thread finds a palindrome product, it lowers the smallest $m$ found so far with CUDA.@atomic, so that two threads can't overwrite each other's results:

function search_kernel!(best_m, first_sum::UInt64, num_groups::UInt64, search::GPUSearch)
    thread_index = UInt64((blockIdx().x - 1) * blockDim().x + threadIdx().x - 1)  # from 0
    num_threads = UInt64(gridDim().x) * UInt64(blockDim().x)
    group = thread_index
    while group < num_groups
        check_group!(best_m, first_sum + 120group, search)
        group += num_threads
    end
    return nothing
end

@inline function check_group!(best_m, group_start::UInt64, search::GPUSearch)
    group_end = group_start + 119
    if group_end < search.only_40_below
        check_sum!(best_m, group_start + 40, 40, search)
    elseif group_end < search.no_carry_below
        check_sum!(best_m, group_start + 16, 16, search)
        check_sum!(best_m, group_start + 40, 40, search)
        check_sum!(best_m, group_start + 64, 64, search)
    else
        for residue in SUM_RESIDUES
            check_sum!(best_m, group_start + residue, residue, search)
        end
    end
    return nothing
end

@inline function check_sum!(best_m, s::UInt64, residue::Int, search::GPUSearch)
    s <= search.max_sum || return nothing  # only the last group can run past max_sum
    s_squared = widemul(s, s)
    four_B = 4 * search.B
    c = 0
    four_cB = UInt128(0)
    while four_cB <= s_squared  # cB ≤ uv ≤ s²/4
        if passes_mod_3_rule(residue, c) && passes_middle_digit_rule_gpu(s_squared, four_cB, residue, c, search)
            m = s - UInt64(c)
            four_uv = four_cB + 4 * lower_half_gpu(m, search)
            if four_uv <= s_squared && is_perfect_square(s_squared - four_uv)
                CUDA.@atomic best_m[1] = min(best_m[1], m)
            end
        end
        c += 1
        four_cB += four_B
    end
    return nothing
end

@inline function passes_middle_digit_rule_gpu(s_squared, four_cB, residue, c, search::GPUSearch)
    d = mod(c - residue, 10)
    return four_cB + UInt128(4d) * search.leading_digit_place <= s_squared
end

Each thread loops over every num_threads-th group, so the kernel can be launched on any number of groups. Rather than testing every residue against the rules, check_group! uses two thresholds: below the first only $s \equiv 40$ can pass, and below the second there's no carry yet, so only $16$, $40$ and $64$ can. check_sum! also keeps $4cB$ as a running total instead of multiplying for every carry, which is the only change in passes_middle_digit_rule_gpu.

The CPU decides how far to search. It launches the kernel on batches of groups from $s = 0$ upwards, doubling the batch each time, so an answer near the top only takes a few small launches while a deep one gets launches big enough to keep the GPU busy. After each launch it copies back the smallest $m$ and stops once no larger $s$ can beat it, using no_smaller_m_from again:

function largest_palindrome_product_gpu(n)
    MIN_GPU_DIGITS <= n <= MAX_GPU_DIGITS ||
        throw(ArgumentError("expected $MIN_GPU_DIGITS ≤ n ≤ $MAX_GPU_DIGITS, but got n = $n"))

    search = GPUSearch(n)
    no_hit = typemax(UInt64)  # the smallest m until the GPU finds a palindrome product
    best_m = CUDA.fill(no_hit, 1)
    num_blocks = BLOCKS_PER_MULTIPROCESSOR * CUDA.attribute(CUDA.device(), CUDA.DEVICE_ATTRIBUTE_MULTIPROCESSOR_COUNT)

    # Launch batches of groups from s = 0 upwards, doubling the batch size so that answers near the top are found
    # quickly. Groups start at multiples of 120, and the last one can run past max_sum.
    first_sum = UInt64(0)
    groups_per_launch = UInt64(FIRST_LAUNCH_GROUPS)
    while first_sum <= search.max_sum
        num_groups = min(groups_per_launch, (search.max_sum - first_sum) ÷ 120 + 1)
        @cuda threads=THREADS_PER_BLOCK blocks=num_blocks search_kernel!(best_m, first_sum, num_groups, search)
        first_sum += 120 * num_groups
        next_sum = min(first_sum, search.max_sum + 1)  # the first sum that hasn't been checked

        m = Array(best_m)[1]  # waits for the kernel to finish
        if m != no_hit && no_smaller_m_from(UInt128(next_sum), m, search)  # UInt128 so that s² doesn't overflow
            return palindrome_and_factors(n, m)
        end
        groups_per_launch = min(2groups_per_launch, UInt64(MAX_LAUNCH_GROUPS))
    end

    error("found no palindrome product with u + v ≤ $(search.max_sum) for n = $n")
end

The GPU only reports $m$, so palindrome_and_factors then works out the carry and the two factors on the CPU. There's no fallback to the plain search, so the GPU version needs at least 6 digits.

The arithmetic takes more care on a GPU. Its integer instructions work on 32 bits, so 64-bit arithmetic takes a few instructions and 128-bit arithmetic many more, and division, which the hardware can't do directly, is the slowest of all.

Base.@kwdef struct GPUSearch
    B::UInt128                    # 10^n
    max_sum::UInt64               # min(10^18, 10^(n-3)), the largest s to check
    leading_digit_place::UInt128  # 10^(n-1), the place value of L's first digit
    tail_digits::Int              # k = min(n, 18)
    ten_to_tail_digits::UInt64    # 10^k
    nines_shift::UInt64           # 10^(n-k), which shifts the reversed tail past the n - k nines that end L
    only_40_below::UInt64         # groups ending below this can only pass with s ≡ 40 (mod 120) and no carry
    no_carry_below::UInt64        # groups ending below this have s² < 4B, so no carry
end

function GPUSearch(n::Integer)
    B = big(10)^n
    k = min(n, 18)
    ceil_sqrt(x) = isqrt(x - 1) + 1

    # The search stops at s = min(10^18, 10^(n-3)). The 10^18 keeps s in a UInt64 and s² in a UInt128, and it's deep
    # enough for every n that's feasible. The 10^(n-3) keeps m ≤ 10^(n-3), since m ≤ s, so the upper half starts with
    # 999 as the residue rules need. Both factors are then within 10^(n-3) of B, so they always have n digits.
    return GPUSearch(
        B=UInt128(B),
        max_sum=UInt64(min(big(10)^18, big(10)^(n - 3))),
        leading_digit_place=UInt128(10)^(n - 1),
        tail_digits=k,
        ten_to_tail_digits=UInt64(10)^k,
        nines_shift=UInt64(10)^(n - k),
        only_40_below=UInt64(ceil_sqrt(16 * big(10)^(n - 1))),  # s ≡ 16 has d = 4, so the middle digit rule needs this
        no_carry_below=UInt64(ceil_sqrt(4B)))                   # a carry needs s² ≥ 4B
end

Two pieces of the CPU code are also rewritten to avoid slow arithmetic. Reversing the tail takes a division by $10$ per digit, so reverse_tail_gpu splits the tail into its last $9$ digits and the rest, which both fit in 32 bits. That makes the whole search about 1.5× faster than dividing a UInt64. And Base's isqrt corrects its floating-point estimate with a Newton step that divides one UInt128 by another, which makes the whole search about 6× slower. So isqrt_uint128 corrects the estimate by comparing squares instead, which only needs multiplications:

@inline function lower_half_gpu(m::UInt64, search::GPUSearch)
    tail = search.ten_to_tail_digits - m
    # From n = 18 on every tail has 18 digits. Passing that as a constant lets the compiler unroll the digit loops,
    # which makes large n about 10% faster.
    reversed_tail = search.tail_digits == 18 ? reverse_tail_gpu(tail, 18) : reverse_tail_gpu(tail, search.tail_digits)
    return widemul(reversed_tail, search.nines_shift) + (search.nines_shift - 1)
end

@inline function reverse_tail_gpu(x::UInt64, num_digits::Int)
    high, low = divrem(x, UInt64(10)^9)
    reversed_low, _ = reverse_last_digits(low % UInt32, min(num_digits, 9))
    reversed_high, high_place = reverse_last_digits(high % UInt32, num_digits - 9)
    return UInt64(reversed_low) * high_place + reversed_high
end

@inline function reverse_last_digits(x::UInt32, num_digits::Int)
    reversed = UInt32(0)
    place = UInt32(1)
    for _ in 1:num_digits
        x, digit = divrem(x, UInt32(10))
        reversed = UInt32(10) * reversed + digit
        place *= UInt32(10)
    end
    return reversed, place
end

@inline function is_perfect_square(x::UInt128)
    isodd(SQUARES_MOD_64 >> ((x % UInt64) & 63)) || return false
    root = isqrt_uint128(x)
    return widemul(root, root) == x
end

@inline function isqrt_uint128(x::UInt128)
    estimate = Float64((x >> 64) % UInt64) * 2.0^64 + Float64(x % UInt64)
    root = unsafe_trunc(UInt64, sqrt(estimate))
    while widemul(root, root) > x
        root -= 1
    end
    while widemul(root + 1, root + 1) <= x
        root += 1
    end
    return root
end

The residue rules don't change, so the GPU code reuses SUM_RESIDUES, passes_mod_3_rule and SQUARES_MOD_64 from the CPU version.

On a GPU the 20-digit case takes 1.236 ms, which is 895× faster than the filtered search on one CPU core, and the 24-digit case takes 245.480 ms.

That's fast enough to keep going past 24 digits!

The UInt128 arithmetic works up to 37 digits, where $4B$ stops fitting, but each digit makes the search about $\sqrt{10}$ times longer. With an even number of digits the answer is at most $10^{n/2}$ palindromes from the top, for a reason we'll see below, so 32 digits takes up to about an hour, while an odd number of digits can take longer, depending on how far down the answer is. The 33- to 36-digit cases took between 4 and 78 hours each. That's about as far as this code goes, since it only searches up to $s = 10^{18}$, and the 37-digit answer is probably further down than that.

Every answer up to 36 digits

Here's the largest palindrome that's a product of two $n$-digit numbers for every $n$ from $2$ to $36$. These are OEIS A327897, and each one has just one factorization into two $n$-digit numbers:

OEIS A308803 asks the same question by the number of digits in the palindrome instead: the largest $k$-digit palindrome that's a product of two numbers with the same number of digits. For even $k$ that's the table above, since a $2n$-digit product needs two $n$-digit factors. For odd $k$ the two factors can be much further apart, which makes those terms easy to find, and they're known up to $95$ digits.

Why so many answers are all nines and zeros

Most of the even-$n$ answers in the table are $n/2$ nines, then $n$ zeros, then $n/2$ nines again. That's because for every even $n$ there's a palindrome product of that shape:

One factor is all nines, and the other is $n/2$ nines followed by zeros and a final $1$. In the notation from earlier, $u = 1$ and $v = 10^{n/2} - 1$, so $u + v = 10^{n/2}$ and $uv = 10^{n/2} - 1$, which is exactly the upper half $10^n - 10^{n/2}$ reversed. It's the palindrome $m = 10^{n/2}$ from the top, so with an even number of digits the answer is never further down than that. The answer is this palindrome unless one of the $10^{n/2} - 1$ palindromes above it is also a product.

How far each method gets

Here's every benchmark from this post in one table. Each method is only benchmarked at sizes it can do in a few seconds, and it shares at least two sizes with the method before it. Under each time is how many times faster that method is than the one to its left. The CPU searches all run on one core.

Searching palindromes instead of products changes how quickly the time grows. From 6 to 8 digits, the pruned search gets 140× slower, but the Fermat search only gets 13.1× slower, since the answer is only around $10^{n/2}$ palindromes from the top. So the palindrome searches only get about $\sqrt{10} \approx 3$ times slower with each extra digit, which makes every later speedup worth a few more digits in the same time: the filters' 25.5× is worth about 3, and the GPU's 895× about 6.

The filters skip about the same share of the work at every size, so they stay around 20.1× to 25.5× faster. The GPU instead pulls further ahead as $n$ grows, since the fixed cost of each kernel launch matters less, going from 84.2× faster than one CPU core at 15 digits to 895× at 20.