Problem 560: Coprime Nim

Project Euler Problem 560 Solution

EulerSolve provides an optimized solution for Project Euler Problem 560, Coprime Nim, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each pile size $n$, one move chooses a number $a$ with $1 \le a \le n$ and $\gcd(a,n)=1$, then replaces the pile by $n-a$. For fixed $N$ and $K$, we consider ordered $K$-tuples of pile sizes drawn from $\{1,2,\dots,N-1\}$. The task is to count how many such $K$-pile positions are losing positions, meaning that the xor of their Grundy values is $0$, and to return the answer modulo $10^9+7$. A direct dynamic program over all $K$-pile states is hopeless for the intended input size, so the solution first classifies the one-pile game completely and then turns the multi-pile count into an xor convolution problem. Mathematical Approach Let $g(n)$ be the Grundy value of one pile of size $n$. The move rule gives $$g(n)=\operatorname{mex}\left\{g(n-a): 1 \le a \le n,\ \gcd(a,n)=1\right\}.$$ If we know the distribution of these one-pile Grundy values on $\{1,\dots,N-1\}$, then Sprague-Grundy theory reduces the full problem to counting xor sums. Step 1: Express the Game in Sprague-Grundy Form The position with $K$ piles $(x_1,\dots,x_K)$ is losing exactly when $$g(x_1)\oplus g(x_2)\oplus \cdots \oplus g(x_K)=0.$$ So the main difficulty is to understand $g(n)$ for one pile. Once that is known, each pile contributes only its Grundy label, and the counting problem becomes purely combinatorial....

Detailed mathematical approach

Problem Summary

For each pile size $n$, one move chooses a number $a$ with $1 \le a \le n$ and $\gcd(a,n)=1$, then replaces the pile by $n-a$. For fixed $N$ and $K$, we consider ordered $K$-tuples of pile sizes drawn from $\{1,2,\dots,N-1\}$. The task is to count how many such $K$-pile positions are losing positions, meaning that the xor of their Grundy values is $0$, and to return the answer modulo $10^9+7$.

A direct dynamic program over all $K$-pile states is hopeless for the intended input size, so the solution first classifies the one-pile game completely and then turns the multi-pile count into an xor convolution problem.

Mathematical Approach

Let $g(n)$ be the Grundy value of one pile of size $n$. The move rule gives

$$g(n)=\operatorname{mex}\left\{g(n-a): 1 \le a \le n,\ \gcd(a,n)=1\right\}.$$

If we know the distribution of these one-pile Grundy values on $\{1,\dots,N-1\}$, then Sprague-Grundy theory reduces the full problem to counting xor sums.

Step 1: Express the Game in Sprague-Grundy Form

The position with $K$ piles $(x_1,\dots,x_K)$ is losing exactly when

$$g(x_1)\oplus g(x_2)\oplus \cdots \oplus g(x_K)=0.$$

So the main difficulty is to understand $g(n)$ for one pile. Once that is known, each pile contributes only its Grundy label, and the counting problem becomes purely combinatorial.

Step 2: Classify the One-Pile Grundy Values

The implementations use the closed form

$$g(1)=1,$$

$$g(n)=0 \quad \text{for even } n,$$

$$g(n)=\pi(\operatorname{spf}(n)) \quad \text{for odd } n>1,$$

where $\operatorname{spf}(n)$ is the smallest prime factor of $n$, and $\pi(t)$ is the number of primes at most $t$, so for a prime $p$, $\pi(p)$ is its $1$-based prime index.

This formula follows from the move structure:

For $n=1$, the only legal move is to $0$, so the mex is $1$.

If $n$ is even, every legal $a$ must be odd, hence every reachable value $n-a$ is a positive odd number. Every positive odd number has positive Grundy value, so $0$ is missing from the move set and therefore $g(n)=0$.

Now let $n>1$ be odd and write $p=\operatorname{spf}(n)$. The move $a=1$ reaches $n-1$, which is even, so Grundy value $0$ is present. The move $a=n-1$ reaches $1$, so Grundy value $1$ is present. More generally, for every odd prime $q<p$, the move to the pile $q$ is legal because

$$\gcd(n,n-q)=\gcd(n,q)=1,$$

since $q$ does not divide $n$. That target has Grundy value $\pi(q)$. Hence every value below $\pi(p)$ occurs among the legal moves.

On the other hand, no legal move can land on an odd number whose smallest prime factor is exactly $p$. If such a target $m$ existed, then $p\mid n$ and $p\mid m$, so

$$\gcd(n,n-m)=\gcd(n,m)\ge p,$$

contradicting the coprimality condition on the move. Therefore Grundy value $\pi(p)$ is absent, and the mex is exactly $\pi(p)$.

Larger Grundy values may also appear among some moves, but they do not matter once all smaller values are present and $\pi(p)$ is missing.

Step 3: Build the Frequency Table of Grundy Values

Define

$$c_t=\#\left\{x\in\{1,\dots,N-1\}: g(x)=t\right\}.$$

Then $c_t$ is the number of pile sizes carrying Grundy label $t$. The classification above immediately implies

$$c_0=\left\lfloor\frac{N-1}{2}\right\rfloor,$$

because exactly the even pile sizes have Grundy value $0$, and

$$c_1=1,$$

because the only pile with Grundy value $1$ is $1$ itself. Every odd $x>1$ contributes to the bucket indexed by the prime index of its smallest prime factor.

The implementations obtain these buckets with a linear sieve that stores the smallest prime factor of every integer up to $N-1$, together with the running prime-count function $\pi$.

Step 4: Turn $K$ Piles into an XOR Convolution

Let $L(N,K)$ denote the number of losing ordered $K$-tuples. If one pile contributes Grundy value $u$ in $c_u$ ways and another contributes $v$ in $c_v$ ways, then the pair contributes xor value $u\oplus v$. This is exactly xor convolution.

For two arrays $f$ and $h$, define

$$\left(f *_{\oplus} h\right)(s)=\sum_{u\oplus v=s} f(u)h(v).$$

Therefore the distribution of xor sums for $K$ piles is the $K$-fold xor convolution of $c$ with itself, and the desired answer is

$$L(N,K)=c^{(*K)}(0).$$

This identity is just Sprague-Grundy theory translated into counting language: pile choices are independent, and a position is losing exactly when the xor sum is $0$.

Step 5: Evaluate the XOR Convolution with the Walsh-Hadamard Transform

Directly forming $K$ xor convolutions would still be too slow, so the implementations switch to the Walsh-Hadamard transform for xor convolution. If $\widehat{c}$ denotes the xor transform of $c$, then

$$\widehat{c^{(*K)}}(i)=\widehat{c}(i)^K.$$

So the whole computation becomes:

Pad the frequency vector to a power-of-two length $m$.

Apply the xor Walsh-Hadamard transform.

Raise each transformed component to the $K$-th power modulo $10^9+7$.

Apply the same transform again.

Multiply every entry by $m^{-1}\pmod{10^9+7}$, because for xor convolution the transform is self-inverse up to the factor $m$.

The coefficient at index $0$ is then exactly $L(N,K)$.

Worked Example: $L(5,2)$

For $N=5$, the allowed pile sizes are $1,2,3,4$. Their Grundy values are

$$g(1)=1,\qquad g(2)=0,\qquad g(3)=2,\qquad g(4)=0.$$

Hence the frequency vector is

$$c_0=2,\qquad c_1=1,\qquad c_2=1.$$

With two piles, xor is $0$ exactly when both piles have the same Grundy value, so

$$L(5,2)=c_0^2+c_1^2+c_2^2=2^2+1^2+1^2=6.$$

This matches the checkpoint used by the implementation and illustrates why the entire problem is governed by the Grundy frequency distribution rather than by the raw pile values.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. First they build a linear sieve up to $N-1$, storing for every number its smallest prime factor and also the prime-count prefix needed to convert a smallest prime factor into its prime index. Then they construct the Grundy frequency vector directly from the closed form above: all even sizes contribute to Grundy $0$, the pile $1$ contributes to Grundy $1$, and every odd size greater than $1$ is placed into the bucket determined by the prime index of its smallest prime factor.

Next the implementation copies that frequency vector into a power-of-two buffer and applies the xor Walsh-Hadamard transform. Each transformed coefficient is raised to the $K$-th power modulo $10^9+7$. After the inverse scaling step, the entry at index $0$ is the number of losing positions. The C++ implementation also includes small brute-force cross-checks and fixed checkpoints to verify the closed form and the convolution result on manageable inputs, while keeping the main large-instance computation entirely in the fast pipeline.

Complexity Analysis

Let $L=N-1$, and let $m$ be the smallest power of two with $m\ge \pi(L)+1$. The linear sieve and the frequency construction take $O(L)$ time. The Walsh-Hadamard transform takes $O(m\log m)$ time, and the pointwise exponentiation takes $O(m\log K)$ time because each coefficient is raised by binary exponentiation. The total memory usage is $O(L+m)$, dominated in practice by the sieve arrays.

Footnotes and References

Problem page: https://projecteuler.net/problem=560
Sprague-Grundy theorem: Wikipedia - Sprague-Grundy theorem
Mex: Wikipedia - Mex (mathematics)
Hadamard transform: Wikipedia - Hadamard transform
Linear sieve: cp-algorithms - Linear Sieve

Problem 560 source code

C++

#include <algorithm>
#include <chrono>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <thread>
#include <vector>

namespace {

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

constexpr u32 kMod = 1'000'000'007U;
constexpr u32 kDefaultN = 10'000'000U;
constexpr u32 kDefaultK = 10'000'000U;

constexpr u32 kCheckpointN1 = 5U;
constexpr u32 kCheckpointK1 = 2U;
constexpr u32 kCheckpointExpected1 = 6U;

constexpr u32 kCheckpointN2 = 10U;
constexpr u32 kCheckpointK2 = 5U;
constexpr u32 kCheckpointExpected2 = 9'964U;

constexpr u32 kCheckpointN3 = 10U;
constexpr u32 kCheckpointK3 = 10U;
constexpr u32 kCheckpointExpected3 = 472'400'303U;

constexpr u32 kCheckpointN4 = 1'000U;
constexpr u32 kCheckpointK4 = 1'000U;
constexpr u32 kCheckpointExpected4 = 954'021'836U;

constexpr u32 kFormulaValidationLimit = 180U;
constexpr u32 kBruteCrossN = 30U;
constexpr u32 kBruteCrossK = 4U;
constexpr u32 kThreadConsistencyN = 100'000U;
constexpr u32 kThreadConsistencyK = 100'000U;

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

struct SieveData {
    u32 limit = 0U;
    std::vector<u32> spf;
    std::vector<u32> pi;
};

inline u32 add_mod(const u32 a, const u32 b) {
    const u32 s = a + b;
    return (s >= kMod) ? (s - kMod) : s;
}

inline u32 sub_mod(const u32 a, const u32 b) {
    return (a >= b) ? (a - b) : (a + kMod - b);
}

inline u32 mul_mod(const u32 a, const u32 b) {
    return static_cast<u32>((static_cast<u64>(a) * static_cast<u64>(b)) % kMod);
}

u32 mod_pow(u32 base, u32 exponent) {
    u32 result = 1U;
    u32 b = base;
    u32 e = exponent;

    while (e > 0U) {
        if ((e & 1U) != 0U) {
            result = mul_mod(result, b);
        }
        b = mul_mod(b, b);
        e >>= 1U;
    }

    return result;
}

bool parse_u32_after_prefix(const std::string& arg, const char* prefix, u32& value) {
    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 u32 digit = static_cast<u32>(c - '0');
        parsed = parsed * 10ULL + static_cast<u64>(digit);
        if (parsed > static_cast<u64>(std::numeric_limits<u32>::max())) {
            return false;
        }
    }

    value = static_cast<u32>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u32 parsed = 0U;
    if (!parse_u32_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    value = static_cast<unsigned>(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 == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u32 parsed_u32 = 0U;
        if (parse_u32_after_prefix(arg, "--n=", parsed_u32)) {
            options.n = parsed_u32;
            continue;
        }
        if (parse_u32_after_prefix(arg, "--k=", parsed_u32)) {
            options.k = 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;
    }
    if (options.k == 0U) {
        std::cerr << "--k must be at least 1.\n";
        return false;
    }

    return true;
}

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

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

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}

SieveData build_linear_sieve(const u32 limit) {
    SieveData data;
    data.limit = limit;
    data.spf.assign(static_cast<std::size_t>(limit) + 1ULL, 0U);
    data.pi.assign(static_cast<std::size_t>(limit) + 1ULL, 0U);

    if (limit < 2U) {
        return data;
    }

    std::vector<u32> primes;
    primes.reserve(static_cast<std::size_t>(limit / 10U) + 32ULL);

    for (u32 i = 2U; i <= limit; ++i) {
        if (data.spf[i] == 0U) {
            data.spf[i] = i;
            primes.push_back(i);
        }

        for (const u32 p : primes) {
            const u64 v = static_cast<u64>(i) * static_cast<u64>(p);
            if (v > static_cast<u64>(limit)) {
                break;
            }
            data.spf[static_cast<std::size_t>(v)] = p;
            if (p == data.spf[i]) {
                break;
            }
        }
    }

    u32 prime_count = 0U;
    for (u32 x = 1U; x <= limit; ++x) {
        if (x >= 2U && data.spf[x] == x) {
            ++prime_count;
        }
        data.pi[x] = prime_count;
    }

    return data;
}

u32 grundy_formula(const u32 x, const SieveData& sieve) {
    if (x == 1U) {
        return 1U;
    }
    if ((x & 1U) == 0U) {
        return 0U;
    }
    return sieve.pi[sieve.spf[x]];
}

std::vector<u32> build_grundy_frequencies(const u32 n) {
    const u32 limit = n - 1U;
    const SieveData sieve = build_linear_sieve(limit);

    const u32 max_index = std::max(1U, sieve.pi[limit]);
    std::vector<u32> counts(static_cast<std::size_t>(max_index) + 1ULL, 0U);

    counts[0] = limit / 2U;
    counts[1] = 1U;

    for (u32 x = 3U; x <= limit; x += 2U) {
        const u32 p = sieve.spf[x];
        const u32 idx = sieve.pi[p];
        ++counts[idx];
    }

    return counts;
}

void fwht_xor(std::vector<u32>& a) {
    const std::size_t n = a.size();

    for (std::size_t len = 1ULL; len < n; len <<= 1ULL) {
        const std::size_t step = len << 1ULL;
        for (std::size_t i = 0ULL; i < n; i += step) {
            for (std::size_t j = 0ULL; j < len; ++j) {
                const u32 u = a[i + j];
                const u32 v = a[i + j + len];
                a[i + j] = add_mod(u, v);
                a[i + j + len] = sub_mod(u, v);
            }
        }
    }
}

void pow_coefficients(std::vector<u32>& coeffs,
                      const u32 exponent,
                      const bool allow_multithreading,
                      const unsigned requested_threads) {
    const std::size_t n = coeffs.size();
    const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, n);

    if (threads == 1U) {
        for (u32& v : coeffs) {
            v = mod_pow(v, exponent);
        }
        return;
    }

    std::vector<std::thread> workers;
    workers.reserve(threads);

    for (unsigned t = 0U; t < threads; ++t) {
        const std::size_t begin = (n * static_cast<std::size_t>(t)) / static_cast<std::size_t>(threads);
        const std::size_t end = (n * static_cast<std::size_t>(t + 1U)) /
                                static_cast<std::size_t>(threads);

        workers.emplace_back([&coeffs, exponent, begin, end]() {
            for (std::size_t i = begin; i < end; ++i) {
                coeffs[i] = mod_pow(coeffs[i], exponent);
            }
        });
    }

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

u32 solve_losing_positions(const u32 n,
                          const u32 k,
                          const bool allow_multithreading,
                          const unsigned requested_threads) {
    const std::vector<u32> counts = build_grundy_frequencies(n);

    std::size_t size = 1ULL;
    while (size < counts.size()) {
        size <<= 1ULL;
    }

    std::vector<u32> transformed(size, 0U);
    for (std::size_t i = 0ULL; i < counts.size(); ++i) {
        transformed[i] = counts[i] % kMod;
    }

    fwht_xor(transformed);
    pow_coefficients(transformed, k, allow_multithreading, requested_threads);
    fwht_xor(transformed);

    const u32 inv_size = mod_pow(static_cast<u32>(size % static_cast<std::size_t>(kMod)), kMod - 2U);
    return mul_mod(transformed[0], inv_size);
}

std::vector<u32> brute_grundy_values(const u32 limit) {
    std::vector<u32> grundy(static_cast<std::size_t>(limit) + 1ULL, 0U);

    for (u32 n = 1U; n <= limit; ++n) {
        std::vector<unsigned char> seen(static_cast<std::size_t>(n) + 8ULL, 0U);

        for (u32 take = 1U; take <= n; ++take) {
            if (std::gcd(take, n) != 1U) {
                continue;
            }
            const u32 next_value = n - take;
            const u32 g = grundy[next_value];
            if (g >= seen.size()) {
                seen.resize(static_cast<std::size_t>(g) + 1ULL, 0U);
            }
            seen[g] = 1U;
        }

        u32 mex = 0U;
        while (mex < seen.size() && seen[mex] != 0U) {
            ++mex;
        }
        grundy[n] = mex;
    }

    return grundy;
}

u32 brute_losing_positions_small(const u32 n, const u32 k) {
    const u32 limit = n - 1U;
    const std::vector<u32> grundy = brute_grundy_values(limit);

    u32 max_grundy = 0U;
    for (u32 x = 1U; x <= limit; ++x) {
        max_grundy = std::max(max_grundy, grundy[x]);
    }

    std::size_t size = 1ULL;
    while (size <= static_cast<std::size_t>(max_grundy)) {
        size <<= 1ULL;
    }

    std::vector<u32> counts(size, 0U);
    for (u32 x = 1U; x <= limit; ++x) {
        ++counts[grundy[x]];
    }

    std::vector<u32> dp(size, 0U);
    std::vector<u32> next(size, 0U);
    dp[0] = 1U;

    for (u32 pile = 0U; pile < k; ++pile) {
        std::fill(next.begin(), next.end(), 0U);
        for (std::size_t xor_value = 0ULL; xor_value < size; ++xor_value) {
            const u32 ways = dp[xor_value];
            if (ways == 0U) {
                continue;
            }
            for (std::size_t g = 0ULL; g < size; ++g) {
                const u32 freq = counts[g];
                if (freq == 0U) {
                    continue;
                }
                const std::size_t target = xor_value ^ g;
                const u32 contribution = mul_mod(ways, freq % kMod);
                next[target] = add_mod(next[target], contribution);
            }
        }
        dp.swap(next);
    }

    return dp[0];
}

bool expect_equal(const char* label, const u32 actual, const u32 expected) {
    if (actual == expected) {
        return true;
    }

    std::cerr << "Checkpoint failed for " << label << ": got " << actual
              << ", expected " << expected << '\n';
    return false;
}

bool run_checkpoints() {
    bool ok = true;

    {
        const SieveData sieve = build_linear_sieve(kFormulaValidationLimit);
        const std::vector<u32> brute = brute_grundy_values(kFormulaValidationLimit);

        for (u32 x = 1U; x <= kFormulaValidationLimit; ++x) {
            const u32 predicted = grundy_formula(x, sieve);
            if (predicted != brute[x]) {
                std::cerr << "Grundy formula mismatch at n=" << x << ": got " << predicted
                          << ", expected " << brute[x] << '\n';
                ok = false;
                break;
            }
        }
    }

    ok = ok && expect_equal("L(5,2)",
                            solve_losing_positions(kCheckpointN1, kCheckpointK1, false, 1U),
                            kCheckpointExpected1);
    ok = ok && expect_equal("L(10,5)",
                            solve_losing_positions(kCheckpointN2, kCheckpointK2, false, 1U),
                            kCheckpointExpected2);
    ok = ok && expect_equal("L(10,10)",
                            solve_losing_positions(kCheckpointN3, kCheckpointK3, false, 1U),
                            kCheckpointExpected3);
    ok = ok && expect_equal("L(10^3,10^3)",
                            solve_losing_positions(kCheckpointN4, kCheckpointK4, false, 1U),
                            kCheckpointExpected4);

    const u32 brute_cross = brute_losing_positions_small(kBruteCrossN, kBruteCrossK);
    const u32 fast_cross = solve_losing_positions(kBruteCrossN, kBruteCrossK, false, 1U);
    ok = ok && expect_equal("brute-vs-fast", fast_cross, brute_cross);

    unsigned hw_threads = std::thread::hardware_concurrency();
    if (hw_threads == 0U) {
        hw_threads = 1U;
    }
    if (hw_threads >= 2U) {
        const u32 single =
            solve_losing_positions(kThreadConsistencyN, kThreadConsistencyK, false, 1U);
        const u32 threaded =
            solve_losing_positions(kThreadConsistencyN, kThreadConsistencyK, true, 2U);
        ok = ok && expect_equal("thread-consistency", threaded, single);
    }

    if (!ok) {
        std::cerr << "At least one checkpoint failed.\n";
    }
    return ok;
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

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

    const auto start = std::chrono::steady_clock::now();
    const u32 answer = solve_losing_positions(
        options.n, options.k, options.allow_multithreading, options.requested_threads);
    const auto end = std::chrono::steady_clock::now();

    const std::chrono::duration<long double> elapsed = end - start;
    std::cout << "L(" << options.n << ", " << options.k << ") mod " << kMod << " = "
              << answer << '\n';
    std::cout << "Elapsed: " << elapsed.count() << " seconds\n";

    return 0;
}

Python

kMod = 1000000007

def build_linear_sieve(limit):
    spf = [0] * (limit + 1)
    pi = [0] * (limit + 1)
    if limit < 2: return spf, pi
    
    primes = []
    
    for i in range(2, limit + 1):
        if spf[i] == 0:
            spf[i] = i
            primes.append(i)
        
        for p in primes:
            v = i * p
            if v > limit: break
            spf[v] = p
            if p == spf[i]: break
            
    prime_count = 0
    for x in range(1, limit + 1):
        if x >= 2 and spf[x] == x:
            prime_count += 1
        pi[x] = prime_count
        
    return spf, pi

def fwht_xor(a):
    n = len(a)
    length = 1
    while length < n:
        step = length << 1
        for i in range(0, n, step):
            for j in range(length):
                u = a[i + j]
                v = a[i + j + length]
                a[i + j] = (u + v) % kMod
                a[i + j + length] = (u - v + kMod) % kMod
        length <<= 1

def solve_losing_positions(n, k):
    limit = n - 1
    spf, pi = build_linear_sieve(limit)
    
    max_index = max(1, pi[limit])
    counts = [0] * (max_index + 1)
    
    counts[0] = limit // 2
    counts[1] = 1
    
    for x in range(3, limit + 1, 2):
        p = spf[x]
        idx = pi[p]
        counts[idx] += 1
        
    size = 1
    while size < len(counts):
        size <<= 1
        
    transformed = [0] * size
    for i in range(len(counts)):
        transformed[i] = counts[i] % kMod
        
    fwht_xor(transformed)
    
    for i in range(size):
        transformed[i] = pow(transformed[i], k, kMod)
        
    fwht_xor(transformed)
    
    inv_size = pow(size % kMod, kMod - 2, kMod)
    return (transformed[0] * inv_size) % kMod

def solve():
    return str(solve_losing_positions(10000000, 10000000))

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

Java

public class Euler560 {
    static final int kMod = 1000000007;

    static int modPow(int base, int exp) {
        long result = 1;
        long b = base;
        int e = exp;
        while (e > 0) {
            if ((e & 1) != 0)
                result = (result * b) % kMod;
            b = (b * b) % kMod;
            e >>= 1;
        }
        return (int) result;
    }

    static class SieveData {
        int[] spf;
        int[] pi;

        SieveData(int limit) {
            spf = new int[limit + 1];
            pi = new int[limit + 1];
            if (limit < 2)
                return;

            int[] primes = new int[limit / 10 + 32];
            int numPrimes = 0;

            for (int i = 2; i <= limit; i++) {
                if (spf[i] == 0) {
                    spf[i] = i;
                    if (numPrimes < primes.length) {
                        primes[numPrimes++] = i;
                    } else {
                        int[] newPrimes = new int[primes.length * 2];
                        System.arraycopy(primes, 0, newPrimes, 0, numPrimes);
                        primes = newPrimes;
                        primes[numPrimes++] = i;
                    }
                }

                for (int j = 0; j < numPrimes; j++) {
                    int p = primes[j];
                    long v = (long) i * p;
                    if (v > limit)
                        break;
                    spf[(int) v] = p;
                    if (p == spf[i])
                        break;
                }
            }

            int primeCount = 0;
            for (int x = 1; x <= limit; x++) {
                if (x >= 2 && spf[x] == x) {
                    primeCount++;
                }
                pi[x] = primeCount;
            }
        }
    }

    static void fwhtXor(int[] a) {
        int n = a.length;
        for (int len = 1; len < n; len <<= 1) {
            int step = len << 1;
            for (int i = 0; i < n; i += step) {
                for (int j = 0; j < len; j++) {
                    int u = a[i + j];
                    int v = a[i + j + len];
                    a[i + j] = (u + v >= kMod) ? (u + v - kMod) : (u + v);
                    a[i + j + len] = (u >= v) ? (u - v) : (u + kMod - v);
                }
            }
        }
    }

    static int solveLosingPositions(int n, int k) {
        int limit = n - 1;
        SieveData sieve = new SieveData(limit);

        int maxIndex = Math.max(1, sieve.pi[limit]);
        int[] counts = new int[maxIndex + 1];

        counts[0] = limit / 2;
        counts[1] = 1;

        for (int x = 3; x <= limit; x += 2) {
            int p = sieve.spf[x];
            int idx = sieve.pi[p];
            counts[idx]++;
        }

        int size = 1;
        while (size < counts.length) {
            size <<= 1;
        }

        int[] transformed = new int[size];
        for (int i = 0; i < counts.length; i++) {
            transformed[i] = counts[i] % kMod;
        }

        fwhtXor(transformed);

        for (int i = 0; i < size; i++) {
            transformed[i] = modPow(transformed[i], k);
        }

        fwhtXor(transformed);

        int invSize = modPow(size % kMod, kMod - 2);
        long answer = ((long) transformed[0] * invSize) % kMod;
        return (int) answer;
    }

    public static String solve() {
        return Integer.toString(solveLosingPositions(10000000, 10000000));
    }

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