2025-03-13 更新： 新增由 Fortran 转换的 Python 代码。
2024-07-20 应网友建议： 对广义水仙花数进行了搜寻，见链接：
搜寻广义水仙花数（包含首位为0的情形）
2024-07-18 更新： 对关键代码做了进一步优化。
搜寻水仙花数的算法很多，代码也不少。
很少有 Fortran 版本的，快到 10 秒以内的算法还没见到。
给出一个搜寻所有水仙花数（阿姆斯特朗数）的快速算法，Fortran90 代码。
计算耗时：
在 32 位 win7+3.2GHz 主频系统上，CVF 编译器，1～39 位 1.6 秒，1～60 位 2.9 秒。
在 64 位 Linux+3.4GHz 主频系统上，gfortran 编译器，1～39 位 1.0 秒，1～60 位 1.8 秒。
算法要点：
A、自定义数据结构 big_integer 表示大整数，构建加法和乘法运算子程序。
B、定义数组 s，预存数字 1 到 9 的 1 到 60 次方及其 0 到 60 倍数值，减少重复计算。
C、采用递归算法，简化代码，将递归层次 n 之外的所有递归参数都定义为公共变量，提高效率；递归枚举数字 9 到 1 的个数组合，总和不超过 nn 位；0 的个数不需要递归求解， 递归结束时自然确定；n 为 1 时，检验是否为水仙花数。
D、提前剪枝策略：
a、递归从 9 开始，到 1 结束；每一层的个数循环从大到小；剩余可选数量记为 nb；
b、检查各位数字幂的和值 tt 的位数，大于 nn 则 cycle；小于 nn，且增加剩余可选数最大幂和值仍不到 nn 位，则退出本层递归；
c、n 大于等于 2 时，根据后续 nb 个数字 s 的最大值，确定 tt 不受影响的高位部分，统计其中已选数字和未选数字的个数；已选数字个数大于选定个数的， 或者 0 到 (n-1) 未选数字个数大于 nb 的，提前剪枝，退出后续递归；根据高位已出现数字调整递归循环的上下限。此为提高效率的关键。
附 Fortran 代码：
! 搜寻所有水仙花数的快速算法
! 2024-07-18
! szw_sh@163.com
! fortran 90
! CVF编译项：DF /Ox
! FPS编译项：FL32 /Ox
module abc ! 公共变量
type big_integer ! 定义大整数类型
integer*1 a(63)
integer*1 digit
end type
type(big_integer) s(9,60,0:60)
type(big_integer) tt(10),t1,t2
integer*1 b(0:9),b1
integer*1 nn,mm
integer*1 nb(10)
integer*1 a(9)
character*60 h
end module
use abc ! 主程序
integer*1 n
call cpu_time(x) ! 启动计时，初始化s
call init
do nn=1,60 ! 循环nn=1到60，调用递归子程序搜寻
write(*,'(2x,a,i<(nn+110)/60>)') 'n = ',nn
call equ(tt(10),0_1)
nb(10)=nn
n=9
b=0
b1=0
call cc(n)
if(nn.eq.39) then ! 输出39位水仙花数后，输出计时
call cpu_time(x)
nt=x*1000
write(*,'(/2x,a,i<log10(nt+0.9)+1>,a/)') 'time = ',nt,' ms'
end if
end do
call cpu_time(x)
nt=x*1000 ! 输出总时耗
write(*,'(/2x,a,i<log10(nt+0.9)+1>,a)') 'time = ',nt,' ms'
write(*,'(2x,a,i2)') 'total = ',mm
end
recursive subroutine cc(n) ! 查找水仙花数的递归子程序
use abc
integer*1 n,i,j,k,ii
do i=nb(n+1)-b1+b(n),b(n),-1 ! 循环+递归，构造9到1的数字个数组合
a(n)=i
nb(n)=nb(n+1)-i
call add(tt(n+1),s(n,nn,i),tt(n)) ! 计算数字nn次幂的和，只计算变化的位
if(tt(n)%digit.gt.nn) then ! 位数超出nn，循环递减a(n)
cycle
else if(tt(n)%digit.lt.nn) then ! 位数不达nn，根据后续nb预测仍不达nn的，退出递归
call add(tt(n),s(n-1,nn,nb(n)),t2)
if(t2%digit.lt.nn) exit
end if
if(n.gt.1) then
k=s(n-1,nn,nb(n))%digit ! 未选数字构成nn次幂和值的最大值产生进位的位置
if(s(n-1,nn,nb(n))%a(k)+tt(n)%a(k).ge.9) k=k+1 ! 位置k数字之和大于9，一定进位；等于9，可能进位
do ii=k,nn ! 向高位检查是否存在连续9，避免出现进位漏查情形
if(tt(n)%a(ii).ne.9) exit ! 经测试，有数百万个
end do
b=0
b1=0
do j=ii+1,nn ! 统计tt高位不变化位置的数字
k=tt(n)%a(j)
b(k)=b(k)+1
if(k.ge.n) then ! 已选数字，若tt中的大于已选，提前退出统计
if(b(k).gt.a(k)) exit
else ! 0到(n-1)的未选数字，若tt中的大于剩余未选个数的，提前退出
b1=b1+1
if(b1.gt.nb(n)) exit
end if
end do
if(j.le.nn) cycle ! 判定tt存在已选、未选数字的矛盾的情形，提前剪枝
call cc(n-1_1) ! 继续递归调用
else ! 递归结束
b=0 ! 统计tt各数字的数量
b1=0
do j=1,nn
k=tt(n)%a(j)
b(k)=b(k)+1
if(k.gt.0) then
if(b(k).gt.a(k)) exit ! 若超出已选的数量，矛盾，提前跳出
else
b1=b1+1
if(b1.gt.nb(1)) exit
end if
end do
if(j.le.nn) cycle ! 若有矛盾，跳过
call chr(tt(n),h) ! tt转为字符串
if(h.eq.'0') cycle ! 结果为0，非正整数，不计入水仙花数
write(*,'(2x,a)') trim(h) ! 输出找到的水仙花数
mm=mm+1 ! 水仙花数计数
end if
end do
end
subroutine equ(aa,n) ! 大整数赋值子程序，aa=n，n<10
use abc
type(big_integer) aa
integer*1 n
aa%a=0
aa%digit=1
aa%a(1)=n
end
subroutine mul(aa,bb,n) ! 大整数乘法子程序，bb=aa*n，n<100
use abc
type(big_integer) aa,bb
integer*2 k,k1,kn
integer*1 i,ii,n
if(n.eq.0.or.aa%digit.eq.1.and.aa%a(1).eq.0) then ! n=0 或者 aa=0
call equ(bb,0_1)
return
else if(n.eq.1) then ! n=1
bb=aa
return
end if
bb%a=0 ! 初始化bb，进位数k1
k1=0
ii=aa%digit
kn=n
do i=1,ii+3 ! 分段乘法，进位
k=aa%a(i)*kn+k1
bb%a(i)=mod(k,10)
k1=k/10
end do
bb%digit=ii ! 确定位数
if(bb%a(ii+1).ne.0) bb%digit=ii+1
if(bb%a(ii+2).ne.0) bb%digit=ii+2 ! +2是针对9^60*60的最大情形
end
subroutine add(aa,bb,cc) ! 大整数加法子程序，cc=aa+bb
use abc
type(big_integer) aa,bb,cc
integer*1 k,k1,i,ii
cc%a=0 ! 初始化cc，进位数k1
k1=0
ii=max(aa%digit,bb%digit)
do i=1,ii+1 ! 分段加法，进位
k=aa%a(i)+bb%a(i)+k1
cc%a(i)=k
k1=0
if(k.ge.10) then
cc%a(i)=k-10
k1=1
end if
end do
cc%digit=ii ! 确定位数
if(cc%a(ii+1).gt.0) cc%digit=ii+1
end
subroutine chr(aa,hh) ! 大整数转字符串子程序，hh=aa
use abc
character*(*) hh
type(big_integer) aa
integer*1 j
hh=' '
do j=1,nn
hh(nn-j+1:nn-j+1)=char(48+aa%a(j))
end do
end
subroutine init
use abc ! 初始化数组s
type(big_integer) s1,s2
integer*1 i,j,k
do i=1,9
call equ(s1,1_1)
do j=1,60
call mul(s1,s2,i)
s1=s2
do k=0,60
call mul(s1,s(i,j,k),k)
end do
end do
end do
end
附 Python 代码：
import time
import copy
# claude3.7 转换自 Fortran
# 定义大整数类型
class BigInteger:
def __init__(self):
self.a = [0] * 63 # 存储各位数字，索引从0开始但对应原Fortran代码从1开始
self.digit = 0 # 位数
# 全局变量定义
s = None # 将在init函数中初始化
tt = [BigInteger() for _ in range(11)]
t1 = BigInteger()
t2 = BigInteger()
b = [0] * 10 # 0-9数字的计数
b1 = 0
nn = 0 # 当前搜索的位数
mm = 0 # 找到的水仙花数的总数
nb = [0] * 11
a = [0] * 10
h = ""
def equ(aa, n):
"""大整数赋值函数，aa=n，n<10"""
aa.a = [0] * 63
aa.digit = 1
aa.a[0] = n
def mul(aa, bb, n):
"""大整数乘法函数，bb=aa*n，n<100"""
if n == 0 or (aa.digit == 1 and aa.a[0] == 0): # n=0 或者 aa=0
equ(bb, 0)
return
elif n == 1: # n=1
bb.a = aa.a.copy()
bb.digit = aa.digit
return
bb.a = [0] * 63 # 初始化bb
k1 = 0 # 进位数
ii = aa.digit
for i in range(ii + 3): # 分段乘法，进位
k = aa.a[i] * n + k1
bb.a[i] = k % 10
k1 = k // 10
bb.digit = ii # 确定位数
if bb.a[ii] != 0:
bb.digit = ii + 1
if bb.a[ii + 1] != 0:
bb.digit = ii + 2 # +2针对9^60*60的最大情形
def add(aa, bb, cc):
"""大整数加法函数，cc=aa+bb"""
cc.a = [0] * 63 # 初始化cc
k1 = 0 # 进位数
ii = max(aa.digit, bb.digit)
for i in range(ii + 1): # 分段加法，进位
k = aa.a[i] + bb.a[i] + k1
cc.a[i] = k
k1 = 0
if k >= 10:
cc.a[i] = k - 10
k1 = 1
cc.digit = ii # 确定位数
if cc.a[ii] > 0:
cc.digit = ii + 1
def chr_to_str(aa):
"""大整数转字符串函数"""
result = ""
for j in range(nn, 0, -1):
result += chr(48 + aa.a[j-1])
return result
def init():
"""初始化数组s"""
global s
# 创建三维数组 s[1-9][1-60][0-60]
s = [[[BigInteger() for k in range(61)] for j in range(61)] for i in range(10)]
s1 = BigInteger()
s2 = BigInteger()
for i in range(1, 10):
equ(s1, 1)
for j in range(1, 61):
mul(s1, s2, i)
s1 = copy.deepcopy(s2)
for k in range(61):
mul(s1, s[i][j][k], k)
def cc(n):
"""查找水仙花数的递归子程序"""
global a, b, b1, tt, nn, mm, nb, s, t2
for i in range(nb[n+1] - b1 + b[n], b[n] - 1, -1): # 循环+递归，构造9到1的数字个数组合
a[n] = i
nb[n] = nb[n+1] - i
add(tt[n+1], s[n][nn][i], tt[n]) # 计算数字nn次幂的和，只计算变化的位
if tt[n].digit > nn: # 位数超出nn，循环递减a(n)
continue
elif tt[n].digit < nn: # 位数不达nn，根据后续nb预测仍不达nn的，退出递归
add(tt[n], s[n-1][nn][nb[n]], t2)
if t2.digit < nn:
break
if n > 1:
k = s[n-1][nn][nb[n]].digit # 未选数字构成nn次幂和值的最大值产生进位的位置
if s[n-1][nn][nb[n]].a[k-1] + tt[n].a[k-1] >= 9: # 位置k数字之和大于9，一定进位；等于9，可能进位
k += 1
ii = k
while ii <= nn and tt[n].a[ii-1] == 9: # 向高位检查是否存在连续9，避免出现进位漏查情形
ii += 1
b = [0] * 10
b1 = 0
j = ii + 1
while j <= nn: # 统计tt高位不变化位置的数字
digit_val = tt[n].a[j-1]
b[digit_val] += 1
if digit_val >= n: # 已选数字，若tt中的大于已选，提前退出统计
if b[digit_val] > a[digit_val]:
break
else: # 0到(n-1)的未选数字，若tt中的大于剩余未选个数的，提前退出
b1 += 1
if b1 > nb[n]:
break
j += 1
if j <= nn: # 判定tt存在已选、未选数字的矛盾的情形，提前剪枝
continue
cc(n-1) # 继续递归调用
else: # 递归结束
b = [0] * 10 # 统计tt各数字的数量
b1 = 0
j = 1
while j <= nn:
digit_val = tt[n].a[j-1]
b[digit_val] += 1
if digit_val > 0:
if b[digit_val] > a[digit_val]: # 若超出已选的数量，矛盾，提前跳出
break
else:
b1 += 1
if b1 > nb[1]:
break
j += 1
if j <= nn: # 若有矛盾，跳过
continue
result = chr_to_str(tt[n]) # tt转为字符串
if result == '0': # 结果为0，非正整数，不计入水仙花数
continue
print(f" {result}") # 输出找到的水仙花数
mm += 1 # 水仙花数计数
# 主程序
def main():
global nn, mm, nb, tt, b, b1
start_time = time.time()
init() # 初始化s
for nn_val in range(1, 61): # 循环nn=1到60，调用递归子程序搜寻
nn = nn_val
print(f" n = {nn}")
equ(tt[10], 0)
nb[10] = nn
n = 9
b = [0] * 10
b1 = 0
cc(n)
if nn == 39: # 输出39位水仙花数后，输出计时
elapsed = (time.time() - start_time) * 1000
print(f"\n time = {int(elapsed)} ms\n")
elapsed = (time.time() - start_time) * 1000 # 输出总时耗
print(f"\n time = {int(elapsed)} ms")
print(f" total = {mm}")
if __name__ == "__main__":
main()
附计算结果：
n = 1
9
8
7
6
5
4
3
2
1
n = 2
n = 3
407
371
370
153
n = 4
9474
8208
1634
n = 5
93084
92727
54748
n = 6
548834
n = 7
9926315
9800817
4210818
1741725
n = 8
88593477
24678051
24678050
n = 9
912985153
534494836
472335975
146511208
n = 10
4679307774
n = 11
94204591914
82693916578
49388550606
44708635679
42678290603
40028394225
32164049651
32164049650
n = 12
n = 13
n = 14
28116440335967
n = 15
n = 16
4338281769391371
4338281769391370
n = 17
35875699062250035
35641594208964132
21897142587612075
n = 18
n = 19
4929273885928088826
4498128791164624869
3289582984443187032
1517841543307505039
n = 20
63105425988599693916
n = 21
449177399146038697307
128468643043731391252
n = 22
n = 23
35452590104031691935943
28361281321319229463398
27907865009977052567814
27879694893054074471405
21887696841122916288858
n = 24
239313664430041569350093
188451485447897896036875
174088005938065293023722
n = 25
4422095118095899619457938
3706907995955475988644381
3706907995955475988644380
1553242162893771850669378
1550475334214501539088894
n = 26
n = 27
177265453171792792366489765
174650464499531377631639254
128851796696487777842012787
121270696006801314328439376
121204998563613372405438066
n = 28
n = 29
23866716435523975980390369295
19008174136254279995012734741
19008174136254279995012734740
14607640612971980372614873089
n = 30
n = 31
2309092682616190307509695338915
1927890457142960697580636236639
1145037275765491025924292050346
n = 32
17333509997782249308725103962772
n = 33
186709961001538790100634132976991
186709961001538790100634132976990
n = 34
1122763285329372541592822900204593
n = 35
12679937780272278566303885594196922
12639369517103790328947807201478392
n = 36
n = 37
1219167219625434121569735803609966019
n = 38
12815792078366059955099770545296129367
n = 39
115132219018763992565095597973971522401
115132219018763992565095597973971522400
time = 1630 ms
(gfortran 1020 ms)
(Python 30122 ms)
n = 40
n = 41
n = 42
n = 43
n = 44
n = 45
n = 46
n = 47
n = 48
n = 49
n = 50
n = 51
n = 52
n = 53
n = 54
n = 55
n = 56
n = 57
n = 58
n = 59
n = 60
time = 2850 ms
(gfortran 1750 ms)
(Python 51376 ms)
total = 88
扫一扫
787
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
