Detailed mathematical approach

Problem Summary

For each ordered coprime pair \((a,b)\) with \(2 \le a,b \le n\), the corresponding Lissajous-curve contribution collapses to a rational expression determined by \(a\), \(b\), and by how the factors \(2\) and \(5\) of \(10\) are distributed between the two frequencies. The goal is to compute the total

$s(n)=\sum_{\substack{2 \le a,b \le n\\ \gcd(a,b)=1}} w(a,b)$

exactly for very large \(n\). The implementations therefore keep the integer numerator \(4s(n)\) throughout the calculation and divide by \(4\) only at the end.

Mathematical Approach

Write

$C_n=\{(a,b):2 \le a,b \le n,\ \gcd(a,b)=1\},\qquad T(k)=\frac{k(k+1)}{2}.$

The full answer is built from two global coprime sums over \(C_n\), plus a correction on a smaller exceptional set determined only by divisibility with respect to \(10\).

Step 1: Rewrite Coprimality with Möbius Inversion

The starting point is

$\mathbf 1_{\gcd(a,b)=1}=\sum_{d \mid a,\ d \mid b}\mu(d)=\sum_{d \mid \gcd(a,b)}\mu(d).$

This converts coprime pair sums into divisor sums. Over the full square \(1 \le a,b \le n\), the two basic aggregates become

$S_{ab}^{\ast}(n)=\sum_{\substack{1 \le a,b \le n\\ \gcd(a,b)=1}}ab=\sum_{d=1}^{n}\mu(d)d^2T\left(\left\lfloor\frac{n}{d}\right\rfloor\right)^2,$

$S_{a}^{\ast}(n)=\sum_{\substack{1 \le a,b \le n\\ \gcd(a,b)=1}}a=\sum_{d=1}^{n}\mu(d)d\,T\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\left\lfloor\frac{n}{d}\right\rfloor.$

These are exactly the two large Möbius sums computed first.

Step 2: Remove the Boundary Cases \(a=1\) or \(b=1\)

The target sum starts at \(2\), not \(1\), so the boundary contributions must be subtracted. Since

$\sum_{m=1}^{n} m = T(n),$

we obtain

$S_{ab}(n)=S_{ab}^{\ast}(n)-2T(n)+1,$

$S_{a}(n)=S_{a}^{\ast}(n)-n-T(n)+1.$

Hence

$S_{ab}(n)=\sum_{(a,b)\in C_n}ab,\qquad S_{a}(n)=\sum_{(a,b)\in C_n}a.$

Because \(C_n\) is symmetric under swapping the coordinates, the same \(S_a(n)\) also equals \(\sum_{(a,b)\in C_n} b\).

Step 3: Isolate the Base-10 Exceptional Set

Define the bucket map

$g(m)=\gcd(m,10)\in\{1,2,5,10\}.$

The correction layer depends only on whether one frequency contributes the factor \(2\) of \(10\) and the other contributes the factor \(5\). That condition is

$E_n=\{(a,b)\in C_n:g(a)g(b)=10\}.$

The possible bucket matches are therefore

$1 \leftrightarrow 10,\qquad 2 \leftrightarrow 5.$

So the entire non-generic part of the problem reduces to counting coprime pairs in four complementary residue-type buckets.

Step 4: Count the Exceptional Pairs Without Enumerating Them

For any subset \(U \subseteq \{1,\dots,n\}\), define for fixed \(a\)

$N_U(a)=\sum_{d \mid a}\mu(d)\,\#\{b \in U:d \mid b\},$

$W_U(a)=\sum_{d \mid a}\mu(d)\sum_{\substack{b \in U\\ d \mid b}} b.$

These give the count and weighted sum of \(b\)-values in \(U\) that are coprime to \(a\). The four needed subsets are \(g(b)=10\), \(g(b)=5\), \(g(b)=2\), and \(g(b)=1\), and each one has a closed-form arithmetic progression formula.

If \(g(b)=10\), then \(b\) is a multiple of \(10\). With \(L=\operatorname{lcm}(d,10)\),

$C_{10}(d)=\left\lfloor\frac{n}{L}\right\rfloor,\qquad Q_{10}(d)=L\,T\left(\left\lfloor\frac{n}{L}\right\rfloor\right).$

If \(g(b)=5\), then \(b\) is an odd multiple of \(5\). When \(2 \mid d\) there are no such multiples. Otherwise, with \(L=\operatorname{lcm}(d,5)\), \(t=\left\lfloor n/L \right\rfloor\), and \(k=\left\lceil t/2 \right\rceil\),

$C_{5}(d)=k,\qquad Q_{5}(d)=Lk^2,$

because \(1+3+\cdots+(2k-1)=k^2\).

If \(g(b)=2\), then \(b\) is even but not divisible by \(5\). When \(5 \mid d\) the count is zero. Otherwise, with \(L=\operatorname{lcm}(d,2)\), \(t=\left\lfloor n/L \right\rfloor\), and \(t_5=\left\lfloor n/(5L) \right\rfloor\),

$C_{2}(d)=t-t_5,\qquad Q_{2}(d)=L\,T(t)-5L\,T(t_5).$

If \(g(b)=1\), then \(b\) is coprime to \(10\). When \(\gcd(d,10)\ne 1\) the count is zero. Otherwise, with \(x=\left\lfloor n/d \right\rfloor\), inclusion-exclusion on divisibility by \(2\) and \(5\) gives

$C_{1}(d)=x-\left\lfloor\frac{x}{2}\right\rfloor-\left\lfloor\frac{x}{5}\right\rfloor+\left\lfloor\frac{x}{10}\right\rfloor,$

$Q_{1}(d)=d\left(T(x)-2T\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-5T\left(\left\lfloor\frac{x}{5}\right\rfloor\right)+10T\left(\left\lfloor\frac{x}{10}\right\rfloor\right)\right).$

For \(g(a)=10\), the value \(b=1\) must be removed afterward, because the global sum begins at \(2\).

Step 5: Assemble the Exact Total

Summing the exceptional counts and weighted sums over all \(a\) gives

$N_{10}(n)=|E_n|,\qquad S_{a,10}(n)=\sum_{(a,b)\in E_n} a,\qquad S_{ab,10}(n)=\sum_{(a,b)\in E_n} ab.$

The implementations then combine everything as

$4s(n)=8S_{ab}(n)-12S_a(n)+6S_{a,10}(n)+4N_{10}(n)-6S_{ab,10}(n).$

Equivalently,

$s(n)=2S_{ab}(n)-3S_a(n)+\frac{6S_{a,10}(n)+4N_{10}(n)-6S_{ab,10}(n)}{4}.$

This can also be read pairwise. Every \((a,b)\in C_n\) contributes

$w(a,b)= \begin{cases} 2ab-\dfrac{3}{2}(a+b), & g(a)g(b)\ne 10,\\[6pt] \dfrac{2ab-3a-3b+4}{4}, & g(a)g(b)=10. \end{cases}$

So the total is exactly

$s(n)=\sum_{(a,b)\in C_n} w(a,b).$

The quarter-factor in the exceptional case is precisely why the programs keep \(4s(n)\) as an integer.

Worked Example: \(n=10\)

For \(n=10\), the Möbius sums give

$S_{ab}(10)=1554,\qquad S_a(10)=261.$

The exceptional set \(E_{10}\) contains \(14\) ordered pairs, including \((2,5)\), \((3,10)\), \((5,8)\), \((10,3)\), and \((10,9)\). Its totals are

$N_{10}(10)=14,\qquad S_{a,10}(10)=89,\qquad S_{ab,10}(10)=580.$

Therefore

$\begin{aligned} 4s(10) &=8\cdot 1554-12\cdot 261+6\cdot 89+4\cdot 14-6\cdot 580\\ &=12432-3132+534+56-3480\\ &=6410, \end{aligned}$

so

$s(10)=\frac{6410}{4}=1602.5.$

This matches the checkpoint built into the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical decomposition. First they build a Möbius table up to \(n\). Next they evaluate the two large divisor sums for \(S_{ab}^{\ast}(n)\) and \(S_a^{\ast}(n)\), then apply the boundary subtraction that removes the cases with \(a=1\) or \(b=1\).

After that, the implementation processes the four base-10 buckets. Instead of factoring every \(a\) separately and enumerating its divisors, it loops over each divisor \(d\), computes the relevant count and weighted-sum formulas once, and distributes the Möbius-weighted contribution to every multiple of \(d\). This produces per-\(a\) tables for the exceptional set in near-harmonic time. The C++ version can run the four independent bucket passes in parallel; the Python and Java versions preserve the same arithmetic structure and final exact value.

Complexity Analysis

The Möbius sieve is linear, so it costs \(O(n)\) time and \(O(n)\) memory. The dominant work is the divisor-to-multiples sweep used by the four bucket tables, whose total size is controlled by

$\sum_{d=1}^{n}\left\lfloor\frac{n}{d}\right\rfloor = O(n\log n).$

Therefore the overall algorithm runs in \(O(n\log n)\) time and uses \(O(n)\) memory. Parallel execution improves wall-clock time in the C++ implementation but does not change the asymptotic bound.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=777
2. Lissajous curves: Wikipedia - Lissajous curve
3. Mobius inversion: Wikipedia - Mobius inversion formula
4. Mobius function: Wikipedia - Mobius function
5. Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
6. Coprime integers: Wikipedia - Coprime integers