({ m CYJian})最近闲的玩起了(gcd)。。他想到了一个非常简单而有意思的式子：
[prod_{i=1}^Nprod_{j=1}^Nfrac{lcm(i,j)}{gcd(i,j)} (mod 104857601)
]
({ m CYJian})已经算出来这个式子的值了。现在请你帮他验算一下吧。({ m CYJian})只给你(0.2s)的时间哦。
(1 leq N leq 1000000)
把式子写成：
[=prod_{i=1}^Nprod_{j=1}^Nfrac{ij}{gcd^2(i,j)}=frac{prod_{i=1}^Nprod_{j=1}^Nij}{(prod_{i=1}^Nprod_{j=1}^Ngcd(i,j))^2}
]
上面的部分就是 ((N!)^{2N}) ，对于下面 (prod_{i=1}^Nprod_{j=1}^Ngcd(i,j)) 部分，考虑枚举 (gcd) 是多少，设 (g(x,n)) 表示 (sum_{i=1}^nsum_{j=1}^n[gcd(i,j)==x]) ，那么这个式子就变成了：
[prod_{i=1}^Ni^{g(i,N)}
]
考虑怎么求所有的 (g(i,N)) ，我们写出式子：
[egin{aligned}g(x,n)&=sum_{i=1}^nsum_{j=1}^n[gcd(i,j)==x]\&=sum_{i=1}^{frac{n}{x}}sum_{j=1}^{frac{n}{x}}[gcd(i,j)==1]\&=g(1,frac{n}{x})end{aligned}
]
然后考虑求 (g(1,n)) ，就有：
[egin{aligned}
g(1,n)&=sum_{i=1}^{n}sum_{j=1}^{n}sum_{d|gcd(i,j)}mu(d)\&=sum_{d=1}^nsum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{n}{d}}mu(d)\&=sum_{d=1}^nlfloorfrac{n}{d}
floorlfloorfrac{n}{d}
floormu(d)end{aligned}]
这个可以整除分块 (sqrt{n}) 完成。
然后现在我们要求：
[prod_{i=1}^Ni^{g(1,frac{N}{i})}
]
对 (frac{N}{i}) 也整除分块，这样就可以做到 (O(N+sqrt Nlog N)) 的复杂度处理出来了。
Code
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
const int N = 1e6;
const int M = 1e5;
const int p = 104857601;
using namespace std;
int n,mul[N + 5],prime[M + 5],pcnt,ans,s;
bool vis[N + 5];
int mypow(int a,long long x){int s = 1;for (;x;x & 1 ? s = 1ll * a * s % p : 0,a = 1ll * a * a % p,x >>= 1);return s;}
long long calc(int n)
{
long long ans = 0;
for (int l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
ans += 1ll * (n / l) * (n / l) * (mul[r] - mul[l - 1]);
}
return ans;
}
void prework()
{
vis[1] = 1;mul[1] = 1;
for (int i = 2;i <= n;i++)
{
if (!vis[i])
{
mul[i] = -1;
prime[++pcnt] = i;
}
for (int j = 1;j <= pcnt && prime[j] * i <= n;j++)
{
vis[prime[j] * i] = 1;
mul[prime[j] * i] = -mul[i];
if (i % prime[j] == 0)
{
mul[prime[j] * i] = 0;
break;
}
}
}
for (int i = 1;i <= n;i++)
mul[i] += mul[i - 1];
}
int main()
{
scanf("%d",&n);
prework();
ans = 1;
for (int i = 1;i <= n;i++)
ans = 1ll * ans * i % p;
ans = mypow(ans,n);
ans = 1ll * ans * ans % p;
s = 1;
for (int l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
int tmp = 1;
for (int i = l;i <= r;i++)
tmp = 1ll * tmp * i % p;
s = 1ll * s * mypow(tmp,calc(n / l)) % p;
}
s = 1ll * s * s % p;
ans = 1ll * ans * mypow(s,p - 2) % p;
cout<<(ans + p) % p<<endl;
return 0;
}