题目
定义 f ( n , 0 ) = 1 f(n,0)=1 f(n,0)=1, f ( n , d ) = ∑ k ∣ n f ( k , d − 1 ) f(n,d)=\sum_{k|n}f(k,d-1) f(n,d)=∑k∣nf(k,d−1)
给出 n , k , d n,k,d n,k,d,你需要求出: ∑ i = 1 n f ( i k , d ) m o d ( 1 0 9 + 7 ) \sum_{i=1}^n f(i^k,d) \ mod\ (10^9+7) i=1∑nf(ik,d) mod (109+7)
n ≤ 1 e 9 , k , d ≤ 1 e 5 n\leq 1e9,k,d\leq1e5 n≤1e9,k,d≤1e5
原题链接
思路
因数?数据范围这么大?那这玩意肯定是积性函数。于是我们通过观察,发现对于固定的 k , d k,d k,d, f ( i k , d ) f(i^k,d) f(ik,d) 就是积性函数。但注意,这不是完全积性函数。
所以我们要想一想 f ( p k , d ) , p i s p r i m e f(p^k,d),p\ is\ prime f(pk,d),p is prime 怎么求。
注意到, p k p^k pk 的因数是: p 0 , p 1 , . . . , p k p^0, p^1,...,p^k p0,p1,...,pk,这相当于一个 k × d k\times d k×d 的网格,你从左上角走到右下角的方案数。于是有:
f ( p k , d ) = C ( d + k , k ) f(p^k,d)=C(d+k,k) f(pk,d)=C(d+k,k)
我们惊喜的发现,如果 k , d k,d k,d 不变,这就是个定值。定值也是多项式的一种,所以可以用 min25 筛。
但是这又和传统的 min25 筛,因为我们要求 f ( p e k , d ) f(p^{ek},d) f(pek,d) 的值。但其实这玩意只有在求 S S S 的时候会发生(毕竟 g 只是求所有质数的前缀和),依旧是很好做的。
一定要注意的是,求 g g g 的时候,我们是在对常数求 min25,注意不要多乘一个 C ( d + k , k ) C(d+k,k) 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();}
}