解法
首先观察一下,可以发现对于a,b,正方形的边长是lcm(a,b)。那么买的地板的数量是
(
a
∗
b
)
g
c
d
(
a
,
b
)
2
\frac{(a*b)}{gcd(a,b)^2}
gcd(a,b)2(a∗b)
那么对于一个n的答案就是:
a
n
s
=
∏
i
=
1
n
∏
j
=
1
n
i
∗
j
g
c
d
(
i
,
j
)
2
ans=\prod_{i=1}^n \prod_{j=1}^n \frac{i*j}{gcd(i,j)^2}
ans=∏i=1n∏j=1ngcd(i,j)2i∗j
因为直接算看起来并不好算,所以考虑分成分子分母两部分分别计算:
分子:
∏
i
=
1
n
∏
j
=
1
n
i
∗
j
=
(
n
!
)
2
n
\prod_{i=1}^n \prod_{j=1}^n i*j=(n!)^{2n}
∏i=1n∏j=1ni∗j=(n!)2n
主要是分母:
∏
i
=
1
n
∏
j
=
1
n
g
c
d
(
i
,
j
)
2
\prod_{i=1}^n \prod_{j=1}^n gcd(i,j)^2
∏i=1n∏j=1ngcd(i,j)2
现在可以计算
n
∗
n
n*n
n∗n的矩阵中每种gcd的个数,然后算出来分母的具体值
直接枚举
g
c
d
(
i
,
j
)
=
d
gcd(i,j)=d
gcd(i,j)=d,然后计算矩阵中有多少个位置的gcd=d,设为cnt(d)
这里需要知道一个矩阵中有多少个位置互质:
∑
i
=
1
n
∑
j
=
1
n
[
g
c
d
(
i
,
j
)
=
1
]
\sum_{i=1}^n\sum_{j=1}^n [gcd(i,j)=1]
∑i=1n∑j=1n[gcd(i,j)=1]
考虑莫比乌斯反演
=
∑
d
=
1
n
∑
i
=
1
n
/
d
∑
j
=
1
n
/
d
μ
(
d
)
=\sum_{d=1}^n\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}\mu(d)
=∑d=1n∑i=1n/d∑j=1n/dμ(d)
=
∑
d
=
1
n
μ
(
d
)
∗
(
n
/
d
)
2
=\sum_{d=1}^n \mu(d)*(n/d)^2
=∑d=1nμ(d)∗(n/d)2
设这个式子的值为
S
(
d
)
S(d)
S(d)
这个式子可以整除分块做到
o
(
n
)
o(\sqrt n)
o(n)
然后回头计算每种gcd的个数:
c
n
t
(
d
)
=
S
(
n
/
d
)
cnt(d)=S(n/d)
cnt(d)=S(n/d)
计算分母:
∏
d
=
1
n
d
c
n
t
(
d
)
∗
2
\prod_{d=1}^n d^{cnt(d)*2}
∏d=1ndcnt(d)∗2
=
∏
d
=
1
n
d
S
(
n
/
d
)
∗
2
=\prod_{d=1}^n d^{S(n/d)*2}
=∏d=1ndS(n/d)∗2
然后又可以对S(n/d)整除分块
细节: c n t ( d ) 或 者 S ( n / d ) cnt(d)或者S(n/d) cnt(d)或者S(n/d)出现在幂的位置,所以对mod-1取模,而不是mod
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int maxn=1e6+5;
const int mod=19260817;
inline int read(){
char c=getchar();int t=0,f=1;
while((!isdigit(c))&&(c!=EOF)){if(c=='-')f=-1;c=getchar();}
while((isdigit(c))&&(c!=EOF)){t=(t<<3)+(t<<1)+(c^48);c=getchar();}
return t*f;
}
int t,n,s[maxn],rec[maxn];
int mu[maxn],vis[maxn],p[maxn],cnt;
inline int am1(int x){
if(x<0)return x+mod-1;
return x>=(mod-1)?x-(mod-1):x;
}
void get(){
mu[1]=1;
for(int i=2;i<=1e6;i++){
if(!vis[i]){
mu[i]=-1;p[++cnt]=i;
}
for(int j=1;j<=cnt&&i*p[j]<=1e6;j++){
mu[i*p[j]]=-mu[i];vis[i*p[j]]=1;
if(i%p[j]==0){
mu[i*p[j]]=0;break;
}
}
}
for(int i=2;i<=1e6;i++)mu[i]=am1(mu[i]+mu[i-1]);
}
int inv[maxn];
inline int ksm(int a,int b){
int ans=1;
while(b){
if(b&1)ans=ans*a%mod;
a=a*a%mod;b>>=1;
}
return ans;
}
inline int mul(int l,int r){
return rec[r]*inv[l-1]%mod*rec[r]%mod*inv[l-1]%mod;
}
inline int sum(int x){
if(s[x])return s[x];
int r,ans=0;
for(int l=1;l<=x;l=r+1){
r=(x/(x/l));
ans=am1(ans+am1(mu[r]-mu[l-1])*(x/l)%(mod-1)*(x/l)%(mod-1));
}
s[x]=ans;
return ans;
}
signed main(){
t=read();
rec[0]=1;inv[0]=1;
for(int i=1;i<=1e6;i++)rec[i]=rec[i-1]*i%mod;
inv[1000000]=ksm(rec[1000000],mod-2);
for(int i=1e6-1;i>=1;i--)inv[i]=inv[i+1]*(i+1)%mod;
get();
while(t--){
n=read();
int r,ans=1;
for(int l=1;l<=n;l=r+1){
r=(n/(n/l));
ans=1ll*ans*(ksm(mul(l,r),(sum(n/l))%(mod-1))%mod)%mod;
//printf("%d %d %d %d\n",l,r,mul(l,r),sum(n/l));
}
printf("%lld\n",ksm(rec[n],(2*n)%(mod-1))*ksm(ans,mod-2)%mod);
}
return 0;
}