39位水仙花数,即39位正整数的各位数字39次方之和等于它本身。
利用fortran语言编程,用递归枚举39个数字构造组合的算法,找到两个39位水仙花数,耗时约400秒。
这两个数是:
115132219018763992565095597973971522401
115132219018763992565095597973971522400
求:更为高效率、更快速的算法。
recursive subroutine cc(a,tt,n) ! 查找水仙花数的递归子程序
use abc
type(big_integer) tt,t1
integer*1 a(60),n
np=9 ! 逆序查找,比较有效率
if(n.gt.1) np=a(n-1)
do i=np,0,-1 ! 循环+递归,构造nn位的数字组合
a(n)=i
t1=tt+s(i,nn,1) ! 计算数字nn次幂的和,只计算变化的位
if(t1.gt.smax(nn)) cycle ! 超出本级构造的最大数,递减a(n)
if(t1+s(i,nn,nn-n).lt.smin(nn)) exit ! 小于nn位数下限值,退出本级(循环)递归
if(n.eq.nn) then ! 符合上下限值的情形
s1=t1
b=a
do k1=0,54,9 ! 拆解出各位数字,与a进行比较
s2=s1/1000000000
s3=s1-s2*1000000000
s1=s2
do k2=1,9
j=k1+k2
if(j.gt.n) exit
s4=s3/10
k=s3-s4*10
s3=s4
do jj=n-j+1,1,-1 ! a中有匹配的,将该位置代换掉,下次匹配减少一个单元
if(k.eq.b(jj)) then
b(jj)=b(n-j+1)
exit
end if
end do
if(jj.ge.1) then ! 找到匹配,继续;无匹配的,跳出
cycle
else
exit
end if
end do
if(jj.lt.1.and.j.le.n) exit ! 两条件标记存在不匹配的数字,提前跳出
end do
if(j.le.n) cycle ! 无匹配,跳过
write(*,'(2x,a)') trim(char(t1)) ! 输出找到的水仙花数
else
call cc(a,t1,n+1) ! 继续递归调用
end if
end do
end