Squarefree Numbers
How many squarefree numbers are there below 2^50? A number is squarefree if it is not divisible by any perfect square other than 1.
Problem Statement
This archive keeps the full statement, math, and original media on the page.
A positive integer \(n\) is called squarefree, if no square of a prime divides \(n\), thus \(1, 2, 3, 5, 6, 7, 10, 11\) are squarefree, but not \(4, 8, 9, 12\).
How many squarefree numbers are there below \(2^{50}\)?
Problem 193: Squarefree Numbers
Mathematical Development
Theorem 1. (Squarefree Indicator Identity.) For all positive integers ,
where is the Mobius function. In particular, if is squarefree and otherwise.
Proof. Both sides are multiplicative functions of , so it suffices to verify equality at prime powers .
- For : . The divisors with are only , so .
- For : . The divisors with are . For , at least and appear, giving . Since for , this equals .
By multiplicativity, the identity holds for all .
Theorem 2. (Squarefree Counting Formula.) The number of squarefree integers not exceeding is
Proof. By Theorem 1:
The exchange of summation is justified since and imply .
Lemma 1. (Mobius Sieve.) The values for can be computed in time and space using a modified Eratosthenes sieve.
Proof. Initialize for all . For each prime : multiply by for all multiples of ; set for all multiples of . Since each composite is processed once per prime factor, the total work is by Mertens’ theorem.
Editorial
The direct definition of squarefree numbers is awkward for counting, but the Mobius inversion formula turns the indicator into a divisor sum over squares. After that transformation, the problem stops being about testing individual integers and becomes a single summation
So the whole task reduces to computing the Mobius function up to sqrt(2^50) = 2^25.
That is exactly what the modified sieve does. Each prime flips the sign of mu on its multiples, and each square of a prime wipes the value back to zero. Once the sieve is finished, one linear pass through k = 1 .. 2^25 evaluates the formula above and accumulates the squarefree count.
Pseudocode
Let N = 2^50 - 1 and M = floor(sqrt(N)).
Initialize mu(1..M) for a Mobius sieve.
For each prime p up to M:
flip the sign of mu on every multiple of p;
set mu to 0 on every multiple of p^2.
Set the answer to 0.
For each k from 1 to M:
add mu(k) * floor(N / k^2) to the answer.
Return the answer.
Complexity Analysis
- Time: for the Mobius sieve, plus for the summation. Total: .
- Space: for storing the values. For , .
Answer
Code
Each problem page includes the exact C++ and Python source files from the local archive.
#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())