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
检查了一下大整数模块，发现效率不高。因为它需要满足所有情形的计算。而搜索水仙花数，可以通过预先保存各数字的1到60次幂，搜索过程只需要用到加法。这几天抽空改写了大整数子程序，效率大大提高，搜索39位水仙花数只需要100秒左右。
module abc ! 公共变量
type big_integer ! 定义大整数类型
integer1 a(63) ! 因为60连续进位最大结果是599，进位为59，所以取60+2位
integer*1 digit ! digit 位数
end type
type(big_integer) s(0:9,60,0:60) ! 0到9的1到60次方及其0到60倍数值表
integer1 b(60) ! 各位数字的临时数组
integer4 nn
character*60 h ! 用于转换和输出大整数
end module
subroutine add(aa,bb,cc) ! 大整数加法子程序，cc=aa+bb
use abc
type(big_integer) aa,bb,cc
integer*1 k,k1
if(aa%digit.eq.1.and.aa%a(1).eq.0) then ! aa=0
cc=bb
return
end if
if(bb%digit.eq.1.and.bb%a(1).eq.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.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
115132219018763992565095597973971522401