Problem 1001: Connections I

View on Project Euler

Project Euler Problem 1001 Solution

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

Problem Summary We are given an array of \(2n\) elements in which every value occurs exactly twice. The two positions of a value form a chord (an interval) over the row. The array is called connectable when all its chords can be drawn above the row without any two of them crossing. From one array we may form \(2^n\) new arrays by independently keeping or deleting both occurrences of each value. The connectivity number is the number of those \(2^n\) sub-arrays that are connectable. Equivalently, it counts the subsets of the \(n\) chords that are pairwise non-crossing. The required answer is this count modulo $$M=1\,003\,443\,221.$$ For example \([0,1,0,1]\) has connectivity number \(3\), and the \(20\)-element array \([0,1,2,3,1,4,0,5,4,2,6,7,3,8,6,5,9,8,9,7]\) has connectivity number \(86\). The target instance is a \(40\,000\)-element array, i.e. \(n=20\,000\) chords. The number of non-crossing subsets is astronomically large (up to \(2^n\)), so it is accumulated modulo \(M\); brute force over all \(2^n\) subsets is only used to validate the fast method on tiny inputs. Mathematical Approach Chords, crossings, and what is being counted If a value sits at positions \(p \lt q\), its chord is the interval \([p,q]\)....

Detailed mathematical approach

Problem Summary

We are given an array of \(2n\) elements in which every value occurs exactly twice. The two positions of a value form a chord (an interval) over the row. The array is called connectable when all its chords can be drawn above the row without any two of them crossing.

From one array we may form \(2^n\) new arrays by independently keeping or deleting both occurrences of each value. The connectivity number is the number of those \(2^n\) sub-arrays that are connectable. Equivalently, it counts the subsets of the \(n\) chords that are pairwise non-crossing. The required answer is this count modulo

$$M=1\,003\,443\,221.$$

For example \([0,1,0,1]\) has connectivity number \(3\), and the \(20\)-element array \([0,1,2,3,1,4,0,5,4,2,6,7,3,8,6,5,9,8,9,7]\) has connectivity number \(86\). The target instance is a \(40\,000\)-element array, i.e. \(n=20\,000\) chords.

The number of non-crossing subsets is astronomically large (up to \(2^n\)), so it is accumulated modulo \(M\); brute force over all \(2^n\) subsets is only used to validate the fast method on tiny inputs.

Mathematical Approach

Chords, crossings, and what is being counted

If a value sits at positions \(p \lt q\), its chord is the interval \([p,q]\). Two chords \([a,b]\) and \([c,d]\) with \(a \lt c\) cross exactly when they interleave,

$$a \lt c \lt b \lt d.$$

If they do not interleave they are either nested (one interval contains the other) or disjoint. A set of chords is drawable above the line without crossings precisely when no two of them cross, so the connectivity number equals the number of crossing-free subsets, including the empty subset. In graph terms, build the crossing graph whose vertices are the chords and whose edges join crossing pairs; the connectivity number is the number of independent sets of that graph.

Non-crossing families are laminar

Any pairwise non-crossing collection of chords is a laminar family: any two members are nested or disjoint, never interleaved. Laminar families are exactly the forests of the containment order, and it is this nesting structure — absent in a general graph — that turns the count into a polynomial-time dynamic program rather than an intractable independent-set count.

Why closing order is enough

The fast method does not explicitly build the whole crossing graph. A chord that opens inside \([\ell_i,r_i]\) and closes after \(r_i\) crosses chord \(i\), so it cannot appear in a subset that also contains \(i\). At the instant chord \(i\) closes, every chord fully nested inside it has already closed, while every chord that crosses it is still open. Thus the value already accumulated at slot \(i+1\) is exactly the count of admissible inner choices that can be combined with chord \(i\).

A dynamic program over the nesting structure

Sort the chords by left endpoint, \(\ell_0 \lt \ell_1 \lt \dots \lt \ell_{n-1}\) (all endpoints are distinct), and write \(r_i\) for the right endpoint of chord \(i\). For each chord define

$$\operatorname{next}(i)=\min\{\,j : \ell_j \gt r_i\,\},$$

the first chord lying entirely to the right of chord \(i\). The chords with index in \((i,\operatorname{next}(i))\) are precisely those that open inside chord \(i\); among them the nested ones close before \(r_i\) while the crossing ones close after \(r_i\).

The algorithm processes chords in order of increasing right endpoint and maintains an array \(\textit{ways}\) (initialised to all ones) together with a companion array \(\textit{delta}\). When chord \(i\) closes, every chord nested inside it has already closed, so the value

$$\textit{inside}(i)=\textit{ways}[i+1]$$

is the number of non-crossing subsets that live strictly inside chord \(i\). The subsets that include chord \(i\) are then obtained by freely combining one such inner configuration with any non-crossing configuration drawn from the chords entirely to its right:

$$\textit{inc}=\textit{inside}(i)\cdot \textit{ways}[\operatorname{next}(i)].$$

This increment is recorded in \(\textit{delta}[i]\) and added into \(\textit{ways}[i]\). It must also reach every earlier chord \(p \lt i\) that is still open and whose nested region has already been passed (\(\operatorname{next}(p)\le i\)): for such a \(p\), the freshly closed chord \(i\) extends the subsets that include \(p\) through the term \(\textit{inside}(p)\cdot \textit{delta}[\operatorname{next}(p)]\). Propagating this contribution down to all earlier slots takes \(O(n)\) work per closing chord, and after the last chord closes the accumulated total is

$$\text{connectivity number}=\textit{ways}[0]\bmod M.$$

The two array entries \(\textit{ways}[i]\) and \(\textit{delta}[i]\) thus encode, respectively, the running count of valid subsets anchored at slot \(i\) and the most recent increment to be forwarded — the bookkeeping that lets a single linear sweep per chord account for both the "nested inside" and "continues to the right" choices.

The invariant behind ways and delta

Think of slot \(s\) as the boundary just before chord \(s\) in left-endpoint order. After each closing event, \(\textit{ways}[s]\) contains all valid selections whose next available left endpoint is at or after slot \(s\), restricted to chords already closed by the sweep. The companion value \(\textit{delta}[s]\) stores only the newest contribution that still has to be forwarded to earlier open ancestors. This separation prevents double-counting: older contributions have already been absorbed into \(\textit{ways}\), while the loop over \(p=i-1,\dots,0\) propagates only the fresh \(\textit{delta}\) term created by the chord that just closed.

Worked example: \([0,1,0,1]\)

Value \(0\) sits at positions \(0,2\) and value \(1\) at positions \(1,3\), giving chords \([0,2]\) and \([1,3]\). Since \(0 \lt 1 \lt 2 \lt 3\) they interleave, so they cross. The crossing-free subsets are therefore \(\varnothing\), \(\{[0,2]\}\) and \(\{[1,3]\}\) — the full pair is excluded — for a connectivity number of \(3\), matching the stated value.

How the Code Works

The C++, Python, and Java implementations share the same pipeline. load_csv reads the comma-separated array. build_intervals records the two positions of each value, forms the chords \([p,q]\), sorts them by left endpoint, and asserts that the left endpoints are strictly increasing.

connectivity_number builds \(\operatorname{next}(i)\) with a binary search (upper_bound over the sorted left endpoints), orders the chords by right endpoint, and runs the sweep above using the \(\textit{ways}\), \(\textit{delta}\), \(\textit{inside}\), and \(\textit{active}\) arrays, with all arithmetic carried out modulo \(M\) through add_mod and mul_mod (the latter widening to \(128\) bits in C++; \(64\) bits already suffice in Java and Python since the factors are below \(M\)).

crosses and brute_connectivity form the validation path: the latter simply enumerates all \(2^n\) subsets and counts those with no crossing pair. run_checkpoints checks the published values \(3,8,5,8,86\) on five small arrays and confirms that the fast result equals the brute-force count reduced modulo \(M\); main then evaluates the \(40\,000\)-element instance.

Complexity Analysis

Building the chords costs \(O(n)\) plus \(O(n\log n)\) for the two sorts and the binary searches. The closing sweep performs an \(O(n)\) propagation for each of the \(n\) chords, so the dominant cost is

$$O(n^2)$$

time with \(O(n)\) memory. For \(n=20\,000\) this is a few hundred million constant-time modular operations, comfortably fast in C++. The brute-force checker is exponential, \(O(2^n\cdot n^2)\), and is restricted to the small checkpoint arrays only.

Footnotes and References

  1. Problem page: Project Euler 1001
  2. Circle graph (chord intersection graph): Wikipedia - Circle graph
  3. Laminar set family: Wikipedia - Laminar set family
  4. Independent set (graph theory): Wikipedia - Independent set
  5. Dynamic programming: Wikipedia - Dynamic programming
  6. Modular arithmetic: Wikipedia - Modular arithmetic

Problem 1001 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <fstream>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using i64 = std::int64_t;

constexpr i64 MOD = 1'003'443'221LL;

struct Interval {
    int left;
    int right;
};

std::vector<int> load_csv(const std::string& path) {
    std::ifstream fin(path);
    if (!fin) {
        throw std::runtime_error("cannot open input file");
    }

    std::string data;
    std::getline(fin, data);
    std::stringstream ss(data);
    std::string token;
    std::vector<int> values;

    while (std::getline(ss, token, ',')) {
        if (!token.empty()) {
            values.push_back(std::stoi(token));
        }
    }

    return values;
}

std::vector<Interval> build_intervals(const std::vector<int>& values) {
    std::unordered_map<int, std::vector<int>> positions;
    positions.reserve(values.size() / 2U + 1U);

    for (int i = 0; i < static_cast<int>(values.size()); ++i) {
        positions[values[static_cast<std::size_t>(i)]].push_back(i);
    }

    std::vector<Interval> intervals;
    intervals.reserve(positions.size());
    for (const auto& entry : positions) {
        assert(entry.second.size() == 2U);
        intervals.push_back({entry.second[0], entry.second[1]});
    }

    std::sort(intervals.begin(), intervals.end(), [](const Interval& a, const Interval& b) {
        return a.left < b.left;
    });

    for (std::size_t i = 1; i < intervals.size(); ++i) {
        assert(intervals[i - 1].left < intervals[i].left);
    }

    return intervals;
}

i64 add_mod(const i64 a, const i64 b, const i64 mod) {
    const i64 s = a + b;
    return s >= mod ? s - mod : s;
}

i64 mul_mod(const i64 a, const i64 b, const i64 mod) {
    return static_cast<i64>((static_cast<__int128>(a) * b) % mod);
}

i64 connectivity_number(const std::vector<int>& values, const i64 mod) {
    const std::vector<Interval> intervals = build_intervals(values);
    const int n = static_cast<int>(intervals.size());

    std::vector<int> left(static_cast<std::size_t>(n));
    std::vector<int> right(static_cast<std::size_t>(n));
    for (int i = 0; i < n; ++i) {
        left[static_cast<std::size_t>(i)] = intervals[static_cast<std::size_t>(i)].left;
        right[static_cast<std::size_t>(i)] = intervals[static_cast<std::size_t>(i)].right;
    }

    std::vector<int> next(static_cast<std::size_t>(n));
    for (int i = 0; i < n; ++i) {
        next[static_cast<std::size_t>(i)] = static_cast<int>(
            std::upper_bound(left.begin(), left.end(), right[static_cast<std::size_t>(i)]) - left.begin());
    }

    std::vector<int> by_right(static_cast<std::size_t>(n));
    for (int i = 0; i < n; ++i) {
        by_right[static_cast<std::size_t>(i)] = i;
    }
    std::sort(by_right.begin(), by_right.end(), [&](const int a, const int b) {
        return right[static_cast<std::size_t>(a)] < right[static_cast<std::size_t>(b)];
    });

    std::vector<i64> inside(static_cast<std::size_t>(n), 0);
    std::vector<i64> ways(static_cast<std::size_t>(n + 1), 1);
    std::vector<i64> delta(static_cast<std::size_t>(n + 1), 0);
    std::vector<unsigned char> active(static_cast<std::size_t>(n), 0);

    for (const int i : by_right) {
        inside[static_cast<std::size_t>(i)] = ways[static_cast<std::size_t>(i + 1)];

        const i64 inc = mul_mod(inside[static_cast<std::size_t>(i)],
                                ways[static_cast<std::size_t>(next[static_cast<std::size_t>(i)])],
                                mod);
        delta[static_cast<std::size_t>(i)] = inc;
        ways[static_cast<std::size_t>(i)] = add_mod(ways[static_cast<std::size_t>(i)], inc, mod);

        for (int p = i - 1; p >= 0; --p) {
            i64 d = delta[static_cast<std::size_t>(p + 1)];
            if (active[static_cast<std::size_t>(p)] != 0 &&
                next[static_cast<std::size_t>(p)] <= i) {
                d = add_mod(d,
                            mul_mod(inside[static_cast<std::size_t>(p)],
                                    delta[static_cast<std::size_t>(next[static_cast<std::size_t>(p)])],
                                    mod),
                            mod);
            }
            delta[static_cast<std::size_t>(p)] = d;
            ways[static_cast<std::size_t>(p)] = add_mod(ways[static_cast<std::size_t>(p)], d, mod);
        }

        active[static_cast<std::size_t>(i)] = 1;
    }

    return ways[0];
}

bool crosses(const Interval& a, const Interval& b) {
    Interval x = a;
    Interval y = b;
    if (y.left < x.left) {
        std::swap(x, y);
    }
    return x.left < y.left && y.left < x.right && x.right < y.right;
}

i64 brute_connectivity(const std::vector<int>& values) {
    const std::vector<Interval> intervals = build_intervals(values);
    const int n = static_cast<int>(intervals.size());
    assert(n <= 20);

    i64 total = 0;
    for (std::uint64_t mask = 0; mask < (1ULL << n); ++mask) {
        bool ok = true;
        for (int i = 0; i < n && ok; ++i) {
            if (((mask >> i) & 1ULL) == 0ULL) {
                continue;
            }
            for (int j = i + 1; j < n; ++j) {
                if (((mask >> j) & 1ULL) != 0ULL &&
                    crosses(intervals[static_cast<std::size_t>(i)],
                            intervals[static_cast<std::size_t>(j)])) {
                    ok = false;
                    break;
                }
            }
        }
        if (ok) {
            ++total;
        }
    }
    return total;
}

void run_checkpoints() {
    const std::vector<std::vector<int>> cases = {
        {0, 1, 0, 1},
        {0, 0, 1, 2, 2, 1},
        {0, 1, 2, 1, 0, 2},
        {0, 1, 2, 2, 1, 0},
        {0, 1, 2, 3, 1, 4, 0, 5, 4, 2, 6, 7, 3, 8, 6, 5, 9, 8, 9, 7},
    };

    const std::vector<i64> expected = {3, 8, 5, 8, 86};
    for (std::size_t i = 0; i < cases.size(); ++i) {
        const i64 brute = brute_connectivity(cases[i]);
        assert(brute == expected[i]);
        assert(connectivity_number(cases[i], MOD) == brute % MOD);
    }
}

}  // namespace

int main() {
    run_checkpoints();

    const std::vector<int> input = load_csv("resources/documents/1001_input.txt");
    assert(input.size() == 40'000U);

    std::cout << connectivity_number(input, MOD) << '\n';
    return 0;
}

Python

import bisect

MOD = 1_003_443_221


def load_csv(path):
    with open(path) as fin:
        data = fin.readline()
    return [int(token) for token in data.split(",") if token.strip() != ""]


def build_intervals(values):
    positions = {}
    for i, value in enumerate(values):
        positions.setdefault(value, []).append(i)

    intervals = []
    for occ in positions.values():
        assert len(occ) == 2
        intervals.append((occ[0], occ[1]))  # (left, right)

    intervals.sort(key=lambda iv: iv[0])
    for i in range(1, len(intervals)):
        assert intervals[i - 1][0] < intervals[i][0]
    return intervals


def connectivity_number(values, mod):
    intervals = build_intervals(values)
    n = len(intervals)

    left = [iv[0] for iv in intervals]
    right = [iv[1] for iv in intervals]

    # next[i] = first interval whose left endpoint lies strictly after right[i]
    nxt = [bisect.bisect_right(left, right[i]) for i in range(n)]

    by_right = sorted(range(n), key=lambda i: right[i])

    inside = [0] * n
    ways = [1] * (n + 1)
    delta = [0] * (n + 1)
    active = [0] * n

    for i in by_right:
        inside[i] = ways[i + 1]

        inc = inside[i] * ways[nxt[i]] % mod
        delta[i] = inc
        ways[i] = (ways[i] + inc) % mod

        for p in range(i - 1, -1, -1):
            d = delta[p + 1]
            if active[p] and nxt[p] <= i:
                d = (d + inside[p] * delta[nxt[p]]) % mod
            delta[p] = d
            ways[p] = (ways[p] + d) % mod

        active[i] = 1

    return ways[0]


def crosses(a, b):
    if b[0] < a[0]:
        a, b = b, a
    return a[0] < b[0] and b[0] < a[1] and a[1] < b[1]


def brute_connectivity(values):
    intervals = build_intervals(values)
    n = len(intervals)
    assert n <= 20

    total = 0
    for mask in range(1 << n):
        ok = True
        for i in range(n):
            if not (mask >> i) & 1:
                continue
            for j in range(i + 1, n):
                if (mask >> j) & 1 and crosses(intervals[i], intervals[j]):
                    ok = False
                    break
            if not ok:
                break
        if ok:
            total += 1
    return total


def run_checkpoints():
    cases = [
        [0, 1, 0, 1],
        [0, 0, 1, 2, 2, 1],
        [0, 1, 2, 1, 0, 2],
        [0, 1, 2, 2, 1, 0],
        [0, 1, 2, 3, 1, 4, 0, 5, 4, 2, 6, 7, 3, 8, 6, 5, 9, 8, 9, 7],
    ]
    expected = [3, 8, 5, 8, 86]
    for case, want in zip(cases, expected):
        brute = brute_connectivity(case)
        assert brute == want
        assert connectivity_number(case, MOD) == brute % MOD


if __name__ == "__main__":
    run_checkpoints()

    values = load_csv("resources/documents/1001_input.txt")
    assert len(values) == 40_000

    print(connectivity_number(values, MOD))

Java

import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Euler1001 {
    static final long MOD = 1_003_443_221L;

    static int[] loadCsv(String path) throws IOException {
        String data;
        try (BufferedReader reader = new BufferedReader(new FileReader(path))) {
            data = reader.readLine();
        }
        List<Integer> values = new ArrayList<>();
        for (String token : data.split(",")) {
            if (!token.trim().isEmpty()) {
                values.add(Integer.parseInt(token.trim()));
            }
        }
        int[] result = new int[values.size()];
        for (int i = 0; i < result.length; ++i) {
            result[i] = values.get(i);
        }
        return result;
    }

    // intervals[i] = { left, right }, sorted by left
    static int[][] buildIntervals(int[] values) {
        Map<Integer, int[]> positions = new HashMap<>();
        Map<Integer, Integer> seen = new HashMap<>();
        for (int i = 0; i < values.length; ++i) {
            int v = values[i];
            int c = seen.getOrDefault(v, 0);
            if (c == 0) {
                positions.put(v, new int[] { i, -1 });
            } else {
                assert c == 1;
                positions.get(v)[1] = i;
            }
            seen.put(v, c + 1);
        }

        List<int[]> intervals = new ArrayList<>();
        for (int[] pair : positions.values()) {
            assert pair[1] != -1;
            intervals.add(pair);
        }
        intervals.sort((a, b) -> Integer.compare(a[0], b[0]));
        for (int i = 1; i < intervals.size(); ++i) {
            assert intervals.get(i - 1)[0] < intervals.get(i)[0];
        }
        return intervals.toArray(new int[0][]);
    }

    // first index whose left endpoint is strictly greater than key (upper_bound)
    static int upperBound(int[] sortedLeft, int key) {
        int lo = 0;
        int hi = sortedLeft.length;
        while (lo < hi) {
            int mid = (lo + hi) >>> 1;
            if (sortedLeft[mid] <= key) {
                lo = mid + 1;
            } else {
                hi = mid;
            }
        }
        return lo;
    }

    static long mulMod(long a, long b, long mod) {
        return a * b % mod;
    }

    static long addMod(long a, long b, long mod) {
        long s = a + b;
        return s >= mod ? s - mod : s;
    }

    static long connectivityNumber(int[] values, long mod) {
        int[][] intervals = buildIntervals(values);
        int n = intervals.length;

        int[] left = new int[n];
        int[] right = new int[n];
        for (int i = 0; i < n; ++i) {
            left[i] = intervals[i][0];
            right[i] = intervals[i][1];
        }

        int[] next = new int[n];
        for (int i = 0; i < n; ++i) {
            next[i] = upperBound(left, right[i]);
        }

        Integer[] byRight = new Integer[n];
        for (int i = 0; i < n; ++i) {
            byRight[i] = i;
        }
        Arrays.sort(byRight, (a, b) -> Integer.compare(right[a], right[b]));

        long[] inside = new long[n];
        long[] ways = new long[n + 1];
        Arrays.fill(ways, 1L);
        long[] delta = new long[n + 1];
        boolean[] active = new boolean[n];

        for (int i : byRight) {
            inside[i] = ways[i + 1];

            long inc = mulMod(inside[i], ways[next[i]], mod);
            delta[i] = inc;
            ways[i] = addMod(ways[i], inc, mod);

            for (int p = i - 1; p >= 0; --p) {
                long d = delta[p + 1];
                if (active[p] && next[p] <= i) {
                    d = addMod(d, mulMod(inside[p], delta[next[p]], mod), mod);
                }
                delta[p] = d;
                ways[p] = addMod(ways[p], d, mod);
            }

            active[i] = true;
        }

        return ways[0];
    }

    static boolean crosses(int[] a, int[] b) {
        int al = a[0], ar = a[1], bl = b[0], br = b[1];
        if (bl < al) {
            int tl = al, tr = ar;
            al = bl; ar = br; bl = tl; br = tr;
        }
        return al < bl && bl < ar && ar < br;
    }

    static long bruteConnectivity(int[] values) {
        int[][] intervals = buildIntervals(values);
        int n = intervals.length;
        assert n <= 20;

        long total = 0;
        for (long mask = 0; mask < (1L << n); ++mask) {
            boolean ok = true;
            for (int i = 0; i < n && ok; ++i) {
                if (((mask >> i) & 1L) == 0L) {
                    continue;
                }
                for (int j = i + 1; j < n; ++j) {
                    if (((mask >> j) & 1L) != 0L && crosses(intervals[i], intervals[j])) {
                        ok = false;
                        break;
                    }
                }
            }
            if (ok) {
                ++total;
            }
        }
        return total;
    }

    static void runCheckpoints() {
        int[][] cases = {
            { 0, 1, 0, 1 },
            { 0, 0, 1, 2, 2, 1 },
            { 0, 1, 2, 1, 0, 2 },
            { 0, 1, 2, 2, 1, 0 },
            { 0, 1, 2, 3, 1, 4, 0, 5, 4, 2, 6, 7, 3, 8, 6, 5, 9, 8, 9, 7 },
        };
        long[] expected = { 3, 8, 5, 8, 86 };
        for (int i = 0; i < cases.length; ++i) {
            long brute = bruteConnectivity(cases[i]);
            assert brute == expected[i];
            assert connectivityNumber(cases[i], MOD) == brute % MOD;
        }
    }

    public static void main(String[] args) throws IOException {
        runCheckpoints();

        int[] values = loadCsv("resources/documents/1001_input.txt");
        assert values.length == 40_000;

        System.out.println(connectivityNumber(values, MOD));
    }
}