Problem 568: Reciprocal Games II

Project Euler Problem 568 Solution

EulerSolve provides an optimized solution for Project Euler Problem 568, Reciprocal Games II, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The quantity from the problem can be reduced to $$D(n)=\frac{H_n}{2^n}, \qquad H_n=\sum_{k=1}^{n}\frac{1}{k},$$ and the task is to return the first seven significant digits of $D(123456789)$. The main difficulty is numerical rather than combinatorial: $2^n$ is unimaginably large, $D(n)$ is tiny, and only the leading digits matter. So the solution never forms $2^n$ directly. Instead it works with logarithms, a sharp asymptotic expansion for the harmonic number, and extra precision in the final subtraction. Mathematical Approach The whole method is built around the observation that significant digits are encoded in the fractional part of a base-10 logarithm. Step 1: Turn Leading Digits into a Logarithm Problem Let $$L=\log_{10} D(n).$$ Write $L=a+f$ where $a=\lfloor L \rfloor$ is the integer part and $$f=L-\lfloor L \rfloor \in [0,1).$$ Then $$D(n)=10^L=10^a\cdot 10^f.$$ The factor $10^a$ only shifts the decimal point, while $10^f\in[1,10)$ is the significand in scientific notation. Therefore the first seven significant digits are $$\left\lfloor 10^{f+6} \right\rfloor.$$ So the problem is reduced to computing the fractional part of $L$ accurately....

Detailed mathematical approach

Problem Summary

The quantity from the problem can be reduced to

$$D(n)=\frac{H_n}{2^n}, \qquad H_n=\sum_{k=1}^{n}\frac{1}{k},$$

and the task is to return the first seven significant digits of $D(123456789)$. The main difficulty is numerical rather than combinatorial: $2^n$ is unimaginably large, $D(n)$ is tiny, and only the leading digits matter. So the solution never forms $2^n$ directly. Instead it works with logarithms, a sharp asymptotic expansion for the harmonic number, and extra precision in the final subtraction.

Mathematical Approach

The whole method is built around the observation that significant digits are encoded in the fractional part of a base-10 logarithm.

Step 1: Turn Leading Digits into a Logarithm Problem

Let

$$L=\log_{10} D(n).$$

Write $L=a+f$ where $a=\lfloor L \rfloor$ is the integer part and

$$f=L-\lfloor L \rfloor \in [0,1).$$

Then

$$D(n)=10^L=10^a\cdot 10^f.$$

The factor $10^a$ only shifts the decimal point, while $10^f\in[1,10)$ is the significand in scientific notation. Therefore the first seven significant digits are

$$\left\lfloor 10^{f+6} \right\rfloor.$$

So the problem is reduced to computing the fractional part of $L$ accurately.

Step 2: Express $L$ Using the Harmonic Number

From the closed form,

$$L=\log_{10} H_n-n\log_{10} 2.$$

The second term is easy in principle, but it is numerically dominant because $n=123456789$ is large. The first term grows only like $\log \log n$. That means the final answer depends on the low-order digits of a subtraction between quantities of very different scales, so ordinary floating-point arithmetic must be handled carefully.

Step 3: Approximate $H_n$ with Euler-Maclaurin

For large $n$, the harmonic number is approximated by the Euler-Maclaurin expansion

$$H_n=\ln n+\gamma+\frac{1}{2n}-\frac{1}{12n^2}+\frac{1}{120n^4}-\frac{1}{252n^6}+\frac{1}{240n^8}-\frac{5}{660n^{10}}+O(n^{-12}),$$

where $\gamma$ is the Euler-Mascheroni constant. For the target input, the neglected tail is astronomically small and cannot affect the first seven significant digits. One implementation also keeps a direct summation branch for modest $n$, but the actual Project Euler input is far into the asymptotic regime.

Step 4: Preserve the Fractional Part During the Critical Subtraction

Once $H_n$ is known, we still need

$$\log_{10} H_n-n\log_{10} 2.$$

The C++ and Java implementations store the important quantity as a high part plus a low correction and carry out compensated addition, subtraction, and multiplication so the fractional information is not lost when $n\log_{10} 2$ is formed. The Python implementation reaches the same goal with high-precision decimal arithmetic. After the subtraction, the fractional part is normalized back into $[0,1)$, because the small correction term can move it slightly below $0$ or above $1$.

Step 5: Convert the Fractional Part Back to Digits

After obtaining $f$, the leading digits come from

$$x=10^{f+6}.$$

Ideally the answer is $\lfloor x \rfloor$. In practice, the implementations add a tiny positive safety margin before flooring so that a value such as $3828124.9999999996$ is not rounded down incorrectly when the true mathematical result is $3828125$. The result is also capped at $9999999$, because $10^f<10$.

Worked Example: $n=6$

The statement example is small enough to do exactly:

$$H_6=1+\frac12+\frac13+\frac14+\frac15+\frac16=\frac{49}{20}.$$

Hence

$$D(6)=\frac{49/20}{2^6}=\frac{49}{1280}=0.03828125=3.828125\times 10^{-2}.$$

So the first seven significant digits are plainly

$$3828125.$$

The logarithmic method gives the same result: if $L=\log_{10}(0.03828125)$, then $f=L-\lfloor L \rfloor=\log_{10}(3.828125)$, and therefore

$$10^{f+6}=3.828125\times 10^6=3828125.$$

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical pipeline. They start from $D(n)=H_n/2^n$, evaluate $H_n$, compute $\log_{10} H_n$, subtract $n\log_{10} 2$, extract the fractional part, and then turn that fraction into seven leading digits with one final power of ten.

The difference is in the precision strategy. The Python implementation uses a high-precision decimal context and evaluates the asymptotic series directly. The C++ and Java implementations treat the critical logarithmic subtraction more carefully by representing the value as two linked floating-point components, which keeps more reliable fractional information than a single double would. One implementation also includes a direct harmonic summation path for smaller inputs and checks the sample $n=6$ before printing the answer for $123456789$.

Complexity Analysis

For the actual target input, the harmonic number is evaluated from a fixed number of asymptotic terms, the logarithmic arithmetic uses a fixed number of operations, and the final digit extraction is constant work. So the running time is $O(1)$ and the memory usage is $O(1)$. If the optional direct summation path is used for smaller $n$, that branch takes $O(n)$ time and still uses $O(1)$ memory.

Footnotes and References

Problem page: Project Euler 568
Harmonic numbers: Wikipedia - Harmonic number
Euler-Maclaurin formula: Wikipedia - Euler-Maclaurin formula
Common logarithm: Wikipedia - Common logarithm
Euler-Mascheroni constant: Wikipedia - Euler-Mascheroni constant
Double-double arithmetic overview: Wikipedia - Double-double arithmetic

Problem 568 source code

C++

#include <cassert>
#include <cstdint>
#include <cmath>
#include <iostream>
#include <limits>

namespace {

using u64 = std::uint64_t;

struct dd {
    double hi;
    double lo;
};

// Error-free transforms for double-double arithmetic.
static inline dd two_sum(double a, double b) {
    double s = a + b;
    double bb = s - a;
    double err = (a - (s - bb)) + (b - bb);
    return {s, err};
}

static inline dd quick_two_sum(double a, double b) {
    // Assumes |a| >= |b|.
    double s = a + b;
    double err = b - (s - a);
    return {s, err};
}

static inline dd two_prod(double a, double b) {
    double p = a * b;
    // Dekker split for IEEE-754 double.
    static constexpr double split = 134217729.0;  // 2^27 + 1
    double ca = split * a;
    double a_hi = ca - (ca - a);
    double a_lo = a - a_hi;
    double cb = split * b;
    double b_hi = cb - (cb - b);
    double b_lo = b - b_hi;
    double err = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;
    return {p, err};
}

static inline dd dd_add(dd x, dd y) {
    dd s = two_sum(x.hi, y.hi);
    double e = x.lo + y.lo;
    dd t = two_sum(s.lo, e);
    dd u = quick_two_sum(s.hi, t.hi);
    double lo = u.lo + t.lo;
    return quick_two_sum(u.hi, lo);
}

static inline dd dd_sub(dd x, dd y) { return dd_add(x, {-y.hi, -y.lo}); }

static inline dd dd_mul_d(dd x, double y) {
    dd p = two_prod(x.hi, y);
    double e = x.lo * y;
    dd s = two_sum(p.lo, e);
    dd u = quick_two_sum(p.hi, s.hi);
    double lo = u.lo + s.lo;
    return quick_two_sum(u.hi, lo);
}

static inline dd to_dd(long double x) {
    double hi = static_cast<double>(x);
    long double rem = x - static_cast<long double>(hi);
    double lo = static_cast<double>(rem);
    return {hi, lo};
}

static long double harmonic_asymptotic(u64 n) {
    // Euler–Maclaurin:
    // H_n = ln(n) + gamma + 1/(2n) - 1/(12n^2) + 1/(120n^4) - 1/(252n^6) + ...
    static constexpr long double gamma = 0.57721566490153286060651209008240243104215933593992L;
    long double nn = static_cast<long double>(n);
    long double inv = 1.0L / nn;
    long double inv2 = inv * inv;
    long double inv4 = inv2 * inv2;
    long double inv6 = inv4 * inv2;
    long double inv8 = inv4 * inv4;
    long double inv10 = inv8 * inv2;
    return logl(nn) + gamma + 0.5L * inv - (1.0L / 12.0L) * inv2 + (1.0L / 120.0L) * inv4 -
           (1.0L / 252.0L) * inv6 + (1.0L / 240.0L) * inv8 - (5.0L / 660.0L) * inv10;
}

static u64 leading7_D(u64 n) {
    // D(n) = J_B(n) - J_A(n) = H_n / 2^n.
    long double Hn = 0;
    if (n <= 2000000ULL) {
        for (u64 k = 1; k <= n; ++k) Hn += 1.0L / static_cast<long double>(k);
    } else {
        Hn = harmonic_asymptotic(n);
    }

    // log10(2) split into hi+lo so n*log10(2) keeps enough fractional precision.
    static constexpr dd log10_2 = {0.3010299956639812, -4.786261105275507e-18};

    dd log10_H = to_dd(log10l(Hn));
    dd nlog10_2 = dd_mul_d(log10_2, static_cast<double>(n));
    dd L = dd_sub(log10_H, nlog10_2);

    // frac(L) in [0,1).
    double frac = L.hi - std::floor(L.hi);
    frac += L.lo;
    frac -= std::floor(frac);

    long double x = powl(10.0L, static_cast<long double>(frac) + 6.0L);
    long double eps =
        8.0L * std::numeric_limits<long double>::epsilon() * fabsl(x);  // n=6 needs ~1e-8.
    u64 ans = static_cast<u64>(floorl(x + eps));
    if (ans >= 10000000ULL) ans = 9999999ULL;  // mantissa is in [1,10)
    return ans;
}

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    // Statement example: D(6)=0.03828125 -> leading 7 digits are 3828125.
    assert(leading7_D(6) == 3828125ULL);

    std::cout << leading7_D(123456789ULL) << '\n';
    return 0;
}

Python

import decimal

def solve():
    decimal.getcontext().prec = 60
    n = 123456789
    
    gamma = decimal.Decimal("0.57721566490153286060651209008240243104215933593992")
    nn = decimal.Decimal(n)
    inv = decimal.Decimal('1') / nn
    inv2 = inv * inv
    inv4 = inv2 * inv2
    inv6 = inv4 * inv2
    inv8 = inv4 * inv4
    inv10 = inv8 * inv2
    
    Hn = nn.ln() + gamma + decimal.Decimal("0.5") * inv \
         - (decimal.Decimal("1") / decimal.Decimal("12")) * inv2 \
         + (decimal.Decimal("1") / decimal.Decimal("120")) * inv4 \
         - (decimal.Decimal("1") / decimal.Decimal("252")) * inv6 \
         + (decimal.Decimal("1") / decimal.Decimal("240")) * inv8 \
         - (decimal.Decimal("5") / decimal.Decimal("660")) * inv10
         
    log10_Hn = Hn.log10()
    log10_2 = decimal.Decimal("2").log10()
    n_log10_2 = decimal.Decimal(n) * log10_2
    
    L = log10_Hn - n_log10_2
    
    frac = L - L.to_integral_value(rounding=decimal.ROUND_FLOOR)
    
    x = (decimal.Decimal(10) ** frac) * decimal.Decimal(1000000)
    ans = int(x + decimal.Decimal("1e-15"))
    if ans >= 10000000:
        ans = 9999999
        
    return str(ans)

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

Java

public class Euler568 {

    static class DD {
        double hi;
        double lo;

        DD(double hi, double lo) {
            this.hi = hi;
            this.lo = lo;
        }
    }

    static DD twoSum(double a, double b) {
        double s = a + b;
        double bb = s - a;
        double err = (a - (s - bb)) + (b - bb);
        return new DD(s, err);
    }

    static DD quickTwoSum(double a, double b) {
        double s = a + b;
        double err = b - (s - a);
        return new DD(s, err);
    }

    static DD twoProd(double a, double b) {
        double p = a * b;
        double split = 134217729.0;
        double ca = split * a;
        double a_hi = ca - (ca - a);
        double a_lo = a - a_hi;
        double cb = split * b;
        double b_hi = cb - (cb - b);
        double b_lo = b - b_hi;
        double err = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;
        return new DD(p, err);
    }

    static DD ddAdd(DD x, DD y) {
        DD s = twoSum(x.hi, y.hi);
        double e = x.lo + y.lo;
        DD t = twoSum(s.lo, e);
        DD u = quickTwoSum(s.hi, t.hi);
        double lo = u.lo + t.lo;
        return quickTwoSum(u.hi, lo);
    }

    static DD ddSub(DD x, DD y) {
        return ddAdd(x, new DD(-y.hi, -y.lo));
    }

    static DD ddMulD(DD x, double y) {
        DD p = twoProd(x.hi, y);
        double e = x.lo * y;
        DD s = twoSum(p.lo, e);
        DD u = quickTwoSum(p.hi, s.hi);
        double lo = u.lo + s.lo;
        return quickTwoSum(u.hi, lo);
    }

    static double harmonicAsymptotic(long n) {
        double gamma = 0.577215664901532860606512090082402431042159;
        double nn = (double) n;
        double inv = 1.0 / nn;
        double inv2 = inv * inv;
        double inv4 = inv2 * inv2;
        double inv6 = inv4 * inv2;
        double inv8 = inv4 * inv4;
        double inv10 = inv8 * inv2;
        return Math.log(nn) + gamma + 0.5 * inv - (1.0 / 12.0) * inv2 + (1.0 / 120.0) * inv4 - (1.0 / 252.0) * inv6
                + (1.0 / 240.0) * inv8 - (5.0 / 660.0) * inv10;
    }

    public static String solve() {
        long n = 123456789L;
        double Hn = harmonicAsymptotic(n);

        DD log10_2 = new DD(0.3010299956639812, -4.786261105275507e-18);
        DD log10_H = new DD(Math.log10(Hn), 0.0);
        DD nlog10_2 = ddMulD(log10_2, (double) n);
        DD L = ddSub(log10_H, nlog10_2);

        double frac = L.hi - Math.floor(L.hi);
        frac += L.lo;
        frac -= Math.floor(frac);

        double x = Math.pow(10.0, frac + 6.0);
        double eps = 8.0 * Math.ulp(1.0) * Math.abs(x);
        long ans = (long) Math.floor(x + eps);
        if (ans >= 10000000L)
            ans = 9999999L;

        return String.valueOf(ans);
    }

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