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 个难点：

1. 一个数的 21 次幂极为庞大，即使 unsigned long 也显得无能为力；
2. 即使是 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 的开发效率非常高，而且强制的缩进，使得不管是写代码的人还是看代码的人都非常清楚。人生只有短短几十年，开发效率低无疑浪费生命。如果不是对运行效率有更为严苛的要求，试着让自己轻松一些吧！