Problem 787: Bézout's Game

Project Euler Problem 787 Solution

EulerSolve provides an optimized solution for Project Euler Problem 787, Bézout's Game, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a fixed $n$, the game positions are the coprime ordered pairs $(a,b)$ of positive integers with $a+b\le n$. A legal move goes to a smaller positive pair $(u,v)$ satisfying $$1\le u<a,\qquad 1\le v<b,\qquad |av-bu|=1.$$ The game convention used here treats every position with $a=1$ or $b=1$ as an immediate win. The goal is to count how many admissible positions are winning when $n=10^9$. Mathematical Approach Let $W(n)$ be the number of winning positions and $T(n)$ the number of all primitive positions. The fast solution first counts every primitive pair, then subtracts the pairs that are losing. Step 1: Count All Primitive Positions Fix the sum $s=a+b$. Then $1\le a\le s-1$ and $b=s-a$, so $$\gcd(a,b)=\gcd(a,s).$$ Therefore the number of coprime ordered pairs with $a+b=s$ is exactly $\varphi(s)$. Summing over every possible total gives $$T(n)=\sum_{s=2}^{n}\varphi(s)=\Phi(n)-1,$$ where $$\Phi(n)=\sum_{m\le n}\varphi(m).$$ The implementations evaluate $\Phi(n)$ through the standard Möbius identity $$\Phi(n)=\frac{1+\sum_{d=1}^{n}\mu(d)\left\lfloor\frac{n}{d}\right\rfloor^2}{2},$$ which turns the problem into prefix queries for the Möbius function. Step 2: Characterize the Losing Positions Because the move rule is symmetric in $a$ and $b$, it is enough to understand the ordered half $1<a<b$....

Detailed mathematical approach

Problem Summary

For a fixed $n$, the game positions are the coprime ordered pairs $(a,b)$ of positive integers with $a+b\le n$. A legal move goes to a smaller positive pair $(u,v)$ satisfying

$$1\le u<a,\qquad 1\le v<b,\qquad |av-bu|=1.$$

The game convention used here treats every position with $a=1$ or $b=1$ as an immediate win. The goal is to count how many admissible positions are winning when $n=10^9$.

Mathematical Approach

Let $W(n)$ be the number of winning positions and $T(n)$ the number of all primitive positions. The fast solution first counts every primitive pair, then subtracts the pairs that are losing.

Step 1: Count All Primitive Positions

Fix the sum $s=a+b$. Then $1\le a\le s-1$ and $b=s-a$, so

$$\gcd(a,b)=\gcd(a,s).$$

Therefore the number of coprime ordered pairs with $a+b=s$ is exactly $\varphi(s)$. Summing over every possible total gives

$$T(n)=\sum_{s=2}^{n}\varphi(s)=\Phi(n)-1,$$

where

$$\Phi(n)=\sum_{m\le n}\varphi(m).$$

The implementations evaluate $\Phi(n)$ through the standard Möbius identity

$$\Phi(n)=\frac{1+\sum_{d=1}^{n}\mu(d)\left\lfloor\frac{n}{d}\right\rfloor^2}{2},$$

which turns the problem into prefix queries for the Möbius function.

Step 2: Characterize the Losing Positions

Because the move rule is symmetric in $a$ and $b$, it is enough to understand the ordered half $1<a<b$. Consider such a state and one of its legal children $(u,v)$. From

$$av-bu=\pm1$$

we obtain

$$v-u=\frac{(b-a)u\pm1}{a}\ge 0.$$

So every child also lies in the ordered half, except for the boundary case $(u,v)=(1,1)$, which is already a base win.

Now argue by induction on $a+b$:

If $a$ is even, then $b$ is odd because $\gcd(a,b)=1$. Reducing $av-bu=\pm1$ modulo $2$ gives $u\equiv1\pmod 2$. Hence every child has odd smaller coordinate, so no child belongs to the losing family described below. Therefore the state is losing.

If $a$ is odd and $a>1$, the congruences $bu\equiv \pm1\pmod a$ give two residues $u$ and $a-u$ in $\{1,\dots,a-1\}$. Since $a$ is odd, exactly one of them is even. Choosing that even residue produces a legal child with $u<v$ and, by parity, $v$ odd. Its smaller coordinate is even, so by induction that child is losing. Therefore the original state is winning.

Thus the losing positions in the ordered half are exactly the primitive pairs

$$ (a,b)=(2x,\,2y+1),\qquad 1\le x\le y,\qquad \gcd(2x,2y+1)=1. $$

Step 3: Count the Ordered Losing Shape Without the Coprime Condition

Define $A(q)$ to be the number of pairs of the form $(2x,2y+1)$ with $1\le x\le y$ and $2x+2y+1\le q$, ignoring coprimality for the moment. Set

$$s=\left\lfloor\frac{q-1}{2}\right\rfloor.$$

Then the constraints become

$$1\le x\le y,\qquad x+y\le s.$$

For fixed $x$, the variable $y$ runs from $x$ to $s-x$, so

$$A(q)=\sum_{x=1}^{\lfloor s/2\rfloor}(s-2x+1).$$

This simplifies to the closed form

$$A(q)= \begin{cases} p^2, & s=2p,\\ p(p+1), & s=2p+1. \end{cases}$$

An equivalent version, closer to the implementation, is obtained from $k=\left\lceil q/2\right\rceil$ and $p=\left\lfloor k/2\right\rfloor$:

$$A(q)= \begin{cases} p^2, & k\text{ is odd},\\ p(p-1), & k\text{ is even}. \end{cases}$$

Step 4: Enforce Coprimality with Möbius Inversion Over Odd Divisors

Every pair counted by $A(q)$ has one even coordinate and one odd coordinate, so any common divisor must itself be odd. That makes the Möbius inversion especially clean:

$$L_{1/2}(n)=\sum_{\substack{d\le n\\ d\text{ odd}}}\mu(d)\,A\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right),$$

where $L_{1/2}(n)$ is the number of losing positions in the ordered half $a<b$.

To evaluate this quickly, define the ordinary Mertens function

$$M(x)=\sum_{m\le x}\mu(m)$$

and the odd-only prefix

$$M_{\mathrm{odd}}(x)=\sum_{\substack{d\le x\\ d\text{ odd}}}\mu(d).$$

The implementations use the identity

$$M_{\mathrm{odd}}(x)=\sum_{j\ge0} M\!\left(\left\lfloor\frac{x}{2^j}\right\rfloor\right),$$

which follows by expanding the right-hand side and observing that the contributions of even integers cancel in pairs, leaving only odd divisors.

Step 5: Use Symmetry and Assemble the Final Formula

Swapping coordinates preserves both the move rule and the winning/losing status. Since primitive pairs never lie on the diagonal except for $(1,1)$, which is already winning, the total number of losing positions is

$$L(n)=2L_{1/2}(n).$$

Therefore

$$\boxed{W(n)=\left(\sum_{s=2}^{n}\varphi(s)\right)-2\sum_{\substack{d\le n\\ d\text{ odd}}}\mu(d)\,A\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right).}$$

This is exactly the arithmetic formula implemented by the fast solution.

Worked Example: $n=20$

First count all primitive states:

$$T(20)=\sum_{s=2}^{20}\varphi(s)=127.$$

Now count the ordered losing shape without enforcing coprimality. Here

$$A(20)=20,$$

coming from the twenty ordered pairs $(2x,2y+1)$ with $2x<2y+1$ and $2x+2y+1\le20$.

Among these, the only non-primitive family with odd gcd comes from $d=3$, because

$$A\!\left(\left\lfloor\frac{20}{3}\right\rfloor\right)=A(6)=1,$$

while $A(\lfloor20/d\rfloor)=0$ for every larger odd $d$. Hence

$$L_{1/2}(20)=A(20)-A(6)=20-1=19,$$

$$L(20)=2\cdot19=38,$$

and finally

$$W(20)=127-38=89.$$

This matches direct recursive play on the small instance.

How the Code Works

The C++, Python, and Java implementations all follow the same numerical plan. They first build a presieved Möbius table up to a fixed bound and store its prefix sums. That gives instant access to $M(x)$ for small $x$.

For larger arguments, the implementation uses the classical divisor-block recurrence

$$M(x)=1-\sum_{l=2}^{x}(r-l+1)\,M\!\left(\left\lfloor\frac{x}{l}\right\rfloor\right),$$

where each block $[l,r]$ shares the same quotient $\left\lfloor x/l\right\rfloor$. Memoization ensures that every large prefix value is computed only once.

The total number of primitive states is then obtained from the summatory-totient identity by grouping equal values of $\left\lfloor n/d\right\rfloor$. The losing half is evaluated with the same quotient blocking, but using prefix differences of $M_{\mathrm{odd}}$ and the explicit closed form for $A(q)$.

After doubling the ordered-half losing count, the implementation subtracts it from the total primitive count. One version also checks the formula on small inputs with direct recursion before printing the large answer, but the actual $n=10^9$ computation is purely arithmetic.

Complexity Analysis

Let $B$ be the presieve limit. Building the Möbius sieve and its prefix sums costs $O(B)$ time and $O(B)$ memory. The outer summations for both the totient total and the losing count run over the distinct values of $\left\lfloor n/d\right\rfloor$, which is $O(\sqrt n)$.

The expensive part is answering large Mertens-prefix queries. Because those queries are memoized and each one is decomposed into quotient blocks as well, the practical running time stays far below linear in $n$. Memory usage is dominated by the presieve arrays plus the cache of already-computed prefix values.

Footnotes and References

Problem page: https://projecteuler.net/problem=787
Euler's totient function: Wikipedia — Euler's totient function
Möbius function: Wikipedia — Möbius function
Mertens function: Wikipedia — Mertens function
Möbius inversion formula: Wikipedia — Möbius inversion formula
Farey sequence: Wikipedia — Farey sequence

Problem 787 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <iostream>
#include <numeric>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

using i64 = long long;
using i128 = __int128_t;

static std::string to_string_i128(i128 x) {
    if (x == 0) {
        return "0";
    }
    bool neg = x < 0;
    if (neg) {
        x = -x;
    }
    std::string s;
    while (x > 0) {
        const int d = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + d));
        x /= 10;
    }
    if (neg) {
        s.push_back('-');
    }
    std::reverse(s.begin(), s.end());
    return s;
}

struct MertensPrefix {
    int limit;
    std::vector<int> pref_mu;
    std::unordered_map<i64, i64> memo_mu;
    std::unordered_map<i64, i64> memo_odd_mu;

    explicit MertensPrefix(int lim) : limit(lim), pref_mu(lim + 1, 0) {
        std::vector<int> mu(lim + 1, 0);
        std::vector<int> primes;
        std::vector<bool> composite(lim + 1, false);

        mu[1] = 1;
        for (int i = 2; i <= lim; ++i) {
            if (!composite[i]) {
                primes.push_back(i);
                mu[i] = -1;
            }
            for (int p : primes) {
                const i64 v = 1LL * i * p;
                if (v > lim) {
                    break;
                }
                composite[static_cast<std::size_t>(v)] = true;
                if (i % p == 0) {
                    mu[static_cast<std::size_t>(v)] = 0;
                    break;
                }
                mu[static_cast<std::size_t>(v)] = -mu[i];
            }
        }

        for (int i = 1; i <= lim; ++i) {
            pref_mu[i] = pref_mu[i - 1] + mu[i];
        }
        memo_mu.reserve(1 << 16);
        memo_odd_mu.reserve(1 << 16);
    }

    i64 mu_prefix(i64 n) {
        if (n <= limit) {
            return pref_mu[static_cast<std::size_t>(n)];
        }
        const auto it = memo_mu.find(n);
        if (it != memo_mu.end()) {
            return it->second;
        }

        i64 ans = 1;
        i64 l = 2;
        while (l <= n) {
            const i64 q = n / l;
            const i64 r = n / q;
            ans -= (r - l + 1) * mu_prefix(q);
            l = r + 1;
        }

        memo_mu[n] = ans;
        return ans;
    }

    i64 odd_mu_prefix(i64 n) {
        if (n <= 0) {
            return 0;
        }
        const auto it = memo_odd_mu.find(n);
        if (it != memo_odd_mu.end()) {
            return it->second;
        }

        i64 ans = 0;
        i64 cur = n;
        while (cur > 0) {
            ans += mu_prefix(cur);
            cur >>= 1;
        }

        memo_odd_mu[n] = ans;
        return ans;
    }
};

static i128 sum_totients(i64 n, MertensPrefix& mp) {
    i128 s = 0;
    i64 l = 1;
    while (l <= n) {
        const i64 q = n / l;
        const i64 r = n / q;
        const i64 mu_seg = mp.mu_prefix(r) - mp.mu_prefix(l - 1);
        s += static_cast<i128>(mu_seg) * static_cast<i128>(q) * static_cast<i128>(q);
        l = r + 1;
    }
    return (s + 1) / 2;
}

static i64 f_odd_floor_sum(i64 m) {
    const i64 k = (m + 1) / 2;
    const i64 p = k / 2;
    if (k & 1LL) {
        return p * p;
    }
    return p * (p - 1);
}

static i128 solve_fast(i64 n) {
    const int sieve_limit = 2'000'000;
    MertensPrefix mp(sieve_limit);

    const i128 total = sum_totients(n, mp) - 1;

    i128 losing_half = 0;
    i64 l = 1;
    while (l <= n) {
        const i64 q = n / l;
        const i64 r = n / q;
        const i64 odd_mu_seg = mp.odd_mu_prefix(r) - mp.odd_mu_prefix(l - 1);
        losing_half += static_cast<i128>(odd_mu_seg) * static_cast<i128>(f_odd_floor_sum(q));
        l = r + 1;
    }

    const i128 losing = 2 * losing_half;
    return total - losing;
}

static bool has_losing_child(
    int a,
    int b,
    std::unordered_map<i64, bool>& memo,
    const std::function<bool(int, int)>& win
) {
    for (int u = 1; u < a; ++u) {
        const int m1 = b * u - 1;
        if (m1 % a == 0) {
            const int v = m1 / a;
            if (v > 0 && v < b) {
                if (!win(u, v)) {
                    return true;
                }
            }
        }

        const int m2 = b * u + 1;
        if (m2 % a == 0) {
            const int v = m2 / a;
            if (v > 0 && v < b) {
                if (!win(u, v)) {
                    return true;
                }
            }
        }
    }
    return false;
}

static i128 solve_bruteforce_game(int n) {
    std::unordered_map<i64, bool> memo;
    memo.reserve(1 << 16);

    std::function<bool(int, int)> win = [&](int a, int b) -> bool {
        const i64 key = (static_cast<i64>(a) << 32) | static_cast<unsigned int>(b);
        const auto it = memo.find(key);
        if (it != memo.end()) {
            return it->second;
        }

        if (a == 1 || b == 1) {
            memo[key] = true;
            return true;
        }

        const bool w = has_losing_child(a, b, memo, win);
        memo[key] = w;
        return w;
    };

    i128 ans = 0;
    for (int a = 1; a < n; ++a) {
        for (int b = 1; a + b <= n; ++b) {
            if (std::gcd(a, b) != 1) {
                continue;
            }
            if (win(a, b)) {
                ++ans;
            }
        }
    }
    return ans;
}

int main() {
    assert(solve_fast(4) == 5);
    assert(solve_fast(100) == 2043);
    assert(solve_fast(40) == solve_bruteforce_game(40));

    const i128 ans = solve_fast(1'000'000'000LL);
    std::cout << to_string_i128(ans) << '\n';
    return 0;
}

Python

class MertensPrefix:
    def __init__(self, lim):
        self.limit = lim
        self.pref_mu = [0] * (lim + 1)
        
        mu = [0] * (lim + 1)
        primes = []
        composite = [False] * (lim + 1)
        
        mu[1] = 1
        for i in range(2, lim + 1):
            if not composite[i]:
                primes.append(i)
                mu[i] = -1
            for p in primes:
                v = i * p
                if v > lim:
                    break
                composite[v] = True
                if i % p == 0:
                    mu[v] = 0
                    break
                mu[v] = -mu[i]
                
        for i in range(1, lim + 1):
            self.pref_mu[i] = self.pref_mu[i - 1] + mu[i]
            
        self.memo_mu = {}
        self.memo_odd_mu = {}

    def mu_prefix(self, n):
        if n <= self.limit:
            return self.pref_mu[n]
        if n in self.memo_mu:
            return self.memo_mu[n]

        ans = 1
        l = 2
        while l <= n:
            q = n // l
            r = n // q
            ans -= (r - l + 1) * self.mu_prefix(q)
            l = r + 1

        self.memo_mu[n] = ans
        return ans

    def odd_mu_prefix(self, n):
        if n <= 0:
            return 0
        if n in self.memo_odd_mu:
            return self.memo_odd_mu[n]

        ans = 0
        cur = n
        while cur > 0:
            ans += self.mu_prefix(cur)
            cur >>= 1

        self.memo_odd_mu[n] = ans
        return ans

def sum_totients(n, mp):
    s = 0
    l = 1
    while l <= n:
        q = n // l
        r = n // q
        mu_seg = mp.mu_prefix(r) - mp.mu_prefix(l - 1)
        s += mu_seg * q * q
        l = r + 1
    return (s + 1) // 2

def f_odd_floor_sum(m):
    k = (m + 1) // 2
    p = k // 2
    if k % 2 == 1:
        return p * p
    return p * (p - 1)

def solve_fast(n):
    sieve_limit = 2000000
    mp = MertensPrefix(sieve_limit)

    total = sum_totients(n, mp) - 1

    losing_half = 0
    l = 1
    while l <= n:
        q = n // l
        r = n // q
        odd_mu_seg = mp.odd_mu_prefix(r) - mp.odd_mu_prefix(l - 1)
        losing_half += odd_mu_seg * f_odd_floor_sum(q)
        l = r + 1

    losing = 2 * losing_half
    return total - losing

def solve():
    return str(solve_fast(10**9))

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

Java

import java.util.ArrayList;
import java.util.HashMap;

public class Euler787 {

    static class MertensPrefix {
        int limit;
        int[] prefMu;
        HashMap<Long, Long> memoMu;
        HashMap<Long, Long> memoOddMu;

        MertensPrefix(int lim) {
            limit = lim;
            prefMu = new int[lim + 1];

            int[] mu = new int[lim + 1];
            ArrayList<Integer> primes = new ArrayList<>();
            boolean[] composite = new boolean[lim + 1];

            mu[1] = 1;
            for (int i = 2; i <= lim; ++i) {
                if (!composite[i]) {
                    primes.add(i);
                    mu[i] = -1;
                }
                for (int p : primes) {
                    long v = (long) i * p;
                    if (v > lim)
                        break;
                    composite[(int) v] = true;
                    if (i % p == 0) {
                        mu[(int) v] = 0;
                        break;
                    }
                    mu[(int) v] = -mu[i];
                }
            }

            for (int i = 1; i <= lim; ++i) {
                prefMu[i] = prefMu[i - 1] + mu[i];
            }

            memoMu = new HashMap<>(); // Primitives mappings would be faster but this requires less external libraries
            memoOddMu = new HashMap<>();
        }

        long muPrefix(long n) {
            if (n <= limit) {
                return prefMu[(int) n];
            }
            Long cached = memoMu.get(n);
            if (cached != null) {
                return cached;
            }

            long ans = 1;
            long l = 2;
            while (l <= n) {
                long q = n / l;
                long r = n / q;
                ans -= (r - l + 1) * muPrefix(q);
                l = r + 1;
            }

            memoMu.put(n, ans);
            return ans;
        }

        long oddMuPrefix(long n) {
            if (n <= 0)
                return 0;
            Long cached = memoOddMu.get(n);
            if (cached != null)
                return cached;

            long ans = 0;
            long cur = n;
            while (cur > 0) {
                ans += muPrefix(cur);
                cur >>= 1;
            }

            memoOddMu.put(n, ans);
            return ans;
        }
    }

    static long sumTotients(long n, MertensPrefix mp) {
        long s = 0;
        long l = 1;
        while (l <= n) {
            long q = n / l;
            long r = n / q;
            long muSeg = mp.muPrefix(r) - mp.muPrefix(l - 1);
            s += muSeg * q * q;
            l = r + 1;
        }
        return (s + 1) / 2;
    }

    static long fOddFloorSum(long m) {
        long k = (m + 1) / 2;
        long p = k / 2;
        if ((k & 1) == 1) {
            return p * p;
        }
        return p * (p - 1);
    }

    static long solveFast(long n) {
        int sieveLimit = 2000000;
        MertensPrefix mp = new MertensPrefix(sieveLimit);

        long total = sumTotients(n, mp) - 1;

        long losingHalf = 0;
        long l = 1;
        while (l <= n) {
            long q = n / l;
            long r = n / q;
            long oddMuSeg = mp.oddMuPrefix(r) - mp.oddMuPrefix(l - 1);
            losingHalf += oddMuSeg * fOddFloorSum(q);
            l = r + 1;
        }

        long losing = 2 * losingHalf;
        return total - losing;
    }

    public static String solve() {
        return Long.toString(solveFast(1000000000L));
    }

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