Project Euler Problem 954

Solution

Answer: 736463823

Let a number $n$ have digits $d_i$ at decimal positions $i$ (least significant digit = position $0$). Suppose we swap digits $d_i=a$ and $d_j=b$. The new number becomes

$$n' = n + (b-a)(10^i - 10^j).$$

We want $n'$ to never be divisible by $7$, while $n$ itself is also not divisible by $7$.

Working modulo $7$, since $10^k \bmod 7$ is periodic with period $6$,

$$10^k \equiv (1,3,2,6,4,5,\ldots)\pmod 7.$$

For a fixed target remainder $r=n\bmod 7\neq 0$, a swap between positions with residues $x=10^i\bmod 7$ and $y=10^j\bmod 7$ is forbidden iff

$$r + (b-a)(x-y)\equiv 0 \pmod 7,$$

equivalently

$$b \equiv a - r(y-x)^{-1}\pmod 7.$$

Thus, for each final remainder $r\in{1,\dots,6}$, we can build a digit-DP:

process digits left to right,
track current remainder mod $7$,
track which digit residues mod $7$ have already appeared in each of the six position residue classes,
reject any digit that would create a swap yielding divisibility by $7$,
enforce the “no leading zero after swap” rule.

This DP reproduces the checks:

$C(100)=74$
$C(10^4)=3737$

as required.

Clean Python implementation:

from functools import lru_cache

mods = [pow(10, i, 7) for i in range(13)]

def count_length(L):
    posmods = list(reversed(mods[:L]))  # MSD -> LSD
    total = 0

    for r in range(1, 7):

        # forbidden digit residues lookup
        block = [[None] * 7 for _ in range(7)]

        for x in range(1, 7):
            for y in range(1, 7):
                arr = [0] * 128

                diff = (y - x) % 7
                single = [0] * 7

                if diff:
                    inv = pow(diff, -1, 7)
                    k = (-r * inv) % 7
                    for m in range(7):
                        single[m] = 1 << ((m + k) % 7)

                for mask in range(128):
                    b = 0
                    mm = mask
                    while mm:
                        lsb = mm & -mm
                        m = lsb.bit_length() - 1
                        b |= single[m]
                        mm ^= lsb
                    arr[mask] = b

                block[x][y] = arr

        x0 = posmods[0]

        @lru_cache(None)
        def dp(i, rem, masks, msd_res):
            if i == L:
                return rem == r

            x = posmods[i]

            forbidden = 0
            tmp = masks

            for y in range(1, 7):
                forbidden |= block[x][y][tmp & 127]
                tmp >>= 7

            msd_forbidden = 0
            if i > 0 and (x0 - x) % 7:
                diff = (x0 - x) % 7
                inv = pow(diff, -1, 7)
                k = (-r * inv) % 7
                msd_forbidden = 1 << ((msd_res + k) % 7)

            ans = 0

            for d in range(10):
                if i == 0 and d == 0:
                    continue

                f = forbidden
                if i > 0 and d != 0:
                    f |= msd_forbidden

                dr = d % 7
                if (f >> dr) & 1:
                    continue

                nmasks = masks
                if i > 0:
                    shift = (x - 1) * 7
                    nmasks |= (1 << dr) << shift

                ans += dp(i + 1, (rem * 10 + d) % 7,
                          nmasks, msd_res)

            return ans

        for first_digit in range(1, 10):
            total += dp(1, first_digit % 7, 0,
                        first_digit % 7)

    return total

answer = sum(count_length(L) for L in range(1, 14))
print(answer)

The computation gives:

$$C(10^{13}) = 736463823.$$

Answer: 736463823