powerful number筛

心血来潮跑来实现以下这个东西

我们应该知道杜教筛的理论是 \(f * g=h\)，那么问题在于如何找 \(g\)。

之前的blog应该提到过可以令 \(g(p)=-f(p)\)，这样一来 \(h\) 就只会在PN处有值。于是可以大力爆搜 \(h\)，而 \(g\) 的块筛又很好处理。

但是这样复杂度会有一个下限为 \(O(n^{\frac 2 3})\)，有没有办法去除呢？

办法是有的，反过来，设 \(h * g=f\)。

此时我们构造 \(g(p)=f(p)\) 即可得到和上面相同的结论，但此时只需处理 \(g\) 的块筛即可，复杂度下降至 \(O(\frac{n^{\frac 3 4}}{\log n})\) 或者更低。

问题来了，当 \(k>2\) 时，\(g(p^k)\) 应该是多少？

实际上是多少否无所谓，因为有 \(\sum_{i=0}^k h(p^i)g(p^{k-i})=f(p^k)\)，一般情况令 \(g(p^k)=0\)。

但是我在实现的时候 推 错 了 \(g\) 的 块 筛 柿 子，懒得重新推。又发现我的柿子实际上是令 \(g(p^k)=f(p)^k\)，如果直接暴力做卷积那么复杂度会变成 \(O(\sqrt n\log n)\)，于是来优化一下：

\[\sum_{i=0}^kh(p^i)f(p)^{k-i}=f(p^k) \]

\[f(p)^k\sum_{i=0}^kh(p^i)f(p)^{-i}=f(p^k) \]

于是处理出 \(\frac{f(p^k)}{f(p)^k}\)，然后来个差分求逆元就好，复杂度变回了 \(O(\sqrt n)\)，不过需要存在逆元才行。

但是？

\[h(p^k)f(p)^{-k}=\frac{f(p^k)}{f(p)^k}-\frac{f(p^{k-1})}{f(p)^{k-1}} \]

\[h(p^k)=f(p^k)-f(p^{k-1})f(p) \]

不需要逆元也可以。

这里丢一下 DIVCNTK 的实现，目前是 lgrk2，spojrk7：

#include
#include
typedef unsigned ui;
typedef __uint128_t L;
typedef unsigned long long ull;
const ui M=1e5+5;
ull n,K,h[M],F[M],G[M],f0[M],g0[M];L B[M];ui S,top,pri[M];
ull DFS(const ull&n,const ui&k){
	ull ans=n<=S?G[n]+1:F[::n/n]+1;
	for(ui i=k+1;i<=top&&1ull*pri[i]*pri[i]<=n;++i){
		ull*H=h+2;
		for(ull N=(n*B[pri[i]]>>64)*B[pri[i]]>>64;N;N=N*B[pri[i]]>>64){
			ans+=*H++*DFS(N,i);
		}
	}
	return ans;
}
inline ull PN(const ull&n){
	ui i,j,k,tp=0,sqr;ull w,lim;top=0;
	for(S=1;1ull*S*S<=n;++S)f0[S]=(n*B[S]>>64)-1,g0[S]=S-1;
	sqr=sqrt(--S);
	for(i=2;i<=S;++i)if(g0[i]^g0[i-1]){
		w=n*B[i]>>64;lim=w*B[i]>>64;if(lim>S)lim=S;k=S*B[i]>>64;++tp;
		const ull&S0=g0[i-1];
		for(j=1;j<=k;++j)f0[j]-=f0[i*j]-S0;
		for(;j<=lim;++j)f0[j]-=g0[w*B[j]>>64]-S0;
		if(i<=sqr){
			for(lim=k*i,j=S;k>=i;lim-=i,--k){
				for(const ull&V0=g0[k]-S0;j>=lim;--j)g0[j]-=V0;
			}
		}
	}
	for(i=1;i<=S;++i)F[i]=f0[i]*K,G[i]=g0[i]*K;top=tp++;
	for(i=S;i>1;--i)if(g0[i]^g0[i-1]){
		const ull&g=g0[i-1]*K;
		if(i<=sqr){
			for(k=i,lim=(k+1)*i,j=lim-i;lim<=S;lim+=i,++k){
				for(const ull&V=K*(G[k]-g);j>64;lim=w*B[i]>>64;if(lim>S)lim=S;k=S*B[i]>>64;pri[--tp]=i;
		for(j=lim;j>k;--j)F[j]+=K*(G[w*B[j]>>64]-g);
		for(j=k;j>=1;--j)F[j]+=K*(F[i*j]-g);
	}
	lim=log(n)/log(2);K=(K-1)*(K-1);
	for(i=1;i<=lim;++i)h[i]=-K*(i-1);
	return DFS(n,0);
}
signed main(){
	ui T;for(T=1;T