luogu4464：莫比乌斯反演，积性函数和伯努利数

题目链接：https://www.luogu.com.cn/problem/P4464

 

简记$gcd(x,y)=(x,y)$。

推式子：

$\sum_{i=1}^{n}{(i,n)^xlcm(i,n)^y}$

$=\sum_{i=1}^{n}{(i,n)^{x-y}(in)^y}$

$=n^y\sum_{d|n}d^{x-y}\sum_{i}i^y[(i,n)=d]$

$=n^y\sum_{d|n}{d^{x-y}\sum_{i=1}^{\frac{n}{d}}{(id)^y[(i,\frac{n}{d})=1]}}$

$=n^y\sum_{d|n}{d^x\sum_{i=1}^{\frac{n}{d}}{i^y[(i,\frac{n}{d})=1]}}$

$=n^y\sum_{d|n}{d^x\sum_{i=1}^{\frac{n}{d}}{i^y\sum_{k|(i,\frac{n}{d})}{\mu(k)}}}$

改变最后的求和式的顺序，注意k的取值要整除$\frac{n}{d}$。

$=n^y\sum_{d|n}{d^x\sum_{k|\frac{n}{d}}{\mu(k)}\sum_{i=1}^{\frac{n}{dk}}{(ik)^y}}$

$=n^y\sum_{d|n}{d^x\sum_{k|\frac{n}{d}}{\mu(k)k^y}\sum_{i=1}^{\frac{n}{dk}}{i^y}}$

若我们把最后的和式看成是一个y+1次的多项式，它的第i项系数为$t_i$，那么

$=n^y\sum_{j=1}^{y+1}{t_i}\sum_{d|n}{d^x\sum_{k|\frac{n}{d}}{\mu(k)k^y}(\frac{n}{dk})^j}$

（注意最后一项的变化）

注意到后面所有的形式都是积性函数的卷积，因此我们尝试一个质数的若干次方的贡献，最后将其相乘就得到某一个对应的j的答案。

于是对于固定的$d=p^q,n=p^w$，最后一个和式的贡献为$(p^{w-q})^j-p^y(p^{w-q-1})^j$（分别考虑$d=1$和$d=p$）。特别地，若$d=n$，那么贡献为1。如果我们能快速得到所有$t_i$，那么就能在$O(n^{0.25}logn+ylog^3n)$的复杂度完成一次询问。

考虑如何得到$t_i$。一种方法是使用插值，这里讨论使用伯努利数的方法（应用，证明看https://www.luogu.com.cn/blog/chihik/bo-nu-li-shu）。

 

递归定义：$\sum_{i=0}^{n}{\tbinom{n+1}{i}B_i=[n=0]}$

令$S_m(n)=\sum_{i=0}^{n-1}{i^m}$，那么$S_m(n)=\frac{1}{m+1}\sum_{i=0}^{m}{\tbinom{m+1}{i}B_in^{m+1-i}}$

于是对于上面的y+1次多项式，有$S_y(\frac{n}{dk})=\sum_{i=0}^{\frac{n}{dk}-1}{i^y}=\frac{1}{y+1}\sum_{i=0}^{y}{\tbinom{y+1}{i}{B_i(\frac{n}{dk})^{y+1-i}}}$，所以$t_{y+1-i}=\frac{1}{y+1}\tbinom{y+1}{i}B_i$。（它的零次项系数为0！（这是感叹号，不是阶乘））

但我们注意到少了一项$(\frac{n}{dk})^y$，因此$t_y$需要增加1（对应了$i=1$）。

 

 

#include<bits/stdc++.h>
using namespace std;
typedef long long int ll;
typedef long double ld;
namespace MATH
{
    inline ll mul(ll x,ll y,ll mod)
    {
        return ((x*y-(ll)((ld)x*y/mod)*mod)%mod+mod)%mod;
    }
    inline ll qpow(ll x,ll y,ll mod)
    {
        ll ans=1,base=x;
        while(y)
        {
            if(y&1)
                ans=mul(ans,base,mod);
            base=mul(base,base,mod);
            y>>=1; 
        }
        return ans;
    }
    const int prime[10]={2,3,5,7,11,13,17,19,23,29};
    const int LEN=10;
    inline bool miller(ll x)
    {
        for(int i=0;i<LEN;++i)
        {
            if(x<prime[i])
                return false;
            else if(x==prime[i])
                return true;
            ll now=x-1,y=qpow(prime[i],x-1,x);
            if(y!=1)
                return false;
            while(now%2==0&&y==1)
            {
                now>>=1;
                y=qpow(prime[i],now,x);
                if(y!=1&&y!=x-1)
                    return false;
            }
        }
        return true;
    }
    ll gcd(ll x,ll y)
    {
        return x%y==0?y:gcd(y,x%y);
    }
    inline ll f(ll x,ll y,ll mod)
    {
        return (mul(x,x,mod)+y)%mod;
    }
    inline ll get(ll n,ll c,int steps)
    {
        if(!(n&1))
            return 2;
        ll x=2,y=2,d=1;
        while(true)
        {
            ll tmpx=x,tmpy=y;
            for(int i=1;i<=steps;++i)
            {
                x=f(x,c,n);
                y=f(y,c,n);
                y=f(y,c,n);
                d=mul(d,((y-x)%n+n)%n,n);
            }
            d=gcd(d,n);
            if(d==1)
                continue;
            if(d!=n)
                return d;
            x=tmpx,y=tmpy;
            for(int i=1;i<=steps;++i)
            {
                x=f(x,c,n);
                y=f(y,c,n);
                y=f(y,c,n);
                d=gcd(((y-x)%n+n)%n,n);
                if(d!=1&&d!=n)
                    return d;
            }
            return 0;
        }
    }
    vector<ll>wait;
    void pollard(ll n)
    {
        if(miller(n))
        {
            wait.push_back(n);
            return;
        }
        ll now=0,c=1,g=pow(n,0.1)+1;
        while(!now)
            now=get(n,++c,g);
        pollard(now),pollard(n/now);
    }
}
const ll mod=1000000007;
ll C[3005][3005],B[3005];
inline ll qpow(ll x,ll y)
{
    ll ans=1,base=x%mod;
    while(y)
    {
        if(y&1)
            ans=ans*base%mod;
        base=base*base%mod;
        y>>=1;
    }
    return ans;
}
inline void init()
{
    for(int i=0;i<=3004;++i)
    {
        C[i][0]=1;
        for(int j=1;j<=i;++j)
            C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
    }
    B[0]=1;
    for(int i=1;i<=3002;++i)
    {
        for(int j=0;j<i;++j)
            B[i]=(B[i]+B[j]*C[i+1][j])%mod;
        B[i]=(B[i]*qpow(-i-1,mod-2)%mod+mod)%mod;
    }
}
ll tmp[3005];
inline void solve()
{
    ll n,x,y;
    cin>>n>>x>>y;
    MATH::wait.clear();
    MATH::pollard(n);
    vector<ll>wait=MATH::wait;
    map<ll,int>vis;
    for(ll x:wait)
        ++vis[x];
    ll ans=0;
    for(int i=0;i<=y;++i)
        tmp[y+1-i]=C[y+1][i]*B[i]%mod*qpow(y+1,mod-2)%mod;
    tmp[y]+=1;
    for(int i=1;i<=y+1;++i)
    {
        ll s=1;
        for(auto A:vis)
        {
            ll p=A.first,k=A.second;
            ll g=0;
            ll now=1,delta=qpow(p,x);
            ll G=qpow(p,y);
            for(int j=0;j<k;++j,now=now*delta%mod)
                g=(g+now*(qpow(qpow(p,k-j),i)-G*qpow(qpow(p,k-j-1),i)%mod)%mod)%mod;
            g=(g+now%mod)%mod;
            s=s*g%mod;
        }
        ans=(ans+tmp[i]*s)%mod;
    }
    ans=ans*qpow(n%mod,y)%mod;
    ans=(ans%mod+mod)%mod;
    cout<<ans<<endl;
}
int main()
{
    ios::sync_with_stdio(false);
    init();
    int T;
    cin>>T;
    while(T--)
        solve();
    return 0;
}