【min25筛】【CF2020F】Count Leaves

题目

定义 $f(n,0)=1$，$f(n,d)={\sum }_{k|n}f(k,d-1)$
给出 $n,k,d$，你需要求出: $$\sum _{i=1}^{n}f(i^k,d)mod(10^9+7)$$
$n\le 1e9，k,d\le 1e5$
原题链接

思路

因数？数据范围这么大？那这玩意肯定是积性函数。于是我们通过观察，发现对于固定的 $k,d$，$f(i^k,d)$ 就是积性函数。但注意，这不是完全积性函数。
所以我们要想一想 $f(p^k,d)，pisprime$ 怎么求。
注意到，$p^k$ 的因数是：$p^0,p^1,...,p^k$，这相当于一个 $k\times d$ 的网格，你从左上角走到右下角的方案数。于是有：
$$f(p^k,d)=C(d+k,k)$$
我们惊喜的发现，如果 $k,d$ 不变，这就是个定值。定值也是多项式的一种，所以可以用 min25 筛。
但是这又和传统的 min25 筛，因为我们要求 $f(p^{ek},d)$ 的值。但其实这玩意只有在求 $S$ 的时候会发生（毕竟 g 只是求所有质数的前缀和），依旧是很好做的。

一定要注意的是，求 $g$ 的时候，我们是在对常数求 min25，注意不要多乘一个 $C(d+k,k)$

代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=5e6+7,inf=1e18,mod=1e9+7;
vector<int> p,sp,g,id1,id2,w;
int n,K,d;
int power(int x,int t)
{
	int b=1;
	while(t)
	{
		if(t&1) b=b*x%mod;
		x=x*x%mod; t>>=1;
	}
	return b;
}
vector<int> fac,unfac;
void initf(int n)
{
	fac.assign(n+1,0);
	unfac.assign(n+1,0);
	fac[0]=1;
	for(int i=1; i<=n; i++)
		fac[i]=fac[i-1]*i%mod;
	unfac[n]=power(fac[n],mod-2);
	for(int i=n-1; i>=0; i--)
		unfac[i]=unfac[i+1]*(i+1)%mod;
}
int C(int x,int y)
{
	if(x<y) return 0;
	return fac[x]*unfac[y]%mod*unfac[x-y]%mod;
}
int tot,sqr,cp;
void init(int n)
{
	p.clear();
	p.push_back(0);
	sp.assign(2*n+7,0);
	w.assign(2*n+7,0);
	g.assign(2*n+7,0);
	id1.assign(2*n+7,0);
	id2.assign(2*n+7,0);
	tot=0;
	
	vector<bool> vis(n+1);
	for(int i=2; i<=n; i++)
	{
		if(!vis[i])
		{
			p.push_back(i);
			int now=p.size()-1;
			sp[now]=(cp*now)%mod;
		}
		for(auto j:p)
		{
			if(!j) continue;
			if(i*j>n) break;
			vis[i*j]=1;
			if(i%j==0) break;
		}
	}
}
int S(int x,int y)
{
	if(x<=1||p[y]>=x) return 0;
	int k=(x<=sqr)?id1[x]:id2[n/x];
	int ans=(mod+g[k]-sp[y])%mod;
	for(int k=y+1; k<p.size()&&p[k]*p[k]<=x; k++)
	{
		int t=p[k];
		for(int e=1; t<=x; e++,t*=p[k])
		{
//			int p=t%mod;
			(ans+=C(d+e*K,d)*(S(x/t,k)+(e!=1))%mod)%=mod;
		}
	}
	return ans;
}
void O_o()
{
	cin>>n>>K>>d;
	cp=C(K+d,d);
	sqr=sqrt(n);
	init(sqr);
	for(int i=1,j; i<=n; i=j+1)
	{
		j=n/(n/i);
		w[++tot]=n/i;
		int now=w[tot]%mod;
		g[tot]=cp*(now-1)%mod;
		if(w[tot]<=sqr)
			id1[w[tot]]=tot;
		else
			id2[n/w[tot]]=tot;
	}
	for(int i=1; i<p.size(); i++)
	{
		for(int j=1; j<=tot&&p[i]*p[i]<=w[j]; j++)
		{
			int k=w[j]/p[i]<=sqr?id1[w[j]/p[i]]:id2[n/(w[j]/p[i])];
			(g[j]+=mod-(g[k]-sp[i-1]+mod)%mod)%=mod;
		}
	}
	int ans=S(n,0)+1;
	cout<<ans%mod<<"\n";
}
signed main()
{
	ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
	cout<<fixed<<setprecision(12);
	int T=1;
	initf(N);
	cin>>T;
	while(T--)
	{
		O_o();
	}
}