zcmimi's blog
查看原题

点击跳转

\begin{aligned} &\quad\sum_{i=1}^n \gcd(i,n)^x lcm(i,n)^y \\ &=\sum_{i=1}^n \frac{in}{lcm(i,n)}^x lcm(i,n)^y \\ &=\sum_{i=1}^n (in)^x lcm(i,n)^{y-x} \\ &=\sum_{i=1}^n (in)^x \frac{in}{\gcd(i,n)}^{y-x}\\ &=\sum_{i=1}^n (in)^y \gcd(i,n)^{x-y} \\ &=n^y\sum_{d|n} d^{x-y}\sum_{i=1}^n i^y [\gcd(i,n)=d] \\ &=n^y\sum_{d|n} d^{x-y}\sum_{i=1}^{n/d} (id)^y [\gcd(i,n)=1] \\ &=n^y\sum_{d|n} d^{x-y}\sum_{i=1}^{n/d} (id)^y \sum_{k|i,k|n}\mu(k) \\ &=n^y\sum_{d|n} d^{x-y}\sum_{k|\frac nd}\mu(k)\sum_{i=1}^{n/dk} (idk)^y \\ &=n^y\sum_{d|n} d^x\sum_{k|\frac nd}\mu(k)k^y\sum_{i=1}^{n/dk} i^y \end{aligned}

\displaystyle S_k(n)=\sum_{i=1}^{n-1}i^k,可以看做一个y+1次多项式,假设求出后每一项系数为t_i

=n^y\sum_{d|n} d^x\sum_{k|\frac nd}\mu(k)k^yS_y(\frac n{dk}) \\ =n^y\sum_{d|n} d^x\sum_{k|\frac nd}\mu(k)k^y\sum_{i=0}^{y+1}t_i(\frac n{dk})^i \\ =n^y\sum_{i=0}^{y+1}t_i \sum_{d|n} d^x\sum_{k|\frac nd}\mu(k)k^y (\frac n{dk})^i

其中,设\displaystyle Z(n)=\sum_{d|n} d^x\sum_{k|\frac nd}\mu(k)k^y (\frac n{dk})^i ,

可以看成f(d)=d^x,g(d)=\mu(d)d^y,h^i(d)=d^i这三个积性函数的狄利克雷卷积

卷积后还是积性函数,且狄利克雷卷积满足交换律,只有当x=1x=p^k\mu(x)\ne 0

\begin{cases} Z(1)=1\\ Z(p)=p^x-p^y+p^i\\ Z(p^k)=(f * h^i)(p^k)-(f * h^i)(p^{k-1}) \cdot p^y \end{cases}

其中\displaystyle (f * h^i)(p^k)=\sum_{i+j=k} f(p^i)\cdot h^i(p^j)

使用Pollard-Rhon分解质因数,枚举每个质因子(包括次幂),然后狄利克雷卷积计算答案

#include<bits/stdc++.h>
const int N=3011,P=1000000007;
typedef __int128 ll;
typedef __int128 i128;
template<typename T>void rd(T&x){x=0;char c;bool f=0;for(c=getchar();c<'0'||'9'<c;c=getchar())f^=c=='-';for(x=c-48,c=getchar();'0'<=c&&c<='9';x=x*10+c-48,c=getchar());x=f?-x:x;}
ll rnd(){
    static unsigned int seed=233;
    return seed=seed*482711;
}
ll pw(i128 x,ll b,ll p=P){
    i128 res=1;
    while(b){
        if(b&1)(res*=x)%=p;
        b>>=1;(x*=x)%=p;
    }
    return res;
}
void mod(ll&x){if(x>=P)x-=P;}
bool chk(const ll&p){
    if(p==2||p==3||p==5)return 1;
    if(!(p&1)||p<2)return 0;
    ll t=p-1;int k=0;
    while(!(t&1))t>>=1,++k;
    static int test[10]={2,3,5,7,11,13,17,19,23,29};
    for(int i=0;i<10&&test[i]<p;++i){
        ll a=pw(test[i],t,p),nxt;
        for(int j=1;j<=k;++j){
            nxt=(i128)a*a%p;
            if(nxt==1&&a!=1&&a!=p-1)return 0;
            a=nxt;
        }
        if(a!=1)return 0;
    }
    return 1;
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll abs(const ll&x){return x<0?-x:x;}
ll pr(ll p){
    ll k=2,x=rnd()%(p-1)+1,y=x,d=1,c=rnd()%99997+1;
    for(int i=1;d==1;++i){
        x=((i128)x*x+c)%p;
        d=gcd(abs(x-y),p);
        if(i==k)y=x,k<<=1;
    }
    return d;
}
ll a[N];int cnt,w[N];
void fact(ll x){
    if(x<2)return;
    if(chk(x)){a[++cnt]=x;return;}
    ll p=x;
    while(p>=x)p=pr(x);
    while(x%p==0)x/=p;
    fact(p),fact(x);
}
ll fac[N],ifac[N],B[N],b[N];
void init(int n){
    fac[0]=1;
    for(int i=1;i<=n;++i)
        fac[i]=fac[i-1]*i%P;
    ifac[n]=pw(fac[n],P-2);
    for(int i=n;i;--i)
        ifac[i-1]=ifac[i]*i%P;
    B[0]=1;
    for(int i=1;i<=n;++i)
        for(int j=0;j<i;++j)
            mod(B[i]+=P-fac[i]*ifac[j]%P*ifac[i+1-j]%P*B[j]%P);
    ++B[1];
}
ll f[N],g[N],F[N],H[N],h[N];
void solve(const ll&n,const int&x,const int&y){
    b[0]=0;
    for(int i=0;i<=y;++i)
        b[y+1-i]=fac[y]*ifac[i]%P*ifac[y+1-i]%P*B[i]%P;
    cnt=0,fact(n),std::sort(a+1,a+cnt+1);
    ll t=n;
    for(int j=1;j<=cnt;++j)
        while(t%a[j]==0)t/=a[j],++w[j];
    ll ans=0;
    for(int j=1;j<=cnt;++j)
        f[j]=pw(a[j],x),g[j]=pw(a[j],y);
    for(int i=0;i<=y+1;++i){
        ll res=1;
        for(int j=1;j<=cnt;++j){
            h[j]=i?h[j]*a[j]%P:1;
            F[0]=H[0]=1;
            for(int k=1;k<=w[j];++k)
                F[k]=F[k-1]*f[j]%P,
                H[k]=H[k-1]*h[j]%P;
            ll FH=0,fh=0;
            for(int k=0;k<=w[j];++k)
                mod(FH+=F[k]*H[w[j]-k]%P);
            for(int k=0;k<w[j];++k)
                mod(fh+=F[k]*H[w[j]-1-k]%P);
            res=res*(FH-fh*g[j]%P+P)%P;
        }
        mod(ans+=res*b[i]%P);
    }
    for(int j=1;j<=cnt;++j)a[j]=w[j]=0;
    printf("%d\n",pw(n,y)*ans%P);
}
int main(){
    int T,x,y;ll n;
    rd(T);
    init(3002);
    while(T--)
        rd(n),rd(x),rd(y),
        solve(n,x,y);
}
LG 4464 [国家集训队]JZPKIL
comment评论
Search
search