Problem 664: An Infinite Game

View on Project Euler

Project Euler Problem 664 Solution

EulerSolve provides an optimized solution for Project Euler Problem 664, An Infinite Game, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For Problem 664, the quantity coming from the infinite game is not evaluated by simulating positions directly. The C++, Python, and Java implementations all use an equivalent closed form built from the golden ratio $$\varphi=\frac{1+\sqrt5}{2}, \qquad A_n=\sum_{d=1}^{\infty}\frac{d^n}{\varphi^d}.$$ Once this series is known, the required answer is $$F(n)=3+\left\lceil \log_{\varphi}(A_n)\right\rceil.$$ The challenge is numerical rather than combinatorial. For large \(n\), the dominant terms of the series are enormous, the tails are tiny, and a naive evaluation would overflow or waste work on irrelevant terms. The solution therefore works in the logarithmic domain and keeps only the narrow peak region that materially contributes to the sum. Mathematical Approach Step 1: Rewrite the series in logarithmic form Let $$a_d=\frac{d^n}{\varphi^d}.$$ Instead of forming \(a_d\) directly, take logarithms: $$L(d)=\log a_d=n\log d-d\log\varphi.$$ Then the infinite series becomes $$A_n=\sum_{d\ge 1} e^{L(d)}.$$ This removes the overflow risk from the dominant terms and turns the problem into computing a stable logarithm of a sum. Step 2: Locate the dominant index The mass of the series is concentrated near the maximum of \(L(d)\)....

Detailed mathematical approach

Problem Summary

For Problem 664, the quantity coming from the infinite game is not evaluated by simulating positions directly. The C++, Python, and Java implementations all use an equivalent closed form built from the golden ratio

$$\varphi=\frac{1+\sqrt5}{2}, \qquad A_n=\sum_{d=1}^{\infty}\frac{d^n}{\varphi^d}.$$

Once this series is known, the required answer is

$$F(n)=3+\left\lceil \log_{\varphi}(A_n)\right\rceil.$$

The challenge is numerical rather than combinatorial. For large \(n\), the dominant terms of the series are enormous, the tails are tiny, and a naive evaluation would overflow or waste work on irrelevant terms. The solution therefore works in the logarithmic domain and keeps only the narrow peak region that materially contributes to the sum.

Mathematical Approach

Step 1: Rewrite the series in logarithmic form

Let

$$a_d=\frac{d^n}{\varphi^d}.$$

Instead of forming \(a_d\) directly, take logarithms:

$$L(d)=\log a_d=n\log d-d\log\varphi.$$

Then the infinite series becomes

$$A_n=\sum_{d\ge 1} e^{L(d)}.$$

This removes the overflow risk from the dominant terms and turns the problem into computing a stable logarithm of a sum.

Step 2: Locate the dominant index

The mass of the series is concentrated near the maximum of \(L(d)\). Treating the index as a real variable \(x\), define

$$f(x)=n\log x-x\log\varphi.$$

Its derivative is

$$f'(x)=\frac{n}{x}-\log\varphi,$$

so the stationary point satisfies

$$x_*=\frac{n}{\log\varphi}.$$

Because

$$f''(x)=-\frac{n}{x^2}\lt 0,$$

this point is a strict maximum. The best integer \(d\) must therefore lie very close to \(n/\log\varphi\), which is why the implementation only needs a tiny local scan around that estimate to find the true peak term.

Step 3: Truncate the infinite sum to a finite window

Once the peak value \(L_{\max}\) is known, terms far away from it are numerically irrelevant. The implementations expand left and right from the peak until

$$L_{\max}-L(d)\gt 120.$$

At that point the term is smaller than the peak term by a factor of \(e^{-120}\), and the remaining tail only shrinks further because the profile is already descending on that side. This gives a practical finite window

$$d\in[L,R]$$

that contains all numerically significant contributions.

Step 4: Apply the log-sum-exp identity

Directly summing \(e^{L(d)}\) would still be unsafe near the peak. The stable reformulation is

$$\log A_n=L_{\max}+\log\left(\sum_{d=L}^{R} e^{L(d)-L_{\max}}\right).$$

Now every exponent \(L(d)-L_{\max}\) is non-positive, so the largest rescaled term is exactly \(1\) and all others lie in \((0,1]\). This is the central numerical device used by all three implementations.

Step 5: Recover the required integer answer

After \(\log A_n\) has been computed, convert from base \(e\) to base \(\varphi\):

$$\log_{\varphi}(A_n)=\frac{\log A_n}{\log\varphi}.$$

The answer is then

$$F(n)=3+\left\lceil \frac{\log A_n}{\log\varphi}\right\rceil.$$

A tiny negative epsilon is subtracted before the ceiling is taken. This protects against a floating-point boundary case where the true value is an integer but the computed value lands infinitesimally above it.

Worked Example: \(n=2\)

For \(n=2\), the series can be checked exactly. Set

$$x=\frac{1}{\varphi}.$$

Then

$$A_2=\sum_{d=1}^{\infty} d^2 x^d.$$

The classical generating-function identity gives

$$\sum_{d=1}^{\infty} d^2 x^d=\frac{x(x+1)}{(1-x)^3}.$$

Since \(x=1/\varphi\), we have \(x+1=\varphi\) and \(1-x=1/\varphi^2\). Therefore

$$A_2=\frac{(1/\varphi)\varphi}{(1/\varphi^2)^3}=\varphi^6.$$

So

$$F(2)=3+\left\lceil \log_{\varphi}(\varphi^6)\right\rceil=3+6=9,$$

which matches the checkpoint verified by the implementation.

How the Code Works

The implementation first computes \(\log\varphi\) once and reuses it everywhere. It estimates the peak position by \(n/\log\varphi\), inspects a very small neighborhood to choose the best integer index, and then expands outward until the 120-log-unit cutoff is reached on both sides.

Next it accumulates the rescaled sum

$$\sum_{d=L}^{R} e^{L(d)-L_{\max}}.$$

The C++ version can split this interval into chunks and sum them in parallel. The Python and Java versions follow the same mathematics in a sequential pass, with Python using high-precision decimal logarithms and exponentials and Java using floating-point arithmetic. All three then take the logarithm of the scaled total, divide by \(\log\varphi\), apply the guarded ceiling, and add \(3\).

Complexity Analysis

If the retained window has width

$$W=R-L+1,$$

then building the window and summing its terms both cost \(O(W)\) time. The extra memory is \(O(1)\) in the sequential versions and \(O(T)\) partial accumulators in the threaded C++ version, where \(T\) is the number of worker threads.

Near the maximum, a quadratic expansion gives

$$f(x)\approx f(x_*)-\frac{(\log\varphi)^2}{2n}(x-x_*)^2,$$

so a fixed drop threshold produces a window whose width is on the order of \(\sqrt{n}\). In practice the retained band is therefore far smaller than the peak index itself.

Footnotes and References

  1. Project Euler problem page: https://projecteuler.net/problem=664
  2. Golden ratio: Wikipedia — Golden ratio
  3. LogSumExp: Wikipedia — LogSumExp
  4. Polylogarithm: Wikipedia — Polylogarithm
  5. Generating function: Wikipedia — Generating function

Problem 664 source code

C++

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

namespace {

using u64 = std::uint64_t;

constexpr u64 kTargetN = 1'234'567ULL;
constexpr long double kSqrt5 = 2.23606797749978969640917366873127623544L;
constexpr long double kPhi = (1.0L + kSqrt5) / 2.0L;
constexpr long double kLogPhi = 0.48121182505960344749775891342436842314L;
constexpr long double kTailDrop = 120.0L;

struct Options {
    u64 n = kTargetN;
    bool run_checkpoints = true;
    unsigned requested_threads = 0;
};

struct Window {
    u64 left = 1;
    u64 right = 1;
    u64 peak = 1;
    long double peak_log_term = 0.0L;
};

long double log_term(u64 n, u64 d) {
    return static_cast<long double>(n) * std::log(static_cast<long double>(d)) -
           static_cast<long double>(d) * kLogPhi;
}

unsigned resolve_thread_count(unsigned requested_threads) {
    if (requested_threads > 0) {
        return requested_threads;
    }
    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 1;
    }
    return threads;
}

Window build_window(u64 n) {
    Window out;

    u64 approx = 1;
    if (n > 0) {
        const long double center = static_cast<long double>(n) / kLogPhi;
        if (center > static_cast<long double>(std::numeric_limits<u64>::max() - 10ULL)) {
            approx = std::numeric_limits<u64>::max() - 10ULL;
        } else {
            approx = static_cast<u64>(center);
            if (approx == 0) {
                approx = 1;
            }
        }
    }

    out.peak = approx;
    out.peak_log_term = log_term(n, out.peak);

    const u64 scan_lo = (approx > 6 ? approx - 6 : 1ULL);
    const u64 scan_hi = approx + 6;
    for (u64 d = scan_lo; d <= scan_hi; ++d) {
        const long double value = log_term(n, d);
        if (value > out.peak_log_term) {
            out.peak_log_term = value;
            out.peak = d;
        }
    }

    out.left = out.peak;
    while (out.left > 1) {
        const long double next_value = log_term(n, out.left - 1);
        if (out.peak_log_term - next_value > kTailDrop) {
            break;
        }
        --out.left;
    }

    out.right = out.peak;
    while (out.right < std::numeric_limits<u64>::max() - 1ULL) {
        const long double next_value = log_term(n, out.right + 1);
        if (out.peak_log_term - next_value > kTailDrop) {
            break;
        }
        ++out.right;
    }

    return out;
}

long double scaled_sum_range(u64 n,
                             u64 left,
                             u64 right,
                             long double peak_log_term,
                             unsigned threads) {
    if (left > right) {
        return 0.0L;
    }

    const u64 total_terms = right - left + 1ULL;
    if (threads <= 1 || total_terms < 50'000ULL) {
        long double total = 0.0L;
        for (u64 d = left; d <= right; ++d) {
            total += std::exp(log_term(n, d) - peak_log_term);
        }
        return total;
    }

    threads = std::min<unsigned>(threads, static_cast<unsigned>(total_terms));

    std::atomic<u64> next{left};
    std::vector<long double> partial(threads, 0.0L);
    std::vector<std::thread> workers;
    workers.reserve(threads);

    constexpr u64 kChunk = 1024ULL;

    for (unsigned tid = 0; tid < threads; ++tid) {
        workers.emplace_back([&, tid]() {
            long double local = 0.0L;
            while (true) {
                const u64 begin = next.fetch_add(kChunk, std::memory_order_relaxed);
                if (begin > right) {
                    break;
                }
                const u64 end = std::min(right, begin + kChunk - 1ULL);
                for (u64 d = begin; d <= end; ++d) {
                    local += std::exp(log_term(n, d) - peak_log_term);
                }
            }
            partial[tid] = local;
        });
    }

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

    long double total = 0.0L;
    for (long double value : partial) {
        total += value;
    }
    return total;
}

long double compute_log_series(u64 n, unsigned threads) {
    const Window window = build_window(n);

    const long double scaled =
        scaled_sum_range(n, window.left, window.right, window.peak_log_term, threads);

    return window.peak_log_term + std::log(scaled);
}

u64 compute_f(u64 n, unsigned threads) {
    const long double log_series = compute_log_series(n, threads);
    const long double raw = log_series / kLogPhi;

    // The closed form is F(n) = 3 + ceil(log_phi(sum_{d>=1} d^n / phi^d)).
    const long double rounded = std::ceil(raw - 1e-18L) + 3.0L;
    return static_cast<u64>(rounded);
}

bool run_validation_checkpoints(unsigned threads) {
    struct Checkpoint {
        u64 n;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {0ULL, 4ULL},
        {1ULL, 6ULL},
        {2ULL, 9ULL},
        {3ULL, 13ULL},
        {11ULL, 58ULL},
        {123ULL, 1173ULL},
    };

    for (const Checkpoint cp : checkpoints) {
        const u64 got = compute_f(cp.n, threads);
        if (got != cp.expected) {
            std::cerr << "Checkpoint failed: F(" << cp.n << ") expected " << cp.expected
                      << ", got " << got << "\n";
            return false;
        }
    }

    if (threads > 1) {
        const u64 thread_check_n = 200'000ULL;
        const u64 single = compute_f(thread_check_n, 1);
        const u64 multi = compute_f(thread_check_n, threads);
        if (single != multi) {
            std::cerr << "Thread consistency failed for n=" << thread_check_n << ": single="
                      << single << ", multi=" << multi << "\n";
            return false;
        }
    }

    std::cout << "Validation checkpoints passed.\n";
    return true;
}

bool parse_u64(const std::string& s, u64& out) {
    try {
        std::size_t pos = 0;
        const unsigned long long value = std::stoull(s, &pos, 10);
        if (pos != s.size()) {
            return false;
        }
        out = static_cast<u64>(value);
        return true;
    } catch (...) {
        return false;
    }
}

Options parse_options(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 (arg.rfind("--n=", 0) == 0) {
            u64 value = 0;
            if (!parse_u64(arg.substr(4), value)) {
                std::cerr << "Invalid --n value: " << arg << "\n";
                std::exit(1);
            }
            options.n = value;
            continue;
        }

        if (arg.rfind("--threads=", 0) == 0) {
            u64 value = 0;
            if (!parse_u64(arg.substr(10), value) || value == 0ULL ||
                value > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
                std::cerr << "Invalid --threads value: " << arg << "\n";
                std::exit(1);
            }
            options.requested_threads = static_cast<unsigned>(value);
            continue;
        }

        std::cerr << "Unknown argument: " << arg << "\n";
        std::exit(1);
    }

    return options;
}

}  // namespace

int main(int argc, char** argv) {
    const Options options = parse_options(argc, argv);
    const unsigned threads = resolve_thread_count(options.requested_threads);

    if (options.run_checkpoints) {
        if (!run_validation_checkpoints(threads)) {
            return 1;
        }
    }

    const u64 answer = compute_f(options.n, threads);
    std::cout << answer << '\n';
    return 0;
}

Python

import math
from decimal import Decimal, getcontext

getcontext().prec = 50

kLogPhi = Decimal('1.61803398874989484820458683436563811772030917980576').ln()
kTailDrop = Decimal('120.0')

def log_term(n, d):
    dn = Decimal(n)
    dd = Decimal(d)
    return dn * dd.ln() - dd * kLogPhi

def build_window(n):
    approx = int(Decimal(n) / kLogPhi)
    if approx == 0:
        approx = 1

    scan_lo = max(1, approx - 6)
    scan_hi = approx + 6

    peak = scan_lo
    peak_log_term = log_term(n, peak)

    for d in range(scan_lo + 1, scan_hi + 1):
        val = log_term(n, d)
        if val > peak_log_term:
            peak_log_term = val
            peak = d

    left = peak
    while left > 1:
        next_val = log_term(n, left - 1)
        if (peak_log_term - next_val) > kTailDrop:
            break
        left -= 1

    right = peak
    while True:
        next_val = log_term(n, right + 1)
        if (peak_log_term - next_val) > kTailDrop:
            break
        right += 1

    return left, right, peak_log_term

def compute_f(n):
    if n == 0: return 4
    if n == 1: return 6
    left, right, peak_log_term = build_window(n)

    scaled_sum = Decimal('0')
    for d in range(left, right + 1):
        diff = log_term(n, d) - peak_log_term
        scaled_sum += diff.exp()

    log_series = peak_log_term + scaled_sum.ln()
    raw = log_series / kLogPhi

    # Using Python's math.ceil on float conversion for raw
    # Since 50 precision is used, converting to float and taking ceil matches C++ logic
    val = float(raw)
    return int(math.ceil(val - 1e-15)) + 3

def solve():
    ans = compute_f(1234567)
    return str(ans)

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

Java

public class Euler664 {

    static final double kSqrt5 = 2.23606797749978969640917366873127623544;
    static final double kPhi = (1.0 + kSqrt5) / 2.0;
    static final double kLogPhi = Math.log(kPhi);
    static final double kTailDrop = 120.0;

    static double logTerm(long n, long d) {
        return (double) n * Math.log((double) d) - (double) d * kLogPhi;
    }

    static class Window {
        long left, right, peak;
        double peakLogTerm;
    }

    static Window buildWindow(long n) {
        Window out = new Window();

        long approx = 1;
        if (n > 0) {
            double center = (double) n / kLogPhi;
            approx = (long) center;
            if (approx == 0)
                approx = 1;
        }

        out.peak = approx;
        out.peakLogTerm = logTerm(n, out.peak);

        long scanLo = Math.max(1L, approx - 6);
        long scanHi = approx + 6;

        for (long d = scanLo; d <= scanHi; ++d) {
            double value = logTerm(n, d);
            if (value > out.peakLogTerm) {
                out.peakLogTerm = value;
                out.peak = d;
            }
        }

        out.left = out.peak;
        while (out.left > 1) {
            double nextVal = logTerm(n, out.left - 1);
            if (out.peakLogTerm - nextVal > kTailDrop)
                break;
            out.left--;
        }

        out.right = out.peak;
        while (true) {
            double nextVal = logTerm(n, out.right + 1);
            if (out.peakLogTerm - nextVal > kTailDrop)
                break;
            out.right++;
        }

        return out;
    }

    static long computeF(long n) {
        Window w = buildWindow(n);

        double total = 0.0;
        for (long d = w.left; d <= w.right; ++d) {
            total += Math.exp(logTerm(n, d) - w.peakLogTerm);
        }

        double logSeries = w.peakLogTerm + Math.log(total);
        double raw = logSeries / kLogPhi;

        long rounded = (long) Math.ceil(raw - 1e-15) + 3;
        return rounded;
    }

    public static String solve() {
        long ans = computeF(1234567L);
        return Long.toString(ans);
    }

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