Problem 981: The Quaternion Group II

Project Euler Problem 981 Solution

EulerSolve provides an optimized solution for Project Euler Problem 981, The Quaternion Group II, with C++, Python, Java, and a step-by-step mathematical explanation.

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$....

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

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

Problem 981 source code

C++

#include <cstdint>
#include <iostream>
#include <vector>

using namespace std;

static constexpr uint64_t MOD = 888888883ULL;

static uint64_t mod_pow(uint64_t base, uint64_t exp) {
    uint64_t result = 1 % MOD;
    base %= MOD;
    while (exp > 0) {
        if (exp & 1ULL) result = (result * base) % MOD;
        base = (base * base) % MOD;
        exp >>= 1ULL;
    }
    return result;
}

static vector<uint64_t> build_factorials(int max_n, vector<uint64_t>& invfac) {
    vector<uint64_t> fac(max_n + 1, 1);
    for (int i = 1; i <= max_n; ++i) {
        fac[i] = (fac[i - 1] * static_cast<uint64_t>(i)) % MOD;
    }
    invfac.assign(max_n + 1, 1);
    invfac[max_n] = mod_pow(fac[max_n], MOD - 2);
    for (int i = max_n; i >= 1; --i) {
        invfac[i - 1] = (invfac[i] * static_cast<uint64_t>(i)) % MOD;
    }
    return fac;
}

static uint64_t multinomial_mod(int n, int a, int b, int c,
                                const vector<uint64_t>& fac,
                                const vector<uint64_t>& invfac) {
    uint64_t res = fac[n];
    res = (res * invfac[a]) % MOD;
    res = (res * invfac[b]) % MOD;
    res = (res * invfac[c]) % MOD;
    return res;
}

static uint64_t compute_N_mod(int X, int Y, int Z,
                              const vector<uint64_t>& fac,
                              const vector<uint64_t>& invfac) {
    int px = X & 1;
    int py = Y & 1;
    int pz = Z & 1;
    if (px != py || py != pz) return 0;

    int S = X + Y + Z;
    // Neutral strings satisfy inversion parity == ceil(S/2) mod 2.
    // Let total be the multinomial count and diff = even - odd.
    // For all even counts, diff equals the multinomial of half counts.
    // For all odd counts, diff is zero so both parities have total/2 strings.
    uint64_t total = multinomial_mod(S, X, Y, Z, fac, invfac);
    uint64_t inv2 = (MOD + 1) / 2;

    if (px == 1) {
        // All odd counts -> even/odd inversion counts are equal.
        return (total * inv2) % MOD;
    }

    // All even counts -> difference equals multinomial of half counts.
    uint64_t half = multinomial_mod(S / 2, X / 2, Y / 2, Z / 2, fac, invfac);
    if ((S / 2) & 1) {
        uint64_t diff = (total + MOD - half) % MOD;
        return (diff * inv2) % MOD;
    }
    uint64_t sum = total + half;
    if (sum >= MOD) sum -= MOD;
    return (sum * inv2) % MOD;
}

static __int128 factorial128(int n) {
    __int128 res = 1;
    for (int i = 2; i <= n; ++i) res *= i;
    return res;
}

static unsigned long long multinomial128(int n, int a, int b, int c) {
    __int128 res = factorial128(n);
    res /= factorial128(a);
    res /= factorial128(b);
    res /= factorial128(c);
    return static_cast<unsigned long long>(res);
}

static unsigned long long compute_N_exact_small(int X, int Y, int Z) {
    int px = X & 1;
    int py = Y & 1;
    int pz = Z & 1;
    if (px != py || py != pz) return 0;

    int S = X + Y + Z;
    unsigned long long total = multinomial128(S, X, Y, Z);

    if (px == 1) {
        return total / 2;
    }

    unsigned long long half = multinomial128(S / 2, X / 2, Y / 2, Z / 2);
    if ((S / 2) & 1) {
        return (total - half) / 2;
    }
    return (total + half) / 2;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    const int max_i = 87;
    vector<int> cubes(max_i + 1, 0);
    for (int i = 0; i <= max_i; ++i) {
        cubes[i] = i * i * i;
    }
    int max_sum = 3 * cubes[max_i];

    vector<uint64_t> invfac;
    vector<uint64_t> fac = build_factorials(max_sum, invfac);

    // Validation checkpoints from the statement.
    {
        unsigned long long n222 = compute_N_exact_small(2, 2, 2);
        unsigned long long n888 = compute_N_exact_small(8, 8, 8);
        if (n222 != 42ULL || n888 != 4732773210ULL) {
            cerr << "Validation failed: N(2,2,2)=" << n222
                 << ", N(8,8,8)=" << n888 << "\n";
            return 1;
        }
    }

    uint64_t answer = 0;
    for (int i = 0; i <= max_i; ++i) {
        int X = cubes[i];
        for (int j = 0; j <= max_i; ++j) {
            int Y = cubes[j];
            for (int k = 0; k <= max_i; ++k) {
                int Z = cubes[k];
                answer += compute_N_mod(X, Y, Z, fac, invfac);
                if (answer >= MOD) answer -= MOD;
            }
        }
    }

    cout << answer % MOD << "\n";
    return 0;
}

Python

import math

MOD = 888888883

def mod_pow(base, exp):
    return pow(base, exp, MOD)

def build_factorials(max_n):
    fac = [1] * (max_n + 1)
    for i in range(1, max_n + 1):
        fac[i] = (fac[i - 1] * i) % MOD
        
    invfac = [1] * (max_n + 1)
    invfac[max_n] = mod_pow(fac[max_n], MOD - 2)
    for i in range(max_n, 0, -1):
        invfac[i - 1] = (invfac[i] * i) % MOD
        
    return fac, invfac

def multinomial_mod(n, a, b, c, fac, invfac):
    res = fac[n]
    res = (res * invfac[a]) % MOD
    res = (res * invfac[b]) % MOD
    res = (res * invfac[c]) % MOD
    return res

def compute_N_mod(X, Y, Z, fac, invfac):
    px = X & 1
    py = Y & 1
    pz = Z & 1
    if px != py or py != pz:
        return 0
        
    S = X + Y + Z
    total = multinomial_mod(S, X, Y, Z, fac, invfac)
    inv2 = (MOD + 1) // 2
    
    if px == 1:
        return (total * inv2) % MOD
        
    half = multinomial_mod(S // 2, X // 2, Y // 2, Z // 2, fac, invfac)
    if (S // 2) & 1:
        diff = (total + MOD - half) % MOD
        return (diff * inv2) % MOD
        
    sum_val = total + half
    return (sum_val * inv2) % MOD

def compute_N_exact_small(X, Y, Z):
    px = X & 1
    py = Y & 1
    pz = Z & 1
    if px != py or py != pz:
        return 0
        
    S = X + Y + Z
    total = math.factorial(S) // (math.factorial(X) * math.factorial(Y) * math.factorial(Z))
    
    if px == 1:
        return total // 2
        
    half = math.factorial(S // 2) // (math.factorial(X // 2) * math.factorial(Y // 2) * math.factorial(Z // 2))
    if (S // 2) & 1:
        return (total - half) // 2
    return (total + half) // 2

def solve():
    max_i = 87
    cubes = [i * i * i for i in range(max_i + 1)]
    max_sum = 3 * cubes[-1]
    
    fac, invfac = build_factorials(max_sum)
    
    answer = 0
    for X in cubes:
        for Y in cubes:
            for Z in cubes:
                answer = (answer + compute_N_mod(X, Y, Z, fac, invfac)) % MOD
                
    return str(answer)

def run_checkpoints():
    assert compute_N_exact_small(2, 2, 2) == 42
    assert compute_N_exact_small(8, 8, 8) == 4732773210

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

Java

import java.math.BigInteger;

public class Euler981 {

    static final long MOD = 888888883L;

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

    static long[] fac;
    static long[] invfac;

    static void buildFactorials(int maxN) {
        fac = new long[maxN + 1];
        fac[0] = 1;
        for (int i = 1; i <= maxN; ++i) {
            fac[i] = (fac[i - 1] * i) % MOD;
        }
        invfac = new long[maxN + 1];
        invfac[maxN] = modPow(fac[maxN], MOD - 2);
        for (int i = maxN; i >= 1; --i) {
            invfac[i - 1] = (invfac[i] * i) % MOD;
        }
    }

    static long multinomialMod(int n, int a, int b, int c) {
        long res = fac[n];
        res = (res * invfac[a]) % MOD;
        res = (res * invfac[b]) % MOD;
        res = (res * invfac[c]) % MOD;
        return res;
    }

    static long computeNMod(int X, int Y, int Z) {
        int px = X & 1;
        int py = Y & 1;
        int pz = Z & 1;
        if (px != py || py != pz)
            return 0;

        int S = X + Y + Z;
        long total = multinomialMod(S, X, Y, Z);
        long inv2 = (MOD + 1) / 2;

        if (px == 1) {
            return (total * inv2) % MOD;
        }

        long half = multinomialMod(S / 2, X / 2, Y / 2, Z / 2);
        if ((S / 2) % 2 != 0) {
            long diff = (total + MOD - half) % MOD;
            return (diff * inv2) % MOD;
        }
        long sum = (total + half) % MOD;
        return (sum * inv2) % MOD;
    }

    static BigInteger factorialBig(int n) {
        BigInteger res = BigInteger.ONE;
        for (int i = 2; i <= n; ++i) {
            res = res.multiply(BigInteger.valueOf(i));
        }
        return res;
    }

    static BigInteger multinomialBig(int n, int a, int b, int c) {
        BigInteger res = factorialBig(n);
        res = res.divide(factorialBig(a));
        res = res.divide(factorialBig(b));
        res = res.divide(factorialBig(c));
        return res;
    }

    static long computeNExactSmall(int X, int Y, int Z) {
        int px = X & 1;
        int py = Y & 1;
        int pz = Z & 1;
        if (px != py || py != pz)
            return 0;

        int S = X + Y + Z;
        BigInteger total = multinomialBig(S, X, Y, Z);

        if (px == 1) {
            return total.divide(BigInteger.valueOf(2)).longValue();
        }

        BigInteger half = multinomialBig(S / 2, X / 2, Y / 2, Z / 2);
        if ((S / 2) % 2 != 0) {
            return total.subtract(half).divide(BigInteger.valueOf(2)).longValue();
        }
        return total.add(half).divide(BigInteger.valueOf(2)).longValue();
    }

    public static String solve() {
        int maxI = 87;
        int[] cubes = new int[maxI + 1];
        for (int i = 0; i <= maxI; ++i) {
            cubes[i] = i * i * i;
        }
        int maxSum = 3 * cubes[maxI];

        buildFactorials(maxSum);

        long answer = 0;
        for (int i = 0; i <= maxI; ++i) {
            int X = cubes[i];
            for (int j = 0; j <= maxI; ++j) {
                int Y = cubes[j];
                for (int k = 0; k <= maxI; ++k) {
                    int Z = cubes[k];
                    answer += computeNMod(X, Y, Z);
                    if (answer >= MOD)
                        answer -= MOD;
                }
            }
        }

        return Long.toString(answer % MOD);
    }

    public static void main(String[] args) {
        long n222 = computeNExactSmall(2, 2, 2);
        long n888 = computeNExactSmall(8, 8, 8);
        if (n222 != 42L || n888 != 4732773210L) {
            System.err.println("Validation failed: N(2,2,2)=" + n222 + ", N(8,8,8)=" + n888);
            return;
        }
        System.out.println(solve());
    }
}