Detailed mathematical approach

Problem Summary

Let \(F_1=F_2=1\). For each integer \(t\) from \(4\) to \(40\), we consider an \(F_{t-1}\times F_t\) array of cells indexed by two periodic coordinates. A symmetry may independently shift or reverse either coordinate, so two colorings are equivalent if one can be transformed into the other by those operations. Only colorings that use all \(t\) colors at least once are valid, and Problem 651 asks for the sum of all such orbit counts modulo \(10^9+7\).

If \(A_q(r,s)\) denotes the number of valid equivalence classes with \(q\) colors and side lengths \(r\) and \(s\), then the target quantity is

$\sum_{t=4}^{40} A_t(F_{t-1},F_t)\pmod{10^9+7}.$

A brute-force search over all colorings is hopeless, so the solution counts fixed colorings under symmetry classes and averages them with Burnside's lemma.

Mathematical Approach

The implementation solves a general instance first and only then applies it to consecutive Fibonacci sizes. The key object is the periodic grid

$X_{r,s}=\mathbb{Z}/r\mathbb{Z}\times \mathbb{Z}/s\mathbb{Z}.$

Step 1: Describe the Symmetries on One Axis

For one periodic axis of length \(n\), the allowed symmetries are the affine maps

$x\mapsto \sigma x+t \pmod n,\qquad \sigma\in\{1,-1\},\ t\in \mathbb{Z}/n\mathbb{Z}.$

When \(n>2\), this set has \(2n\) elements; for \(n=1\) and \(n=2\) some maps coincide, so those small cases are handled separately by the implementation.

The maps with \(\sigma=1\) are pure shifts. If such a shift has order \(u\), then every orbit on that axis is a cycle of length \(u\), and there are \(n/u\) such cycles. The number of shifts of order \(u\) is \(\varphi(u)\), because

$u=\frac{n}{\gcd(n,t)}$

and exactly \(\varphi(u)\) offsets produce that order. So every divisor \(u\mid n\) contributes one shift class with multiplicity \(\varphi(u)\) and cycle profile

$u^{\,n/u}.$

The maps with \(\sigma=-1\) are reversals. Their cycle structure depends only on the parity of \(n\):

$\begin{aligned} n\ \text{odd}&:\quad 1^1\,2^{(n-1)/2},\\ n\ \text{even, one half of the reversals}&:\quad 1^2\,2^{(n-2)/2},\\ n\ \text{even, the other half}&:\quad 2^{n/2}. \end{aligned}$

This follows from the congruence \(2x\equiv t\pmod n\): for odd \(n\) it has one solution, for even \(n\) it has either two or zero solutions depending on the parity of \(t\).

Step 2: Apply Burnside's Lemma

Let \(\Gamma_r\) and \(\Gamma_s\) be the symmetry sets on the two axes. A coloring is a function

$f:X_{r,s}\to \{1,2,\dots,q\}.$

For \((g,h)\in \Gamma_r\times \Gamma_s\), let \(c(g,h)\) be the number of cycles of the induced permutation on the cell set \(X_{r,s}\). A coloring fixed by \((g,h)\) must be constant on each cycle, so the fixed colorings that use every color are exactly the surjections from a \(c(g,h)\)-element cycle set onto the \(q\) colors.

By inclusion-exclusion, that number is

$\operatorname{Surj}(c,q)=\sum_{j=0}^{q}(-1)^j\binom{q}{j}(q-j)^c.$

Therefore Burnside gives

$A_q(r,s)=\frac{1}{|\Gamma_r|\,|\Gamma_s|}\sum_{g\in\Gamma_r}\sum_{h\in\Gamma_s}\operatorname{Surj}(c(g,h),q).$

This formula is exact, but it becomes fast only after we compress group elements by cycle type.

Step 3: Count Cycles on the Product Grid

Suppose one symmetry on the first axis contains \(m_i\) cycles of length \(\ell_i\), and one symmetry on the second axis contains \(n_j\) cycles of length \(r_j\). On the Cartesian product of one \(\ell_i\)-cycle and one \(r_j\)-cycle, the action advances simultaneously around both cycles.

The orbit length is \(\operatorname{lcm}(\ell_i,r_j)\), so the \(\ell_i r_j\) points split into

$\frac{\ell_i r_j}{\operatorname{lcm}(\ell_i,r_j)}=\gcd(\ell_i,r_j)$

product cycles. Summing over all profile pairs yields

$c(g,h)=\sum_{i}\sum_{j} m_i n_j\,\gcd(\ell_i,r_j).$

This is the central formula used to turn two one-dimensional cycle profiles into the number of cycles on the full two-dimensional pattern.

Step 4: Compress Burnside by Symmetry Classes

All group elements with the same cycle profile contribute the same fixed-coloring count, so the Burnside sum can be rewritten over classes instead of individual symmetries.

For one axis of length \(n\), the implementation needs only

$\tau(n)+O(1)$

different profiles: one shift profile for each divisor of \(n\), plus one or two reversal profiles according to parity. The multiplicities of those profiles are \(\varphi(u)\) for shifts and the obvious parity-dependent counts for reversals.

This compression is why huge symmetry groups never need to be enumerated explicitly. In the actual Fibonacci inputs, the number of class pairs is tiny compared with the raw number of group elements.

Step 5: Worked Example \(\boldsymbol{A_2(2,3)=11}\)

This checkpoint appears in the implementation and is small enough to compute by hand.

For length \(2\), the axis has two classes:

$1^2\ \text{with multiplicity }1,\qquad 2^1\ \text{with multiplicity }1.$

For length \(3\), the axis has three classes:

$1^3\ \text{with multiplicity }1,\qquad 3^1\ \text{with multiplicity }2,\qquad 1^1 2^1\ \text{with multiplicity }3.$

Because \(q=2\), the surjection count simplifies to

$\operatorname{Surj}(c,2)=2^c-2.$

Using the product-cycle formula, the six class pairs give

$\begin{aligned} c(1^2,1^3)&=6, & \operatorname{mult}&=1,\\ c(1^2,3^1)&=2, & \operatorname{mult}&=2,\\ c(1^2,1^1 2^1)&=4, & \operatorname{mult}&=3,\\ c(2^1,1^3)&=3, & \operatorname{mult}&=1,\\ c(2^1,3^1)&=1, & \operatorname{mult}&=2,\\ c(2^1,1^1 2^1)&=3, & \operatorname{mult}&=3. \end{aligned}$

So the weighted Burnside numerator is

$1\cdot 62+2\cdot 2+3\cdot 14+1\cdot 6+2\cdot 0+3\cdot 6=132.$

The group sizes are \(|\Gamma_2|=2\) and \(|\Gamma_3|=6\), hence

$A_2(2,3)=\frac{132}{2\cdot 6}=11.$

Step 6: Apply the General Formula to Consecutive Fibonacci Sizes

Once the class-compressed Burnside count \(A_q(r,s)\) is available, Problem 651 simply evaluates it at

$q=t,\qquad r=F_{t-1},\qquad s=F_t,\qquad 4\le t\le 40,$

and accumulates the results modulo \(10^9+7\):

$\boxed{\sum_{t=4}^{40} A_t(F_{t-1},F_t)\pmod{10^9+7}.}$

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical pipeline. First they build the binomial coefficients needed by the inclusion-exclusion formula and use fast modular exponentiation for the powers \((q-j)^c\).

For each of the two periodic lengths, the implementation factorizes the length, enumerates its divisors, converts divisor data into shift classes with multiplicities \(\varphi(u)\), and then appends the parity-dependent reversal classes. The tiny degenerate lengths \(1\) and \(2\) are treated explicitly.

Next the implementation loops over every pair of classes from the two axes. For each pair it computes the total cycle count by summing \(\gcd(\ell_i,r_j)\) across the profile entries. Many different class pairs lead to the same cycle count, so the value of \(\operatorname{Surj}(c,q)\) is cached and reused.

Each Burnside contribution is the class multiplicity product times the cached surjection count. After summing all contributions, the implementation divides by the group size using a modular inverse and adds the resulting orbit count to the running Fibonacci total.

Complexity Analysis

For one general instance \(A_q(r,s)\), factorizing the two side lengths by trial division costs \(O(\sqrt{r}+\sqrt{s})\) time in the implementations shown here. If \(K(r)\) and \(K(s)\) denote the numbers of symmetry classes on the two axes, then the Burnside loop uses \(O(K(r)K(s))\) class pairs.

Each class profile has only one or two cycle lengths, so the product-cycle computation for one pair is \(O(1)\). Every distinct cycle count triggers one inclusion-exclusion evaluation in \(O(q)\), but \(q\le 40\) in Problem 651 and many class pairs share the same cycle count. Memory usage is \(O(K(r)+K(s)+U+q^2)\), where \(U\) is the number of distinct cached cycle counts. In practice the method is extremely small compared with the combinatorial size of the raw coloring space.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=651
2. Burnside's lemma: Wikipedia — Burnside's lemma
3. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
4. Euler's totient function: Wikipedia — Euler's totient function
5. Dihedral group and cycle symmetries: Wikipedia — Dihedral group
6. Stirling numbers of the second kind: Wikipedia — Stirling numbers of the second kind