线性筛的魔法改造:用因子个数统计破解球盒难题
在算法竞赛中,我们常常会遇到一些看似是组合数学问题,实则暗藏数论玄机的题目。今天要探讨的这个"球盒问题"就是典型代表——将n个球放入n个盒子,要求每个盒子里的球与其编号的因子个数相同。初看像是排列组合,实则核心在于高效计算因子个数和巧妙利用模数特性。而解决这个问题的关键钥匙,竟是我们熟悉的线性筛法。
1. 问题本质与数学洞察
让我们先抛开代码,从数学角度理解这个问题。题目要求每个盒子的编号和放入球的编号具有相同数量的因子。这意味着我们需要:
- 对所有1到n的数字计算其因子个数
- 将相同因子个数的数字分为一组
- 每组内部的排列方案数为该组大小的阶乘
- 最终答案为各组阶乘的乘积
关键观察点:
- 因子个数相同的数字之间可以任意排列
- 当模数较小时(本题为500009),存在临界点使得某个因子个数的出现次数超过模数,此时整个乘积必然为0
- 通过打表发现,当n≥2,250,000时必然存在这样的因子个数
# 简单示例:n=3时的计算过程 因子个数统计: 1: 1 (因子数1) 2: 2 (因子数2) 3: 2 (因子数2) 分组情况: 因子数为1的数字:[1] 因子数为2的数字:[2,3] 方案计算: 1! * 2! = 1 * 2 = 22. 线性筛法的深度改造
标准的线性筛法用于高效找出素数,但我们可以改造它,使其在筛素数的同时计算每个数的因子个数。这需要深入理解线性筛的工作原理并进行巧妙扩展。
2.1 标准线性筛回顾
标准线性筛的核心思想是:
- 维护一个素数列表
- 对每个数i,用已知素数p去筛i*p
- 当i能被p整除时停止,确保每个合数只被最小素因子筛掉
// 标准线性筛伪代码 vector<bool> is_prime(n+1, true); vector<int> primes; for (int i = 2; i <= n; i++) { if (is_prime[i]) primes.push_back(i); for (int p : primes) { if (i*p > n) break; is_prime[i*p] = false; if (i % p == 0) break; } }2.2 扩展为因子个数筛
要计算因子个数,我们需要记录:
- 每个数的最小素因子及其幂次
- 利用数论公式:若n = p₁^a₁ * p₂^a₂...pₖ^aₖ,则因子个数d(n) = (a₁+1)(a₂+1)...*(aₖ+1)
改造后的筛法分为三个阶段处理:
阶段一:处理p³ ≤ n的素数
- 对每个素数p,处理其所有幂次p^k
- 对每个p^k的倍数,乘以(k+1)的贡献
阶段二:处理p² ≤ n的剩余素数
- 类似阶段一,但只需处理k=1和k=2的情况
阶段三:处理剩余的大素数
- 这些素数在平方以上,每个只能贡献一个素因子,即乘以2
vector<int> ndivisors(MAXN + 1, 1); // 初始化为1,因为1没有素因子 // 阶段一:处理p³ ≤ MAXN的素数 for (int p = 2; p*p*p <= MAXN; ++p) { if (ndivisors[p] == 1) { // p是素数 int pk = p; for (int k = 1; ; ++k) { int mm = 0; for (int q = pk; q <= MAXN; q += pk) { ++mm; if (mm == p) mm = 0; // 避免重复计算 else ndivisors[q] *= (k + 1); } LL pk2 = LL(pk) * p; if (pk2 > MAXN) break; pk = int(pk2); } } }3. 模数特性的巧妙利用
本题模数为500009,这个"小模数"带来了重要的优化机会。根据威尔逊定理的扩展思想,当某个因子个数的出现次数超过模数时,其阶乘必然包含模数的倍数,因此整个乘积模500009的结果为0。
关键优化步骤:
- 预处理阶段设定上限MAXN=2,250,000
- 当n ≥ MAXN时直接返回0
- 否则使用预处理的因子个数统计结果
vector<unsigned> cnt(MAXN + 1, 0); vector<unsigned> res(MAXN + 1, 0); res[0] = 1; for (int n = 1; n <= MAXN; ++n) { int d = ndivisors[n]; ++cnt[d]; res[n] = res[n-1] * (cnt[d] % 1024u) % unsigned(MOD); if (cnt[d] >= 1024u) { unsigned a = res[n-1] * (cnt[d] >> 10) % unsigned(MOD); res[n] = (res[n] + (a << 10)) % unsigned(MOD); } if (res[n] == 0) break; // 提前终止 }4. 性能优化与实现细节
为了处理T≤1e5的大规模查询,我们需要:
- 预处理所有结果:预先计算1到MAXN的所有答案
- 高效查询:O(1)时间响应每个查询
- 内存优化:使用紧凑的数据结构存储中间结果
阶乘计算的优化技巧:
- 将计数器cnt[d]拆分为低10位和高位部分
- 分别计算贡献后再合并,避免大数运算
- 一旦结果变为0即可提前终止后续计算
注意:在实际编码中,位运算技巧可以显著提升模运算效率,但需要确保不会导致数值溢出。
5. 从数学到代码的完整实现
将上述洞察转化为完整解决方案,我们需要:
- 初始化因子个数数组ndivisors
- 三阶段筛法计算每个数的因子个数
- 统计每个因子个数的出现次数cnt[d]
- 预处理答案数组res[n]
- 处理查询
完整解决方案框架:
#include <vector> using namespace std; const int MOD = 500009; const int MAXN = 2250000; vector<int> precompute() { vector<int> ndivisors(MAXN + 1, 1); // 三阶段筛法计算因子个数 // ... (省略具体实现) vector<unsigned> cnt(MAXN + 1, 0); vector<unsigned> res(MAXN + 1, 0); res[0] = 1; for (int n = 1; n <= MAXN; ++n) { int d = ndivisors[n]; ++cnt[d]; // ... (阶乘计算优化) if (res[n] == 0) break; } return res; } int main() { auto res = precompute(); // 处理查询 int T, n; scanf("%d", &T); while (T--) { scanf("%d", &n); printf("%d\n", n <= MAXN ? res[n] : 0); } return 0; }6. 同类问题的扩展思考
这种基于线性筛的扩展技术可以应用于多种数论问题:
- 因子和计算:类似因子个数,但改为计算(p^(a+1)-1)/(p-1)的乘积
- 欧拉函数计算:记录每个数的最小素因子来递推计算
- 莫比乌斯函数:在筛法过程中标记平方因子
性能对比表:
| 算法类型 | 预处理时间 | 单次查询时间 | 适用场景 |
|---|---|---|---|
| 暴力计算 | O(1) | O(√n) | 少量查询 |
| 普通筛法 | O(n logn) | O(1) | 中等规模n |
| 改造线性筛 | O(n) | O(1) | 大规模n |
在实际比赛中,这种对经典算法的深度理解和灵活改造能力,往往是解决难题的关键。当遇到模数特殊的问题时,不妨多思考模数本身的性质可能带来的优化机会。