题目链接
原始题面
题目背景
\({\rm CYJian}\): " 听说 \(gcd\) 和 \(\sum\) 套起来比较好玩？？那我就......"
题目描述
\({\rm CYJian}\) 最近闲的玩起了 \(gcd\). 他想到了一个非常简单而有意思的式子:
\[ \prod_{i=1}^N\prod_{j=1}^N\frac{lcm(i,j)}{gcd(i,j)}\ (\bmod\ 104857601) \]
\({\rm CYJian}\) 已经算出来这个式子的值了。现在请你帮他验算一下吧. \({\rm CYJian}\) 只给你 \(0.2s\) 的时间哦
输入输出格式
输入格式
一行一个正整数 \(N\)
输出格式
一行一个正整数，表示答案模 \(104857601\) 的值
输入输出样例
输入样例 #1
输出样例 #1
说明
样例解释:
| 1 | 1 | 2 | 3 | 4 | 5 |
| 2 | 2 | 1 | 6 | 2 | 10 |
| 3 | 3 | 6 | 1 | 12 | 15 |
| 4 | 4 | 2 | 12 | 1 | 20 |
| 5 | 5 | 10 | 15 | 20 | 1 |
对于 \(30\%\) 的数据: \(1 \leq N \leq 5000\)
对于 \(100\%\) 的数据: \(1 \leq N \leq 1000000\)
解题思路
推式子
\[ \begin{aligned} \prod_{i=1}^n\prod_{j=1}^n \frac{ij}{(i,j)^2} & =\frac{(n!)^{2n}}{\prod_{i=1}^n\prod_{j=1}^n(i,j)^2} \\ & =\frac{(n!)^{2n}}{\prod_{d=1}^nd^{2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[(i,j)=1]}} \\ & =\frac{(n!)^{2n}}{\prod_{d=1}^nd^{2\sum_{e=1}^{\lfloor\frac{n}{d}\rfloor}\mu(e)\lfloor\frac{n}{de}\rfloor^2}} \\ \end{aligned} \]
两层整除分块即可
时间复杂度
\(O(n)\)
代码参考
Show code
Luogu_P5221view raw1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
|
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 2, mod = 104857601;
i64 qpow(i64 a, i64 b) {
i64 res = 1;
for (; b; b >>= 1, a = a * a % mod)
if (b & 1) res = res * a % mod;
return res;
}
i64 inv(i64 a) { return qpow(a, mod - 2); }
int prime[N / 10], cnt;
bool vis[N];
int smu[N];
void seive(int n = N - 2) {
smu[1] = 1;
for (int i = 2; i <= n; ++i) {
if (!vis[i]) smu[prime[++cnt] = i] = -1;
for (int j = 1; j <= cnt && i * prime[j] <= n; ++j) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
smu[i * prime[j]] = -smu[i];
}
}
for (int i = 2; i <= n; ++i) smu[i] += smu[i - 1];
}
int main() {
seive();
int n;
cin >> n;
i64 ans = 1;
for (int i = 2; i <= n; ++i) (ans *= i) %= mod;
ans = qpow(ans, 2 * n);
i64 fl = 1, fr = 1;
int il = 1, ir = 1;
for (int l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
i64 exp = 0;
for (int l2 = 1, r2, n2 = n / l; l2 <= n2; l2 = r2 + 1) {
r2 = n2 / (n2 / l2);
exp += 1ll * (n2 / l2) * (n2 / l2) % (mod - 1) *
(smu[r2] - smu[l2 - 1] + mod - 1) % (mod - 1);
exp %= mod - 1;
}
while (ir <= r) (fr *= ir++) %= mod;
while (il < l) (fl *= il++) %= mod;
ans *= qpow(inv(fr) * fl % mod, exp * 2 % (mod - 1));
ans %= mod;
}
cout << ans << endl;
return 0;
}
|