| lemonoil

blog源地址 感谢__debug大佬。

洲阁筛

一种在 \( O(\frac{n^{\frac{3}{4}}}{\log n})\) 的时间中求出大多数积性函数的前缀和的方法.

引入 求 \( \sum_{i = 1}^{n} F(i) \) 其中 \( F(x) \) 是一个积性函数. 要求在低于线性的时间内求出.

转化 如果 \( F(x) \) 是一个十分特殊的积性函数, 比如 \( \mu^2(x) \), 那么可以非常快地求出; 就算它不那么特殊, 如果我们能找到两个合适的函数 \( G, H \) 满足\( F * G = H \), 并且\( G, H \) 的前缀和容易求出, 那么可以利用杜教筛在 \( O(n^{\frac{2}{3}}) \)的时间内求出. 但是当\( F(x)\) 不具备这些性质的时候呢?

当然了, 还是有一个要求: 当 \( p\) 为质数的时候, \( F(p^c) \) 是一个关于 (p) 的低阶多项式.

我们将 \( [1, n]\)的所有数按照是否有 \( > \sqrt{n}\) 的质因子分为两类, 那么显然有: \( \sum_{i = 1}^{n} F(i) = \sum_{\substack{1 \le i \le n \ i \text{ have no prime factors } > \sqrt{n} }} F(i) \left( 1 + \sum_{\substack{\sqrt{n} < j \le \lfloor \frac{n}{i} \rfloor \ j \text{ is prime}}} F(j) \right) \)

由于 \( i \ge \sqrt{n}\) 时后面一部分等于 \(1\), 所以需要计算以下两个东西:

对于每个 \(1 \le i < \sqrt{n}\), 计算 \( \sum_{\substack{\sqrt{n} < j \le \lfloor \frac{n}{i} \rfloor \ j \text{ is prime}}} F(j) \sum_{\substack{\sqrt{n} \le i \le n \ i \text{ have no prime factors } > \sqrt{n} }} F(i) \) 当然, 最后还要用线性筛求出 \([1, \sqrt{n})\)的 \(F(x)\) 乘上相应的系数对答案的贡献, 这一部显然不会成为瓶颈.

Part 1

设 \(g_k(i, j)\) 表示 \([1, j]\) 中与前 \(i\) 个质数互质的数的 \(k\) 次幂和. 显然有转移 \(g_k(i, j) = g_k(i - 1, j) - p_i^k g_k(i - 1, \lfloor \frac{j}{p_i} \rfloor) \)

观察到 \(j\) 的取值只有 \(\sqrt{n}\) 种, 于是直接暴力计算的复杂度为 \(O(\frac{n}{\log n})\).

考虑优化. 首先注意到当 \(p_{i + 1} > j\) 时 \(g_k(i, j) = 1\), 但是这个观察并不能带来什么优化; 紧接着我们发现如果 \(p_i > \lfloor \frac{j}{p_i} \rfloor\) 即 \(p_i^2 > j\) 时, \(g_k(i, j)\) 的转移变为: \(g_k(i, j) = g_k(i - 1, j) - p_i^k\)

我们从小到大枚举 \(i\), 对于某个 \(j\) 一旦 \(p_{i_0}^2 > j\) 便可以不再转移, 之后如果其他的值需要使用到它在 \(i_1\) 时的值, 直接用 \(g_k(i_0, j) - \sum_{l = i_0}^{i_1 - 1} p_l^k\) 即可. (其实这个地方还有一些细节, 可以自己思考)

此时的复杂度可以简单地用积分近似为 \(O(\frac{n\frac{3}{4}}{\log n})\).

Part 2

设 \(f(i, j)\) 表示 \([1, j]\) 中仅由前 \(i\) 个质数组成的数的 \(F(x)\) 之和. 显然有转移 \(f(i, j) = f(i - 1, j) + \sum_{c \ge 1} F(p_i^c) f(i - 1, \lfloor \frac{j}{p_i^c} \rfloor)\)

虽然多出来一个 \(\sum_{c \ge 1}\), 但是直接暴力计算的复杂度仍然是 \(O(\frac{n}{\log n})\). 如何优化? 似乎不太能沿用之前的方法了.

但是如果将状态定义中的"前 \(i\) 个质数"改为"(小于 \( \sqrt{n}\)的后 \(i\) 个质数", 此时当 \(p_i > j\) 时, 一定有 \(f(i, j) = 1\). 类似地, 当 \(p_i^2 > j\) 时转移变为: \(f(i, j) = f(i - 1, j) + F(p_i) \)

所以可以从大到小枚举 \(i\), 如果对于某个 \(j\) 有 \(p_i^2 > j\), 可以不转移, 每次用的时候加入 \([p_i, \min(j, \sqrt{n})]\) 这一段的质数的 \(F(p)\) 就可以了.

类似上面的分析, 复杂度为 \(O(\frac{n^{\frac{3}{4}}}{\log n})\).

小结 从上面的推导可以看出, 如果 \(F(p^c)\) 是一个关于 \(p\) 的常数阶多项式, 那么我们可以在 \(O(\frac{n^{\frac{3}{4}}}{\log n})\) 的时间内求出 \(F(x)\) 的前 \(n\) 项和. 利用滚动数组, 空间复杂度是 \(O(\sqrt{n})\) 的.

主要用到了以下几个想法:

将所有数分按照有无大于 \(\sqrt{n}\) 的质因子为两类 整除的结果只有 \(O(\sqrt{n})\) 种 一个关于递推的优化技巧 这几个想法相对独立, 都挺有启发性的, 感觉甚至可以优化一些 DP.

例题: SPOJ DIVCNT3 Description 求 \( \sum_{i = 1}^{n} d(i^3)\)

其中 \(d(x)\) 表示 \(x\) 的约数个数.

\( n \le 10^{11} \) , 时限 20s.

Solution 直接上洲阁筛, 其中 \(F(p^c) = 3c + 1\).

提供主要代码以供参考. 其中 \(sump[i]\) 表示小于等于 \(i\) 的质数之和, \(d3[i]\) 表示 \(d(i^3)\).

LL N;
int pbnd;
int vbnd;
int l0[SIZE], l[SIZE];
LL g0[SIZE], g[SIZE], f0[SIZE], f[SIZE];

void calcG()
{
    for (int i = 1; i < vbnd; ++i) {
	g0[i] = i;
    }
    for (int i = vbnd - 1; i >= 1; --i) {
	g[i] = N / i;
    }
    for (int i = 0; i < pbnd; ++i) {
	int p = primes[i];
	for (int j = 1; j < vbnd && i < l[j]; ++j) {
	    LL y = (N / j) / p;
	    g[j] -= (y < vbnd ? g0[y] - std::max(0, i - l0[y]) : g[N / y] - std::max(0, i - l[N / y]));
	}
	for (int j = vbnd - 1; j >= 1 && i < l0[j]; --j) {
	    LL y = j / p;
	    g0[j] -= g0[y] - std::max(0, i - l0[y]);
	}
    }
    for (int i = 1; i < vbnd; ++i) {
	g[i] -= pbnd - l[i];
    }
}

void calcF()
{
    std::fill(f0 + 1, f0 + vbnd, 1);
    std::fill(f + 1, f + vbnd, 1);
    for (int i = pbnd - 1; i >= 0; --i) {
	int p = primes[i];
	for (int j = 1; j < vbnd && i < l[j]; ++j) {
	    LL y = (N / j) / p;
	    for (int z = 1; y; ++z, y /= p) {
		f[j] += (3 * z + 1) * (y < vbnd ?
				       f0[y] + 4 * std::max(0, sump[y] - std::max(l0[y], i + 1)) :
				       f[N / y] + 4 * (pbnd - std::max(l[N / y], i + 1)));
	    }
	}
	for (int j = vbnd - 1; j >= 1 && i < l0[j]; --j) {
	    int y = j / p;
	    for (int z = 1; y; ++z, y /= p) {
		f0[j] += (3 * z + 1) * (f0[y] + 4 * std::max(0, sump[y] - std::max(l0[y], i + 1)));
	    }
	}
    }
    for (int i = 1; i < vbnd; ++i) {
	f[i] += 4 * (pbnd - l[i]);
    }
}

LL calcSumS3()
{
    for (vbnd = 1; (LL)vbnd * vbnd <= N; ++vbnd) { }
    for (pbnd = 0; (LL)primes[pbnd] * primes[pbnd] <= N; ++pbnd) { }
    for (int i = 1; i < vbnd; ++i) {
	for (l0[i] = l0[i - 1]; (LL)primes[l0[i]] * primes[l0[i]] <= i; ++l0[i]) { }
    }
    l[vbnd] = 0;
    for (int i = vbnd - 1; i >= 1; --i) {
	LL x = N / i;
	for (l[i] = l[i + 1]; (LL)primes[l[i]] * primes[l[i]] <= x; ++l[i]) { }
    }

    calcG();
    calcF();

    LL ret = f[1];
    for (int i = 1; i < vbnd; ++i) {
	ret += d3[i] * 4 * (g[i] - 1);
    }
    return ret;
}