3 分钟内找到全部的 21 位水仙花数
一、问题描述
一个 N 位的十进制正整数，如果它的每个位上的数字的 N 次方的和等于这个数本身，则称其为水仙花数。
例如：当 N=3 时，153 就满足条件，因为 \(1^3 + 5^3 + 3^3 = 153\)；当 N=4 时，1634 满足条件，因为 \(1^4 + 6^4 + 3^4 + 4^4 = 1634\)；当 N=5 时，92727 满足条件。实际上，对 N 的每个取值，可能有多个数字满足条件。
程序的任务是：求 N=21 时，所有满足条件的花朵数。
注意：这个整数有 21 位，它的各个位数字的 21 次方之和正好等于这个数本身。如果满足条件的数字不只有一个，请从小到大输出所有符合条件的数字，每个数字占一行。因为这个数字很大，请注意解法时间上的可行性。要求程序在 3 分钟内运行完毕。
二、难点与思路
这个问题用 C 语言实现主要面临 2 个难点：
- 一个数的 21 次幂极为庞大，即使 unsigned long 也显得无能为力；
- 即使是 C 语言这样运行效率极高的语言，穷举判断也很有可能超出 3 分钟的时限。
1、借助数组逐位运算
为了进一步减少程序运行时间，我们可以预存 0 到 9 的 21 次幂的结果，之后采用查表的方法求各位数的 N 次幂之和。
#include"stdio.h"
void pow21(int b) // h 高 7 位数，m 中 7 位数，l 低 7 位数
{
unsigned long h=0, m=0, l=b, i = 1;
for(; i < 21 ; i++)
{
l = l * b;
m = m * b + l/10000000; // 中位加低位的进位
h = h * b + m/10000000; // 高位加中位的进位
l %= 10000000; // 低位保留 7 位数据
m %= 10000000; // 中位保留 7 位数据
}
printf("%07ld%07ld%07ld\n",h,m,l);
}
void main(void)
{
int i = 0;
for(; i < 10 ; i++)
pow21(i); // 求 i 的 21 次幂
}
000000000000000000000
000000000000000000001
000000000000002097152
000000000010460353203
000000004398046511104
000000476837158203125
000021936950640377856
000558545864083284007
009223372036854775808
109418989131512359209
2、搜索水仙花数的方法
从 \(x_0=9\times10^{20}\) 开始，通过 \(y \cdot V\) 之间的关系判定所属情况（如下表所列），然后搜寻下一个 \(x\)，具体步骤为：
- 从高位开始依次寻找 \(x\) 中首个为 0 的位 \(b_k\)：若为情况 2，则将高一位赋值于它 \(b_k= b_{k+1}\)；若为情况 3，则将高一位减 1，即 \(b_{k+1}=b_{k+1}-1\)；
- 如果 \(x\) 不存在 0，则从高位开始依次寻找 \(x\) 中首个为 1 的位 \(b_k\)，将高一位减 1，即 \(b_{k+1}=b_{k+1}-1\)，同时清零 \(b_i,\,i<=k\)；
- 如果 \(x\) 不存在 0 和 1，则将 x 的末位 \(b_0\) 减小 1。
显然，按上述方法操作 \(x\) 低位的数总不能大于高位的数，考虑下表情况 1 可知：程序终止条件为 \(b_{20}<8\)。
| 情况 | \(x=\displaystyle\sum_{i=0}^{20}b_i\times10^i\) | \(y=f(x)=\displaystyle\sum_{i=0}^{20}b_i^{21}\) | \(V=(10^{20},10^{21})\) | 说明 |
|---|---|---|---|---|
| 1 | \(\forall\, b_i < 8\) | \(\max y = 21\times7^{21}\) | \(y \notin V,\, y<10^{20}\) | 若 \(x\) 不含数字 8 或 9，则 \(y\) 一定不足 21 位 |
| 2 | \(9\times10^{20}\) | \(109418989131512359209\) | \(y \in V\) | 若 \(x,y\) 只是数字的排列顺序不同，则 \(y\) 为水仙花数 |
| 3 | \((10^{10}-1)\times10^{10}\) | \(1094189891315123592090\) | \(y \notin V ,\, y>10^{21}\) | 若 \(y\) 超出 \(V\) 的上界，则改变 \(x\) 使 \(y\) 减小 |
再详细解释下为何这样搜索 \(x\)（以 3 位数为例）：
- 譬如 \(x=900\) 开始，\(y\) 也是三位数，但不是水仙花数，就应该找下一个 \(x=990\)，即尽可能保证改变一个数位后 \(y\) 增加最大，这样可以保证覆盖到所有可能的数；
- 由于 \(f(990)\) 超出了 \(V\) 的范围，因此缩小它到 \(x=980\)，继续，直到 \(x=960\) 又回到 \(V\) 的范围；
- 因为 990，980，970 都在 960 之前判断过，因此改为判断 \(x=966\)，继续，直到 \(x=961\)，还在 \(V\) 的范围内；
- 因为 960 也已经判断过，所以接下来该判断 \(x=950\)，还在 \(V\) 的范围内，继续，直到 \(x=740\)，求得 \(y=407\)；
- 显然 \(x\) 与 \(y\) 只是数字的排列顺序不同，因此 \(y\) 是水仙花数，继续，直至 \(x\) 的最高位为 3；
- 由于 \(x\) 低位的数不会大于高位，因此 \(y\) 的最大值也不超过三位数，判断结束。
三、C 程序源代码及其运行结果
#include<stdio.h>
#include<time.h>
#define N 21
typedef unsigned int INT;
INT tab[10][N]={{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,9,7,1,5,2},{0,0,0,0,0,0,0,0,0,0,1,0,4,6,0,3,5,3,2,0,3},{0,0,0,0,0,0,0,0,4,3,9,8,0,4,6,5,1,1,1,0,4},{0,0,0,0,0,0,4,7,6,8,3,7,1,5,8,2,0,3,1,2,5},{0,0,0,0,2,1,9,3,6,9,5,0,6,4,0,3,7,7,8,5,6},{0,0,0,5,5,8,5,4,5,8,6,4,0,8,3,2,8,4,0,0,7},{0,0,9,2,2,3,3,7,2,0,3,6,8,5,4,7,7,5,8,0,8},{1,0,9,4,1,8,9,8,9,1,3,1,5,1,2,3,5,9,2,0,9}};
void main(void)
{
int i=0, j=0, c=0, v=0, flag=0;
INT x[N]={9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, y[N]={0}, temp[N]={0};
clock_t begin = clock();
double cost;
while(x[0]>7)
{
for(i=0; i<N; i++)
y[i] = 0; // 清零数组 y
for(i=0; i<N; i++)
for(j=0; j<N; j++)
y[i] += tab[x[j]][i]; // 对数组 x 中的每一位分别求和
for(i=N-1; i>0; i--)
{
c = y[i]/10;
y[i] %= 10;
y[i-1] += c;
} // 对数据进行进位处理
v = y[0] > 0 && y[0] < 10; // y 是否是一个 21 位的数
if(v)
{
for(i=0; i<N; i++)
temp[i]=x[i]; // 把 y 临时存储于 temp
for(i=0; i<N; i++)
for(j=0; j<N; j++)
if(y[i]==temp[j])
{ // 在数组 temp 中搜寻数组 y 的每一位
temp[j]=10;
break;
} // 将数组 temp 中首个匹配的位标记为 10
flag=0;
for(i=0; i<N; i++)
flag += temp[i]; // 计算数组 temp 各位的和
if(flag==N*10) // 和为 210 则证明数组 y 代表着一个水仙花数
{
for(i=0;i<N;i++)
printf("%d",y[i]);
printf("\n"); // 这里直接输出结果，否则保存起来排序后再输出
}
}
flag = N;
for(i=0; i<N; i++) // 寻找数组 x 中是否有 0
if(x[i]==0)
{
flag = i;
break;
} // flag 为数组 x 中第一个 0 的索引
if(flag < N)
{
if(v) x[flag] = x[flag-1]; // 数组 x 中 flag 之前的一位赋值给 flag 位
else x[flag-1] -= 1; // 数组 x 中 flag 之前的一位自减 1
}
else
{
flag = N;
for(i=0; i<N; i++) // 寻找数组 x 中是否有 1
if(x[i]==1)
{
flag = i;
break;
} // flag 为数组 x 中第一个 1 的索引
if(flag < N)
{
x[flag-1] -= 1;
for(i=flag; i<N; i++)
x[i] = 0;
} // 数组 x 中 flag 之前的一位自减 1，之后的每一位都置零
else x[N-1] -= 1; // 如果数组 x 中不含 0 或 1，则将其末尾自减 1
}
}
cost = (double)(clock() - begin) / CLOCKS_PER_SEC;
printf("%lf seconds\n", cost);
}
449177399146038697307
128468643043731391252
24.890000 seconds
四、Python 实现方案及其运行结果
import time
time_count = time.clock()
p, x, flower= [i**21 for i in range(10)], [int(i) for i in str(int(9*1e20))], []
f = lambda x: sum(p[b] for b in x)
while(x[0] > 7):
y = f(x)
v = 1e20 < y < 1e21
if v and sorted(x) == sorted(int(i) for i in str(y)): flower.append(y)
try: k = x.index(0)
except:
try: k = x.index(1)
except: x[-1] -= 1
else:
x[k-1] -= 1
for j in range(k,len(x)): x[j] = 0
else:
if v: x[k] = x[k-1]
else: x[k-1] -= 1
print(sorted(flower),'\nTime : ',int(time.clock()-time_count),'s')
[128468643043731391252, 449177399146038697307]
Time : 167 s
比 C 语言多耗费了 5.7 倍的执行时间，不过好在勉强完成了任务。
人生苦短，我用 Python。Pythoner 都以这句话引以为傲，因为 Python 的开发效率非常高，而且强制的缩进，使得不管是写代码的人还是看代码的人都非常清楚。人生只有短短几十年，开发效率低无疑浪费生命。如果不是对运行效率有更为严苛的要求，试着让自己轻松一些吧！
