Problem 728: Circle of Coins
View on Project EulerProject Euler Problem 728 Solution
EulerSolve provides an optimized solution for Project Euler Problem 728, Circle of Coins, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each pair \((n,k)\) with \(1 \le k \le n\), the circle-of-coins problem has a counting function \(F(n,k)\) for the admissible binary coin configurations associated with that pair. The overall task is to evaluate $$S(N)=\sum_{n=1}^{N}\sum_{k=1}^{n} F(n,k)\pmod{10^9+7}.$$ The key point in the implementations is that the combinatorics collapse to arithmetic data: only \(\gcd(n,k)\), the exponent of \(2\) dividing \(n\) and \(k\), and a totient-based regrouping are needed. That removes any need to enumerate coin states directly. Mathematical Approach The C++, Python, and Java implementations all use the same arithmetic decomposition. Once the fixed-pair formula is known, the remaining work is to reorganize the double sum so that equal contributions are collected together efficiently. Step 1: Closed Form for a Fixed Pair Let $$g=\gcd(n,k),$$ and let \(v_2(x)\) denote the exponent of \(2\) in \(x\). The implementations use the closed form $$F(n,k)=\begin{cases} 2^{n-g+1}, & v_2(k)\le v_2(n),\\ 2^{n-g}, & v_2(k)>v_2(n). \end{cases}$$ So the entire dependence on the original circle is compressed into two quantities: the common divisor \(g\), which controls the main exponent \(n-g\), and a single parity-sensitive comparison of 2-adic valuations, which decides whether there is one extra factor of \(2\)....
Detailed mathematical approach
Problem Summary
For each pair \((n,k)\) with \(1 \le k \le n\), the circle-of-coins problem has a counting function \(F(n,k)\) for the admissible binary coin configurations associated with that pair. The overall task is to evaluate
$$S(N)=\sum_{n=1}^{N}\sum_{k=1}^{n} F(n,k)\pmod{10^9+7}.$$
The key point in the implementations is that the combinatorics collapse to arithmetic data: only \(\gcd(n,k)\), the exponent of \(2\) dividing \(n\) and \(k\), and a totient-based regrouping are needed. That removes any need to enumerate coin states directly.
Mathematical Approach
The C++, Python, and Java implementations all use the same arithmetic decomposition. Once the fixed-pair formula is known, the remaining work is to reorganize the double sum so that equal contributions are collected together efficiently.
Step 1: Closed Form for a Fixed Pair
Let
$$g=\gcd(n,k),$$
and let \(v_2(x)\) denote the exponent of \(2\) in \(x\). The implementations use the closed form
$$F(n,k)=\begin{cases} 2^{n-g+1}, & v_2(k)\le v_2(n),\\ 2^{n-g}, & v_2(k)>v_2(n). \end{cases}$$
So the entire dependence on the original circle is compressed into two quantities: the common divisor \(g\), which controls the main exponent \(n-g\), and a single parity-sensitive comparison of 2-adic valuations, which decides whether there is one extra factor of \(2\).
Step 2: Separate the GCD from the Reduced Step
Write
$$n=t h,\qquad k=t r,\qquad t=\gcd(n,k),\qquad \gcd(r,h)=1.$$
Then \(t\) is exactly the gcd that appears in the fixed-pair formula, so
$$n-g=t h-t=t(h-1).$$
After removing the common factor \(t\), every pair \((n,k)\) is described by a reduced denominator \(h\) and a reduced residue \(r\) coprime to \(h\). The formula becomes
$$F(t h,t r)=\begin{cases} 2^{t(h-1)+1}, & v_2(r)\le v_2(h),\\ 2^{t(h-1)}, & v_2(r)>v_2(h). \end{cases}$$
This shows that, for fixed \(h\) and \(t\), all dependence on the reduced step is packed into how many coprime residues \(r\) satisfy the valuation test.
Step 3: Count the Reduced Residues for Fixed \(h\)
Now fix \(h \ge 2\) and sum over all \(r\) with \(1 \le r \le h\) and \(\gcd(r,h)=1\).
If \(h\) is even, every residue coprime to \(h\) must be odd. Hence \(v_2(r)=0\le v_2(h)\), so every reduced residue contributes the larger value \(2^{t(h-1)+1}\). Since there are \(\varphi(h)\) such residues, their total contribution is
$$\varphi(h)\cdot 2^{t(h-1)+1}=2\varphi(h)\,2^{t(h-1)}.$$
If \(h\) is odd, the reduced residues split evenly into odd and even values. Indeed, for every reduced residue \(r\), the paired residue \(h-r\) is also reduced and has opposite parity. Therefore exactly half of the \(\varphi(h)\) residues satisfy \(v_2(r)=0=v_2(h)\), while the other half fail the inequality.
Step 4: Derive the Coefficient \(c(h)\)
From the parity split above, the total contribution for one fixed \(h \ge 2\) and one fixed \(t\) is
$$c(h)\,2^{t(h-1)},$$
where
$$c(h)=\begin{cases} 2\varphi(h), & h \text{ even},\\ \dfrac{3\varphi(h)}{2}, & h \text{ odd}. \end{cases}$$
The odd case is just
$$\frac{\varphi(h)}{2}\cdot 2^{t(h-1)+1}+\frac{\varphi(h)}{2}\cdot 2^{t(h-1)} =\frac{3\varphi(h)}{2}\,2^{t(h-1)}.$$
This is the crucial compression step: once pairs are grouped by \(h\), all the messy dependence on \(k\) is replaced by a single totient-based coefficient.
Step 5: Rebuild the Whole Sum
The case \(h=1\) is special. Then \(k=n\), so \(g=n\) and the fixed-pair formula gives \(F(n,n)=2\). Summed over \(n=1,\dots,N\), this contributes
$$2N.$$
For every \(h \ge 2\), the remaining parameter is \(t\), and the condition \(n=t h \le N\) means
$$1\le t\le \left\lfloor\frac{N}{h}\right\rfloor.$$
Therefore the full sum becomes
$$\boxed{S(N)=2N+\sum_{h=2}^{N} c(h)\sum_{t=1}^{\lfloor N/h\rfloor} 2^{t(h-1)} \pmod{10^9+7}.}$$
Mathematically the inner sum is a finite geometric progression, but the implementations simply step through its exponents directly after precomputing powers of \(2\).
Worked Example: \(N=3\)
The special term is
$$2N=6.$$
For \(h=2\), we have \(\varphi(2)=1\) and \(c(2)=2\). Also \(\lfloor 3/2\rfloor=1\), so this block contributes
$$2\cdot 2^{1}=4.$$
For \(h=3\), we have \(\varphi(3)=2\) and \(c(3)=3\). Also \(\lfloor 3/3\rfloor=1\), so this block contributes
$$3\cdot 2^{2}=12.$$
Hence
$$S(3)=6+4+12=22,$$
which matches the checkpoint used by the implementation. A single-pair checkpoint is \(F(9,3)=2^{9-3+1}=128\), because \(\gcd(9,3)=3\) and \(v_2(3)\le v_2(9)\).
How the Code Works
The implementations first build Euler's totient values for all integers up to \(N\) with a linear sieve. They also precompute the table
$$2^0,2^1,\dots,2^N \pmod{10^9+7},$$
so every power lookup in the main sum is constant time.
After that, the algorithm starts with the special contribution \(2N\). It then iterates through \(h=2,3,\dots,N\), computes the coefficient \(c(h)\) from the parity of \(h\) and the totient value, and adds
$$c(h)\,2^{t(h-1)}$$
for each \(t=1,\dots,\lfloor N/h\rfloor\). All arithmetic is reduced modulo \(10^9+7\) after each addition. The C++, Python, and Java implementations follow exactly the same mathematical plan; the C++ version also includes small checkpoint assertions before producing the final large-input result.
Complexity Analysis
The totient sieve runs in \(O(N)\) time and uses \(O(N)\) memory. The double summation over blocks costs
$$\sum_{h=2}^{N}\left\lfloor\frac{N}{h}\right\rfloor=O(N\log N).$$
Therefore the total running time is \(O(N\log N)\), while the memory usage remains \(O(N)\).
Footnotes and References
- Problem page: Project Euler 728 — Circle of Coins
- Greatest common divisor: Wikipedia — Greatest common divisor
- Euler's totient function: Wikipedia — Euler's totient function
- \(p\)-adic valuation, including the special case \(v_2\): Wikipedia — \(p\)-adic valuation
- Geometric progression: Wikipedia — Geometric progression
Problem 728 source code
C++
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <vector>
namespace {
using i64 = std::int64_t;
using u64 = std::uint64_t;
constexpr i64 kMod = 1'000'000'007LL;
int v2(const int x) {
return __builtin_ctz(static_cast<unsigned>(x));
}
i64 mod_pow2(i64 exp) {
i64 base = 2;
i64 result = 1;
while (exp > 0) {
if (exp & 1LL) {
result = static_cast<i64>((__int128)result * base % kMod);
}
base = static_cast<i64>((__int128)base * base % kMod);
exp >>= 1LL;
}
return result;
}
i64 F(const int n, const int k) {
const int g = std::gcd(n, k);
if (v2(k) <= v2(n)) {
return mod_pow2(static_cast<i64>(n - g + 1));
}
return mod_pow2(static_cast<i64>(n - g));
}
i64 solve(const int N) {
std::vector<int> phi(static_cast<std::size_t>(N + 1), 0);
std::vector<int> primes;
primes.reserve(static_cast<std::size_t>(N / 10));
phi[1] = 1;
for (int i = 2; i <= N; ++i) {
if (phi[static_cast<std::size_t>(i)] == 0) {
phi[static_cast<std::size_t>(i)] = i - 1;
primes.push_back(i);
}
for (const int p : primes) {
const i64 v = static_cast<i64>(i) * p;
if (v > N) {
break;
}
if (i % p == 0) {
phi[static_cast<std::size_t>(v)] = phi[static_cast<std::size_t>(i)] * p;
break;
}
phi[static_cast<std::size_t>(v)] = phi[static_cast<std::size_t>(i)] * (p - 1);
}
}
std::vector<i64> pow2(static_cast<std::size_t>(N + 1), 1);
for (int i = 1; i <= N; ++i) {
pow2[static_cast<std::size_t>(i)] = (pow2[static_cast<std::size_t>(i - 1)] * 2) % kMod;
}
i64 ans = (2LL * N) % kMod; // h = 1 contribution for each n
for (int h = 2; h <= N; ++h) {
i64 coeff;
if ((h & 1) == 0) {
coeff = (2LL * phi[static_cast<std::size_t>(h)]) % kMod;
} else {
coeff = (3LL * (phi[static_cast<std::size_t>(h)] / 2LL)) % kMod;
}
const int step = h - 1;
const int count = N / h;
int exp = step;
for (int i = 0; i < count; ++i) {
const i64 term = static_cast<i64>((__int128)coeff * pow2[static_cast<std::size_t>(exp)] % kMod);
ans += term;
if (ans >= kMod) {
ans -= kMod;
}
exp += step;
}
}
return ans;
}
} // namespace
int main() {
assert(F(3, 2) == 4);
assert(F(8, 3) == 256);
assert(F(9, 3) == 128);
assert(solve(3) == 22);
assert(solve(10) == 10'444);
assert(solve(1'000) == 853'837'042);
std::cout << solve(10'000'000) << '\n';
return 0;
}
Python
import math
def solve():
MOD = 1000000007
N = 10000000
def mod_pow(base, exp):
r = 1; base %= MOD
while exp > 0:
if exp & 1: r = r * base % MOD
base = base * base % MOD; exp >>= 1
return r
phi = list(range(N + 1))
primes = []
for i in range(2, N + 1):
if phi[i] == i: phi[i] = i - 1; primes.append(i)
for p in primes:
if i * p > N: break
if i % p == 0: phi[i * p] = phi[i] * p; break
phi[i * p] = phi[i] * (p - 1)
pow2 = [1] * (N + 1)
for i in range(1, N + 1): pow2[i] = pow2[i-1] * 2 % MOD
ans = 2 * N % MOD # h=1
def v2(x):
c = 0
while x & 1 == 0: x >>= 1; c += 1
return c
for h in range(2, N + 1):
if h & 1:
coeff = 3 * (phi[h] // 2) % MOD
else:
coeff = 2 * phi[h] % MOD
step = h - 1; count = N // h
exp = step
for i in range(count):
ans = (ans + coeff * pow2[exp]) % MOD
exp += step
return str(ans)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.List;
public class Euler728 {
static final long kMod = 1000000007L;
public static String solve() {
int N = 10000000;
int[] phi = new int[N + 1];
List<Integer> primes = new ArrayList<>(N / 10);
phi[1] = 1;
for (int i = 2; i <= N; ++i) {
if (phi[i] == 0) {
phi[i] = i - 1;
primes.add(i);
}
for (int p : primes) {
long v = (long) i * p;
if (v > N)
break;
if (i % p == 0) {
phi[(int) v] = phi[i] * p;
break;
}
phi[(int) v] = phi[i] * (p - 1);
}
}
int[] pow2 = new int[N + 1];
pow2[0] = 1;
for (int i = 1; i <= N; ++i) {
pow2[i] = (int) (((long) pow2[i - 1] * 2) % kMod);
}
long ans = (2L * N) % kMod;
for (int h = 2; h <= N; ++h) {
long coeff;
if ((h & 1) == 0) {
coeff = (2L * phi[h]) % kMod;
} else {
coeff = (3L * (phi[h] / 2L)) % kMod;
}
int step = h - 1;
int count = N / h;
int exp = step;
for (int i = 0; i < count; ++i) {
long term = (coeff * pow2[exp]) % kMod;
ans += term;
if (ans >= kMod) {
ans -= kMod;
}
exp += step;
}
}
return Long.toString(ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}