10毫秒找出16位以内43个水仙花数

代码为 Fortran90 格式。

之前发过大整数的完整版本，有网友希望能给出 Java 或者 Python 代码。遗憾的是，借助 AI 也没能完整转换。于是把简化的 Fortran 代码发上来。

使用 integer*8 （64 位整数类型），最大可以计算到 16 位水仙花数。代码很短，不到 90 行。在 win7 / CPU 3.2GHz （睿频 3.6GHz） 配置的电脑上，运行耗时不到 10 毫秒。

具体算法思路可以参见：
搜寻所有水仙花数（阿姆斯特朗数）的快速算法（搜完1到60位数字组合耗时2秒）

---

程序代码：

! 2024-07-18
! szw_sh@163.com

module abc
integer*8 s(9,16,0:16),tt(10),t1,t2,smin(16),smax(16)
integer*1 b(0:9),b1,nn,nb(10),a(9),sd(9,16,0:16),mm
end module

use abc
call cpu_time(x)
do i=1,9
   do j=1,16
      do k=0,16
         s(i,j,k)=int8(i)**j*k
         sd(i,j,k)=log10(max(s(i,j,k)*1d0,1d0))+1
      end do
   end do
end do
do i=1,16
   smax(i)=int8(10)**i-1
   smin(i)=int8(10)**(i-1)
end do
do nn=1,16
   write(*,'(2x,a,g0)') 'n = ',nn
   tt(10)=0
   nb(10)=nn
   b=0 
   b1=0
   call cc(9_1)
end do
call cpu_time(x)
write(*,'(/2x,a,g0,a)') 'time = ',int(x*1000),' ms'
write(*,'(2x,a,g0)') 'total = ',mm
end

recursive subroutine cc(n)
use abc
integer*1 i,j,k,n,ii
do i=nb(n+1)-b1+b(n),b(n),-1
   a(n)=i
   nb(n)=nb(n+1)-i
   tt(n)=tt(n+1)+s(n,nn,i)
   if(tt(n).gt.smax(nn)) then
      cycle
   else if(tt(n).lt.smin(nn)) then
      t2=tt(n)+s(n-1,nn,nb(n))
      if(t2.lt.smin(nn)) exit
   end if
   b=0
   b1=0
   if(n.gt.1) then
      kk=sd(n-1,nn,nb(n))+1
      t1=tt(n)/smin(kk)
      do ii=kk,nn
         k=mod(t1,10)
         t1=t1/10
         if(k.ne.9) exit
      end do
      t2=tt(n)/smin(min(ii+1,nn))
      ! 遇连续9进位时，避免smin越界
      do j=ii+1,nn
         k=mod(t2,10)
         t2=t2/10
         b(k)=b(k)+1
         if(k.ge.n) 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 cc(n-1_1)
   else
      t1=tt(n)
      do j=1,nn
         k=mod(t1,10)
         t1=t1/10
         b(k)=b(k)+1
         if(k.ge.n) 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
      write(*,'(2x,g0)') tt(n)
      mm=mm+1
   end if
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

  time = 9 ms
  total = 43

---

附图：www.onlinegdb.com 在线编译运行的截图