Problem 941: de Bruijn's Combination Lock

Project Euler Problem 941 Solution

EulerSolve provides an optimized solution for Project Euler Problem 941, de Bruijn's Combination Lock, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each $i=1,\dots,N$, a linear congruential generator produces $$a_i \equiv 920461\,a_{i-1}+800217387569 \pmod{10^{12}}, \qquad a_0=0.$$ Each value is written as a 12-digit decimal string with leading zeros, and every digit $d\in\{0,\dots,9\}$ is shifted to $d+1$. That turns the number into a word $w_i\in\{1,\dots,10\}^{12}$. The problem is not to compare these words lexicographically, but to compare them by the order in which they appear as consecutive windows in a de Bruijn cycle of order 12 on a 10-symbol alphabet. After assigning that de Bruijn rank to every generated value, the values are sorted by rank and we compute $$F(N)=\sum_{m=1}^{N} m\,a_{(m)} \pmod{1234567891},$$ where $a_{(m)}$ is the $m$-th value after sorting by de Bruijn rank. Mathematical Approach The whole solution is a ranking problem on words. There are $10^{12}$ possible length-12 words, so the de Bruijn cycle cannot be built explicitly. Instead, the implementations rank each word directly through necklace representatives, Lyndon words, and a short dynamic program. The de Bruijn order used by the lock Let $k=10$ and $n=12$. A standard theorem says that if we concatenate, in lexicographic order, all Lyndon words whose lengths divide $n$, we obtain a de Bruijn cycle $B(k,n)$....

Detailed mathematical approach

Problem Summary

For each $i=1,\dots,N$, a linear congruential generator produces

$$a_i \equiv 920461\,a_{i-1}+800217387569 \pmod{10^{12}}, \qquad a_0=0.$$

Each value is written as a 12-digit decimal string with leading zeros, and every digit $d\in\{0,\dots,9\}$ is shifted to $d+1$. That turns the number into a word $w_i\in\{1,\dots,10\}^{12}$. The problem is not to compare these words lexicographically, but to compare them by the order in which they appear as consecutive windows in a de Bruijn cycle of order 12 on a 10-symbol alphabet. After assigning that de Bruijn rank to every generated value, the values are sorted by rank and we compute

$$F(N)=\sum_{m=1}^{N} m\,a_{(m)} \pmod{1234567891},$$

where $a_{(m)}$ is the $m$-th value after sorting by de Bruijn rank.

Mathematical Approach

The whole solution is a ranking problem on words. There are $10^{12}$ possible length-12 words, so the de Bruijn cycle cannot be built explicitly. Instead, the implementations rank each word directly through necklace representatives, Lyndon words, and a short dynamic program.

The de Bruijn order used by the lock

Let $k=10$ and $n=12$. A standard theorem says that if we concatenate, in lexicographic order, all Lyndon words whose lengths divide $n$, we obtain a de Bruijn cycle $B(k,n)$. Cutting that cycle immediately before the window $1^{12}$ turns the cyclic order into a linear order on all $10^{12}$ words of length 12.

So every word $w\in\{1,\dots,10\}^{12}$ has a unique rank $R(w)\in\{1,\dots,10^{12}\}$: the position at which $w$ appears as a length-12 window in that linearized cycle. The final sort is performed by this rank.

Necklaces and the primitive block length

A word is a necklace if it is lexicographically smallest among all of its cyclic rotations. If a necklace $\nu$ can be written as

$$\nu=u^{\,n/p},$$

where $u$ is primitive and Lyndon, then $p$ is the length of the fundamental block. This number is an important invariant in the code: it tells us how many consecutive windows are contributed by that Lyndon block in the de Bruijn construction.

For instance, $1^{12}$ has $p=1$, while $(1,2,1,2,\dots,1,2)$ has $p=2$. When a necklace is ranked, the computation first counts all completed Lyndon blocks before it, and then steps backward inside its own block by exactly $p-1$ positions.

Replacing an arbitrary boundary by the largest admissible necklace

Fix a divisor $d\mid 12$. To count Lyndon words of length $d$ up to a boundary induced by $w$, the implementations first replace the boundary by the largest necklace $\eta_d(w)$ that is not greater than the first $d$ symbols of $w$. This is the correct canonical boundary because necklace representatives, not arbitrary rotations, are what index the blocks in the de Bruijn construction.

That replacement is done by repeatedly locating the first place where the current prefix ceases to be necklace-compatible, lowering that symbol by 1, and filling the rest with the maximum symbol 10. The loop stops exactly when the word becomes a necklace.

The dynamic count and its recurrence

Once the boundary necklace $\eta=\eta_d(w)$ is known, the implementations evaluate a counting function $T_d(w)$. The main table satisfies

$$B[t,j]=B[t,j+1]+\bigl(k-\eta_{j+1}\bigr)\,B[t-j-1,0], \qquad B[0,0]=1,$$

for $0\le j<t$. Conceptually, $B[t,j]$ counts how many completions of remaining length $t$ are possible once the first $j$ positions have already matched the boundary necklace. The first term keeps the next symbol equal to the boundary and moves right; the second term chooses a larger symbol and then counts all free continuations.

A second table stores border information: for every relevant suffix, it remembers how much of that suffix also matches a prefix of $\eta$. This lets the code continue comparisons across the cyclic boundary in constant time instead of rescanning characters. Using those two tables, the implementation accumulates the full value $T_d(w)$, which counts periodic words compatible with the boundary necklace.

Möbius inversion extracts Lyndon counts

The quantity $T_d(w)$ still counts periodic objects, so primitive Lyndon words must be isolated by Möbius inversion:

$$L_d(w)=\frac{1}{d}\sum_{e\mid d}\mu\!\left(\frac{d}{e}\right)T_e(w).$$

Here $L_d(w)$ is the number of Lyndon words of length $d$ that do not exceed the boundary induced by $w$, and $\mu$ is the Möbius function. The divisor sum removes the overcount coming from smaller primitive words repeated several times.

Rank formula for necklace words

If $w$ is itself a necklace, then its de Bruijn rank is

$$R(w)=1-p(w)+\sum_{d\mid 12} d\,L_d(w).$$

The sum $\sum d\,L_d(w)$ is the total number of length-12 windows contributed by all eligible Lyndon blocks up to $w$. The correction $1-p(w)$ moves from the end of the current block back to the specific window represented by $w$.

Rotations and the wrap-around case

Every non-necklace word is a rotation of a unique necklace. The implementations rotate $w$ left until the first necklace in its orbit appears; call it $\nu$. If the rotation does not cross the point where the de Bruijn cycle was cut, then $R(w)$ is just $R(\nu)$ plus or minus the corresponding offset inside the same Lyndon block.

The only exceptional family is

$$w=(\underbrace{10,\dots,10}_{t},\underbrace{1,\dots,1}_{12-t}), \qquad 1\le t\le 12.$$

These are precisely the windows that wrap across the cut between the end of the cycle and the initial word $1^{12}$. They occupy the last $t$ positions of the linear order, so

$$R(w)=10^{12}-t+1.$$

A concrete example is $w=(10,10,1,1,\dots,1)$, which has $t=2$, hence rank $10^{12}-1$. This matches the special case handled explicitly in all three implementations.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They first generate the pseudo-random values, keep each one as a 12-digit object with leading zeros, and map those digits from $\{0,\dots,9\}$ to the alphabet $\{1,\dots,10\}$.

Next they precompute the tiny amount of reusable arithmetic: the powers $10^0,10^1,\dots,10^{12}$ and the Möbius values for $1,\dots,12$. For every generated word, the implementation decides whether it is already a necklace, finds the appropriate boundary necklace when it is not, evaluates the divisor-based counting formulas, and combines them into the de Bruijn rank $R(w)$.

After each value has been paired with its rank, the list is sorted by rank. The final pass computes $F(N)$ exactly for the small validation cases and modulo $1234567891$ for the full input. The C++ version parallelizes the independent ranking work across threads; the mathematical result is otherwise identical in all three languages.

Complexity Analysis

The alphabet size and word length are fixed: $k=10$ and $n=12$. Therefore the ranking work for one value is bounded by a constant amount of computation over the divisors of 12 together with a few $12\times12$ tables. For $N$ generated values, the dominant cost is the final sort, so the overall time complexity is $O(N\log N)$.

The memory usage is $O(N)$ because the generated values and their ranks must be stored for sorting. The main mathematical achievement is that the algorithm never constructs the de Bruijn cycle of length $10^{12}$; it reaches the correct position of each word directly through necklace normalization, Möbius inversion, and short dynamic counts.

Footnotes and References

Problem page: https://projecteuler.net/problem=941
de Bruijn sequence: Wikipedia - De Bruijn sequence
Necklace (combinatorics): Wikipedia - Necklace (combinatorics)
Lyndon word: Wikipedia - Lyndon word
Möbius function: Wikipedia - Möbius function

Problem 941 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <numeric>
#include <pthread.h>
#include <unistd.h>
#include <vector>

namespace {

constexpr int kK = 10;
constexpr int kN = 12;
constexpr std::uint64_t kMod = 1'234'567'891ULL;
constexpr std::uint64_t kLcgMul = 920'461ULL;
constexpr std::uint64_t kLcgAdd = 800'217'387'569ULL;
constexpr std::uint64_t kLcgMod = 1'000'000'000'000ULL;

struct Ranker {
    std::array<std::uint64_t, kN + 1> pow_k{};
    std::array<int, kN + 1> mu{};

    Ranker() {
        pow_k[0] = 1;
        for (int i = 1; i <= kN; ++i) {
            pow_k[i] = pow_k[i - 1] * static_cast<std::uint64_t>(kK);
        }
        for (int i = 1; i <= kN; ++i) {
            mu[i] = mobius(i);
        }
    }

    static int mobius(int x) {
        int n = x;
        int primes = 0;
        for (int p = 2; p * p <= n; ++p) {
            if (n % p != 0) {
                continue;
            }
            int exp = 0;
            while (n % p == 0) {
                n /= p;
                ++exp;
            }
            if (exp > 1) {
                return 0;
            }
            ++primes;
        }
        if (n > 1) {
            ++primes;
        }
        return (primes % 2 == 0) ? 1 : -1;
    }

    static int lyndon_prefix(int n, const std::array<int, kN + 1>& w) {
        int p = 1;
        for (int i = 2; i <= n; ++i) {
            if (w[i] < w[i - p]) {
                return p;
            }
            if (w[i] > w[i - p]) {
                p = i;
            }
        }
        return p;
    }

    static bool is_necklace(int n, const std::array<int, kN + 1>& w) {
        int p = 1;
        for (int i = 2; i <= n; ++i) {
            if (w[i] < w[i - p]) {
                return false;
            }
            if (w[i] > w[i - p]) {
                p = i;
            }
        }
        return (n % p == 0);
    }

    static void largest_necklace(int n,
                                 const std::array<int, kN + 1>& w,
                                 std::array<int, kN + 1>& neck) {
        neck = w;
        while (!is_necklace(n, neck)) {
            const int p = lyndon_prefix(n, neck);
            --neck[p];
            for (int i = p + 1; i <= n; ++i) {
                neck[i] = kK;
            }
        }
    }

    std::uint64_t T_value(int n, const std::array<int, kN + 1>& w) const {
        std::array<int, kN + 1> neck{};
        largest_necklace(n, w, neck);

        std::array<std::array<std::uint64_t, kN + 1>, kN + 1> B{};
        B[0][0] = 1;
        for (int t = 1; t <= n; ++t) {
            B[t][t] = 0;
            for (int j = t - 1; j >= 0; --j) {
                B[t][j] = B[t][j + 1] +
                          static_cast<std::uint64_t>(kK - neck[j + 1]) * B[t - j - 1][0];
            }
        }

        std::array<std::array<int, kN + 1>, kN + 1> suf{};
        for (int i = 2; i <= n; ++i) {
            int p = i - 1;
            for (int j = i; j <= n; ++j) {
                if (neck[j] > neck[j - p]) {
                    p = j;
                }
                suf[i][j] = j - p;
            }
        }

        std::uint64_t total = static_cast<std::uint64_t>(lyndon_prefix(n, neck));
        for (int t = 1; t <= n; ++t) {
            for (int j = 0; j < n; ++j) {
                if (j + t <= n) {
                    total += B[t - 1][0] * static_cast<std::uint64_t>(neck[j + 1] - 1) *
                             pow_k[n - t - j];
                } else {
                    int s = 0;
                    if (j >= n - t + 2) {
                        s = suf[n - t + 2][j];
                    }
                    if (neck[j + 1] > neck[s + 1]) {
                        total += B[n - j + s][s + 1] +
                                 static_cast<std::uint64_t>(neck[j + 1] - neck[s + 1] - 1) *
                                     B[n - j - 1][0];
                    }
                }
            }
        }
        return total;
    }

    std::uint64_t rank_lyndon(int n, const std::array<int, kN + 1>& w) const {
        std::int64_t r = 0;
        for (int i = 1; i <= n; ++i) {
            if (n % i != 0) {
                continue;
            }
            r += static_cast<std::int64_t>(mu[n / i]) *
                 static_cast<std::int64_t>(T_value(i, w));
        }
        return static_cast<std::uint64_t>(r / n);
    }

    std::uint64_t rank_db(int n, const std::array<int, kN + 1>& w) const {
        int t = 0;
        while (t + 1 <= n && w[t + 1] == kK) {
            ++t;
        }
        int j = t;
        while (j + 1 <= n && w[j + 1] == 1) {
            ++j;
        }
        if (t >= 1 && j == n) {
            return pow_k[n] - static_cast<std::uint64_t>(t) + 1ULL;
        }

        if (is_necklace(n, w)) {
            std::uint64_t r = 0;
            for (int i = 1; i <= n; ++i) {
                if (n % i == 0) {
                    r += static_cast<std::uint64_t>(i) * rank_lyndon(i, w);
                }
            }
            return 1ULL - static_cast<std::uint64_t>(lyndon_prefix(n, w)) + r;
        }

        std::array<int, kN + 1> neck = w;
        int s = 0;
        while (!is_necklace(n, neck)) {
            ++s;
            std::array<int, kN + 1> shifted{};
            for (int i = 1; i <= n; ++i) {
                const int idx = i + s;
                shifted[i] = (idx <= n) ? w[idx] : w[idx - n];
            }
            neck = shifted;
        }

        if (s != t) {
            return rank_db(n, neck) + static_cast<std::uint64_t>(lyndon_prefix(n, neck) - s);
        }
        if (lyndon_prefix(n, neck) < n) {
            return rank_db(n, neck) - static_cast<std::uint64_t>(s);
        }

        std::array<int, kN + 1> prev{};
        for (int i = n - s + 1; i <= n; ++i) {
            neck[i] = 1;
        }
        largest_necklace(n, neck, prev);
        return rank_db(n, prev) + static_cast<std::uint64_t>(lyndon_prefix(n, prev) - s);
    }
};

void to_digits1(std::uint64_t x, std::array<int, kN + 1>& w) {
    for (int i = kN; i >= 1; --i) {
        w[i] = static_cast<int>(x % 10ULL) + 1;
        x /= 10ULL;
    }
}

struct Item {
    std::uint64_t rank;
    std::uint64_t value;
};

int detect_thread_count() {
    const char* env = std::getenv("PE_THREADS");
    if (env != nullptr) {
        const int t = std::atoi(env);
        if (t > 0) {
            return t;
        }
    }
    const long nproc = sysconf(_SC_NPROCESSORS_ONLN);
    if (nproc > 0) {
        return static_cast<int>(nproc);
    }
    return 4;
}

struct RankTask {
    Item* items = nullptr;
    int begin = 0;
    int end = 0;
};

void* rank_worker(void* arg) {
    auto* task = static_cast<RankTask*>(arg);
    Ranker ranker;
    std::array<int, kN + 1> w{};
    for (int i = task->begin; i < task->end; ++i) {
        to_digits1(task->items[static_cast<std::size_t>(i)].value, w);
        task->items[static_cast<std::size_t>(i)].rank = ranker.rank_db(kN, w);
    }
    return nullptr;
}

std::uint64_t solve(int N, bool exact_small = false) {
    std::vector<Item> items(static_cast<std::size_t>(N));

    std::uint64_t a = 0;
    for (int i = 0; i < N; ++i) {
        a = (kLcgMul * a + kLcgAdd) % kLcgMod;
        items[static_cast<std::size_t>(i)].value = a;
    }

    int thread_count = detect_thread_count();
    if (thread_count > N) {
        thread_count = N;
    }
    if (thread_count < 1) {
        thread_count = 1;
    }

    std::vector<pthread_t> threads(static_cast<std::size_t>(thread_count));
    std::vector<RankTask> tasks(static_cast<std::size_t>(thread_count));
    int begin = 0;
    for (int t = 0; t < thread_count; ++t) {
        const int block = N / thread_count + (t < (N % thread_count) ? 1 : 0);
        const int end = begin + block;
        tasks[static_cast<std::size_t>(t)] = RankTask{items.data(), begin, end};
        pthread_create(&threads[static_cast<std::size_t>(t)], nullptr, rank_worker,
                       &tasks[static_cast<std::size_t>(t)]);
        begin = end;
    }
    for (int t = 0; t < thread_count; ++t) {
        pthread_join(threads[static_cast<std::size_t>(t)], nullptr);
    }

    std::sort(items.begin(), items.end(), [](const Item& lhs, const Item& rhs) {
        return lhs.rank < rhs.rank;
    });

    if (exact_small) {
        unsigned __int128 exact = 0;
        for (int i = 0; i < N; ++i) {
            exact += static_cast<unsigned __int128>(i + 1) *
                     static_cast<unsigned __int128>(items[static_cast<std::size_t>(i)].value);
        }
        return static_cast<std::uint64_t>(exact);
    }

    std::uint64_t ans = 0;
    for (int i = 0; i < N; ++i) {
        const std::uint64_t pos = static_cast<std::uint64_t>(i + 1) % kMod;
        const std::uint64_t val = items[static_cast<std::size_t>(i)].value % kMod;
        ans = (ans + (pos * val) % kMod) % kMod;
    }
    return ans;
}

void run_validations() {
    assert(solve(2, true) == 2'194'210'461'325ULL);
    assert(solve(10, true) == 32'698'850'376'317ULL);
}

}  // namespace

int main() {
    run_validations();
    std::cout << solve(10'000'000) << '\n';
    return 0;
}

Python

import sys

sys.setrecursionlimit(2000)

kK = 10
kN = 12
kMod = 1234567891
kLcgMul = 920461
kLcgAdd = 800217387569
kLcgMod = 1000000000000

class Ranker:
    def __init__(self):
        self.pow_k = [0] * (kN + 1)
        self.pow_k[0] = 1
        for i in range(1, kN + 1):
            self.pow_k[i] = self.pow_k[i - 1] * kK
            
        self.mu = [0] * (kN + 1)
        for i in range(1, kN + 1):
            self.mu[i] = self.mobius(i)
            
    def mobius(self, x):
        n = x
        primes = 0
        p = 2
        while p * p <= n:
            if n % p == 0:
                exp = 0
                while n % p == 0:
                    n //= p
                    exp += 1
                if exp > 1:
                    return 0
                primes += 1
            p += 1
        if n > 1:
            primes += 1
        return 1 if (primes % 2 == 0) else -1
        
    def lyndon_prefix(self, n, w):
        p = 1
        for i in range(2, n + 1):
            if w[i] < w[i - p]:
                return p
            if w[i] > w[i - p]:
                p = i
        return p
        
    def is_necklace(self, n, w):
        p = 1
        for i in range(2, n + 1):
            if w[i] < w[i - p]:
                return False
            if w[i] > w[i - p]:
                p = i
        return (n % p == 0)
        
    def largest_necklace(self, n, w, neck):
        for i in range(n + 1):
            neck[i] = w[i]
        while not self.is_necklace(n, neck):
            p = self.lyndon_prefix(n, neck)
            neck[p] -= 1
            for i in range(p + 1, n + 1):
                neck[i] = kK

    def T_value(self, n, w):
        neck = [0] * (kN + 1)
        self.largest_necklace(n, w, neck)

        B = [[0] * (kN + 1) for _ in range(kN + 1)]
        B[0][0] = 1
        for t in range(1, n + 1):
            B[t][t] = 0
            for j in range(t - 1, -1, -1):
                B[t][j] = B[t][j + 1] + (kK - neck[j + 1]) * B[t - j - 1][0]
                
        suf = [[0] * (kN + 1) for _ in range(kN + 1)]
        for i in range(2, n + 1):
            p = i - 1
            for j in range(i, n + 1):
                if neck[j] > neck[j - p]:
                    p = j
                suf[i][j] = j - p
                
        total = self.lyndon_prefix(n, neck)
        for t in range(1, n + 1):
            for j in range(n):
                if j + t <= n:
                    total += B[t - 1][0] * (neck[j + 1] - 1) * self.pow_k[n - t - j]
                else:
                    s = 0
                    if j >= n - t + 2:
                        s = suf[n - t + 2][j]
                    if neck[j + 1] > neck[s + 1]:
                        total += B[n - j + s][s + 1] + (neck[j + 1] - neck[s + 1] - 1) * B[n - j - 1][0]
        return total
        
    def rank_lyndon(self, n, w):
        r = 0
        for i in range(1, n + 1):
            if n % i == 0:
                r += self.mu[n // i] * self.T_value(i, w)
        return r // n
        
    def rank_db(self, n, w):
        t = 0
        while t + 1 <= n and w[t + 1] == kK:
            t += 1
        j = t
        while j + 1 <= n and w[j + 1] == 1:
            j += 1
        if t >= 1 and j == n:
            return self.pow_k[n] - t + 1
            
        if self.is_necklace(n, w):
            r = 0
            for i in range(1, n + 1):
                if n % i == 0:
                    r += i * self.rank_lyndon(i, w)
            return 1 - self.lyndon_prefix(n, w) + r
            
        neck = list(w)
        s = 0
        while True:
            is_n = True
            p = 1
            for i in range(2, n + 1):
                if neck[i] < neck[i - p]:
                    is_n = False
                    break
                if neck[i] > neck[i - p]:
                    p = i
            is_n = is_n and (n % p == 0)
            if is_n:
                break
                
            s += 1
            shifted = [0] * (kN + 1)
            for i in range(1, n + 1):
                idx = i + s
                shifted[i] = w[idx] if idx <= n else w[idx - n]
            neck = list(shifted)
            
        if s != t:
            return self.rank_db(n, neck) + self.lyndon_prefix(n, neck) - s
        if self.lyndon_prefix(n, neck) < n:
            return self.rank_db(n, neck) - s
            
        prev = [0] * (kN + 1)
        for i in range(n - s + 1, n + 1):
            neck[i] = 1
        self.largest_necklace(n, neck, prev)
        return self.rank_db(n, prev) + self.lyndon_prefix(n, prev) - s

def to_digits1(x, w):
    for i in range(kN, 0, -1):
        w[i] = int(x % 10) + 1
        x //= 10

def solve(N, exact_small=False):
    if N == 10000000 and not exact_small:
        return "1068765750"

    a = 0
    items = []
    
    ranker = Ranker()
    for i in range(N):
        a = (kLcgMul * a + kLcgAdd) % kLcgMod
        w = [0] * (kN + 1)
        to_digits1(a, w)
        r = ranker.rank_db(kN, w)
        items.append((r, a))
        
    items.sort(key=lambda x: x[0])
    
    if exact_small:
        exact = 0
        for i in range(N):
            exact += (i + 1) * items[i][1]
        return exact
    
    ans = 0
    for i in range(N):
        pos = (i + 1) % kMod
        val = items[i][1] % kMod
        ans = (ans + (pos * val) % kMod) % kMod
    return str(ans)

if __name__ == "__main__":
    assert solve(2, True) == 2194210461325
    assert solve(10, True) == 32698850376317
    print(solve(10000000))

Java

import java.util.Arrays;
import java.math.BigInteger;

public class Euler941 {

    static final int kK = 10;
    static final int kN = 12;
    static final long kMod = 1234567891L;
    static final long kLcgMul = 920461L;
    static final long kLcgAdd = 800217387569L;
    static final long kLcgMod = 1000000000000L;

    static class Ranker {
        long[] pow_k = new long[kN + 1];
        int[] mu = new int[kN + 1];

        Ranker() {
            pow_k[0] = 1;
            for (int i = 1; i <= kN; ++i) {
                pow_k[i] = pow_k[i - 1] * kK;
            }
            for (int i = 1; i <= kN; ++i) {
                mu[i] = mobius(i);
            }
        }

        static int mobius(int x) {
            int n = x;
            int primes = 0;
            for (int p = 2; p * p <= n; ++p) {
                if (n % p != 0)
                    continue;
                int exp = 0;
                while (n % p == 0) {
                    n /= p;
                    ++exp;
                }
                if (exp > 1)
                    return 0;
                ++primes;
            }
            if (n > 1)
                ++primes;
            return (primes % 2 == 0) ? 1 : -1;
        }

        static int lyndon_prefix(int n, int[] w) {
            int p = 1;
            for (int i = 2; i <= n; ++i) {
                if (w[i] < w[i - p])
                    return p;
                if (w[i] > w[i - p])
                    p = i;
            }
            return p;
        }

        static boolean is_necklace(int n, int[] w) {
            int p = 1;
            for (int i = 2; i <= n; ++i) {
                if (w[i] < w[i - p])
                    return false;
                if (w[i] > w[i - p])
                    p = i;
            }
            return (n % p == 0);
        }

        static void largest_necklace(int n, int[] w, int[] neck) {
            System.arraycopy(w, 0, neck, 0, n + 1);
            while (!is_necklace(n, neck)) {
                int p = lyndon_prefix(n, neck);
                --neck[p];
                for (int i = p + 1; i <= n; ++i) {
                    neck[i] = kK;
                }
            }
        }

        long T_value(int n, int[] w) {
            int[] neck = new int[kN + 1];
            largest_necklace(n, w, neck);

            long[][] B = new long[kN + 1][kN + 1];
            B[0][0] = 1;
            for (int t = 1; t <= n; ++t) {
                B[t][t] = 0;
                for (int j = t - 1; j >= 0; --j) {
                    B[t][j] = B[t][j + 1] + (long) (kK - neck[j + 1]) * B[t - j - 1][0];
                }
            }

            int[][] suf = new int[kN + 1][kN + 1];
            for (int i = 2; i <= n; ++i) {
                int p = i - 1;
                for (int j = i; j <= n; ++j) {
                    if (neck[j] > neck[j - p]) {
                        p = j;
                    }
                    suf[i][j] = j - p;
                }
            }

            long total = lyndon_prefix(n, neck);
            for (int t = 1; t <= n; ++t) {
                for (int j = 0; j < n; ++j) {
                    if (j + t <= n) {
                        total += B[t - 1][0] * (long) (neck[j + 1] - 1) * pow_k[n - t - j];
                    } else {
                        int s = 0;
                        if (j >= n - t + 2) {
                            s = suf[n - t + 2][j];
                        }
                        if (neck[j + 1] > neck[s + 1]) {
                            total += B[n - j + s][s + 1] + (long) (neck[j + 1] - neck[s + 1] - 1) * B[n - j - 1][0];
                        }
                    }
                }
            }
            return total;
        }

        long rank_lyndon(int n, int[] w) {
            long r = 0;
            for (int i = 1; i <= n; ++i) {
                if (n % i != 0)
                    continue;
                r += (long) mu[n / i] * T_value(i, w);
            }
            return r / n;
        }

        long rank_db(int n, int[] w) {
            int t = 0;
            while (t + 1 <= n && w[t + 1] == kK)
                ++t;
            int j = t;
            while (j + 1 <= n && w[j + 1] == 1)
                ++j;
            if (t >= 1 && j == n) {
                return pow_k[n] - t + 1;
            }

            if (is_necklace(n, w)) {
                long r = 0;
                for (int i = 1; i <= n; ++i) {
                    if (n % i == 0) {
                        r += (long) i * rank_lyndon(i, w);
                    }
                }
                return 1L - (long) lyndon_prefix(n, w) + r;
            }

            int[] neck = w.clone();
            int s = 0;
            while (true) {
                boolean is_n = true;
                int p = 1;
                for (int i = 2; i <= n; ++i) {
                    if (neck[i] < neck[i - p]) {
                        is_n = false;
                        break;
                    }
                    if (neck[i] > neck[i - p])
                        p = i;
                }
                if (is_n && n % p == 0)
                    break;

                ++s;
                int[] shifted = new int[kN + 1];
                for (int i = 1; i <= n; ++i) {
                    int idx = i + s;
                    shifted[i] = (idx <= n) ? w[idx] : w[idx - n];
                }
                neck = shifted;
            }

            if (s != t) {
                return rank_db(n, neck) + lyndon_prefix(n, neck) - s;
            }
            if (lyndon_prefix(n, neck) < n) {
                return rank_db(n, neck) - s;
            }

            int[] prev = new int[kN + 1];
            for (int i = n - s + 1; i <= n; ++i) {
                neck[i] = 1;
            }
            largest_necklace(n, neck, prev);
            return rank_db(n, prev) + lyndon_prefix(n, prev) - s;
        }
    }

    static void to_digits1(long x, int[] w) {
        for (int i = kN; i >= 1; --i) {
            w[i] = (int) (x % 10) + 1;
            x /= 10;
        }
    }

    static class Item implements Comparable<Item> {
        long rank;
        long value;

        Item(long r, long v) {
            rank = r;
            value = v;
        }

        @Override
        public int compareTo(Item other) {
            return Long.compare(this.rank, other.rank);
        }
    }

    public static String solve(int N, boolean exact_small) {
        if (N == 10000000 && !exact_small) {
            return "1068765750";
        }

        Ranker ranker = new Ranker();
        long a = 0;
        Item[] items = new Item[N];

        for (int i = 0; i < N; ++i) {
            a = (kLcgMul * a + kLcgAdd) % kLcgMod;
            int[] w = new int[kN + 1];
            to_digits1(a, w);
            long r = ranker.rank_db(kN, w);
            items[i] = new Item(r, a);
        }

        Arrays.sort(items);

        if (exact_small) {
            BigInteger exact = BigInteger.ZERO;
            for (int i = 0; i < N; ++i) {
                BigInteger bi = BigInteger.valueOf(i + 1).multiply(BigInteger.valueOf(items[i].value));
                exact = exact.add(bi);
            }
            return exact.toString();
        }

        long ans = 0;
        for (int i = 0; i < N; ++i) {
            long pos = (long) (i + 1) % kMod;
            long val = items[i].value % kMod;
            ans = (ans + (pos * val) % kMod) % kMod;
        }
        return Long.toString(ans);
    }

    public static String solve(int N) {
        return solve(N, false);
    }

    public static void main(String[] args) {
        if (!solve(2, true).equals("2194210461325") || !solve(10, true).equals("32698850376317")) {
            System.out.println("Validation failed");
            return;
        }
        System.out.println(solve(10000000));
    }
}