Problem 186: Connectedness of a Network

Project Euler Problem 186 Solution

EulerSolve provides an optimized solution for Project Euler Problem 186, Connectedness of a Network, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary There are $N=1000000$ subscribers labeled $0,1,\dots,999999$. The distinguished subscriber is the Prime Minister, whose number is $524287$. Call attempts are generated by the deterministic sequence $$S_k=(100003-200003k+300007k^3)\bmod 1000000 \qquad (1\le k\le 55),$$ $$S_k=(S_{k-24}+S_{k-55})\bmod 1000000 \qquad (k\ge 56).$$ The $n$-th call attempt uses the ordered pair $(S_{2n-1},S_{2n})$. If the two numbers are equal, the call is a misdial and does not count. Otherwise it is a successful call and creates an undirected connection between those two subscribers. The goal is to find the first successful-call count $t$ for which the Prime Minister belongs to a connected component of size at least $99\%$ of the network, namely $990000$ subscribers. Mathematical Approach The graph only grows: successful calls add edges, and nothing ever removes one. Because of that, the problem is really about tracking a partition of the subscriber set into connected components, not about storing the whole history of calls. The call stream is deterministic Although the story describes a pseudo-random telephone log, the data are completely determined by the recurrence for $S_k$. Once the first 55 values are fixed, every later value follows from the lagged relation $S_{k-24}+S_{k-55}$....

Detailed mathematical approach

Problem Summary

A telephone system has one million users with the numbers $0,1,\dots,999999$. Record $n$ of its log says that user $S_{2n-1}$ called user $S_{2n}$, where a lagged Fibonacci generator gives, for $1\le k\le55$,

$$\begin{aligned}S_k=(&100003-200003k\\&{}+300007k^3)\bmod 10^6,\end{aligned}$$

and for $k\ge56$,

$$S_k=(S_{k-24}+S_{k-55})\bmod 10^6.$$

A record with two equal numbers is a misdial and fails; every other record is a successful call. Two users are friends if one has called the other, and the relation extends along chains: a friend of a friend counts, and so on. We must find after how many successful calls, misdials not counted, at least $99\%$ of the users, that is $990000$ users including the Prime Minister himself, are linked by such chains to the Prime Minister, whose number is $524287$. The statement lists the first three records: $200007$ calls $100053$, $600183$ calls $500439$, and $600863$ calls $701497$. ^[1]

The whole solution in one sentence: successful calls are lines between users, and groups joined by lines only merge, so we replay the calls in order, keep the groups in a union–find structure that merges groups and reports their sizes in a few steps, and stop at the first call after which the Prime Minister's group has at least $990000$ members.

Mathematical Approach

We turn friendship into graph language, prove that calls only merge groups, replay the log with union–find, and compare the replay with random graph theory.

1. The call records

The first $55$ numbers come from the cubic. For $k=2$ it gives $100003-400006+2400056$${}=2100053$, so $S_2=100053$; in the same way $S_1=200007$ and $S_3=600183$, and $S_4$, $S_5$, $S_6$ complete the three records of the statement. At $k=55$ the cubic reaches $49902764463$, more than $2^{31}$, so it needs $64$-bit integers before the reduction. From $k=56$ on, each number is the sum of the numbers $24$ and $55$ places back, for example

$$\begin{aligned}S_{56}&=S_{32}+S_1\\&=329283+200007\\&=529290,\end{aligned}$$

and the reduction of the sum first matters at $S_{58}=S_{34}+S_3-10^6=375212$. This is an additive lagged Fibonacci generator with the classical lags $24$ and $55$. ^[2] Only the last $55$ numbers are ever needed. If $S_k$ is stored in slot $(k-1)\bmod 55$, this slot still holds $S_{k-55}$ when $S_k$ is due, and $S_{k-55}$ is then used for the last time; $S_{k-24}$ lies $31$ slots further on, because $-24\equiv31\pmod{55}$.

2. Friends are connected components

Draw each user as a point and each successful call as a line between its two users. After $t$ successful calls this gives a graph $G_t$ with $n=10^6$ vertices and $t$ edges, in which a repeated call is one more parallel line, and a chain of friends is a path in it. Being joined by a path, with every user joined to himself, is an equivalence relation, since paths can be walked backwards and joined end to end; its classes are the connected components of $G_t$. ^[3] With $C_t(x)$ the component of user $x$, the task asks for the smallest $t$ with

$$|C_t(524287)|\ge990000,$$

where the Prime Minister counts himself. A misdial, a line from a user to himself, would link nobody new, and it is not counted.

3. A toy log

Take $8$ users $0,1,\dots,7$, let user $0$ be the Prime Minister and stop when his component has $6$ users, $75\%$ of all. The records are $(0,1)$, $(2,3)$, $(4,4)$, $(1,2)$, $(3,0)$, $(5,6)$, $(6,5)$ and $(2,5)$. Calls $1$ and $2$ form the groups $\{0,1\}$ and $\{2,3\}$, and the record $(4,4)$ is a misdial. Call $3$, the record $(1,2)$, merges the two groups into $\{0,1,2,3\}$. Call $4$, $(3,0)$, joins two users who are already linked and changes nothing. Call $5$ forms $\{5,6\}$, and call $6$, $(6,5)$, again changes nothing. Call $7$, $(2,5)$, gives the Prime Minister a component of $6$ users. The answer for the toy log is $7$: the misdial is not counted, but the two calls that changed nothing are, because the question counts successful calls, not merges.

4. What one call does

Claim. Let the $t$-th successful call join $u$ and $v$. If $u$ and $v$ lie in one component of $G_{t-1}$, then $G_t$ has the same components. Otherwise the components of $G_t$ are those of $G_{t-1}$, except that $C_{t-1}(u)$ and $C_{t-1}(v)$ are replaced by their union.

Proof. Every path of $G_{t-1}$ is a path of $G_t$, and the new edge links $C_{t-1}(u)$ to $C_{t-1}(v)$. Conversely, take a path of $G_t$ from $x$ to $y$ that uses the new edge. The part before its first use links $x$ to $u$ or $v$ inside $G_{t-1}$, and the part after its last use does the same for $y$, so both lie in $C_{t-1}(u)\cup C_{t-1}(v)$. Hence two users are linked in $G_t$ exactly when they were linked before or both lie in this union, which is the claim, since for linked $u$ and $v$ the union is one old component. $\square$

So the partition into components is the whole state; the calls that created it no longer matter.

5. The Prime Minister's component only grows

By Section 4, $C_{t-1}(x)\subseteq C_t(x)$ for every user $x$, so $|C_t(524287)|$ never decreases. Once $990000$ users are linked to the Prime Minister they stay linked, and the requested number is

$$\begin{aligned}T=\min\{\,t:{}&|C_t(524287)|\\&\quad{}\ge990000\,\}.\end{aligned}$$

Replaying the calls in order and testing the size after each one finds $T$; since the size never decreases, a binary search over $t$ would find it too, which the verification below uses.

6. Union–find

A union–find structure, or disjoint-set forest, stores every component as a rooted tree on its users. ^[4] Each user $x$ has a pointer $\mathrm{parent}(x)$; a root points to itself and represents its component, and $\mathrm{size}(r)$ at a root $r$ counts the users of its tree. Find follows the pointers from $x$ to its root, so two users are linked exactly when they have the same root. Unite handles a successful call $(a,b)$: it finds both roots, does nothing if they are equal, and otherwise points the root of the smaller tree to the root of the larger one and adds the sizes, which is called union by size. At the start every user is a root of size $1$, as in $G_0$, so by Section 4 and induction on $t$ the trees are always the components of $G_t$, and $|C_t(x)|=\mathrm{size}(\mathrm{find}(x))$.

In the toy log, calls $1$ and $2$ set $\mathrm{parent}(1)=0$ and $\mathrm{parent}(3)=2$. At call $3$ both trees have size $2$, and on a tie the programs keep the caller's root: $\mathrm{parent}(2)=0$ and $\mathrm{size}(0)=4$, although $2$ and $0$ never called each other, since a pointer records membership, not a call. Call $4$ walks from $3$ via $2$ to $0$, the root of $0$ itself, and merges nothing. Call $5$ sets $\mathrm{parent}(6)=5$, call $6$ finds the common root $5$, and call $7$ points the root $5$ of the smaller tree to the root $0$, so that $\mathrm{size}(0)=6$.

7. Why the trees stay shallow

Claim. With union by size no user is ever more than $\log_2 n$ steps below its root; for $n=10^6$ the depth is at most $19$.

Proof. The depth of $x$ grows only when the root of its tree is pointed to another root, and by the size rule that other tree is at least as large, so the tree containing $x$ at least doubles. It starts with one user and never has more than $n$, so this happens at most $\lfloor\log_2 n\rfloor$ times, and $\lfloor\log_2 10^6\rfloor=19$. $\square$

In the toy log, user $6$ reaches depth $2$: its tree grows from $1$ to $2$ users at call $5$ and from $2$ to $6$ at call $7$. The programs also shorten the paths they walk: path compression (C++) points every user on the walk straight to the root, path halving (Python, Java) points every other user on the walk to its grandparent, and in call $4$ of the toy log both let $3$ point straight to $0$. Such changes only move pointers to ancestors and never increase a depth, so the claim stays true. Tarjan proved, and with van Leeuwen extended to path halving, that $m$ operations then take time $O(m\,\alpha(n))$, where the inverse Ackermann function $\alpha(n)$ is at most $4$ for every $n$ that could ever be written down. ^[4]

8. Random calls and the giant component

If the calls were independent random pairs, $G_t$ with its rare repeated pairs merged would be the random graph $G(n,M)$ of Erdős and Rényi with $M\approx t$ edges, in which a user has $c=2t/n$ partners on average. For fixed $c$, with probability tending to $1$ as $n$ grows, the components change sharply at $c=1$, after $n/2$ calls: for $c\lt1$ every component has $O(\log n)$ users, at $c=1$ the largest has the order $n^{2/3}$, here $10^4$, and for $c\gt1$ a single giant component holds a fraction close to $y$ of the users, the positive solution of ^[5] ^[6]

$$y=1-e^{-cy}.$$

Heuristically, the number $D$ of partners of a user is nearly Poisson with mean $c$, that is $D=d$ with probability $e^{-c}c^d/d!$, since each call involves the user with probability $2/n$. ^[7] Treat the partners as independent branches, each leading into the giant component with probability about $y$; the user stays outside when no branch does, so

$$1-y=\sum_{d\ge0}e^{-c}\frac{c^d}{d!}(1-y)^d=e^{-cy}.$$

For $y$ near $1$ this is close to $e^{-c}$, the chance of having no partner at all, so late in the process most outsiders have never made or received a call.

9. What the real log looks like

The generated calls behave like random ones: the first of the rare misdials is record $442959$, in which $977223$ calls itself. Up to call $461514$ every successful call merges two components; call $461515$ is the first between two users who are already linked. The largest component has $150$ users after $400000$ calls and $5040$ after $500000$, where $c=1$. After $600000$, $700000$ and $10^6$ calls it has $315111$, $512114$ and $796599$ users, while the equation predicts the fractions $0.3137$, $0.5110$ and $0.7968$.

The Prime Minister stays out of all this for a long time: his number is absent from the first $1840580$ records, $3681160$ numbers, which for random numbers has probability $(1-10^{-6})^{3681160}\approx e^{-3.68}\approx0.025$. In record $1840581$, the successful call $1840579$, user $910544$ of the giant component calls him, and his component jumps from $1$ to $972133$ users, close to the predicted fraction $0.9721$ for $c\approx3.68$; every threshold from $1\%$ to $97\%$ is reached at this call, and the C++ program with --threshold-percent=50 prints $1840579$. Of the $27867$ users then outside, $25193$ have never been in a successful call, against $n\,e^{-c}\approx25194$ predicted, so the $99\%$ mark, at most $10000$ users outside, is reached mainly through first calls of such users.

How the Code Works

The C++ program ^[8] starts in main, where parse_arguments reads the options --subscribers=, --pm= and --threshold-percent=, with the defaults $10^6$, $524287$ and $99$, and --skip-checkpoints, and run_checkpoints compares the first six numbers with the statement. solve creates a LaggedFibonacci generator and a DSU, sets the target subscribers * threshold_percent / 100, here $990000$, and runs Section 5 literally: while the size at the root of the Prime Minister is below the target, it draws a caller and a called number, skips the record if they are equal, and otherwise counts the call and calls unite. The count is printed. The integer division rounds down; this is exact here, but for $3$ users and $99\%$ it gives the target $2$ instead of $3$.

LaggedFibonacci::next evaluates the cubic in std::int64_t for its first $55$ values and then uses the window of Section 1: idx_ is the slot of $S_{k-55}$, (idx_ + 31) % 55 that of $S_{k-24}$, and the new value overwrites slot idx_. In DSU, find is recursive with full path compression, safe because Section 7 limits the recursion to $20$ levels, and unite attaches the smaller tree below the larger one, keeping the caller's root on a tie.

The Python and Java programs run the same loop for the fixed parameters, without options or checkpoint test. They keep every generated number in a growing list and form a new number from the entries $24$ and $55$ places from its end, buf[-24] and buf[-55] in Python; find is a loop with path halving. Python evaluates the cubic with unbounded integers, Java with long constants such as 300007L. Java unboxes the entries of its ArrayList<Integer> into int variables, so caller == called compares numbers, not Integer objects.

Complexity and Verification

Let $A$ be the number of records read before the stop. They cost $2A$ generator steps and $A+1$ finds for the loop test, the last one included, and each successful call adds two finds in unite. By Section 7 the time is $O(n+A\log n)$, including the set-up of the arrays, and Tarjan's bound improves it to $O(n+A\,\alpha(n))$. C++ needs $O(n)$ memory, two arrays of $10^6$ entries of $32$ bits ($8$ MB) and a window of $55$ numbers; Java and Python also keep every generated number and need $O(n+A)$, and the Python lists hold references to integer objects. Measured once on a desktop computer, the C++ program compiled with -O2 takes about $0.05$ seconds, Java about $0.15$ and Python about $1.3$.

Besides the C++ checkpoint test, we checked the following for this article, outside the published solutions:

the $55$-slot window agrees with a plain array built from the definition for the first $10^7$ numbers;
a brute force that recomputes the Prime Minister's component by breadth-first search after every call gives $7$ for the toy log and matches the C++ program, whose generator then works modulo $n$, in $64$ cases: $n=100$, $300$, $1000$ or $2000$ users, the Prime Minister at $0$, $7$, $n-1$ or $524287\bmod n$, and the thresholds $10$, $50$, $90$ and $99\%$;
the three programs print the same number, and so does an independent program that builds the numbers from the definition, computes the components of the first $t$ calls with a standard routine and finds the first $t$ with $990000$ users around the Prime Minister by binary search, valid by Section 5; the first call after which any component has $990000$ users is the same;
union by size without path shortening reaches depth at most $9$ in the whole run, and the C++ forest ends with depth at most $3$;
only $3$ misdials occur before the stop, so counting records instead of successful calls would be off by exactly $3$.

Footnotes and References

The references supply the background facts; the merge rule, the depth bound and all numbers about the real log were derived or computed for this article.

Project Euler 186 — Connectedness of a Network. The official statement defines the records through the lagged Fibonacci generator, misdials and friendship through chains of calls, and asks after how many successful calls 99% of the one million users, the Prime Minister 524287 included, are linked to him.
Wikipedia — Lagged Fibonacci generator. Additive generators Sₙ ≡ Sₙ₋ⱼ + Sₙ₋ₖ (mod m) with 0 < j < k; (24, 55) is among the popular lags and heads Knuth's list in volume 2 of The Art of Computer Programming.
Wikipedia — Component (graph theory). Components are the equivalence classes of reachability, and a disjoint-set structure tracks them as edges are added.
Wikipedia — Disjoint-set data structure. Parent-pointer forests with union by size, path compression and path halving; Tarjan first proved the bound O(m α(n)) in 1975, path halving keeps this complexity, and α(n) ≤ 4 for any n that can be written in the physical universe.
Wikipedia — Giant component. For fixed ε > 0 and with probability tending to 1 as n grows, all components have size O(log n) when the edge probability is p ≤ (1 − ε)/n, and there is a single giant component when p ≥ (1 + ε)/n; at p = 1/n the largest component has size of order n^(2/3).
Wikipedia — Erdős–Rényi model. G(n, M) and G(n, p) behave alike when M ≈ (n choose 2)·p; the percolation section gives the relative size P∞ = p′(1 − exp(−⟨k⟩P∞)) of the giant component after a random fraction 1 − p′ of the nodes is removed, which is y = 1 − e^(−cy) for p′ = 1.
Wikipedia — Poisson distribution. B(n, p) is close to the Poisson distribution with mean np for large n and small p, whose generating function is exp(np(z − 1)).
C++ — Euler186.cpp. The linked immutable revision contains LaggedFibonacci, DSU, solve and run_checkpoints; it is the implementation record for this article.

Problem 186 source code

C++

#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>

namespace {

struct Options {
    int subscribers = 1000000;
    int pm = 524287;
    int threshold_percent = 99;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    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, "--subscribers=", options.subscribers) ||
            parse_int_after_prefix(arg, "--pm=", options.pm) ||
            parse_int_after_prefix(arg, "--threshold-percent=", options.threshold_percent)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.subscribers >= 2 && options.pm >= 0 && options.pm < options.subscribers &&
           options.threshold_percent > 0 && options.threshold_percent <= 100;
}

class LaggedFibonacci {
public:
    explicit LaggedFibonacci(const int mod) : mod_(mod), idx_(0), generated_(0) {
        values_.assign(55, 0);
    }

    int next() {
        ++generated_;
        if (generated_ <= 55) {
            const std::int64_t k = static_cast<std::int64_t>(generated_);
            const std::int64_t value = 100003LL - 200003LL * k + 300007LL * k * k * k;
            const int result = static_cast<int>((value % mod_ + mod_) % mod_);
            values_[static_cast<std::size_t>(generated_ - 1)] = result;
            return result;
        }

        const int a = values_[idx_];
        const int b = values_[(idx_ + 31) % 55];  // k-24 in a 55-slot ring
        const int result = (a + b) % mod_;
        values_[idx_] = result;
        idx_ = (idx_ + 1) % 55;
        return result;
    }

private:
    int mod_;
    int idx_;
    int generated_;
    std::vector<int> values_;
};

struct DSU {
    std::vector<int> parent;
    std::vector<int> size;

    explicit DSU(const int n) : parent(static_cast<std::size_t>(n)), size(static_cast<std::size_t>(n), 1) {
        std::iota(parent.begin(), parent.end(), 0);
    }

    int find(const int x) {
        if (parent[static_cast<std::size_t>(x)] == x) {
            return x;
        }
        parent[static_cast<std::size_t>(x)] = find(parent[static_cast<std::size_t>(x)]);
        return parent[static_cast<std::size_t>(x)];
    }

    void unite(const int a, const int b) {
        int ra = find(a);
        int rb = find(b);
        if (ra == rb) {
            return;
        }
        if (size[static_cast<std::size_t>(ra)] < size[static_cast<std::size_t>(rb)]) {
            std::swap(ra, rb);
        }
        parent[static_cast<std::size_t>(rb)] = ra;
        size[static_cast<std::size_t>(ra)] += size[static_cast<std::size_t>(rb)];
    }
};

int solve(const int subscribers, const int pm, const int threshold_percent) {
    LaggedFibonacci gen(subscribers);
    DSU dsu(subscribers);

    const int target_size = subscribers * threshold_percent / 100;
    int successful_calls = 0;
    while (dsu.size[static_cast<std::size_t>(dsu.find(pm))] < target_size) {
        const int caller = gen.next();
        const int called = gen.next();
        if (caller == called) {
            continue;
        }
        ++successful_calls;
        dsu.unite(caller, called);
    }
    return successful_calls;
}

bool run_checkpoints() {
    LaggedFibonacci gen(1000000);
    const int s1 = gen.next();
    const int s2 = gen.next();
    const int s3 = gen.next();
    const int s4 = gen.next();
    const int s5 = gen.next();
    const int s6 = gen.next();
    if (s1 != 200007 || s2 != 100053 || s3 != 600183 ||
        s4 != 500439 || s5 != 600863 || s6 != 701497) {
        std::cerr << "Checkpoint failed for Lagged Fibonacci prefix" << '\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 2;
    }
    std::cout << solve(options.subscribers, options.pm, options.threshold_percent) << '\n';
    return 0;
}

Python

# Problem 186: Connectedness of a network.
# Union-Find with Lagged Fibonacci generator.

def solve():
    subscribers = 1000000
    pm = 524287
    threshold = 99  # percent

    # Lagged Fibonacci
    s = [0]*(56)
    for k in range(1, 56):
        s[k] = (100003 - 200003*k + 300007*k*k*k) % 1000000
    buf = s[1:]
    idx = 0
    def next_val():
        nonlocal idx, buf
        if idx < 55:
            v = buf[idx]
        else:
            buf.append((buf[-24] + buf[-55]) % 1000000)
            v = buf[idx]
        idx += 1
        return v

    # Union-Find
    parent = list(range(subscribers))
    size = [1]*subscribers
    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x
    def unite(a, b):
        a, b = find(a), find(b)
        if a == b: return
        if size[a] < size[b]: a, b = b, a
        parent[b] = a
        size[a] += size[b]

    target_size = subscribers * threshold // 100
    successful = 0
    while size[find(pm)] < target_size:
        caller = next_val()
        called = next_val()
        if caller == called: continue
        successful += 1
        unite(caller, called)
    print(successful)

solve()

Java

public class Euler186 {
    static int[] parent, size;

    static int find(int x) {
        while (parent[x] != x) {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        return x;
    }

    static void unite(int a, int b) {
        a = find(a);
        b = find(b);
        if (a == b)
            return;
        if (size[a] < size[b]) {
            int t = a;
            a = b;
            b = t;
        }
        parent[b] = a;
        size[a] += size[b];
    }

    public static void main(String[] args) {
        int N = 1000000, pm = 524287;
        parent = new int[N];
        size = new int[N];
        for (int i = 0; i < N; i++) {
            parent[i] = i;
            size[i] = 1;
        }
        int[] s = new int[56];
        for (int k = 1; k <= 55; k++)
            s[k] = (int) ((100003L - 200003L * k + 300007L * k * k * k) % 1000000);
        int[] buf = new int[55];
        for (int i = 0; i < 55; i++)
            buf[i] = s[i + 1];
        int idx = 0;
        int bufLen = 55;
        int target = N * 99 / 100;
        int succ = 0;
        java.util.ArrayList<Integer> bufList = new java.util.ArrayList<>();
        for (int b : buf)
            bufList.add(b);
        while (size[find(pm)] < target) {
            int caller, called;
            if (idx < bufList.size()) {
                caller = bufList.get(idx);
            } else {
                bufList.add((bufList.get(bufList.size() - 24) + bufList.get(bufList.size() - 55)) % 1000000);
                caller = bufList.get(idx);
            }
            idx++;
            if (idx - 1 < bufList.size() - 1 && idx < bufList.size()) {
                called = bufList.get(idx);
            } else {
                bufList.add((bufList.get(bufList.size() - 24) + bufList.get(bufList.size() - 55)) % 1000000);
                called = bufList.get(idx);
            }
            idx++;
            if (caller == called)
                continue;
            succ++;
            unite(caller, called);
        }
        System.out.println(succ);
    }
}