很少有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