简要题意
给出 \(N\)，计算：
\[\prod_{i=1}^N\prod_{j=1}^N\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}\pmod{104857601}
\]
\]
\(1 \leq N \leq 10^6\)。时间限制 \(200\operatorname{ms}\)，空间限制 \(7.81\operatorname{MB}\)。
思路
又是一道没有看题解做出来的反演题。
消灭莫比乌斯反演暴政！世界属于欧拉反演！
首先给出 \(\operatorname{lcm}\) 的定义：
\[\operatorname{lcm}(i,j)=\dfrac{ij}{\gcd(i,j)}
\]
\]
这个式子非常好用。
然后就可以愉快地推式子了：
\[\begin{aligned}
&\prod_{i=1}^N\prod_{j=1}^N\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}\\
&=\prod_{i=1}^N\prod_{j=1}^N\frac{ij}{\gcd(i,j)^2}\\
&=\dfrac{\prod_{i=1}^N\prod_{j=1}^N{ij}}{\prod_{i=1}^N\prod_{j=1}^N\gcd(i,j)^2}\\
&=\dfrac{\prod_{i=1}^N\prod_{j=1}^N{ij}}{(\prod_{i=1}^N\prod_{j=1}^N\gcd(i,j))^2}
\end{aligned}
\]
&\prod_{i=1}^N\prod_{j=1}^N\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}\\
&=\prod_{i=1}^N\prod_{j=1}^N\frac{ij}{\gcd(i,j)^2}\\
&=\dfrac{\prod_{i=1}^N\prod_{j=1}^N{ij}}{\prod_{i=1}^N\prod_{j=1}^N\gcd(i,j)^2}\\
&=\dfrac{\prod_{i=1}^N\prod_{j=1}^N{ij}}{(\prod_{i=1}^N\prod_{j=1}^N\gcd(i,j))^2}
\end{aligned}
\]
把分子分母拆开计算。先算分子：
\[\begin{aligned}
&\prod_{i=1}^N\prod_{j=1}^N{ij}\\
&=\prod_{i=1}^N i^n(n!)\\
&=(n!)^n\prod_{i=1}^N i^n\\
&=(n!)^{2n}
\end{aligned}
\]
&\prod_{i=1}^N\prod_{j=1}^N{ij}\\
&=\prod_{i=1}^N i^n(n!)\\
&=(n!)^n\prod_{i=1}^N i^n\\
&=(n!)^{2n}
\end{aligned}
\]
再算分母中的底数：
\[\begin{aligned}
&\prod_{i=1}^N\prod_{j=1}^N\gcd(i,j)\\
&=\prod_{d=1}^N d^{\sum_{i=1}^N\sum_{j=1}^N [\gcd(i,j)=d]}
\end{aligned}
\]
&\prod_{i=1}^N\prod_{j=1}^N\gcd(i,j)\\
&=\prod_{d=1}^N d^{\sum_{i=1}^N\sum_{j=1}^N [\gcd(i,j)=d]}
\end{aligned}
\]
再看指数：
\[\begin{aligned}
&\sum_{i=1}^N\sum_{j=1}^N [\gcd(i,j)=d]\\
&=\sum_{i=1}^N\sum_{j=1}^N [\gcd(\lfloor\frac{i}{d}\rfloor,\lfloor\frac{j}{d}\rfloor)=1]\\
&=\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{d}\rfloor}[\gcd(i,j)=1]
\end{aligned}
\]
&\sum_{i=1}^N\sum_{j=1}^N [\gcd(i,j)=d]\\
&=\sum_{i=1}^N\sum_{j=1}^N [\gcd(\lfloor\frac{i}{d}\rfloor,\lfloor\frac{j}{d}\rfloor)=1]\\
&=\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{d}\rfloor}[\gcd(i,j)=1]
\end{aligned}
\]
然后你可以用莫比乌斯反演：
\[\begin{aligned}
&\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{d}\rfloor}[\gcd(i,j)=1]\\
&=\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{k}\rfloor}\sum_{k|i,k|j}\mu(d)\\
&=\sum_{k=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{d}\rfloor}[k|i][k|j]\mu(d)\\
&=\sum_{k=1}^{\lfloor\frac{N}{d}\rfloor}\mu(k)\lfloor\frac{N}{kd}\rfloor^2
\end{aligned}
\]
&\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{d}\rfloor}[\gcd(i,j)=1]\\
&=\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{k}\rfloor}\sum_{k|i,k|j}\mu(d)\\
&=\sum_{k=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{d}\rfloor}[k|i][k|j]\mu(d)\\
&=\sum_{k=1}^{\lfloor\frac{N}{d}\rfloor}\mu(k)\lfloor\frac{N}{kd}\rfloor^2
\end{aligned}
\]
但这是何必呢？我们可以使用欧拉函数的性质：
\[\begin{aligned}
&\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{d}\rfloor}[\gcd(i,j)=1]\\
&=2(\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\varphi(i))-1
\end{aligned}
\]
&\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{d}\rfloor}[\gcd(i,j)=1]\\
&=2(\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\varphi(i))-1
\end{aligned}
\]
然后发现指数有一点点大，我们可以使用欧拉定理 \(\gcd(a,M)=1,a^b\equiv a^{b\bmod \varphi(M)}\pmod{M}\) 进行缩小。
时间复杂度 \(O(n\log n)\) 空间复杂度 \(O(n)\)。
代码
#include <bits/stdc++.h>
using namespace std;
const int MOD = 104857601;
inline int M(const long long x,const int mod=MOD){return (x%mod+mod)%mod;}
inline int fact(int x){
long long ret=1;
for(int i=1;i<=x;i++) ret=M(ret*i);
return ret;
}
const int N = 1000005;
int n;
int phi[N];
int pri[N],tot;
void sieve(){
phi[1]=1;
for(int i=2;i<=n;i++){
if(!phi[i]){
pri[++tot]=i;phi[i]=i-1;
}
for(int j=1;j<=tot&&(i*pri[j]<=n);j++){
if(!(i%pri[j])){
phi[i*pri[j]]=phi[i]*pri[j];
break;
}
phi[i*pri[j]]=phi[i]*phi[pri[j]];
}
}
for(int i=2;i<=n;i++) phi[i]=M(phi[i]+phi[i-1],MOD-1);
}
long long fastpow(long long a,int b,const int mod){
if(b==0) return 1;
if(b==1) return M(a,mod);
long long ret=fastpow(a,b>>1,mod);
ret=M(ret*ret,mod);
if(b&1) ret=M(ret*a);
return ret;
}
signed main(){
cin>>n;
sieve();
long long a=1;
for(int d=1;d<=n;d++) a=M(a*fastpow(d,phi[(n/d)]*2-1,MOD));
cout<<(M(fastpow(fact(n),n<<1,MOD)*fastpow(M(a*a),MOD-2,MOD)));
return 0;
}
