Mobius反演与积性函数前缀和演学习笔记 BZOJ 4176 Lucas的数论 SDOI 2015 约数个数和
下文中所有讨论都在数论函数范围内开展.
数论函数指的是定义域为正整数域, 且值域为复数域的函数.
数论意义下的和式处理技巧
-
因子
\[\sum_{d | n} a_d = \sum_{d | n} a_{\frac n d} \] -
双重因子
- 遍历 + 因子
- \(\mu\)函数相关(\(\mu\)函数后面会讲到, 这里先放公式)
莫比乌斯反演
引入
现给定一个函数
根据\(F\)的定义, 我们有
假设我们已知所有\(F\)值, 考虑如何倒推出\(f\)值:
倒推的过程是否有什么通法呢?
实际上是有的, 这里就用到了Mobius反演:
Mobius函数
我们注意到, 这个式子中出现了函数\(\mu\). \(\mu(d)\)为Mobius函数, 其定义如下:
我们尝试列举出\(\mu\)的前几个值
| \(n\) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| \(\mu (n)\) | 1 | -1 | -1 | 0 | -1 | 1 | -1 | 0 | 0 | 1 | -1 | 0 |
Mobius反演的证明
回到那个公式, 考虑怎么定理证明:
首先证明怎么在已知左边等式的情况下得到右边等式:
然后我们还需要反过来证明右边怎么得到左边:
莫比乌斯反演还有一种扩展:
证明:
证明左推右:
右推左:
证毕.
积性函数
定义
设\(f(x)\)是一个数论函数, 若\(\forall a \perp b\)都有\(f(ab) = f(a) f(b)\), 则称\(f(x)\)是积性函数.
设\(f(x)\)是一个数论函数, 若\(\forall a, b\)都有\(f(ab) = f(a)f(b)\), 则称\(f(x)\)为完全积性函数.
常见的积性函数
除数函数: \(\sigma_k(n)\)表示所有约数的\(k\)次幂之和.
Euler函数: \(\phi(n)\)表示不超过\(n\)且与\(n\)互质的正整数个数.
Mobius函数: 就是前文讲述的\(\mu\).
常见的完全积性函数
幂函数: \(Id_k(n) = n^k\). 特别地, 我们令\(Id(n) = 1\).
单位函数:
线性筛
线性筛是计算积性函数的常用方法, 可以\(O(n)\)计算出2~n的最小质因子及其次数, 利用这些信息对前\(n\)项的数值完成递推.
例: 计算前\(n\)个Euler函数
首先, Euler函数的计算公式如下:
感性的理解是这样的: 假设\(m \perp n\), 则\(m\)不含有\(n\)的任意质因子, 因此我们通过\(\frac{p - 1}p\)将含有\(p\)这个质因子的所有数筛掉即可. 又由于质数两两互质, 因此不会出现筛重复的情况.
至于为什么Euler函数是积性函数, 通过公式显而易见.
考虑如何用线性筛求\(\phi(1) \cdots \phi(n)\): 假设我们已经求出了\(\phi(1) \cdots \phi(n - 1)\), 且\(p\)为\(n\)的最小质因子. 我们令\(n' = \frac n p\), 则分以下两种情况讨论:
- \(p | n'\), 则\(n'\)含有\(n\)中所有质因子, 因此\(为的质因子为质因子\phi(n) = n \prod_{p为n的质因子} \frac{p - 1} p = p \times n' \prod_{p为n'质因子}\frac {p - 1} p = p \times \phi(n - 1)\).
- \(p \nmid n'\), 则由于\(\phi\)为积性函数, 因而\(\phi(n) = \phi(p) \times \phi(n') = (p - 1) \phi(n')\).
因而我们只要通过线性筛得到每个数的最小质因子即可.
下面的代码可以线性求所有常见积性函数
/*
minDivisor 表示最小质因子^其次方数
minDivisorDegree 最小质因子次数
mu
phi
sgm[0] sigma_0函数
sgm[1] sigma_1函数
*/
#include <cstdio>
#include <cstring>
const int N = (int)1e5;
typedef int arr[N + 7];
int tot;
arr isNotPrime, prm, minDivisor, minDivisorDegree, phi, mu, sgm[2];
inline void initialize()
{
memset(isNotPrime, 0, sizeof(isNotPrime)); tot = 0;
mu[1] = 1; minDivisor[1] = minDivisorDegree[1] = -1; phi[1] = 0; sgm[0][1] = sgm[1][1] = 1;
for (int i = 2; i <= N; ++ i)
{
if (! isNotPrime[i])
{
prm[tot ++] = i;
minDivisor[i] = i;
minDivisorDegree[i] = 1;
mu[i] = -1;
phi[i] = i - 1;
sgm[0][i] = 2; sgm[1][i] = i + 1;
}
for (int j = 0; j < tot && i * prm[j] <= N; ++ j)
{
int x = i * prm[j];
isNotPrime[x] = 1;
if (i % prm[j])
{
minDivisor[x] = prm[j];
minDivisorDegree[x] = 1;
mu[x] = - mu[i];
phi[x] = (prm[j] - 1) * phi[i];
sgm[0][x] = sgm[0][i] << 1;
sgm[1][x] = sgm[1][i] * (prm[j] + 1);
}
else
{
minDivisor[x] = prm[j] * minDivisor[i];
minDivisorDegree[x] = minDivisorDegree[i] + 1;
mu[x] = 0;
phi[x] = prm[j] * phi[i];
sgm[0][x] = sgm[0][i / minDivisor[i]] * (minDivisorDegree[x] + 1);
sgm[1][x] = sgm[1][i / minDivisor[i]] * minDivisor[x] + sgm[1][i];
}
}
}
}
int main()
{
initialize();
}
Dirichlet卷积
我们定义两个数论函数\(f(x)\), \(g(x)\)的Dirichlet卷积为
Dirichlet卷积的运算性质:
- 交换律: \(f \times g = g \times f\).
- 结合律: \(f \times g \times h = f \times (g \times h)\).
- 分配律: \(f \times (g + h) = f \times g + f \times h\).
- 单位元: \(f \times \epsilon = f\).
- 封闭性: 假如\(f\), \(g\)均为积性函数, 则\(f \times g\)也是积性函数.
主要形式
设\(f(n)\)是一个数论函数, 我们需要计算\(F(n) = \sum_{i = 1}^n f(i)\).
我们取另一个函数\(g\), 则有
因此
假如我们可以\(O(1)\)求出\(\sum_{k = 1}^n (f \times g) (k)\)以及\(\sum_{k =1}^n g(k)\), 则我们的主要时间复杂度在于递归计算\(F(\lfloor \frac n k \rfloor)\). 递归计算每一个\(F(\lfloor \frac n k \rfloor)\), 并进行记忆化.
对于一个\(n\), 考虑\(\lfloor \frac n k \rfloor\)有多少个值:
- \(1 \le k \le \sqrt{n}\), 最多\(\sqrt{n}\)个值.
- \(\sqrt n < k \le n\), 则\(\lfloor \frac n k \rfloor \le \sqrt n\), 因而最多\(\sqrt n\)个值.
因此对于每一个\(n\), 在计算\(F(n)\)时我们要作\(O(\sqrt n)\)次递归运算.
考虑
因此对于某个\(n\), 在计算\(F(n)\)时, 衍生出的所有需要计算的不同\(F(n')\)个数仍然为\(\lfloor \frac n k \rfloor\)的不同个数个, 即\(O(\sqrt n)\)个. 记忆化递归即可.
时间复杂度的计算有一点麻烦:
我们把\(\lfloor \frac n k \rfloor\)分为两部分来看:
- \(k \le \sqrt{n}\), 我们假设每个\(\lfloor \frac n k \rfloor\)都是不同的, 则这一部分的计算次数可以近似地看作是
因此时间复杂度不超过\(O(n^\frac 3 4)\)
-
\(k > \sqrt{n}\), 则有\(\lfloor \frac n k \rfloor \le \sqrt{n}\). 我们假设这\(\sqrt{n}\)个值都出现一次, 则计算次数可以近似地看作是
-
\[\int_0^{n^\frac 1 2} \sqrt{x} dx = \frac 2 3 n^{\frac 3 4} \]
因此时间复杂度不超过\(O(n^\frac 3 4)\)
因此总时间复杂度为\(O(n^\frac 3 4)\). 当然, 注意到上面的两个"假设", 这注定了计算出的时间复杂度上界是较为宽松的.
还没完~~
我们还可以通过预处理前面一部分的\(F\)值加快速度. 我们线性预处理前\(n^{\frac 2 3}\)的\(F\)值, 剩下的通过记忆化搜索来实现. 这样一来, 我们就只需要计算\(\lfloor \frac n k \rfloor\)中\(k \le n^\frac 1 3\)的部分, 最大计算次数可以近似地看作是
时间复杂度: \(O(n^\frac 2 3)\).
总结: 我们通过把\(F(n) = \sum_{i = 1}^n f(i)\)转化为\(g(1) F(n) = \sum_{k = 1}^n (f \times g) (k) - \sum_{k = 2}^n g(k) F( \lfloor \frac n k \rfloor )\), 通过构造一个合适的\(g\), 使得\(\sum_{k = 1}^n (f \times g) (k)\)和\(\sum_{k =1}^n g(k)\)都能快速得到, 从而用预处理前缀 + 记忆化搜索的方法将求前缀和的方法优化到\(\frac 2 3\)次的时间复杂度.
至于有些资料上说这种求和方法只能用于积性函数, 我不太认可. 实际上, 对于所有能线性预处理出前缀和的函数, 均可能可以采用以上方法; 即便不能线性预处理前缀和, 也可以用\(O(n^\frac 3 4)\)的方法优化时间复杂度.
应用
理论部分讲完了, 当然就到应用部分啦
SDOI 2015 约数个数和
Description:
共\(T \le 5 \times 10^4\)组询问, 每组询问给定\(n\)和\(m\), 请你求出
Solution:
先给出一个结论:
证明: 我们令\(i = p_1^{a_1} p_2^{a_2} \cdots\), \(j = p_1^{b_1} p_2^{b_2} \cdots\), \(d | ij\)且\(d = p_1^{c_1} p_2^{c_2} \cdots\), 则\(c_n \le a_n + b_n\).
考虑如何不重复地统计每一个\(d\): 令\(c_n = A_n + B_n\), 其中\(A_n\)和\(B_n\)分别为\(i\)和\(j\)对\(c_n\)的贡献, 则我们要求
这样一来, \(c_n\)的表示形式就变成唯一的了, 因而不会被重复统计. 我们再考虑如何统计这样的\(A_n\)和\(B_n\): 我们令\(A_n' = a_n - A_n\), 则约束条件变为
等价于\(\gcd(A_n', B_n) = 1\).
因此得证.
好吧, 假如看不懂上面的这一些证明, 就这么想吧: \(i\)表示\(a\)中不取多少, \(j\)表示\(b\)中取多少, 只要保证\(\gcd(a, b) = 1\), 就不会重复统计.
因此我们考虑原题的式子
分块处理后半部分即可.
时间复杂度: 预处理\(O(n)\), 单次询问\(O(n^\frac 1 2)\)
#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>
using namespace std;
namespace Zeonfai
{
inline int getInt()
{
int a = 0, sgn = 1; char c;
while (! isdigit(c = getchar())) if (c == '-') sgn *= -1;
while (isdigit(c)) a = a * 10 + c - '0', c = getchar();
return a * sgn;
}
}
const int N = (int)5e4, MOD = (int)1e9;
typedef int arr[N + 7];
typedef long long Larr[N + 7];
int tot;
arr isNotPrime, prm, mu, minDivisor, minDivisorDegree, sgm;
Larr a, b, c;
inline void initialize()
{
memset(isNotPrime, 0, sizeof(isNotPrime));
tot = 0;
sgm[1] = mu[1] = 1;
for (int i = 2; i <= N; ++ i)
{
if (! isNotPrime[i])
{
prm[tot ++] = i;
mu[i] = -1;
minDivisor[i] = i;
minDivisorDegree[i] = 1;
sgm[i] = 2;
}
for (int j = 0; j < tot && i * prm[j] <= N; ++ j)
{
int x = i * prm[j]; isNotPrime[x] = 1;
if (i % prm[j])
{
mu[x] = - mu[i];
minDivisor[x] = prm[j];
minDivisorDegree[x] = 1;
sgm[x] = sgm[i] * 2;
}
else
{
mu[x] = 0;
minDivisor[x] = minDivisor[i] * prm[j];
minDivisorDegree[x] = minDivisorDegree[i] + 1;
sgm[x] = sgm[i / minDivisor[i]] * (minDivisorDegree[x] + 1);
}
}
}
a[0] = b[0] = c[0] = 0;
for (int i = 1; i <= N; ++ i) a[i] = a[i - 1] + sgm[i], b[i] = a[i] * a[i], c[i] = c[i - 1] + mu[i];
}
int main()
{
using namespace Zeonfai;
initialize();
int T = getInt();
for (int cs = 0; cs < T; ++ cs)
{
int n = getInt(), m = getInt();
long long ans = 0;
int L = 1;
while (L <= min(n, m))
{
int R = min(n / (n / L), m / (m / L));
ans = ans + a[n / L] * a[m / L] * (c[R] - c[L - 1]);
L = R + 1;
}
printf("%lld\n", ans);
}
}
BZOJ 4176 Lucas的数论
Description: 给定正整数\(n \le 10^9\), 要你求
Solution: 乍一看这题跟上一题真的非常相像... 所以我们就照搬上一题的公式吧
我们令
则
注意到\(n\)的范围比较大, 考虑怎么处理\(\mu\)的前缀和以及\(s\)值.
考虑怎么求\(s(n)\):
按照\(\lfloor \frac n d \rfloor\)分块即可.
再考虑怎么求\(\sum_{k = 1}^n \mu(k)\): 按照杜教筛的一般形式
我们令\(t(n) = sum_{k = 1}^n \mu(k)\), 有
预处理前缀和 + 记忆化 + 递归即可.
时间复杂度: 计算\(t(n)\)是\(O(n^\frac 2 3)\), 计算\(s(n)\)是\(O(n^\frac 3 4)\), 因此总时间复杂度为\(O(n^\frac 3 4)\)
#include <cstdio>
#include <cstring>
#include <tr1/unordered_map>
#include <algorithm>
#include <cmath>
using namespace std;
using namespace std::tr1;
typedef long long LL;
const int BND = (int)1e6, MOD = (int)1e9 + 7;
LL n;
unordered_map<int, int> mp;
inline void initialize()
{
static int isNotPrime[BND + 7], prm[BND + 7], mu[BND + 7];
memset(isNotPrime, 0, sizeof isNotPrime);
int tot = 0; mu[1] = 1;
int bnd = round(pow(n, (double)2 / 3));
for (int i = 2; i <= bnd; ++ i)
{
if (! isNotPrime[i])
{
mu[i] = -1;
prm[tot ++] = i;
}
for (int j = 0; j < tot && i * prm[j] <= BND; ++ j)
{
isNotPrime[i * prm[j]] = 1;
if (i % prm[j] == 0)
{
mu[i * prm[j]] = 0;
break;
}
mu[i * prm[j]] = - mu[i];
}
}
int sum = 0; mp.clear();
for (int i = 1; i <= bnd; ++ i) mp[i] = (sum = sum + mu[i]);
}
inline int power(LL a) { return a * a % MOD; }
inline LL s(int n)
{
LL res = 0;
for (int i = 1; i <= n; ++ i)
{
LL L = i, R = n / (n / i);
res = (res + (n / L) * (R - L + 1)) % MOD;
i = R;
}
return res;
}
inline LL t(int n)
{
if (n == 0) return 0;
if (mp.find(n) != mp.end()) return mp[n];
LL res = 1;
for (int i = 2; i <= n; ++ i)
{
LL L = i, R = n / (n / i);
res = (res - t(n / L) * (R - L + 1) % MOD + MOD) % MOD;
i = R;
}
return mp[n] = res;
}
int main()
{
scanf("%lld", &n);
initialize();
LL ans = 0;
for (int i = 1; i <= n; ++ i)
{
LL L = i, R = n / (n / i);
ans = (ans + (t(R) - t(L - 1) + MOD) * power(s(n / L))) % MOD;
i = R;
}
printf("%lld\n", ans);
}
