那 些 回 忆 丶 挥 之 不 去 Cnblogs Pandaman Dashboard

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 的开发效率非常高,而且强制的缩进,使得不管是写代码的人还是看代码的人都非常清楚。人生只有短短几十年,开发效率低无疑浪费生命。如果不是对运行效率有更为严苛的要求,试着让自己轻松一些吧!