【题解】Luogu-P5221 Product

P5221 Product

Description

• 给定整数 $n$，请求出
  $$\left[\prod _{i=1}^n\prod _{j=1}^n\frac{\mathrm{lcm}(i,j)}{\gcd (i,j)}\right]\mathrm{m}\mathrm{o}\mathrm{d}104857601$$

• 对于 $100\%$ 的数据，$1\le n\le 10^6$。

Solution

$$\begin{array}{rl}\prod _{i=1}^n\prod _{j=1}^n\frac{\mathrm{lcm}(i,j)}{\gcd (i,j)}& =\prod _{i=1}^n\prod _{j=1}^n\frac{ij}{\gcd (i,j{)}^2}\\ & =\frac{\prod _{i=1}^n\prod _{j=1}^{n}ij}{{\left[\prod _{i=1}^n\prod _{j=1}^n\gcd (i,j)\right]}^2}\end{array}$$
对于分子
$$\begin{array}{rl}\prod _{i=1}^n\prod _{j=1}^{n}ij& =\prod _{i=1}^n{i}^n\cdot n!\\ & =(n!{)}^n\cdot (n!{)}^n\\ & =(n!{)}^{2n}\end{array}$$
对于分母（忽略 $2$ 次方）
$$\begin{array}{rl}\prod _{i=1}^n\prod _{j=1}^n\gcd (i,j)& =\prod _{d=1}^n\prod _{i=1}^n\prod _{j=1}^{n}d[\gcd (i,j)=d]\\ & =\prod _{d=1}^n\prod _{i=1}^{\lfloor \frac{n}{d}\rfloor }\prod _{j=1}^{\lfloor \frac{n}{d}\rfloor }d[\gcd (i,j)=1]\\ & =\prod _{d=1}^n{d}^{{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{n}{d}\rfloor }[\gcd (i,j)=1]}\end{array}$$
拆指数
$$\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }[\gcd (i,j)=1]=\sum _{k=1}^{\lfloor \frac{n}{d}\rfloor }\mu (k){\lfloor \frac{n}{dk}\rfloor }^2$$

代回去
$$\prod _{d=1}^n{d}^{{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{n}{d}\rfloor }[\gcd (i,j)=1]}=\prod _{d=1}^n{d}^{{\sum }_{k=1}^{\lfloor \frac{n}{d}\rfloor }\mu (k){\lfloor \frac{n}{dk}\rfloor }^2}$$
指数用欧拉定理 $\mathrm{m}\mathrm{o}\mathrm{d}(p-1)$ 转成正数并降至 $\mathrm{O}(p)$ 级别。

暴力枚举 + 整除分块。

时间复杂度为 $\mathrm{O}\left({\sum }_{i=1}^n\sqrt{\frac{n}{i}}\cdot \log p\right)=\mathrm{O}(\log p{\sum }_{i=1}^n\sqrt{\frac{n}{i}})$。

当积分算
$${\int }_1^n\sqrt{\frac{n}{x}}dx=2n-2\sqrt{n}$$

所以整个是 $\mathrm{O}(n\log p)$ 的。

注意本题卡空间，$\mu$ 数组的前缀和直接覆盖到 $\mu$ 数组上即可。

Code

// 18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;

const int MAXN = 1e6 + 5;
const int MOD = 104857601;

int p[MAXN / 10], mu[MAXN];
bool vis[MAXN];

void pre(int n)
{
	mu[1] = 1;
	for (int i = 2; i <= n; i++)
	{
		if (!vis[i])
		{
			p[++p[0]] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= p[0] && i * p[j] <= n; j++)
		{
			vis[i * p[j]] = true;
			if (i % p[j] == 0)
			{
				mu[i * p[j]] = 0;
				break;
			}
			mu[i * p[j]] = mu[i] * mu[p[j]];
		}
	}
	for (int i = 1; i <= n; i++)
	{
		mu[i] = (mu[i - 1] + mu[i] + MOD - 1) % (MOD - 1);
	}
}

int GetSum(int l, int r)
{
	return (mu[r] - mu[l - 1] + MOD - 1) % (MOD - 1);
}

int block(int n)
{
	int res = 0;
	for (int l = 1, r; l <= n; l = r + 1)
	{
		int k = n / l;
		r = n / k;
		res = (res + (ll)GetSum(l, r) * k % (MOD - 1) * k % (MOD - 1)) % (MOD - 1);
	}
	return res;
}

int qpow(int a, int b)
{
	int base = a, ans = 1;
	while (b)
	{
		if (b & 1)
		{
			ans = (ll)ans * base % MOD;
		}
		base = (ll)base * base % MOD;
		b >>= 1;
	}
	return ans;
}

int inv(int a)
{
	return qpow(a, MOD - 2);
}

int main()
{
	int n;
	scanf("%d", &n);
	pre(n);
	int fac = 1;
	for (int i = 1; i <= n; i++)
	{
		fac = (ll)fac * i % MOD;
	}
	fac = qpow(fac, 2 * n);
	int ans = 1;
	for (int i = 1; i <= n; i++)
	{
		ans = (ll)ans * qpow(i, block(n / i)) % MOD;
	}
	printf("%d\n", (ll)fac * inv((ll)ans * ans % MOD) % MOD);
	return 0;
}