什么是完全数?

完全数(Perfect number), 又称完美数或完备数, 是一些特殊的自然数。它所有的真因子（即除了自身以外的约数）的和,
恰好等于它本身。如果一个数恰好等于它的因子之和, 则称该数为“完全数”。

以上来自百度百科.

完全数的定义理解起来并不难. 让我们试着做一两道题吧.

一、HDU 1406 完数

Problem Description

完数的定义：如果一个大于1的正整数的所有因子之和等于它的本身，则称这个数是完数，比如6，28都是完数：6=1+2+3；
28=1+2+4+7+14。

本题的任务是判断两个正整数之间完数的个数。

Input

输入数据包含多行，第一行是一个正整数n，表示测试实例的个数，然后就是n个测试实例，每个实例占一行，由两个正整数
num1和num2组成,(1<num1,num2<10000) 。

Output

对于每组测试数据，请输出num1和num2之间（包括num1和num2）存在的完数个数。

Sample Input

2
2 5
5 7

Sample Output

0
1

题目解析

    1. 多次访问[num1, num2]之间的数有多少个完全数, 所以一开始应该先求出2 ~ 9999的完全数.
    2. 范围: 1 < num1, num2 < 10000, 而num1不一定小于num2, 计算数量时需注意num1, num2的大小.

首先解决最简单的输入输出.

#include <stdio.h>
int perfect[6];             //存储完全数
int num;                    //存储完全数的个数
void swap(int* a, int* b) {
    int t = *a; *a = *b; *b = t;
}
//求2 ~ N的完全数的函数
void solve(int N) { ... }
int main() {
    solve(9999);            //调用该函数
    int n, a, b;
    scanf("%d", &n);        //读入数据的组数
    while (n--) {
        scanf("%d%d", &a, &b);  //读入a, b
        if (a > b)              //若a > b, 交换它们
            swap(&a, &b);
        int ans = 0;            //记录[a, b]的完全数的个数
        //寻找在[a, b]之间的完全数
        for (int i = 0; i < num && perfect[i] <= b; ++i)
            if (perfect[i] >= a) ++ans;
        printf("%d\n", ans);    //输出结果
    }
    return 0;
}

现在考虑如何求解2 ~ N 的完全数:

对于一个数$m \geq 2$, 若有两因子$p, q(p \leq q)$, 使$m = p * q$, 那么必然有:
$m = p \bullet q \geq p^{2}$
$\therefore p \leq \sqrt{m}$
因此, 当m的因子$p < \sqrt{m}$时, 表示$p, m / p$为$m$其中两个因子, 当m的因子$p = \sqrt{m}$时, 表示$p$为$m$的因子.
所以, 可以写一个暴力算法求解任意正整数是否为完全数:

//判断N是否为完全数
bool isPerfect(int N) {
    if (N == 1) return false;   //1不是完全数
    int i;
    int s = 1;                  //存储因子和, 初始值为1(N > 1时, 1为该数的因子)
    for (i = 2; i * i < N; ++i)
        if (N % i == 0)         //若i < sqrt(N)且i | N, 那么i, N / i为N的因子
            s += i + N / i;
    if (i * i == N) s += i;     //若i = sqrt(N)且i | N, 那么i为N的因子
    return s == N;              //当且仅当因子和等于N时, N为完全数
}
//求2 ~ N的完全数的函数
void solve(int N) {
    num = 0;
    for (int i = 2; i <= N; ++i)
        if (isPerfect(i)) perfect[num++] = i;   //当i为完全数, 存储进perfect[]里
}

$isPerfect()$判断$N$是否为完全数, 变量$i$从$2$遍历至$\sqrt{N}$, 所以它的时间复杂度为: $O(\sqrt{N})$
由于要判断$2~N$共$N-1$个数是否为完全数, 每次判断需要$O(\sqrt{N}$)的时间. 因此总的时间复杂度为: $O(N\sqrt{N})$
HDU1406的$N$的范围不大, 还能接受.
若N为$10^{18}$, 即使电脑每秒能处理$10^{14}$的运算量, 求出$10^{18}$内的所有完全数至少也需要30万年……所以, 必须得想一个更为高效的算法.


古希腊数学家欧几里得在名著《几何原本》中论述完全数时提出：
如果$2^p-1$是素数, 则$2^{p - 1}\bullet(2^{p}-1)$是完全数.

最终, 欧拉证明欧几里得的公式可以给出所有的偶完全数, 即:
如果n是偶完全数, 则n为 $n=2^{p - 1}\bullet(2^p-1)$ 的形式,
其中$2^{p}-1$为梅森素数(梅森素数, 是指形如$2^p-1$的素数, 其中指数p也是素数).

也就是说, 只要$2^{p}-1$为梅森素数, 那么$2^{p - 1}\bullet(2^p-1)$必然是偶完全数. 剩下的问题就是判断p为素数时$2^{p}-1$是否为梅森素数了.
判断一个数是否是素数有很多种方法, 如试除法(时间复杂度为: $O(\sqrt{N})$), 拉宾-米勒素性检验(时间复杂度为: $O(k \bullet log{N})$, $k$为试验次数, 通过一次测试后, n为合数的概率少于$1/4$), 然而试除法太慢了, 拉宾-米勒素性检验算法是挺高效的, 然而它是通过概率去判断该素数是否为素数, 有没有更好的方法呢？
对于形如$2^{p}-1$梅森数, 我们可以通过一种高效且确定的算法:

卢卡斯-莱默检验法(Lucas–Lehmer primality test).

利用卢卡斯-莱默检验法可以得到下面的算法: (只要N在long long范围都能算出答案)

long long perfect[8];       //存储完全数
int num;                    //存储完全数的个数
//试除法判断小整数p是否为素数
bool isPrime(int N) {
    if (N < 4) return N > 1;
    if (!(N & 1)) return false;
    for (int i = 3; i * i <= N; i += 2)
        if (!(N % i)) return false;
    return true;
}
//Lucas–Lehmer primality test
bool primality(int N, long long M) {
    if (N == 2) return true;    //Lucas–Lehmer只能判断N>2的情况, 而N为2时2^2-1确实为梅森素数
    long long s = 4;
    for (int i = 0; i < N - 2; ++i)
        s = (s * s - 2) % M;
    return s == 0;              //N>2时, 当且仅当s为0, 该梅森数为素数
}
//求2 ~ N的完全数的函数
void solve(long long N) {
    num = 0;
    for (int i = 2; i < 32; ++i) {
        long long M = (1ll << i) - 1;   //梅森数: 2^p-1
        long long t = M << (i - 1);     //2^(p-1)*(2^p-1)
        if (t > N) break;   //t超过N的范围, 停止计算
        //当p为素数且梅森数2^p-1为素数时2^(p-1)*(2^p-1)为完全数
        if (isPrime(i) && primality(i, M))
            perfect[num++] = t;         //储存完全数t
    }
}

该算法的时间复杂度如何呢?
首先, 每个$N=2^{p}-1$需要判断p是否为素数(时间复杂度为$O(\sqrt p)$), 若$p$为素数, 再判断$2^{p}-1$是否为梅森素数(时间复杂度为$O(p)$). 对于每个数N时间复杂度为$O(p)$即$O(log{N})$, 所以对于求解2 ~ N的完全数, 该算法的复杂度为$O((log{N})^{2})$. 因此, 该算法能快速得到答案, 即使N为$10^{18}$.