水仙花数（阿姆斯特朗数）并行算法，快至0.4秒

2024-07-18 更新：对关键代码做了进一步优化。

---

得空将此前写的水仙花数（阿姆斯特朗数）Fortran 代码改为并行计算，在 12核CPU（16线程）/ 3.4GHz 机器上，耗时约 0.4 秒；在 6核CPU / 3.0GHz 机器上，耗时约 0.8 秒；在 4核CPU / 3.2GHz 机器上，耗时约 1.2 秒。在目前普通配置的个人电脑上，找到全部88个水仙花数（阿姆斯特朗数），搜完 1 位到 60 位（包括大整数）范围内所有可能的正整数组合，多线程耗时减少到 1 秒。

注： nn = 35～42 时，相应长度正整数组合的搜索耗时最长可达 200 ms，加上部署线程和汇总结果所耗 100～200 ms， 400 ms 已是该算法的极限。

---

这应该是目前网络可见的个人电脑水仙花数（阿姆斯特朗数）的最快算法了。

---

将主程序做了修改，通过外部递归实现并行计算。该并行算法无需依赖专门并行环境，代码具有很好的可移植性。

除此之外，只有递归子程序 cc 中有一行修改：
write(*,‘(2x,a)’) trim(h)
修改为：
write(1,‘(2x,a)’) trim(h)
其余子程序代码均未做修改。

---

完整代码和算法注解参见：
搜寻所有水仙花数（阿姆斯特朗数）的快速算法（搜完1到60位数字组合耗时 2 秒）

---

附修改为并行计算的主程序代码：

$freeform

! 搜寻所有水仙花数的快速算法，多线程并行版本
! 2024-07-18

! szw_sh@163.com

! 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)                     ! 1到9的1到60次方及其0到60倍数值表
   type(big_integer) tt(10),t1,t2

   integer*1 b(0:9),b1                                ! 统计数字个数的临时数组
   integer*1 nn                                       ! 水仙花数的位数
   integer*1 nb(10)                                   ! 未选数字的个数
   integer*1 a(9)                                     ! 数字个数组合

   character*60 h                                     ! 用于转换和输出大整数

end module



use msflib                                            ! 主程序
use abc
parameter(ms=64)
integer*1 n
integer ck(ms)                                        ! 记录数据文件冲突测试的错误码
character*60 hh(60,10)                                ! 记录从合并文件读取的水仙花数结果
character*10 fn(ms),ff                                ! 多线程对应的数据文件名
character*200 pp                                      ! 本程序名
character*1 qq


key=nargs()                                           ! 作为是否带参数调用程序的标记
ii=getenvqq('number_of_processors',ff)                ! 获取CPU环境变量
read(ff,'(i10)') mm                                   ! 读取CPU数量，作为开启并行计算线程的数量

call getarg(0,pp)                                     ! 获取本程序名，双击启动时包含路径
if(index(pp,'\').gt.0) key_stop=1
np=len_trim(pp)+1

if(key.gt.1) then                                     ! 带参数调用本程序，启动新线程

   call getarg(1,ff)                                  ! 获取参数包含的线程序号
   read(ff(5:6),'(i2)') ifn

   open(1,file=ff,share='denyrw')                     ! 以禁止其它进程读写的方式打开数据记录文件

   call init

   do nn=ifn,60,mm                                    ! 循环nn=1到60，调用递归子程序搜寻
      write(1,'(2x,a,i<(nn+110)/60>)') 'n = ',nn
      call equ(tt(10),0_1)
      nb(10)=nn
      n=9
      call cc(n)
   end do                                             ! 每个并行线程起点不同，步长均为mm，时间基本一致

   close(1)                                           ! 关闭数据文件，作为主程序对并行线程是否结束的测试依据

else                                                  ! 主程序，不带参数启动时，进行线程分配和启动

   x=timef()                                          ! 启动计时，初始化s

   if(mm.gt.1) then
      write(*,'(/2x,a,i<mm/10+1>,a/)') 'CPU = ',mm,'，启动多线程计算 ...'
   else
      write(*,'(/2x,a/)') '单核 CPU，无法进行多线程计算 ...'
      !stop                                           ! 在linux/wine环境下，因系统缺陷会导致文件访问冲突
   end if

   fn='~sxh00.t$t'                                    ! 初始化数据文件名

   do i=1,mm                                          ! 根据序号，用start命令启动并行线程，参数realtime指定最高级CUP占用
      write(fn(i)(5:6),'(i2.2)') i
      ii=systemqq('start /b /min /realtime '//pp(1:np+1)//fn(i))
   end do                                             ! 参数b，不创建新控制台；参数min，最小化窗口，在wine下参数b不起作用

   ck=1                                               ! 初始化ck
   do while(1)
      call sleepqq(5)                                 ! 短延时，避免过多的文件打开动作影响运行效率
      do i=1,mm                                       ! 若已经成功打开过，则跳过，不在测试是否占用
         if(ck(i).ne.0) open(i,file=fn(i),iostat=ck(i))
         if(ck(i).eq.0) close(i)                      ! 若成功打开，则关闭文件
         if(ck(i).ne.0) exit                          ! 若打开失败，退出循环
      end do
      if(i.gt.mm) exit                                ! 所有并行线程的数据文件均能打开，即所有并行线程均已结束，退出测试
   end do

   ii=systemqq('copy ~sxh??.t$t out.t$t > nul')       ! 合并计算结果，删除线程的数据文件
   ii=systemqq('del ~sxh??.t$t > nul')

   open(1,file='out.t$t')

   hh=' '
   do i=1,1000                                        ! 读取合并文件，同步完成排序

      read(1,'(a)',iostat=io) h
      if(io.ne.0) exit
      if(h(1:1).ne.' ') exit                          ! 在FPS编译程序下，合并文件最后一行为一个字节，char(26)，需要忽略

      if(h(3:3).eq.'n') then
         read(h(7:8),'(i2)') n
         k=1
         hh(n,k)=h
      else
         k=k+1
         hh(n,k)=h
      end if
   end do

   close(1,status='delete')                           ! 关闭并清除合并文件

   do i=1,60                                          ! 输出水仙花数的搜索结果
      do j=1,10
         if(hh(i,j).eq.' ') exit
         write(*,'(a)') trim(hh(i,j))
      end do
   end do

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

   if(key_stop.eq.1) then                             ! 双击启动的情形，暂停等候按键退出
      write(*,'(/2x,a\)') '按任意键退出'
      qq=getcharqq()
   end if

end if

end



! 省略其它子程序代码