39位水仙花数的优化算法

Question

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

Answer 1

检查了一下大整数模块，发现效率不高。因为它需要满足所有情形的计算。而搜索水仙花数，可以通过预先保存各数字的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