[DarkBZOJ2440]完全平方数

题目

传送门 to DarkBZOJ

题目概要
求这个集合中第 $k$ 小的值：
$$\{x\in {\mathbb{N}}^{+}|\mu (x)\ne 0\}$$

数据范围与约定
$k⩽10^9$ 。数据组数 $T⩽50$ 。事实上这个数据范围非常宽松。

思路

其实蛮简单一道题？

我一开始想，每个质因数的指数不超过 $1$，那么最多只需要 $30$ 个质数。于是我把最小的 $30$ 个质数乘起来，发现爆 $\mathtt{l}\mathtt{o}\mathtt{n}\mathtt{g}\mathtt{l}\mathtt{o}\mathtt{n}\mathtt{g}$ 了，以为不可做……

仔细想想，只用前 $30$ 个质数固然可以凑出 $k$ 个数，但是使用其他的质数，也可以得到一些数字，并且比只用前 $30$ 个质数的最大数更小。一个非常粗略的估计是，$x$ 以内的质数个数约为 $\frac{x}{\ln x}$，假如这个数大于 $k$，就一定可行了。$x=10^{11}$ 时，已经满足，显然真正的第 $10^9$ 个数会小很多。

既然上界是有限的，就可以进行二分。既然不能含有平方因子，就直接做容斥，
$$\sum _{p_1\in P}\lfloor \frac{n}{p_1^2}\rfloor -\sum _{p_1,p_2\in P}\lfloor \frac{n}{p_1^2{p}_2^2}\rfloor +\cdots$$

就可以用到 $\mu$ 了，可以写成
$$\sum _{x=2}^{+\infty }\mu (x)\lfloor \frac{n}{x^2}\rfloor$$

直接枚举 $x$，复杂度是 $\mathcal{O}(\sqrt{n})$ 的。是否能做到更快呢？显然这里有很明显的整除分块倾向。

然而我一估算：$x\in [n^{1/4},n^{1/2}]$ 时，$\lfloor \frac{n}{x^2}\rfloor \in [1,\sqrt{n}]$，两边都是 $\mathcal{O}(\sqrt{n})$ 的范围，那么约莫是一一对应，故复杂度始终为 $\sqrt{n}$ 。可是仿照整除分块的复杂度证明：

• 当 $x⩽n^{1/3}$ 时，最多 $\mathcal{O}(\sqrt[3]{n})$ 个 $x$ 。
• 当 $x>\sqrt[3]{n}$，有 $\lfloor \frac{n}{x^2}\rfloor ⩽\sqrt[3]{n}$ 。

所以整除分块的复杂度是 $\mathcal{O}(\sqrt[3]{n})$ 的，是个大优化！求 $x$ 的范围，即求 $\sqrt{\frac{n}{v}}(v⩽\sqrt[3]{n})$，可以预处理实数值 $\sqrt{v}$，每组询问求一次 $\sqrt{n}$ 的实数值，用 $x=\lfloor \frac{\sqrt{n}}{\sqrt{v}}\rfloor$ 作为估计值，然后牛顿迭代一次、微调一次。

总时间复杂度 $\mathcal{O}(\sqrt{n}+T\sqrt[3]{n}\log n)$ 。

代码

#include <cstdio>
#include <iostream>
#include <vector>
#include <cstring>
#include <algorithm>
#include <cctype>
#include <cmath>
using namespace std;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
inline int readint(){
	int a = 0, c = getchar(), f = 1;
	for(; !isdigit(c); c=getchar())
		if(c == '-') f = -f;
	for(; isdigit(c); c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}
typedef long long int_;
inline void writeUnsigned(const unsigned &x){
	if(x > 9) writeUnsigned(x/10);
	putchar(char((x-x/10*10)^48));
}

const int MAXN = 40600;
int primes[MAXN], primes_size;
bool isPrime[MAXN]; int mu[MAXN];
int ddg[MAXN]; ///< prefix sum of @a mu
void sieve(const int &n){
	memset(isPrime+2,true,n-1);
	for(int i=2; i<=n; ++i){
		if(isPrime[i]) mu[primes[++ primes_size] = i] = -1;
		for(int j=1; j<=primes_size&&primes[j]<=n/i; ++j){
			isPrime[i*primes[j]] = false;
			if(i%primes[j]) mu[i*primes[j]] = -mu[i];
			else { mu[i*primes[j]] = 0; break; }
		}
	}
	rep(i,mu[1]=1,n) ddg[i] = ddg[i-1]+mu[i];
}

const int BOUND = 4641;
double pfg[BOUND+1]; // square root
int_ solve(const int_ &n){
	int_ res = 0; int bound;
	for(int i=1; true; ++i){
		res += ((n/i)/i)*mu[i];
		if(n/i/i < i){ bound = i; break; }
	}
	if(bound > (n-1)/bound) return res;
	double sqrtn = sqrt(n);
	for(int v=1,l,r=int(sqrtn); true; ++v,r=l){
		l = int(sqrtn/pfg[v+1]);
		/// do Newton Iteration once
		l = int((l+(n/(v+1)/l))>>1);
		if(l > n/(v+1)/l) -- l; // floor
		if(l <= bound){ // exit
			res += int_(ddg[r]-ddg[bound])*v;
			break; // the last shot
		}
		res += int_(ddg[r]-ddg[l])*v;
	}
	return res;
}

int main(){
	sieve(MAXN-5); // get mu
	rep(i,1,BOUND) pfg[i] = sqrt(i);
	for(int T=readint(),k; T; --T){
		k = readint();
		unsigned l = 0, r = 1644934081, x;
		for(x=r>>1; l!=r; x=((l+r)>>1))
			(solve(x) >= k) ? (r = x) : (l = x+1);
		writeUnsigned(l), putchar('\n');
	}
	return 0;
}