Counting Square-free Numbers
National Olympiad — National round (2023 season)
Submission name: square_free.cpp
Problem summary
For each n, count positive integers at most n that are not divisible by the square of any integer greater than one.
Input
q, followed by q values n.
Output
The square-free count for each n.
Constraints
1 ≤ q ≤ 500; 1 ≤ n ≤ 10^14.
Worked example
An independently prepared example for this explanation.
Input
2 4 10
Output
3 7
Solution analysis
- Inclusion–exclusion over prime squares gives count(n)=Σ mu(d)*floor(n/d²), for d≤sqrt(n), where mu is the Möbius function. Precompute mu and its prefix sums up to sqrt(maximum query) with a linear sieve.
- For small d, evaluate the terms directly. Beyond a cutoff of 90000, group consecutive d values having the same quotient v=floor(n/d²). The final index of a group is floor(sqrt(n/v)), and its contribution is v*(prefixMu[end]−prefixMu[start−1]). Correct floating-point square-root estimates with integer comparisons.
- For n=10, the nonzero terms are 10−floor(10/4)−floor(10/9)=10−2−1=7. The numbers are 1,2,3,5,6,7,10. Grouping affects only evaluation speed and leaves this exact inclusion–exclusion sum unchanged.
Why it works
For an integer x, summing mu(d) over all d whose square divides x equals one if x has no repeated prime factor, and zero otherwise: it is the product of (1−1) over its squared prime divisors. Summing this indicator over x≤n yields the formula. Every quotient group has a constant floor(n/d²), so replacing its individual Möbius terms by a prefix-sum difference is exact.
Watch out. Use 64-bit n and d*d. Include d=1 so that 1 is counted. Correct square-root rounding at perfect squares and reuse preprocessing for all queries.
Commented C++17 solution
An independently written complete program with standard headers.
#include<iostream>
#include<vector>
#include<cmath>
#include<algorithm>
// Algorithm: Inclusion–exclusion over prime squares gives count(n)=Σ mu(d)*floor(n/d²), for
// d≤sqrt(n), where mu is the Möbius function.
using namespace std;
long long integerRoot(long long value){
long long root=static_cast<long long>(sqrt(static_cast<long double>(value)));
while((root+1)*(root+1)<=value){
++root;
}
while(root*root>value){
--root;
}
return root;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(nullptr);
int q;
cin>>q;
vector<long long>queries(q);
long long maximum=1;
for(auto&value:queries){
cin>>value;
maximum=max(maximum,value);
}
int limit=static_cast<int>(integerRoot(maximum));
vector<signed char>mobius(limit+1,0);
vector<bool>composite(limit+1,false);
vector<int>primes,prefix(limit+1,0);
mobius[1]=1;
for(int value=2;value<=limit;++value){
if(!composite[value]){
primes.push_back(value);
mobius[value]=-1;
}
for(int prime:primes){
if(1LL*value*prime>limit){
break;
}
int product=value*prime;
composite[product]=true;
if(value%prime==0){
mobius[product]=0;
break;
}
mobius[product]=-mobius[value];
}
}
for(int i=1;i<=limit;++i){
prefix[i]=prefix[i-1]+mobius[i];
}
for(long long n:queries){
int end=static_cast<int>(integerRoot(n));
int cutoff=min(end,90000);
long long answer=0;
for(long long d=1;d<=cutoff;++d){
answer+=mobius[d]*(n/(d*d));
}
long long start=cutoff+1;
while(start<=end){
long long quotient=n/(start*start);
long long last=integerRoot(n/quotient);
answer+=quotient*(prefix[last]-prefix[start-1]);
start=last+1;
}
cout<<answer<<'\n';
}
}