搜寻所有水仙花数（阿姆斯特朗数）的快速算法（搜完1到60位数字组合耗时2秒）

置顶 f90boy

已于 2024-11-13 17:16:15 修改

阅读量2.4w

点赞数 18

文章标签：算法数据结构剪枝 Fortran 水仙花数最快算法递归算法

于 2023-02-10 00:00:00 首次发布

本文链接：https://blog.csdn.net/fortranboy_sh/article/details/128950838

版权

2024-07-18 更新： 对关键代码做了进一步优化。
2024-07-20 应网友建议： 对广义水仙花数进行了搜寻，见链接：
搜寻广义水仙花数（包含首位为0的情形）

搜寻水仙花数的算法很多，代码也不少。
很少有 fortran 版本的，快到 10 秒以内的算法还没见到。
给出一个搜寻所有水仙花数（阿姆斯特朗数）的快速算法，fortran90 代码。

计算耗时：
在 32 位 win7+3.2GHz 主频系统上，CVF 编译器，1～39 位 1.6 秒，1～60 位 2.9 秒。
在 64 位 Linux+3.4GHz 主频系统上，gfortran 编译器，1～39 位 1.0 秒，1～60 位 1.8 秒。

算法要点：
A、自定义数据结构 big_integer 表示大整数，构建加法和乘法运算子程序。
B、定义数组 s，预存数字 1 到 9 的 1 到 60 次方及其 0 到 60 倍数值，减少重复计算。
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 的，提前剪枝，退出后续递归；根据高位已出现数字调整递归循环的上下限。此为提高效率的关键。

附代码：

! 搜寻所有水仙花数的快速算法
! 2024-07-18

! szw_sh@163.com

! fortran 90
! CVF编译项：DF /Ox
! FPS编译项：FL32 /Ox


module abc                                            ! 公共变量

   type big_integer                                   ! 定义大整数类型
      integer*1 a(63)
      integer*1 digit
   end type

   type(big_integer) s(9,60,0:60)
   type(big_integer) tt(10),t1,t2

   integer*1 b(0:9),b1
   integer*1 nn,mm
   integer*1 nb(10)
   integer*1 a(9)

   character*60 h

end module



use abc                                               ! 主程序
integer*1 n

call cpu_time(x)                                      ! 启动计时，初始化s
call init

do nn=1,60                                            ! 循环nn=1到60，调用递归子程序搜寻

   write(*,'(2x,a,i<(nn+110)/60>)') 'n = ',nn

   call equ(tt(10),0_1)
   nb(10)=nn
   n=9
   b=0
   b1=0
   call cc(n)

   if(nn.eq.39) then                                  ! 输出39位水仙花数后，输出计时
      call cpu_time(x)
      nt=x*1000
      write(*,'(/2x,a,i<log10(nt+0.9)+1>,a/)') 'time = ',nt,' ms'
   end if

end do

call cpu_time(x)
nt=x*1000                                             ! 输出总时耗
write(*,'(/2x,a,i<log10(nt+0.9)+1>,a)') 'time = ',nt,' ms'
write(*,'(2x,a,i2)') 'total = ',mm

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,nn,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,nn,nb(n)),t2)
      if(t2%digit.lt.nn) exit
   end if

   if(n.gt.1) then

      k=s(n-1,nn,nb(n))%digit                         ! 未选数字构成nn次幂和值的最大值产生进位的位置
      if(s(n-1,nn,nb(n))%a(k)+tt(n)%a(k).ge.9) k=k+1  ! 位置k数字之和大于9，一定进位；等于9，可能进位

      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(1)) 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)                       ! 输出找到的水仙花数
      mm=mm+1                                         ! 水仙花数计数

   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
if(bb%a(ii+2).ne.0) bb%digit=ii+2                     ! +2是针对9^60*60的最大情形

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
      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*(*) 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,60
      call mul(s1,s2,i)
      s1=s2
      do k=0,60
         call mul(s1,s(i,j,k),k)
      end do
   end do
end do

end

附计算结果：

  n = 1
  9
  8
  7
  6
  5
  4
  3
  2
  1
  n = 2
  n = 3
  407
  371
  370
  153
  n = 4
  9474
  8208
  1634
  n = 5
  93084
  92727
  54748
  n = 6
  548834
  n = 7
  9926315
  9800817
  4210818
  1741725
  n = 8
  88593477
  24678051
  24678050
  n = 9
  912985153
  534494836
  472335975
  146511208
  n = 10
  4679307774
  n = 11
  94204591914
  82693916578
  49388550606
  44708635679
  42678290603
  40028394225
  32164049651
  32164049650
  n = 12
  n = 13
  n = 14
  28116440335967
  n = 15
  n = 16
  4338281769391371
  4338281769391370
  n = 17
  35875699062250035
  35641594208964132
  21897142587612075
  n = 18
  n = 19
  4929273885928088826
  4498128791164624869
  3289582984443187032
  1517841543307505039
  n = 20
  63105425988599693916
  n = 21
  449177399146038697307
  128468643043731391252
  n = 22
  n = 23
  35452590104031691935943
  28361281321319229463398
  27907865009977052567814
  27879694893054074471405
  21887696841122916288858
  n = 24
  239313664430041569350093
  188451485447897896036875
  174088005938065293023722
  n = 25
  4422095118095899619457938
  3706907995955475988644381
  3706907995955475988644380
  1553242162893771850669378
  1550475334214501539088894
  n = 26
  n = 27
  177265453171792792366489765
  174650464499531377631639254
  128851796696487777842012787
  121270696006801314328439376
  121204998563613372405438066
  n = 28
  n = 29
  23866716435523975980390369295
  19008174136254279995012734741
  19008174136254279995012734740
  14607640612971980372614873089
  n = 30
  n = 31
  2309092682616190307509695338915
  1927890457142960697580636236639
  1145037275765491025924292050346
  n = 32
  17333509997782249308725103962772
  n = 33
  186709961001538790100634132976991
  186709961001538790100634132976990
  n = 34
  1122763285329372541592822900204593
  n = 35
  12679937780272278566303885594196922
  12639369517103790328947807201478392
  n = 36
  n = 37
  1219167219625434121569735803609966019
  n = 38
  12815792078366059955099770545296129367
  n = 39
  115132219018763992565095597973971522401
  115132219018763992565095597973971522400

  time = 1630 ms
  (gfortran 1020 ms)

  n = 40
  n = 41
  n = 42
  n = 43
  n = 44
  n = 45
  n = 46
  n = 47
  n = 48
  n = 49
  n = 50
  n = 51
  n = 52
  n = 53
  n = 54
  n = 55
  n = 56
  n = 57
  n = 58
  n = 59
  n = 60

  time = 2850 ms
  (gfortran 1750 ms)
  total = 88