【转】如何用 C/C++ 求 1 到 1000 内的所有完全数

作者：刘海洋
链接：https://www.zhihu.com/question/26494835/answer/33451991
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

这是个小题目，按定义穷举算是基本的练习，并太多特别值得说的地方。选个合适的语言，也就是一两行代码按定义表示出来，比如 Python：

[n for n in range(1, 1000)
    if sum(k for k in range(1, n) if n%k == 0) == n]

另一方面这又是个数学题目。可以想见会有人想出一些办法去做「优化」，却轻视编程码代码之前的数学与算法上的分析。实际上，如果未经深思熟虑，这类优化能起到的效果甚微，在本质上不能减少多少运算。因此，倒不如仔细从数学角度好好看看这个问题，给出一个完全不一样的算法，或许会让这个原本不起眼的小题目吃起来更有嚼头一些。

---

由 Euclid-Euler 定理，一个数 n 是偶完全数，当且仅当 ，其中 是一个 梅森素数。

奇完全数是否存在，是当今数论中未解决的问题。但 已有文献指出，奇完全数即使存在，则不小于 ，因此在合理范围内无须考虑。

于是现在问题归结到求一组较小的 k 使 是梅森素数，然后套用 Euclid-Euler 定理公式计算。

验证梅森素数的方法，较为有效的是 Lucas–Lehmer 检验法。即定义

从而 是素数当且仅当

其中 k 是奇数，由梅森素数的性质，等价于 k 是奇素数。k 为偶数 2 则是特殊情形， 是素数。

为方便计算，事实上并不需要求出所有的 ，只要求得余数即可（否则中间结果会很大）。而这个条件用计算机非常容易验证。

由于只要求 1000 以内的完全数，因此用 n 的表达式粗略地估计，最多只要求出 1000 以内的梅森素数。注意到 2 的 10 次方即超出 1000，因此只需要考虑 k 小于 10 的情形进行验证。因此这也是很有效的算法，远好于用定义试除计算。

下面是代码：

#include <iostream>

inline long mersenne(int k)
{
	return (1L << k) - 1L;
}

// euclid_euler(k) is perfect if mersenne(k) is prime
// http://en.wikipedia.org/wiki/Euclid%E2%80%93Euler_theorem
inline long euclid_euler(int k)
{
	return ((1L<<k) - 1L) << (k-1);
}

// return true if mersenne(k) is prime
bool lucas_lehmer_test(int k)
{
	// special case: merssene(2) = 3 is prime
	if (k <= 2)
		return k == 2;
	// Lucas-Lehmer test
	// http://mathworld.wolfram.com/Lucas-LehmerTest.html
	long s = 4;
	long m = mersenne(k);
	for (int i = 0; i < k - 2; ++i) {
		s = s * s - 2;
		if (m != 0)
			s %= m;
	}
	return s == 0;
}

int main()
{
	const long perfect_number_bound = 1000;
	for (int k = 0; euclid_euler(k) < perfect_number_bound; ++k) {
		if (lucas_lehmer_test(k))
			std::cout << euclid_euler(k) << '\n';
	}
	return 0;
}
// Result:
//   6
//   28
//   496