首先，一个排列对应的权值为其中所有循环长度的最小公倍数（(operatorname{lcm})）。
那么对每个质因数分别考虑：如某个 (p^x le N)，我们希望统计 (p^x) 对答案贡献了多少个幂次。
例如当 (N = 5) 时样例解释给出了 ([1, 25, 20, 30, 24, 20])。
所以 (2^1) 贡献了 (25 + 20 = 45) 次，(2^2) 贡献了 (30) 次，(3^1) 贡献了 (20 + 20 = 40) 次，(5^1) 贡献了 (24) 次。
答案就是 (2^{45} cdot 4^{30} cdot 3^{40} cdot 5^{24})，容易验证是正确的。
我们发现这样还是有点怪，因为统计的是“恰好”，例如 (2^2) 是无法贡献给 (2^1) 的。
所以我们这样考虑：如果 (operatorname{lcm}) 质因数分解中的一项为 (p^x)，则会给所有 (p^y)（(1 le y le x)）都贡献一次。
也就是现在统计的是 (operatorname{lcm}) 中的 (p) 的幂次超过（而不是恰好等于）(x) 的排列的个数。
那么现在 (2^1) 贡献了 (25 + 30 + 20 = 75) 次，而 (2^2) 贡献了 (30) 次。
但是现在算贡献，因为改变了统计的东西，所以是 (2^{75} cdot 2^{30} cdot 3^{40} cdot 5^{24})（注意是 (2^{30}) 而不是原先的 (4^{30})）。
综上，我们现在计算的目标是
[egin{aligned}
mathbf{Ans} &= prod_{A in S_{[N]}} operatorname{step}(A) \
&= prod_{ egin{subarray}{c} p in mathbb{P} \ p le N end{subarray} } exp_p sum_{p^x le N} sum_{A in S_{[N]}} [p^x mid operatorname{step}(A)]
end{aligned}
]
也就是对于每个 (p^x le n) 计算循环长度 (operatorname{lcm}) 是 (p^x) 的倍数的排列个数，然后计算 (p) 的那么多次方再相乘。
循环长度 (operatorname{lcm}) 要是 (p^x) 的倍数，则其中至少有一个循环的长度是 (p^x) 的倍数。
容斥，转化为所有排列的个数减去每个循环长度都不是 (p^x) 的倍数的排列个数。
所有排列的个数即为 (N!)，而每个循环长度都不是 (p^x) 的倍数，就是强制不选长度为 (p^x) 的倍数的循环的有序组合。
如果你没看懂前一句话的后半部分，这时就是生成函数出场的时候了。为了方便表示，令 (k = p^x)。
我们知道一个排列就是若干个循环的有序组合（循环也就是圆排列）。所以应该有此式：
[exp mathbf{EGF}( extrm{cyclic permutation}) = mathbf{EGF}( extrm{permutation})
]
事实上这确实是成立的，因为
[egin{aligned}
mathbf{EGF}( extrm{permutation}) &= sum_{i = 0}^{infty} frac{i!}{i!} z^i \
&= frac{1}{1 - z} \
mathbf{EGF}( extrm{cyclic permutation}) &= sum_{i = 1}^{infty} frac{(i - 1)!}{i!} z^i \
&= sum_{i = 1}^{infty} frac{z^i}{i} \
&= ln frac{1}{1 - z}
end{aligned}
]
那么每个循环长度都不是 (k = p^x) 的倍数的排列个数，应该是去除长度为 (k) 的倍数的圆排列的有序组合。
[egin{aligned}
mathbf{EGF}( extrm{required sequence}_k) &= exp left[ sum_{i = 1}^{infty} frac{z^i}{i} - sum_{i = 1}^{infty} frac{z^{i k}}{i k}
ight] \
&= frac{1}{1 - z} div exp frac{1}{k} ln frac{1}{1 - z^k} \
&= frac{1}{1 - z} div {left( frac{1}{1 - z^k}
ight)}^{1 / k} \
&= frac{{(1 - z^k)}^{1 / k}}{1 - z}
end{aligned}
]
而我们要求的就是 ( extrm{required sequence}_k) 的第 (N) 项系数，也就是 (displaystyle N! [z^N] frac{{(1 - z^k)}^{1 / k}}{1 - z})。
注意到 ({(1 - z^k)}^{1 / k}) 只有在 (k) 的倍数次项有值，而 (displaystyle frac{1}{1 - z}) 实际上是一个前缀和操作，所以做如下变换：
令 (l = lfloor N / k
floor)，要求的即是（这一步由 EntropyIncreaser 提供）
[egin{aligned}
N! [z^N] frac{{(1 - z^k)}^{1 / k}}{1 - z} &= N! [z^{l k}] frac{{(1 - z^k)}^{1 / k}}{1 - z} \
&= N! [z^{l k}] frac{{(1 - z^k)}^{1 / k}}{1 - z^k} \
&= N! [z^{l k}] {(1 - z^k)}^{1 / k - 1} \
&= N! {(-1)}^l inom{1 / k - 1}{l} \
&= N! inom{l - 1 / k}{l}
end{aligned}
]
那么最终结果也就是 (displaystyle N! - N! inom{l - 1 / k}{l})。
然而组合数的上指标不是整数（因为位置在指数上，要对 (M - 1) 取模，必须避免除法）必须进一步化简：
[egin{align*}
N! inom{l - 1 / k}{l} &= frac{N!}{l!} prod_{i = 1}^{l} (i - 1 / k) \
&= frac{N!}{k^l l!} prod_{i = 1}^l (i k - 1) ag{2} \
&= frac{1 imes 2 imes cdots imes N}{k imes 2 k imes cdots imes l k} prod_{i = 1}^l (i k - 1) \
&= N^{underline{N - lk}} imes prod_{i = 1}^{l} (i k - 1) {(i k - 1)}^{underline{k - 1}}
end{align*}
]
至此，就化为了只有乘法而没有除法的式子了。而且其中涉及的数的范围都在 ([1, N]) 内。
可以发现每个要求的值都是下降幂形式，通过一些方式预处理下降幂即可。
例如如果暴力预处理所有可能的下降幂，或者用暴力预处理组合数乘阶乘的方式，复杂度为 (mathcal O (N^2))。
如果使用猫树预处理，复杂度为 (mathcal O (N log N))。
如果对 (M - 1) 分解质因数并预处理阶乘在每个质因数上的幂次，也可以做到单次询问 (mathcal O (omega(M - 1)))。
（此时直接对 ((2)) 式进行计算即可，对于每个 (k) 需要花费 (mathcal O (l + omega(M - 1))) 的时间）
总之不论使用何种预处理方式，最后答案的表达式就是：
[egin{aligned}mathbf{Ans} &= prod_{ egin{subarray}{c} p in mathbb{P} \ p le N end{subarray} } exp_p sum_{p^x le N} sum_{A in S_{[N]}} [p^x mid operatorname{step}(A)] \&= prod_{ egin{subarray}{c} p in mathbb{P} \ p le N end{subarray} } exp_p sum_{ egin{subarray}{c} k = p^x le N \ l = lfloor N / k
floor end{subarray} } left[ N! - frac{N!}{k^l l!} prod_{i = 1}^l (i k - 1)
ight] \&= prod_{ egin{subarray}{c} p in mathbb{P} \ p le N end{subarray} } exp_p sum_{ egin{subarray}{c} k = p^x le N \ l = lfloor N / k
floor end{subarray} } left[ N! - N^{underline{N - lk}} imes prod_{i = 1}^{l} (i k - 1) {(i k - 1)}^{underline{k - 1}}
ight]end{aligned}
]
然后使用你喜欢的预处理方式计算下降幂，最后再用快速幂计算答案即可。
预处理复杂度在 (mathcal O (N log log N)) 至 (mathcal O (N^2)) 不等，计算复杂度为 (displaystyle mathcal O left( N log log N + frac{N log M}{log N}
ight))。