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

瞎扯

建议在阅读题解之前欣赏这首由普莉兹姆利巴姐妹带来的的合奏。

Q:你参加省选吗?不是说好了考完 NOIP 就退吗。
A:对啊。
Q:那你学这玩意干啥?
A:对啊,我学这玩意干啥?

写这题的动机?
一是一直很喜欢的曲子,感觉快退役了,圆个梦。
二是写了很多题解了,之前认为最优秀的是 NOI嘉年华的题解,但被叉掉之后不知道该怎么改了,于是删了。其他的都太不精致,都不满意。想在退役之前留下一篇最优秀的题解,于是瞅准了这题。
再有,就是想争口气吧。

最后扯一句,题面里将露娜萨(Lunasa)误写成了 (Lusana)。

简述

原题面:Luogu

给定参数 p,有 T 组数据,每次给定参数 A,B,C,求:

i=1Aj=1Bk=1C(lcm(i,j)gcd(i,k))f(type)

其中 f(type) 的取值如下:

f(type)={1type=0i×j×ktype=1gcd(i,j,k)type=2

1A,B,C105107p1.05×109pPT=70
2.5S,128MB。

分析

原式

先化下原式,原式等于:

i=1Aj=1Bk=1C(i×jgcd(i,j)×gcd(i,k))f(type)

发现每一项仅与两个变量有关,设:

f1(a,b,c)=i=1aj=1bk=1cif(type)f2(a,b,c)=i=1aj=1bk=1cgcd(i,j)f(type)

发现 可以随意交换,则原式等价于:

f1(a,b,c)×f1(b,a,c)f2(a,b,c)×f2(a,c,b)

考虑在 type 取值不同时,如何快速求得 f1f2
一共有 6 个需要推导的式子,不妨就叫它们 16 面了(


type = 0

f1(a,b,c)=i=1aj=1bk=1cif2(a,b,c)=i=1aj=1bk=1cgcd(i,j)

对于 1 面,显然有:

i=1aj=1bk=1ci=(i=1ai)b×c

预处理阶乘 + 快速幂即可,单次计算时间复杂度 O(logn)


再考虑 2 面,套路地枚举 gcd,显然有:

i=1aj=1bk=1cgcd(i,j)=(i=1aj=1bgcd(i,j))c=(d=1d(i=1aj=1b[gcd(i,j)=d]))c

指数是个套路,可以看这里:P3455 [POI2007]ZAP-Queries。于是有:

i=1aj=1b[gcd(i,j)=d]=i=1adj=1bd[gcd(i,j)=1]=i=1adj=1bdkgcd(i,j)μ(k)=k=1μ(k)akdbkd

代回原式,略做处理,则原式等于:

(d=1d(k=1μ(k)akdbkd))c=(d=1(dk=1μ(k))akdbkd)c=d=1(k=1ndd(μ(k)akdbkd))c

「SDOI2017」数字表格 一样,考虑枚举 t=kdd,则原式等于:

t=1n((d|tdμ(td))atbt)c

设:

g0(t)=d|tdμ(td)

线性筛预处理 μ 后,g0(t) 可以用埃氏筛预处理,时间复杂度 O(nlogn)。再代回原式,原式等于:

t=1a(g0(t)atbt)c

预处理 g0(t) 的前缀积和前缀积的逆元,复杂度 O(nlogn)
数论分块 + 快速幂计算即可,单次时间复杂度 O(nlogn)


type = 1

f1(a,b,c)=i=1aj=1bk=1cii×j×kf2(a,b,c)=i=1aj=1bk=1cgcd(i,j)i×j×k

考虑 3 面,把 k 扔到指数位置,有:

i=1aj=1bk=1cii×j×k=i=1aj=1bi(i×j×k=1ck)

再把 j 也扔到指数位置,引入 sum(n)=i=1ni=n(n+1)2,原式等于:

(i=1aii)sum(b)×sum(c)

预处理 ii 的前缀积,复杂度 O(nlogn)
指数可以 O(1) 算出,然后快速幂,单次时间复杂度 O(logn)

根据费马小定理,指数需要对 p1 取模。注意 p1 不是质数,计算 sum 时不能用逆元,但乘不爆 LL,直接算就行。


再考虑 4 面,发现 kgcd 无关,则同样把 k 扔到指数位置,则有:

i=1aj=1bk=1cgcd(i,j)i×j×k=(i=1aj=1bgcd(i,j)i×j)sum(c)

套路地枚举 gcd,原式等于:

(d=1d(i=1aj=1bi×j[gcd(i,j)=d]))sum(c)

大力化简指数,有:

i=1aj=1bi×j[gcd(i,j)=d]=d2i=1adj=1bdi×j[gcd(i,j)=1=d2i=1adij=1bdjt|gcd(i,j)μ(t)=d2i=1adij=1bdjk|gcd(i,j)μ(k)=d2k=1μ(k)i=1adi[k|i]j=1bdj[k|j]=d2k=1k2μ(k)i=1akdij=1bkdj=d2k=1k2μ(k)sum(akd)sum(bkd)

指数化不动了,代回原式,原式等于:

(d=1d(d2k=1k2μ(k)sum(akd)sum(bkd)))sum(c)

同 2 面的情况,先展开一下,再枚举 t=kdd,原式等于:

(d=1(k=1ndd(d2k2μ(k)))(sum(akd)sum(bkd)))sum(c)=t=1((d|td(d2(td)2μ(td)))sum(at)sum(bt))sum(c)=t=1((d|td(t2μ(td)))sum(at)sum(bt))sum(c)

与二面相同,设:

g1(t)=d|td(t2μ(td))

g1(t) 可以用埃氏筛套快速幂筛出,时间复杂度 O(nlog2n)。再代回原式,原式等于:

t=1(g1(t)sum(at)sum(bt))sum(c)

同样预处理 g1(t) 的前缀积及其逆元,时间复杂度 O(nlogn)
整除分块 + 快速幂即可,单次时间复杂度 O(nlogn)

注意指数的取模。


type = 2

f1(a,b,c)=i=1aj=1bk=1cigcd(i,j,k)f2(a,b,c)=i=1aj=1bk=1cgcd(i,j)gcd(i,j,k)

考虑 5 面,手段同上,大力反演化简一波,再调换枚举对象,则有:

i=1aj=1bk=1cigcd(i,j,k)=d=1i=1ai(dj=1bk=1c[gcd(i,j,k)=d])=d=1i=1ai(dj=1bdk=1cd[gcd(id,j,k)=1])=d=1i=1ad(id)(dj=1bdk=1cd[gcd(i,j,k)=1])=d=1i=1ad(id)(dj=1bdk=1cdx|gcd(i,j,k)μ(x))=d=1i=1ad(id)(dx=1μ(x)[x|i]j=1bd[x|j]k=1cd[x|k])=d=1x=1i=1ad(id)(d×μ(x)[x|i]j=1bd[x|j]k=1cd[x|k])=d=1x=1i=1axd(ixd)(d×μ(x)bxdcxd)=t=1d|Ti=1at(it)(d×μ(td)btct)=t=1d|T(tati=1ati)d×μ(td)btct

引入 fac(n)=i=1ni,再根据枚举对象调整一下指数,原式等于:

t=1d|t(tat×fac(at))(d×μ(td)btct)=t=1(d|t(tat×fac(at))d×μ(td))btct=t=1((tat×fac(at))d|td×μ(td))btct

指数中出现了一个经典的狄利克雷卷积的形式,对其进行反演。
(Idμ)(n)=φ(n) 代入原式,原式等于:

t=1(tat×fac(at))φ(t)btct=t=1(tφ(t)atbtct×fac(at)φ(t)btct)

预处理 tφ(t) 的前缀积及逆元,阶乘的前缀积及阶乘逆元,(modp1) 下的 φ(t) 的前缀和(指数
),时间复杂度 O(nlogn)
同样整除分块 + 快速幂即可,单次时间复杂度 O(nlogn)


然后是最掉 sans 的 6 面。有 gcd(i,j,k)=gcd(gcd(i,j),k),考虑先枚举 gcd(i,j),然后套路化式子,则有:

i=1aj=1bk=1cgcd(i,j)gcd(i,j,k)=d=1i=1aj=1bk=1c[gcd(i,j)=d]dgcd(d,k)=d=1(d(i=1aj=1b[gcd(i,j)=d]))k=1cgcd(d,k)

先考虑最外面的指数,这也是个套路,可以参考 一个例子。用 Id=φ1 反演,显然有:

k=1cgcd(d,k)=k=1cx|gcd(d,k)φ(x)=x=1φ(x)[x|d]k=1c[x|k]=x|dφ(x)cx

再考虑里面的指数,发现这式子在 2 面已经推了一遍了,于是直接拿过来用,有:

i=1aj=1b[gcd(i,j)=d]=y=1μ(y)aydbyd

将化简后的两个指数代入原式,原式等于:

d=1(d(y=1μ(y)aydbyd))x|dφ(x)cx=d=1(y=1d(μ(y)aydbyd))x|dφ(x)cx

与 2、4 面同样套路地,考虑枚举 t=ydd,再略作调整,原式等于:

d=1(y=1d(μ(y)aydbyd))x|dφ(x)cx=t=1d|td(μ(td)atbtx|dφ(x)cx)=t=1(d|td(μ(td)x|dφ(x)cx))atbt=t=1(d|tx|dd(μ(td)φ(x)cx))atbt

发现要同时枚举 dx,化不动了。
从题解里学到一个比较神的技巧,考虑把 d 拆成 xdx 分别计算贡献再相乘,即分别计算下两式:

t=1(d|tx|dx(μ(td)φ(x)cx))atbtt=1(d|tx|d(dx)(μ(td)φ(x)cx))atbt


先考虑 x 的情况,首先把枚举 x 调整到最外层。设 lim=max(a,b,c),则原式等于:

x=1t=1lim[x|t](d|t[x|d]x(μ(td)φ(x)cx))atbt=x=1t=1limx(d|tx(μ(txdx)φ(x)cx))atxbtx=x=1t=1limxd|tx(φ(x)cxatxbtxμ(td))

t 挪到指数位置,原式等于:

x=1x(t=1limxd|tφ(x)cxatxbtxμ(td))=x=1x(φ(x)cxt=1limxatxbtxd|tμ(td))

指数中又出现了一个经典的狄利克雷卷积的形式,对其进行反演。
(μ1)(n)=ϵ(n)=[n=1] 代入原式,原式等于:

x=1x(φ(x)cxt=1limxatxbtx[t=1])=x=1x(φ(x)axbxcx)

得到了一个非常优美的式子,而且发现这个式子是 5 面最终答案的一部分。同 5 面的做法,直接整除分块即可。


再考虑 dx 的情况,同上先把枚举 x 放到最外层,并调整一下指数,则原式等于:

x=1t=1lim[x|t](d|t[x|d](dx)(μ(td)φ(x)cx))atbt=x=1t=1limx(d|tx[x|d](dx)(μ(txd)φ(x)cx))atxbtx=x=1(t=1limx(d|tx[x|d](dx)μ(txd))atxbtx)φ(x)cx

考虑枚举 dx,替换原来的 d,注意一下这里的倍数关系。原式等于:

x=1(t=1limx(d|tdμ(td))atxbtx)φ(x)cx

发现最内层的式子 d|tdμ(td),即为二面处理过的 g0(t)。直接代入,原式等于:

x=1(t=1limxg0(t)atxbtx)φ(x)cx

一个小结论,证明可以看 这里

a,b,cZ,abc=abc

则原式等于:

x=1(t=1limxg0(t)axtbxt)φ(x)cx

于是可以先对外层整除分块,再对内层整除分块。

然后就做完了,哈哈哈。

实现

一些实现上的小技巧:

  • 逆元能预处理就预处理。
  • 注意对指数取模,模数为 p1
复制复制
//知识点:莫比乌斯反演
/*
By:Luckyblock
用了比较清晰易懂的变量名,阅读体验应该会比较好。
vsc 的自动补全真是太好啦!
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
using std::min;
using std::max;
#define LL long long
const int Lim = 1e5;
const int kN = 1e5 + 10;
//=============================================================
LL A, B, C, mod, ans;
int T, p_num, p[kN];
bool vis[kN];
LL mu[kN], phi[kN], fac[kN], g[2][kN];
LL sumphi[kN], prodt_phi[kN], prodi_i[kN], prodg[2][kN];
LL inv[kN], inv_fac[kN], inv_prodt_phi[kN], inv_prodg[2][kN];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) {
w = (w << 3) + (w << 1) + (ch ^ '0');
}
return f * w;
}
void Chkmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
x_ %= mod;
y_ %= mod - 1;
LL ret = 1;
for (; y_; y_ >>= 1ll) {
if (y_ & 1) ret = ret * x_ % mod;
x_ = x_ * x_ % mod;
}
return ret;
}
LL Inv(LL x_) {
return QPow(x_, mod - 2);
}
LL Sum(LL n_) {
return ((n_ * (n_ + 1ll)) / 2ll) % (mod - 1);
}
void Euler() {
vis[1] = true, mu[1] = phi[1] = 1; //初值
for (int i = 2; i <= Lim; ++ i) {
if (! vis[i]) {
p[++ p_num] = i;
mu[i] = -1;
phi[i] = i - 1;
}
for (int j = 1; j <= p_num && i * p[j] <= Lim; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
phi[i * p[j]] = phi[i] * p[j];
break;
}
mu[i * p[j]] = -mu[i];
phi[i * p[j]] = phi[i] * (p[j] - 1);
}
}
}
void Prepare() {
Euler();
inv[1] = fac[0] = prodt_phi[0] = prodi_i[0] = 1;
for (int i = 1; i <= Lim; ++ i) {
g[0][i] = g[1][i] = 1;
fac[i] = 1ll * fac[i - 1] * i % mod;
sumphi[i] = (sumphi[i - 1] + phi[i]) % (mod - 1);
prodi_i[i] = prodi_i[i - 1] * QPow(i, i) % mod;
if (i > 1) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
prodt_phi[i] = prodt_phi[i - 1] * QPow(i, phi[i]) % mod;
inv_prodt_phi[i] = Inv(prodt_phi[i]);
}
for (int d = 1; d <= Lim; ++ d) {
for (int j = 1; d * j <= Lim; ++ j) {
int t = d * j;
if (mu[j] == 1) {
g[0][t] = g[0][t] * d % mod;
g[1][t] = g[1][t] * QPow(1ll * d, 1ll * t * t) % mod;
} else if (mu[j] == -1) {
g[0][t] = g[0][t] * inv[d] % mod;
g[1][t] = g[1][t] * Inv(QPow(1ll * d, 1ll * t * t)) % mod;
}
}
}
inv_prodg[0][0] = prodg[0][0] = 1;
inv_prodg[1][0] = prodg[1][0] = 1;
inv_prodt_phi[0] = 1;
for (int i = 1; i <= Lim; ++ i) {
for (int j = 0; j <= 1; ++ j) {
prodg[j][i] = prodg[j][i - 1] * g[j][i] % mod;
inv_prodg[j][i] = Inv(prodg[j][i]);
}
}
}
LL f1(LL a_, LL b_, LL c_, int type) {
if (! type) return QPow(fac[a_], b_ * c_);
if (type == 1) return QPow(prodi_i[a_], Sum(b_) * Sum(c_));
LL ret = 1, lim = min(min(a_, b_), c_);
for (LL l = 1, r = 1; l <= lim; l = r + 1) {
r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
ret = ret * QPow(fac[a_ / l], (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
}
return ret;
}
LL f2_2(LL a_, LL b_) {
LL ret = 1;
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
ret = ret * QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l)) % mod;
}
return ret;
}
LL f2(LL a_, LL b_, LL c_, int type) {
LL ret = 1;
if (! type) {
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
LL val = QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l));
ret = (ret * QPow(val, c_)) % mod;
}
} else if (type == 1) {
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
LL val = QPow(prodg[1][r] * inv_prodg[1][l - 1], Sum(a_ / l) * Sum(b_ / l));
ret = (ret * QPow(val, Sum(c_))) % mod;
}
} else {
LL lim = min(min(a_, b_), c_);
for (LL l = 1, r = 1; l <= lim; l = r + 1) {
r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
ret = ret * QPow(f2_2(a_ / l, b_ / l), (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (c_ / l)) % mod;
ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
}
}
return ret;
}
//=============================================================
int main() {
T = read(), mod = read();
Prepare();
while (T -- ) {
A = read(), B = read(), C = read();
for (int i = 0; i <= 2; ++ i) {
ans = f1(A, B, C, i) * f1(B, A, C, i) % mod;
ans = ans * Inv(f2(A, B, C, i)) % mod * Inv(f2(A, C, B, i)) % mod;
printf("%lld ", ans);
}
printf("\n");
}
return 0;
}
posted @   Luckyblock  阅读(540)  评论(12编辑  收藏  举报
编辑推荐:
· 开发者必知的日志记录最佳实践
· SQL Server 2025 AI相关能力初探
· Linux系列:如何用 C#调用 C方法造成内存泄露
· AI与.NET技术实操系列(二):开始使用ML.NET
· 记一次.NET内存居高不下排查解决与启示
阅读排行:
· Manus重磅发布:全球首款通用AI代理技术深度解析与实战指南
· 被坑几百块钱后,我竟然真的恢复了删除的微信聊天记录!
· 没有Manus邀请码?试试免邀请码的MGX或者开源的OpenManus吧
· 园子的第一款AI主题卫衣上架——"HELLO! HOW CAN I ASSIST YOU TODAY
· 【自荐】一款简洁、开源的在线白板工具 Drawnix
点击右上角即可分享
微信分享提示