「学习笔记」莫比乌斯反演
点击查看目录
目录
- 「学习笔记」莫比乌斯反演
- 前置知识
- 整除分块
- 积性函数
- 线性筛任意积性函数
- 莫比乌斯反演
- 莫比乌斯函数
- 莫比乌斯反演公式
- 例题
- [HAOI2011] Problem b
- YY 的 GCD
- [SDOI2014] 数表
- DZY Loves Math
- [SDOI2015] 约数个数和
- [SDOI2017] 数字表格
- 于神之怒加强版
- [国家集训队] Crash的数字表格 / JZPTAB
- [湖北省队互测2014] 一个人的数论
前置知识
整除分块
考虑快速求：
[sum_{i = 1}^{n}leftlfloordfrac{n}{i}rightrfloor
]
]
发现连续的一段 (leftlfloordfrac{n}{i}rightrfloor) 取值是一样的，而且取值最多只有 (2sqrt{n}) 种，考虑从这里入手，把连续的一段统一处理 .
当前段左端点 (l) 可以通过上一段右端点加一得到，如何快速求右端点 (r) ？给出结论是 (r = leftlfloordfrac{n}{leftlfloorfrac{n}{l}rightrfloor}rightrfloor)，证明如下：
对于 (leftlfloorfrac{n}{x}rightrfloor = k)，有 (n = xk + r(0le rle x))，可推导出不等式 (nge xk)，移项得 (xleleftlfloorfrac{n}{k}rightrfloor)，此时 (x) 最大为 (leftlfloorfrac{n}{k}rightrfloor) .
当前块 (k = leftlfloorfrac{n}{l}rightrfloor)，所以右端点（最大值）取 (r = leftlfloordfrac{n}{leftlfloorfrac{n}{l}rightrfloor}rightrfloor) .
(O(1)) 跑到下一个块，最多 (2sqrt{n}) 个块，所以时间复杂度 (O(sqrt{n}))，比较厉害 .
积性函数
给定数论函数 (f(x))，如果对于任意一组互质的整数 (a, b) 存在 (f(ab) = f(a)f(b))，则称 (f(x)) 为「积性函数」 .
特别的，如果对于任何一组（不要求互质）整数 (a, b) 都存在 (f(ab) = f(a)f(b))，则称 (f(x)) 为「完全积性函数」 .
比较常见的积性函数：
- (varphi(n))：欧拉函数 .
- (sigma_{k}(x))：约数函数，公式为 (sigma_{k}(x) = sum_{dmid x}d^k) . 为方便一般把 (sigma_{0}(x)) 简记为 (tau)，把 (sigma_{1}(x)) 简记为 (sigma(x)) .
- (mu(x))：莫比乌斯函数，本文核心内容，放在下文写 .
比较常见的完全积性函数：
- (epsilon(x) = [n = 1]) .
- (operatorname{id}_k(x)=x^k)，(operatorname{id}_{1}(x)) 通常简记作 (operatorname{id}(x)) .
- (1(x) = 1) .
好像狄利克雷卷积会用，这里挂个名 .
线性筛任意积性函数
看名字就感觉很厉害！
所有积性函数都可以被线性筛，如果只求积性函数前缀和还有「杜教筛」这种复杂度低于线性的高级筛法 . 但没学，目前题也不用 .
注意到线性筛中所有数都只会被它的最小质因子筛到，那么所有只有一个质因子的数的函数值都可以基于积性函数「(f(ab) = f(a)f(b))」这一性质进行处理 .
考虑不止一个质因子的数怎么办 . 设当前筛的函数为 (f(x))，最小质因子为 (j)，处理到的数为 (itimes j)，其中 (i = j^ktimes p_1^{c_1}times p_2^{c_2}timescdotstimes p_n^{c_n}) . 把 (itimes j) 分解为 (dfrac{i}{j^k}) 和 (j^{k + 1}) 这两个互质的数即可算出 (f(itimes j) = fleft(dfrac{i}{j^k}right)fleft(j^{k + 1}right)) . 那你存一下对于每个数的 (k)（最小质因数的指数）就行了 .
如果完全积性就随便了 .
一份筛 (sigma(x)) 的实现，更具一般性：
d[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) {
prime.push_back (i);
low[i] = i, d[i] = i + 1;
}
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) {
low[i * j] = low[i] * j;
if (low[i] == i) d[i * j] = d[i] + i * j;
else d[i * j] = d[i / low[i]] * d[low[i] * j];
break;
}
low[i * j] = j, d[i * j] = d[i] * d[j];
}
}
一份筛 (mu(x)) 的实现，由于定义了有完全平方因子的数的函数值所以很简单（不清楚定义先往下看）：
mu[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) prime.push_back (i), mu[i] = -1;
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
}
莫比乌斯反演
莫比乌斯函数
[mu(x) =
begin{cases}
1 &x = 1\
0 &x 含有完全平方因子\
(-1)^k &x 不含有完全平方因子，不同质因数个数为 k
end{cases}
]
begin{cases}
1 &x = 1\
0 &x 含有完全平方因子\
(-1)^k &x 不含有完全平方因子，不同质因数个数为 k
end{cases}
]
性质：
- (sum_{dmid n}mu(d) = [n = 1])：最重要的性质 . 常用于转化 ([gcd(x, y) = 1])，在狄利克雷卷积中也会出现 .
- 积性函数：意味着可以被快速筛出来 .
莫比乌斯反演公式
形式一：
[f(n) = sum_{dmid n}g(d) Longleftrightarrow g(n) = sum_{dmid n}f(d)muleft(frac{n}{d}right)
]
]
证明直接大力代入，同时使用莫比乌斯函数性质：
[begin{aligned}
sum_{dmid n}mu(d)fleft(frac{n}{d}right)
&= sum_{dmid n}mu(d)sum_{kmid frac{n}{d}}g(k)\
&= sum_{kmid n}g(k)sum_{dmid frac{n}{k}}mu(d)\
&= sum_{kmid n}left[frac{n}{k} = 1right]g(k)\
&= g(n)\
end{aligned}
]
sum_{dmid n}mu(d)fleft(frac{n}{d}right)
&= sum_{dmid n}mu(d)sum_{kmid frac{n}{d}}g(k)\
&= sum_{kmid n}g(k)sum_{dmid frac{n}{k}}mu(d)\
&= sum_{kmid n}left[frac{n}{k} = 1right]g(k)\
&= g(n)\
end{aligned}
]
[tag*{$square$}
]
]
形式二：
[f(n) = sum_{nmid d}g(d) Longleftrightarrow g(n) = sum_{nmid d}f(d)muleft(frac{d}{n}right)
]
]
这个形式好像很不常用啊 .
证明：
令 (k = frac{d}{n}) .
[begin{aligned}
sum_{nmid d}f(d)muleft(frac{d}{n}right)
&= sum_{k = 1}^{+infty}mu(k)g(nk)\
&= sum_{k = 1}^{+infty}mu(k)sum_{nkmid t}f(t)\
&= sum_{nmid t}f(t)sum_{kmid frac{t}{n}}mu(k)\
&= sum_{nmid t}left[frac{t}{n} = 1right]f(k)\
&= g(n)\
end{aligned}
]
sum_{nmid d}f(d)muleft(frac{d}{n}right)
&= sum_{k = 1}^{+infty}mu(k)g(nk)\
&= sum_{k = 1}^{+infty}mu(k)sum_{nkmid t}f(t)\
&= sum_{nmid t}f(t)sum_{kmid frac{t}{n}}mu(k)\
&= sum_{nmid t}left[frac{t}{n} = 1right]f(k)\
&= g(n)\
end{aligned}
]
[tag*{$square$}
]
]
（事实上你会发现，这两个式子完全不会也能做题，因为一般来说你都可以用直接推的方式代替）
例题
会把本来想把常见套路单开一个部分的，想了想还是放在例题中比较好 .
都标的十分明显 .
[HAOI2011] Problem b
这么典的么！
可以容斥，所以只考虑上限为 (n, m) 怎么做 . 下文为了方便钦定 (n .
有式子：
[sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = k]
]
]
「套路一」：变换上界：
[sum_{i = 1}^{leftlfloor n/krightrfloor}sum_{j = 1}^{leftlfloor m/krightrfloor}[gcd(i, j) = 1]
]
]
「套路二」：用莫比乌斯函数的性质把 ([gcd(i, j) = 1]) 换掉：
[sum_{i = 1}^{leftlfloor n/krightrfloor}sum_{j = 1}^{leftlfloor m/krightrfloor}sum_{dmidgcd(i, j)}mu(d)
]
]
「套路三」：变换求和顺序，枚举 (d)：
[sum_{d = 1}^{n}mu(d)sum_{i = 1}^{leftlfloor n/krightrfloor}sum_{j = 1}^{leftlfloor m/krightrfloor}[dmidgcd(i, j)]
]
]
化简成：
[sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{dk}rightrfloorleftlfloorfrac{m}{dk}rightrfloor
]
]
复杂度是 (O(n)) 的，考虑使用整除分块，预处理出 (mu) 函数前缀和，即可 (O(sqrt{n})) 通过 .
点击查看代码
const int N = 5e4 + 10, MX = 5e4;
namespace SOLVE {
int a, b, c, d, k, mu[N], smu[N]; ll ans;
bool vis[N]; std::vector prime;
inline int rnt () {
int x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
smu[i] = smu[i - 1] + mu[i];
}
return;
}
inline int F (int a, int b) {
ll sum = 0;
for (int l = 1, r = 0; l * k YY 的 GCD
考虑枚举质数：
[sum_{pinmathbb{P}}sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = p]
]
]
使用套路一：
[sum_{pinmathbb{P}}sum_{i = 1}^{leftlfloor n / prightrfloor}sum_{j = 1}^{leftlfloor m / prightrfloor}[gcd(i, j) = 1]
]
]
使用套路二：
[sum_{pinmathbb{P}}sum_{i = 1}^{leftlfloor n / prightrfloor}sum_{j = 1}^{leftlfloor m / prightrfloor}sum_{dmidgcd(i, j)}mu(d)
]
]
使用套路三：
[sum_{pinmathbb{P}}sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{pd}rightrfloorleftlfloorfrac{m}{pd}rightrfloor
]
]
「套路四」：为了方便整除分块，枚举分母：
令 (T = pd)，则：
[begin{aligned}
&sum_{pinmathbb{P}}sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor\
=&sum_{t = 1}^{n}sum_{dmid T}mu(d)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor\
=&sum_{t = 1}^{n}leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloorsum_{dmid T}mu(d)\
end{aligned}
]
&sum_{pinmathbb{P}}sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor\
=&sum_{t = 1}^{n}sum_{dmid T}mu(d)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor\
=&sum_{t = 1}^{n}leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloorsum_{dmid T}mu(d)\
end{aligned}
]
预处理出每个数的所有因数的 (mu) 函数之和即可 . 这个题可以直接用调和级数预处理搞，比较厉害 .
点击查看代码
const int N = 1e7 + 10, MX = 1e7;
namespace SOLVE {
int n, m, mu[N], mus[N]; ll ans;
bool vis[N]; std::vector prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
}
far (p, prime) _for (i, 1, MX / p) mus[i * p] += mu[i];
_for (i, 1, MX) mus[i] += mus[i - 1];
return;
}
inline void In () {
n = rnt (), m = rnt ();
return;
}
inline void Solve () {
ans = 0;
for (int l = 1, r = 0; l [SDOI2014] 数表
入门题是不是过了，那不写什么套路几了 .
首先我们不考虑 (sigma(k)le a) 的限制写出式子然后推：
[begin{aligned}
sum_{k = 1}^{n}sigma(k)sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = k]
&=sum_{k = 1}^{n}sigma(k)sum_{i = 1}^{leftlfloor n / prightrfloor}sum_{j = 1}^{leftlfloor m / prightrfloor}[gcd(i, j) = 1]\
&=sum_{k = 1}^{n}sigma(k)sum_{i = 1}^{leftlfloor n / prightrfloor}sum_{j = 1}^{leftlfloor m / prightrfloor}sum_{dmidgcd(i, j)}mu(d)\
&=sum_{k = 1}^{n}sigma(k)sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{dk}rightrfloorleftlfloorfrac{m}{dk}rightrfloor\
end{aligned}
]
sum_{k = 1}^{n}sigma(k)sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = k]
&=sum_{k = 1}^{n}sigma(k)sum_{i = 1}^{leftlfloor n / prightrfloor}sum_{j = 1}^{leftlfloor m / prightrfloor}[gcd(i, j) = 1]\
&=sum_{k = 1}^{n}sigma(k)sum_{i = 1}^{leftlfloor n / prightrfloor}sum_{j = 1}^{leftlfloor m / prightrfloor}sum_{dmidgcd(i, j)}mu(d)\
&=sum_{k = 1}^{n}sigma(k)sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{dk}rightrfloorleftlfloorfrac{m}{dk}rightrfloor\
end{aligned}
]
令 (T = dk)：
[begin{aligned}
sum_{k = 1}^{n}sigma(k)sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{dk}rightrfloorleftlfloorfrac{m}{dk}rightrfloor
&=sum_{T = 1}^{n}leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloorsum_{kmid T}sigma(k)muleft(frac{T}{k}right)
end{aligned}
]
sum_{k = 1}^{n}sigma(k)sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{dk}rightrfloorleftlfloorfrac{m}{dk}rightrfloor
&=sum_{T = 1}^{n}leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloorsum_{kmid T}sigma(k)muleft(frac{T}{k}right)
end{aligned}
]
发现后面是个狄利克雷卷积的形式，然后发现卷不了，因为还有个 (sigma(k)le a) 的限制 .
那么考虑对每个询问的 (a) 升序排序，每次询问时加入满足 (sigma(k)le a) 的 (k)，由于查询的是前缀和可以直接用树状数组维护 .
模数比较特殊是 (2^{31})，直接自然溢出可以跑得飞快，实测正常取模 2.08s，自然溢出 1.15s .
点击查看代码
const ll N = 1e5 + 10, MX = 1e5, Q = 2e4 + 10, P = 1ll prime; bool vis[N]; BIT::BIT tr;
class QU {
public:
int n, m, a, id;
friend inline bool operator MX) break;
vis[i * j] = true;
if (!(i % j)) {
mu[i * j] = 0;
low[i * j] = low[i] * j;
if (low[i] == i) si[i * j] = si[i] + i * j;
else si[i * j] = si[i / low[i]] * si[low[i] * j];
break;
}
mu[i * j] = -mu[i];
low[i * j] = j, si[i * j] = si[i] * si[j];
}
od[i] = i;
}
std::sort (od + 1, od + MX + 1, [](int i, int j) { return si[i] DZY Loves Math
[begin{aligned}
sum_{i = 1}^{n}sum_{j = 1}^{m}f(gcd(i, j))
&= sum_{d = 1}^{n}f(d)sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = d]\
&= sum_{d = 1}^{n}f(d)sum_{i = 1}^{leftlfloor n/drightrfloor}sum_{j = 1}^{leftlfloor m/drightrfloor}sum_{xmid gcd(i, j)}mu(x)\
&= sum_{d = 1}^{n}f(d)sum_{x = 1}^{n}mu(x)leftlfloorfrac{n}{dx}rightrfloorleftlfloorfrac{m}{dx}rightrfloor\
end{aligned}
]
sum_{i = 1}^{n}sum_{j = 1}^{m}f(gcd(i, j))
&= sum_{d = 1}^{n}f(d)sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = d]\
&= sum_{d = 1}^{n}f(d)sum_{i = 1}^{leftlfloor n/drightrfloor}sum_{j = 1}^{leftlfloor m/drightrfloor}sum_{xmid gcd(i, j)}mu(x)\
&= sum_{d = 1}^{n}f(d)sum_{x = 1}^{n}mu(x)leftlfloorfrac{n}{dx}rightrfloorleftlfloorfrac{m}{dx}rightrfloor\
end{aligned}
]
令 (T = dx)：
[begin{aligned}
sum_{d = 1}^{n}f(d)sum_{x = 1}^{n}mu(x)leftlfloorfrac{n}{dx}rightrfloorleftlfloorfrac{m}{dx}rightrfloor
&= sum_{d = 1}^{n}f(d)sum_{x = 1}^{n}muleft(frac{T}{d}right)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor\
&= sum_{T = 1}^{n}leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloorsum_{dmid T}f(d)muleft(frac{T}{d}right)\
end{aligned}
]
sum_{d = 1}^{n}f(d)sum_{x = 1}^{n}mu(x)leftlfloorfrac{n}{dx}rightrfloorleftlfloorfrac{m}{dx}rightrfloor
&= sum_{d = 1}^{n}f(d)sum_{x = 1}^{n}muleft(frac{T}{d}right)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor\
&= sum_{T = 1}^{n}leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloorsum_{dmid T}f(d)muleft(frac{T}{d}right)\
end{aligned}
]
调和级数预处理 .
点击查看代码
const ll N = 1e7 + 10, MX = 1e7;
namespace SOLVE {
ll n, m, mu[N], f[N], low[N], cnt[N], sum[N], ans;
bool vis[N]; std::vector prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x MX) break;
vis[i * j] = true;
if (!(i % j)) {
mu[i * j] = 0, low[i * j] = low[i] * j;
f[i * j] = std::max (f[i], f[low[i]] + 1);
break;
}
mu[i * j] = -mu[i], f[i * j] = f[i], low[i * j] = j;
}
}
_for (i, 1, MX) _for (j, 1, MX / i) sum[i * j] += mu[i] * f[j];
_for (i, 1, MX) sum[i] += sum[i - 1];
return;
}
inline void In () {
n = rnt (), m = rnt ();
return;
}
inline void Solve () {
ans = 0;
for (ll l = 1, r = 0; l [SDOI2015] 约数个数和
给出一个人类想不出来的结论：
[d(ij) = sum_{xmid i}sum_{ymid j}[gcd(x, y) = 1]
]
]
简单解释：假设对于 (ij) 的质因数 (p)，其在 (i) 中的指数为 (a)，在 (j) 中的指数为 (b)，则其在 (ij) 中的指数为 (a + b)，考虑如何选出来所有可能的 (a + b + 1) 种指数 . 我们钦定要选的指数 (kle a) 时从 (i) 取，否则把 (i) 取光后从 (j) 里取剩下的 . 为了方便，「(i) 全选后在 (j) 里选 (k – a) 个」变为「(i) 里不选直接在 (j) 里选剩下的 (k – a) 个」，那就能保证选出来的两个数是互质的，就成了上面那个式子 .
呃呃，写的好抽象，感性理解罢 .
然后就直接推式子：
[begin{aligned}
sum_{i = 1}^{n}sum_{j = 1}^{m}d(ij)
&= sum_{i = 1}^{n}sum_{j = 1}^{m}sum_{xmid i}sum_{ymid j}[gcd(x, y) = 1]\
&= sum_{i = 1}^{n}sum_{j = 1}^{m}sum_{xmid i}sum_{ymid j}sum_{kmidgcd(x, y)}mu(k)\
&= sum_{x = 1}^{n}sum_{y = 1}^{m}bigglfloorfrac{n}{x}biggrfloorbigglfloorfrac{m}{y}biggrfloorsum_{kmidgcd(x, y)}mu(k)\
&= sum_{k = 1}^{n}mu(k)sum_{x = 1}^{lfloor n/krfloor}sum_{y = 1}^{lfloor m/krfloor}bigglfloorfrac{n}{kx}biggrfloorbigglfloorfrac{m}{ky}biggrfloor\
&= sum_{k = 1}^{n}mu(k)(sum_{x = 1}^{lfloor n/krfloor}bigglfloorfrac{n}{kx}biggrfloor)(sum_{y = 1}^{lfloor m/krfloor}bigglfloorfrac{m}{ky}biggrfloor)\
end{aligned}
]
sum_{i = 1}^{n}sum_{j = 1}^{m}d(ij)
&= sum_{i = 1}^{n}sum_{j = 1}^{m}sum_{xmid i}sum_{ymid j}[gcd(x, y) = 1]\
&= sum_{i = 1}^{n}sum_{j = 1}^{m}sum_{xmid i}sum_{ymid j}sum_{kmidgcd(x, y)}mu(k)\
&= sum_{x = 1}^{n}sum_{y = 1}^{m}bigglfloorfrac{n}{x}biggrfloorbigglfloorfrac{m}{y}biggrfloorsum_{kmidgcd(x, y)}mu(k)\
&= sum_{k = 1}^{n}mu(k)sum_{x = 1}^{lfloor n/krfloor}sum_{y = 1}^{lfloor m/krfloor}bigglfloorfrac{n}{kx}biggrfloorbigglfloorfrac{m}{ky}biggrfloor\
&= sum_{k = 1}^{n}mu(k)(sum_{x = 1}^{lfloor n/krfloor}bigglfloorfrac{n}{kx}biggrfloor)(sum_{y = 1}^{lfloor m/krfloor}bigglfloorfrac{m}{ky}biggrfloor)\
end{aligned}
]
考虑预处理出所有 (sum_{i = 1}^{n}leftlfloorfrac{n}{i}rightrfloor)，可以用整除分块做到单次询问 (O(sqrt{n})) .
点击查看代码
const ll N = 5e4 + 10, MX = 5e4;
namespace SOLVE {
ll n, m, mu[N], mus[N], sum[N], ans;
bool vis[N]; std::vector prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
mus[i] = mus[i - 1] + mu[i];
}
_for (i, 1, MX) {
for (ll l = 1, r = 0; l [SDOI2017] 数字表格
这个是乘积，交换运算顺序要化成指数.
[begin{aligned}
prod_{i = 1}^{n}prod_{j = 1}^{m}f(gcd(i, j))
&= prod_{k = 1}^{n}f(k)^{sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = k]}\
&= prod_{k = 1}^{n}f(k)^{sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{dk}rightrfloorleftlfloorfrac{m}{dk}rightrfloor}\
end{aligned}
]
prod_{i = 1}^{n}prod_{j = 1}^{m}f(gcd(i, j))
&= prod_{k = 1}^{n}f(k)^{sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = k]}\
&= prod_{k = 1}^{n}f(k)^{sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{dk}rightrfloorleftlfloorfrac{m}{dk}rightrfloor}\
end{aligned}
]
令 (T = dk)：
[begin{aligned}
prod_{k = 1}^{n}f(k)^{sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{dk}rightrfloorleftlfloorfrac{m}{dk}rightrfloor}
&= prod_{T = 1}^{n}prod_{dmid T}fleft(frac{T}{d}right)^{mu(d)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor}\
&= prod_{T = 1}^{n}left(prod_{dmid T}fleft(frac{T}{d}right)^{mu(d)}right)^{leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor}\
end{aligned}
]
prod_{k = 1}^{n}f(k)^{sum_{d = 1}^{n}mu(d)leftlfloorfrac{n}{dk}rightrfloorleftlfloorfrac{m}{dk}rightrfloor}
&= prod_{T = 1}^{n}prod_{dmid T}fleft(frac{T}{d}right)^{mu(d)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor}\
&= prod_{T = 1}^{n}left(prod_{dmid T}fleft(frac{T}{d}right)^{mu(d)}right)^{leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor}\
end{aligned}
]
(prod_{dmid T}fleft(frac{T}{d}right)^{mu(d)}) 是可以调和级数预处理出来的，剩下部分使用整除分块解决 .
点击查看代码
const ll N = 1e6 + 10, MX = 1e6, P = 1e9 + 7;
namespace SOLVE {
ll n, m, f[N], mu[N], pro[N], inv_f[N], inv_pro[N], ans;
bool vis[N]; std::vector prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x >= 1;
}
return ans;
}
inline void Pre () {
f[1] = 1, inv_f[1] = 1, mu[1] = 1;
_for (i, 2, MX) {
f[i] = (f[i - 1] + f[i - 2]) % P;
inv_f[i] = FastPow (f[i], P - 2);
if (!vis[i]) prime.push_back (i), mu[i] = -1;
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
}
pro[0] = 1, inv_pro[0] = 1;
_for (i, 1, MX) pro[i] = 1;
_for (i, 1, MX) {
_for (j, 1, MX / i) {
if (mu[j] > 0) pro[i * j] = pro[i * j] * f[i] % P;
else if (mu[j] 于神之怒加强版
推导部分比较平凡：
[begin{aligned}
sum_{i = 1}^{n}sum_{j = 1}^{m}gcd(i, j)^k
&= sum_{d = 1}^{n}d^ksum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = d]\
&= sum_{d = 1}^{n}d^ksum_{i = 1}^{leftlfloor n/drightrfloor}sum_{j = 1}^{leftlfloor m/drightrfloor}sum_{xmid gcd(i, j)}mu(x)\
&= sum_{d = 1}^{n}d^ksum_{x = 1}^{n}mu(x)leftlfloorfrac{n}{dx}rightrfloorleftlfloorfrac{m}{dx}rightrfloor\
end{aligned}
]
sum_{i = 1}^{n}sum_{j = 1}^{m}gcd(i, j)^k
&= sum_{d = 1}^{n}d^ksum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) = d]\
&= sum_{d = 1}^{n}d^ksum_{i = 1}^{leftlfloor n/drightrfloor}sum_{j = 1}^{leftlfloor m/drightrfloor}sum_{xmid gcd(i, j)}mu(x)\
&= sum_{d = 1}^{n}d^ksum_{x = 1}^{n}mu(x)leftlfloorfrac{n}{dx}rightrfloorleftlfloorfrac{m}{dx}rightrfloor\
end{aligned}
]
令 (T = dx)：
[begin{aligned}
sum_{d = 1}^{n}d^ksum_{x = 1}^{n}mu(x)leftlfloorfrac{n}{dx}rightrfloorleftlfloorfrac{m}{dx}rightrfloor
&= sum_{d = 1}^{n}d^ksum_{x = 1}^{n}muleft(frac{T}{d}right)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor\
&= sum_{T = 1}^{n}leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloorsum_{dmid T}d^kmuleft(frac{T}{d}right)\
end{aligned}
]
sum_{d = 1}^{n}d^ksum_{x = 1}^{n}mu(x)leftlfloorfrac{n}{dx}rightrfloorleftlfloorfrac{m}{dx}rightrfloor
&= sum_{d = 1}^{n}d^ksum_{x = 1}^{n}muleft(frac{T}{d}right)leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloor\
&= sum_{T = 1}^{n}leftlfloorfrac{n}{T}rightrfloorleftlfloorfrac{m}{T}rightrfloorsum_{dmid T}d^kmuleft(frac{T}{d}right)\
end{aligned}
]
令 (g(T) = sum_{dmid T}^{n}d^kmuleft(frac{T}{d}right))，这部分可以埃筛处理，虽然能过但是交一发就会喜提最劣解 .
观察发现是个卷积形式，而且卷的两个函数都是积性函数，说明 (g) 也是个积性函数 . 考虑如何线性筛筛出来 .
如果 (T) 是一个质数则有 (g(T) = T^k – 1)，否则根据积性函数定义把它分解质因数然后推一下：
[begin{aligned}
g(T)
&=prod_{i = 1}^{t}gleft(p_i^{c_i}right)\
&=prod_{i = 1}^{t}((p_i^{c_i – 1})^kmu(p_i) + (p_i^{c_i})^kmu(1))\
&=prod_{i = 1}^{t}p_i^{k(c_i – 1)}(p_i^k – 1)\
end{aligned}
]
g(T)
&=prod_{i = 1}^{t}gleft(p_i^{c_i}right)\
&=prod_{i = 1}^{t}((p_i^{c_i – 1})^kmu(p_i) + (p_i^{c_i})^kmu(1))\
&=prod_{i = 1}^{t}p_i^{k(c_i – 1)}(p_i^k – 1)\
end{aligned}
]
就可以线性筛预处理了 .
点击查看代码
const ll N = 5e6 + 10, MX = 5e6, P = 1e9 + 7;
namespace SOLVE {
ll k, n, m, mu[N], dk[N], sum[N], ans;
bool vis[N]; std::vector prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x >= 1;
}
return ans;
}
inline void Pre () {
mu[1] = 1, dk[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) prime.push_back (i), mu[i] = -1;
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
dk[i] = FastPow (i, k);
}
_for (i, 1, MX) _for (j, 1, MX / i) sum[i * j] += dk[i] * mu[j] % P;
_for (i, 1, MX) sum[i] = ((sum[i - 1] + sum[i]) % P + P) % P;
return;
}
inline void In () {
n = rnt (), m = rnt ();
return;
}
inline void Solve () {
ans = 0;
for (ll l = 1, r = 0; l [国家集训队] Crash的数字表格 / JZPTAB
硬推：
[begin{aligned}
sum_{i = 1}^{n}sum_{j = 1}^{m}operatorname{lcm}(i, j)
&= sum_{i = 1}^{n}sum_{j = 1}^{m}frac{ij}{gcd(i, j)}\
&= sum_{k = 1}^{n}frac{1}{k}sum_{i = 1}^{n}sum_{j = 1}^{m}ij[gcd(i, j) = k]\
&= sum_{k = 1}^{n}frac{1}{k}sum_{i = 1}^{leftlfloor n/krightrfloor}sum_{j = 1}^{leftlfloor m/krightrfloor}ij[gcd(i, j) = 1]\
&= sum_{k = 1}^{n}frac{1}{k}sum_{i = 1}^{leftlfloor n/krightrfloor}sum_{j = 1}^{leftlfloor m/krightrfloor}k^2ij[gcd(i, j) = 1]\
&= sum_{k = 1}^{n}ksum_{i = 1}^{leftlfloor n/krightrfloor}sum_{j = 1}^{leftlfloor m/krightrfloor}ijsum_{dmidgcd(i, j)}mu(d)\
&= sum_{k = 1}^{n}ksum_{d = 1}^{n}mu(d)sum_{i = 1}^{leftlfloor n/dkrightrfloor}sum_{j = 1}^{leftlfloor m/dkrightrfloor}d^2ij\
&= sum_{k = 1}^{n}ksum_{d = 1}^{n}d^2mu(d)left(sum_{i = 1}^{leftlfloor n/dkrightrfloor}iright)left(sum_{j = 1}^{leftlfloor m/dkrightrfloor}jright)\
end{aligned}
]
sum_{i = 1}^{n}sum_{j = 1}^{m}operatorname{lcm}(i, j)
&= sum_{i = 1}^{n}sum_{j = 1}^{m}frac{ij}{gcd(i, j)}\
&= sum_{k = 1}^{n}frac{1}{k}sum_{i = 1}^{n}sum_{j = 1}^{m}ij[gcd(i, j) = k]\
&= sum_{k = 1}^{n}frac{1}{k}sum_{i = 1}^{leftlfloor n/krightrfloor}sum_{j = 1}^{leftlfloor m/krightrfloor}ij[gcd(i, j) = 1]\
&= sum_{k = 1}^{n}frac{1}{k}sum_{i = 1}^{leftlfloor n/krightrfloor}sum_{j = 1}^{leftlfloor m/krightrfloor}k^2ij[gcd(i, j) = 1]\
&= sum_{k = 1}^{n}ksum_{i = 1}^{leftlfloor n/krightrfloor}sum_{j = 1}^{leftlfloor m/krightrfloor}ijsum_{dmidgcd(i, j)}mu(d)\
&= sum_{k = 1}^{n}ksum_{d = 1}^{n}mu(d)sum_{i = 1}^{leftlfloor n/dkrightrfloor}sum_{j = 1}^{leftlfloor m/dkrightrfloor}d^2ij\
&= sum_{k = 1}^{n}ksum_{d = 1}^{n}d^2mu(d)left(sum_{i = 1}^{leftlfloor n/dkrightrfloor}iright)left(sum_{j = 1}^{leftlfloor m/dkrightrfloor}jright)\
end{aligned}
]
令 (T = dk, S(n) = dbinom{n}{2})：
[begin{aligned}
sum_{k = 1}^{n}ksum_{d = 1}^{n}d^2mu(d)left(sum_{i = 1}^{leftlfloor n/dkrightrfloor}iright)left(sum_{j = 1}^{leftlfloor m/dkrightrfloor}jright)
&= sum_{k = 1}^{n}ksum_{d = 1}^{n}d^2mu(d)Sleft(leftlfloorfrac{n}{T}rightrfloorright)Sleft(leftlfloorfrac{m}{T}rightrfloorright)\
&= sum_{T = 1}^{n}Sleft(leftlfloorfrac{n}{T}rightrfloorright)Sleft(leftlfloorfrac{m}{T}rightrfloorright)sum_{dmid T}d^2mu(d)frac{T}{d}\
&= sum_{T = 1}^{n}Sleft(leftlfloorfrac{n}{T}rightrfloorright)Sleft(leftlfloorfrac{m}{T}rightrfloorright)Tsum_{dmid T}dmu(d)\
end{aligned}
]
sum_{k = 1}^{n}ksum_{d = 1}^{n}d^2mu(d)left(sum_{i = 1}^{leftlfloor n/dkrightrfloor}iright)left(sum_{j = 1}^{leftlfloor m/dkrightrfloor}jright)
&= sum_{k = 1}^{n}ksum_{d = 1}^{n}d^2mu(d)Sleft(leftlfloorfrac{n}{T}rightrfloorright)Sleft(leftlfloorfrac{m}{T}rightrfloorright)\
&= sum_{T = 1}^{n}Sleft(leftlfloorfrac{n}{T}rightrfloorright)Sleft(leftlfloorfrac{m}{T}rightrfloorright)sum_{dmid T}d^2mu(d)frac{T}{d}\
&= sum_{T = 1}^{n}Sleft(leftlfloorfrac{n}{T}rightrfloorright)Sleft(leftlfloorfrac{m}{T}rightrfloorright)Tsum_{dmid T}dmu(d)\
end{aligned}
]
考虑线性筛出 (g(T) = sum_{dmid T}dmu(d)) 后整除分块 .
卷的这两个都是积性函数所以 (g(T)) 也是积性函数，对于 (T = ip(pinmathbb{P})
) 进行分类讨论：
- (i = 1(Tinmathbb{P}))：(g(T) = 1mu(1) + Tmu(T) = 1 – T) .
- (p) 为 (i) 的质因子：(mu(p^2) = 0) 所以产生不了任何贡献 .
- (pperp i)：由于是积性函数所以直接算 (g(T) = g(i)g(p)) .
点击查看代码
const int N = 1e7 + 10, MX = 1e7, P = 1e8 + 9;
namespace SOLVE {
int n, m, g[N], sum[N], s[N], ans;
bool vis[N]; std::vector prime;
inline int rnt () {
int x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x MX) break;
vis[i * j] = true;
if (!(i % j)) { g[i * j] = g[i]; break; }
g[i * j] = 1ll * g[i] * g[j] % P;
}
}
_for (i, 1, MX) {
sum[i] = (sum[i - 1] + 1ll * g[i] * i % P) % P;
s[i] = (s[i - 1] + i) % P;
}
return;
}
inline void In () {
n = rnt (), m = rnt ();
if (n > m) std::swap (n, m);
return;
}
inline void Solve () {
ans = 0;
for (int l = 1, r = 0; l [湖北省队互测2014] 一个人的数论
推式子：
[begin{aligned}
sum_{i = 1}^{n}i^k[gcd(i, j) = 1]
&= sum_{i = 1}^{n}i^ksum_{dmid i, dmid n}mu(d)\
&= sum_{dmid n}mu(d)sum_{dmid i}i^k\
&= sum_{dmid n}mu(d)d^ksum_{i = 1}^{lfloor n/drfloor}i^k\
end{aligned}
]
sum_{i = 1}^{n}i^k[gcd(i, j) = 1]
&= sum_{i = 1}^{n}i^ksum_{dmid i, dmid n}mu(d)\
&= sum_{dmid n}mu(d)sum_{dmid i}i^k\
&= sum_{dmid n}mu(d)d^ksum_{i = 1}^{lfloor n/drfloor}i^k\
end{aligned}
]
发现后面自然数幂和 .
喜报：我不会伯努利数！
但是之前扰动法的闲话中写过（原来闲话能有这种用途）：
[S_{k}(n) = sum_{i = 0}^{n}i^k = frac{(n + 1) ^ {k + 1} – sum_{j = 0}^{k – 1}dbinom{k + 1}{j}S_j(n)}{(k + 1)}
]
还是可以看出来它是一个 (k+1) 次多项式的 .
设 (f(x) = sum_{i = 0}^{x}i^k = sum_{i = 0}^{k + 1}a_ix^i) .
那么继续推式子：
[begin{aligned}
sum_{dmid n}mu(d)d^ksum_{i = 1}^{lfloor n/drfloor}i^k
&= sum_{dmid n}mu(d)d^kfleft(frac{n}{d}right)\
&= sum_{dmid n}mu(d)d^ksum_{i = 0}^{k + 1}a_{i}left(frac{n}{d}right)^{i}\
&= sum_{i = 0}^{k + 1}a_{i}n^{i}sum_{dmid n}mu(d)d^{k – i}\
end{aligned}
]
sum_{dmid n}mu(d)d^ksum_{i = 1}^{lfloor n/drfloor}i^k
&= sum_{dmid n}mu(d)d^kfleft(frac{n}{d}right)\
&= sum_{dmid n}mu(d)d^ksum_{i = 0}^{k + 1}a_{i}left(frac{n}{d}right)^{i}\
&= sum_{i = 0}^{k + 1}a_{i}n^{i}sum_{dmid n}mu(d)d^{k – i}\
end{aligned}
]
令 (g(n) = sum_{dmid n}mu(d)d^{k – i})，是个积性函数，所以 (g(n) = prod_{i = 1}^{w}g(p_i^{alpha_i})) . 考虑计算 (g(p_i^{alpha_i}))，(mu) 只有 (d = 1) 和 (d = p_i) 时非零，易得 (g(p_i^{alpha_i}) = mu(1)1^{k – i} + mu(p_i)p_i^{k – i} = 1 – p_i^{k – i})，即：
[sum_{i = 0}^{k + 1}a_{i}n^{i}prod_{i = 1}^{w}(1 – p_i^{k – i})
]
]
唯一的问题是 (a_i) 怎么求，数据范围允许 (O(d^3))，所以直接考虑直接高斯消元 .
点击查看代码
const ll W = 1e3 + 10, N = 110, P = 1e9 + 7;
namespace SOLVE {
ll k, w, n, p[W], alpha[W], g[N][N], a[N], ans;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x >= 1;
}
return ans;
}
inline ll Inv (ll a) { return FastPow (a, P - 2); }
inline void Build () {
_for (i, 0, k + 1) {
_for (j, 0, k + 1) g[i][j] = FastPow (i + 1, j);
_for (j, 1, i + 1) g[i][k + 2] = (g[i][k + 2] + FastPow (j, k)) % P;
}
return;
}
inline void Gauss () {
_for (i, 0, k + 1) {
ll l = i;
_for (j, i + 1, k + 1) if (g[j][i] > g[l][i]) l = j;
std::swap (g[i], g[l]);
ll fm = Inv (g[i][i]);
_for (j, i, k + 2) g[i][j] = g[i][j] * fm % P;
_for (j, 0, k + 1) {
if (i == j) continue;
_for (l, i + 1, k + 2) g[j][l] = (g[j][l] + P - g[j][i] * g[i][l] % P) % P;
}
}
_for (i, 0, k + 1) a[i] = g[i][k + 2];
return;
}
inline void In () {
k = rnt (), w = rnt (), n = 1;
_for (i, 1, w) {
p[i] = rnt (), alpha[i] = rnt ();
n = n * FastPow (p[i], alpha[i]) % P;
}
return;
}
inline void Solve () {
Build (), Gauss ();
ans = 0;
_for (i, 0, k + 1) {
ll prod = 1;
_for (j, 1, w) prod = prod * (1 - FastPow (p[j], k - i + P - 1) + P) % P;
ans = (ans + a[i] * FastPow (n, i) % P * prod % P) % P;
}
return;
}
inline void Out () {
printf ("%lldn", ans);
return;
}
}
机房租用，北京机房托管，大带宽租用，IDC机房服务器主机租用托管-价格及服务咨询 www.e1idc.net
