莫比乌斯反演3

题目戳这里（洛谷P5221）

原式

$$ans=\prod _{i=1}^n\prod _{j=1}^n\frac{lcm(i,j)}{\gcd (i,j)}(mod104857601)$$

解法：

由于 $lcm(i,j)=\frac{i\ast j}{gcd(i,j)}$，原式可化为：
$$ans=\prod _{i=1}^n\prod _{j=1}^n\frac{i\ast j}{{\gcd }^2(i,j)}(mod104857601)$$

$$ans=(\prod _{i=1}^n\prod _{j=1}^{n}i\ast j)\ast (\prod _{i=1}^n\prod _{j=1}^n\frac{1}{{\gcd }^2(i,j)})(mod104857601)$$

前者容易看出：
$$pre=\prod _{i=1}^n\prod _{j=1}^{n}i\ast j=(n!{)}^{2n}$$

后者由于最后结果取模，得到答案 $last=inv(T)\ast inv(T)$ ，其中：
$$last=\prod _{i=1}^n\prod _{j=1}^n\frac{1}{{\gcd }^2(i,j)}$$
$$T=\prod _{i=1}^n\prod _{j=1}^{n}gcd(i,j)$$

对 $T$进行化简，设 $gcd(i,j)=d$ ，对 $d$ 进行枚举得：

$$\prod _{d=1}^{min(n,m)}{d}^{\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)==d]}$$

由于：
$①\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)==x]=\sum _{d=1}^{n/x}\mu (d)\lfloor \frac{n}{dx}{\rfloor }^2$

$②\phi (n)=\sum _{d|n}\mu (d)\ast \lfloor \frac{n}{d}\rfloor$

$得：\sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j)=\sum _{d=1}^n\sum _{d^{\prime }=1}^n\phi (d^{\prime })=2\sum _{d=1}^{n/d}\phi (d)-1$

于是有：已确保 $(n<m)$
$$T=\prod _{d=1}^n{d}^{2\sum _{d^{\prime }=1}^{n/d}\phi (d^{\prime })-1}$$

将上述式子结合得：

$$ans=pre\ast last=pre\ast in{v}^2(T)$$
$$ans=(n!{)}^{2n}\ast in{v}^2(\prod _{d=1}^n{d}^{2\sum _{d^{\prime }=1}^{n/d}\phi (d^{\prime })-1})$$

$$=(n!{)}^{2n}\ast inv(\prod _{d=1}^n{d}^{2\ast (2\sum _{d^{\prime }=1}^{n/d}\phi (d^{\prime })-1)})$$
$因为2\ast (2\sum _{d^{\prime }=1}^{n/d}\phi (d^{\prime })-1)=4\sum _{d^{\prime }=1}^{n/d}\phi (d^{\prime })-2$
所以可以对 $\sum _{d^{\prime }=1}^{n/d}\phi (d^{\prime })$记录前缀和.

代码：

//Siberian Squirrel
//C'est le vie
#include<bits/stdc++.h>
using namespace std;
//#define ACM_LOCAL
typedef long long ll;
#define printd(x) printf("%d\n", x)
#define printlld(x) printf("%lld\n", x)
#define printlf(x) printf("%.3f\n", x)
#define rep(i, l, r) for(int i = l, i##end = r; i <= i##end; ++i)
#define per(i, l, r) for(int i = l, i##end = r; i >= i##end; --i)
const int inf = 0x3f3f3f3f, MAXN = 1000001, MOD = 104857601;
//
int prime[100000], phi[MAXN], cnt = 0;
bool vis[MAXN];
//
int quick_pow(int ans, int p, int res = 1) {
	for (; p; p >>= 1, ans = 1ll * ans * ans % MOD) {
		if (p & 1) res = 1ll * res * ans % MOD;
	}
	return res % MOD;
}
int inv(int ans) {
	return quick_pow(ans, MOD - 2);
}
void get_phi(int n) {//Acquisition de phi, fonction
	phi[1] = 1;
	for(int i = 2; i <= n; ++i) {
		if(!vis[i]) prime[++cnt] = i, phi[i] = i - 1;
	 for(int j = 1; j <= cnt && i * prime[j] <= n; ++j) {
			vis[i * prime[j]] = true;
			if(i % prime[j] == 0) {
				phi[i * prime[j]] = phi[i] * prime[j];
				break;
			} phi[i * prime[j]] = phi[i] * (prime[j] - 1);
		}
	}
	for(int i = 2; i <= n; ++i) phi[i] = (phi[i - 1] + phi[i]) % (MOD - 1);
}
int main() {
	int n, res = 1;
	scanf("%d", &n);
	get_phi(n);//Acquisition de phi, sum
	int fac = 1, pre, T = 1;
	for(int i = 2; i <= n; ++i) fac = 1ll * fac * i % MOD;
	//fac = n!;
	pre = quick_pow(fac, 1ll * 2 * n % (MOD - 1));
	for(int d = 1; d <= n; ++d) {
		int sum = 0;
		T = 1ll * T * quick_pow(d ,1ll * 2 * (2 * phi[n / d] - 1) % (MOD - 1)) % MOD;
	}
	int invT = inv(T) % MOD;
	res = 1ll * pre * invT % MOD;
	printlld(res);
    return 0;
}