Problem 905: Lattice Point Visibility

Mathematical Analysis

Mobius Inversion

The indicator $[\gcd(a,b) = 1]$ can be expressed via the Mobius function:

This is the fundamental identity of Mobius inversion, following from $\sum_{d \mid n} \mu(d) = [n = 1]$.

Closed-Form Sum

Proof. Substituting (1) and swapping summation order:

Asymptotic Density

Theorem. $\lim_{N \to \infty} \frac{V(N)}{N^2} = \frac{6}{\pi^2} = \frac{1}{\zeta(2)}$.

Proof. $\frac{V(N)}{N^2} = \sum_{d=1}^{N} \mu(d) \left(\frac{\lfloor N/d \rfloor}{N}\right)^2 \to \sum_{d=1}^{\infty} \frac{\mu(d)}{d^2} = \frac{1}{\zeta(2)} = \frac{6}{\pi^2}$. The interchange of limit and sum is justified by absolute convergence. The identity $\sum_{d=1}^{\infty} \mu(d)/d^s = 1/\zeta(s)$ for $\text{Re}(s) > 1$ is a standard result from the Euler product $\zeta(s) = \prod_p (1 - p^{-s})^{-1}$. $\square$

Relation to Euler’s Totient

Note that $V(N)$ is closely related to the totient summatory function:

More precisely, $V(N) = -1 + 2\sum_{k=1}^{N} \varphi(k)$ counts visible points $(a,b)$ with $1 \le a, b \le N$ (the $-1$ adjusts for the point $(1,1)$ counted once vs. twice in the Farey-based counting).

Actually, the exact relation is: counting coprime pairs $(a,b)$ with $1 \le a, b \le N$ equals $\sum_{d=1}^{N} \mu(d) \lfloor N/d \rfloor^2$. This also equals $2\left(\sum_{k=1}^{N} \varphi(k)\right) - 1$ since each Farey fraction $a/b$ with $\gcd(a,b) = 1$ and $1 \le a < b \le N$ corresponds to a visible point, and we count both $(a,b)$ and $(b,a)$ plus $(1,1)$.

Hyperbola Method Optimization

For large $N$, the sum (2) can be accelerated using the fact that $\lfloor N/d \rfloor$ takes only $O(\sqrt{N})$ distinct values. Grouping terms by equal values of $\lfloor N/d \rfloor$ and precomputing prefix sums of $\mu(d)$ (the Mertens function) yields an $O(\sqrt{N})$-time evaluation after an $O(N)$ sieve.

Concrete Examples

Editorial

Count visible lattice points (a,b) with 1 <= a,b <= N from the origin, where visibility means gcd(a,b) = 1. Find V(10^6) mod 10^9+7. Key formula: V(N) = sum_{d=1}^{N} mu(d) * floor(N/d)^2. We sieve** $\mu(d)$ for $d = 1, \ldots, N$ using a linear sieve in $O(N)$ time. Finally, evaluate** $V(N) = \sum_{d=1}^{N} \mu(d) \lfloor N/d \rfloor^2 \bmod (10^9+7)$.

Pseudocode

Sieve** $\mu(d)$ for $d = 1, \ldots, N$ using a linear sieve in $O(N)$ time
Evaluate** $V(N) = \sum_{d=1}^{N} \mu(d) \lfloor N/d \rfloor^2 \bmod (10^9+7)$

Proof of Correctness

Theorem. Formula (2) correctly counts $V(N)$.

Proof. The Mobius identity (1) is equivalent to the Dirichlet series identity $1/\zeta(s) = \sum \mu(n)/n^s$. Substituting into the double sum and rearranging is a standard application of Mobius inversion on the divisor lattice. Each coprime pair is counted exactly once. $\square$

Complexity Analysis

• Sieve: $O(N \log \log N)$ for Eratosthenes-based, or $O(N)$ for linear sieve.
• Summation: $O(N)$ for direct evaluation, $O(\sqrt{N})$ with hyperbola method.
• Space: $O(N)$ for the sieve array.

For $N = 10^6$, the sieve and sum complete in well under a second.

Answer

#include <bits/stdc++.h>
using namespace std;

/*
 * Problem 905: Lattice Point Visibility
 *
 * Count lattice points (a,b) with 1 <= a,b <= N and gcd(a,b) = 1.
 * Formula: V(N) = sum_{d=1}^{N} mu(d) * floor(N/d)^2
 *
 * Two methods:
 *   1. Direct Mobius summation: O(N) after O(N) sieve
 *   2. Hyperbola method: O(sqrt(N)) summation after sieve
 *
 * Asymptotic: V(N) ~ 6N^2/pi^2
 */

const long long MOD = 1e9 + 7;

/*
 * Linear sieve for Mobius function
 * mu[n] = 0 if n has a squared prime factor
 * mu[n] = (-1)^k if n is a product of k distinct primes
 */
vector<int> mobius_sieve(int n) {
    vector<int> mu(n + 1, 0);
    mu[1] = 1;
    vector<bool> is_prime(n + 1, true);
    vector<int> primes;
    for (int i = 2; i <= n; i++) {
        if (is_prime[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            if ((long long)i * p > n) break;
            is_prime[i * p] = false;
            if (i % p == 0) {
                mu[i * p] = 0;
                break;
            }
            mu[i * p] = -mu[i];
        }
    }
    return mu;
}

/*
 * Method 1: Direct summation
 * V(N) = sum_{d=1}^{N} mu(d) * floor(N/d)^2 mod p
 */
long long solve_direct(int N, const vector<int>& mu) {
    long long result = 0;
    for (int d = 1; d <= N; d++) {
        if (mu[d] != 0) {
            long long q = (N / d) % MOD;
            long long term = q * q % MOD;
            if (mu[d] == 1) {
                result = (result + term) % MOD;
            } else {
                result = (result - term + MOD) % MOD;
            }
        }
    }
    return result;
}

/*
 * Method 2: Hyperbola method with prefix sums of mu
 * Group terms by the value of floor(N/d), which takes O(sqrt(N)) values.
 * Requires prefix sums of mu (Mertens function).
 */
long long solve_hyperbola(int N, const vector<int>& mu) {
    // Mertens function: M(k) = sum_{d=1}^{k} mu(d)
    vector<long long> M(N + 1, 0);
    for (int i = 1; i <= N; i++)
        M[i] = M[i - 1] + mu[i];

    long long result = 0;
    int d = 1;
    while (d <= N) {
        int q = N / d;
        int d2 = N / q;  // largest d' with floor(N/d') = q
        // Sum mu[d..d2] = M(d2) - M(d-1)
        long long mu_sum = ((M[d2] - M[d - 1]) % MOD + MOD) % MOD;
        long long qm = q % MOD;
        result = (result + mu_sum % MOD * (qm * qm % MOD)) % MOD;
        d = d2 + 1;
    }
    return result;
}

int main() {
    const int N = 1000000;

    auto mu = mobius_sieve(N);

    long long ans1 = solve_direct(N, mu);
    long long ans2 = solve_hyperbola(N, mu);

    assert(ans1 == ans2);

    // Verify small case: V(10) = 63
    auto mu10 = mobius_sieve(10);
    long long v10 = 0;
    for (int d = 1; d <= 10; d++) {
        if (mu10[d] != 0) {
            long long q = 10 / d;
            v10 += mu10[d] * q * q;
        }
    }
    assert(v10 == 63);

    cout << ans1 << endl;
    return 0;
}

"""
Problem 905: Lattice Point Visibility

Count visible lattice points (a,b) with 1 <= a,b <= N from the origin,
where visibility means gcd(a,b) = 1.  Find V(10^6) mod 10^9+7.

Key formula: V(N) = sum_{d=1}^{N} mu(d) * floor(N/d)^2

Methods:
    1. Mobius sieve + direct summation - O(N log log N)
    2. Euler totient summation (cross-check for small N)
    3. Brute force gcd check (verification for small N)
"""

from math import gcd

MOD = 10**9 + 7

# Mobius sieve (linear sieve)
def mobius_sieve(n: int) -> list:
    """Compute mu(k) for k = 0, 1, ..., n via linear sieve."""
    mu = [0] * (n + 1)
    mu[1] = 1
    is_prime = [True] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if is_prime[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            if i * p > n:
                break
            is_prime[i * p] = False
            if i % p == 0:
                mu[i * p] = 0
                break
            mu[i * p] = -mu[i]
    return mu

def visible_count_mobius(N: int, mod: int):
    """V(N) = sum_{d=1}^{N} mu(d) * floor(N/d)^2 mod p."""
    mu = mobius_sieve(N)
    total = 0
    for d in range(1, N + 1):
        if mu[d] != 0:
            q = N // d
            total = (total + mu[d] * q * q) % mod
    return total % mod

def visible_count_totient(N: int):
    """V(N) = 2 * sum_{k=1}^{N} phi(k) - 1."""
    phi = list(range(N + 1))
    for p in range(2, N + 1):
        if phi[p] == p:  # p is prime
            for multiple in range(p, N + 1, p):
                phi[multiple] -= phi[multiple] // p
    return 2 * sum(phi[1:]) - 1

def visible_count_brute(N: int):
    """Count pairs (a,b) with gcd(a,b) = 1 by direct computation."""
    count = 0
    for a in range(1, N + 1):
        for b in range(1, N + 1):
            if gcd(a, b) == 1:
                count += 1
    return count

# Solve
N = 10**6

ans = visible_count_mobius(N, MOD)

# Verify on small cases with totient
for test_n in [5, 10, 50, 100]:
    v_mob = visible_count_mobius(test_n, 10**18)
    v_tot = visible_count_totient(test_n)
    assert v_mob == v_tot, f"N={test_n}: Mobius={v_mob}, Totient={v_tot}"

# Small brute force
for test_n in [5, 10]:
    v_brute = visible_count_brute(test_n)
    v_mob = visible_count_mobius(test_n, 10**18)
    assert v_brute == v_mob, f"N={test_n}: brute={v_brute}, Mobius={v_mob}"

print(ans)