题解「Luogu5221 Product」
题意
求这个东西:
\[\prod_{i=1}^N\prod_{j=1}^N\frac{{\rm{lcm}}(i,j)}{{\rm{gcd}}(i,j)} \ ({\rm{mod}} \ 104857601)
\]
题解
根据
\[{\rm{lcm}}(i,j)=\frac{i \times j}{{\rm{gcd}}(i,j)}
\]
化简式子:
\[\begin{aligned}
\text{Ans}&=\prod_{i=1}^N\prod_{j=1}^N\frac{i\times j}{{\rm{gcd}}^2(i,j)}\\
&=\frac{\prod_{i=1}^N\prod_{j=1}^Nij}{\prod_{i=1}^N\prod_{j=1}^N{\rm{gcd}}^2(i,j)}
\end{aligned}
\]
先看位于分子的式子:
\[\begin{aligned}
\prod_{i=1}^N\prod_{j=1}^Nij&=\prod_{i=1}^Ni^N (N!)\\
&=(N!)^N(N!)^N\\
&=(N!)^{2N}
\end{aligned}
\]
快速幂求即可。
然后是分母的式子:
先不考虑 \({\rm{gcd}}(i,j)\) 的次数。枚举 \({\rm{gcd}}(i,j)\) ,若有数对 \(i,j\) 满足 \({\rm{gcd}}(i,j)=d\) ,那么 \(d\) 对答案会产生贡献。有:
\[\prod_{i=1}^N\prod_{j=1}^N{\rm{gcd}}(i,j)=\prod_{d=1}^Nd^{f(d)}
\]
\[\begin{aligned}
f(d)&=\sum_{i=1}^N\sum_{j=1}^N[{\rm{gcd}}(i,j)=d]\\
&=\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{k|i,k|j}\mu(k)\\
&=\sum_{k=1}^{\lfloor\frac{N}{d}\rfloor}\mu(k)\lfloor\frac{N}{kd}\rfloor^2
\end{aligned}
\]
枚举 \(d\) ,整除分块求 \(f(d)\) 即可。
\(\text{Code}:\)
#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define maxn 1000005
#define Rint register int
#define INF 0x3f3f3f3f
using namespace std;
typedef long long lxl;
const lxl mod=104857601;
template <typename T>
inline T read()
{
T x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}
while(ch>='0'&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return x*f;
}
int prime[maxn],cnt;
bool flag[maxn];
int mu[maxn];
inline void sieve()
{
mu[1]=1;
for(int i=2;i<maxn;++i)
{
if(!flag[i]) prime[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*prime[j]<maxn;++j)
{
flag[i*prime[j]]=true;
if(!(i%prime[j])) break;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<maxn;++i)
mu[i]+=mu[i-1];
}
inline lxl fmi(lxl a,lxl b)
{
lxl ans=1;
a%=mod;
while(b>0)
{
if(b&1) ans=(ans*a)%mod;
a=(a*a)%mod;
b>>=1;
}
return ans;
}
inline lxl calcu(lxl N)
{
lxl res=0;
for(lxl l=1,r=0;l<=N;l=r+1)
{
r=N/(N/l);
res+=(mu[r]-mu[l-1])*(N/l)*(N/l);
}
return res;
}
inline lxl solve(int N)
{
lxl res=1;
for(lxl i=1;i<=N;++i) res=(res*i)%mod;
res=fmi(res,2*N);
lxl flr=1;
for(lxl d=1;d<=N;++d)
flr=flr*fmi(d,calcu(N/d))%mod;
flr=fmi(flr*flr%mod,mod-2);
return res*flr%mod;
}
int main()
{
// freopen("P5221.in","r",stdin);
sieve();
int N=read<int >();
printf("%lld\n",solve(N));
return 0;
}