【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;elseid2[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();}
}