Problem 612: Friend Numbers

View on Project Euler

Project Euler Problem 612 Solution

EulerSolve provides an optimized solution for Project Euler Problem 612, Friend Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let \(f(10^K)\) denote the number of pairs \((p,q)\) with \(1 \le p \lt q \lt 10^K\) such that the decimal representations of \(p\) and \(q\) share at least one digit. For Problem 612 we need $$f(10^{18}) \bmod 1000267129.$$ A direct scan is impossible: there are \(10^{18}-1\) candidate numbers and about \(\binom{10^{18}-1}{2}\) pairs. The workable idea is to forget the exact values and remember only which decimal digits occur in each number. Mathematical Approach The solution counts numbers by their exact digit sets, then subtracts the pairs whose digit sets are disjoint. Step 1: Encode Each Number by Its Digit Set For a positive integer \(n\), let \(S(n)\subseteq\{0,1,\dots,9\}\) be the set of digits appearing in its usual decimal expansion. Two numbers are friends exactly when $$S(p)\cap S(q)\neq \varnothing.$$ For each mask \(A\subseteq\{0,1,\dots,9\}\), define \(C_K(A)\) as the number of integers \(n\) with \(1\le n\lt 10^K\) and \(S(n)=A\). Once all values \(C_K(A)\) are known, the original problem becomes a counting problem over only \(2^{10}=1024\) masks. Step 2: Count Exact Masks with a Length DP Let \(D_\ell(A)\) be the number of \(\ell\)-digit positive integers whose digit set is exactly \(A\)....

Detailed mathematical approach

Problem Summary

Let \(f(10^K)\) denote the number of pairs \((p,q)\) with \(1 \le p \lt q \lt 10^K\) such that the decimal representations of \(p\) and \(q\) share at least one digit. For Problem 612 we need

$$f(10^{18}) \bmod 1000267129.$$

A direct scan is impossible: there are \(10^{18}-1\) candidate numbers and about \(\binom{10^{18}-1}{2}\) pairs. The workable idea is to forget the exact values and remember only which decimal digits occur in each number.

Mathematical Approach

The solution counts numbers by their exact digit sets, then subtracts the pairs whose digit sets are disjoint.

Step 1: Encode Each Number by Its Digit Set

For a positive integer \(n\), let \(S(n)\subseteq\{0,1,\dots,9\}\) be the set of digits appearing in its usual decimal expansion. Two numbers are friends exactly when

$$S(p)\cap S(q)\neq \varnothing.$$

For each mask \(A\subseteq\{0,1,\dots,9\}\), define \(C_K(A)\) as the number of integers \(n\) with \(1\le n\lt 10^K\) and \(S(n)=A\).

Once all values \(C_K(A)\) are known, the original problem becomes a counting problem over only \(2^{10}=1024\) masks.

Step 2: Count Exact Masks with a Length DP

Let \(D_\ell(A)\) be the number of \(\ell\)-digit positive integers whose digit set is exactly \(A\). Because a decimal representation cannot start with zero, the base layer is

$$D_1(\{d\})=1 \quad \text{for } d\in\{1,\dots,9\},$$

and all other one-digit states are zero.

If an \(\ell\)-digit number currently uses digit set \(A\), then appending a digit \(x\in\{0,\dots,9\}\) produces an \((\ell+1)\)-digit number with digit set \(A\cup\{x\}\). So the transition is simply

$$A \longrightarrow A\cup\{x\}\qquad (x=0,1,\dots,9).$$

Summing over all permitted lengths gives

$$C_K(A)=\sum_{\ell=1}^{K} D_\ell(A).$$

This matches the range \(1\le n\lt 10^K\): every such number has between \(1\) and \(K\) digits, and leading zeros are never introduced.

Step 3: Count the Complement First

Let

$$M=10^K-1.$$

Then the total number of unordered pairs is

$$\binom{M}{2}.$$

A pair is not friendly precisely when the two digit sets are disjoint. If \(N_K\) denotes the number of unordered pairs with

$$A\cap B=\varnothing,$$

then the desired answer is

$$f(10^K)=\binom{M}{2}-N_K.$$

So the task is reduced to counting disjoint mask pairs.

Step 4: Use Subset Sums to Count Disjoint Masks Quickly

For every mask \(T\subseteq\{0,\dots,9\}\), define

$$Z_K(T)=\sum_{B\subseteq T} C_K(B).$$

Now fix a mask \(A\). A second number is disjoint from it exactly when its mask \(B\) is contained in the complement

$$\overline{A}=\{0,1,\dots,9\}\setminus A.$$

Therefore the number of choices for the second number is \(Z_K(\overline{A})\), and the number of ordered disjoint pairs is

$$N_K^{\mathrm{ord}}=\sum_{A} C_K(A)\,Z_K(\overline{A}).$$

The empty mask never occurs for a positive integer, so a disjoint pair can never be a self-pair. Hence every unordered disjoint pair appears exactly twice in the ordered count, which gives

$$N_K=\frac{N_K^{\mathrm{ord}}}{2}.$$

The values \(Z_K(T)\) are computed by the standard subset-zeta transform in \(O(10\cdot 2^{10})\) time.

Step 5: Worked Example for \(K=2\)

Now the numbers run from \(1\) to \(99\), so \(M=99\) and

$$\binom{99}{2}=4851.$$

The exact mask counts are easy to list:

\(\{d\}\) with \(d\in\{1,\dots,9\}\): two numbers, namely \(d\) and \(dd\).

\(\{0,d\}\): one number, namely \(d0\).

\(\{a,b\}\) with \(1\le a\lt b\le 9\): two numbers, namely \(ab\) and \(ba\).

There are no other masks below \(100\).

Now count disjoint unordered pairs by type:

$$\begin{aligned} \{a\},\{b\}:&\quad \binom{9}{2}\cdot 2\cdot 2 = 144,\\ \{a\},\{0,b\},\ a\neq b:&\quad 9\cdot 8\cdot 2\cdot 1 = 144,\\ \{a\},\{b,c\},\ a\notin\{b,c\}:&\quad \binom{9}{2}\cdot 7\cdot 2\cdot 2 = 1008,\\ \{0,a\},\{b,c\},\ a\notin\{b,c\}:&\quad \binom{9}{2}\cdot 7\cdot 1\cdot 2 = 504,\\ \{a,b\},\{c,d\},\ |\{a,b,c,d\}|=4:&\quad \binom{9}{4}\cdot 3\cdot 2\cdot 2 = 1512. \end{aligned}$$

Therefore

$$N_2=144+144+1008+504+1512=3312,$$

and the number of friend pairs is

$$f(10^2)=4851-3312=1539,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations keep small arrays indexed by the 1024 digit masks. One array stores the current length layer, another stores cumulative counts over all lengths up to \(K\). The process starts from the nine one-digit states, repeatedly appends each decimal digit \(0\) through \(9\), and adds every finished layer into the totals.

Next, the implementation copies the exact-mask counts into another 1024-entry array and performs ten subset-zeta passes, one pass per digit. After that transform, each mask stores the sum of all counts over its submasks, so the number of masks disjoint from a given mask can be read directly from its complement.

Finally, the implementation accumulates all ordered disjoint pairs, multiplies by the modular inverse of \(2\) to obtain unordered disjoint pairs, computes \(\binom{10^K-1}{2}\) modulo \(1000267129\), and subtracts. For the actual Project Euler input it sets \(K=18\).

Complexity Analysis

The length DP uses \(K\) layers, \(2^{10}\) masks, and \(10\) appended digits, so it costs \(O(K\cdot 10 \cdot 2^{10})\) time. The subset-zeta transform costs \(O(10\cdot 2^{10})\) time. Memory usage is \(O(2^{10})\), because only a few arrays of length \(1024\) are stored.

With decimal digits fixed, the whole method is tiny. For \(K=18\), the runtime is effectively constant relative to the enormous size of the original search space.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=612
  2. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
  3. Dynamic programming: Wikipedia — Dynamic programming
  4. Bit array / bitmask representation: Wikipedia — Bit array
  5. Subset iteration and SOS-style preprocessing: CP-Algorithms — Enumerating submasks of a bitmask

Problem 612 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

// Project Euler 612: count pairs of integers < 10^K whose decimal digit-sets intersect.

using i64 = long long;
using i128 = __int128_t;

static constexpr i64 MOD = 1000267129LL;
static constexpr int DIGS = 10;
static constexpr int FULL = (1 << DIGS) - 1;

static inline i64 mod_norm(i64 x) {
    x %= MOD;
    if (x < 0) x += MOD;
    return x;
}

static inline i64 mul_mod(i64 a, i64 b) { return (i64)((i128)a * b % MOD); }

static i64 pow_mod(i64 a, std::uint64_t e) {
    i64 r = 1 % MOD;
    a %= MOD;
    while (e) {
        if (e & 1) r = mul_mod(r, a);
        a = mul_mod(a, a);
        e >>= 1;
    }
    return r;
}

static std::vector<i64> count_by_mask(int max_len) {
    std::vector<i64> cnt(1 << DIGS, 0);
    std::vector<i64> dp(1 << DIGS, 0), ndp(1 << DIGS, 0);

    for (int d = 1; d <= 9; ++d) {
        dp[1 << d] += 1;
        if (dp[1 << d] >= MOD) dp[1 << d] -= MOD;
    }

    for (int len = 1; len <= max_len; ++len) {
        for (int m = 0; m <= FULL; ++m) {
            cnt[m] += dp[m];
            if (cnt[m] >= MOD) cnt[m] -= MOD;
        }
        if (len == max_len) break;

        std::fill(ndp.begin(), ndp.end(), 0);
        for (int m = 0; m <= FULL; ++m) {
            const i64 v = dp[m];
            if (!v) continue;
            for (int d = 0; d <= 9; ++d) {
                const int nm = m | (1 << d);
                ndp[nm] += v;
                if (ndp[nm] >= MOD) ndp[nm] -= MOD;
            }
        }
        dp.swap(ndp);
    }
    return cnt;
}

static i64 f_power10(int K) {
    const i64 inv2 = (MOD + 1) / 2;

    const std::vector<i64> cnt = count_by_mask(K);

    std::vector<i64> sum_sub = cnt; // SOS DP: sum_sub[mask] = sum_{s subset mask} cnt[s]
    for (int b = 0; b < DIGS; ++b) {
        for (int mask = 0; mask <= FULL; ++mask) {
            if (mask & (1 << b)) {
                sum_sub[mask] += sum_sub[mask ^ (1 << b)];
                if (sum_sub[mask] >= MOD) sum_sub[mask] -= MOD;
            }
        }
    }

    i64 ordered_disjoint = 0; // ordered pairs (p,q) with disjoint digit-sets
    for (int a = 0; a <= FULL; ++a) {
        ordered_disjoint += mul_mod(cnt[a], sum_sub[FULL ^ a]);
        ordered_disjoint %= MOD;
    }
    const i64 disjoint_pairs = mul_mod(ordered_disjoint, inv2);

    const i64 M = mod_norm(pow_mod(10, (std::uint64_t)K) - 1); // count of integers in [1,10^K-1]
    const i64 total_pairs = mul_mod(mul_mod(M, mod_norm(M - 1)), inv2);

    return mod_norm(total_pairs - disjoint_pairs);
}

static int digit_mask_u64(std::uint64_t x) {
    int m = 0;
    while (x) {
        m |= 1 << (int)(x % 10);
        x /= 10;
    }
    return m;
}

static std::uint64_t brute_f(int n) {
    std::vector<int> masks((std::size_t)n, 0);
    for (int i = 1; i < n; ++i) masks[i] = digit_mask_u64((std::uint64_t)i);

    std::uint64_t ans = 0;
    for (int p = 1; p < n; ++p) {
        for (int q = p + 1; q < n; ++q) {
            if (masks[p] & masks[q]) ++ans;
        }
    }
    return ans;
}

int main() {
    assert(f_power10(2) == 1539);
    assert(f_power10(3) == (i64)(brute_f(1000) % (std::uint64_t)MOD));

    std::cout << f_power10(18) << "\n";
    return 0;
}

Python

MOD = 1000267129
DIGS = 10
FULL = (1 << DIGS) - 1

def count_by_mask(max_len):
    cnt = [0] * (FULL + 1)
    dp = [0] * (FULL + 1)
    ndp = [0] * (FULL + 1)
    
    for d in range(1, 10):
        dp[1 << d] += 1
        
    for length in range(1, max_len + 1):
        for m in range(FULL + 1):
            cnt[m] = (cnt[m] + dp[m]) % MOD
            
        if length == max_len:
            break
            
        ndp = [0] * (FULL + 1)
        for m in range(FULL + 1):
            v = dp[m]
            if not v: continue
            for d in range(10):
                nm = m | (1 << d)
                ndp[nm] = (ndp[nm] + v) % MOD
                
        dp = ndp
        
    return cnt

def solve_f(K):
    inv2 = (MOD + 1) // 2
    cnt = count_by_mask(K)
    
    sum_sub = list(cnt)
    for b in range(DIGS):
        bit = 1 << b
        for mask in range(FULL + 1):
            if mask & bit:
                sum_sub[mask] = (sum_sub[mask] + sum_sub[mask ^ bit]) % MOD
                
    ordered_disjoint = 0
    for a in range(FULL + 1):
        ordered_disjoint = (ordered_disjoint + cnt[a] * sum_sub[FULL ^ a]) % MOD
        
    disjoint_pairs = (ordered_disjoint * inv2) % MOD
    
    M = (pow(10, K, MOD) - 1) % MOD
    total_pairs = (M * (M - 1) % MOD) * inv2 % MOD
    
    return (total_pairs - disjoint_pairs + MOD) % MOD

def solve():
    return str(solve_f(18))

if __name__ == '__main__':
    print(solve())

Java

public class Euler612 {
    static final long MOD = 1000267129L;
    static final int DIGS = 10;
    static final int FULL = (1 << DIGS) - 1;

    static long modPow(long a, long e) {
        long r = 1;
        a %= MOD;
        while (e > 0) {
            if ((e & 1) == 1)
                r = (r * a) % MOD;
            a = (a * a) % MOD;
            e >>= 1;
        }
        return r;
    }

    static long[] countByMask(int maxLen) {
        long[] cnt = new long[FULL + 1];
        long[] dp = new long[FULL + 1];
        long[] ndp = new long[FULL + 1];

        for (int d = 1; d <= 9; d++) {
            dp[1 << d]++;
        }

        for (int len = 1; len <= maxLen; len++) {
            for (int m = 0; m <= FULL; m++) {
                cnt[m] += dp[m];
                if (cnt[m] >= MOD)
                    cnt[m] -= MOD;
            }
            if (len == maxLen)
                break;

            for (int m = 0; m <= FULL; m++)
                ndp[m] = 0;
            for (int m = 0; m <= FULL; m++) {
                long v = dp[m];
                if (v == 0)
                    continue;
                for (int d = 0; d <= 9; d++) {
                    int nm = m | (1 << d);
                    ndp[nm] += v;
                    if (ndp[nm] >= MOD)
                        ndp[nm] -= MOD;
                }
            }
            long[] tmp = dp;
            dp = ndp;
            ndp = tmp;
        }
        return cnt;
    }

    static String solveF(int K) {
        long inv2 = (MOD + 1) / 2;
        long[] cnt = countByMask(K);

        long[] sumSub = new long[FULL + 1];
        System.arraycopy(cnt, 0, sumSub, 0, FULL + 1);

        for (int b = 0; b < DIGS; b++) {
            int bit = 1 << b;
            for (int mask = 0; mask <= FULL; mask++) {
                if ((mask & bit) != 0) {
                    sumSub[mask] += sumSub[mask ^ bit];
                    if (sumSub[mask] >= MOD)
                        sumSub[mask] -= MOD;
                }
            }
        }

        long orderedDisjoint = 0;
        for (int a = 0; a <= FULL; a++) {
            long term = (cnt[a] * sumSub[FULL ^ a]) % MOD;
            orderedDisjoint = (orderedDisjoint + term) % MOD;
        }
        long disjointPairs = (orderedDisjoint * inv2) % MOD;

        long M = (modPow(10, K) - 1 + MOD) % MOD;
        long totalPairs = (((M * (M - 1 + MOD) % MOD) % MOD) * inv2) % MOD;

        long result = (totalPairs - disjointPairs + MOD) % MOD;
        return Long.toString(result);
    }

    public static String solve() {
        return solveF(18);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}