题解 P5221 【Product】
Nemlit
2019-02-17 19:51:53
## [原文地址](https://www.cnblogs.com/bcoier/p/10392413.html)
~~先吐槽一波：凉心出题人又卡时间又卡空间~~
先来化简一波柿子
$$\prod_{i=1}^{n}\prod_{j=1}^{n}\frac{lcm(i,j)}{gcd(i,j)}$$
$$=\prod_{i=1}^{n}\prod_{j=1}^{n}\frac{i*j}{gcd(i,j)^2}$$
$$=(\prod_{i=1}^{n}\prod_{j=1}^{n}i*j)*(\prod_{i=1}^{n}\prod_{j=1}^{n}gcd(i,j))^{-2}$$
先看前面的那一坨：
$$\prod_{i=1}^{n}\prod_{j=1}^{n}i*j$$
$$=\prod_{i=1}^{n}(i^n*n!)$$（不理解可以把上述式子打开，就可以发现了）
$$=(n!)^{n}*\prod_{i=1}^{n}i^n$$
$$=(n!)^n*(n!)^n$$
$$=(n!)^{2n}$$
然后我们来看后面那一坨（先不看-2次方）：
$$\prod_{i=1}^{n}\prod_{j=1}^{n}gcd(i,j)$$
$$=\prod_{d=1}^{n}\prod_{i=1}^{n}\prod_{j=1}^{n}[gcd(i,j)==d]$$
$$=\prod_{d=1}^{n}d^{\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)==d]}$$
$$=\prod_{d=1}^{n}d^{\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[gcd(i,j)==1]}$$
先只看指数：
$$\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[gcd(i,j)==1]$$
这不就是[仪仗队](https://tbr-blog.blog.luogu.org/ou-la-han-shuo-zong-jie)吗？
所以我们只要求出啦欧拉函数前缀和，就可以用欧拉函数$O(1)$算出指数了~
所以我们把柿子综合一下：
令$$sum[x]=\sum_{i=1}^x\phi(i)$$
原式化为：$$(n!)^{2n}*(\Pi_{d=1}^{n}d^{2*sum[\frac{n}{d}]-1})^{-2}$$
然后我们就可以$O(nlogn)$求出答案了
最后注意一下，因为欧拉函数前缀和会爆int，所以要用longlong，但是这样就会MLE，所以我们要考虑优化
根据欧拉定理，因为模数为质数，所以$\phi(mod)=mod-1$，所以原式我们可以进一步化为：
$$(n!)^{2n}*(\prod_{d=1}^{n}d^{(2*sum[\frac{n}{d}]-1)\%(mod-1)})^{-2}$$
这样就可以不需要longlong了
下面给出代码：
```
#include<bits/stdc++.h>
using namespace std;
#define il inline
#define re register
#define debug printf("Now is Line : %d\n",__LINE__)
#define file(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
#define mod 104857601
il int read()
{
re int x=0,f=1;re char c=getchar();
while(c<'0'||c>'9') {if(c=='-') f=-1;c=getchar();}
while(c>='0'&&c<='9') x=x*10+c-48,c=getchar();
return x*f;
}
#define maxn 1000000+5
int n,cnt,ans1=1,prim[80000],ans2=1,pai[maxn];
bool vis[maxn];
il int qpow(int a,ll b)
{
int r=1;
while(b)
{
if(b&1ll) r=1ll*r*a%mod;
b>>=1ll;
a=1ll*a*a%mod;
}
return r;
}
int main()
{
n=read();
pai[1]=1;
for(re int i=2;i<=n;++i)
{
ans1=1ll*ans1*i%mod;
if(!vis[i]) prim[++cnt]=i,pai[i]=i-1;
for(re int j=1;j<=cnt;++j)
{
if(prim[j]*i>n) break;
vis[prim[j]*i]=1;
if(i%prim[j]==0) {pai[i*prim[j]]=pai[i]*prim[j];break;}
pai[i*prim[j]]=pai[prim[j]]*pai[i];
}
}
for(re int i=1;i<=n;++i) pai[i]=pai[i]*2+pai[i-1]%(mod-1);
ans1=qpow(ans1,2*n);
for(re int i=2;i<=n;++i) ans2=1ll*ans2*qpow(i,pai[n/i]-1)%mod;
printf("%d",(1ll*ans1*qpow(1ll*ans2*ans2%mod,mod-2))%mod);
return 0;
}
```