下文中所有讨论都在数论函数范围内开展.

数论函数指的是定义域为正整数域, 且值域为复数域的函数.

数论意义下的和式处理技巧

1. 因子

   \[\sum_{d | n} a_d = \sum_{d | n} a_{\frac n d}
   \]

2. 双重因子

1. 遍历 + 因子

1. \(\mu\)函数相关(\(\mu\)函数后面会讲到, 这里先放公式)

莫比乌斯反演

引入

现给定一个函数

根据\(F\)的定义, 我们有

假设我们已知所有\(F\)值, 考虑如何倒推出\(f\)值:

倒推的过程是否有什么通法呢?

实际上是有的, 这里就用到了Mobius反演:

Mobius函数

我们注意到, 这个式子中出现了函数\(\mu\). \(\mu(d)\)为Mobius函数, 其定义如下:

我们尝试列举出\(\mu\)的前几个值

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\)分为两部分来看:

1. \(k \le \sqrt{n}\), 我们假设每个\(\lfloor \frac n k \rfloor\)都是不同的, 则这一部分的计算次数可以近似地看作是

​ 因此时间复杂度不超过\(O(n^\frac 3 4)\)

1. \(k > \sqrt{n}\), 则有\(\lfloor \frac n k \rfloor \le \sqrt{n}\). 我们假设这\(\sqrt{n}\)个值都出现一次, 则计算次数可以近似地看作是

2. \[\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);

}