powerful number筛

update on 25.08.15 之前写的太难懂了，重构一下

这个 PN 筛的做法和一般的 PN 筛不太一样，做法类似于将叶筛 DP 正反各做一次。

问题：对于积性函数 \(f\) 满足 \(f(p^k)\) 在 \(k\geq 2\) 时能快速求单点且在 \(k=1\) 时为关于 \(p\) 的多项式，求 \(f\) 的块筛。

我们设积性函数 \(g,h\) 满足 \(g*h=f\)，且 \(g\) 为完全积性函数，\(h\) 为满足对于任意 \(p\) 有 \(h(p)=0\) 的积性函数。

我们令 \(g(p^k)=f(p)^k\)，推导 \(h\) 的通项：

\[\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) \]

\[\sum_{i=0}^k\frac{h(p^i)}{f(p)^i}=\frac{f(p^k)}{f(p)^k} \]

\[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) \]

由于 \(f(p^k)\) 可能为 \(0\)，所以验证一下：

\[\sum_{i=1}^k(f(p^i)-f(p^{i-1})f(p))f(p)^{k-i}+f(p)^k \]

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

\[f(p^k)-f(p)^k+f(p)^k=f(p^k) \]

于是在获得了 \(g\) 的块筛后，我们就可以 \(O(\sqrt n\log n)\) 或 \(O(n^{\frac{2}{3}})\) 获得 \(f\) 的块筛。

考虑如何获得 \(g\) 的块筛。我们注意到叶筛 DP 实际上是将完全积性函数前缀和 \(\sum_{i=1}^n f(i)\)筛成了完全积性函数在质数位置的前缀和 \(\sum_{p\leq n}f(p)\)，因此我们可以反着做这个过程。

具体地，先求出 \(\sum_{p\leq n}g(p)=\sum_{p\leq n}f(p)\)，然后做叶筛的转置即可得到 \(\sum_{i=1}^n g(i)\) 的块筛。

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

#include<cstdio>
#include<cmath>
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<lim;++j)G[j]+=V;
			}
			for(const ull&V=K*(G[k]-g);j<=S;++j)G[j]+=V;
		}
		w=n*B[i]>>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<M;++T)B[T]=((L(1)<<64)+T-1)/T;scanf("%u",&T);
	while(T--)scanf("%llu%llu",&n,&K),++K,printf("%llu\n",PN(n));
}

upd:写了个人能看的 P5235 代码：

#include<cstdio>
#include<cmath>
typedef long long ll;
const int M=1e5+5,mod=1e9+7;
ll n;int S,F[M],G[M],f1[M],f2[M],g1[M],g2[M];int top,*h[M],pri[M],buf[M*20],*now;
inline ll min(const ll&a,const ll&b){
	return a>b?b:a;
}
inline int DFS(ll N,int k){
	int ans=N<=S?G[N]+1:F[n/N]+1;
	for(int i=k+1;i<=top&&1ll*pri[i]*pri[i]<=N;++i){
		int k(2);
		for(ll PK=1ll*pri[i]*pri[i];PK<=N;PK*=pri[i]){
			ans=(ans+1ll*h[i][k++]*DFS(N/PK,i))%mod;
		}
	}
	return ans;
}
inline int PN(const ll&n){
	int i,j,k,tp=0,sqr,lim;ll w;now=buf;top=0;
	for(S=1;1ll*S*S<=n;++S){
		w=n/S%mod;
		f1[S]=w*(w+1ll)/2%mod;f2[S]=f1[S]*(2*w+1ll)%mod*333333336%mod;--f1[S];--f2[S];
		g1[S]=S*(S+1ll)/2%mod;g2[S]=g1[S]*(2*S+1ll)%mod*333333336%mod;--g1[S];--g2[S];
	}
	sqr=sqrt(--S);
	for(i=2;i<=S;++i)if(g1[i]^g1[i-1]){
		lim=min(S,n/i/i);++tp;
		int P1=i,P2=1ll*i*i%mod,S1=g1[i-1],S2=g2[i-1];
		for(j=1;j<=lim;++j){
			if(i*j<=S){
				f1[j]=(f1[j]+mod-1ll*P1*(f1[i*j]+mod-S1)%mod)%mod;
				f2[j]=(f2[j]+mod-1ll*P2*(f2[i*j]+mod-S2)%mod)%mod;
			}
			else{
				f1[j]=(f1[j]+mod-1ll*P1*(g1[n/i/j]+mod-S1)%mod)%mod;
				f2[j]=(f2[j]+mod-1ll*P2*(g2[n/i/j]+mod-S2)%mod)%mod;
			}
		}
		if(i<=sqr){
			for(j=S;j>=i*i;--j){
				g1[j]=(g1[j]+mod-1ll*P1*(g1[j/i]+mod-S1)%mod)%mod;
				g2[j]=(g2[j]+mod-1ll*P2*(g2[j/i]+mod-S2)%mod)%mod;
			}
		}
	}
	for(i=1;i<=S;++i)F[i]=(mod+f2[i]-f1[i])%mod,G[i]=(mod+g2[i]-g1[i])%mod;top=tp++;
	for(i=S;i>1;--i)if(g1[i]^g1[i-1]){
		const int&f=i*(i-1ll)%mod,&g=(mod+g2[i-1]-g1[i-1])%mod;
		if(i<=sqr){
			for(j=i*i;j<=S;++j){
				G[j]=(G[j]+1ll*f*(G[j/i]+mod-g))%mod;
			}
		}
		lim=min(S,n/i/i);
		for(int j=lim;j>=1;--j){
			if(i*j<=S)F[j]=(F[j]+1ll*f*(F[i*j]+mod-g))%mod;
			else F[j]=(F[j]+1ll*f*(G[n/i/j]+mod-g))%mod;
		}
		pri[--tp]=i;h[tp]=now;lim=log(n)/log(i);
		for(w=i,j=1;j<=lim;++j,w=w*i%mod)now[j]=w*(w-1)%mod;
		for(j=lim;j>1;--j)now[j]=(now[j]+1ll*(mod-f)*now[j-1])%mod;now+=lim+1;
	}
	return DFS(n,0);
}
signed main(){
	scanf("%lld",&n);printf("%d",PN(n));
}