P4464 [国家集训队]JZPKIL
看到题面底部有个 by gyz ，查了一下，这个人是当年NOIrank1，AK了提高组/fad，这种神仙出的必然是神仙题吧
发现这个 \(t|n\) 好像怎么都换不掉了，可是 \(n\le 10^{18}\) 啊。
忽然，脑子里闪过一个恐怖的东西（真不知道当时场上的选手是什么感觉）：这必然是 Pollard-Rho 吧！
那么假设我们已经 \(O(n^{\frac{1}{4}})\) 分解了质因数（这个就是PR干的事情）。
喂喂喂你2012年出PR是什么毒瘤啊，这东西现在都不普及吧
果然，集训队互测就是集训队选手往死里出题/fad
有一个非常重要的结论（被我特意收集到这里的）：\(10^{18}\) 以内，因数最多的数有 \(103680\) 个因数。
接下去这段是我的翻车现场，以为推完了，可以跳过qaq，只能拿部分分。
第一个 \(\sum\) ，直接枚举。
第二个 \(\sum\) ，考虑用伯努利数计算自然数幂的前缀和。这个可以直接用伯努利数推通项，然后 \(O(y)\) 计算。貌似拉格朗日插值也可以。
第三个 \(\sum\) ，其实我们就是要对于每一个 \(d|n\) 算出 \(\sum_{i|d}d^x\) ，可以 \(O(\sigma(n))\) 在分解完质因数之后枚举因数。
但是怎么统计呢？我想到了狄利克雷前缀和这个东西，之前看过模板题的描述，感觉一模一样啊。。。 不过不用貌似也没关系，复杂度一样。
复杂度是 \(O(T(n^{\frac{1}{4}}+\sigma(n)\log \sigma(n)+y\sigma(n)+\sigma(n)\log\log\sigma(n)))\)
诶，等等，里面怎么还有个 \(O(Ty\sigma(n))\) ，草草草，跑不过去的！！！
看了半天没思路，然后，就暴力把伯努利数求自然数幂和通项搞出来带进去吧。。。
自然数 \(k\) 次幂和的多项式用伯努利数写就是
伯努利数直接 \(O(y^2)\) 递推即可，相关内容可以看看 多项式笔记（二）
就是一个 \(k+1\) 次多项式，换一种写法，好看一点。
然后带进去
嗯？三个积性函数狄利克雷卷积起来是什么情况？
暴毙了啊，数论没学好啊，推了半天一点进展没有。
……
woc我就是个sb，积性函数不是在互质的时候可以直接乘起来的嘛。
然后把 \(n\) 分解质因数，得到 \(n=p_1^{k_1}p_2^{k_2}\cdots\) ，对于每一个质因子 \(p_i^{k_i}\) 算出答案然后乘起来不就完事了？
可能怎么算积性函数还要再分析一下。
分解完质因数只需要计算 \(ans(p^k)\) 即可。
注意到 \(g(n)\) 里面那个莫比乌斯函数在 \(n=p^k\) 的时候有很好的性质，当且仅当 \(k=0\) 或 \(k=1\) 时有值。
那么就枚举 \(g\) 的值，然后计算 \(f*h\) ，即 \(ans(n)=g(1)(f*h)(n)+g(p)(f*h)(\dfrac{n}{p})\)
设质因子 \(p\) 的最大次幂为 \(k\)
前半部分：
后半部分：
质因子次数很小，个数也很少，暴力卷就可以了。
复杂度：\(O(y^2+T(n^{\frac{1}{4}}+y\log n))\) 。
这个题好难码啊/kk，主要是一看这个时间限制以及使用的算法就感觉很卡常，还要码PR。边写边卡的后果就是代码自己都读不懂了qaq
调了好几个小时，原来是老坑了：伯努利求和的上界是 \(n-1\) ，所以要单独给 \(a_k\) 加 \(1\) ！
一发最优解第二。或许换个PR板子可以更快？
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp(x,y) make_pair(x,y)
#define pb(x) push_back(x)
#define sz(v) (int)v.size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
LL mul(LL x,LL y,LL zky){
LL res=x*y-(long long)((long double)x/zky*y+0.5)*zky;
return res<0?res+zky:res;
}
LL qpow(LL n,LL k,LL zky){
LL res=1;for(;k;k>>=1,n=mul(n,n,zky))if(k&1)res=mul(res,n,zky);
return res;
}
LL gcd(LL x,LL y){
while(y){LL t=x%y;x=y,y=t;}
return x;
}
namespace MR{
static const int P[8]={2,3,5,7,11,19,61,233};
bool mr(LL x,LL p){
if(x%p==0||qpow(p,x-1,x)!=1)return 0;
LL k=x-1;
while(!(k&1)){
LL t=qpow(p,k>>=1,x);
if(t!=1&&t!=x-1)return 0;
if(t==x-1)break;
}
return 1;
}
bool test(LL x){
for(int i=0;i<8;++i){
if(x==P[i])return 1;
if(!mr(x,P[i]))return 0;
}
return 1;
}
}
namespace PR{
LL s[150],d[150];
int tot;
mt19937_64 rnd(679);
#define rad(l,r) (rnd()%((r)-(l)+1)+(l))
LL pr(LL x,LL c){
if(x==4)return 2;
LL v0=c,v1=mul(v0,v0,x)+c,g=1;int len=0;
while(1){
s[++len]=v1-v0,g=mul(g,v1-v0,x);
if(len==127){if(gcd(g,x)>1)break;len=0;}
v0=mul(v0,v0,x)+c,v1=mul(v1,v1,x)+c,v1=mul(v1,v1,x)+c;
}
for(int i=1;i<=len;++i)if((g=gcd(s[i],x))>1)return g;
return x;
}
void solve(LL x){
if(x==1)return;
if(MR::test(x))return d[++tot]=x,void();
LL y=x;while(y==x)y=pr(x,rad(1,x-1));
solve(x/y),solve(y);
}
void work(LL x){tot=0,solve(x),sort(d+1,d+tot+1);}
}
#define zky 1000000007
void fmod(int&x){x-=zky,x+=x>>31&zky;}
inline int qpow(int n,int k){int res=1;for(;k;k>>=1,n=1ll*n*n%zky)if(k&1)res=1ll*n*res%zky;return res;}
int B[3005],fac[3005],ifc[3005],a[3005],inv[3005];
int comb(int n,int m){return n<m?0:1ll*fac[n]*ifc[m]%zky*ifc[n-m]%zky;}
LL n;
int x,y,ans;
void Main(){
ans=0;
scanf("%lld%d%d",&n,&x,&y);
rep(i,0,y)a[y-i+1]=1ll*B[i]*comb(y+1,i)%zky*inv[y+1]%zky;
++a[y];
PR::work(n);
rep(t,1,y+1){
int l=1,r=0,res=1;
while(l<=PR::tot){
while(r<PR::tot&&PR::d[r+1]==PR::d[l])++r;
LL p=PR::d[l];
int k=r-l+1,tmp=0,iv,now,bs,tem;
tem=p%zky,now=qpow(tem,t*(k-1)),iv=qpow(qpow(tem,t),zky-2),bs=qpow(tem,x);
rep(i,0,k-1)fmod(tmp+=now),now=1ll*now*bs%zky*iv%zky;
tmp=1ll*tmp*(zky-qpow(tem,y))%zky;
now=qpow(tem,t*k);
rep(i,0,k)fmod(tmp+=now),now=1ll*now*bs%zky*iv%zky;
res=1ll*res*tmp%zky,l=r+1;
}
fmod(ans+=1ll*a[t]*res%zky);
}
ans=1ll*ans*qpow(n%zky,y)%zky;
printf("%d\n",ans);
}
signed main(){
B[0]=fac[0]=ifc[0]=inv[1]=1;
rep(i,1,3003)fac[i]=1ll*fac[i-1]*i%zky,ifc[i]=qpow(fac[i],zky-2);
rep(i,2,3003)inv[i]=1ll*inv[zky%i]*(zky-zky/i)%zky;
rep(i,1,3003){
int res=0;
rep(j,0,i-1)fmod(res+=1ll*B[j]*comb(i+1,j)%zky);
B[i]=1ll*(zky-res)*qpow(i+1,zky-2)%zky;
}
for(int T=read();T;--T)Main();
}
模数设的好啊！改了个模数就调出来了。
