首页 >

搜寻所有水仙花数的快速算法

作者:fortranboy  日期:02-19
来源:Fcode研讨团队
搜寻水仙花数的算法很多,代码也不少。
很少有fortran版本的,快到一分钟内的算法还没见到。
给出一个搜寻所有水仙花数的快速算法,fortran90的代码,搜完1到60位组合耗时仅20秒。

附代码:
! 搜寻所有水仙花数的快速算法
! 2023-02-12

! 作者:szw_sh@163.com

! 算法要点:

! 一、自定义数据结构big_integer表示大整数,构建加法和乘法运算子程序。
! 二、定义数组s,预存数字1到9的1到60次方及其0到60倍数值,减少重复计算。
! 三、采用递归算法,简化代码,将递归层次n之外的所有递归参数都定义为公共变量,提高效率;
!     递归枚举数字9到1的个数组合,总和不超过nn位;0的个数不需要递归求解, 递归结束时自
!     然确定;n为1时,检验是否为水仙花数。
! 四、提前剪枝策略:
!     1、递归从9开始,到1结束;每一层的个数循环从大到小;剩余可选数量记为nb;
!     2、检查各位数字幂的和值tt的位数,大于nn则cycle;小于nn,且增加剩余可选数最大幂和
!        值仍不到nn位,则退出本层递归;
!     3、n大于等于2时,根据后续nb个数字s的最大值,确定tt不受影响的高位部分,统计其中已
!        选数字和未选数字的个数;已选数字个数大于选定个数的, 或者未选数字个数大于nb及
!        受影响位数之和的,提前剪枝,退出后续递归;此为提高效率的关键。

! 计算耗时:39位约10秒,全部60位约20秒

Module BigInteger
  implicit none
  private
  Type , public :: big_integer
    Integer(1) a(63)
    Integer(1) digit
  contains
    Procedure :: equ , mul , add , toString
    Generic :: assignment(=) => equ
    Generic :: operator(+) => add
    Generic :: operator(*) => mul
  End Type big_integer
  
contains

  Subroutine equ(this, n) ! 大整数赋值子程序,aa=n,n<10
    class(big_integer) , intent(OUT) :: this
    Integer(1) , intent(IN) :: n
    this%a = 0
    this%digit = 1
    this%a(1) = n
  End Subroutine equ

  Type(big_integer) Function mul(aa, n) result(bb) ! 大整数乘法子程序,bb=aa*n,n<100
    class(big_integer) , intent(IN) :: aa
    Integer(1) , intent(IN) :: n
    Integer(2) k, k1, kn
    Integer(1) i, ii
    If (n==0 .Or. aa%digit==1 .And. aa%a(1)==0) Then ! n=0 或者 aa=0
      bb = 0_1
      Return
    Else If (n==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_2)
      k1 = k/10
    End Do
    bb%digit = ii ! 确定位数
    If (bb%a(ii+1)/=0) bb%digit = ii + 1
    If (bb%a(ii+2)/=0) bb%digit = ii + 2 ! +2是针对9^60*60的最大情形
  End Function mul

  Type(big_integer) Function add(aa, bb) result(cc) ! 大整数加法子程序,cc=aa+bb
    class(big_integer) , intent(IN) :: aa, bb
    Integer(1) k, k1, i, ii
    If (aa%digit==1 .And. aa%a(1)==0) Then ! aa=0
      cc = bb
      Return
    End If
    If (bb%digit==1 .And. bb%a(1)==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>=10) Then
        cc%a(i) = k - 10
        k1 = 1
      End If
    End Do
    cc%digit = ii ! 确定位数
    If (cc%a(ii+1)>0) cc%digit = ii + 1
  End Function add

  Character(len=60) Function toString(aa) result(hh) ! 大整数转字符串子程序,hh=aa
    class(big_integer) , intent(IN) :: aa ! hh长度按照7*9=63设定
    Integer(1) ii, j
    hh = ' '
    ii = aa%digit
    Do j = 1, ii
      Write (hh(60-j+1:60-j+1), '(i1.1)') aa%a(j)
    End Do
    hh(1:60-aa%digit) = ' ' ! 删除左端的0
    hh = adjustl(hh) ! 左对齐
  End Function toString
  
End Module BigInteger
  
Module abc ! 公共变量
  use BigInteger
  implicit none
  Type (big_integer) s(9, 60, 0:60)
  Type (big_integer) tt(10), t1, t2

  Integer(1) b(9)
  Integer(1) nn
  Integer(1) nb(10)
  Integer(1) a(9)
  Character(len=60) h
  
contains

  Recursive Subroutine cc(n) ! 查找水仙花数的递归子程序
    Integer(1) n, i, j, k, ii
    Do i = nb(n+1), 0, -1 ! 循环+递归,构造9到1的数字个数组合
      a(n) = i
      nb(n) = nb(n+1) - i
      tt(n) = tt(n+1) + s(n,nn,i) ! 计算数字nn次幂的和,只计算变化的位

      If (tt(n)%digit>nn) Then ! 位数超出nn,循环递减a(n)
        Cycle
      Else If (tt(n)%digit<nn) Then ! 位数不达nn,根据后续nb预测仍不达nn的,退出递归
        if(n>1) t2 = tt(n) + s(n-1,nn,nb(n))
        If (t2%digit<nn) Exit
      End If

      If (n>1) Then
        k = s(n-1, nn, nb(n))%digit + 1 ! 未选数字构成nn次幂和值的最大值可能产生的进位位置
        Do ii = k, nn ! 向高位检查是否存在连续9,避免出现进位漏查情形
          If (tt(n)%a(ii)/=9) Exit ! 经测试,可能有数百万个
        End Do
        b = 0
        Do j = ii + 1, nn ! 统计tt高位不变化位置的数字
          k = tt(n)%a(j)
          If (k==0) Cycle
          b(k) = b(k) + 1
          If (k>=n) Then ! 已选数字,若tt中的大于已选,提前退出统计
            If (b(k)>a(k)) Exit
          Else ! 未选数字,若tt中的大于剩余未选+低位未统计位数
            If (b(k)>nb(n)+ii) Exit ! 按照未统计低位每位都是未选某一数字的最不利情形
          End If
        End Do
        If (j<=nn) Cycle ! 判定tt存在已选、未选数字的矛盾的情形,提前剪枝
        Call cc(n-1_1) ! 继续递归调用
      Else ! 递归结束
        b = 0 ! 统计tt各数字的数量
        Do j = 1, nn
          k = tt(n)%a(j)
          If (k==0) Cycle ! 0不统计
          b(k) = b(k) + 1
          If (b(k)>a(k)) Exit ! 若超出已选的数量,矛盾,提前跳出
        End Do
        If (j<=nn) Cycle ! 若有矛盾,跳过
        Do j = 9, 1, -1 ! 检验各数字的数量,不符合的,提前跳出
          If (a(j)/=b(j)) Exit
        End Do
        If (j>=1) Cycle ! 数字组合不相同的,跳过
        h = tt(n)%toString() ! tt转为字符串
        If (h=='0') Cycle ! 结果为0的,不计入水仙花数
        Write (*, '(2x,a)') trim(h) ! 输出找到的水仙花数
      End If
    End Do
  End Subroutine cc
  
  Subroutine init()
    Type (big_integer) s1, s2
    Integer(1) i, j, k
    Do i = 1, 9
      s1 = 1_1
      Do j = 1, 60
        s2 = s1 * i
        s1 = s2
        Do k = 0, 60
          s(i,j,k) = s1 * k
        End Do
      End Do
    End Do
  End Subroutine init  
  
End Module abc

Program fcode_cn
  Use abc ! 主程序
  use ifport
  implicit none
  Integer(1) n
  integer :: nt
  real(8) :: x
  x = timef() ! 启动计时,初始化s
  Call init()
  Do nn = 1, 60 ! 循环nn=1到60,调用递归子程序搜寻
    Write (*, '(2x,"n = ",g0)') nn
    tt(10) = 0_1
    nb(10) = nn
    n = 9
    Call cc(n)
    If (nn==39) Then ! 输出39位水仙花数后,输出计时
      nt = timef()*1000
      Write (*, '(/2x,"time = ",g0," ms"/)') nt
    End If
  End Do
  nt = timef()*1000 ! 输出总时耗
  Write (*, '(/2x,"time = ",g0," ms")') nt
End Program fcode_cn
常规|工具|专业|读物|
代码|教学|算法|
首页 >
FortranCoder手机版-导航