下文中所有讨论都在数论函数范围内开展.
数论函数指的是定义域为正整数域, 且值域为复数域的函数.
数论意义下的和式处理技巧
-
因子
[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 imes n' prod_{p为n'质因子}frac {p - 1} p = p imes phi(n - 1)).
- (p mid n'), 则由于(phi)为积性函数, 因而(phi(n) = phi(p) imes 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 imes g = g imes f).
- 结合律: (f imes g imes h = f imes (g imes h)).
- 分配律: (f imes (g + h) = f imes g + f imes h).
- 单位元: (f imes epsilon = f).
- 封闭性: 假如(f), (g)均为积性函数, 则(f imes g)也是积性函数.
主要形式
设(f(n))是一个数论函数, 我们需要计算(F(n) = sum_{i = 1}^n f(i)).
我们取另一个函数(g), 则有
因此
假如我们可以(O(1))求出(sum_{k = 1}^n (f imes g) (k))以及(sum_{k =1}^n g(k)), 则我们的主要时间复杂度在于递归计算(F(lfloor frac n k floor)). 递归计算每一个(F(lfloor frac n k floor)), 并进行记忆化.
对于一个(n), 考虑(lfloor frac n k floor)有多少个值:
- (1 le k le sqrt{n}), 最多(sqrt{n})个值.
- (sqrt n < k le n), 则(lfloor frac n k floor le sqrt n), 因而最多(sqrt n)个值.
因此对于每一个(n), 在计算(F(n))时我们要作(O(sqrt n))次递归运算.
考虑
因此对于某个(n), 在计算(F(n))时, 衍生出的所有需要计算的不同(F(n'))个数仍然为(lfloor frac n k floor)的不同个数个, 即(O(sqrt n))个. 记忆化递归即可.
时间复杂度的计算有一点麻烦:
我们把(lfloor frac n k floor)分为两部分来看:
- (k le sqrt{n}), 我们假设每个(lfloor frac n k floor)都是不同的, 则这一部分的计算次数可以近似地看作是
因此时间复杂度不超过(O(n^frac 3 4))
-
(k > sqrt{n}), 则有(lfloor frac n k floor 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 floor)中(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 imes g) (k) - sum_{k = 2}^n g(k) F( lfloor frac n k floor )), 通过构造一个合适的(g), 使得(sum_{k = 1}^n (f imes g) (k))和(sum_{k =1}^n g(k))都能快速得到, 从而用预处理前缀 + 记忆化搜索的方法将求前缀和的方法优化到(frac 2 3)次的时间复杂度.
至于有些资料上说这种求和方法只能用于积性函数, 我不太认可. 实际上, 对于所有能线性预处理出前缀和的函数, 均可能可以采用以上方法; 即便不能线性预处理前缀和, 也可以用(O(n^frac 3 4))的方法优化时间复杂度.
应用
理论部分讲完了, 当然就到应用部分啦
SDOI 2015 约数个数和
Description:
共(T le 5 imes 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
", ans);
}
}
BZOJ 4176 Lucas的数论
Description: 给定正整数(n le 10^9), 要你求
Solution: 乍一看这题跟上一题真的非常相像... 所以我们就照搬上一题的公式吧
我们令
则
注意到(n)的范围比较大, 考虑怎么处理(mu)的前缀和以及(s)值.
考虑怎么求(s(n)):
按照(lfloor frac n d floor)分块即可.
再考虑怎么求(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
", ans);
}