Problem 281: Pizza Toppings

Project Euler Problem 281 Solution

EulerSolve provides an optimized solution for Project Euler Problem 281, Pizza Toppings, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For integers $m\ge2$ and $n\ge1$, we place $mn$ toppings around a circle, with exactly $n$ copies of each of the $m$ colors. Two arrangements are considered the same if one is a rotation of the other. Let $f(m,n)$ be the number of distinct rotational classes. The Project Euler task is to sum all values satisfying $$f(m,n)\le 10^{15}.$$ The final numeric sum is intentionally omitted here; the goal is to explain the counting formula and the search bounds used by the C++ solution. Mathematical Approach 1) Burnside's lemma is the natural tool The acting group is the cyclic rotation group $$C_{mn}=\{0,1,\dots,mn-1\},$$ where rotation by $r$ means shifting every position by $r$ places around the circle. Burnside's lemma says that the number of rotational classes is the average number of arrangements fixed by each rotation: $$f(m,n)=\frac{1}{mn}\sum_{r=0}^{mn-1}\mathrm{Fix}(r).$$ So the whole problem is to understand $\mathrm{Fix}(r)$ for one rotation $r$. 2) Cycle structure of one rotation Set $$N=mn,\qquad d=\gcd(N,r).$$ The rotation by $r$ decomposes the $N$ positions into $$d$$ disjoint cycles, each of length $$\ell=\frac{N}{d}.$$ This is the standard cycle decomposition of a cyclic shift....

Detailed mathematical approach

Problem Summary

For integers $m\ge2$ and $n\ge1$, we place $mn$ toppings around a circle, with exactly $n$ copies of each of the $m$ colors. Two arrangements are considered the same if one is a rotation of the other. Let $f(m,n)$ be the number of distinct rotational classes. The Project Euler task is to sum all values satisfying

$$f(m,n)\le 10^{15}.$$

The final numeric sum is intentionally omitted here; the goal is to explain the counting formula and the search bounds used by the C++ solution.

Mathematical Approach

1) Burnside's lemma is the natural tool

The acting group is the cyclic rotation group

$$C_{mn}=\{0,1,\dots,mn-1\},$$

where rotation by $r$ means shifting every position by $r$ places around the circle.

Burnside's lemma says that the number of rotational classes is the average number of arrangements fixed by each rotation:

$$f(m,n)=\frac{1}{mn}\sum_{r=0}^{mn-1}\mathrm{Fix}(r).$$

So the whole problem is to understand $\mathrm{Fix}(r)$ for one rotation $r$.

2) Cycle structure of one rotation

Set

$$N=mn,\qquad d=\gcd(N,r).$$

The rotation by $r$ decomposes the $N$ positions into

$$d$$

disjoint cycles, each of length

$$\ell=\frac{N}{d}.$$

This is the standard cycle decomposition of a cyclic shift.

An arrangement is fixed by this rotation if and only if every cycle is monochromatic: once one position in a cycle is chosen, all other positions in that cycle are forced to have the same color.

3) When can a rotation fix any valid arrangement?

Each color appears exactly $n$ times. If every cycle has length $\ell$, then every color must occupy a whole number of cycles. Therefore

$$\ell \mid n$$

is necessary.

It is also sufficient. If $\ell\mid n$, then each color must occupy exactly

$$\frac{n}{\ell}$$

cycles, because each such cycle contributes $\ell$ positions of that color.

Since there are

$$d=\frac{N}{\ell}=\frac{mn}{\ell}$$

cycles in total, a fixed arrangement exists precisely when $\ell\mid n$.

4) Counting the fixed arrangements for one cycle length

Assume $\ell\mid n$. Then the $d=\frac{mn}{\ell}$ cycles are distinguishable cycle slots, and we must color them so that each of the $m$ colors is used exactly $\frac{n}{\ell}$ times.

So the number of fixed arrangements is the multinomial coefficient

$$\mathrm{Fix}(\ell)=\frac{\left(\frac{mn}{\ell}\right)!}{\left(\frac{n}{\ell}!\right)^m}.$$

If $\ell\nmid n$, then

$$\mathrm{Fix}(\ell)=0.$$

5) How many rotations have the same cycle length?

Now fix a divisor $\ell$ of $N$. We want the number of rotations $r$ whose cycle length is exactly $\ell$, that is, whose gcd with $N$ equals $N/\ell$.

Write

$$r=\frac{N}{\ell}k.$$

Then

$$\gcd(N,r)=\frac{N}{\ell}\gcd(\ell,k).$$

So the condition $\gcd(N,r)=N/\ell$ is equivalent to

$$\gcd(\ell,k)=1.$$

The number of such residues $k$ modulo $\ell$ is exactly Euler's totient:

$$\varphi(\ell).$$

Therefore there are $\varphi(\ell)$ rotations with cycle length $\ell$.

6) The divisor formula used by the code

Only divisors $\ell\mid n$ contribute, so Burnside becomes

$$f(m,n)=\frac{1}{mn}\sum_{\ell\mid n}\varphi(\ell)\, \frac{(mn/\ell)!}{\bigl((n/\ell)!\bigr)^m}.$$

This is exactly the formula implemented by f_value() in the C++ solution.

7) Worked example: $f(3,2)=16$

Here $m=3$, $n=2$, so $N=6$. The divisors of $n$ are $1$ and $2$. Hence

$$f(3,2)=\frac{1}{6}\left[ \varphi(1)\frac{6!}{(2!)^3}+ \varphi(2)\frac{3!}{(1!)^3} \right].$$

Now

$$\varphi(1)=1,\qquad \varphi(2)=1,$$

$$f(3,2)=\frac{1}{6}\left[\frac{720}{8}+6\right] =\frac{1}{6}(90+6)=16.$$

This is one of the explicit checkpoints in the code.

8) Another small example: $f(2,3)=4$

With two colors and three copies of each,

$$f(2,3)=\frac{1}{6}\left[ \varphi(1)\frac{6!}{(3!)^2}+ \varphi(3)\frac{2!}{(1!)^2} \right].$$

Since $\varphi(1)=1$ and $\varphi(3)=2$, we get

$$f(2,3)=\frac{1}{6}(20+4)=4.$$

The brute-force checkpoint in the program confirms this value by explicit necklace enumeration.

9) Why the search over $n$ is finite

The Burnside sum is positive-term, so just the identity rotation $\ell=1$ already gives a lower bound:

$$f(m,n)\ge \frac{1}{mn}\frac{(mn)!}{(n!)^m}=:g(m,n).$$

For fixed $n$, this lower bound increases with $m$. Indeed,

$$\frac{g(m+1,n)}{g(m,n)} =\frac{m}{m+1}\cdot \frac{\prod_{k=1}^{n}(mn+k)}{n!} >1.$$

So among all $m\ge2$, the smallest possible value occurs at $m=2$. Therefore

$$f(m,n)\ge g(2,n)=\frac{(2n)!}{2n\,(n!)^2} =\frac{1}{2n}\binom{2n}{n}.$$

If this lower bound already exceeds $10^{15}$, then no value with that $n$ can contribute, for any $m\ge2$.

This is exactly why the code computes

$$\frac{1}{2n}\binom{2n}{n}$$

in compute_n_upper_bound().

10) Why the search over $m$ is finite

Take $n=1$. Then we place $m$ distinct colors once each around the circle, modulo rotation. The number of circular orders is

$$f(m,1)=(m-1)!.$$

So if

$$ (m-1)! > 10^{15}, $$

then even the smallest possible $n$ already overshoots the limit, and this $m$ can never contribute.

That is why compute_m_upper_bound() effectively finds the largest $m$ with

$$ (m-1)! \le 10^{15}. $$

Code Logic

1) Totients. compute_totients() sieves all $\varphi(\ell)$ values up to the maximum relevant $n$.

2) Factorials and factorial powers. factorials_up_to() precomputes $k!$, and precompute_fact_powers() stores $(n!)^m$. This avoids repeating huge multiplications in every Burnside term.

3) Burnside evaluator. f_value(m,n,...) loops over all divisors $\ell\mid n$, computes the multinomial fixed-count term, multiplies by $\varphi(\ell)$, sums, and divides by $mn$.

4) Checkpoints. The code explicitly verifies

$$f(2,1)=1,\qquad f(2,2)=2,\qquad f(3,1)=2,\qquad f(3,2)=16,$$

and then compares the formula with brute-force necklace counting for

$$ (m,n)\in\{(2,3),(2,4),(3,2),(4,2)\}. $$

5) Final scan. After the safe upper bounds are known, the program checks every $(m,n)$ in that rectangle and adds only values $\le10^{15}$.

Complexity Analysis

For each pair $(m,n)$, the work is proportional to the number of divisors of $n$, plus big-integer arithmetic on factorial ratios. The precomputed totients, factorials, and factorial powers greatly reduce repeated work. Memory usage is dominated by those precomputed big-integer tables.

Problem 281 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <set>
#include <string>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using boost::multiprecision::cpp_int;

constexpr u64 kLimit = 1000000000000000ULL;

struct Options {
    bool run_checkpoints = true;
};

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

std::vector<int> compute_totients(const int n_max) {
    std::vector<int> phi(static_cast<std::size_t>(n_max + 1));
    for (int i = 0; i <= n_max; ++i) {
        phi[static_cast<std::size_t>(i)] = i;
    }

    for (int p = 2; p <= n_max; ++p) {
        if (phi[static_cast<std::size_t>(p)] != p) {
            continue;
        }
        for (int multiple = p; multiple <= n_max; multiple += p) {
            phi[static_cast<std::size_t>(multiple)] -= phi[static_cast<std::size_t>(multiple)] / p;
        }
    }

    return phi;
}

cpp_int central_binomial(const int n) {
    cpp_int c = 1;
    for (int i = 1; i <= n; ++i) {
        c *= (n + i);
        c /= i;
    }
    return c;
}

int compute_n_upper_bound(const u64 limit) {
    const cpp_int limit_big = limit;
    int n = 1;
    int best = 0;
    while (true) {
        const cpp_int bound = central_binomial(n) / (2 * n);
        if (bound > limit_big) {
            break;
        }
        best = n;
        ++n;
    }
    return best;
}

int compute_m_upper_bound(const u64 limit) {
    const cpp_int limit_big = limit;

    cpp_int factorial = 1;  // 1! = 1
    int best = 1;

    for (int m = 2;; ++m) {
        if (factorial > limit_big) {
            break;
        }
        best = m;
        factorial *= m;  // factorial now becomes m!
    }

    return best;
}

std::vector<cpp_int> factorials_up_to(const int n_max) {
    std::vector<cpp_int> fact(static_cast<std::size_t>(n_max + 1));
    fact[0] = 1;
    for (int i = 1; i <= n_max; ++i) {
        fact[static_cast<std::size_t>(i)] = fact[static_cast<std::size_t>(i - 1)] * i;
    }
    return fact;
}

std::vector<std::vector<cpp_int>> precompute_fact_powers(
    const std::vector<cpp_int>& fact,
    const int n_max,
    const int m_max
) {
    std::vector<std::vector<cpp_int>> powers(
        static_cast<std::size_t>(n_max + 1),
        std::vector<cpp_int>(static_cast<std::size_t>(m_max + 1), cpp_int(1))
    );

    for (int n = 0; n <= n_max; ++n) {
        powers[static_cast<std::size_t>(n)][0] = 1;
        for (int m = 1; m <= m_max; ++m) {
            powers[static_cast<std::size_t>(n)][static_cast<std::size_t>(m)] =
                powers[static_cast<std::size_t>(n)][static_cast<std::size_t>(m - 1)] *
                fact[static_cast<std::size_t>(n)];
        }
    }

    return powers;
}

cpp_int f_value(
    const int m,
    const int n,
    const std::vector<int>& phi,
    const std::vector<cpp_int>& fact,
    const std::vector<std::vector<cpp_int>>& fact_powers
) {
    const int total = m * n;
    cpp_int sum = 0;

    for (int l = 1; l <= n; ++l) {
        if ((n % l) != 0) {
            continue;
        }

        const int slice = n / l;
        const cpp_int numerator = fact[static_cast<std::size_t>(total / l)];
        const cpp_int denominator = fact_powers[static_cast<std::size_t>(slice)][static_cast<std::size_t>(m)];
        const cpp_int multinomial = numerator / denominator;

        sum += static_cast<cpp_int>(phi[static_cast<std::size_t>(l)]) * multinomial;
    }

    return sum / total;
}

std::vector<int> rotate_left(const std::vector<int>& v, const int shift) {
    const int n = static_cast<int>(v.size());
    std::vector<int> out(static_cast<std::size_t>(n));
    for (int i = 0; i < n; ++i) {
        out[static_cast<std::size_t>(i)] = v[static_cast<std::size_t>((i + shift) % n)];
    }
    return out;
}

u64 brute_count_necklaces(const int m, const int n) {
    std::vector<int> arrangement;
    arrangement.reserve(static_cast<std::size_t>(m * n));
    for (int color = 0; color < m; ++color) {
        for (int cnt = 0; cnt < n; ++cnt) {
            arrangement.push_back(color);
        }
    }

    std::set<std::vector<int>> canonical;

    std::sort(arrangement.begin(), arrangement.end());
    do {
        std::vector<int> best = arrangement;
        for (int shift = 1; shift < m * n; ++shift) {
            std::vector<int> r = rotate_left(arrangement, shift);
            if (r < best) {
                best = std::move(r);
            }
        }
        canonical.insert(std::move(best));
    } while (std::next_permutation(arrangement.begin(), arrangement.end()));

    return static_cast<u64>(canonical.size());
}

bool run_checkpoints(
    const std::vector<int>& phi,
    const std::vector<cpp_int>& fact,
    const std::vector<std::vector<cpp_int>>& fact_powers
) {
    struct Case {
        int m;
        int n;
        u64 expected;
    };

    const std::vector<Case> explicit_cases = {
        {2, 1, 1},
        {2, 2, 2},
        {3, 1, 2},
        {3, 2, 16},
    };

    for (const Case& c : explicit_cases) {
        const cpp_int value = f_value(c.m, c.n, phi, fact, fact_powers);
        if (value != c.expected) {
            std::cerr << "Explicit checkpoint failed for f(" << c.m << ',' << c.n
                      << "): got " << value << ", expected " << c.expected << '\n';
            return false;
        }
    }

    const std::vector<std::pair<int, int>> brute_cases = {
        {2, 3},
        {2, 4},
        {3, 2},
        {4, 2},
    };

    for (const auto& [m, n] : brute_cases) {
        const cpp_int formula = f_value(m, n, phi, fact, fact_powers);
        const u64 brute = brute_count_necklaces(m, n);
        if (formula != brute) {
            std::cerr << "Brute checkpoint failed for f(" << m << ',' << n
                      << "): formula=" << formula << ", brute=" << brute << '\n';
            return false;
        }
    }

    return true;
}

u64 solve() {
    const int n_upper = compute_n_upper_bound(kLimit);
    const int m_upper = compute_m_upper_bound(kLimit);

    const int max_factorial_needed = m_upper * n_upper;
    const std::vector<cpp_int> fact = factorials_up_to(max_factorial_needed);
    const std::vector<int> phi = compute_totients(n_upper);
    const std::vector<std::vector<cpp_int>> fact_powers = precompute_fact_powers(fact, n_upper, m_upper);

    if (!run_checkpoints(phi, fact, fact_powers)) {
        return 0;
    }

    const cpp_int limit_big = kLimit;
    u64 answer = 0;

    for (int m = 2; m <= m_upper; ++m) {
        for (int n = 1; n <= n_upper; ++n) {
            const cpp_int value = f_value(m, n, phi, fact, fact_powers);
            if (value <= limit_big) {
                answer += value.convert_to<u64>();
            }
        }
    }

    return answer;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    const int n_upper = compute_n_upper_bound(kLimit);
    const int m_upper = compute_m_upper_bound(kLimit);
    const int max_factorial_needed = m_upper * n_upper;

    const std::vector<cpp_int> fact = factorials_up_to(max_factorial_needed);
    const std::vector<int> phi = compute_totients(n_upper);
    const std::vector<std::vector<cpp_int>> fact_powers = precompute_fact_powers(fact, n_upper, m_upper);

    if (options.run_checkpoints && !run_checkpoints(phi, fact, fact_powers)) {
        return 2;
    }

    const cpp_int limit_big = kLimit;
    u64 answer = 0;

    for (int m = 2; m <= m_upper; ++m) {
        for (int n = 1; n <= n_upper; ++n) {
            const cpp_int value = f_value(m, n, phi, fact, fact_powers);
            if (value <= limit_big) {
                answer += value.convert_to<u64>();
            }
        }
    }

    std::cout << answer << '\n';
    return 0;
}

Python

import math

def compute_totients(n_max):
    phi = list(range(n_max + 1))
    for p in range(2, n_max + 1):
        if phi[p] == p:
            for multiple in range(p, n_max + 1, p):
                phi[multiple] -= phi[multiple] // p
    return phi

def central_binomial(n):
    return math.comb(2 * n, n)

def compute_n_upper_bound(limit):
    n = 1
    best = 0
    while True:
        bound = central_binomial(n) // (2 * n)
        if bound > limit:
            break
        best = n
        n += 1
    return best

def compute_m_upper_bound(limit):
    factorial = 1
    best = 1
    m = 2
    while True:
        if factorial > limit:
            break
        best = m
        factorial *= m
        m += 1
    return best

def f_value(m, n, phi, fact, fact_powers):
    total = m * n
    sum_val = 0
    
    for l in range(1, n + 1):
        if n % l != 0:
            continue
            
        slice_len = n // l
        numerator = fact[total // l]
        denominator = fact_powers[slice_len][m]
        multinomial = numerator // denominator
        
        sum_val += phi[l] * multinomial
        
    return sum_val // total

def solve(limit=1000000000000000):
    n_upper = compute_n_upper_bound(limit)
    m_upper = compute_m_upper_bound(limit)
    
    max_factorial_needed = m_upper * n_upper
    fact = [1] * (max_factorial_needed + 1)
    for i in range(1, max_factorial_needed + 1):
        fact[i] = fact[i - 1] * i
        
    phi = compute_totients(n_upper)
    
    fact_powers = [[1] * (m_upper + 1) for _ in range(n_upper + 1)]
    for n in range(n_upper + 1):
        fact_powers[n][0] = 1
        for m in range(1, m_upper + 1):
            fact_powers[n][m] = fact_powers[n][m - 1] * fact[n]
            
    answer = 0
    for m in range(2, m_upper + 1):
        for n in range(1, n_upper + 1):
            val = f_value(m, n, phi, fact, fact_powers)
            if val <= limit:
                answer += val
                
    return str(answer)

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

Java

import java.math.BigInteger;

public class Euler281 {
    static final long LIMIT = 1000000000000000L;

    static int[] computeTotients(int nMax) {
        int[] phi = new int[nMax + 1];
        for (int i = 0; i <= nMax; ++i)
            phi[i] = i;

        for (int p = 2; p <= nMax; ++p) {
            if (phi[p] == p) {
                for (int multiple = p; multiple <= nMax; multiple += p) {
                    phi[multiple] -= phi[multiple] / p;
                }
            }
        }
        return phi;
    }

    static BigInteger centralBinomial(int n) {
        BigInteger c = BigInteger.ONE;
        for (int i = 1; i <= n; ++i) {
            c = c.multiply(BigInteger.valueOf(n + i))
                    .divide(BigInteger.valueOf(i));
        }
        return c;
    }

    static int computeNUpperBound(long limit) {
        BigInteger limitBig = BigInteger.valueOf(limit);
        int n = 1;
        int best = 0;
        while (true) {
            BigInteger bound = centralBinomial(n).divide(BigInteger.valueOf(2L * n));
            if (bound.compareTo(limitBig) > 0)
                break;
            best = n;
            ++n;
        }
        return best;
    }

    static int computeMUpperBound(long limit) {
        BigInteger limitBig = BigInteger.valueOf(limit);
        BigInteger factorial = BigInteger.ONE;
        int best = 1;
        for (int m = 2;; ++m) {
            if (factorial.compareTo(limitBig) > 0)
                break;
            best = m;
            factorial = factorial.multiply(BigInteger.valueOf(m));
        }
        return best;
    }

    static BigInteger[] factorialsUpTo(int nMax) {
        BigInteger[] fact = new BigInteger[nMax + 1];
        fact[0] = BigInteger.ONE;
        for (int i = 1; i <= nMax; ++i) {
            fact[i] = fact[i - 1].multiply(BigInteger.valueOf(i));
        }
        return fact;
    }

    static BigInteger[][] precomputeFactPowers(BigInteger[] fact, int nMax, int mMax) {
        BigInteger[][] powers = new BigInteger[nMax + 1][mMax + 1];
        for (int n = 0; n <= nMax; ++n) {
            powers[n][0] = BigInteger.ONE;
            for (int m = 1; m <= mMax; ++m) {
                powers[n][m] = powers[n][m - 1].multiply(fact[n]);
            }
        }
        return powers;
    }

    static BigInteger fValue(int m, int n, int[] phi, BigInteger[] fact, BigInteger[][] factPowers) {
        int total = m * n;
        BigInteger sum = BigInteger.ZERO;

        for (int l = 1; l <= n; ++l) {
            if (n % l != 0)
                continue;

            int slice = n / l;
            BigInteger numerator = fact[total / l];
            BigInteger denominator = factPowers[slice][m];
            BigInteger multinomial = numerator.divide(denominator);

            sum = sum.add(multinomial.multiply(BigInteger.valueOf(phi[l])));
        }

        return sum.divide(BigInteger.valueOf(total));
    }

    public static String solve() {
        int nUpper = computeNUpperBound(LIMIT);
        int mUpper = computeMUpperBound(LIMIT);

        int maxFactorialNeeded = mUpper * nUpper;
        BigInteger[] fact = factorialsUpTo(maxFactorialNeeded);
        int[] phi = computeTotients(nUpper);
        BigInteger[][] factPowers = precomputeFactPowers(fact, nUpper, mUpper);

        BigInteger limitBig = BigInteger.valueOf(LIMIT);
        long answer = 0;

        for (int m = 2; m <= mUpper; ++m) {
            for (int n = 1; n <= nUpper; ++n) {
                BigInteger value = fValue(m, n, phi, fact, factPowers);
                if (value.compareTo(limitBig) <= 0) {
                    answer += value.longValue();
                }
            }
        }

        return String.valueOf(answer);
    }

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