Problem 643: $2$-Friendly
View on Project EulerProject Euler Problem 643 Solution
EulerSolve provides an optimized solution for Project Euler Problem 643, $2$-Friendly, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary A pair \((a,b)\) with \(1 \le a \lt b \le N\) is called 2-friendly when \(\gcd(a,b)\) is a power of two strictly larger than \(1\). If \(f(N)\) denotes the number of such pairs, the task is to compute $$f(10^{11}) \pmod{10^9+7}.$$ A direct scan over all pairs would be quadratic in \(N\), so the solution has to reorganize the counting problem around the exact gcd. Mathematical Approach We want to count all pairs whose gcd has the form \(2^t\) with \(t \ge 1\). The decisive observation is that once the exact gcd is fixed, what remains is a coprime pair-counting problem. Step 1: Separate the Exact Power-of-Two gcd Suppose \(\gcd(a,b)=2^t\) for some \(t \ge 1\). Then we can write $$a=2^t x,\qquad b=2^t y,$$ with $$1 \le x \lt y \le \left\lfloor \frac{N}{2^t} \right\rfloor.$$ Because the gcd was assumed to be exactly \(2^t\), the reduced pair must satisfy $$\gcd(x,y)=1.$$ Conversely, every coprime pair \((x,y)\) in that range produces exactly one original pair \((2^t x,2^t y)\) whose gcd is exactly \(2^t\). So there is no overcounting between different values of \(t\). Step 2: Count Coprime Pairs with Euler's Totient Function For a fixed upper bound \(m\), define $$C(m)=\#\{(x,y):1 \le x \lt y \le m,\ \gcd(x,y)=1\}.$$ If we fix the larger coordinate \(y\), then the valid values of \(x\) are precisely the integers \(1 \le x \lt y\) that are coprime to \(y\)....
Detailed mathematical approach
Problem Summary
A pair \((a,b)\) with \(1 \le a \lt b \le N\) is called 2-friendly when \(\gcd(a,b)\) is a power of two strictly larger than \(1\). If \(f(N)\) denotes the number of such pairs, the task is to compute
$$f(10^{11}) \pmod{10^9+7}.$$
A direct scan over all pairs would be quadratic in \(N\), so the solution has to reorganize the counting problem around the exact gcd.
Mathematical Approach
We want to count all pairs whose gcd has the form \(2^t\) with \(t \ge 1\). The decisive observation is that once the exact gcd is fixed, what remains is a coprime pair-counting problem.
Step 1: Separate the Exact Power-of-Two gcd
Suppose \(\gcd(a,b)=2^t\) for some \(t \ge 1\). Then we can write
$$a=2^t x,\qquad b=2^t y,$$
with
$$1 \le x \lt y \le \left\lfloor \frac{N}{2^t} \right\rfloor.$$
Because the gcd was assumed to be exactly \(2^t\), the reduced pair must satisfy
$$\gcd(x,y)=1.$$
Conversely, every coprime pair \((x,y)\) in that range produces exactly one original pair \((2^t x,2^t y)\) whose gcd is exactly \(2^t\). So there is no overcounting between different values of \(t\).
Step 2: Count Coprime Pairs with Euler's Totient Function
For a fixed upper bound \(m\), define
$$C(m)=\#\{(x,y):1 \le x \lt y \le m,\ \gcd(x,y)=1\}.$$
If we fix the larger coordinate \(y\), then the valid values of \(x\) are precisely the integers \(1 \le x \lt y\) that are coprime to \(y\). Their number is \(\varphi(y)\). Therefore
$$C(m)=\sum_{y=2}^{m}\varphi(y).$$
Now introduce the summatory totient function
$$\Phi(m)=\sum_{k=1}^{m}\varphi(k).$$
Since \(\varphi(1)=1\), we get the compact identity
$$C(m)=\Phi(m)-1.$$
Step 3: Sum over All Powers of Two
For a fixed \(t\), the admissible pairs with gcd \(2^t\) are counted by
$$C\left(\left\lfloor \frac{N}{2^t} \right\rfloor\right)=\Phi\left(\left\lfloor \frac{N}{2^t} \right\rfloor\right)-1.$$
Hence
$$f(N)=\sum_{t \ge 1}\left(\Phi\left(\left\lfloor \frac{N}{2^t} \right\rfloor\right)-1\right).$$
Only finitely many terms are nonzero, because once \(\left\lfloor N/2^t \right\rfloor \le 1\), there is no room for a pair with \(x \lt y\).
Step 4: Derive a Fast Recurrence for \(\Phi(n)\)
The implementations do not sum \(\varphi(k)\) up to \(n\) from scratch for every query. Instead they use the classical divisor identity
$$\sum_{d \mid m}\varphi(d)=m.$$
Summing this from \(m=1\) to \(m=n\) gives
$$\sum_{m=1}^{n} m=\sum_{m=1}^{n}\sum_{d \mid m}\varphi(d).$$
After exchanging the order of summation, this becomes
$$\frac{n(n+1)}{2}=\sum_{q=1}^{n}\Phi\left(\left\lfloor \frac{n}{q} \right\rfloor\right).$$
Separating the \(q=1\) term yields the recurrence
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{q=2}^{n}\Phi\left(\left\lfloor \frac{n}{q} \right\rfloor\right).$$
The floor value \(\left\lfloor n/q \right\rfloor\) is constant on intervals, so the sum is grouped by quotient blocks rather than processed one index at a time.
Step 5: Use Quotient Grouping and Memoization
If a block starts at \(l\), let
$$v=\left\lfloor \frac{n}{l} \right\rfloor,\qquad r=\left\lfloor \frac{n}{v} \right\rfloor.$$
Then every \(q\) in the interval \(l \le q \le r\) has the same quotient \(v\), so their total contribution is
$$ (r-l+1)\,\Phi(v). $$
This reduces the amount of work dramatically. Memoization then stores each large \(\Phi(v)\) value after its first computation, which is important because the outer formula asks for \(\Phi\) at the related arguments \(N/2,N/4,N/8,\dots\).
Worked Example: \(N=10\)
For \(N=10\), the relevant halvings are
$$\left\lfloor \frac{10}{2} \right\rfloor=5,\qquad \left\lfloor \frac{10}{4} \right\rfloor=2,\qquad \left\lfloor \frac{10}{8} \right\rfloor=1.$$
The last term contributes nothing because \(\Phi(1)-1=0\). So
$$f(10)=\left(\Phi(5)-1\right)+\left(\Phi(2)-1\right).$$
Now
$$\Phi(5)=1+1+2+2+4=10,\qquad \Phi(2)=1+1=2,$$
hence
$$f(10)=9+1=10.$$
The ten pairs are
$$(2,4),(2,6),(2,8),(2,10),(4,6),(4,10),(6,8),(6,10),(8,10),(4,8).$$
The first nine have gcd \(2\), and the last one has gcd \(4\).
How the Code Works
The C++, Python, and Java implementations follow the same structure. They first precompute \(\varphi(n)\) up to a fixed cutoff of five million with a linear sieve and store the prefix sums modulo \(10^9+7\). That makes every small \(\Phi(n)\) query an \(O(1)\) table lookup.
For larger \(n\), the implementation evaluates the recurrence
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{q=2}^{n}\Phi\left(\left\lfloor \frac{n}{q} \right\rfloor\right)$$
using quotient grouping, so each interval of equal floor value is handled in one step. The result is memoized, which prevents the same large summatory-totient value from being recomputed.
Finally, the main solve loop starts with \(m=\lfloor N/2 \rfloor\), repeatedly halves \(m\), and accumulates \(\Phi(m)-1\) until \(m \lt 2\). Every arithmetic step is reduced modulo \(10^9+7\).
Complexity Analysis
Let \(B=5{,}000{,}000\) be the precomputation cutoff used by the implementations. The linear sieve and prefix table take \(O(B)\) time and \(O(B)\) memory. The outer summation over \(N/2,N/4,N/8,\dots\) has only \(O(\log N)\) terms.
The expensive part is evaluating large \(\Phi(n)\) values, but quotient grouping compresses each summation into blocks of equal floor quotient, and memoization ensures repeated subproblems are solved once. That hybrid strategy is what makes \(N=10^{11}\) practical.
Footnotes and References
- Problem page: https://projecteuler.net/problem=643
- Euler's totient function: Wikipedia — Euler's totient function
- Coprime integers: Wikipedia — Coprime integers
- Greatest common divisor: Wikipedia — Greatest common divisor
- Dirichlet hyperbola method: Wikipedia — Dirichlet hyperbola method
Problem 643 source code
C++
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <unordered_map>
#include <vector>
namespace {
using u32 = std::uint32_t;
using u64 = std::uint64_t;
constexpr u64 kMod = 1'000'000'007ULL;
constexpr u64 kInv2 = 500'000'004ULL;
struct PhiSummatory {
int limit = 0;
std::vector<int> primes;
std::vector<int> lp;
std::vector<u32> phi;
std::vector<u64> pref;
std::unordered_map<u64, u64> memo;
explicit PhiSummatory(const int n) : limit(n), lp(n + 1, 0), phi(n + 1, 0), pref(n + 1, 0) {
primes.reserve(static_cast<std::size_t>(n / 10));
phi[1] = 1;
for (int i = 2; i <= n; ++i) {
if (lp[i] == 0) {
lp[i] = i;
primes.push_back(i);
phi[i] = static_cast<u32>(i - 1);
}
for (int p : primes) {
if (p > lp[i] || (u64)i * (u64)p > (u64)n) break;
lp[i * p] = p;
if (p == lp[i]) {
phi[i * p] = phi[i] * static_cast<u32>(p);
break;
}
phi[i * p] = phi[i] * static_cast<u32>(p - 1);
}
}
for (int i = 1; i <= n; ++i) {
pref[i] = (pref[i - 1] + static_cast<u64>(phi[i])) % kMod;
}
memo.reserve(1 << 20);
}
u64 sum_phi(const u64 n) {
if (n <= static_cast<u64>(limit)) return pref[static_cast<std::size_t>(n)];
const auto it = memo.find(n);
if (it != memo.end()) return it->second;
const u64 nn = n % kMod;
u64 res = nn * ((n + 1) % kMod) % kMod;
res = res * kInv2 % kMod;
for (u64 l = 2; l <= n;) {
const u64 q = n / l;
const u64 r = n / q;
const u64 cnt = (r - l + 1) % kMod;
res = (res + kMod - cnt * sum_phi(q) % kMod) % kMod;
l = r + 1;
}
memo.emplace(n, res);
return res;
}
};
u64 solve(const u64 n, PhiSummatory& ph) {
u64 ans = 0;
for (u64 m = n / 2; m >= 2; m /= 2) {
const u64 add = (ph.sum_phi(m) + kMod - 1) % kMod;
ans += add;
ans %= kMod;
}
return ans;
}
u64 brute(const u32 n) {
u64 cnt = 0;
for (u32 a = 1; a <= n; ++a) {
for (u32 b = a + 1; b <= n; ++b) {
const u32 g = std::gcd(a, b);
if (g > 1 && (g & (g - 1)) == 0) ++cnt;
}
}
return cnt % kMod;
}
} // namespace
int main() {
PhiSummatory ph(5'000'000);
assert(solve(100, ph) == 1031);
assert(solve(1'000'000, ph) == 321'418'433ULL);
assert(solve(2000, ph) == brute(2000));
std::cout << solve(100'000'000'000ULL, ph) << "\n";
return 0;
}
Python
import math
import sys
sys.setrecursionlimit(2000)
MOD = 1000000007
INV2 = 500000004
class PhiSummatory:
def __init__(self, limit):
self.limit = limit
self.primes = []
lp = [0] * (limit + 1)
self.phi = [0] * (limit + 1)
self.pref = [0] * (limit + 1)
self.phi[1] = 1
for i in range(2, limit + 1):
if lp[i] == 0:
lp[i] = i
self.primes.append(i)
self.phi[i] = i - 1
for p in self.primes:
if p > lp[i] or i * p > limit: break
lp[i * p] = p
if p == lp[i]:
self.phi[i * p] = self.phi[i] * p
break
self.phi[i * p] = self.phi[i] * (p - 1)
for i in range(1, limit + 1):
self.pref[i] = (self.pref[i - 1] + self.phi[i]) % MOD
self.memo = {}
def sum_phi(self, n):
if n <= self.limit: return self.pref[n]
if n in self.memo: return self.memo[n]
nn = n % MOD
res = (nn * ((n + 1) % MOD)) % MOD
res = (res * INV2) % MOD
l = 2
while l <= n:
q = n // l
r = n // q
cnt = (r - l + 1) % MOD
res = (res - cnt * self.sum_phi(q)) % MOD
l = r + 1
res = (res + MOD) % MOD
self.memo[n] = res
return res
def solve_impl(n, ph):
ans = 0
m = n // 2
while m >= 2:
add_val = (ph.sum_phi(m) - 1) % MOD
ans = (ans + add_val) % MOD
m //= 2
return ans
def solve():
ph = PhiSummatory(5000000)
ans = solve_impl(100000000000, ph)
return str(ans)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.HashMap;
public class Euler643 {
static final long kMod = 1000000007L;
static final long kInv2 = 500000004L;
static class PhiSummatory {
int limit;
ArrayList<Integer> primes;
int[] lp;
int[] phi;
long[] pref;
HashMap<Long, Long> memo;
PhiSummatory(int n) {
limit = n;
lp = new int[n + 1];
phi = new int[n + 1];
pref = new long[n + 1];
primes = new ArrayList<>(n / 10);
memo = new HashMap<>();
phi[1] = 1;
for (int i = 2; i <= n; ++i) {
if (lp[i] == 0) {
lp[i] = i;
primes.add(i);
phi[i] = i - 1;
}
for (int p : primes) {
if (p > lp[i] || (long) i * p > n)
break;
lp[i * p] = p;
if (p == lp[i]) {
phi[i * p] = phi[i] * p;
break;
}
phi[i * p] = phi[i] * (p - 1);
}
}
for (int i = 1; i <= n; ++i) {
pref[i] = (pref[i - 1] + phi[i]) % kMod;
}
}
long sum_phi(long n) {
if (n <= limit)
return pref[(int) n];
Long val = memo.get(n);
if (val != null)
return val;
long nn = n % kMod;
long res = (nn * ((n + 1) % kMod)) % kMod;
res = (res * kInv2) % kMod;
for (long l = 2; l <= n;) {
long q = n / l;
long r = n / q;
long cnt = (r - l + 1) % kMod;
res = (res + kMod - (cnt * sum_phi(q)) % kMod) % kMod;
l = r + 1;
}
memo.put(n, res);
return res;
}
}
static long solveImpl(long n, PhiSummatory ph) {
long ans = 0;
for (long m = n / 2; m >= 2; m /= 2) {
long add = (ph.sum_phi(m) + kMod - 1) % kMod;
ans = (ans + add) % kMod;
}
return ans;
}
public static String solve() {
PhiSummatory ph = new PhiSummatory(5000000);
long ans = solveImpl(100000000000L, ph);
return Long.toString(ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}