Problem 467: Superinteger

Project Euler Problem 467 Solution

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

Problem Summary Let $p_1,p_2,\dots$ be the primes and $c_1,c_2,\dots$ be the composite numbers. Define the digital-root sequences $$P_n=(\operatorname{dr}(p_1),\dots,\operatorname{dr}(p_n)),\qquad C_n=(\operatorname{dr}(c_1),\dots,\operatorname{dr}(c_n)),$$ where $$\operatorname{dr}(x)=1+((x-1)\bmod 9).$$ The task is to form the shortest decimal digit string that contains both $P_n$ and $C_n$ as subsequences. If several shortest strings exist, we take the lexicographically smallest one; because every digit lies in $\{1,\dots,9\}$, this is also the numerically smallest shortest superinteger. The required output is that enormous integer modulo $10^9+7$, and the implementations target $n=10000$. Mathematical Approach Write $$A=(a_1,\dots,a_n)=P_n,\qquad B=(b_1,\dots,b_n)=C_n.$$ We must compute the lexicographically least shortest common supersequence of these two digit sequences. Step 1: Build the Two Digit Sequences The first stage is purely arithmetic. We scan the integers from $2$ upward, separate primes from composites, and replace each accepted value by its digital root. Because $$\operatorname{dr}(x)\in\{1,2,\dots,9\},$$ both sequences consist only of ordinary decimal digits and contain no zeros. This matters later, because equal-length digit strings are ordered lexicographically exactly as their decimal integers are ordered....

Detailed mathematical approach

Problem Summary

Let $p_1,p_2,\dots$ be the primes and $c_1,c_2,\dots$ be the composite numbers. Define the digital-root sequences

$$P_n=(\operatorname{dr}(p_1),\dots,\operatorname{dr}(p_n)),\qquad C_n=(\operatorname{dr}(c_1),\dots,\operatorname{dr}(c_n)),$$

where

$$\operatorname{dr}(x)=1+((x-1)\bmod 9).$$

The task is to form the shortest decimal digit string that contains both $P_n$ and $C_n$ as subsequences. If several shortest strings exist, we take the lexicographically smallest one; because every digit lies in $\{1,\dots,9\}$, this is also the numerically smallest shortest superinteger. The required output is that enormous integer modulo $10^9+7$, and the implementations target $n=10000$.

Mathematical Approach

Write

$$A=(a_1,\dots,a_n)=P_n,\qquad B=(b_1,\dots,b_n)=C_n.$$

We must compute the lexicographically least shortest common supersequence of these two digit sequences.

Step 1: Build the Two Digit Sequences

The first stage is purely arithmetic. We scan the integers from $2$ upward, separate primes from composites, and replace each accepted value by its digital root. Because

$$\operatorname{dr}(x)\in\{1,2,\dots,9\},$$

both sequences consist only of ordinary decimal digits and contain no zeros. This matters later, because equal-length digit strings are ordered lexicographically exactly as their decimal integers are ordered.

Step 2: Reformulate the Problem as SCS

A shortest common supersequence (SCS) of $A$ and $B$ is a shortest string that contains both sequences in order, not necessarily contiguously. Every valid answer must preserve the internal order of $A$ and the internal order of $B$. If $a_i=b_j$, one output digit can serve both sequences simultaneously; if they differ, the next output digit must come from either $A$ or $B$.

So the problem is not about decimal arithmetic first. It is a sequence-merging problem whose final merged digit string is interpreted as an integer only at the end.

Step 3: Dynamic Programming for the Optimal Length

Let $L(i,j)$ be the length of the shortest common supersequence of the suffixes

$$a_i a_{i+1}\dots a_n\qquad\text{and}\qquad b_j b_{j+1}\dots b_n,$$

with the convention that $L(n+1,n+1)=0$. The boundary cases are immediate:

$$L(n+1,j)=n-j+1,\qquad L(i,n+1)=n-i+1.$$

For interior states we have

$$L(i,j)= \begin{cases} 1+L(i+1,j+1), & a_i=b_j,\\ 1+\min\{L(i+1,j),L(i,j+1)\}, & a_i\ne b_j. \end{cases}$$

The recurrence is standard: matching digits are consumed together, while different digits force a choice of which sequence contributes the next output digit.

Step 4: Recover the Lexicographically Smallest Optimal Answer

The length table tells us which moves keep the answer shortest. Reconstruction then follows these rules:

1. If one sequence is exhausted, append the rest of the other sequence.

2. If $a_i=b_j$, output that digit once and advance in both sequences.

3. If $a_i\ne b_j$, compare $L(i+1,j)$ and $L(i,j+1)$.

If one option gives a smaller remaining length, that branch is forced. If both options preserve the optimal length, then both produce shortest supersequences, so the first digit decides lexicographic order immediately. Therefore we output $\min(a_i,b_j)$. By induction on the remaining suffix length, this greedy tie-break yields the lexicographically smallest shortest supersequence.

Step 5: Evaluate the Huge Integer Modulo $10^9+7$

If the final digit string is $d_1d_2\dots d_t$, its value modulo $10^9+7$ can be accumulated online by

$$V_0=0,\qquad V_{k+1}\equiv 10V_k+d_{k+1}\pmod{10^9+7}.$$

This avoids constructing a huge big integer. The algorithm only needs the next chosen digit, so reconstruction and modular evaluation can happen in the same pass.

Worked Example: $n=10$

The first ten prime digital roots are

$$P_{10}=(2,3,5,7,2,4,8,1,5,2),$$

and the first ten composite digital roots are

$$C_{10}=(4,6,8,9,1,3,5,6,7,9).$$

The common subsequence $(4,8,1,5)$ already shows that the answer can be shorter than length $20$, and the dynamic-programming table proves that the true optimum length is $16$. After lexicographic tie-breaking, the shortest supersequence becomes

$$2357246891352679.$$

This is the checkpoint used by the implementations for $f(10)$. They also verify that $f(100)\equiv 771661825 \pmod{10^9+7}$.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. First they sieve enough integers to collect the first $n$ primes and the first $n$ composites, converting each accepted number to its digital root immediately. Next they fill an $(n+1)\times(n+1)$ dynamic-programming table from bottom right to top left so that every state already knows the optimal remaining length of its suffixes.

After that, the implementation reconstructs one digit at a time. Equal digits are merged, unequal digits consult the two neighboring table entries, and ties are broken by the smaller next digit. At the same moment, the answer is updated modulo $10^9+7$. For small checkpoint cases the full digit string can also be materialized, but that is not needed for the final $n=10000$ computation.

Complexity Analysis

Let $M$ be the largest integer that must be sieved in order to collect $n$ primes and $n$ composites. Building primality information up to $M$ costs $O(M\log\log M)$ time and $O(M)$ memory. Since the prime side is the sparse one, $M$ is on the order of the $n$-th prime, so $M=\Theta(n\log n)$.

The SCS stage dominates: filling the length table takes $O(n^2)$ time and $O(n^2)$ memory. Reconstruction adds only $O(t)$ time, where $t\le 2n$. For the target $n=10000$, the quadratic dynamic programming is the main cost.

Footnotes and References

Problem page: https://projecteuler.net/problem=467
Digital root: Wikipedia — Digital root
Shortest common supersequence: Wikipedia — Shortest common supersequence problem
Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes
Dynamic programming: Wikipedia — Dynamic programming

Problem 467 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#include <functional>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;

constexpr u32 kMod = 1'000'000'007U;

struct Options {
    int n = 10'000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        out = std::stoi(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    if (options.n <= 0) {
        std::cerr << "--n must be positive.\n";
        return false;
    }
    return true;
}

int estimate_limit_for_prime_count(const int n) {
    if (n < 6) {
        return 32;
    }
    const double dn = static_cast<double>(n);
    const double estimate = dn * (std::log(dn) + std::log(std::log(dn)));
    return static_cast<int>(estimate) + 64;
}

std::vector<bool> sieve(const int limit) {
    std::vector<bool> is_prime(static_cast<std::size_t>(limit + 1), true);
    is_prime[0] = false;
    is_prime[1] = false;
    for (int p = 2; static_cast<int64_t>(p) * p <= limit; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int x = p * p; x <= limit; x += p) {
            is_prime[static_cast<std::size_t>(x)] = false;
        }
    }
    return is_prime;
}

std::uint8_t digital_root(const int x) {
    return static_cast<std::uint8_t>(1 + (x - 1) % 9);
}

void build_sequences(const int n, std::vector<std::uint8_t>& prime_digits,
                     std::vector<std::uint8_t>& composite_digits) {
    int limit = std::max(64, estimate_limit_for_prime_count(n));
    while (true) {
        const std::vector<bool> is_prime = sieve(limit);
        prime_digits.clear();
        composite_digits.clear();
        prime_digits.reserve(static_cast<std::size_t>(n));
        composite_digits.reserve(static_cast<std::size_t>(n));

        for (int x = 2; x <= limit; ++x) {
            if (is_prime[static_cast<std::size_t>(x)]) {
                if (static_cast<int>(prime_digits.size()) < n) {
                    prime_digits.push_back(digital_root(x));
                }
            } else {
                if (static_cast<int>(composite_digits.size()) < n) {
                    composite_digits.push_back(digital_root(x));
                }
            }
            if (static_cast<int>(prime_digits.size()) == n &&
                static_cast<int>(composite_digits.size()) == n) {
                return;
            }
        }
        limit *= 2;
    }
}

inline std::size_t idx(const int i, const int j, const int width) {
    return static_cast<std::size_t>(i) * static_cast<std::size_t>(width) +
           static_cast<std::size_t>(j);
}

u32 compute_scs_mod(const std::vector<std::uint8_t>& a, const std::vector<std::uint8_t>& b,
                    std::string* out_string = nullptr) {
    const int n = static_cast<int>(a.size());
    const int width = n + 1;
    std::vector<std::uint16_t> dp(static_cast<std::size_t>(width) * static_cast<std::size_t>(width),
                                  0U);

    for (int i = n; i >= 0; --i) {
        for (int j = n; j >= 0; --j) {
            const std::size_t at = idx(i, j, width);
            if (i == n) {
                dp[at] = static_cast<std::uint16_t>(n - j);
                continue;
            }
            if (j == n) {
                dp[at] = static_cast<std::uint16_t>(n - i);
                continue;
            }
            if (a[static_cast<std::size_t>(i)] == b[static_cast<std::size_t>(j)]) {
                dp[at] = static_cast<std::uint16_t>(1 + dp[idx(i + 1, j + 1, width)]);
            } else {
                const std::uint16_t da = dp[idx(i + 1, j, width)];
                const std::uint16_t db = dp[idx(i, j + 1, width)];
                dp[at] = static_cast<std::uint16_t>(1 + (da < db ? da : db));
            }
        }
    }

    int i = 0;
    int j = 0;
    u32 result_mod = 0U;
    if (out_string != nullptr) {
        out_string->clear();
        out_string->reserve(static_cast<std::size_t>(dp[0]));
    }

    while (i < n || j < n) {
        std::uint8_t digit = 0U;
        if (i == n) {
            digit = b[static_cast<std::size_t>(j++)];
        } else if (j == n) {
            digit = a[static_cast<std::size_t>(i++)];
        } else if (a[static_cast<std::size_t>(i)] == b[static_cast<std::size_t>(j)]) {
            digit = a[static_cast<std::size_t>(i)];
            ++i;
            ++j;
        } else {
            const std::uint16_t da = dp[idx(i + 1, j, width)];
            const std::uint16_t db = dp[idx(i, j + 1, width)];
            if (da < db || (da == db && a[static_cast<std::size_t>(i)] < b[static_cast<std::size_t>(j)])) {
                digit = a[static_cast<std::size_t>(i++)];
            } else {
                digit = b[static_cast<std::size_t>(j++)];
            }
        }

        result_mod = static_cast<u32>((static_cast<u64>(result_mod) * 10ULL + digit) % kMod);
        if (out_string != nullptr) {
            out_string->push_back(static_cast<char>('0' + digit));
        }
    }

    return result_mod;
}

u32 solve_mod(const int n, std::string* out_string = nullptr) {
    std::vector<std::uint8_t> prime_digits;
    std::vector<std::uint8_t> composite_digits;
    build_sequences(n, prime_digits, composite_digits);
    return compute_scs_mod(prime_digits, composite_digits, out_string);
}

bool run_checkpoints() {
    std::string f10;
    solve_mod(10, &f10);
    if (f10 != "2357246891352679") {
        std::cerr << "Checkpoint failed: f(10)\n";
        return false;
    }
    if (solve_mod(100, nullptr) != 771'661'825U) {
        std::cerr << "Checkpoint failed: f(100) mod 1e9+7\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 1;
    }

    std::cout << solve_mod(options.n, nullptr) << '\n';
    return 0;
}

Python

import math

def solve():
    MOD = 1000000007
    n = 10000

    def digital_root(x): return 1 + (x - 1) % 9

    limit = max(64, int(n * (math.log(n) + math.log(math.log(n))))) + 64
    sieve = bytearray(b'\x01') * (limit + 1); sieve[0] = sieve[1] = 0
    for p in range(2, int(limit**0.5)+1):
        if sieve[p]:
            for j in range(p*p, limit+1, p): sieve[j] = 0

    pd, cd = [], []
    for x in range(2, limit + 1):
        if sieve[x]:
            if len(pd) < n: pd.append(digital_root(x))
        else:
            if len(cd) < n: cd.append(digital_root(x))
        if len(pd) == n and len(cd) == n: break

    # SCS DP
    w = n + 1
    dp = [[0]*w for _ in range(w)]
    for i in range(n, -1, -1):
        for j in range(n, -1, -1):
            if i == n: dp[i][j] = n - j
            elif j == n: dp[i][j] = n - i
            elif pd[i] == cd[j]: dp[i][j] = 1 + dp[i+1][j+1]
            else: dp[i][j] = 1 + min(dp[i+1][j], dp[i][j+1])

    i, j, result = 0, 0, 0
    while i < n or j < n:
        if i == n: d = cd[j]; j += 1
        elif j == n: d = pd[i]; i += 1
        elif pd[i] == cd[j]: d = pd[i]; i += 1; j += 1
        else:
            da, db = dp[i+1][j], dp[i][j+1]
            if da < db or (da == db and pd[i] < cd[j]):
                d = pd[i]; i += 1
            else: d = cd[j]; j += 1
        result = (result * 10 + d) % MOD

    return str(result)

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

Java

public class Euler467 {
    public static String solve() {
        int kMod = 1000000007;
        int n = 10000;

        int limit = 150000;
        boolean[] isPrime = new boolean[limit + 1];
        for (int i = 2; i <= limit; i++)
            isPrime[i] = true;

        for (int p = 2; p * p <= limit; p++) {
            if (isPrime[p]) {
                for (int i = p * p; i <= limit; i += p) {
                    isPrime[i] = false;
                }
            }
        }

        byte[] primeDigits = new byte[n];
        byte[] compositeDigits = new byte[n];
        int pCount = 0;
        int cCount = 0;

        for (int x = 2; x <= limit; x++) {
            if (isPrime[x]) {
                if (pCount < n) {
                    primeDigits[pCount++] = (byte) (1 + (x - 1) % 9);
                }
            } else {
                if (cCount < n) {
                    compositeDigits[cCount++] = (byte) (1 + (x - 1) % 9);
                }
            }
            if (pCount == n && cCount == n)
                break;
        }

        int width = n + 1;
        short[] dp = new short[width * width];

        for (int i = n; i >= 0; i--) {
            int rowOffset = i * width;
            int nextRowOffset = (i + 1) * width;
            byte pi = i < n ? primeDigits[i] : 0;

            for (int j = n; j >= 0; j--) {
                int at = rowOffset + j;
                if (i == n) {
                    dp[at] = (short) (n - j);
                    continue;
                }
                if (j == n) {
                    dp[at] = (short) (n - i);
                    continue;
                }

                if (pi == compositeDigits[j]) {
                    dp[at] = (short) (1 + dp[nextRowOffset + j + 1]);
                } else {
                    short da = dp[nextRowOffset + j];
                    short db = dp[rowOffset + j + 1];
                    dp[at] = (short) (1 + (da < db ? da : db));
                }
            }
        }

        int i = 0;
        int j = 0;
        long resultMod = 0;

        while (i < n || j < n) {
            byte digit;
            if (i == n) {
                digit = compositeDigits[j++];
            } else if (j == n) {
                digit = primeDigits[i++];
            } else if (primeDigits[i] == compositeDigits[j]) {
                digit = primeDigits[i];
                i++;
                j++;
            } else {
                short da = dp[(i + 1) * width + j];
                short db = dp[i * width + j + 1];
                if (da < db || (da == db && primeDigits[i] < compositeDigits[j])) {
                    digit = primeDigits[i++];
                } else {
                    digit = compositeDigits[j++];
                }
            }
            resultMod = (resultMod * 10 + digit) % kMod;
        }

        return Long.toString(resultMod);
    }

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