This is a problem on a three-dimensional grid, but we will first consider the simplified version on a one-dimensional grid, then extend it to two-dimensional, and the original three-dimensional version.

[1] One-dimensional version

We consider the following problem.

At first he is in square \(0\) on an infinite one-dimensional grid. When he is in square \(x\), he can move to square \(x-1\) or square \(x+1\). How many ways are there to reach square \(X\) after \(N\) moves?

If \(N < |X|\), then it is impossible to reach square \(X\) after \(N\) moves, so the answer is \(0\).
Also, after moving odd number of times, the resulting position \(x\) is always odd, and after moving even number of times, \(x\) is always even, so if \(N \bmod 2 \neq X \bmod 2\), the answer is \(0\) too.

We consider other cases. If \(X \geq 0\) (if \(X < 0\)), “he is in square \(X\) after exactly \(N\) moves” if and only if “in the \(N\) moves, he moves \((N+|X|)/2\) times in the positive (negative) direction and \((N-|X|)/2\) times in the negative (positive) direction of \(x\) coordinate.” Therefore, the desired number of ways can be expressed with a binomial coefficient as \(\binom{N}{(N+|X|)/2}\).

Therefore, the number of ways \(f_1(N, X)\) to reach square \(X\) after exactly \(N\) moves is

\[ f_1(N, X) = \begin{cases} \displaystyle\binom{N}{(N+|X|)/2} &(\text{if }N \geq |X|\text{ and }N \bmod 2 = X \bmod 2)\\ 0&(\text{otherwise}). \end{cases} \]

By precalculating the factorials, \(f_1(N, X)\) can be calculated for given \(N\) and \(X\) in an \(\mathrm{O}(1)\) time.

[2] Two-dimensional version

Consider the following problem.

At first he is in square \((0, 0)\) on an infinite two-dimensional grid. When he is in square \((x, y)\), he can move to square \((x-1, y), (x+1, y), (x, y-1)\) or square \((x, y+1)\). How many ways are there to reach square \(X\) after \(N\) moves?

In order to solve this problem, we introduce the following coordinate transformation on \(xy\)-coordinate:

\[ \begin{cases} u = x+y\\ v = x-y\\ \end{cases} \]

and we will rephrase the problem as a problem on \(uv\)-coordinate. Also, let \((U, V) := (X+Y, X-Y)\).

Then, “the number of ways to reach square \((X, Y)\) after \(N\) moves” is equal to “the number of ways to reach \((u, v) = (U, V)\) after \(N\) moves.”

Also, the counterpart of “when he is in square \((x, y)\), he can move to square \((x-1, y), (x+1, y), (x, y-1)\) or square \((x, y+1)\)” on \(uv\)-coordinate system is “when he is at \((u, v)\), he can move to \((u-1, v-1), (u-1, v+1), (u+1, v-1)\), or \((u+1, v+1)\). Therefore, on each move, “whether to increment or decrement \(u\)” and “whether to increment or decrement \(v\)” can be determined independently.

Hence, in order to determine a way to reach \((u, v) = (U, V)\) after exactly \(N\) moves, we can choose the following two routes independently:

• the route along the \(u\) axis to reach \(u = U\) after \(N\) moves ( \(f_1(N, U)\) choices) and
• the route along the \(v\) axis to reach \(v = V\) after \(N\) moves ( \(f_1(N, V)\) choices),

so the number of ways \(f_2(N, X, Y)\) to reach square \((X, Y)\) after exactly \(N\) moves can be expressed as

\[ f_2(N, X, Y) = f_1(N, U) f_1(N, V) = f_1(N, X+Y) f_1(N, X-Y). \]

[3] Three-dimensional problem

We now consider the original three-dimensional problem.

We will count the number of ways to reach square \((X, Y, Z)\) such that he moves in the direction of \(z\) axis \(k\) times. Such a way is determined by the following three factors:

• there are \(\binom{N}{k}\) ways of choosing which \(k\) moves of the \(N\) moves are in the direction of \(z\) axis;
• the \(k\) moves along \(z\) axis are chosen so that he reaches \(z = Z\) after \(k\) moves along \(z\) axis, so there are \(f_1(k, Z)\) ways to do so;
• the other \((N-k)\) moves in the direction other than \(x\) axis are chosen so that he reaches \((x, y) = (X, Y)\) after \((N-k)\) moves on the two-dimensional \(xy\)-grid, so there are \(f_2(N-k, X, Y)\) ways to do so.

Therefore, the number of ways to reach square \((X, Y, Z)\) such that he moves in the direction of \(z\) axis \(k\) times is

\[ \binom{N}{k} f_1(k, Z) f_2(N-k, X, Y). \]

By finding the sum over \(k = 0, 1, \ldots, N\), the answer for this problem \(f_3(N, X, Y, Z)\) is obtained as

\[ f_3(N, X, Y, Z) = \sum_{k = 0}^N \binom{N}{k} f_1(k, Z) f_2(N-k, X, Y) = \sum_{k = 0}^N \binom{N}{k} f_1(k, Z) f_1(N-k, X+Y) f_1(N-k, X-Y). \]

With proper precalculations, necessary binomial coefficients can be computed in an \(\mathrm{O}(1)\) time each, so \(f_3(N, X, Y, Z)\) can be computed in a total of \(\mathrm{O}(N)\) time.

Hence, the problem can be solved in a total of \(\mathrm{O}(N)\) time.