Problem 430: Range Flips

Project Euler Problem 430 Solution

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

Problem Summary We have $N$ disks, all initially white. One move chooses two endpoints $A,B \in \{1,\dots,N\}$ independently and uniformly, then flips every disk in the closed interval between them. If a disk is flipped an odd number of times it ends black; if it is flipped an even number of times it ends white. The quantity of interest is $E(N,M)$, the expected number of white disks after $M$ moves. For the target input $N=10^{10}$ and $M=4000$, a direct simulation or a direct $O(N)$ summation is far too slow, so the solution reduces the whole problem to a one-disk probability formula plus a carefully bounded truncation of the final sum. Mathematical Approach 1. One-Turn Behavior of a Fixed Disk Fix a disk position $i$. In one move, disk $i$ is not flipped exactly when both chosen endpoints lie strictly to its left or both lie strictly to its right. Because ordered pairs $(A,B)$ are chosen uniformly from $N^2$ possibilities, the one-turn no-flip probability is $$q_i=\frac{(i-1)^2+(N-i)^2}{N^2}.$$ Therefore the one-turn flip probability is $$p_i=1-q_i.$$ It is convenient to combine these into the signed coefficient $$r_i=q_i-p_i=2q_i-1=\frac{2\big((i-1)^2+(N-i)^2\big)-N^2}{N^2}.$$ This is exactly the quantity raised to the $M$-th power by the implementations. 2....

Detailed mathematical approach

Problem Summary

We have $N$ disks, all initially white. One move chooses two endpoints $A,B \in \{1,\dots,N\}$ independently and uniformly, then flips every disk in the closed interval between them. If a disk is flipped an odd number of times it ends black; if it is flipped an even number of times it ends white.

The quantity of interest is $E(N,M)$, the expected number of white disks after $M$ moves. For the target input $N=10^{10}$ and $M=4000$, a direct simulation or a direct $O(N)$ summation is far too slow, so the solution reduces the whole problem to a one-disk probability formula plus a carefully bounded truncation of the final sum.

Mathematical Approach

1. One-Turn Behavior of a Fixed Disk

Fix a disk position $i$. In one move, disk $i$ is not flipped exactly when both chosen endpoints lie strictly to its left or both lie strictly to its right. Because ordered pairs $(A,B)$ are chosen uniformly from $N^2$ possibilities, the one-turn no-flip probability is

$$q_i=\frac{(i-1)^2+(N-i)^2}{N^2}.$$

Therefore the one-turn flip probability is

$$p_i=1-q_i.$$

It is convenient to combine these into the signed coefficient

$$r_i=q_i-p_i=2q_i-1=\frac{2\big((i-1)^2+(N-i)^2\big)-N^2}{N^2}.$$

This is exactly the quantity raised to the $M$-th power by the implementations.

2. From Flip Parity to the Probability of White

For disk $i$, introduce a sign variable $S_t \in \{+1,-1\}$, where $+1$ means white after $t$ moves and $-1$ means black. Initially $S_0=+1$. Each move multiplies the sign by $+1$ with probability $q_i$ and by $-1$ with probability $p_i$, so

$$\mathbb{E}[S_{t+1}\mid S_t]=r_i\,S_t.$$

Taking expectations repeatedly gives the simple recurrence

$$\mathbb{E}[S_t]=r_i^t,$$

because $\mathbb{E}[S_0]=1$. Now the white-indicator of disk $i$ after $M$ moves is

$$\mathbf{1}_{i,\text{white}}=\frac{1+S_M}{2},$$

hence

$$\Pr(\text{disk } i \text{ is white after } M \text{ moves})=\frac{1+r_i^M}{2}.$$

This is the core closed form. It converts a random sequence of interval flips into a deterministic per-position contribution.

3. Summing All Disks by Linearity of Expectation

Expected values add even when the disk colors are not independent. Therefore

$$E(N,M)=\sum_{i=1}^{N}\Pr(\text{disk } i \text{ is white after } M \text{ moves})$$

and the previous formula yields

$$E(N,M)=\sum_{i=1}^{N}\frac{1+r_i^M}{2}=\frac{N}{2}+\frac{1}{2}\sum_{i=1}^{N} r_i^M.$$

So the entire problem is reduced to evaluating the sum of $r_i^M$ over all positions.

4. Symmetry and Monotonicity Near the Edges

The coefficients are symmetric:

$$r_i=r_{N+1-i}.$$

This follows directly from the formula, since replacing $i$ by $N+1-i$ swaps the two squared terms.

On the left half of the board, the sequence is monotone decreasing. A direct subtraction gives

$$r_{i+1}-r_i=\frac{4(2i-N)}{N^2}.$$

Hence $r_{i+1} \lt r_i$ whenever $i \lt N/2$. The importance of this monotonicity is algorithmic: once a threshold is chosen, the last significant edge position can be found by binary search instead of by scanning all the way from the edge to the center.

For large $M$, any value with $\lvert r_i\rvert \lt 1$ becomes tiny after taking the $M$-th power. Since the disks in the middle have coefficients much farther from $1$ than the disks near the edges, almost the entire contribution comes from the two boundaries.

5. Truncating the Middle with a Rigorous Error Bound

Choose a threshold $0 \lt \rho \lt 1$, and let $L$ be the largest index on the left side with $r_L \gt \rho$. By symmetry, the same number of disks is retained on the right side. The approximation is then

$$E(N,M)\approx \frac{N}{2}+\sum_{i=1}^{L} r_i^M,$$

because the factor $1/2$ in the exact formula is canceled by the two symmetric edge blocks.

The omitted middle block contains $N-2L$ disks. In the large-$N$ regime where the accelerated method is used, those terms satisfy $\lvert r_i\rvert \le \rho$, so the neglected contribution is bounded by

$$\left|E(N,M)-\left(\frac{N}{2}+\sum_{i=1}^{L} r_i^M\right)\right|\le \frac{1}{2}(N-2L)\rho^M.$$

The implementations start from a target tail size $\varepsilon=10^{-4}$ and choose

$$\rho=\left(\frac{2\varepsilon}{N}\right)^{1/M},$$

so that the theoretical omitted contribution is tiny. They then compute the explicit upper bound above and compare it with $0.005$, the threshold relevant for safe rounding to two decimal places. If needed, the threshold is tightened slightly and the edge sum is recomputed.

6. Small Worked Example

The checkpoint $E(3,1)=10/9$ follows immediately from the formula. For $N=3$, the edge disks have

$$q_1=q_3=\frac{4}{9},\qquad r_1=r_3=-\frac{1}{9},$$

while the middle disk has

$$q_2=\frac{2}{9},\qquad r_2=-\frac{5}{9}.$$

Therefore

$$E(3,1)=\frac{3}{2}+\frac{1}{2}\left(-\frac{1}{9}-\frac{5}{9}-\frac{1}{9}\right)=\frac{10}{9}.$$

Squaring the same coefficients gives

$$E(3,2)=\frac{3}{2}+\frac{1}{2}\left(\frac{1}{81}+\frac{25}{81}+\frac{1}{81}\right)=\frac{5}{3}.$$

These are exactly the small exact checkpoints used by the implementations before the large-input calculation.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. For small inputs they evaluate the exact formula

$$E(N,M)=\frac{N}{2}+\frac{1}{2}\sum_{i=1}^{N} r_i^M$$

directly, which is useful for validation and checkpoint testing. For very large inputs they switch to the accelerated method: compute the coefficient for a single disk, find the edge cutoff by binary search, sum only the retained edge terms, and use symmetry to reconstruct the total contribution.

The expensive part is the edge summation, so the implementation partitions that retained range across multiple workers and combines the partial sums at the end. The arithmetic itself remains simple: each retained disk contributes one power $r_i^M$, and the final estimate is protected by the tail bound described above.

Complexity Analysis

The exact method is $O(N)$ time and $O(1)$ memory, which is fine for checkpoints but impossible for $N=10^{10}$. The accelerated method uses $O(\log N)$ work to find the cutoff and $O(L)$ work to sum the retained edge terms, where $L \ll N$. Memory usage stays $O(1)$ apart from a small number of partial accumulators used for parallel execution.

Footnotes and References

Problem page: https://projecteuler.net/problem=430
Expected value and linearity: Wikipedia - Expected value
Two-state Markov chains: Wikipedia - Markov chain
Binomial parity identity and repeated Bernoulli trials: Wikipedia - Binomial theorem

Problem 430 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    u64 n = 10000000000ULL;
    int m = 4000;
    int threads = 0;
    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;
    }
    try {
        value = std::stoi(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = static_cast<u64>(std::stoull(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_u64_after_prefix(arg, "--n=", options.n) ||
            parse_int_after_prefix(arg, "--m=", options.m) ||
            parse_int_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1ULL && options.m >= 0 && options.threads >= 0;
}

int choose_thread_count(int requested, u64 work_items) {
    if (work_items <= 1ULL) {
        return 1;
    }
    int threads = requested;
    if (threads <= 0) {
        threads = static_cast<int>(std::thread::hardware_concurrency());
    }
    if (threads <= 0) {
        threads = 4;
    }
    if (threads > static_cast<int>(work_items)) {
        threads = static_cast<int>(work_items);
    }
    if (threads < 1) {
        threads = 1;
    }
    return threads;
}

long double r_value(u64 n, u64 i) {
    // r_i = 2*P(not flipped in one turn) - 1.
    const long double nf = static_cast<long double>(n);
    const long double n2 = nf * nf;
    const long double a = static_cast<long double>(i - 1ULL);
    const long double b = static_cast<long double>(n - i);
    return (2.0L * (a * a + b * b) - n2) / n2;
}

long double expected_exact(u64 n, int m) {
    long double ans = 0.0L;
    for (u64 i = 1; i <= n; ++i) {
        const long double r = r_value(n, i);
        ans += 0.5L * (1.0L + std::powl(r, static_cast<long double>(m)));
    }
    return ans;
}

u64 edge_limit_by_rho(u64 n, long double rho) {
    u64 lo = 0;
    u64 hi = n / 2ULL;
    while (lo < hi) {
        const u64 mid = lo + (hi - lo + 1ULL) / 2ULL;
        if (r_value(n, mid) > rho) {
            lo = mid;
        } else {
            hi = mid - 1ULL;
        }
    }
    return lo;
}

long double edge_sum_parallel(u64 n, int m, u64 edge_count, int requested_threads) {
    if (edge_count == 0ULL) {
        return 0.0L;
    }
    const int thread_count = choose_thread_count(requested_threads, edge_count);
    std::vector<long double> partial(static_cast<std::size_t>(thread_count), 0.0L);
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(thread_count));

    for (int t = 0; t < thread_count; ++t) {
        const u64 begin = edge_count * static_cast<u64>(t) / static_cast<u64>(thread_count) + 1ULL;
        const u64 end =
            edge_count * static_cast<u64>(t + 1) / static_cast<u64>(thread_count);
        workers.emplace_back([&, begin, end, t]() {
            long double local = 0.0L;
            for (u64 i = begin; i <= end; ++i) {
                const long double r = r_value(n, i);
                local += std::powl(r, static_cast<long double>(m));
            }
            partial[static_cast<std::size_t>(t)] = local;
        });
    }
    for (auto& w : workers) {
        w.join();
    }

    long double sum = 0.0L;
    for (long double x : partial) {
        sum += x;
    }
    return sum;
}

long double expected_fast(u64 n, int m, int threads) {
    if (m == 0) {
        return static_cast<long double>(n);
    }
    if (n <= 2000000ULL) {
        return expected_exact(n, m);
    }

    // Tail bound target: omitted contribution to E is < 1e-4.
    const long double eps = 1e-4L;
    const long double rho = std::expl(std::log((2.0L * eps) / static_cast<long double>(n)) /
                                      static_cast<long double>(m));

    const u64 l = edge_limit_by_rho(n, rho);
    const long double edge = edge_sum_parallel(n, m, l, threads);

    // E = N/2 + 1/2 * sum_i r_i^M.
    // Approximation keeps the 2*L edge terms exactly (by symmetry) and drops middle terms.
    long double estimate = static_cast<long double>(n) * 0.5L + edge;

    // Safety check for rounding to 2 decimals.
    const u64 omitted = n - 2ULL * l;
    const long double err_bound =
        0.5L * static_cast<long double>(omitted) * std::powl(rho, static_cast<long double>(m));
    if (err_bound >= 0.005L) {
        // Fallback: tighter threshold if needed.
        const long double rho2 = rho * 0.99L;
        const u64 l2 = edge_limit_by_rho(n, rho2);
        const long double edge2 = edge_sum_parallel(n, m, l2, threads);
        estimate = static_cast<long double>(n) * 0.5L + edge2;
    }

    return estimate;
}

bool close_to(long double a, long double b, long double tol) {
    return std::fabsl(a - b) <= tol;
}

bool run_checkpoints() {
    if (!close_to(expected_exact(3, 1), 10.0L / 9.0L, 1e-15L)) {
        std::cerr << "Checkpoint failed: E(3,1)\n";
        return false;
    }
    if (!close_to(expected_exact(3, 2), 5.0L / 3.0L, 1e-15L)) {
        std::cerr << "Checkpoint failed: E(3,2)\n";
        return false;
    }
    if (!close_to(expected_exact(10, 4), 5.15702608L, 1e-12L)) {
        std::cerr << "Checkpoint failed: E(10,4)\n";
        return false;
    }
    if (!close_to(expected_exact(100, 10), 51.8928010382861L, 1e-10L)) {
        std::cerr << "Checkpoint failed: E(100,10)\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;
    }

    const long double answer = expected_fast(options.n, options.m, options.threads);
    std::cout << std::fixed << std::setprecision(2) << static_cast<double>(answer) << '\n';
    return 0;
}

Python

import sys
import math
from multiprocessing import Pool, cpu_count

def r_value(n, i):
    nf = float(n)
    n2 = nf * nf
    a = float(i - 1)
    b = float(n - i)
    return (2.0 * (a * a + b * b) - n2) / n2

def expected_exact(n, m):
    ans = 0.0
    for i in range(1, n + 1):
        r = r_value(n, i)
        ans += 0.5 * (1.0 + math.pow(r, m))
    return ans

def edge_limit_by_rho(n, rho):
    lo = 0
    hi = n // 2
    while lo < hi:
        mid = lo + (hi - lo + 1) // 2
        if r_value(n, mid) > rho:
            lo = mid
        else:
            hi = mid - 1
    return lo

def edge_sum_chunk(n, m, begin, end):
    local = 0.0
    for i in range(begin, end + 1):
        r = r_value(n, i)
        local += math.pow(r, m)
    return local

def edge_sum_parallel(n, m, edge_count, num_threads):
    if edge_count == 0:
        return 0.0
    
    threads = min(num_threads, edge_count)
    if threads < 1: threads = 1
    
    tasks = []
    for t in range(threads):
        begin = edge_count * t // threads + 1
        end = edge_count * (t + 1) // threads
        tasks.append((n, m, begin, end))
        
    with Pool(threads) as pool:
        results = pool.starmap(edge_sum_chunk, tasks)
        
    return sum(results)

def expected_fast(n, m, threads=0):
    if m == 0:
        return float(n)
    if n <= 2000000:
        return expected_exact(n, m)
        
    if threads <= 0:
        threads = cpu_count()
        if threads <= 0: threads = 4
        
    eps = 1e-4
    rho = math.exp(math.log((2.0 * eps) / float(n)) / float(m))
    
    l = edge_limit_by_rho(n, rho)
    edge = edge_sum_parallel(n, m, l, threads)
    
    estimate = float(n) * 0.5 + edge
    
    omitted = n - 2 * l
    err_bound = 0.5 * float(omitted) * math.pow(rho, m)
    if err_bound >= 0.005:
        rho2 = rho * 0.99
        l2 = edge_limit_by_rho(n, rho2)
        edge2 = edge_sum_parallel(n, m, l2, threads)
        estimate = float(n) * 0.5 + edge2
        
    return estimate

def solve():
    n = 10000000000
    m = 4000
    ans = expected_fast(n, m)
    return f"{ans:.2f}"

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

Java

import java.util.stream.IntStream;
import java.util.stream.LongStream;

public class Euler430 {

    static double rValue(long n, long i) {
        double nf = (double) n;
        double n2 = nf * nf;
        double a = (double) (i - 1);
        double b = (double) (n - i);
        return (2.0 * (a * a + b * b) - n2) / n2;
    }

    static double expectedExact(long n, int m) {
        double ans = 0.0;
        for (long i = 1; i <= n; i++) {
            double r = rValue(n, i);
            ans += 0.5 * (1.0 + Math.pow(r, m));
        }
        return ans;
    }

    static long edgeLimitByRho(long n, double rho) {
        long lo = 0;
        long hi = n / 2;
        while (lo < hi) {
            long mid = lo + (hi - lo + 1) / 2;
            if (rValue(n, mid) > rho) {
                lo = mid;
            } else {
                hi = mid - 1;
            }
        }
        return lo;
    }

    static double edgeSumParallel(long n, int m, long edgeCount) {
        if (edgeCount == 0)
            return 0.0;
        int threads = Math.min(Runtime.getRuntime().availableProcessors(), (int) Math.min(edgeCount, 1000));
        if (threads < 1)
            threads = 1;

        final int numThreads = threads;
        return IntStream.range(0, numThreads)
                .parallel()
                .mapToDouble(t -> {
                    long begin = edgeCount * t / numThreads + 1;
                    long end = edgeCount * (t + 1) / numThreads;
                    double local = 0.0;
                    for (long i = begin; i <= end; i++) {
                        double r = rValue(n, i);
                        local += Math.pow(r, m);
                    }
                    return local;
                })
                .sum();
    }

    static double expectedFast(long n, int m) {
        if (m == 0)
            return (double) n;
        if (n <= 2000000)
            return expectedExact(n, m);

        double eps = 1e-4;
        double rho = Math.exp(Math.log((2.0 * eps) / (double) n) / (double) m);

        long l = edgeLimitByRho(n, rho);
        double edge = edgeSumParallel(n, m, l);

        double estimate = (double) n * 0.5 + edge;

        long omitted = n - 2 * l;
        double errBound = 0.5 * (double) omitted * Math.pow(rho, m);
        if (errBound >= 0.005) {
            double rho2 = rho * 0.99;
            long l2 = edgeLimitByRho(n, rho2);
            double edge2 = edgeSumParallel(n, m, l2);
            estimate = (double) n * 0.5 + edge2;
        }

        return estimate;
    }

    public static String solve() {
        long n = 10000000000L;
        int m = 4000;
        double ans = expectedFast(n, m);
        return String.format(java.util.Locale.US, "%.2f", ans);
    }

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