题目链接
洛谷：https://www.luogu.org/problemnew/show/P4464
Solution
这题是真的毒....数论大杂烩，窝断断续续写了两天。
众所周知：
[{
m lcm}(x,y)=frac{xy}{gcd(x,y)}
]
带进去，顺便枚举(gcd)值套(mobius)反演公式然后换一下求和符号：
[egin{align}
ans=&n^ysum_{d|n}d^{x-y}sum_{i=1}^{n}i^y[gcd(i,n)=1]\
=&n^ysum_{d|n}d^xsum_{i=1}^{n/d}i^y[gcd(id,n)=1]\
=&n^ysum_{d|n}d^xsum_{i=1}^{n/d}i^ysum_{t|i,t|frac{n}{d}}mu(t)\
=&n^ysum_{d|n}d^xsum_{t|frac{n}{d}}mu(t)t^ysum_{i=1}^{n/dt}i^y\
end{align}
]
众所周知，自然数(k)次幂和可以表示为一个(k+1)次多项式，设多项式系数为(s)，即：
[sum_{i=1}^{n}i^y=sum_{i=0}^{y+1}s_in^i
]
那么带进去顺便把(s_i)丢最前面：
[egin{align}
ans=&n^ysum_{d|n}d^xsum_{t|frac{n}{d}}mu(t)t^ysum_{i=1}^{y+1}s_i(frac{n}{dt})^i\
=&n^ysum_{i=1}^{y+1}s_isum_{d|n}d^{x}sum_{t|frac{n}{d}}mu(t)t^{y}(frac{n}{dt})^i\
end{align}
]
可以注意到后面是一个狄利克雷卷积的形式，设：
[f_i(n)=sum_{d|n}d^{x}sum_{t|frac{n}{d}}mu(t)t^{y}(frac{n}{dt})^i\
a(n)=d^n,b(n)=mu(n)n^y,c_i(n)=n^i
]
显然：
[f_i=a*b*c_i
]
也就是说(f_i)是一个积性函数，我们只需要关系(f(p^w))的值就好了，然后用( m pollard rho)分解(n)，暴力算就好了。
显然根据(mu)函数的性质，(f_i)的定义式中枚举的(t)只有等于(1 or p)的时候才有贡献，我们分类讨论：
-
(t=1)：
[f(n)=n^isum_{d|n}d^{x-i}\ f(p^w)=p^{wi}sum_{j=0}^{w}p^{j(x-i)} ] -
(t=p)：
[f(n)=-n^isum_{d|n}d^{x-i}p^{y-i}\f(p^w)=-p^{wi+y-i}sum_{j=0}^{w-1}p^{j(x-i)} ]
这里可以用等比数列求和，当然暴力枚举复杂度也是对的，本题复杂度瓶颈在质因数分解上，但是这题卡pollard rho常就很不友好。
那个自然数幂和的系数可以用伯努利数求，具体上百度，伯努利数直接大力预处理就好了。
然后就做完了，复杂度(O(Tsqrt[4]{n}+Tylog n+y^2))。
// #pragma GCC optimize(3)
#include<bits/stdc++.h>
using namespace std;
template<typename T> void read(T &x) {
x=0;T f=1;char ch=getchar();
for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
void print(int x) {
if(x<0) putchar('-'),x=-x;
if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}
#define lf double
#define ll long long
#define pii pair<int,int >
#define vec vector<int >
#define pb push_back
#define mp make_pair
#define fr first
#define sc second
#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++)
const int maxn = 4e3+10;
const int N = 3e3+10;
const int inf = 1e9;
const lf eps = 1e-8;
ll r[maxn],cnt;
namespace Pollard_rho {
const int pri[] = {0,2,3,5,7,11,13,23,31};
ll gcd(ll a,ll b) {return !b?a:gcd(b,a%b);}
ll mul(ll x,ll y,ll p) {
ll t=(long double)x*y/p;
ll res=x*y-t*p;res%=p;if(res<0) res+=p;
return res;
}
ll qpow(ll a,ll x,ll p) {
ll res=1;
for(;x;x>>=1,a=mul(a,a,p)) if(x&1) res=mul(res,a,p);
return res;
}
bool MR(ll n) {
if(n<2) return 1;
ll e=(n-1)>>__builtin_ctz(n-1);
for(int i=1;i<=8;i++) {
if(pri[i]>=n) break;
for(ll d=e,lst=qpow(pri[i],d,n),now;d!=n-1;d<<=1,lst=now) {
now=mul(lst,lst,n);
if(lst!=1&&lst!=n-1&&now==1) return 0;
}if(qpow(pri[i],n-1,n)!=1) return 0;
}return 1;
}
ll calc(ll n,ll v) {
ll a=rand()%n,b=a,g=1;
while(g==1) {
a=(mul(a,a,n)+v)%n;
b=(mul(b,b,n)+v)%n;
b=(mul(b,b,n)+v)%n;
g=gcd(abs(a-b),n);
}return g;
}
void PR(ll n) {
if(MR(n)) return n<2?0:r[++cnt]=n,void();
ll d=n;
while(d==n) d=calc(n,1ull*rand()*rand()*rand()*rand()%(n-1)+1);
PR(d),PR(n/d);
}
void divide(ll n) {
for(int i=1;i<=8;i++) while(n%pri[i]==0) n/=pri[i],r[++cnt]=pri[i];
PR(n);
}
}
const int mod = 1e9+7;
int fac[maxn],ifac[maxn],inv[maxn],b[maxn];
int add(int x,int y) {return x+y>=mod?x+y-mod:x+y;}
int del(int x,int y) {return x-y<0?x-y+mod:x-y;}
int mul(int x,int y) {return 1ll*x*y-1ll*x*y/mod*mod;}
int c(int x,int y) {return mul(fac[x],mul(ifac[y],ifac[x-y]));}
int qpow(int a,int x) {
int res=1,f=0;if(x<0) x=-x,f=1;
for(;x;x>>=1,a=mul(a,a)) if(x&1) res=mul(res,a);
return f?qpow(res,mod-2):res;
}
namespace Bernoulli {
void sieve() {
fac[0]=ifac[0]=inv[0]=inv[1]=1;
for(int i=1;i<maxn;i++) fac[i]=mul(fac[i-1],i);
for(int i=2;i<maxn;i++) inv[i]=mul(mod-mod/i,inv[mod%i]);
for(int i=1;i<maxn;i++) ifac[i]=mul(ifac[i-1],inv[i]);
}
void get_B() {
b[0]=1;
for(int i=1;i<N;i++) {
for(int j=0;j<i;j++) b[i]=add(b[i],mul(b[j],c(i+1,j)));
b[i]=del(0,mul(b[i],inv[i+1]));
}b[1]++;
}
}
int x,y,s[maxn],tot;ll n;
struct data {ll p;int x;}a[maxn];
int calc(int p,int w,int i) {
int res=0;
for(int j=0;j<=w;j++) res=add(res,qpow(p,j*(x-i)+w*i));
for(int j=0;j<=w-1;j++) res=del(res,qpow(p,j*(x-i)+w*i-i+y));
return res;
}
void solve() {
read(n),read(x),read(y);cnt=0;s[0]=0;tot=0;
if(!n) return puts("0"),void();
for(int i=0;i<=y;i++) s[y+1-i]=mul(inv[y+1],mul(b[i],c(y+1,i)));
Pollard_rho :: divide(n);
sort(r+1,r+cnt+1);
for(int i=1;i<=cnt;i++)
if(r[i]!=r[i-1]) a[++tot].p=r[i],a[tot].x=1;
else a[tot].x++;
int ans=0;
for(int i=0;i<=y+1;i++) {
int res=s[i];
for(int j=1;j<=tot;j++) {
int p=a[j].p%mod,w=a[j].x;
res=mul(res,calc(p,w,i));
}ans=add(ans,res);
}ans=mul(ans,qpow(n%mod,y));
write(ans);
for(int i=0;i<=cnt+5;i++) a[i].p=a[i].x=r[i]=0;
for(int i=0;i<=y+1+5;i++) s[i]=0;
}
int main() {
srand(19280817);
Bernoulli :: sieve(),Bernoulli :: get_B();
int t;read(t);while(t--) solve();
return 0;
}