Problem 193: Squarefree Numbers

Mathematical Foundation

Theorem 1. (Squarefree Indicator Identity.) For all positive integers $n$,

where $\mu$ is the Mobius function. In particular, $|\mu(n)|^2 = 1$ if $n$ is squarefree and $0$ otherwise.

Proof. Both sides are multiplicative functions of $n$, so it suffices to verify equality at prime powers $n = p^k$.

• For $k = 1$: $|\mu(p)|^2 = 1$. The divisors $d$ with $d^2 \mid p$ are only $d = 1$, so $\sum_{d^2 \mid p} \mu(d) = \mu(1) = 1$. $\square_{\text{case}}$
• For $k \geq 2$: $|\mu(p^k)|^2 = 0$. The divisors $d$ with $d^2 \mid p^k$ are $d \in \{1, p, \ldots, p^{\lfloor k/2 \rfloor}\}$. For $k \geq 2$, at least $d = 1$ and $d = p$ appear, giving $\mu(1) + \mu(p) + \cdots$. Since $\mu(p^j) = 0$ for $j \geq 2$, this equals $\mu(1) + \mu(p) = 1 + (-1) = 0$. $\square_{\text{case}}$

By multiplicativity, the identity holds for all $n$. $\square$

Theorem 2. (Squarefree Counting Formula.) The number of squarefree integers not exceeding $N$ is

Proof. By Theorem 1:

The exchange of summation is justified since $d^2 \mid n$ and $n \leq N$ imply $d \leq \sqrt{N}$. $\square$

Lemma 1. (Mobius Sieve.) The values $\mu(k)$ for $1 \leq k \leq M$ can be computed in $O(M \log \log M)$ time and $O(M)$ space using a modified Eratosthenes sieve.

Proof. Initialize $\mu(k) = 1$ for all $k$. For each prime $p \leq M$: multiply $\mu(k)$ by $-1$ for all multiples $k$ of $p$; set $\mu(k) = 0$ for all multiples $k$ of $p^2$. Since each composite is processed once per prime factor, the total work is $\sum_{p \leq M} M/p = O(M \log \log M)$ by Mertens’ theorem. $\square$

Editorial

Count squarefree numbers below 2^50. Uses: Q(N) = sum_{k=1}^{sqrt(N)} mu(k) * floor(N/k^2). We sieve Mobius function. We then iterate over p from 2 to M. Finally, iterate over k from p to M step p.

Pseudocode

Sieve Mobius function
for p from 2 to M
for k from p to M step p
Compute Q(N)
for k from 1 to M

Complexity Analysis

• Time: $O(\sqrt{N} \log \log \sqrt{N})$ for the Mobius sieve, plus $O(\sqrt{N})$ for the summation. Total: $O(\sqrt{N} \log \log \sqrt{N})$.
• Space: $O(\sqrt{N})$ for storing the $\mu$ values. For $N = 2^{50}$, $\sqrt{N} = 2^{25} \approx 3.4 \times 10^7$.

Answer

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

int main() {
    // Count squarefree numbers below 2^50
    // Q(N) = sum_{k=1}^{sqrt(N)} mu(k) * floor(N/k^2)
    const long long N = (1LL << 50) - 1;
    const int SQ = 1 << 25; // sqrt(2^50) = 2^25

    // Sieve Mobius function
    vector<int> mu(SQ + 1, 1);
    vector<bool> is_prime(SQ + 1, true);
    vector<int> primes;

    // We'll use a sieve that tracks the smallest prime factor
    // and whether a number is divisible by p^2 for any prime p
    // Simple approach: for each prime p, multiply mu by -1 for multiples of p,
    // set mu=0 for multiples of p^2

    for (int i = 2; i <= SQ; i++) {
        if (is_prime[i]) {
            primes.push_back(i);
            // i is prime
            for (long long j = i; j <= SQ; j += i) {
                is_prime[j] = (j == i);
                mu[j] *= -1;
            }
            long long i2 = (long long)i * i;
            for (long long j = i2; j <= SQ; j += i2) {
                mu[j] = 0;
            }
        }
    }

    long long ans = 0;
    for (int k = 1; k <= SQ; k++) {
        if (mu[k] != 0) {
            ans += (long long)mu[k] * (N / ((long long)k * k));
        }
    }

    cout << ans << endl;
    return 0;
}

"""
Problem 193: Squarefree Numbers

Count squarefree numbers below 2^50.
Uses: Q(N) = sum_{k=1}^{sqrt(N)} mu(k) * floor(N/k^2)
"""

def solve():
    N = (1 << 50) - 1
    SQ = 1 << 25  # 2^25 = 33554432

    # Sieve Mobius function up to SQ
    mu = [0] * (SQ + 1)
    mu[1] = 1

    # Use a sieve: for each prime p, multiply mu by -1 for multiples,
    # set mu=0 for multiples of p^2
    is_prime = bytearray(b'\x01') * (SQ + 1)
    is_prime[0] = is_prime[1] = 0
    mu = [1] * (SQ + 1)

    for i in range(2, SQ + 1):
        if is_prime[i]:
            for j in range(i, SQ + 1, i):
                if j != i:
                    is_prime[j] = 0
                mu[j] *= -1
            i2 = i * i
            for j in range(i2, SQ + 1, i2):
                mu[j] = 0

    ans = 0
    for k in range(1, SQ + 1):
        if mu[k] != 0:
            ans += mu[k] * (N // (k * k))

    return ans

print(solve())