Luogu P4464 [国家集训队]JZPKIL

Link

\[\begin{aligned} ans&=\sum\limits_{i=1}^n\gcd(i,n)^x\operatorname{lcm}(i,n)^y\\ &=\sum\limits_{i=1}^n(in)^y\gcd(i,n)^{x-y}\\ &=n^y\sum\limits_{d|n}d^x\sum\limits_{t|\frac nd}\mu(t)t^y\sum\limits_{i=1}^{\frac n{dt}}i^y\\ &=\frac{n^y}{y+1}\sum\limits_{d|n}d^x\sum\limits_{t|\frac nd}\mu(t)t^y\sum\limits_{k=0}^y{y+1\choose k}B^+_k(\frac n{dt})^{y-k+1}\\ &=\sum\limits_{k=0}^y{y+1\choose k}B^+_k(id_x\times(\mu\cdot id_y)\times id_{y-k+1})(n) \end{aligned} \]

\(id_x\times(\mu\cdot id_y)\times id_{y-k+1}\)是一个积性函数，因此我们只需要关心\(id_x\times(\mu\cdot id_y)\times id_{y-k+1}(p^e)\)的值。
我们知道\(\forall e\in[2,+\infty),(\mu\cdot id_y)(p^e)=0\)，那么只需要对\((\mu\cdot id_y)(1),(\mu\cdot id_y)(p)\)的情况分类讨论，剩下的暴力卷积即可。

#include<map>
#include<cctype>
#include<cstdio>
#include<random>
#include<cstring>
#include<algorithm>
using i64=long long;
const int N=3007,P=1000000007;
i64 fac[N],ifac[N],B[N];
i64 read(){i64 x=0;char c=getchar();while(isspace(c))c=getchar();while(isdigit(c))(x*=10)+=c&15,c=getchar();return x;}
void inc(i64&a,i64 b){a+=b-P,a+=a>>63&P;}
void dec(i64&a,i64 b){a-=b,a+=a>>63&P;}
i64 pow(i64 a,i64 b){i64 r=1;for(a%=P;b;b>>=1,a=a*a%P)if(b&1)r=r*a%P;return r;}
i64 C(int n,int m){return fac[n]*ifac[m]%P*ifac[n-m]%P;}
void init(int n)
{
    fac[0]=ifac[0]=B[0]=1;
    for(int i=1;i<=n;++i) ifac[i]=pow(fac[i]=i*fac[i-1]%P,P-2);
    for(int i=1;i<=n;++i)
	if(i==1||~i&1)
	{
	    for(int k=0;k<i;++k) dec(B[i],C(i+1,k)*B[k]%P);
	    B[i]=B[i]*pow(C(i+1,i),P-2)%P;
	}
    B[1]=P-B[1];
}
namespace Factoring
{
    std::map<i64,int>cnt;std::mt19937 rng(19260817);
    i64 mul(i64 a,i64 b,i64 p){return (__int128)a*b%p;}
    i64 add(i64 a,i64 b,i64 p){return a+=b-p,a+=a>>63&p;}
    i64 pow(i64 a,i64 b,i64 p){i64 r=1;for(;b;b>>=1,a=mul(a,a,p))if(b&1)r=mul(a,r,p);return r;}
    int Judge(i64 x)
    {
	static int p[7]={2,3,5,7,11,61,24151};
	for(int i=0;i<7;++i) if(x==p[i]) return 1; else if(!(x%p[i])||pow(p[i],x-1,x)!=1) return 0;
	i64 phi=x-1;while(~phi&1) phi>>=1;
	for(int i=0;i<7;++i)
	{
	i64 y=pow(p[i],phi,x);
	if(y==1) continue;
	while(y!=1&&y!=x-1) y=mul(y,y,x);
	if(y==1) return 0;
	}
	return 1;
    }
    i64 Find(i64 x)
    {
	i64 c=rng()%(x-1)+1,las=0,now=0;
	for(int len=1;;las=now,len<<=1)
	{
	    i64 prod=1,d;
	    for(int i=1;i<=len;++i) if(now=mul(now,now,x)+c,prod=mul(prod,llabs(now-las),x),!(i&127)&&(d=std::__gcd(prod,x))!=1) return d;
	    if((d=std::__gcd(prod,x))!=1) return d;
	}
    }
    void Factor(i64 x,int f)
    {
	if(x==1) return ;
	if(Judge(x)) return cnt[x]+=f,void();
	i64 p=Find(x);int c=0;
	while(!(x%p))x/=p,++c;
	Factor(p,f*c),Factor(x,f);
    }
}using Factoring::cnt;
void solve()
{
    static i64 f[N],g[N];
    i64 n=read(),x=read(),y=read(),ans=0;
    cnt.clear(),Factoring::Factor(n,1);
    for(int k=0;k<=y;++k)
    {
	if(!B[k]) continue;
	i64 t=C(y+1,k)*B[k]%P,z=y+1-k;
	for(auto[p,e]:cnt)
	{
	    i64 t1=0,t2=0,u=f[1]=pow(p,x),v=g[1]=pow(p,z);f[0]=g[0]=1;
	    for(int j=2;j<=e;++j) f[j]=u*f[j-1]%P,g[j]=v*g[j-1]%P;
	    for(int j=0;j<=e;++j) (t1+=f[j]*g[e-j])%=P;
	    for(int j=0;j<e;++j) (t2+=f[j]*g[e-j-1])%=P;
	    (t2*=P-pow(p,y))%=P,(t*=(t1+t2))%=P;
	}
	inc(ans,t);
    }
    printf("%lld\n",ans*pow(n,y)%P*pow(y+1,P-2)%P);
}
int main()
{
    init(3001);
    for(int t=read();t;--t) solve();
}