很少有fortran版本的，快到一分钟内的算法还没见到。
给出一个搜寻所有水仙花数的快速算法，fortran90的代码，搜完1到60位组合耗时仅20秒。
附代码：
! 搜寻所有水仙花数的快速算法
! 2023-02-12
! 作者：szw_sh@163.com
! 算法要点：
! 一、自定义数据结构big_integer表示大整数，构建加法和乘法运算子程序。
! 二、定义数组s，预存数字1到9的1到60次方及其0到60倍数值，减少重复计算。
! 三、采用递归算法，简化代码，将递归层次n之外的所有递归参数都定义为公共变量，提高效率；
! 递归枚举数字9到1的个数组合，总和不超过nn位；0的个数不需要递归求解， 递归结束时自
! 然确定；n为1时，检验是否为水仙花数。
! 四、提前剪枝策略：
! 1、递归从9开始，到1结束；每一层的个数循环从大到小；剩余可选数量记为nb；
! 2、检查各位数字幂的和值tt的位数，大于nn则cycle；小于nn，且增加剩余可选数最大幂和
! 值仍不到nn位，则退出本层递归；
! 3、n大于等于2时，根据后续nb个数字s的最大值，确定tt不受影响的高位部分，统计其中已
! 选数字和未选数字的个数；已选数字个数大于选定个数的， 或者未选数字个数大于nb及
! 受影响位数之和的，提前剪枝，退出后续递归；此为提高效率的关键。
! 计算耗时：39位约10秒，全部60位约20秒
Module BigInteger
implicit none
private
Type , public :: big_integer
Integer(1) a(63)
Integer(1) digit
contains
Procedure :: equ , mul , add , toString
Generic :: assignment(=) => equ
Generic :: operator(+) => add
Generic :: operator(*) => mul
End Type big_integer
contains
Subroutine equ(this, n) ! 大整数赋值子程序，aa=n，n<10
class(big_integer) , intent(OUT) :: this
Integer(1) , intent(IN) :: n
this%a = 0
this%digit = 1
this%a(1) = n
End Subroutine equ
Type(big_integer) Function mul(aa, n) result(bb) ! 大整数乘法子程序，bb=aa*n，n<100
class(big_integer) , intent(IN) :: aa
Integer(1) , intent(IN) :: n
Integer(2) k, k1, kn
Integer(1) i, ii
If (n==0 .Or. aa%digit==1 .And. aa%a(1)==0) Then ! n=0 或者 aa=0
bb = 0_1
Return
Else If (n==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_2)
k1 = k/10
End Do
bb%digit = ii ! 确定位数
If (bb%a(ii+1)/=0) bb%digit = ii + 1
If (bb%a(ii+2)/=0) bb%digit = ii + 2 ! +2是针对9^60*60的最大情形
End Function mul
Type(big_integer) Function add(aa, bb) result(cc) ! 大整数加法子程序，cc=aa+bb
class(big_integer) , intent(IN) :: aa, bb
Integer(1) k, k1, i, ii
If (aa%digit==1 .And. aa%a(1)==0) Then ! aa=0
cc = bb
Return
End If
If (bb%digit==1 .And. bb%a(1)==0) Then ! bb=0
cc = aa
Return
End If
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>=10) Then
cc%a(i) = k - 10
k1 = 1
End If
End Do
cc%digit = ii ! 确定位数
If (cc%a(ii+1)>0) cc%digit = ii + 1
End Function add
Character(len=60) Function toString(aa) result(hh) ! 大整数转字符串子程序，hh=aa
class(big_integer) , intent(IN) :: aa ! hh长度按照7*9=63设定
Integer(1) ii, j
hh = ' '
ii = aa%digit
Do j = 1, ii
Write (hh(60-j+1:60-j+1), '(i1.1)') aa%a(j)
End Do
hh(1:60-aa%digit) = ' ' ! 删除左端的0
hh = adjustl(hh) ! 左对齐
End Function toString
End Module BigInteger
Module abc ! 公共变量
use BigInteger
implicit none
Type (big_integer) s(9, 60, 0:60)
Type (big_integer) tt(10), t1, t2
Integer(1) b(9)
Integer(1) nn
Integer(1) nb(10)
Integer(1) a(9)
Character(len=60) h
contains
Recursive Subroutine cc(n) ! 查找水仙花数的递归子程序
Integer(1) n, i, j, k, ii
Do i = nb(n+1), 0, -1 ! 循环+递归，构造9到1的数字个数组合
a(n) = i
nb(n) = nb(n+1) - i
tt(n) = tt(n+1) + s(n,nn,i) ! 计算数字nn次幂的和，只计算变化的位
If (tt(n)%digit>nn) Then ! 位数超出nn，循环递减a(n)
Cycle
Else If (tt(n)%digit<nn) Then ! 位数不达nn，根据后续nb预测仍不达nn的，退出递归
if(n>1) t2 = tt(n) + s(n-1,nn,nb(n))
If (t2%digit<nn) Exit
End If
If (n>1) Then
k = s(n-1, nn, nb(n))%digit + 1 ! 未选数字构成nn次幂和值的最大值可能产生的进位位置
Do ii = k, nn ! 向高位检查是否存在连续9，避免出现进位漏查情形
If (tt(n)%a(ii)/=9) Exit ! 经测试，可能有数百万个
End Do
b = 0
Do j = ii + 1, nn ! 统计tt高位不变化位置的数字
k = tt(n)%a(j)
If (k==0) Cycle
b(k) = b(k) + 1
If (k>=n) Then ! 已选数字，若tt中的大于已选，提前退出统计
If (b(k)>a(k)) Exit
Else ! 未选数字，若tt中的大于剩余未选+低位未统计位数
If (b(k)>nb(n)+ii) Exit ! 按照未统计低位每位都是未选某一数字的最不利情形
End If
End Do
If (j<=nn) Cycle ! 判定tt存在已选、未选数字的矛盾的情形，提前剪枝
Call cc(n-1_1) ! 继续递归调用
Else ! 递归结束
b = 0 ! 统计tt各数字的数量
Do j = 1, nn
k = tt(n)%a(j)
If (k==0) Cycle ! 0不统计
b(k) = b(k) + 1
If (b(k)>a(k)) Exit ! 若超出已选的数量，矛盾，提前跳出
End Do
If (j<=nn) Cycle ! 若有矛盾，跳过
Do j = 9, 1, -1 ! 检验各数字的数量，不符合的，提前跳出
If (a(j)/=b(j)) Exit
End Do
If (j>=1) Cycle ! 数字组合不相同的，跳过
h = tt(n)%toString() ! tt转为字符串
If (h=='0') Cycle ! 结果为0的，不计入水仙花数
Write (*, '(2x,a)') trim(h) ! 输出找到的水仙花数
End If
End Do
End Subroutine cc
Subroutine init()
Type (big_integer) s1, s2
Integer(1) i, j, k
Do i = 1, 9
s1 = 1_1
Do j = 1, 60
s2 = s1 * i
s1 = s2
Do k = 0, 60
s(i,j,k) = s1 * k
End Do
End Do
End Do
End Subroutine init
End Module abc
Program fcode_cn
Use abc ! 主程序
use ifport
implicit none
Integer(1) n
integer :: nt
real(8) :: x
x = timef() ! 启动计时，初始化s
Call init()
Do nn = 1, 60 ! 循环nn=1到60，调用递归子程序搜寻
Write (*, '(2x,"n = ",g0)') nn
tt(10) = 0_1
nb(10) = nn
n = 9
Call cc(n)
If (nn==39) Then ! 输出39位水仙花数后，输出计时
nt = timef()*1000
Write (*, '(/2x,"time = ",g0," ms"/)') nt
End If
End Do
nt = timef()*1000 ! 输出总时耗
Write (*, '(/2x,"time = ",g0," ms")') nt
End Program fcode_cn
