ICode9
精准搜索请尝试： 精确搜索
首页
编程语言>
数据库
系统相关
互联网
其他分享
Python
Java
JavaScript
android
PHP
首页 > 其他分享> 文章详细

「题解」Luogu P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题

2022-03-08 14:00:07  阅读：158  来源： 互联网

标签：P5518 lfloor right int 题解 sum MtOI2019 prod left


P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题

Description

  • 多测。

  • 在一开始给定模数 \(p\)。

  • 每次给定 \(3\) 个正整数 \(a, b, c\)，求

    \[\prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \left(\dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right)^{f(type)} \bmod p \]

    其中 \(type\in \{0, 1, 2\}\)：

    \[f(0) = 1 \\ f(1) = ijk \\ f(2) = \gcd(i, j, k) \]

    每次分别回答 \(type = 0/ 1/ 2\) 时的答案。

  • \(1\le a, b, c\le 10^5, 10^7 \le p \le 1.05\times 10^9, p\in \mathbb{P}, T = 70\)。

Solution

以下令 \(S(n) = \sum_{i = 1}^n i\)​。

参考了 peterwuyihong 大佬提供的思路，在 \(\prod\) 较多较复杂时，可以对原式去 \(\ln\)，将其变成 \(\sum\) 的形式，处理完之后再使用 \(\exp\) 转回 \(\prod\)​ 的形式。

\[ans = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \left(\dfrac{ij}{\gcd(i, j)\gcd(i, k)} \right)^{f(type)} \]

发现每一项只与 \(2\) 个 \(\prod\) 有关，于是设：

\[g(a, b, c) = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c i^{f(type)} \\ h(a, b, c) = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \gcd(i, j)^{f(type)} \]

则

\[ans = \dfrac{g(a, b, c) g(b, a, c)}{h(a, b, c) h(a, c, b)} \]

分类讨论。

type = 0

\[\begin{aligned} g(a, b, c) & = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c i^1 \\ & = (a!)^{bc} \end{aligned} \]

预处理阶乘做到单次 \(\Omicron(\log p)\)​​。

\[h(a, b, c) = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \gcd(i, j)^1 \]

不妨设 \(a\le b\)。

取 \(\ln\) 得

\[\begin{aligned} \sum_{i = 1}^a \sum_{j = 1}^b \sum_{k = 1}^c \ln(\gcd(i, j)) & = c \sum_{d = 1}^a \ln(d) \sum_{i = 1}^a \sum_{j = 1}^b [\gcd(i, j) = d] \\ & = c \sum_{d = 1}^a \ln(d) \sum_{k = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \mu(k) \left\lfloor\dfrac{a}{dk}\right\rfloor \left\lfloor\dfrac{b}{dk}\right\rfloor \\ & = c \sum_{T = 1}^a \left\lfloor\dfrac{a}{T}\right\rfloor \left\lfloor\dfrac{b}{T}\right\rfloor \sum_{d\mid T} \ln(d) \mu\left(\dfrac{T}{d}\right) \end{aligned} \]

取 \(\exp\) 得

\[\left[\prod_{T = 1}^a \left(\prod_{d \mid T} d^{\mu\left(\frac{T}{d}\right)} \right)^{\left\lfloor\frac{a}{T}\right\rfloor \left\lfloor\frac{b}{T}\right\rfloor} \right]^c \]

中间那个可以 \(\Omicron(n\ln n\log p)\) 预处理，结合整除分块即可 \(\Omicron(\sqrt{n}\log p)\)。

type = 1

\[\begin{aligned} g(a, b, c) & = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c i^{ijk} \\ & = \left(\prod_{i = 1}^a i^i \right)^{S(b) S(c)} \end{aligned} \]

\(\Omicron(n\log p)\) 预处理，\(\Omicron(\log p)\)​ 回答。

\[h(a, b, c) = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \gcd(i, j)^{ijk} \]

不妨设 \(a\le b\)。

取 \(\ln\) 得

\[\begin{aligned} \sum_{i = 1}^a \sum_{j = 1}^b \sum_{k = 1}^c \ln(\gcd(i, j)) ijk & = \sum_{d = 1}^a \ln(d) \sum_{i = 1}^a i \sum_{j = 1}^b j \sum_{k = 1}^c k [\gcd(i, j) = d] \\ & = S(c) \sum_{d = 1}^a \ln(d) \sum_{i = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} id \sum_{j = 1}^{\left\lfloor\frac{b}{d}\right\rfloor} jd [\gcd(i, j) = 1] \\ & = S(c) \sum_{d = 1}^a \ln(d) \sum_{k = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \mu(k) \sum_{i = 1}^{\left\lfloor\frac{a}{dk}\right\rfloor} idk \sum_{j = 1}^{\left\lfloor\frac{b}{dk}\right\rfloor} jdk \\ & = S(c) \sum_{d = 1}^a \ln(d) d^2 \sum_{k = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \mu(k) k^2 S\left(\left\lfloor\dfrac{a}{dk}\right\rfloor \right) S\left(\left\lfloor\dfrac{b}{dk}\right\rfloor \right) \\ & = S(c) \sum_{T = 1}^a S\left(\left\lfloor\dfrac{a}{T}\right\rfloor \right) S\left(\left\lfloor\dfrac{b}{T}\right\rfloor \right) \sum_{d\mid T} \ln(d) d^2 \mu\left(\dfrac{T}{d}\right) \left(\dfrac{T}{d}\right)^2 \\ & = S(c) \sum_{T = 1}^a T^2 S\left(\left\lfloor\dfrac{a}{T}\right\rfloor\right) S\left(\left\lfloor\dfrac{b}{T}\right\rfloor\right) \sum_{d\mid T} \ln(d) \mu\left(\dfrac{T}{d}\right) \end{aligned} \]

取 \(\exp\) 得

\[\left[\prod_{T = 1}^a \left(\prod_{d\mid T} d^{\mu\left(\frac{T}{d}\right) T^2} \right)^{S\left(\left\lfloor\frac{a}{T}\right\rfloor \right) S\left(\left\lfloor\frac{b}{T}\right\rfloor \right)} \right]^{S(c)} \]

照样是 \(\Omicron(n\ln n\log p)\) 的预处理和 \(\Omicron(\sqrt{n} \log p)\) 的回答。

type = 2

直接计算原式（注意原式不可调换 \(a, b\)）

\[\prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \left(\dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right)^{\gcd(i, j, k)} \]

取 \(\ln\) 得

\[\begin{aligned} \sum_{i = 1}^a \sum_{j = 1}^b \sum_{k = 1}^c \ln\left(\dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right) \gcd(i, j, k) & = \sum_{d = 1}^{\min(a, b, c)} \varphi(d) \sum_{i = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \sum_{j = 1}^{\left\lfloor\frac{b}{d}\right\rfloor} \sum_{k = 1}^{\left\lfloor\frac{c}{d}\right\rfloor} \ln\left(\dfrac{\operatorname{lcm}(id, jd)}{\gcd(id, kd)} \right) \\ & = \sum_{d = 1}^{\min(a, b, c)} \varphi(d) \sum_{i = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \sum_{j = 1}^{\left\lfloor\frac{b}{d}\right\rfloor} \sum_{k = 1}^{\left\lfloor\frac{c}{d}\right\rfloor} \ln\left(\dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right) \end{aligned} \]

取 \(\exp\) 得E

\[\prod_{d = 1}^{\min(a, b, c)} \left(\prod_{i = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \prod_{j = 1}^{\left\lfloor\frac{b}{d}\right\rfloor} \prod_{k = 1}^{\left\lfloor\frac{c}{d}\right\rfloor} \dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right)^{\varphi(d)} \]

发现没？中间那一块就是 \(type = 0\) 时的答案！

于是整除分块套整除分块实现 \(\Omicron(\sqrt{n} \log p \cdot \sqrt{n}\log p) = \Omicron(n\log^2 p)\)。

Code

tips : 注意哪些是指数，要 \(\bmod p - 1\)。

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

namespace IO
{
}
using IO::read;

int p, q;

int add1(int a, int b) {return (a + b) % p;}
int add2(int a, int b) {return (a + b) % q;}
int sub1(int a, int b) {return (a - b + p) % p;}
int sub2(int a, int b) {return (a - b + q) % q;}
int mul1(int a, int b) {return (ll)a * b % p;}
int mul2(int a, int b) {return (ll)a * b % q;}
int qpow(int a, int b) {if (b < 0) b += q; int base = a, ans = 1; while (b) {if (b & 1) ans = mul1(ans, base); base = mul1(base, base); b >>= 1;} return ans;}
int inv(int a) {return qpow(a, p - 2);}

const int MAXN = 1e5 + 5;
typedef int arr[MAXN];
const int N = 1e5;

arr prime, mu, phi, phiS, fac, dyy, dyyP, dyyPP, dyyPPI, dyyPI, indexP, xhj, xhjP, xhjPP, xhjPPI, xhjPI;
bool vis[MAXN];

void pre()
{
	mu[1] = phi[1] = 1;
	for (int i = 2; i <= N; i++)
	{
		if (!vis[i])
		{
			prime[++prime[0]] = i;
			mu[i] = -1;
			phi[i] = i - 1;
		}
		for (int j = 1; j <= prime[0] && i * prime[j] <= N; j++)
		{
			vis[i * prime[j]] = true;
			if (i % prime[j] == 0)
			{
				mu[i * prime[j]] = 0;
				phi[i * prime[j]] = phi[i] * prime[j];
				break;
			}
			mu[i * prime[j]] = mu[i] * mu[prime[j]];
			phi[i * prime[j]] = phi[i] * phi[prime[j]];
		}
	}
	
	for (int i = 1; i <= N; i++)
	{
		phiS[i] = add2(phiS[i - 1], phi[i]);
	}
	
	fac[0] = 1;
	for (int i = 1; i <= N; i++)
	{
		fac[i] = mul1(fac[i - 1], i);
	}
	
	for (int i = 1; i <= N; i++)
	{
		dyy[i] = 1;
	}
	for (int i = 1; i <= N; i++)
	{
		for (int j = 1; i * j <= N; j++)
		{
			dyy[i * j] = mul1(dyy[i * j], qpow(i, mu[j]));
		}
	}
	dyyP[0] = dyyPP[0] = 1;
	for (int i = 1; i <= N; i++)
	{
		dyyP[i] = mul1(dyyP[i - 1], dyy[i]);
		dyyPP[i] = mul1(dyyPP[i - 1], dyyP[i]);
	}
	dyyPPI[N] = inv(dyyPP[N]);
	dyyPI[0] = 1;
	for (int i = N - 1; i >= 0; i--)
	{
		dyyPPI[i] = mul1(dyyPPI[i + 1], dyyP[i + 1]);
		dyyPI[i + 1] = mul1(dyyPPI[i + 1], dyyPP[i]);
	}
	
	indexP[0] = 1;
	for (int i = 1; i <= N; i++)
	{
		indexP[i] = mul1(indexP[i - 1], qpow(i, i));
	}
	
	for (int i = 1; i <= N; i++)
	{
		xhj[i] = 1;
	}
	for (int i = 1; i <= N; i++)
	{
		for (int j = 1; i * j <= N; j++)
		{
			xhj[i * j] = mul1(xhj[i * j], qpow(i, mul2(mu[j], mul2(i * j, i * j))));
		}
	}
	xhjP[0] = xhjPP[0] = 1;
	for (int i = 1; i <= N; i++)
	{
		xhjP[i] = mul1(xhjP[i - 1], xhj[i]);
		xhjPP[i] = mul1(xhjPP[i - 1], xhjP[i]);
	}
	xhjPPI[N] = inv(xhjPP[N]);
	xhjPI[0] = 1;
	for (int i = N - 1; i >= 0; i--)
	{
		xhjPPI[i] = mul1(xhjPPI[i + 1], xhjP[i + 1]);
		xhjPI[i + 1] = mul1(xhjPPI[i + 1], xhjPP[i]);
	}
}

int S(int n)
{
	return (ll)n * (n + 1) / 2 % q;
}

int GetSum_phi(int l, int r)
{
	return sub2(phiS[r], phiS[l - 1]);
}

int GetPro_dyy(int l, int r)
{
	return mul1(dyyP[r], dyyPI[l - 1]);
}

int GetPro_xhj(int l, int r)
{
	return mul1(xhjP[r], xhjPI[l - 1]);
}

struct Type0
{
	int g(int a, int b, int c)
	{
		return qpow(fac[a], mul2(b, c));
	}
	
	int h(int a, int b, int c)
	{
		if (a > b)
		{
			swap(a, b);
		}
		int res = 1;
		for (int l = 1, r; l <= a; l = r + 1)
		{
			int k1 = a / l, k2 = b / l;
			r = min(a / k1, b / k2);
			res = mul1(res, qpow(GetPro_dyy(l, r), mul2(k1, k2)));
		}
		return qpow(res, c);
	}
	
	int GetAns(int a, int b, int c)
	{
		int x = mul1(g(a, b, c), g(b, a, c)), y = mul1(h(a, b, c), h(a, c, b));
		return mul1(x, inv(y));
	}
}type0;

struct Type1
{
	int g(int a, int b, int c)
	{
		return qpow(indexP[a], mul2(S(b), S(c)));
	}
	
	int h(int a, int b, int c)
	{
		if (a > b)
		{
			swap(a, b);
		}
		int res = 1;
		for (int l = 1, r; l <= a; l = r + 1)
		{
			int k1 = a / l, k2 = b / l;
			r = min(a / k1, b / k2);
			res = mul1(res, qpow(GetPro_xhj(l, r), mul2(S(k1), S(k2))));
		}
		return qpow(res, S(c));
	}
	
	int GetAns(int a, int b, int c)
	{
		int x = mul1(g(a, b, c), g(b, a, c)), y = mul1(h(a, b, c), h(a, c, b));
		return mul1(x, inv(y));
	}
}type1;

struct Type2
{
	int GetAns(int a, int b, int c)
	{
		int res = 1, upbound = min(a, min(b, c));
		for (int l = 1, r; l <= upbound; l = r + 1)
		{
			int k1 = a / l, k2 = b / l, k3 = c / l;
			r = min(a / k1, min(b / k2, c / k3));
			res = mul1(res, qpow(type0.GetAns(k1, k2, k3), GetSum_phi(l, r)));
		}
		return res;
	}
}type2;

signed main()
{
	int t = read();
	p = read(), q = p - 1;
	pre();
	while (t--)
	{
		int a = read(), b = read(), c = read();
		printf("%d %d %d\n", type0.GetAns(a, b, c), type1.GetAns(a, b, c), type2.GetAns(a, b, c));
	}
	return 0;
}

标签：P5518,lfloor,right,int,题解,sum,MtOI2019,prod,left
来源： https://www.cnblogs.com/mangoworld/p/Solution-Luogu-P5518.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

关于我们 | 联系我们 | 留言反馈

专注分享技术，共同学习，共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有