零基础学Python之42：水仙花数

前面41节《零基础学Python》的内容请点【这里】。

今天解决Python初学者经常会遇到的一道练习题：

请打印所有的水仙花数，所谓水仙花数，是指这样的3位数，它的每个位上的数字的3次方之和等于它本身，例如：153就是一个水仙花数，因为：1^3 + 5^3+ 3^3 = 153。另外还有三个，分别是：370、371和407。

先不要看答案，也不要百度，自己试着写一下程序。

直接解法

遇到这样的问题，需要尝试用《怎样解题》中的办法将问题简化和分解成这样一些子问题，当把这些子问题都解决之后，整个问题也就迎刃而解：

求出某个数m在个位、十位和百位上的的数字
求出某个数m在每个位上的数字的3次方之和
列出所有的三位数

第一步，求一个3位数在个位、十位和百位上的数，可以用除法和取余解决。

m  = 153print(m % 10)print(int(m / 10) % 10)print(int(m / 100))

第二步，一个数的3次方可以用**运算符，再求和：

s = (m % 10) ** 3 + int((m / 10) % 10) ** 3 + int(m / 100) ** 3

第三步，所有的三位数就是从100到999，利用range()函数，注意range()函数里有一个差1问题，可以用这个代码验证结果：

print(list(range(100, 1000, 1)))

加上逻辑判断语句，最后的代码：

for m in range(100, 1000, 1):    s = (m % 10) ** 3 + int((m / 10) % 10) ** 3 + int(m / 100) ** 3    if(s == m):        print(m)

可以找到4个水仙花数：

153
370
371
407

更进一步

这个问题解决了，还可以更深入地研究一下，比如学习一下“列表推导”的语法，上面的代码可以精简为一条核心语句（不算print语句）：

flower = [m for m in range(100, 1000, 1) if (m % 10) ** 3 + int((m / 10) % 10) ** 3 + int(m / 100) ** 3 == m]print(flower)

结果不变，只不过得到的是一个列表：

[153, 370, 371, 407]

再再进一步

三位数的情况叫水仙花数，对于n位数，每个位上的数字的 n 次幂之和如果等于它本身，这种数在数学上还有一个统一的称呼：自幂数。国外对于这类数的名称是：Armstrong Number。

中国人给它们起了一些有趣的名字：

一位自幂数：独身数
两位自幂数：无
三位自幂数：水仙花数
四位自幂数：四叶玫瑰数
五位自幂数：五角星数
六位自幂数：六合数
七位自幂数：北斗七星数
八位自幂数：八仙数
九位自幂数：九九重阳数
十位自幂数：十全十美数

对于三位数，上面的程序没问题，但对于n位数，程序需要调整一下，先从4位数入手：

m = 1634print(m % 10)print(int(m / 10) % 10)print(int(m / 100) % 10)print(int(m / 1000))

再改一下，就可以适应n位数。




    
m = 1634print(int(m / 1) % 10)print(int(m / 10) % 10)print(int(m / 100) % 10)print(int(m / 1000) % 10)

再改一下，先用字符串函数求出数字的位数n，再把个位、百位、千位等等都求出来：

m = 1634n = len(str(m))digits = [(int(m / (10 ** i)) % 10) for i in range(n)]

最后的代码是这样，求出所有6位以下的自幂数需要运行几秒钟：

for m in range(1, 1000000):    n = len(str(m))    digits = [(int(m / (10 ** i)) % 10) for i in range(n)]    if(sum([d ** n for d in digits]) == m):        print(m)

结果有20个：

123456789153370371407163482089474547489272793084548834

效率再优化一些

刚才的程序找6位数的自幂数，花不了几秒钟，而要找7位数的自幂数，则需要运行几分钟，我在我的机器上跑了110秒。

我这里给出一个简单优化的思路，程序中求取n次方的数学运算比较花时间，可以把这些数值提前保存在字典中，以后直接查表即可，这样程序从110秒优化到55秒，快了一倍。

import time
start = time.time()
ten_pow = [10 ** i for i in range(20)]   # [1, 10, 100...]pow_n = [i for i in range(10)]
for m in range(1, 10000000):    n = len(str(m))    if(n == len(str(m-1)) + 1 ): # 更新字典        pow_n = [i ** n for i in range(10)]    digits = [(int(m / ten_pow[i]) % 10) for i in range(n)]    if(sum([pow_n[d] for d in digits]) == m):        print(m)
print(time.time() - start)

当然程序仍有非常大的优化空间，比如几个位置上的n次方求和时，如果已经超过数本身，则可以过滤掉一些数。代码会变得更为复杂，留着感兴趣的朋友自行探索吧。

在这个网站（https://oeis.org/A005188）有一段程序，秒求9位以下的所有自幂数，不过算法不容易看懂。

from itertools import combinations_with_replacementlist = []for k in range(1, 10):    a = [i**k for i in range(10


    
)]    for b in combinations_with_replacement(range(10), k):        x = sum(map(lambda y:a[y], b))        if x > 0 and tuple(int(d) for d in sorted(str(x))) == b:            list.append(x)list = sorted(list) # Chai Wah Wu, Aug 25 2015print(list)

运行结果：

[1, 2, 3, 4, 5, 6, 7, 8, 9, 153, 370, 371, 407, 1634, 8208, 9474, 54748, 92727, 93084, 548834, 1741725, 4210818, 9800817, 9926315, 24678050, 24678051, 88593477, 146511208, 472335975, 534494836, 912985153]

欣赏一下最大的一个，即最后一个（第88个）自幂数：

115132219018763992565095597973971522401

长达39位数。

--- END ---

近期文章：