Problem 645: Every Day Is a Holiday

Project Euler Problem 645 Solution

EulerSolve provides an optimized solution for Project Euler Problem 645, Every Day Is a Holiday, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Place the $D$ days on a cycle. One random draw chooses a day uniformly, and that draw covers two consecutive days: the chosen day and its cyclic successor. Let $T_D$ be the number of draws needed until every day has been covered at least once. The task is to compute the expectation $$E(D)=\mathbb{E}[T_D]$$ for the large target value $D=10000$. The C++, Python, and Java implementations do not simulate the process. Instead, they turn the expectation into a one-dimensional integral whose integrand comes from inclusion-exclusion on uncovered days and a $2\times2$ transfer matrix on the cycle. Mathematical Approach Number the days modulo $D$, and let a draw on day $i$ cover the pair $(i,i+1)$. The stopping time is $T_D$, the first moment when all days are covered. Step 1: Rewrite the expectation as a tail sum For any nonnegative integer-valued stopping time, $$E(D)=\sum_{n=0}^{\infty}\Pr(T_D>n).$$ So the whole problem becomes: after $n$ draws, what is the probability that at least one day is still uncovered? A day $j$ is uncovered exactly when no draw landed on $j$ and no draw landed on $j-1$, because those are the only two starting points that could cover $j$....

Detailed mathematical approach

Problem Summary

Place the $D$ days on a cycle. One random draw chooses a day uniformly, and that draw covers two consecutive days: the chosen day and its cyclic successor. Let $T_D$ be the number of draws needed until every day has been covered at least once. The task is to compute the expectation

$$E(D)=\mathbb{E}[T_D]$$

for the large target value $D=10000$.

The C++, Python, and Java implementations do not simulate the process. Instead, they turn the expectation into a one-dimensional integral whose integrand comes from inclusion-exclusion on uncovered days and a $2\times2$ transfer matrix on the cycle.

Mathematical Approach

Number the days modulo $D$, and let a draw on day $i$ cover the pair $(i,i+1)$. The stopping time is $T_D$, the first moment when all days are covered.

Step 1: Rewrite the expectation as a tail sum

For any nonnegative integer-valued stopping time,

$$E(D)=\sum_{n=0}^{\infty}\Pr(T_D>n).$$

So the whole problem becomes: after $n$ draws, what is the probability that at least one day is still uncovered?

A day $j$ is uncovered exactly when no draw landed on $j$ and no draw landed on $j-1$, because those are the only two starting points that could cover $j$.

Step 2: Apply inclusion-exclusion to uncovered sets

For a subset $S$ of days, define

$$N(S)=S\cup(S-1), \qquad S-1=\{i-1 \pmod D : i\in S\}.$$

If every day in $S$ is still uncovered after $n$ draws, then every draw must avoid the set $N(S)$. Since each draw is uniform over the $D$ starting days,

$$\Pr(A_S(n))=\left(1-\frac{|N(S)|}{D}\right)^n,$$

where $A_S(n)$ is the event that all days in $S$ remain uncovered after $n$ draws.

Now inclusion-exclusion over the uncovered days gives

$$\Pr(T_D>n)=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\left(1-\frac{|N(S)|}{D}\right)^n,$$

where $C_D$ denotes the cycle of $D$ days.

Step 3: Sum the geometric series and convert to an integral

Insert the previous formula into the tail sum and interchange the order of summation:

$$\frac{E(D)}{D}=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\frac{1}{|N(S)|}.$$

The reciprocal is converted with

$$\frac{1}{m}=\int_0^1 x^{m-1}\,dx.$$

Therefore

$$\frac{E(D)}{D}=\int_0^1\frac{1-Z_D(x)}{x}\,dx,$$

where

$$Z_D(x)=\sum_{S\subseteq C_D}(-1)^{|S|}x^{|N(S)|}.$$

The empty set contributes the leading $1$, which is why the numerator becomes $1-Z_D(x)$.

Step 4: Encode the subsets with a transfer matrix

Represent a subset $S$ by a cyclic binary word $s_1,\dots,s_D$, where $s_i=1$ means day $i$ is required to remain uncovered. Then

$$|N(S)|=\sum_{i=1}^{D}s_i+\sum_{i=1}^{D}(1-s_i)s_{i+1}, \qquad s_{D+1}=s_1.$$

The first sum counts the chosen days themselves, and the second sum counts one extra predecessor for each run of consecutive $1$'s around the cycle.

This makes the weight factorize locally:

$$(-1)^{|S|}x^{|N(S)|}=\prod_{i=1}^{D} M_{s_i,s_{i+1}}(x),$$

with transfer matrix

$$M(x)=\begin{pmatrix}1 & -x^2 \\ 1 & -x\end{pmatrix}.$$

Summing over all cyclic binary words is exactly the trace of the $D$-th matrix power, so

$$Z_D(x)=\operatorname{tr}(M(x)^D).$$

Step 5: Diagonalize the matrix

The characteristic polynomial of $M(x)$ is

$$\lambda^2-(1-x)\lambda-x(1-x)=0.$$

Its two eigenvalues are

$$\lambda_{1,2}(x)=\frac{1-x\pm\sqrt{(1-x)(1+3x)}}{2}.$$

Hence

$$Z_D(x)=\lambda_1(x)^D+\lambda_2(x)^D,$$

and the expectation becomes

$$\boxed{\frac{E(D)}{D}=\int_0^1\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}\,dx.}$$

This is the exact formula used by the implementation.

Step 6: Endpoint behavior and numerical shape

The integrand looks singular because of the division by $x$, but the singularity at $x=0$ is removable. Expanding the eigenvalues near zero gives

$$\lambda_1(x)=1-x^2+O(x^3), \qquad \lambda_2(x)=-x+O(x^2),$$

$$1-\lambda_1(x)^D-\lambda_2(x)^D = D x^2 + O(x^3),$$

and therefore the integrand is $Dx+O(x^2)$, which tends to $0$. At $x=1$, both eigenvalues are $0$, so the integrand tends to $1$.

Worked Example: $D=5$

For five days, the transfer-matrix trace simplifies to

$$Z_5(x)=1-5x^2+5x^3-x^5.$$

Thus

$$\frac{E(5)}{5}=\int_0^1\frac{5x^2-5x^3+x^5}{x}\,dx=\int_0^1(5x-5x^2+x^4)\,dx.$$

Integrating term by term gives

$$\frac{E(5)}{5}=\frac{5}{2}-\frac{5}{3}+\frac{1}{5}=\frac{31}{30},$$

$$E(5)=\frac{31}{6}.$$

This matches the exact checkpoint used by the implementations. For the smallest nontrivial cycle, $Z_2(x)=1-x^2$, giving $E(2)=1$.

How the Code Works

The C++, Python, and Java implementations all evaluate the same integrand

$$I_D(x)=\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}.$$

They first compute the discriminant $\sqrt{(1-x)(1+3x)}$, then the two eigenvalues, and then the numerator. The endpoints are handled explicitly: the limiting value at $x=0$ is $0$, and the value at $x=1$ is $1$.

The delicate part is the term $1-\lambda_1(x)^D$. When $\lambda_1(x)$ is very close to $1$, direct subtraction loses significant digits. To avoid that, the implementation rewrites it as

$$1-\lambda_1(x)^D=-\operatorname{expm1}\!\bigl(D\log \lambda_1(x)\bigr)$$

whenever $D\log \lambda_1(x)$ is close to $0$. The second eigenvalue is smaller in magnitude and is raised directly to the $D$-th power.

The integral on $[0,1]$ is then evaluated with Simpson's rule using an even number of intervals. Starting from an initial mesh, the code repeatedly doubles the interval count and compares two consecutive Simpson estimates:

$$|S_{2N}-S_N|\le \text{rtol}\cdot \max(1,|S_{2N}|).$$

Once this relative criterion is satisfied, the refined estimate is returned. The C++ implementation also supports parallel accumulation of the interior Simpson nodes when the mesh is large, and it verifies the derivation against the known values $E(2)=1$, $E(5)=31/6$, and $E(365)\approx 1174.3501$.

Complexity Analysis

If the final Simpson mesh uses $N$ subintervals, one pass costs $O(N)$ evaluations of the integrand. Because the mesh doubles geometrically, the total work up to the accepted mesh is still $O(N)$. The single-threaded implementations use $O(1)$ auxiliary memory, while the parallel C++ version uses $O(T)$ extra memory for $T$ partial sums.

Footnotes and References

Problem page: https://projecteuler.net/problem=645
Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
Transfer-matrix method: Wikipedia - Transfer-matrix method
Eigenvalues and eigenvectors: Wikipedia - Eigenvalues and eigenvectors
Simpson's rule: Wikipedia - Simpson's rule

Problem 645 source code

C++

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

namespace {

struct Options {
    int days = 10'000;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0;
    int initial_intervals = 4'096;
    int max_intervals = 1 << 20;
    long double relative_tolerance = 1e-12L;
};

bool parse_int_after_prefix(const std::string& arg, const char* prefix, int& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    std::int64_t parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const int digit = c - '0';
        if (parsed > (std::numeric_limits<std::int64_t>::max() - digit) / 10) {
            return false;
        }
        parsed = parsed * 10 + digit;
    }

    if (parsed > std::numeric_limits<int>::max()) {
        return false;
    }
    value = static_cast<int>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    int parsed = 0;
    if (!parse_int_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed < 0) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_long_double_after_prefix(const std::string& arg,
                                    const char* prefix,
                                    long double& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    char* end_ptr = nullptr;
    const long double parsed = std::strtold(tail.c_str(), &end_ptr);
    if (end_ptr == nullptr || *end_ptr != '\0' || !std::isfinite(parsed)) {
        return false;
    }
    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 == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }

        int parsed_int = 0;
        if (parse_int_after_prefix(arg, "--days=", parsed_int) ||
            parse_int_after_prefix(arg, "--d=", parsed_int)) {
            options.days = parsed_int;
            continue;
        }
        if (parse_int_after_prefix(arg, "--initial-intervals=", parsed_int)) {
            options.initial_intervals = parsed_int;
            continue;
        }
        if (parse_int_after_prefix(arg, "--max-intervals=", parsed_int)) {
            options.max_intervals = parsed_int;
            continue;
        }

        unsigned parsed_unsigned = 0;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        long double parsed_ld = 0.0L;
        if (parse_long_double_after_prefix(arg, "--rtol=", parsed_ld)) {
            options.relative_tolerance = parsed_ld;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.days < 2) {
        std::cerr << "--days must be >= 2.\n";
        return false;
    }
    if (options.initial_intervals < 4 || (options.initial_intervals % 2) != 0) {
        std::cerr << "--initial-intervals must be an even integer >= 4.\n";
        return false;
    }
    if (options.max_intervals < options.initial_intervals || (options.max_intervals % 2) != 0) {
        std::cerr << "--max-intervals must be an even integer >= --initial-intervals.\n";
        return false;
    }
    if (!(options.relative_tolerance > 0.0L) || !std::isfinite(options.relative_tolerance)) {
        std::cerr << "--rtol must be a positive finite number.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(bool allow_multithreading,
                             unsigned requested_threads,
                             std::size_t workload) {
    if (!allow_multithreading || workload < 20'000ULL) {
        return 1;
    }

    unsigned threads = requested_threads;
    if (threads == 0) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0) {
            threads = 1;
        }
    }
    return std::max(1u, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}

long double pow_int_long_double(long double base, int exponent) {
    if (exponent == 0) {
        return 1.0L;
    }

    int e = exponent;
    long double b = base;
    long double result = 1.0L;
    while (e > 0) {
        if ((e & 1) != 0) {
            result *= b;
        }
        b *= b;
        e >>= 1;
    }
    return result;
}

long double integrand_value(int days, long double x) {
    if (x <= 0.0L) {
        return 0.0L;
    }
    if (x >= 1.0L) {
        return 1.0L;
    }

    // Inclusion-exclusion + transfer matrix on the cycle gives:
    // E(D) / D = integral_0^1 (1 - Z_D(x)) / x dx,
    // Z_D(x) = lambda_1(x)^D + lambda_2(x)^D.
    const long double disc = std::sqrt((1.0L - x) * (1.0L + 3.0L * x));
    const long double lambda1 = (1.0L - x + disc) * 0.5L;
    const long double lambda2 = (1.0L - x - disc) * 0.5L;

    long double one_minus_lambda1_pow = 1.0L;
    if (lambda1 > 0.0L) {
        const long double log_term = static_cast<long double>(days) * std::log(lambda1);
        if (log_term > -1e-4L) {
            one_minus_lambda1_pow = -std::expm1(log_term);
        } else {
            one_minus_lambda1_pow = 1.0L - std::exp(log_term);
        }
    }

    const long double lambda2_pow = pow_int_long_double(lambda2, days);
    const long double numerator = one_minus_lambda1_pow - lambda2_pow;
    return numerator / x;
}

long double simpson_integral(int days,
                             int intervals,
                             bool allow_multithreading,
                             unsigned requested_threads) {
    const long double h = 1.0L / static_cast<long double>(intervals);
    const std::size_t workload = static_cast<std::size_t>(intervals - 1);
    const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, workload);

    if (threads == 1) {
        long double sum = integrand_value(days, 0.0L) + integrand_value(days, 1.0L);
        for (int i = 1; i < intervals; ++i) {
            const long double x = static_cast<long double>(i) * h;
            const int weight = (i & 1) != 0 ? 4 : 2;
            sum += static_cast<long double>(weight) * integrand_value(days, x);
        }
        return static_cast<long double>(days) * h * sum / 3.0L;
    }

    std::vector<long double> partial(threads, 0.0L);
    std::vector<std::thread> workers;
    workers.reserve(threads);

    const int block = (intervals - 1 + static_cast<int>(threads) - 1) / static_cast<int>(threads);
    for (unsigned t = 0; t < threads; ++t) {
        const int begin = 1 + static_cast<int>(t) * block;
        const int end = std::min(intervals, begin + block);
        workers.emplace_back([&, t, begin, end]() {
            long double local = 0.0L;
            for (int i = begin; i < end; ++i) {
                const long double x = static_cast<long double>(i) * h;
                const int weight = (i & 1) != 0 ? 4 : 2;
                local += static_cast<long double>(weight) * integrand_value(days, x);
            }
            partial[t] = local;
        });
    }

    for (std::thread& worker : workers) {
        worker.join();
    }

    long double sum = integrand_value(days, 0.0L) + integrand_value(days, 1.0L);
    for (const long double value : partial) {
        sum += value;
    }
    return static_cast<long double>(days) * h * sum / 3.0L;
}

long double compute_expectation(int days,
                                bool allow_multithreading,
                                unsigned requested_threads,
                                int initial_intervals,
                                int max_intervals,
                                long double relative_tolerance) {
    int intervals = initial_intervals;
    long double previous =
        simpson_integral(days, intervals, allow_multithreading, requested_threads);

    while (intervals < max_intervals) {
        intervals *= 2;
        const long double current =
            simpson_integral(days, intervals, allow_multithreading, requested_threads);
        const long double delta = std::fabs(current - previous);
        const long double scale = std::max(1.0L, std::fabs(current));
        if (delta <= relative_tolerance * scale) {
            return current;
        }
        previous = current;
    }

    return previous;
}

bool run_checkpoints(const Options& options) {
    const long double e2 = compute_expectation(
        2, options.allow_multithreading, options.requested_threads, 512, 32'768, 1e-13L);
    const long double e5 = compute_expectation(
        5, options.allow_multithreading, options.requested_threads, 512, 32'768, 1e-13L);
    const long double e365 = compute_expectation(
        365, options.allow_multithreading, options.requested_threads, 1'024, 65'536, 1e-12L);

    const long double expected_e2 = 1.0L;
    const long double expected_e5 = 31.0L / 6.0L;
    const long double expected_e365 = 1174.3501L;

    const bool ok2 = std::fabs(e2 - expected_e2) < 1e-11L;
    const bool ok5 = std::fabs(e5 - expected_e5) < 1e-11L;
    const bool ok365 = std::fabs(e365 - expected_e365) < 5e-5L;

    if (!ok2 || !ok5 || !ok365) {
        std::cerr << std::setprecision(18);
        std::cerr << "Checkpoint failed.\n";
        std::cerr << "E(2):   got " << e2 << ", expected " << expected_e2 << '\n';
        std::cerr << "E(5):   got " << e5 << ", expected " << expected_e5 << '\n';
        std::cerr << "E(365): got " << e365 << ", expected approx " << expected_e365 << '\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(options)) {
        return 1;
    }

    const long double answer = compute_expectation(options.days,
                                                   options.allow_multithreading,
                                                   options.requested_threads,
                                                   options.initial_intervals,
                                                   options.max_intervals,
                                                   options.relative_tolerance);

    std::cout << std::fixed << std::setprecision(4) << answer << '\n';
    return 0;
}

Python

import math

def integrand_value(days, x):
    if x <= 0.0: return 0.0
    if x >= 1.0: return 1.0
    
    disc = math.sqrt((1.0 - x) * (1.0 + 3.0 * x))
    lambda1 = (1.0 - x + disc) * 0.5
    lambda2 = (1.0 - x - disc) * 0.5
    
    one_minus_lambda1_pow = 1.0
    if lambda1 > 0.0:
        log_term = days * math.log(lambda1)
        if log_term > -1e-4:
            one_minus_lambda1_pow = -math.expm1(log_term)
        else:
            one_minus_lambda1_pow = 1.0 - math.exp(log_term)
            
    lambda2_pow = math.pow(lambda2, days) if isinstance(days, int) else lambda2 ** days
    numerator = one_minus_lambda1_pow - lambda2_pow
    return numerator / x

def simpson_integral(days, intervals):
    h = 1.0 / intervals
    
    sum_val = integrand_value(days, 0.0) + integrand_value(days, 1.0)
    for i in range(1, intervals):
        x = i * h
        weight = 4 if (i % 2) != 0 else 2
        sum_val += weight * integrand_value(days, x)
        
    return days * h * sum_val / 3.0

def compute_expectation(days, initial_intervals, max_intervals, relative_tolerance):
    intervals = initial_intervals
    previous = simpson_integral(days, intervals)
    
    while intervals < max_intervals:
        intervals *= 2
        current = simpson_integral(days, intervals)
        delta = abs(current - previous)
        scale = max(1.0, abs(current))
        if delta <= relative_tolerance * scale:
            return current
        previous = current
        
    return previous

def solve():
    ans = compute_expectation(10000, 4096, 1048576, 1e-12)
    return "{:.4f}".format(ans)

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

Java

import java.util.Locale;

public class Euler645 {

    static double integrandValue(int days, double x) {
        if (x <= 0.0)
            return 0.0;
        if (x >= 1.0)
            return 1.0;

        double disc = Math.sqrt((1.0 - x) * (1.0 + 3.0 * x));
        double lambda1 = (1.0 - x + disc) * 0.5;
        double lambda2 = (1.0 - x - disc) * 0.5;

        double oneMinusLambda1Pow = 1.0;
        if (lambda1 > 0.0) {
            double logTerm = days * Math.log(lambda1);
            if (logTerm > -1e-4) {
                oneMinusLambda1Pow = -Math.expm1(logTerm);
            } else {
                oneMinusLambda1Pow = 1.0 - Math.exp(logTerm);
            }
        }

        double lambda2Pow = Math.pow(lambda2, days);
        double numerator = oneMinusLambda1Pow - lambda2Pow;
        return numerator / x;
    }

    static double simpsonIntegral(int days, int intervals) {
        double h = 1.0 / intervals;
        double sum = integrandValue(days, 0.0) + integrandValue(days, 1.0);
        for (int i = 1; i < intervals; ++i) {
            double x = i * h;
            int weight = (i % 2 != 0) ? 4 : 2;
            sum += weight * integrandValue(days, x);
        }
        return days * h * sum / 3.0;
    }

    static double computeExpectation(int days, int initialIntervals, int maxIntervals, double relativeTol) {
        int intervals = initialIntervals;
        double previous = simpsonIntegral(days, intervals);

        while (intervals < maxIntervals) {
            intervals *= 2;
            double current = simpsonIntegral(days, intervals);
            double delta = Math.abs(current - previous);
            double scale = Math.max(1.0, Math.abs(current));
            if (delta <= relativeTol * scale) {
                return current;
            }
            previous = current;
        }
        return previous;
    }

    public static String solve() {
        double ans = computeExpectation(10000, 4096, 1048576, 1e-12);
        return String.format(Locale.US, "%.4f", ans);
    }

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