网上有很多讨论 21 位水仙花数算法 的旧帖子，大概这曾经是一个竞赛题目。所见到的各种语言的代码，运行耗时大多在 10 秒到数分钟，也有个别是在 10 秒以内的，高位水仙花数算法还有不少优化空间。
网上找不到 21位水仙花数的 fortran 版本，可能这门语言太古老了！
本文的这个算法是 fortran90 版的，将运行耗时缩短到了 20 毫秒。
有详细注解，供参考。
如果是 39 位，运行耗时为 150 毫秒。
测试环境：Win7，CPU主频3.2GHz / 3.6GHz(睿频)。
算法要点：
A、自定义数据结构 big_integer 表示大整数，构建加法和乘法运算子程序。
B、定义数组 s，预存数字 1 到 9 的 nn 次方及其 0 到 nn 倍数值，减少重复计算。
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 的，提前剪枝，退出后续递归；根据高位已出现数字调整递归循环的上下限。此为提高效率的关键。
代码以 21 位为例，如需计算 39 位水仙花数，只需修改常量 nn = 39 。
$freeform
! 搜寻nn位水仙花数的快速算法
! 2024-07-18
! szw_sh@163.com
! CVF编译项：DF /Ox
module abc ! 公共变量
integer*1 nn ! 水仙花数的位数
parameter(nn=21)
type big_integer ! 定义大整数类型
integer*1 a(nn+3) ! 考虑进位计算，取nn+3位
integer*1 digit ! a 各位数字，a(1)是最低位；digit 位数
end type
type(big_integer) s(9,0:nn) ! 1到9的nn次方及其0到nn倍数值表
type(big_integer) tt(10),t1,t2
integer*1 b(0:9),b1 ! 统计数字个数的临时数组
integer*1 nb(10) ! 未选数字的个数
integer*1 a(9) ! 数字个数组合
character*(nn) h ! 用于转换和输出大整数
end module
use abc ! 主程序
integer*1 n
call cpu_time(x) ! 启动计时，初始化s
call init
write(*,'(/2x,a,i2)') 'n = ',nn
call equ(tt(10),0_1)
nb(10)=nn
n=9
b=0
b1=0
call cc(n)
call cpu_time(x) ! 输出总时耗
nt=x*1000
write(*,'(/2x,a,g0,a)') 'time = ',nt,' ms'
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,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,nb(n)),t2)
if(t2%digit.lt.nn) exit
end if
if(n.gt.1) then
k=s(n-1,nb(n))%digit+1 ! 未选数字构成nn次幂和值的最大值可能产生的进位位置
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(n)) 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) ! 输出找到的水仙花数
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
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 ! 因为最大进位为1，不用mod，直接用减法
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*(nn) 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,nn
call mul(s1,s2,i)
s1=s2
if(j.lt.nn) cycle
do k=0,nn
call mul(s1,s(i,k),k)
end do
end do
end do
end
21位的运行结果：
n = 21
449177399146038697307
128468643043731391252
time = 18 ms
39位的运行结果：
n = 39
115132219018763992565095597973971522401
115132219018763992565095597973971522400
time = 152 ms
软木酒塞 + 牛皮纸袋
配图无关代码
683
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
