- 给定(n)，求(prod_{i=1}^nprod_{j=1}^nfrac{ij}{gcd(i,j)^2})。
- 数据组数(le10^3,nle10^6)
一道很水的反演题，然而日常智障地推错了一步，结果卡了半天。不过最后还是完全靠自己推出来的。
由于好久没做过莫比乌斯反演的题目了，这篇博客就写得详细一些吧。
第一阶段
套路地枚举(gcd)，但要注意由于这道题中是(prod)而不是(sum)，([gcd(i,j)=1])要放在指数上：
[prod_{d=1}^nprod_{i=1}^{lfloorfrac nd
floor}prod_{j=1}^{lfloorfrac nd
floor}(ij)^{[gcd(i,j)=1]}
]
根据(sum_{d|n}mu(d)=[n=1])这一莫比乌斯反演基本性质，转化得到：
[prod_{d=1}^nprod_{i=1}^{lfloorfrac nd
floor}prod_{j=1}^{lfloorfrac nd
floor}(ij)^{sum_{p|gcd(i,j)}mu(p)}
]
根据指数运算法则(x^{a+b}=x^acdot x^b)，可以把指数中的(sum)移下来变成(prod)：
[prod_{d=1}^nprod_{i=1}^{lfloorfrac nd
floor}prod_{j=1}^{lfloorfrac nd
floor}prod_{p|gcd(i,j)}(ij)^{mu(p)}
]
调整枚举顺序，把(p)放到(i,j)前面：
[prod_{d=1}^nprod_{p=1}^{lfloorfrac nd
floor}(prod_{i=1}^{lfloorfrac n{dp}
floor}prod_{j=1}^{lfloorfrac n{dp}
floor}(ipcdot jp))^{mu(p)}
]
而(prod_{i=1}^{lfloorfrac n{dp} floor}prod_{j=1}^{lfloorfrac n{dp} floor}(ipcdot jp))这一项其实就等于((lfloorfrac n{dp} floor!cdot p^{lfloorfrac n{dp} floor})^{2lfloorfrac n{dp} floor})，于是得到：
[prod_{d=1}^nprod_{p=1}^{lfloorfrac nd
floor}(lfloorfrac n{dp}
floor!cdot p^{lfloorfrac n{dp}
floor})^{2mu(p)lfloorfrac n{dp}
floor}
]
第二阶段
套路地令(D=dp)，然后改为枚举(D,p)，得到：
[prod_{D=1}^nprod_{p|D}(lfloorfrac nD
floor!cdot p^{lfloorfrac nD
floor})^{2mu(p)lfloorfrac n{D}
floor}
]
由于题目中的(n)是每次给定的，所以我们需要尽可能将含(n)的项单独抠出，因此得到：
[prod_{D=1}^n(prod_{p|D}(lfloorfrac nD
floor!)^{2mu(p)lfloorfrac n{D}
floor})cdot((prod_{p|D} p^{2mu(p)})^{lfloorfrac n{D}
floor^2})
]
先看(prod_{p|D}(lfloorfrac nD floor!)^{2mu(p)lfloorfrac n{D} floor})这一项，我们重新把(p)移回指数，得到：
[(lfloorfrac nD
floor!)^{2lfloorfrac n{D}
floorsum_{p|D}mu(p)}=(lfloorfrac nD
floor!)^{2lfloorfrac nD
floor[D=1]}
]
也就是说，这一项只有在(D=1)时值才不为(1)，且此时值为((n!)^{2n})。
然后考虑令(f(D)=prod_{p|D}p^{2mu(p)})，这个可以(O(nlogn))预处理。
那么原式就被化简成了：
[(n!)^{2n}prod_{D=1}^nf(D)^{lfloorfrac nD
floor^2}
]
此时只要再对于每次询问除法分块一下就可以了。
代码：(O(nlogn+Tsqrt nlogn))
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 1000000
#define X 19260817
using namespace std;
int n,f[N+5],g[N+5],Fac[N+5],Pt,P[N+5],mu[N+5];
I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
I void Sieve()//线性筛预处理μ
{
RI i,j;for(mu[1]=1,i=2;i<=N;++i)
for(!P[i]&&(mu[P[++Pt]=i]=X-2),j=1;j<=Pt&&i*P[j]<=N;++j)
if(P[i*P[j]]=1,i%P[j]) mu[i*P[j]]=X-1-mu[i];else break;
}
int main()
{
RI Tt,i,j,t;for(Sieve(),Fac[0]=i=1;i<=N;++i) Fac[i]=1LL*Fac[i-1]*i%X,f[i]=1;//预处理阶乘，并给f数组赋初值1
for(i=1;i<=N;++i) for(t=1LL*QP(i,2*mu[i]%(X-1)),j=i;j<=N;j+=i) f[j]=1LL*f[j]*t%X;//预处理f
for(f[0]=g[0]=i=1;i<=N;++i) f[i]=1LL*f[i-1]*f[i]%X,g[i]=QP(f[i],X-2);//预处理f的前缀积及其逆元
RI l,r;scanf("%d",&Tt);W(Tt--)//处理询问
{
for(scanf("%d",&n),t=QP(Fac[n],n<<1),l=1;l<=n;l=r+1)//除法分块
r=n/(n/l),t=1LL*t*QP(1LL*f[r]*g[l-1]%X,1LL*(n/l)*(n/l)%(X-1))%X;
printf("%d
",t);
}return 0;
}