P10182 题解

· · 题解

对于 n=\prod p_i^{c_i}，有：

\sigma_0^t(n^s)=\prod(1+s\times c_i)^{t}

然后设：

f(n)=\sum_{d\mid n} \mu\left(\dfrac{n}{d}\right)\sigma_0^t(n^s)

所以有：

\sigma_0^t(n^s)=\sum_{d\mid n}f(d)

由于 \sigma_0^t 是积性函数，不妨先假设 n=p^c。则：

\sigma_0^t(p^{cs})=\sum_{i=0}^{c}f(p^i)

也就是说对于 n=p^c，f(p^c) 是 \sigma_0^t 关于 c 的高维差分。

类似的，设：

f(n)=\sum_{d\mid n} g(d)

那个 g 是 f 关于 c 的差分，即 g 是 \sigma_0^t 的二阶高维差分。具体地就有：

g(p^c)=(1+c_is)^t-[c_i\ge 1]2(1+(c_i-1)s)^t+[c_i\ge 2] (1+(c_i-2)s)^t

那么 g 是积性函数，就有：

g(n)=\prod((1+c_is)^t-[c_i\ge 1]2(1+(c_i-1)s)^t+[c_i\ge 2] (1+(c_i-2)s)^t)

矩阵可以拆成：

\begin{aligned} a_{i,j}&=\sum_{d\mid \gcd(i,j)}g(d)=\sum_{d=1}^{n}[d\mid i][d\mid j]g(d)\\\\ b_{i,j}&=[j\mid i]g(j)\\\\ c_{i,j}&=[i\mid j]\\\\ A&=B*C \end{aligned} - $\text{det}(B)=\prod_{i=1}^n b_{i,i}=\prod_{i=1}^ng(i)

故答案为：

\det(A)=\text{det} (B)\times \det(C)=\prod_{m=1}^{n}\prod((1+c_is)^t-[c_i\ge 1]2(1+(c_i-1)s)^t+[c_i\ge 2] (1+(c_i-2)s)^t)

依次考虑每个质数的贡献。

对于 p\le \sqrt{n} 的部分暴力枚举 p 和 c，暴力求快速幂，复杂度 \mathcal{O}(\sqrt{n}\log t)。

对于 p>\sqrt{n}，可知 c\le 1，当 c=0 时贡献为 1 不用考虑，当 c=1 时 g(p)=(1+s)^t-2，故只需要筛素数个数即可，在本题范围内，使用 min25 筛即可。

#include<bits/stdc++.h>

#define ll long long
#define mk make_pair
#define fi first
#define se second

using namespace std;

const int mod=1e9+7;
int ksm(int x,ll y,int p=mod){
    int ans=1;y%=(p-1);
    for(int i=y;i;i>>=1,x=1ll*x*x%p)if(i&1)ans=1ll*ans*x%p;
    return ans%p;
}
int inv(int x,int p=mod){return ksm(x,p-2,p)%p;}
mt19937 rnd(time(0));
int randint(int l,int r){return rnd()%(r-l+1)+l;}
void add(int &x,int v){x+=v;if(x>=mod)x-=mod;}
void Mod(int &x){if(x>=mod)x-=mod;}
int cmod(int x){if(x>=mod)x-=mod;return x;}

const int N=1e6+5;
const int B=320000;
ll n;
int cp,pri[N],s,t;
bool vis[N];
void Euler(int V){
    for(int i=2;i<=V;i++){
        if(!vis[i])pri[++cp]=i;
        for(int j=1;j<=cp&&pri[j]*i<=V;j++){
            vis[i*pri[j]]=1;
            if(i%pri[j]==0)break;
        }
    }
}
ll f[N<<1],num[N<<1];
int id1[N],id2[N];
int getid(ll x){
    if(x>B)return id2[n/x];
    return id1[x];
}
void calcpri(){
    Euler(B);
    int ncnt=0;
    for(ll l=1,r=0;l<=n;l=r+1){
        ll val=(n/l);
        if(val<=B)num[id1[val]=++ncnt]=val;
        else num[id2[n/val]=++ncnt]=val;
        r=n/(n/l);r=min(r,n);
    }
    for(int i=1;i<=ncnt;i++)f[i]=num[i]-1;
    for(int i=1;i<=cp;i++){
        for(int j=1;j<=ncnt&&num[j]>=1ll*pri[i]*pri[i];j++){
            int k=getid(num[j]/pri[i]);
            f[j]-=f[k]-(i-1);
        }
    }
}

int g[100];
void solve(){
    cin>>n>>s>>t;
    if(s==0)return puts(n==1?"1":"0"),void();

    int ans=1;calcpri();
    g[1]=cmod(ksm(s+1,t)+mod-2);
    for(int i=2;i<=40;i++){
        add(g[i],ksm(1ll*i*s%mod+1,t));
        add(g[i],mod-2ll*ksm(1ll*(i-1)*s%mod+1,t)%mod);
        add(g[i],ksm(1ll*(i-2)*s%mod+1,t));
    }
    int lim=n/B;
    for(int i=1;i<=cp&&n/pri[i]>lim;i++){
        ll val=pri[i];int k=1;
        while(val<=n){
            ll cc=n/val-n/(val*pri[i]);
            ans=1ll*ans*ksm(g[k],cc)%mod,val*=pri[i],k++;
        }
    }
    ll cc=0;
    // floor(n/p) >= i <==> i*p <= n <==> p <= floor(n/i)
    // floor(n/p) == i <==> floor(n/(i+1)) < p <= floor(n/i)
    for(int i=1;i<=lim;i++){ // floor(n/p) = i
        ll cnt=f[getid(n/i)]-f[getid(n/(i+1))];
        cc+=1ll*i*cnt;
    }
    ans=1ll*ans*ksm(g[1],cc)%mod;
    cout<<ans<<endl;
}

signed main(void){

    solve();

    return 0;
}