Problem 806: Nim on Towers of Hanoi

Project Euler Problem 806 Solution

EulerSolve provides an optimized solution for Project Euler Problem 806, Nim on Towers of Hanoi, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary In the standard optimal Tower of Hanoi transfer with $n$ disks, after each move we record the peg populations $(a,b,c)$. The position is losing in normal Nim exactly when $$a\oplus b\oplus c=0.$$ The required quantity is the sum of all move numbers $t\in\{1,\dots,2^n-1\}$ for which the Hanoi position is Nim-losing. The C++, Python, and Java implementations compute this value modulo $10^9+7$ without simulating the full exponential move sequence. Mathematical Approach The key reduction is to count how many losing times correspond to each admissible ordered population triple, and only at the very end convert that count into a sum of move indices. Step 1: Turn the index sum into a counting problem The optimal Hanoi path has length $2^n-1$. Reversing time sends move $t$ to move $2^n-1-t$, and this reversal only permutes peg roles, so the Nim condition is preserved. Therefore losing times come in complementary pairs whose indices add to $2^n-1$. If $\Lambda(n)$ denotes the total number of losing times, then $$S(n)=\frac{2^n-1}{2}\,\Lambda(n)\pmod{10^9+7}.$$ So the real task is to compute $\Lambda(n)$, the number of Nim-losing moments along the Hanoi path....

Detailed mathematical approach

Problem Summary

In the standard optimal Tower of Hanoi transfer with $n$ disks, after each move we record the peg populations $(a,b,c)$. The position is losing in normal Nim exactly when

$$a\oplus b\oplus c=0.$$

The required quantity is the sum of all move numbers $t\in\{1,\dots,2^n-1\}$ for which the Hanoi position is Nim-losing. The C++, Python, and Java implementations compute this value modulo $10^9+7$ without simulating the full exponential move sequence.

Mathematical Approach

The key reduction is to count how many losing times correspond to each admissible ordered population triple, and only at the very end convert that count into a sum of move indices.

Step 1: Turn the index sum into a counting problem

The optimal Hanoi path has length $2^n-1$. Reversing time sends move $t$ to move $2^n-1-t$, and this reversal only permutes peg roles, so the Nim condition is preserved. Therefore losing times come in complementary pairs whose indices add to $2^n-1$.

If $\Lambda(n)$ denotes the total number of losing times, then

$$S(n)=\frac{2^n-1}{2}\,\Lambda(n)\pmod{10^9+7}.$$

So the real task is to compute $\Lambda(n)$, the number of Nim-losing moments along the Hanoi path.

Step 2: Describe all admissible losing triples

Any Hanoi state satisfies

$$a+b+c=n.$$

For a losing Nim state we also need

$$a\oplus b\oplus c=0.$$

Bit by bit, xor zero means that at every binary position there are either zero 1s or exactly two 1s among $a,b,c$. So the contribution of bit $2^r$ to the sum $a+b+c$ is either $0$ or $2^{r+1}$. This immediately shows that $n$ must be even; if $n$ is odd, the answer is $0$.

Write $n=2m$, and let the binary expansion of $m$ be

$$m=\sum_{r\in R}2^r.$$

For each $r\in R$, exactly one peg omits the bit $2^r$, so the local contribution is one of

$$\left(0,2^r,2^r\right),\qquad \left(2^r,0,2^r\right),\qquad \left(2^r,2^r,0\right).$$

Adding these choices independently over all set bits of $m$ generates every ordered triple $(a,b,c)$ with $a+b+c=n$ and $a\oplus b\oplus c=0$, and it generates each such triple exactly once. Hence the number of admissible losing triples is

$$3^{s_2(m)}=3^{s_2(n/2)},$$

where $s_2(\cdot)$ is the number of 1s in the binary expansion.

Step 3: Introduce the auxiliary generating function

The implementation evaluates an auxiliary coefficient extracted from

$$\Psi(X,Y,Z)=\frac{1}{1-X^2-Y^2-Z^2-2XYZ}.$$

Define

$$\mathcal{A}(\alpha,\beta,\gamma)=\left[X^\alpha Y^\beta Z^\gamma\right]\Psi(X,Y,Z).$$

Using the geometric-series expansion,

$$\Psi(X,Y,Z)=\sum_{i\ge 0}\left(X^2+Y^2+Z^2+2XYZ\right)^i.$$

So in the $i$-th layer we choose some copies of $X^2$, some of $Y^2$, some of $Z^2$, and some of $2XYZ$. That multinomial expansion is exactly the algebra encoded by the implementations.

Step 4: Closed form for the auxiliary coefficient

Fix nonnegative $\alpha,\beta,\gamma$ and set

$$m=\alpha+\beta+\gamma.$$

In one term of degree $i$, suppose $k$ factors contribute $2XYZ$. Then

$$k=m-2i,$$

and the remaining exponents must come from squares, so

$$r_X=\frac{\alpha-k}{2},\qquad r_Y=\frac{\beta-k}{2},\qquad r_Z=\frac{\gamma-k}{2}.$$

These must be nonnegative integers. Equivalently, with

$$A=\frac{\beta+\gamma}{2},\qquad B=\frac{\alpha+\gamma}{2},\qquad C=\frac{\alpha+\beta}{2},$$

we have

$$r_X=i-A,\qquad r_Y=i-B,\qquad r_Z=i-C.$$

Therefore the $i$-th contribution is

$$Q_i(\alpha,\beta,\gamma)=\frac{i!\,2^{m-2i}}{(m-2i)!\,(i-A)!\,(i-B)!\,(i-C)!},$$

and the full auxiliary term is

$$\mathcal{A}(\alpha,\beta,\gamma)=\sum_{i=\ell}^{\lfloor m/2\rfloor} Q_i(\alpha,\beta,\gamma),$$

where

$$\ell=\max\left(\left\lceil\frac{m}{3}\right\rceil,A,B,C\right).$$

If the parity conditions fail, meaning $\alpha,\beta,\gamma$ do not all have the same parity as $m$, then $\mathcal{A}(\alpha,\beta,\gamma)=0$.

Step 5: Recover the exact multiplicity of an ordered triple

The actual number of losing occurrences of a specific ordered Hanoi population triple is obtained by multiplying the generating function by a correction polynomial. Define

$$M(a,b,c)=\left[X^aY^bZ^c\right]\left(1+X+Z+XY+YZ-Y^2\right)\Psi(X,Y,Z).$$

Extracting coefficients gives the six-term identity

$$\begin{aligned} M(a,b,c)=&\ \mathcal{A}(a,b,c)+\mathcal{A}(a-1,b,c)+\mathcal{A}(a,b,c-1)\\ &+\mathcal{A}(a-1,b-1,c)+\mathcal{A}(a,b-1,c-1)-\mathcal{A}(a,b-2,c). \end{aligned}$$

This expression is intentionally not symmetric in $a,b,c$. The standard Hanoi transfer treats the three pegs as source, auxiliary, and target in a fixed order, so the middle peg receives a different correction term.

Step 6: Sum over all losing triples

Let $\mathcal{T}_n$ be the set of ordered triples satisfying

$$a+b+c=n,\qquad a\oplus b\oplus c=0.$$

Then the total number of losing times is

$$\Lambda(n)=\sum_{(a,b,c)\in\mathcal{T}_n} M(a,b,c).$$

Combining this with Step 1 yields the implemented formula

$$\boxed{S(n)=\frac{2^n-1}{2}\sum_{\substack{a+b+c=n\\a\oplus b\oplus c=0}} M(a,b,c)\pmod{10^9+7}.}$$

Worked Example: $n=4$

Here $n/2=2$ has one set bit, so there are exactly three admissible losing triples:

$$\left(0,2,2\right),\qquad \left(2,0,2\right),\qquad \left(2,2,0\right).$$

For $(2,2,0)$, only one summation level contributes to the auxiliary coefficient, giving

$$\mathcal{A}(2,2,0)=\frac{2!}{0!\,1!\,1!\,0!}=2.$$

Also $\mathcal{A}(2,0,0)=1$, and the other shifted terms in the six-term identity vanish, so

$$M(2,2,0)=2-1=1.$$

Similarly,

$$M(0,2,2)=1,\qquad M(2,0,2)=2.$$

Therefore

$$\Lambda(4)=1+1+2=4.$$

The required sum of move indices is then

$$S(4)=\frac{2^4-1}{2}\cdot 4=\frac{15}{2}\cdot 4=30,$$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical pipeline. First they reject odd $n$, because Step 2 proves that no losing triple can exist in that case. Next they precompute factorials, inverse factorials, powers of two, and modular inverses up to about $n/2$, so each multinomial term can be evaluated quickly modulo $10^9+7$.

Then the implementation enumerates every admissible triple by scanning the set bits of $n/2$ and, for each such bit, choosing which peg does not receive that contribution. For one ordered triple, the auxiliary coefficient is not recomputed from scratch for every $i$; instead the code starts at the smallest feasible $i$ and advances with the exact ratio

$$\frac{Q_{i+1}}{Q_i}=\frac{(i+1)(m-2i)(m-2i-1)}{4(i+1-A)(i+1-B)(i+1-C)}.$$

That turns the inner sum into a linear sweep with $O(1)$ modular work per step. After the auxiliary values are known, the implementation combines six shifted coefficients to obtain $M(a,b,c)$, sums these values over all admissible triples to get $\Lambda(n)$, and finally multiplies by $(2^n-1)/2$. The C++ version parallelizes the outer triple loop; the Python and Java versions keep the same mathematics in serial form.

Complexity Analysis

If $n$ is odd, the answer is returned immediately. For even $n$, the number of admissible ordered triples is exactly

$$T=3^{s_2(n/2)}.$$

For each triple, the auxiliary summation runs over at most $\lfloor n/2\rfloor+1$ indices $i$, and every transition uses only $O(1)$ modular operations because the tables are precomputed. Therefore the worst-case running time is

$$O\!\left(n\cdot 3^{s_2(n/2)}\right),$$

with memory usage

$$O(n).$$

Since $s_2(n/2)\le \lfloor\log_2 n\rfloor$, this is dramatically smaller than brute-force Hanoi simulation, which would require $\Theta(2^n)$ moves and states.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=806
Tower of Hanoi: Wikipedia — Tower of Hanoi
Nim: Wikipedia — Nim
Gray code: Wikipedia — Gray code
Multinomial theorem: Wikipedia — Multinomial theorem
Generating function: Wikipedia — Generating function

Problem 806 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <future>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <string>
#include <thread>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>

using namespace std;

static const uint32_t MOD = 1000000007u;

static inline uint32_t addmod(uint32_t a, uint32_t b){ uint32_t r = a + b; if(r >= MOD) r -= MOD; return r; }
static inline uint32_t submod(uint32_t a, uint32_t b){ return a >= b ? a - b : a + MOD - b; }
static inline uint32_t fastmod(uint64_t x){
    static const uint64_t inv = (uint64_t)((__uint128_t(1) << 64) / MOD);
    uint64_t q = (uint64_t)((__uint128_t)x * inv >> 64);
    uint64_t r = x - q * MOD;
    if(r >= MOD) r -= MOD;
    return (uint32_t)r;
}
static inline uint32_t mulmod(uint32_t a, uint32_t b){ return fastmod((uint64_t)a * b); }

static uint32_t modpow(uint32_t a, uint32_t e){
    uint32_t r = 1u;
    a %= MOD;
    while(e > 0){
        if(e & 1) r = mulmod(r, a);
        a = mulmod(a, a);
        e >>= 1;
    }
    return r;
}

struct Comb {
    int maxI;
    vector<uint32_t> fact, invfact, pow2, inv;
    explicit Comb(int maxI_=0){ init(maxI_); }
    void init(int maxI_){
        maxI = max(0, maxI_);
        fact.assign(maxI+1, 1);
        invfact.assign(maxI+1, 1);
        pow2.assign(maxI+1, 1);
        inv.assign(maxI+1, 0);
        for(int i=1;i<=maxI;i++){
            fact[i] = mulmod(fact[i-1], (uint32_t)i);
            pow2[i] = addmod(pow2[i-1], pow2[i-1]);
        }
        invfact[maxI] = modpow(fact[maxI], MOD-2);
        for(int i=maxI;i>=1;i--){
            invfact[i-1] = mulmod(invfact[i], (uint32_t)i);
        }
        if(maxI >= 1){
            inv[1] = 1;
            for(int i=2;i<=maxI;i++){
                inv[i] = MOD - mulmod((uint32_t)(MOD / i), inv[MOD % i]);
            }
        }
    }
};

static inline uint32_t coefP(const Comb& C, int i, int u, int v, int w){
    int n = u + v + w;
    int k = n - 2*i;
    if(k < 0 || k > i) return 0;
    if(u < k || v < k || w < k) return 0;
    if(((u-k)&1) || ((v-k)&1) || ((w-k)&1)) return 0;
    int p = (u-k)/2, q = (v-k)/2, r = (w-k)/2;
    if(p < 0 || q < 0 || r < 0) return 0;
    uint32_t res = C.fact[i];
    res = mulmod(res, C.invfact[k]);
    res = mulmod(res, C.invfact[p]);
    res = mulmod(res, C.invfact[q]);
    res = mulmod(res, C.invfact[r]);
    res = mulmod(res, C.pow2[k]);
    return res;
}

static uint32_t H(const Comb& C, int u, int v, int w){
    if(u < 0 || v < 0 || w < 0) return 0;
    int n = u + v + w;
    if(((u ^ n) & 1) || ((v ^ n) & 1) || ((w ^ n) & 1)) return 0;
    int A = (v + w) / 2;
    int B = (u + w) / 2;
    int Cc = (u + v) / 2;
    int lo = (n + 2) / 3;
    if(A > lo) lo = A;
    if(B > lo) lo = B;
    if(Cc > lo) lo = Cc;
    int hi = n / 2;
    if(lo > hi) return 0;
    int i0 = lo;
    int k0 = n - 2 * i0;
    int p0 = i0 - A;
    int q0 = i0 - B;
    int r0 = i0 - Cc;
    uint32_t coef = C.fact[i0];
    coef = mulmod(coef, C.invfact[k0]);
    coef = mulmod(coef, C.invfact[p0]);
    coef = mulmod(coef, C.invfact[q0]);
    coef = mulmod(coef, C.invfact[r0]);
    coef = mulmod(coef, C.pow2[k0]);
    uint32_t s = coef;
    const uint32_t inv4 = C.inv[4];
    int k = k0;
    for(int i = i0; i < hi; i++){
        uint32_t num = mulmod((uint32_t)(i + 1), mulmod((uint32_t)k, (uint32_t)(k - 1)));
        uint32_t den = mulmod(C.inv[i + 1 - A], mulmod(C.inv[i + 1 - B], C.inv[i + 1 - Cc]));
        coef = mulmod(coef, num);
        coef = mulmod(coef, den);
        coef = mulmod(coef, inv4);
        s = addmod(s, coef);
        k -= 2;
    }
    return s;
}

static uint32_t gCoeff(const Comb& C, int a, int b, int c){
    uint32_t res = 0;
    res = addmod(res, H(C, a,   b,   c));
    res = addmod(res, H(C, a-1, b,   c));
    res = addmod(res, H(C, a,   b,   c-1));
    res = addmod(res, H(C, a-1, b-1, c));
    res = addmod(res, H(C, a,   b-1, c-1));
    res = submod(res, H(C, a,   b-2, c));
    return res;
}

static vector<array<int,3>> enumerateXorTriples(int n){
    vector<array<int,3>> cur;
    cur.push_back({0,0,0});
    if(n & 1) return {};
    for(int bit = 1; (1<<bit) <= n*2 || bit <= 20; bit++){
        if((n >> bit) & 1){
            int v = 1 << (bit-1);
            vector<array<int,3>> nxt;
            nxt.reserve(cur.size() * 3);
            for(auto t: cur){
                nxt.push_back({t[0], t[1] + v, t[2] + v});
                nxt.push_back({t[0] + v, t[1], t[2] + v});
                nxt.push_back({t[0] + v, t[1] + v, t[2]});
            }
            cur.swap(nxt);
        }
        if((1<<bit) > n && bit > 20) break;
    }
    vector<array<int,3>> out;
    out.reserve(cur.size());
    for(auto t: cur){
        if(t[0] + t[1] + t[2] == n && ((t[0] ^ t[1] ^ t[2]) == 0)) out.push_back(t);
    }
    return out;
}

static uint32_t solveFast(const Comb& C, int n){
    if(n < 0) return 0;
    if(n & 1) return 0;
    auto triples = enumerateXorTriples(n);

    unsigned hw = thread::hardware_concurrency();
    unsigned T = max(1u, hw);
    T = min<unsigned>(T, (unsigned)triples.size());
    if(T == 0) T = 1;

    vector<uint32_t> partial(T, 0);
    vector<thread> threads;
    threads.reserve(T);
    atomic<size_t> done(0);
    atomic<bool> running(true);
    const size_t total = triples.size();
    thread progress;
    if(total > 0){
        progress = thread([&](){
            const int width = 40;
            auto start = chrono::steady_clock::now();
            while(running.load(memory_order_relaxed)){
                size_t d = done.load(memory_order_relaxed);
                double ratio = total ? (double)d / (double)total : 1.0;
                if(ratio > 1.0) ratio = 1.0;
                int filled = (int)(ratio * width);
                auto now = chrono::steady_clock::now();
                double secs = chrono::duration<double>(now - start).count();
                cerr << "\r[";
                for(int i=0;i<width;i++) cerr << (i < filled ? '#' : '-');
                cerr << "] " << fixed << setprecision(2) << (ratio * 100.0) << "% "
                     << d << "/" << total << " " << (uint64_t)secs << "s" << flush;
                this_thread::sleep_for(chrono::milliseconds(200));
            }
        });
    }

    auto worker = [&](unsigned tid){
        size_t L = triples.size();
        size_t from = (L * tid) / T;
        size_t to   = (L * (tid + 1)) / T;
        uint32_t acc = 0;
        for(size_t i = from; i < to; i++){
            auto t = triples[i];
            acc = addmod(acc, gCoeff(C, t[0], t[1], t[2]));
            done.fetch_add(1, memory_order_relaxed);
        }
        partial[tid] = acc;
    };
    for(unsigned t=0; t<T; t++) threads.emplace_back(worker, t);
    for(auto& th: threads) th.join();
    running.store(false, memory_order_relaxed);
    if(progress.joinable()){
        progress.join();
        cerr << "\r";
        cerr << string(80, ' ') << "\r";
    }

    uint32_t K = 0;
    for(unsigned t=0; t<T; t++) K = addmod(K, partial[t]);

    const uint32_t inv2 = (MOD + 1) / 2;
    uint32_t pow2n = modpow(2u, (uint32_t)n);
    uint32_t ans = mulmod(K, submod(pow2n, 1));
    ans = mulmod(ans, inv2);
    return ans;
}

static uint32_t solveBrute(int n){
    if(n < 0) return 0;
    if(n == 0) return 0;
    int N = n;

    vector<uint8_t> loc(N+1, 0);
    int cnt0 = N, cnt1 = 0, cnt2 = 0;
    vector<int8_t> dir(N+1, 0);
    for(int d=1; d<=N; d++){
        dir[d] = (((N-d)&1) ? 1 : -1);
    }

    uint64_t sumIdx = 0;
    auto isLose = [&](int a,int b,int c){ return ((a ^ b ^ c) == 0); };
    if(isLose(cnt0,cnt1,cnt2)) sumIdx = 0;
    uint32_t totalStates = (N >= 31 ? 0u : (1u << N));

    for(uint32_t m = 1; m < totalStates; m++){
        int d = __builtin_ctz(m) + 1;
        int oldp = loc[d];
        int newp = oldp + dir[d];
        if(newp < 0) newp += 3;
        else if(newp >= 3) newp -= 3;
        loc[d] = (uint8_t)newp;
        if(oldp == 0) cnt0--;
        else if(oldp == 1) cnt1--;
        else cnt2--;
        if(newp == 0) cnt0++;
        else if(newp == 1) cnt1++;
        else cnt2++;

        if(isLose(cnt0,cnt1,cnt2)){
            sumIdx += (uint64_t)m;
            sumIdx %= MOD;
        }
    }
    return (uint32_t)sumIdx;
}

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

    int n = 100000;
    if(argc > 1){
        n = atoi(argv[1]);
    }

    int maxN = max(n, 10);
    int maxI = maxN / 2 + 5;
    Comb C(maxI);

    {
        uint32_t f4_fast = solveFast(C, 4);
        uint32_t f10_fast = solveFast(C, 10);
        if(f4_fast != 30 || f10_fast != 67518){
            cerr << "Validation failed (fast): f(4)=" << f4_fast
                 << ", f(10)=" << f10_fast << "\n";
            return 1;
        }
        uint32_t f4_brute = solveBrute(4);
        uint32_t f10_brute = solveBrute(10);
        if(f4_brute != 30 || f10_brute != 67518 || f4_brute != f4_fast || f10_brute != f10_fast){
            cerr << "Validation failed (brute): f(4)=" << f4_brute
                 << ", f(10)=" << f10_brute << "\n";
            return 1;
        }
    }

    cout << solveFast(C, n) << "\n";
    return 0;
}

Python

import sys

sys.set_int_max_str_digits(200000)
MOD = 1000000007

class Comb:
    def __init__(self, maxI):
        self.maxI = max(0, maxI)
        self.fact = [1] * (self.maxI + 1)
        self.invfact = [1] * (self.maxI + 1)
        self.pow2 = [1] * (self.maxI + 1)
        self.inv = [0] * (self.maxI + 1)
        
        for i in range(1, self.maxI + 1):
            self.fact[i] = (self.fact[i - 1] * i) % MOD
            self.pow2[i] = (self.pow2[i - 1] * 2) % MOD
            
        self.invfact[self.maxI] = pow(self.fact[self.maxI], MOD - 2, MOD)
        for i in range(self.maxI, 0, -1):
            self.invfact[i - 1] = (self.invfact[i] * i) % MOD
            
        if self.maxI >= 1:
            self.inv[1] = 1
            for i in range(2, self.maxI + 1):
                self.inv[i] = (MOD - (MOD // i) * self.inv[MOD % i] % MOD) % MOD

def H(C, u, v, w):
    if u < 0 or v < 0 or w < 0:
        return 0
    n = u + v + w
    if ((u ^ n) & 1) or ((v ^ n) & 1) or ((w ^ n) & 1):
        return 0
        
    A = (v + w) // 2
    B = (u + w) // 2
    Cc = (u + v) // 2
    
    lo = (n + 2) // 3
    if A > lo: lo = A
    if B > lo: lo = B
    if Cc > lo: lo = Cc
    
    hi = n // 2
    if lo > hi:
        return 0
        
    i0 = lo
    k0 = n - 2 * i0
    p0 = i0 - A
    q0 = i0 - B
    r0 = i0 - Cc
    
    coef = C.fact[i0]
    coef = (coef * C.invfact[k0]) % MOD
    coef = (coef * C.invfact[p0]) % MOD
    coef = (coef * C.invfact[q0]) % MOD
    coef = (coef * C.invfact[r0]) % MOD
    coef = (coef * C.pow2[k0]) % MOD
    
    s = coef
    inv4 = C.inv[4]
    k = k0
    
    for i in range(i0, hi):
        num = ((i + 1) * k * (k - 1)) % MOD
        den = (C.inv[i + 1 - A] * C.inv[i + 1 - B] % MOD) * C.inv[i + 1 - Cc] % MOD
        coef = (coef * num) % MOD
        coef = (coef * den) % MOD
        coef = (coef * inv4) % MOD
        s = (s + coef) % MOD
        k -= 2
        
    return s

def gCoeff(C, a, b, c):
    res = 0
    res = (res + H(C, a, b, c)) % MOD
    res = (res + H(C, a - 1, b, c)) % MOD
    res = (res + H(C, a, b, c - 1)) % MOD
    res = (res + H(C, a - 1, b - 1, c)) % MOD
    res = (res + H(C, a, b - 1, c - 1)) % MOD
    res = (res - H(C, a, b - 2, c) + MOD) % MOD
    return res

def enumerate_xor_triples(n):
    cur = [(0, 0, 0)]
    if n % 2 == 1:
        return []
        
    bit = 1
    while (1 << bit) <= n * 2 or bit <= 20:
        if (n >> bit) & 1:
            v = 1 << (bit - 1)
            nxt = []
            for t in cur:
                nxt.append((t[0], t[1] + v, t[2] + v))
                nxt.append((t[0] + v, t[1], t[2] + v))
                nxt.append((t[0] + v, t[1] + v, t[2]))
            cur = nxt
        if (1 << bit) > n and bit > 20:
            break
        bit += 1
        
    out = []
    for t in cur:
        if t[0] + t[1] + t[2] == n and (t[0] ^ t[1] ^ t[2]) == 0:
            out.append(t)
    return out

def solve_fast(C, n):
    if n < 0 or n % 2 == 1:
        return 0
        
    triples = enumerate_xor_triples(n)
    K = 0
    for t in triples:
        K = (K + gCoeff(C, t[0], t[1], t[2])) % MOD
        
    inv2 = (MOD + 1) // 2
    pow2n = pow(2, n, MOD)
    ans = (K * (pow2n - 1 + MOD)) % MOD
    ans = (ans * inv2) % MOD
    return ans

def solve():
    n = 100000
    maxN = max(n, 10)
    maxI = maxN // 2 + 5
    C = Comb(maxI)
    return str(solve_fast(C, n))

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

Java

import java.util.ArrayList;

public class Euler806 {

    static final long MOD = 1000000007L;

    static long addmod(long a, long b) {
        long r = a + b;
        if (r >= MOD)
            r -= MOD;
        return r;
    }

    static long submod(long a, long b) {
        return a >= b ? a - b : a + MOD - b;
    }

    static long mulmod(long a, long b) {
        return (a * b) % MOD;
    }

    static long modpow(long a, long e) {
        long r = 1L;
        a %= MOD;
        while (e > 0) {
            if ((e & 1L) == 1L)
                r = mulmod(r, a);
            a = mulmod(a, a);
            e >>= 1L;
        }
        return r;
    }

    static class Comb {
        int maxI;
        long[] fact, invfact, pow2, inv;

        Comb(int maxI_) {
            maxI = Math.max(0, maxI_);
            fact = new long[maxI + 1];
            invfact = new long[maxI + 1];
            pow2 = new long[maxI + 1];
            inv = new long[maxI + 1];

            fact[0] = 1;
            pow2[0] = 1;
            invfact[0] = 1;

            for (int i = 1; i <= maxI; i++) {
                fact[i] = mulmod(fact[i - 1], i);
                pow2[i] = addmod(pow2[i - 1], pow2[i - 1]);
            }
            invfact[maxI] = modpow(fact[maxI], MOD - 2);
            for (int i = maxI; i >= 1; i--) {
                invfact[i - 1] = mulmod(invfact[i], i);
            }
            if (maxI >= 1) {
                inv[1] = 1;
                for (int i = 2; i <= maxI; i++) {
                    inv[i] = MOD - mulmod(MOD / i, inv[(int) (MOD % i)]);
                }
            }
        }
    }

    static long H(Comb C, int u, int v, int w) {
        if (u < 0 || v < 0 || w < 0)
            return 0;
        int n = u + v + w;
        if (((u ^ n) & 1) != 0 || ((v ^ n) & 1) != 0 || ((w ^ n) & 1) != 0)
            return 0;

        int A = (v + w) / 2;
        int B = (u + w) / 2;
        int Cc = (u + v) / 2;

        int lo = (n + 2) / 3;
        if (A > lo)
            lo = A;
        if (B > lo)
            lo = B;
        if (Cc > lo)
            lo = Cc;

        int hi = n / 2;
        if (lo > hi)
            return 0;

        int i0 = lo;
        int k0 = n - 2 * i0;
        int p0 = i0 - A;
        int q0 = i0 - B;
        int r0 = i0 - Cc;

        long coef = C.fact[i0];
        coef = mulmod(coef, C.invfact[k0]);
        coef = mulmod(coef, C.invfact[p0]);
        coef = mulmod(coef, C.invfact[q0]);
        coef = mulmod(coef, C.invfact[r0]);
        coef = mulmod(coef, C.pow2[k0]);

        long s = coef;
        long inv4 = C.inv[4];
        int k = k0;

        for (int i = i0; i < hi; i++) {
            long num = mulmod(i + 1, mulmod(k, k - 1));
            long den = mulmod(C.inv[i + 1 - A], mulmod(C.inv[i + 1 - B], C.inv[i + 1 - Cc]));
            coef = mulmod(coef, num);
            coef = mulmod(coef, den);
            coef = mulmod(coef, inv4);
            s = addmod(s, coef);
            k -= 2;
        }

        return s;
    }

    static long gCoeff(Comb C, int a, int b, int c) {
        long res = 0;
        res = addmod(res, H(C, a, b, c));
        res = addmod(res, H(C, a - 1, b, c));
        res = addmod(res, H(C, a, b, c - 1));
        res = addmod(res, H(C, a - 1, b - 1, c));
        res = addmod(res, H(C, a, b - 1, c - 1));
        res = submod(res, H(C, a, b - 2, c));
        return res;
    }

    static class Triple {
        int a, b, c;

        Triple(int a, int b, int c) {
            this.a = a;
            this.b = b;
            this.c = c;
        }
    }

    static ArrayList<Triple> enumerateXorTriples(int n) {
        ArrayList<Triple> cur = new ArrayList<>();
        cur.add(new Triple(0, 0, 0));
        if ((n & 1) == 1)
            return new ArrayList<>();

        for (int bit = 1; (1 << bit) <= n * 2 || bit <= 20; bit++) {
            if (((n >> bit) & 1) == 1) {
                int v = 1 << (bit - 1);
                ArrayList<Triple> nxt = new ArrayList<>();
                for (Triple t : cur) {
                    nxt.add(new Triple(t.a, t.b + v, t.c + v));
                    nxt.add(new Triple(t.a + v, t.b, t.c + v));
                    nxt.add(new Triple(t.a + v, t.b + v, t.c));
                }
                cur = nxt;
            }
            if ((1 << bit) > n && bit > 20)
                break;
        }

        ArrayList<Triple> out = new ArrayList<>();
        for (Triple t : cur) {
            if (t.a + t.b + t.c == n && ((t.a ^ t.b ^ t.c) == 0)) {
                out.add(t);
            }
        }
        return out;
    }

    static long solveFast(Comb C, int n) {
        if (n < 0)
            return 0;
        if ((n & 1) == 1)
            return 0;
        ArrayList<Triple> triples = enumerateXorTriples(n);

        long K = 0;
        for (Triple t : triples) {
            K = addmod(K, gCoeff(C, t.a, t.b, t.c));
        }

        long inv2 = (MOD + 1) / 2;
        long pow2n = modpow(2, n);
        long ans = mulmod(K, submod((int) pow2n, 1));
        ans = mulmod(ans, (int) inv2);
        return ans;
    }

    public static String solve() {
        int n = 100000;
        int maxI = Math.max(n, 10) / 2 + 5;
        Comb C = new Comb(maxI);
        return Long.toString(solveFast(C, n));
    }

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