Problem 441: The Inverse Summation of Coprime Couples

Project Euler Problem 441 Solution

EulerSolve provides an optimized solution for Project Euler Problem 441, The Inverse Summation of Coprime Couples, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer $M \ge 2$, define $$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$$ The checkpoints are $S(2)=\tfrac12$, $S(10)=6.9146825397\ldots$, and $S(100)=58.2962380621\ldots$. The task is to evaluate $S(10^7)$ to four decimal places. A direct evaluation over all triples $(M,p,q)$ is far too slow, so the implementation reorganizes the sum and replaces the coprimality test by Möbius inclusion-exclusion. Mathematical Approach Step 1: Exchange the Order of Summation Fix a coprime pair $1 \le p \lt q \le N$. Such a pair contributes to $R(M)$ exactly when $q \le M$ and $M \le p+q$. Since $M$ is also bounded by $N$, the admissible values are $$q \le M \le \min(N,p+q).$$ Therefore this pair appears with multiplicity $$w_N(p,q)=\min(N,p+q)-q+1.$$ After switching the order of summation, $$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}.$$ This is the key simplification: the problem becomes a sum of independent contributions indexed by the larger denominator $q$....

Detailed mathematical approach

Problem Summary

For each integer $M \ge 2$, define

$$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$$

The checkpoints are $S(2)=\tfrac12$, $S(10)=6.9146825397\ldots$, and $S(100)=58.2962380621\ldots$. The task is to evaluate $S(10^7)$ to four decimal places. A direct evaluation over all triples $(M,p,q)$ is far too slow, so the implementation reorganizes the sum and replaces the coprimality test by Möbius inclusion-exclusion.

Mathematical Approach

Step 1: Exchange the Order of Summation

Fix a coprime pair $1 \le p \lt q \le N$. Such a pair contributes to $R(M)$ exactly when $q \le M$ and $M \le p+q$. Since $M$ is also bounded by $N$, the admissible values are

$$q \le M \le \min(N,p+q).$$

Therefore this pair appears with multiplicity

$$w_N(p,q)=\min(N,p+q)-q+1.$$

After switching the order of summation,

$$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}.$$

This is the key simplification: the problem becomes a sum of independent contributions indexed by the larger denominator $q$.

Step 2: Encode Coprimality with Möbius Inversion

Define the coprime reciprocal sum and the coprime counting function

$$A_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}\frac{1}{p}, \qquad C_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}1.$$

Using the identity

$$\mathbf{1}_{\gcd(p,q)=1}=\sum_{d \mid \gcd(p,q)}\mu(d),$$

we can rewrite both quantities in terms of divisors of $q$. If $H_t=\sum_{j=1}^{t}\frac{1}{j}$ denotes the $t$-th harmonic number, then

$$A_q(m)=\sum_{d \mid q}\frac{\mu(d)}{d}H_{\lfloor m/d \rfloor}, \qquad C_q(m)=\sum_{d \mid q}\mu(d)\left\lfloor \frac{m}{d}\right\rfloor.$$

Only squarefree divisors matter, because $\mu(d)=0$ whenever $d$ contains a repeated prime factor. That is why the implementation factors $q$ into distinct primes and enumerates all squarefree divisor subsets with their Möbius signs.

Step 3: The Easy Half $q \le \lfloor N/2 \rfloor$

If $q \le \lfloor N/2 \rfloor$, then for every $p \lt q$ we have

$$p+q \le (q-1)+q = 2q-1 \le N-1,$$

so the upper bound $\min(N,p+q)$ is just $p+q$. Hence

$$w_N(p,q)=p+1.$$

The contribution of this $q$ becomes

$$\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{p+1}{p} =\frac{1}{q}\left(\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}1+\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$

The first sum is $\varphi(q)$, because among $1,2,\dots,q-1$ there are exactly $\varphi(q)$ integers coprime to $q$. Therefore

$$T_N(q)=\frac{\varphi(q)+A_q(q-1)}{q}\qquad \text{for } q \le \left\lfloor \frac{N}{2}\right\rfloor.$$

Step 4: The Truncated Half $q > \lfloor N/2 \rfloor$

Now let $q > \lfloor N/2 \rfloor$ and set $a=N-q$. Then $0 \le a \lt q$, and the weight changes at $p=a$:

$$w_N(p,q)= \begin{cases} p+1, & 1 \le p \le a,\\ a+1, & a \lt p \lt q. \end{cases}$$

So the contribution splits into two pieces:

$$T_N(q)=\frac{1}{q}\left(\sum_{\substack{1\le p \le a\\ \gcd(p,q)=1}}\left(1+\frac{1}{p}\right) +(a+1)\sum_{\substack{a \lt p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$

Writing the second reciprocal sum as $A_q(q-1)-A_q(a)$, we obtain

$$T_N(q)=\frac{C_q(a)+A_q(a)+(a+1)\bigl(A_q(q-1)-A_q(a)\bigr)}{q}.$$

This is the exact branch used by the implementation once $q$ passes the midpoint.

Step 5: Final Closed Form

Combining the two ranges gives the complete formula

$$\boxed{ S(N)= \sum_{q=2}^{\lfloor N/2 \rfloor}\frac{\varphi(q)+A_q(q-1)}{q} + \sum_{q=\lfloor N/2 \rfloor+1}^{N} \frac{C_q(N-q)+A_q(N-q)+(N-q+1)\bigl(A_q(q-1)-A_q(N-q)\bigr)}{q} .}$$

All expensive work has now been reduced to three precomputable ingredients: Euler's totient values, harmonic prefixes, and factorizations of the integers $2,\dots,N$.

Worked Example: $N=10$

For $q=5$, we are still in the first branch. Since $\varphi(5)=4$ and

$$A_5(4)=1+\frac12+\frac13+\frac14=\frac{25}{12},$$

we get

$$T_{10}(5)=\frac{4+\frac{25}{12}}{5}=\frac{73}{60}=1.216666\ldots$$

For $q=8$, we are in the second branch with $a=10-8=2$. Here

$$C_8(2)=1,\qquad A_8(2)=1,\qquad A_8(7)-A_8(2)=\frac13+\frac15+\frac17=\frac{71}{105},$$

$$T_{10}(8)=\frac{1+1+3\cdot \frac{71}{105}}{8}=\frac{141}{280}=0.503571\ldots$$

Summing all contributions from $q=2$ to $10$ yields

$$S(10)=\frac{3485}{504}=6.914682539682\ldots,$$

which matches the published checkpoint.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they build a smallest-prime-factor table up to $N$, which makes factorization of each $q$ fast and also supports a linear-time computation of $\varphi(q)$. Next they precompute the harmonic prefix array $H_0,H_1,\dots,H_N$. Then each value of $q$ is processed independently: factor $q$, enumerate the squarefree divisors induced by its distinct prime factors, evaluate the Möbius formulas for $A_q$ and $C_q$, apply the appropriate branch, and add the result to the global sum.

The C++ and Java implementations additionally use compensated summation to reduce floating-point drift, and the C++ version can split the outer loop into parallel chunks. None of those engineering choices changes the mathematics; they only improve numerical stability and runtime for the target size.

Complexity Analysis

Building the smallest-prime-factor table, the totient table, and the harmonic prefix array costs $O(N)$ time and $O(N)$ memory. For each $q$, the factorization is recovered from the precomputed table, and the Möbius expansion enumerates exactly $2^{\omega(q)}$ squarefree divisors, where $\omega(q)$ is the number of distinct prime factors of $q$. Therefore the total running time is

$$O\!\left(N+\sum_{q=2}^{N}2^{\omega(q)}\right),$$

with $O(N)$ memory usage. For the specific target $N=10^7$, every integer has at most eight distinct prime factors, so the subset enumeration remains very small in practice. Parallel execution reduces wall-clock time but does not change the asymptotic work.

Footnotes and References

Problem page: https://projecteuler.net/problem=441
Möbius function and inversion: Wikipedia — Möbius function
Harmonic numbers: Wikipedia — Harmonic number
Euler's totient function: Wikipedia — Euler's totient function
Compensated summation: Wikipedia — Kahan summation algorithm

Problem 441 source code

C++

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

namespace {

using i64 = std::int64_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

constexpr u32 kDefaultN = 10'000'000U;
constexpr u32 kParallelThreshold = 200'000U;
constexpr u32 kChunkSize = 2'048U;

struct Options {
    u32 n = kDefaultN;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
    bool run_checks = true;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value_out) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

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

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

    value_out = parsed;
    return true;
}

bool parse_u32_after_prefix(const std::string& arg, const char* prefix, u32& value_out) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<u32>::max())) {
        return false;
    }
    value_out = static_cast<u32>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value_out) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value_out = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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-checks") {
            options.run_checks = false;
            continue;
        }

        u32 parsed_u32 = 0U;
        if (parse_u32_after_prefix(arg, "--n=", parsed_u32)) {
            options.n = parsed_u32;
            continue;
        }

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

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

    if (options.n < 2U) {
        std::cerr << "--n must be at least 2.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const u64 workload_units) {
    if (!allow_multithreading || workload_units < 2ULL) {
        return 1U;
    }

    const unsigned max_threads = static_cast<unsigned>(
        std::min<u64>(workload_units, static_cast<u64>(std::numeric_limits<unsigned>::max())));

    if (requested_threads > 0U) {
        return std::max(1U, std::min(requested_threads, max_threads));
    }

    if (workload_units < static_cast<u64>(kParallelThreshold)) {
        return 1U;
    }

    unsigned hardware_threads = std::thread::hardware_concurrency();
    if (hardware_threads == 0U) {
        hardware_threads = 1U;
    }

    return std::max(1U, std::min(hardware_threads, max_threads));
}

struct KahanAccumulator {
    long double sum = 0.0L;
    long double compensation = 0.0L;

    void add(const long double value) {
        const long double y = value - compensation;
        const long double t = sum + y;
        compensation = (t - sum) - y;
        sum = t;
    }
};

class Euler441Solver {
public:
    explicit Euler441Solver(const u32 n)
        : n_(n),
          half_(n / 2U),
          spf_(static_cast<std::size_t>(n) + 1U, 0U),
          phi_(static_cast<std::size_t>(n) + 1U, 0U),
          harmonic_(static_cast<std::size_t>(n) + 1U, 0.0L) {
        build_linear_sieve();
        build_harmonic_prefix();
    }

    long double solve(const unsigned thread_count) const {
        if (n_ < 2U) {
            return 0.0L;
        }

        const unsigned clamped_threads =
            std::max(1U, std::min(thread_count, static_cast<unsigned>(n_ - 1U)));

        if (clamped_threads == 1U) {
            KahanAccumulator acc;
            for (u32 q = 2U; q <= n_; ++q) {
                acc.add(compute_contribution(q));
            }
            return acc.sum;
        }

        std::atomic<u32> next_q(2U);
        std::vector<long double> partial(clamped_threads, 0.0L);
        std::vector<std::thread> workers;
        workers.reserve(clamped_threads);

        for (unsigned tid = 0U; tid < clamped_threads; ++tid) {
            workers.emplace_back([&, tid]() {
                KahanAccumulator local;

                while (true) {
                    const u32 q_begin = next_q.fetch_add(kChunkSize, std::memory_order_relaxed);
                    if (q_begin > n_) {
                        break;
                    }

                    const u32 q_end = std::min<u32>(n_, q_begin + kChunkSize - 1U);
                    for (u32 q = q_begin; q <= q_end; ++q) {
                        local.add(compute_contribution(q));
                    }
                }

                partial[static_cast<std::size_t>(tid)] = local.sum;
            });
        }

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

        KahanAccumulator total;
        for (const long double value : partial) {
            total.add(value);
        }
        return total.sum;
    }

private:
    static constexpr int kMaxPrimeFactors = 16;
    static constexpr int kMaxSquarefreeDivisors = 1024;

    u32 n_ = 0U;
    u32 half_ = 0U;
    std::vector<u32> spf_;
    std::vector<u32> phi_;
    std::vector<long double> harmonic_;

    void build_linear_sieve() {
        std::vector<u32> primes;
        primes.reserve(static_cast<std::size_t>(n_) / 10U + 16U);

        phi_[1] = 1U;

        for (u32 i = 2U; i <= n_; ++i) {
            if (spf_[i] == 0U) {
                spf_[i] = i;
                phi_[i] = i - 1U;
                primes.push_back(i);
            }

            for (const u32 p : primes) {
                const u64 v = static_cast<u64>(i) * static_cast<u64>(p);
                if (v > n_) {
                    break;
                }

                spf_[static_cast<std::size_t>(v)] = p;
                if (p == spf_[i]) {
                    phi_[static_cast<std::size_t>(v)] = phi_[i] * p;
                    break;
                }

                phi_[static_cast<std::size_t>(v)] = phi_[i] * (p - 1U);
            }
        }
    }

    void build_harmonic_prefix() {
        harmonic_[0] = 0.0L;
        for (u32 i = 1U; i <= n_; ++i) {
            harmonic_[static_cast<std::size_t>(i)] =
                harmonic_[static_cast<std::size_t>(i - 1U)] + 1.0L / static_cast<long double>(i);
        }
    }

    long double compute_contribution(const u32 q) const {
        u32 prime_factors[kMaxPrimeFactors];
        int factor_count = 0;

        u32 x = q;
        while (x > 1U) {
            const u32 p = spf_[x];
            prime_factors[factor_count++] = p;
            while (x % p == 0U) {
                x /= p;
            }
        }

        u32 divisors[kMaxSquarefreeDivisors];
        int signs[kMaxSquarefreeDivisors];
        int subset_count = 1;

        divisors[0] = 1U;
        signs[0] = 1;

        for (int i = 0; i < factor_count; ++i) {
            const u32 p = prime_factors[i];
            for (int j = 0; j < subset_count; ++j) {
                divisors[subset_count + j] = divisors[j] * p;
                signs[subset_count + j] = -signs[j];
            }
            subset_count <<= 1;
        }

        const u32 qm1 = q - 1U;
        long double a_full = 0.0L;
        for (int i = 0; i < subset_count; ++i) {
            const u32 d = divisors[i];
            const long double signed_inv_d =
                static_cast<long double>(signs[i]) / static_cast<long double>(d);
            a_full += signed_inv_d * harmonic_[static_cast<std::size_t>(qm1 / d)];
        }

        if (q <= half_) {
            return (static_cast<long double>(phi_[q]) + a_full) / static_cast<long double>(q);
        }

        const u32 a = n_ - q;
        long double a_a = 0.0L;
        long double b = 0.0L;
        i64 c = 0;

        for (int i = 0; i < subset_count; ++i) {
            const u32 d = divisors[i];
            const int sign = signs[i];
            const long double signed_inv_d =
                static_cast<long double>(sign) / static_cast<long double>(d);
            const u32 a_over_d = a / d;

            c += static_cast<i64>(sign) * static_cast<i64>(a_over_d);
            const long double ha = harmonic_[static_cast<std::size_t>(a_over_d)];
            const long double hq = harmonic_[static_cast<std::size_t>(qm1 / d)];

            a_a += signed_inv_d * ha;
            b += signed_inv_d * (hq - ha);
        }

        const long double numerator =
            static_cast<long double>(c) + a_a + static_cast<long double>(a + 1U) * b;
        return numerator / static_cast<long double>(q);
    }
};

long double brute_force_s(const u32 n) {
    long double total = 0.0L;

    for (u32 i = 2U; i <= n; ++i) {
        long double r_i = 0.0L;
        for (u32 p = 1U; p < i; ++p) {
            for (u32 q = p + 1U; q <= i; ++q) {
                if (p + q < i) {
                    continue;
                }
                if (std::gcd(p, q) != 1U) {
                    continue;
                }
                r_i += 1.0L / (static_cast<long double>(p) * static_cast<long double>(q));
            }
        }
        total += r_i;
    }

    return total;
}

bool near_abs(const long double a, const long double b, const long double eps) {
    return std::fabsl(a - b) <= eps;
}

bool run_problem_checkpoints() {
    struct Checkpoint {
        u32 n;
        long double expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {2U, 0.5L},
        {10U, 6.9146825396825396825L},
        {100U, 58.296238062166323935L},
    };

    for (const Checkpoint& checkpoint : checkpoints) {
        const Euler441Solver solver(checkpoint.n);
        const long double got = solver.solve(1U);
        if (!near_abs(got, checkpoint.expected, 1e-12L)) {
            std::cerr << "Checkpoint failed for n=" << checkpoint.n << ": got "
                      << std::setprecision(18) << static_cast<double>(got) << ", expected "
                      << static_cast<double>(checkpoint.expected) << '\n';
            return false;
        }
    }

    return true;
}

bool run_bruteforce_crosscheck() {
    for (u32 n = 2U; n <= 40U; ++n) {
        const Euler441Solver solver(n);
        const long double fast = solver.solve(1U);
        const long double brute = brute_force_s(n);
        if (!near_abs(fast, brute, 1e-12L)) {
            std::cerr << "Bruteforce mismatch at n=" << n << ": fast=" << std::setprecision(18)
                      << static_cast<double>(fast) << ", brute=" << static_cast<double>(brute)
                      << '\n';
            return false;
        }
    }
    return true;
}

bool run_parallel_consistency_check() {
    unsigned hardware_threads = std::thread::hardware_concurrency();
    if (hardware_threads < 2U) {
        return true;
    }

    const u32 n = 200'000U;
    const Euler441Solver solver(n);
    const long double single = solver.solve(1U);

    const unsigned threads = std::max(2U, std::min(8U, hardware_threads));
    const long double parallel = solver.solve(threads);

    if (!near_abs(single, parallel, 1e-12L)) {
        std::cerr << "Parallel consistency check failed: single=" << std::setprecision(18)
                  << static_cast<double>(single) << ", parallel="
                  << static_cast<double>(parallel) << ", threads=" << threads << '\n';
        return false;
    }

    return true;
}

bool run_validation_suite(const bool allow_multithreading) {
    if (!run_problem_checkpoints()) {
        return false;
    }
    if (!run_bruteforce_crosscheck()) {
        return false;
    }
    if (allow_multithreading && !run_parallel_consistency_check()) {
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checks && !run_validation_suite(options.allow_multithreading)) {
        return 1;
    }

    const Euler441Solver solver(options.n);
    const unsigned threads = choose_thread_count(options.allow_multithreading,
                                                 options.requested_threads,
                                                 static_cast<u64>(options.n) - 1ULL);
    const long double answer = solver.solve(threads);

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

Python

import math

def solve():
    n = 10000000

    # Linear sieve for SPF and phi
    spf = [0]*(n+1); phi = [0]*(n+1); phi[1] = 1
    primes = []
    for i in range(2, n+1):
        if spf[i] == 0:
            spf[i] = i; phi[i] = i-1; primes.append(i)
        for p in primes:
            if i*p > n: break
            spf[i*p] = p
            if p == spf[i]: phi[i*p] = phi[i]*p; break
            phi[i*p] = phi[i]*(p-1)

    # Harmonic prefix
    H = [0.0]*(n+1)
    for i in range(1, n+1): H[i] = H[i-1] + 1.0/i

    half = n // 2
    total = 0.0
    for q in range(2, n+1):
        # Get distinct prime factors
        pf = []; x = q
        while x > 1:
            p = spf[x]; pf.append(p)
            while x % p == 0: x //= p

        # Generate squarefree divisors with signs
        divs = [1]; signs = [1]
        for p in pf:
            m = len(divs)
            for j in range(m):
                divs.append(divs[j]*p); signs.append(-signs[j])

        qm1 = q - 1
        if q <= half:
            a_full = 0.0
            for i in range(len(divs)):
                a_full += signs[i] / divs[i] * H[qm1 // divs[i]]
            total += (phi[q] + a_full) / q
        else:
            a = n - q
            a_a = b = 0.0; c = 0
            for i in range(len(divs)):
                d = divs[i]; s = signs[i]; inv_d = s / d
                ad = a // d; c += s * ad
                a_a += inv_d * H[ad]; b += inv_d * (H[qm1 // d] - H[ad])
            total += (c + a_a + (a+1)*b) / q

    return f'{total:.4f}'

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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;

public class Euler441 {
    static final int N = 10000000;
    static final int HALF = N / 2;

    static int[] spf;
    static int[] phi;
    static double[] harmonic;

    static void buildSieveAndHarmonic() {
        spf = new int[N + 1];
        phi = new int[N + 1];
        harmonic = new double[N + 1];

        phi[1] = 1;
        List<Integer> primes = new ArrayList<>();

        for (int i = 2; i <= N; i++) {
            if (spf[i] == 0) {
                spf[i] = i;
                phi[i] = i - 1;
                primes.add(i);
            }

            for (int p : primes) {
                long v = (long) i * p;
                if (v > N)
                    break;
                spf[(int) v] = p;
                if (p == spf[i]) {
                    phi[(int) v] = phi[i] * p;
                    break;
                }
                phi[(int) v] = phi[i] * (p - 1);
            }
        }

        harmonic[0] = 0.0;
        for (int i = 1; i <= N; i++) {
            harmonic[i] = harmonic[i - 1] + 1.0 / i;
        }
    }

    static class KahanAccumulator {
        double sum = 0.0;
        double compensation = 0.0;

        void add(double value) {
            double y = value - compensation;
            double t = sum + y;
            compensation = (t - sum) - y;
            sum = t;
        }
    }

    static double computeContribution(int q) {
        int[] primeFactors = new int[16];
        int factorCount = 0;

        int x = q;
        while (x > 1) {
            int p = spf[x];
            primeFactors[factorCount++] = p;
            while (x % p == 0) {
                x /= p;
            }
        }

        int[] divisors = new int[1024];
        int[] signs = new int[1024];
        int subsetCount = 1;

        divisors[0] = 1;
        signs[0] = 1;

        for (int i = 0; i < factorCount; i++) {
            int p = primeFactors[i];
            for (int j = 0; j < subsetCount; j++) {
                divisors[subsetCount + j] = divisors[j] * p;
                signs[subsetCount + j] = -signs[j];
            }
            subsetCount <<= 1;
        }

        int qm1 = q - 1;

        if (q <= HALF) {
            double aFull = 0.0;
            for (int i = 0; i < subsetCount; i++) {
                int d = divisors[i];
                double signedInvD = (double) signs[i] / d;
                aFull += signedInvD * harmonic[qm1 / d];
            }
            return (phi[q] + aFull) / q;
        }

        int a = N - q;
        double aA = 0.0;
        double b = 0.0;
        long c = 0;

        for (int i = 0; i < subsetCount; i++) {
            int d = divisors[i];
            int sign = signs[i];
            double signedInvD = (double) sign / d;
            int aOverD = a / d;

            c += (long) sign * aOverD;
            double ha = harmonic[aOverD];
            double hq = harmonic[qm1 / d];

            aA += signedInvD * ha;
            b += signedInvD * (hq - ha);
        }

        double numerator = c + aA + (a + 1) * b;
        return numerator / q;
    }

    public static String solve() {
        buildSieveAndHarmonic();

        int numThreads = Math.max(1, Runtime.getRuntime().availableProcessors());
        int chunkSize = 2048;
        int totalChunks = (N - 1 + chunkSize - 1) / chunkSize;
        double[] partials = new double[totalChunks];

        IntStream.range(0, totalChunks).parallel().forEach(chunkIdx -> {
            int qBegin = 2 + chunkIdx * chunkSize;
            int qEnd = Math.min(N, qBegin + chunkSize - 1);
            KahanAccumulator local = new KahanAccumulator();
            for (int q = qBegin; q <= qEnd; q++) {
                local.add(computeContribution(q));
            }
            partials[chunkIdx] = local.sum;
        });

        KahanAccumulator total = new KahanAccumulator();
        for (double p : partials) {
            total.add(p);
        }

        return String.format(java.util.Locale.US, "%.4f", total.sum);
    }

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