Detailed mathematical approach

Problem Summary

For each triple \((a,b,c)\) with \(0 \le a,b,c \le 87\), define \(X=a^3\), \(Y=b^3\), and \(Z=c^3\). The quantity \(N(X,Y,Z)\) counts words made from \(X\) copies of \(i\), \(Y\) copies of \(j\), and \(Z\) copies of \(k\) whose product in the quaternion group \(Q_8\) lands on the required central sign. The final task is

$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$

The outer loop over \(88^3\) cube triples is not the hard part. The real issue is to evaluate \(N(X,Y,Z)\) without enumerating the raw multinomial family of words, whose size is

$\binom{X+Y+Z}{X,Y,Z}=\frac{(X+Y+Z)!}{X!\,Y!\,Z!}.$

The implementations solve that counting problem exactly by turning quaternion multiplication into a statement about inversion parity in multiset permutations.

Mathematical Approach

Write \(S=X+Y+Z\). Every word with those letter counts can be viewed as a permutation of the multiset containing \(X\) letters \(i\), \(Y\) letters \(j\), and \(Z\) letters \(k\). The key is to compare such a word with the canonical ordered word \(i^Xj^Yk^Z\).

Canonical Ordering Of A Quaternion Word

Fix the alphabet order \(i \lt j \lt k\). For a word \(w\), let \(\operatorname{inv}(w)\) be the number of pairs that are out of that order, equivalently the total number of occurrences of \(j\) before \(i\), plus \(k\) before \(i\), plus \(k\) before \(j\).

Each adjacent swap of two distinct generators flips the sign, because

$ji=-ij,\qquad kj=-jk,\qquad ik=-ki.$

Therefore reordering any word into \(i^Xj^Yk^Z\) gives

$P(w)=(-1)^{\operatorname{inv}(w)}\,i^Xj^Yk^Z.$

This is the central invariant behind the whole solution: the word value depends only on the fixed ordered product \(i^Xj^Yk^Z\) and on the inversion parity of the word.

When The Ordered Product Can Be Central

The quaternion group has central elements only \(\pm1\). So we first ask when \(i^Xj^Yk^Z\) is already central.

If the parities of \(X\), \(Y\), and \(Z\) are not all the same, then one of the basis directions survives, and the ordered product is one of \(\pm i\), \(\pm j\), or \(\pm k\). In that case no word can contribute to \(N(X,Y,Z)\), so

$N(X,Y,Z)=0.$

If \(X\), \(Y\), and \(Z\) are all even, say \(X=2x\), \(Y=2y\), \(Z=2z\), then

$i^Xj^Yk^Z=(i^2)^x(j^2)^y(k^2)^z=(-1)^{x+y+z}=(-1)^{S/2}.$

If they are all odd, say \(X=2x+1\), \(Y=2y+1\), \(Z=2z+1\), then

$i^Xj^Yk^Z=(-1)^{x+y+z}ijk=(-1)^{x+y+z+1}=(-1)^{\lfloor S/2\rfloor}.$

So whenever the three parities agree, the ordered product simplifies to the compact formula

$i^Xj^Yk^Z=(-1)^{\lfloor S/2\rfloor}.$

The Admissible Inversion Parity

The condition encoded by the implementations is that the quaternion word must land on the central sign \((-1)^S\). Combining that with the previous identity gives

$(-1)^{\operatorname{inv}(w)}(-1)^{\lfloor S/2\rfloor}=(-1)^S.$

Equivalently,

$\operatorname{inv}(w)\equiv \left\lceil\frac{S}{2}\right\rceil \pmod 2.$

So \(N(X,Y,Z)\) is not an arbitrary group-theoretic count any more. It is exactly the number of multiset permutations with a prescribed inversion parity, provided the parity filter from the previous subsection is satisfied.

Splitting The Multinomial Count Into Even And Odd Inversions

Let

$T=\binom{S}{X,Y,Z},$

and let \(E\) and \(O\) denote the numbers of words with even and odd inversion counts. Then

$E+O=T.$

To separate the two halves, introduce the inversion generating function

$\sum_w q^{\operatorname{inv}(w)}=\binom{S}{X,Y,Z}_q,$

the \(q\)-multinomial over the multiset words \(w\). Evaluating at \(q=-1\) gives

$E-O=\binom{S}{X,Y,Z}_{q=-1}.$

For this problem, the only cases that survive the central-sign filter are especially simple:

$E-O=0 \quad \text{when } X,Y,Z \text{ are all odd},$

and if \(X,Y,Z\) are all even, with

$H=\binom{S/2}{X/2,Y/2,Z/2},$

then

$E-O=H.$

That half-scale multinomial \(H\) is exactly the correction term used by the implementations.

Closed Form For \(N(X,Y,Z)\)

Because the required inversion parity is \(\lceil S/2\rceil \bmod 2\), the final count is

$\boxed{ N(X,Y,Z)= \begin{cases} 0, & \text{if } X,Y,Z \text{ do not have the same parity},\\[6pt] \dfrac{T}{2}, & \text{if } X,Y,Z \text{ are all odd},\\[10pt] \dfrac{T+H}{2}, & \text{if } X,Y,Z \text{ are all even and } S/2 \text{ is even},\\[10pt] \dfrac{T-H}{2}, & \text{if } X,Y,Z \text{ are all even and } S/2 \text{ is odd}. \end{cases}}$

This is precisely the branch structure shared by the C++, Python, and Java implementations.

Worked Checkpoints

The small exact values used in the implementations are good sanity checks for the derivation.

For \((X,Y,Z)=(2,2,2)\), we have \(S=6\),

$T=\binom{6}{2,2,2}=90,\qquad H=\binom{3}{1,1,1}=6.$

Since \(S/2=3\) is odd, the admissible words are the odd-inversion ones, so

$N(2,2,2)=\frac{90-6}{2}=42.$

For \((X,Y,Z)=(8,8,8)\), we have \(S=24\),

$T=\binom{24}{8,8,8}=9465511770,\qquad H=\binom{12}{4,4,4}=34650.$

Now \(S/2=12\) is even, so the admissible words are the even-inversion ones, giving

$N(8,8,8)=\frac{9465511770+34650}{2}=4732773210.$

Those are the exact checkpoints verified before the main modular sweep.

The Final Sum Over Cubes

Once the closed form for \(N(X,Y,Z)\) is known, the Project Euler quantity is just

$A=\sum_{a=0}^{87}\sum_{b=0}^{87}\sum_{c=0}^{87} N(a^3,b^3,c^3)\pmod{888888883}.$

Because cubing preserves parity, a triple contributes only when \(a\), \(b\), and \(c\) are all even or all odd. Every surviving term is then obtained by a constant-time multinomial calculation.

How the Code Works

Modular Combinatorial Tables

The implementations first build the list of cubes \(0^3,1^3,\dots,87^3\) and precompute factorials and inverse factorials up to

$3\cdot 87^3=1975509.$

Since the modulus \(M=888888883\) is prime, every inverse factorial is obtained with Fermat's little theorem:

$a^{-1}\equiv a^{M-2}\pmod M.$

This makes every multinomial evaluation \(O(1)\) after preprocessing:

$\binom{n}{a,b,c}\equiv n!\,(a!)^{-1}(b!)^{-1}(c!)^{-1}\pmod M.$

Evaluating One Triple And Sweeping The Whole Grid

For one triple \((X,Y,Z)\), the implementation first applies the parity filter. If the parities do not match, the contribution is zero immediately. Otherwise it computes the full multinomial \(T\), and in the all-even case it also computes the half-scale multinomial \(H\). Dividing by 2 modulo \(M\) is done by multiplying with \(2^{-1}=(M+1)/2\).

Before the main accumulation, the implementations also confirm the exact integer values \(N(2,2,2)=42\) and \(N(8,8,8)=4732773210\) using wide exact arithmetic. After that validation, they run a straightforward serial triple loop over all \(88^3\) cube triples and add the modular contribution of each one.

Complexity Analysis

The factorial and inverse-factorial tables dominate memory, and they are sized up to \(3\cdot 87^3\). That preprocessing phase costs \(O(87^3)\) time and \(O(87^3)\) memory.

The final sweep visits exactly \(88^3\) triples. Each visit performs only constant-time arithmetic once the tables exist, so that phase is \(O(88^3)\). Overall, the algorithm runs in

$O(87^3+88^3)$

time and uses \(O(87^3)\) additional memory. The exact checkpoint calculations are tiny compared with those two main phases.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=981
2. Quaternion group: Wikipedia - Quaternion group
3. Multinomial theorem: Wikipedia - Multinomial theorem
4. Inversion in combinatorics: Wikipedia - Inversion (discrete mathematics)
5. Gaussian binomial coefficient: Wikipedia - Gaussian binomial coefficient
6. Fermat's little theorem: Wikipedia - Fermat's little theorem