神仙数论题（
首先上手推式子，考虑枚举 \(\newcommand{\lcm}{\operatorname{lcm}}\gcd\)，有 \[ \begin{aligned} & \prod\limits_{x_1,x_2,\ldots,x_n \le m} \left(\lcm_{i=1}^n x_i\right)^{\gcd_{i=1}^n x_i} \\ =& \prod\limits_{d=1}^m \left(\prod\limits_{\substack{x_1,x_2,\ldots,x_n \le \left\lfloor \frac md\right\rfloor}} \left(\lcm_{i=1}^n x_i\right)^{[\gcd_{i=1}^n x_i =1]}\right)^d d^{d\sum_{x_1,x_2,\ldots,x_n \le {\scriptsize\left\lfloor\frac md\right\rfloor}}[\gcd_{i=1}^n x_i=1]} \end{aligned} \]
考虑设 \[f(m) = \prod\limits_{\substack{x_1,x_2,\ldots,x_n \le m}} \left(\lcm_{i=1}^n x_i\right)^{[\gcd_{i=1}^n x_i=1]}\]
和 \[g(m) = \sum_{x_1,x_2,\ldots,x_n \le m}[\mathrm{gcd}_{i=1}^n x_i=1]\]
则 \[ \begin{aligned} & \prod\limits_{x_1,x_2,\ldots,x_n \le m} \left(\lcm_{i=1}^n x_i\right)^{\gcd_{i=1}^n x_i} \\ =& \prod\limits_{d=1}^m f^d\left(\left\lfloor\frac md\right\rfloor\right) d^{d \cdot g\left(\left\lfloor\frac md\right\rfloor\right)} \end{aligned} \]
先考虑 \(g\)，设 \[g_d(m) = \sum_{x_1,x_2,\ldots,x_n \le m}[\mathrm{gcd}_{i=1}^n x_i=d]\]
则显然有 \[g_d(m) = g\left(\left\lfloor\frac md\right\rfloor\right)\]
以及 \[\sum\limits_{d=1}^m g_d(m) = m^n\]
于是有 \[g(m) = m^n - \sum\limits_{d=2}^m g\left(\left\lfloor\frac md\right\rfloor\right)\]
故可以使用类似杜教筛的方式在 \(O(n^{3/4})\) 的时间内求出所需的 \(g\) 值。
考虑 \(f\)，设 \[f_d(m) = \prod\limits_{\substack{x_1,x_2,\ldots,x_n \le m}} \left(\lcm_{i=1}^n x_i\right)^{[\gcd_{i=1}^n x_i=d]}\]
则显然有 \[f_d(m) = f\left(\left\lfloor\frac md\right\rfloor\right) d^{g\left(\left\lfloor\frac md\right\rfloor\right)}\]
以及 \[\prod\limits_{d=1}^m f_d(m) = \prod\limits_{x_1,x_2,\ldots,x_n \le m} \lcm_{i=1}^n x_i\]
此时考虑设 \[h(m) = \prod\limits_{x_1,x_2,\ldots,x_n \le m} \lcm_{i=1}^n x_i\]
则有 \[f(m) = h(m)\left(\prod\limits_{d=2}^m f\left(\left\lfloor\frac md\right\rfloor\right)d^{g\left(\left\lfloor\frac md\right\rfloor\right)}\right)^{-1}\]
考虑如何求 \(h\)。
注意 lcm 和 gcd 的本质：对所有数的每个质因子的指数取 max / min 并乘起来。
注意到，此题中同为乘，那么可以考虑对所有质因子同时做 min-max 容斥。
即 \[h(m) = \prod\limits_{x_1,x_2,\ldots,x_n \le m} \prod\limits_{T \subseteq S} \left(\mathrm{gcd}_{i \in T} x_i\right)^{(-1)^{|T|+1}}\ (S = [1,n] \cap \mathbb Z)\]
考虑枚举 gcd 并计算贡献，设 \[s_t(m) = \sum\limits_{x_1,x_2,\ldots,x_n \le m}\sum\limits_{T \subseteq S} [\mathrm{gcd}_{i \in T} x_i = 1](-1)^{|T|+1}\prod\limits_{i \in T} [x_i \le t]\]
则有 \[h(m) = \prod\limits_{d=1}^m d^{s_{\lfloor m/d\rfloor}(m)}\]
设 \[s_{t,d}(m) = \sum\limits_{x_1,x_2,\ldots,x_n \le m}\sum\limits_{T \subseteq S} [\mathrm{gcd}_{i \in T} x_i = d](-1)^{|T|+1}\prod\limits_{i \in T} [x_i \le t]\]
显然有 \[s_{t,d}(m) = s_{\scriptsize\left\lfloor\frac td\right\rfloor}(m)\]
以及 \[\sum\limits_{d=1}^t s_{t,d}(m) = \sum\limits_{k=1}^n \binom nk t^k m^{n-k} (-1)^{k+1}\]
注意到这很像二项式展开的形式，考虑化一化 \[ \begin{aligned} \sum\limits_{k=1}^n \binom nk t^k m^{n-k} (-1)^{k+1} &= \sum\limits_{k=0}^n \binom nk t^k m^{n-k} (-1)^{k+1} - \binom n0 t^0 m^n (-1)^1 \\ &= m^n + \sum\limits_{k=0}^n \binom nk t^k m^{n-k} (-1)^{k+1} \\ &= m^n - \sum\limits_{k=0}^n \binom nk t^k m^{n-k} (-1)^k \\ &= m^n - \sum\limits_{k=0}^n \binom nk (-t)^k m^{n-k} \\ &= m^n - (m - t)^n \end{aligned} \]
故 \[s_t(m) = m^n - (m-t)^n - \sum\limits_{d=2}^t s_{\scriptsize\left\lfloor\frac td\right\rfloor}(m)\]
求 \(f,g,h,s\) 的总复杂度瓶颈在于 \(h\) 和 \(s\) —— 相当于使用了数论分块套数论分块套数论分块，而其他部分带快速幂的 \(\log\)，故复杂度是 \(O(m^{7/8} + m^{3/4} \log n)\) 的。
然而注意到要求得答案，实际上还需要求出 \(\prod\limits_{d=1}^m d^d\)，而其中的 \(m\) 取 \(O(\sqrt n)\) 种值。
考虑先求 \(\prod\limits_{d=1}^m d^{m-d+1}\) 再除 \((m!)^{m+1}\)。
注意到 \(\prod\limits_{d=1}^m d^{m-d+1} = \prod\limits_{d=1}^m \prod\limits_{i=d}^m d = \prod\limits_{i=1}^m \prod\limits_{d=1}^i d = \prod\limits_{i=1}^m i!\)，所以可以预处理阶乘。
代码： 1
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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
using namespace std;
const int N = 1e8;
const int M = 2e5;
const int MX = 450;
const int mod = 998244353;
inline int fpow(int a,int b,int m = mod)
{
int ret = 1;
for(;b;b >>= 1)
(b & 1) && (ret = (long long)ret * a % m),a = (long long)a * a % m;
return ret;
}
int n,m,lim,ans = 1;
int fac[M + 5],ifac[M + 5],prod[M + 5],iprod[M + 5];
int le[MX + 5],ge[MX + 5],tot;
inline int &id(int x)
{
return x <= lim ? le[x] : ge[m / x];
}
int mem_f[2 * MX + 5],mem_g[2 * MX + 5],mem_h[2 * MX + 5],mem_s[2 * MX + 5][2 * MX + 5];
int g(int m)
{
if(~mem_g[id(m)])
return mem_g[id(m)];
int ret = fpow(m,n,mod - 1);
for(register int l = 2,r;l <= m;l = r + 1)
{
r = m / (m / l);
ret = (ret - (long long)(r - l + 1) * g(m / l) % (mod - 1) + mod - 1) % (mod - 1);
}
return mem_g[id(m)] = ret;
}
int s(int t,int m)
{
if(t > m)
return 0;
if(~mem_s[id(t)][id(m)])
return mem_s[id(t)][id(m)];
int ret = (fpow(m,n,mod - 1) - fpow(m - t,n,mod - 1) + mod - 1) % (mod - 1);
for(register int l = 2,r;l <= t;l = r + 1)
{
r = t / (t / l);
ret = (ret - (long long)(r - l + 1) * s(t / l,m) % (mod - 1) + mod - 1) % (mod - 1);
}
return mem_s[id(t)][id(m)] = ret;
}
int h(int m)
{
if(~mem_h[id(m)])
return mem_h[id(m)];
int ret = 1;
for(register int l = 1,r;l <= m;l = r + 1)
{
r = m / (m / l);
ret = (long long)ret * fpow((long long)fac[r] * ifac[l - 1] % mod,s(m / l,m)) % mod;
}
return mem_h[id(m)] = ret;
}
int f(int m)
{
if(~mem_f[id(m)])
return mem_f[id(m)];
int ret = 1;
for(register int l = 2,r;l <= m;l = r + 1)
{
r = m / (m / l);
ret = (long long)ret * fpow(f(m / l),r - l + 1) % mod * fpow((long long)fac[r] * ifac[l - 1] % mod,g(m / l)) % mod;
}
ret = (long long)h(m) * fpow(ret,mod - 2) % mod;
return mem_f[id(m)] = ret;
}
int main()
{
freopen("lg.in","r",stdin),freopen("lg.out","w",stdout);
memset(mem_f,-1,sizeof mem_f),memset(mem_g,-1,sizeof mem_g),memset(mem_h,-1,sizeof mem_h),memset(mem_s,-1,sizeof mem_s);
scanf("%d%d",&n,&m),fac[0] = prod[0] = 1,lim = sqrt(m);
for(register int i = 1;i <= m;++i)
fac[i] = (long long)fac[i - 1] * i % mod,
prod[i] = (long long)prod[i - 1] * fac[i] % mod;
ifac[m] = fpow(fac[m],mod - 2),
iprod[m] = fpow(prod[m],mod - 2);
for(register int i = m;i;--i)
ifac[i - 1] = (long long)ifac[i] * i % mod,
iprod[i - 1] = (long long)iprod[i] * fac[i] % mod;
for(register int l = 1,r;l <= m;l = r + 1)
{
r = m / (m / l);
id(m / l) = ++tot;
}
for(register int l = 1,r;l <= m;l = r + 1)
{
r = m / (m / l);
ans = (long long)ans * fpow(f(m / l),(long long)(l + r) * (r - l + 1) / 2 % (mod - 1)) % mod * fpow((long long)fpow(fac[r],r + 1) * iprod[r] % mod * fpow(ifac[l - 1],l) % mod * prod[l - 1] % mod,g(m / l)) % mod;
}
printf("%d\n",ans);
}