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

A round pizza is cut into $mn$ equal slices, and every slice gets exactly one of $m$ different toppings, each topping on exactly $n$ slices, where $m\ge2$ and $n\ge1$. Two arrangements count as the same exactly when one is a rotation of the other; reflections are not used, so an arrangement and its mirror image count as different unless they are also rotations of each other. With $f(m,n)$ the number of different arrangements, the statement gives $f(2,1)=1$, $f(2,2)=f(3,1)=2$ and $f(3,2)=16$, and asks for the sum of all values $f(m,n)$ that do not exceed $10^{15}$. ^[1] In the language of combinatorics, $f(m,n)$ counts the necklaces of $mn$ beads with $n$ beads of each of $m$ colours, which may be turned but not flipped over. ^[2]

The whole solution in one sentence: the number of arrangements up to rotation is the average, over all rotations, of the number of arrangements that a rotation leaves unchanged (Burnside's lemma), which gives a short formula for $f(m,n)$, and since $f(m,n)$ grows quickly in both $m$ and $n$, only the pairs with $m\le18$ and $n\le29$ have to be checked against the limit.

Mathematical Approach

We first count the arrangements on a straight strip, then correct for rotations, and finally bound the pairs $(m,n)$ that can satisfy the limit.

1. Small cases by hand

Number the slices $0,1,\dots,N-1$ around the pizza, with $N=mn$, and write an arrangement as a word $a_0a_1\cdots a_{N-1}$ whose letters are the toppings. The rotation by $r$ slices, for $0\le r\lt N$, turns this word into $a_ra_{r+1}\cdots a_{N-1}a_0\cdots a_{r-1}$. The words that can be turned into one another form a class, and $f(m,n)$ is the number of classes. With the toppings $\mathtt A$, $\mathtt B$ and $\mathtt C$:

$f(2,2)$: of the six words with two $\mathtt A$ and two $\mathtt B$, the words $\mathtt{AABB}$, $\mathtt{ABBA}$, $\mathtt{BBAA}$ and $\mathtt{BAAB}$ are rotations of one another, and so are $\mathtt{ABAB}$ and $\mathtt{BABA}$: two classes, of sizes $4$ and $2$.
$f(3,1)$: the classes $\{\mathtt{ABC},\mathtt{BCA},\mathtt{CAB}\}$ and $\{\mathtt{ACB},\mathtt{CBA},\mathtt{BAC}\}$ are mirror images of each other; since reflections are distinct, they are two different arrangements.
$f(3,2)$: there are $90$ words with two copies of each topping (Section 2). If every class had six words, there would be $90/6=15$ classes, but the answer is $16$. The words $\mathtt{ABCABC}$, $\mathtt{BCABCA}$ and $\mathtt{CABCAB}$ repeat after three slices, so they form a class of only three words, and so do $\mathtt{ACBACB}$, $\mathtt{CBACBA}$ and $\mathtt{BACBAC}$. The other $84$ words form $14$ classes of six, and $14+2=16$.

The classes have different sizes, and dividing by $N$ fails exactly because of the words that a rotation by some $r\ne0$ leaves unchanged.

2. Words with repeated letters

On a strip that cannot be turned, the number of words with $n$ copies of each of $m$ letters is the multinomial coefficient

$$W(m,n)=\frac{(mn)!}{(n!)^m}.$$

To see this, give the $n$ copies of each letter temporary labels. The $mn$ labelled copies can be put in a row in $(mn)!$ ways, and every word arises exactly $(n!)^m$ times, once for each way to permute the labels within each of the $m$ letters. ^[3] For example $W(2,2)=4!/(2!)^2=6$ and $W(3,2)=6!/(2!)^3=90$, the numbers of words in Section 1.

3. Counting classes by fixed words

For a rotation $r$, let $\mathrm{Fix}(r)$ be the set of words that it leaves unchanged. Burnside's lemma, also called the orbit-counting theorem, says that the number of classes is the average size of these sets: ^[4]

$$f(m,n)=\frac1N\sum_{r=0}^{N-1}\bigl|\mathrm{Fix}(r)\bigr|.$$

For rotations the proof is short. Fix a word $x$ and let $p$ be the smallest $r\ge1$ for which the rotation by $r$ leaves $x$ unchanged; the full turn $r=N$ always does. If two rotations leave $x$ unchanged, so do their sum and difference modulo $N$. Dividing with remainder, $N=qp+s$ with $0\le s\lt p$, shows that the rotation by $s$ leaves $x$ unchanged too, so $s=0$ by the choice of $p$; in the same way, the rotations that fix $x$ are exactly the $N/p$ multiples of $p$. The rotations by $0,1,\dots,p-1$ give $p$ different words, because two equal ones would give a fixing rotation smaller than $p$, and every other rotation gives the same word as the rotation by its remainder modulo $p$; so the class of $x$ has exactly $p$ words. Now count the pairs of a rotation and a word that it fixes. Rotation by rotation, there are $\sum_r|\mathrm{Fix}(r)|$ of them. Word by word, each of the $p$ words of a class of size $p$ is fixed by $N/p$ rotations, so the class contributes $p\cdot N/p=N$ pairs. Hence the sum is $N$ times the number of classes.

For $f(3,2)$, the rotation by $0$ fixes all $90$ words, the half turn $r=3$ fixes the $6$ words that repeat after three slices, and the other four rotations fix none, so $f(3,2)=(90+6)/6=16$.

4. The cycles of one rotation

Under the rotation by $r$, slice $i$ receives the topping of slice $i+r$, counted modulo $N$, and following $i\to i+r\to i+2r\to\cdots$ splits the slices into cycles. Claim: with $d=\gcd(r,N)$, there are $d$ cycles, each of length $\ell=N/d$. Indeed, we come back to $i$ after $k$ steps exactly when $N$ divides $kr$, that is, when $N/d$ divides $k\cdot r/d$. Since $N/d$ and $r/d$ have no common factor, this happens exactly when $N/d$ divides $k$. So every cycle has length $N/d$, and there are $d$ of them. For $N=6$, the rotation by $1$ has the single cycle $\{0,1,2,3,4,5\}$, the rotation by $2$ has $\{0,2,4\}$ and $\{1,3,5\}$, and the rotation by $3$ has $\{0,3\}$, $\{1,4\}$ and $\{2,5\}$; the rotation by $0$ has $\gcd(0,N)=N$ cycles of length $1$.

5. The words fixed by one rotation

A word is fixed by the rotation by $r$ exactly when $a_{i+r}=a_i$ for every $i$, that is, when all slices of each cycle carry the same topping. Then every topping covers whole cycles, so its $n$ slices are a union of cycles of length $\ell$, and $\ell$ must divide $n$. If it does, each topping covers exactly $n/\ell$ of the $N/\ell=m\cdot n/\ell$ cycles, and choosing these cycles is the problem of Section 2 with the cycles in place of the slices:

$$\bigl|\mathrm{Fix}(r)\bigr|=\begin{cases}W(m,n/\ell),&\ell\mid n,\\0,&\ell\nmid n.\end{cases}$$

For $m=3$ and $n=2$, the half turn has $\ell=2$, and its three cycles get the three toppings in $W(3,1)=3!=6$ ways, the six words of Section 1. The rotations by $2$ and $4$ have $\ell=3$, and those by $1$ and $5$ have $\ell=6$; neither length divides $2$.

6. How many rotations share a cycle length

Let $\ell$ divide $N$. The rotation by $r$ has cycle length $\ell$ when $\gcd(r,N)=N/\ell$, so $r$ is a multiple $r=kN/\ell$ with $0\le k\lt\ell$. Then $\gcd(r,N)=(N/\ell)\gcd(k,\ell)$, which equals $N/\ell$ exactly when $\gcd(k,\ell)=1$. Replacing $k=0$ by $k=\ell$ changes no greatest common divisor, so the number of such $k$ is Euler's totient $\varphi(\ell)$, the number of integers among $1,\dots,\ell$ that have no common factor with $\ell$. ^[5] For $N=6$, the cycle lengths $1$, $2$, $3$ and $6$ belong to $r=0$, to $r=3$, to $r=2,4$ and to $r=1,5$, that is, to $\varphi(1)=1$, $\varphi(2)=1$, $\varphi(3)=2$ and $\varphi(6)=2$ rotations.

7. The divisor formula

Sections 3, 5 and 6 together give, since every divisor of $n$ also divides $N=mn$,

$$f(m,n)=\frac1{mn}\sum_{\ell\mid n}\varphi(\ell)\,W\Bigl(m,\frac n\ell\Bigr).$$

The division is exact, because by Section 3 the sum is $mn$ times a count. The same formula, written for arbitrary numbers of beads of each colour, is the standard count of necklaces with prescribed content. ^[2] It reproduces the four values of the statement: $f(2,1)=2/2=1$, $f(2,2)=(6+2)/4=2$, $f(3,1)=6/3=2$ and $f(3,2)=(90+6)/6=16$. A complete example with four divisors is $f(2,6)$, where $N=12$:

$\ell=1$: the rotation by $0$ fixes all $W(2,6)=12!/(6!)^2=924$ words;
$\ell=2$: the rotation by $6$ has six cycles of two slices and fixes $W(2,3)=20$ words;
$\ell=3$: the $\varphi(3)=2$ rotations by $4$ and $8$ have four cycles and fix $W(2,2)=6$ words each;
$\ell=6$: the $\varphi(6)=2$ rotations by $2$ and $10$ have two cycles and fix $W(2,1)=2$ words each.

The other six rotations have the cycle lengths $4$ and $12$, which do not divide $6$, so $f(2,6)=(924+20+12+4)/12$${}=80$. For a prime $n=p$ only $\ell=1$ and $\ell=p$ occur, and $W(m,1)=m!$, so

$$f(m,p)=\frac{W(m,p)+(p-1)\,m!}{mp}.$$

For instance $f(2,3)=(20+2\cdot2)/6=4$, with the classes of $\mathtt{AAABBB}$, $\mathtt{AABABB}$, $\mathtt{AABBAB}$ and $\mathtt{ABABAB}$; the second and the third are mirror images of each other. The values $1,2,4,10,26,80,\dots$ for $m=2$ and $2,16,188,2896,\dots$ for $m=3$ are listed in the OEIS. ^[6]

8. How fast the counts grow

All terms of the divisor formula are nonnegative, so the term $\ell=1$ alone gives the lower bound

$$\begin{aligned}f(m,n)\ge g(m,n)&=\frac{W(m,n)}{mn}\\&=\frac{(mn)!}{mn\,(n!)^m}.\end{aligned}$$

This bound increases in both directions. When $n$ grows by one,

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

Every factor of the product is at least $1$. Since $m\ge2$, the factors for $k=1$ and $k=2$ are at least $\frac{2n+1}{n+1}$ and $2$, so the product is at least $\frac{2(2n+1)}{n+1}\ge3$, while $\frac{n}{n+1}\ge\frac12$; the ratio is therefore at least $3/2$. When $m$ grows by one, the ratio $g(m+1,n)/g(m,n)$ is

$$\frac{m}{m+1}\binom{(m+1)n}{n}\ge mn\ge2.$$

Here the quotient $\binom K{j+1}/\binom Kj=(K-j)/(j+1)$ is at least $1$ for $j\le(K-1)/2$, so the entries of a row of Pascal's triangle increase towards the middle, and by the symmetry $\binom Kj=\binom K{K-j}$ they decrease after it. Hence $\binom{K}{j}\ge\binom K1=K$ for $1\le j\le K-1$, and with $K=(m+1)n$ and $j=n$ the binomial coefficient is at least $(m+1)n$. So $g(m,n)\ge g(m,1)$ and $g(m,n)\ge g(2,n)$, and both ends are explicit, $g(m,1)=m!/m=(m-1)!$ and $g(2,n)=\binom{2n}{n}/(2n)$:

$$\begin{aligned}f(m,n)&\ge(m-1)!,\\f(m,n)&\ge\binom{2n}{n}\Big/(2n).\end{aligned}$$

The first bound is attained at $n=1$: $f(m,1)=(m-1)!$, the number of ways to seat $m$ people at a round table.

9. The search box

Let $L=10^{15}$. Both bounds of Section 8 increase: the step from $m$ to $m+1$ multiplies $(m-1)!$ by $m\ge2$, and $g(2,n)$ grows by the ratio estimate. So if $M$ is the largest $m$ with $(m-1)!\le L$ and $K$ is the largest $n$ with $\binom{2n}{n}/(2n)\le L$, every pair with $f(m,n)\le L$ has $m\le M$ and $n\le K$. Since $17!=355\,687\,428\,096\,000$ is at most $L$ and $18!=6\,402\,373\,705\,728\,000$ is not, $M=18$. The central binomial coefficient grows roughly like $4^n/\sqrt{\pi n}$, ^[7] so $K$ is small: $\binom{58}{29}/58\approx5.18\cdot10^{14}$ is at most $L$, while $\binom{60}{30}/60\approx1.97\cdot10^{15}$ is not, so $K=29$.

The box $2\le m\le18$, $1\le n\le29$ holds $17\cdot29=493$ pairs, and $74$ of them satisfy $f(m,n)\le L$, from all $n\le29$ for $m=2$ and the $n\le12$ for $m=3$ down to $n=1$ alone for $11\le m\le18$. Both corners are used: $f(18,1)=17!$ is below the limit, and so is $f(2,29)=\bigl(\binom{58}{29}+56\bigr)/58$ by the prime case of Section 7. So the box cannot be made smaller.

How the Code Works

The C++ program ^[8] evaluates the divisor formula of Section 7 on the box of Section 9. After parse_arguments, which accepts only --skip-checkpoints, main calls compute_n_upper_bound, which increases $n$ until central_binomial(n) / (2 * n) exceeds kLimit $=10^{15}$ and returns the last $n$ before that, and compute_m_upper_bound, which multiplies $1\cdot2\cdot3\cdots$ and returns the last $m$ with $(m-1)!\le L$. central_binomial multiplies by $n+i$ and divides by $i$ for $i=1,\dots,n$; the value after step $i$ is $\binom{n+i}{i}$, so every division is exact. Rounding the quotient by $2n$ down could only add an $n$ to the box, never lose one, and for this limit it changes nothing.

factorials_up_to stores $k!$ for $k\le MK=522$, and precompute_fact_powers stores $(k!)^j$ for $k\le K$ and $j\le M$. compute_totients starts with every entry equal to its index. An entry $p$ that is still unchanged when the loop reaches it is a prime, since every composite number has been reduced by its smallest prime factor, and then each multiple $k$ of $p$ gets phi[k] -= phi[k] / p, an exact multiplication by $1-1/p$, as only smaller primes were divided out before. This is Euler's product formula $\varphi(k)=k\prod_{p\mid k}(1-1/p)$. ^[5]

f_value(m, n, ...) is the formula of Section 7: for each divisor l of $n$ it divides fact[total / l] by the stored $((n/l)!)^m$, which gives $W(m,n/l)$, adds $\varphi(l)$ times this, and at the end divides by total $=mn$. Unless the option is given, run_checkpoints compares f_value with the four values of the statement and with brute_count_necklaces for $(2,3)$, $(2,4)$, $(3,2)$ and $(4,2)$. That function visits all distinct words with std::next_permutation, starting from the sorted word, ^[9] keeps the smallest rotation of each, built with rotate_left, in a std::set and returns its size; a failed check ends the program with exit code $2$. Finally main evaluates f_value on the whole box, adds the values that do not exceed the limit in a 64-bit unsigned integer and prints the sum. A function solve repeats these steps but is never called.

The formula needs arbitrary-precision integers, cpp_int of Boost.Multiprecision in C++, BigInteger in Java and the built-in integers of Python: ^[10] $522!$ has $1194$ digits, and $f(18,29)$ has $634$. The running sum, u64 in C++ and long in Java, would be safe even if all $493$ pairs qualified, since every added value is at most $10^{15}$ and $493\cdot10^{15}\lt2^{63}\approx9.2\cdot10^{18}$. Python and Java compute the same bounds, tables and formula, without the checkpoints and the option; Python builds the tables inside solve and takes $\binom{2n}{n}$ from math.comb, Java uses camel-case names such as fValue, and both return the sum from solve as a string.

Complexity and Verification

For each pair, f_value tests every $l\le n$ for divisibility, $7395$ cheap tests in all, and computes a term only for the divisors. The numbers $n\le29$ have $103$ divisors together, at most $8$ for $n=24$, so the $493$ pairs need $17\cdot103=1751$ terms, each an exact division, a multiplication and an addition of integers with at most $1194$ digits, and one final division per pair; in terms of the bounds, these are $O(MK^2)$ tests and $O(MK\log K)$ terms. The tables take $522$ multiplications for the factorials and $540$ for the powers, and they are the whole memory: $523$ factorials, $570$ powers and $30$ totients, about $140$ kilobytes of integer data. On a desktop computer all three programs finish within a few hundredths of a second, mostly start-up time.

The C++ program checks the four values of the statement and the brute-force counts $f(2,3)=4$, $f(2,4)=10$, $f(3,2)=16$ and $f(4,2)=318$, from $20$, $70$, $90$ and $2520$ words; the Python and Java versions contain no checks. For this article the following was verified separately:

a separately written brute force agrees with the divisor formula for all $23$ pairs with $mn\le12$, up to the $479\,001\,600$ words behind $f(12,1)=11!$;
the formula gives the OEIS values for $m=2$, $n\le14$ and for $m=3$, $n\le9$; ^[6]
the two ratios of Section 8, computed exactly for $2\le m\le40$ and $1\le n\le80$, are never below $3/2$ and $2$, and no pair of this range outside the box has $f(m,n)\le10^{15}$;
the sieve gives the correct $\varphi(k)$ for every $k\le1000$;
the C++ program built with GCC and with Clang, also with --skip-checkpoints, and the Python and Java programs print the same sum.

Footnotes and References

The references below supply the background facts used above. The proof of Burnside's lemma for rotations, the count of fixed words, the growth bounds and the search box are derived explicitly in this article.

Project Euler 281 — Pizza Toppings. The official statement: a pizza of m·n equal slices gets one topping per slice, m ≥ 2 different toppings, each on exactly n ≥ 1 slices; reflections are distinct, rotations are not. It gives f(2, 1) = 1, f(2, 2) = f(3, 1) = 2 and f(3, 2) = 16 and asks for the sum of all f(m, n) ≤ 10¹⁵.
Wikipedia — Necklace (combinatorics). Necklaces are classes of strings under rotation, bracelets under rotation and reflection. For a multiset with nᵢ beads of colour i, the article counts the necklaces by a sum of φ(d) times multinomial coefficients over the divisors d of gcd(n₁, …, nₖ), divided by the number of beads; for every nᵢ = n it is the formula of Section 7.
Wikipedia — Multinomial theorem. The multinomial coefficient is the number of distinct permutations of a multiset, illustrated with the letters of MISSISSIPPI.
Wikipedia — Burnside's lemma. Also called the orbit-counting theorem or the Cauchy–Frobenius lemma: the number of orbits is the average number of points fixed by a group element, proved by counting fixed pairs in two ways.
Wikipedia — Euler's totient function. φ(n) counts the k with 1 ≤ k ≤ n and gcd(n, k) = 1, and Euler's product formula gives φ(n) = n ∏ (1 − 1/p) over the distinct primes p dividing n.
OEIS A003239 is also the number of necklaces with 2n beads, n white and n black, and OEIS A118644 the number of distinct (n red, n blue, n green)-bead necklaces, 1, 2, 16, 188, 2896, …
Wikipedia — Central binomial coefficient. The article gives C(2n, n) = (4ⁿ / √(πn)) (1 − 1/(8n) + O(n⁻²)).
C++ — Euler281.cpp. The linked immutable C++ revision contains f_value, the two bound functions and run_checkpoints. It is the implementation record for this article.
cppreference — std::next_permutation. It turns a range into the lexicographically next permutation and returns false after the last one; the example on the page lists the six arrangements of AABB.
Boost.Multiprecision — cpp_int, Java — BigInteger and Python — numeric types describe the three arbitrary-precision integer types used by the programs.

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());
    }
}