洛谷 P10182 一径入繁华

题目背景

伴随龙年到来的，还有帆巨很喜欢的九省联考。为了爆踩压轴题。帆巨狠狠地重温了数论。

数论所生，繁华之地！

题目描述

帆巨觉得求 ��xa 在  mod  �mod p 意义下的值太简单了，所以他想求 �0�(��)σ0s​(xt) 在  mod  �mod p 意义下的值。

帆帆不满足于只计算一次，于是他列了一个 �×�n×n 的数表 �A，保证第 �i 行第 �j 列（1≤�,�≤�1≤i,j≤n）中的元素 ��,�ai,j​ 满足：

��,�=∑�∣gcd⁡(�,�)�(gcd⁡(�,�)�)×(�0(��))�ai,j​=d∣gcd(i,j)∑​μ(dgcd(i,j)​)×(σ0​(ds))t

帆帆想知道这个数表长什么样子，但这个数表实在太大了，所以请你告诉他 det⁡�detA 对 109+7109+7 取模后的结果。

注释：

1. 表达式中的各种函数含义在 这里（�μ 表示莫比乌斯函数，�0σ0​ 表示约数个数函数）。
2. det⁡�detA 表示方阵 �A 的 行列式。

输入格式

一行，输入三个整数 �,�,�n,s,t。

输出格式

一行，输出一个整数表示答案。

输入输出样例

输入 #1复制

2 1 2

输出 #1复制

2

输入 #2复制

2 3 4

输出 #2复制

254

输入 #3复制

19 8 10

输出 #3复制

913255725

输入 #4复制

10000000000 1 2

输出 #4复制

880793261

说明/提示

【样例 11 解释】

矩阵 �A 如下：

[1113][11​13​]

行列式为 1×3−1×1=21×3−1×1=2。

【样例 22 解释】

矩阵 �A 如下：

[111255][11​1255​]

行列式为 1×255−1×1=2541×255−1×1=254。

数据范围

本题采用 子任务捆绑测试。

对于 100%100% 的数据，保证 1≤�≤10111≤n≤1011，0≤�,�<109+70≤s,t<109+7。

子任务编号	�n	特殊性质	分值
Subtask #1	≤500≤500	无	88
Subtask #2	≤107≤107	�=1,�=2s=1,t=2	55
Subtask #3	≤107≤107	�=1s=1	1010
Subtask #4	≤1011≤1011	�=1,�=2s=1,t=2	1010
Subtask #5	≤1011≤1011	�=1s=1	1010
Subtask #6	≤1011≤1011	�=1t=1	22
Subtask #7	≤107≤107	�≤9t≤9	1010
Subtask #8	≤1011≤1011	�≤9t≤9	1515
Subtask #9	≤107≤107	无	1010
Subtask #10	≤1011≤1011	无	2020

特殊性质 一栏为空则表示没有特殊性质。子任务中没有规定范围的变量的值均在 [0,109+7)[0,109+7) 范围内生成。

时间限制：2000 ms2000 ms；

空间限制：512 MB512 MB。

对于 �=∏����n=∏pici​​，有：

�0�(��)=∏(1+�×��)�σ0t​(ns)=∏(1+s×ci​)t

然后设：

�(�)=∑�∣��(��)�0�(��)f(n)=d∣n∑​μ(dn​)σ0t​(ns)

所以有：

�0�(��)=∑�∣��(�)σ0t​(ns)=d∣n∑​f(d)

由于 �0�σ0t​ 是积性函数，不妨先假设 �=��n=pc。则：

�0�(���)=∑�=0��(��)σ0t​(pcs)=i=0∑c​f(pi)

也就是说对于 �=��n=pc，�(��)f(pc) 是 �0�σ0t​ 关于 �c 的高维差分。

类似的，设：

�(�)=∑�∣��(�)f(n)=d∣n∑​g(d)

那个 �g 是 �f 关于 �c 的差分，即 �g 是 �0�σ0t​ 的二阶高维差分。具体地就有：

�(��)=(1+���)�−[��≥1]2(1+(��−1)�)�+[��≥2](1+(��−2)�)�g(pc)=(1+ci​s)t−[ci​≥1]2(1+(ci​−1)s)t+[ci​≥2](1+(ci​−2)s)t

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

�(�)=∏((1+���)�−[��≥1]2(1+(��−1)�)�+[��≥2](1+(��−2)�)�)g(n)=∏((1+ci​s)t−[ci​≥1]2(1+(ci​−1)s)t+[ci​≥2](1+(ci​−2)s)t)

矩阵可以拆成：

��,�=∑�∣gcd⁡(�,�)�(�)=∑�=1�[�∣�][�∣�]�(�)��,�=[�∣�]�(�)��,�=[�∣�]�=�∗�ai,j​bi,j​ci,j​A​=d∣gcd(i,j)∑​g(d)=d=1∑n​[d∣i][d∣j]g(d)=[j∣i]g(j)=[i∣j]=B∗C​

∗∗ 是矩阵乘法。�,�B,C 都是上、下三角矩阵，于是

• det(�)=∏�=1���,�=∏�=1��(�)det(B)=∏i=1n​bi,i​=∏i=1n​g(i)
• det(�)=∏�=1���,�=1det(C)=∏i=1n​ci,i​=1

故答案为：

det⁡(�)=det(�)×det⁡(�)=∏�=1�∏((1+���)�−[��≥1]2(1+(��−1)�)�+[��≥2](1+(��−2)�)�)det(A)=det(B)×det(C)=m=1∏n​∏((1+ci​s)t−[ci​≥1]2(1+(ci​−1)s)t+[ci​≥2](1+(ci​−2)s)t)

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

对于 �≤�p≤n​ 的部分暴力枚举 �p 和 �c，暴力求快速幂，复杂度 �(�log⁡�)O(n​logt)。

对于 �>�p>n​，可知 �≤1c≤1，当 �=0c=0 时贡献为 11 不用考虑，当 �=1c=1 时 �(�)=(1+�)�−2g(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;
}

拜拜!