You're correct that this is a nim variant. More generally, all impartial games (those where both players have the same moves available to them) of this sort fall under the broad category of the Sprage-Grundy theorem which says that their values are equivalent to the value of some nim-pile.
But for determining who's to win from a given 'single-position' configuration, you don't need any of the broader information; instead, you just need the concepts of 'next player wins' ($\mathcal{N}$-positions) and 'previous player has won' ($\mathcal{P}$-positions), along with the rules for determining the value of a position from all of its options: if any of the moves available from a given position is to a $\mathcal{P}$-position, then that position is an $\mathcal{N}$-position (since the next player can win by moving to the $\mathcal{P}$-position). Otherwise (i.e., if all of the positions that can be reached from the given position are $\mathcal{N}$-positions), the position is a $\mathcal{P}$-position; there are no good moves available.
This evaluation can then be used in one of two ways; if you're interested in the value of a single position then you can simply build the tree (more accurately, DAG — Directed Acyclic Graph) of moves from that position (being careful to cache 'transpositions' - e.g., since making move $m_1$ followed by $m_2$ leads to the same position $P_{12}$ as making move $m_2$ followed by $m_1$, you don't want to compute the value of $P_{12}$ twice) and evaluate it.
If you're going to be computing the value of many positions then it will generally be better to build a table of values, but in games which have many piles this is likely to be a very large table. Since all of the moves have to decrease one or more coordinates, the table can be built outwards from the origin in 'shells': first fill in the value for the origin (which by definition will be a $\mathcal{P}$-position, since there's no legal move to make). Next, fill in the values for all cells with the sum of their pile-sizes equal to 1 (i.e., $(1,0,0,\ldots),\ (0, 1, 0, \ldots),\ (0, 0, 1,\ldots),\ \ldots$); since any move decreases one or more coordinates it has to decrease the sum of the pile-sizes and so you'll already have processed all of the positions a given position's value can depend on before you process that position itself. Note that for a fixed number of dimensions (i.e., piles), it's easy to write the iteration over each shell as a series of nested loops with appropriate limits. For instance, the four-pile case would look like:
for ( totalSum = 0; totalSum <= 16; totalSum++ ) {
for ( pile0 = 0; pile0 <= totalSum; pile0++ ) {
for ( pile1 = 0; pile1 <= totalSum-pile0; pile1++ ) {
for ( pile2 = 0; pile2 <= totalSum-(pile0+pile1); pile2++ ) {
pile3 = totalSum-(pile0+pile1+pile2);
...
}
}
}
}
If the number of piles is given as part of the input you'll have a bit more work to do, but there are still relatively straightforward ways of iterating over all combinations summing to a given number.