2026-04-08 更新更新由 Fortran 转换的 Python 代码。2024-07-20 应网友建议对广义水仙花数进行了搜寻见链接搜寻广义水仙花数包含首位为0的情形2024-07-18 更新对关键代码做了进一步优化。搜寻水仙花数的算法很多代码也不少。很少有 Fortran 版本的快到 10 秒以内的算法还没见到。给出一个搜寻所有水仙花数阿姆斯特朗数的快速算法Fortran90 代码。计算耗时在 32 位 win73.2GHz 主频系统上CVF 编译器139 位1.6 秒160 位2.9 秒。在 64 位 Linux3.4GHz 主频系统上gfortran 编译器139 位1.0 秒160 位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 结束每一层的个数循环从大到小剩余可选数量记为 nbb、检查各位数字幂的和值 tt 的位数大于 nn 则 cycle小于 nn且增加剩余可选数最大幂和值仍不到 nn 位则退出本层递归c、n 大于等于 2 时根据后续 nb 个数字 s 的最大值确定 tt 不受影响的高位部分统计其中已选数字和未选数字的个数已选数字个数大于选定个数的 或者 0 到 (n-1) 未选数字个数大于 nb 的提前剪枝退出后续递归根据高位已出现数字调整递归循环的上下限。此为提高效率的关键。附 Fortran 代码!搜寻所有水仙花数的快速算法!2024-07-18!szw_sh163.com!fortran90!CVF编译项DF /Ox!FPS编译项FL32 /Ox module abc!公共变量typebig_integer!定义大整数类型 integer*1 a(63)integer*1 digit endtypetype(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 initdonn1,60!循环nn1到60调用递归子程序搜寻 write(*,(2x,a,i(nn110)/60))n ,nn call equ(tt(10),0_1)nb(10)nnn9b0b10call cc(n)if(nn.eq.39)then!输出39位水仙花数后输出计时 call cpu_time(x)ntx*1000 write(*,(/2x,a,ilog10(nt0.9)1,a/))time ,nt, msendifenddocall cpu_time(x)ntx*1000!输出总时耗 write(*,(/2x,a,ilog10(nt0.9)1,a))time ,nt, mswrite(*,(2x,a,i2))total ,mm end recursive subroutine cc(n)!查找水仙花数的递归子程序 use abc integer*1 n,i,j,k,iidoinb(n1)-b1b(n),b(n),-1!循环递归构造9到1的数字个数组合 a(n)i nb(n)nb(n1)-i call add(tt(n1),s(n,nn,i),tt(n))!计算数字nn次幂的和只计算变化的位 if(tt(n)%digit.gt.nn)then!位数超出nn循环递减a(n)cycleelseif(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)exitendifif(n.gt.1)thenks(n-1,nn,nb(n))%digit!未选数字构成nn次幂和值的最大值产生进位的位置 if(s(n-1,nn,nb(n))%a(k)tt(n)%a(k).ge.9)kk1!位置k数字之和大于9一定进位等于9可能进位doiik,nn!向高位检查是否存在连续9避免出现进位漏查情形 if(tt(n)%a(ii).ne.9)exit!经测试有数百万个 enddob0b10dojii1,nn!统计tt高位不变化位置的数字ktt(n)%a(j)b(k)b(k)1 if(k.ge.n)then!已选数字若tt中的大于已选提前退出统计 if(b(k).gt.a(k))exitelse!0到(n-1)的未选数字若tt中的大于剩余未选个数的提前退出b1b11 if(b1.gt.nb(n))exitendifenddoif(j.le.nn)cycle!判定tt存在已选、未选数字的矛盾的情形提前剪枝 call cc(n-1_1)!继续递归调用else!递归结束b0!统计tt各数字的数量b10doj1,nnktt(n)%a(j)b(k)b(k)1 if(k.gt.0)thenif(b(k).gt.a(k))exit!若超出已选的数量矛盾提前跳出elseb1b11 if(b1.gt.nb(1))exitendifenddoif(j.le.nn)cycle!若有矛盾跳过 call chr(tt(n),h)!tt转为字符串 if(h.eq.0)cycle!结果为0非正整数不计入水仙花数 write(*,(2x,a))trim(h)!输出找到的水仙花数mmmm1!水仙花数计数 endifenddoend subroutine equ(aa,n)!大整数赋值子程序aann10use abc type(big_integer)aa integer*1 n aa%a0aa%digit1aa%a(1)n end subroutine mul(aa,bb,n)!大整数乘法子程序bbaa*nn100use 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!n0或者aa0call equ(bb,0_1)returnelseif(n.eq.1)then!n1bbaareturnendifbb%a0!初始化bb进位数k1k10iiaa%digitknndoi1,ii3!分段乘法进位kaa%a(i)*knk1 bb%a(i)mod(k,10)k1k/10 enddobb%digitii!确定位数 if(bb%a(ii1).ne.0)bb%digitii1 if(bb%a(ii2).ne.0)bb%digitii2!2是针对9^60*60的最大情形 end subroutine add(aa,bb,cc)!大整数加法子程序ccaabb use abc type(big_integer)aa,bb,cc integer*1 k,k1,i,ii cc%a0!初始化cc进位数k1k10iimax(aa%digit,bb%digit)doi1,ii1!分段加法进位kaa%a(i)bb%a(i)k1 cc%a(i)kk10if(k.ge.10)thencc%a(i)k-10k11endifenddocc%digitii!确定位数 if(cc%a(ii1).gt.0)cc%digitii1 end subroutine chr(aa,hh)!大整数转字符串子程序hhaa use abc character*(*)hh type(big_integer)aa integer*1 jhh doj1,nn hh(nn-j1:nn-j1)char(48aa%a(j))enddoend subroutine init use abc!初始化数组s type(big_integer)s1,s2 integer*1 i,j,kdoi1,9call equ(s1,1_1)doj1,60call mul(s1,s2,i)s1s2dok0,60call mul(s1,s(i,j,k),k)enddoenddoenddoend附 Python 代码# 搜寻所有水仙花数的快速算法# 2025-06-12# szw_sh163.comimportmathimporttimeclassABC:def__init__(self):存储计算水仙花数所需的全局数据# s[i][n][k]存储k个数字i的n次幂之和self.s[[[0]*61for_inrange(61)]for_inrange(10)]# 存储中间计算结果self.tt[0]*11# 记录n位数的上下界self.smin[0]*61# 10^(n-1)self.smax[0]*61# (10^n)-1# 存储数字出现次数统计self.b[0]*10self.b10# 未选中数字计数self.nn0# 当前搜索的位数# 存储各数字个数分配self.nb[0]*11# 数字i的数量self.a[0]*10# 临时使用# 存储幂和结果的位数self.sd[[[0]*61for_inrange(61)]for_inrange(10)]self.mm0# 水仙花数计数器# 创建全局模块实例modABC()definit():初始化数组s和sd预计算数字幂和及位数foriinrange(1,10):# 数字1-9forjinrange(1,61):# 位数1-60basei**j# 数字i的j次幂forkinrange(0,61):# 数字i出现的次数0-60mod.s[i][j][k]base*k# 计算并存储结果的位数ifmod.s[i][j][k]0:mod.sd[i][j][k]math.floor(math.log10(mod.s[i][j][k]))1else:mod.sd[i][j][k]0# 初始化n位数的范围限制foriinrange(1,61):# 位数1-60mod.smin[i]10**(i-1)# 最小n位数mod.smax[i]10**i-1# 最大n位数defcc(n):递归搜索水仙花数的核心函数 使用分治剪枝策略避免无效搜索# 计算循环范围逆序减少无效计算startmod.nb[n1]-mod.b1mod.b[n]endmod.b[n]step-1foriinrange(start,end-1,step):# 更新数字n的出现次数和总数限制mod.a[n]i mod.nb[n]mod.nb[n1]-i# 累加n次幂和仅计算变化部分mod.tt[n]mod.tt[n1]mod.s[n][mod.nn][i]# 重要剪枝1超出当前位数的范围ifmod.tt[n]mod.smax[mod.nn]:continueelifmod.tt[n]mod.smin[mod.nn]:# 检查后续可能的最大值仍不足t2_valmod.tt[n]mod.s[n-1][mod.nn][mod.nb[n]]ift2_valmod.smin[mod.nn]:break# 重置数字频率统计mod.b[0]*10mod.b10ifn1:# 递归情形# 计算可能的进位位置kkmod.sd[n-1][mod.nn][mod.nb[n]]1divisormod.smin[kk]ifkkmod.nnelsemod.smin[mod.nn]t1_valmod.tt[n]//divisorifdivisor0elsemod.tt[n]# 重要优化查找连续9的位置避免漏掉进位iikk# 默认从进位位置开始foridxinrange(kk,mod.nn1):ift1_val0:# 已处理完所有数字iiidxbreakk_digitt1_val%10# 取当前位数字t1_val//10ifk_digit!9:# 关键优化点iiidx# 记录第一个非9的位置break# 结尾处理ifidxmod.nn:iiidx1# 重要剪枝2检查高位数字分布是否合法divisor2mod.smin[min(ii1,mod.nn)]t2_valmod.tt[n]//divisor2ifdivisor20elsemod.tt[n]validTrueforjinrange(ii1,mod.nn1):ift2_val0:# 高位已处理完毕breakk_digitt2_val%10t2_val//10mod.b[k_digit]1# 判断数字k是否已被选中ifk_digitn:# 已选中数字ifmod.b[k_digit]mod.a[k_digit]:validFalsebreakelse:# 未选中数字mod.b11ifmod.b1mod.nb[n]:validFalsebreakifnotvalid:continue# 存在冲突跳过当前组合# 递归处理下一个数字cc(n-1)else:# 基本情况n1# 最终校验数字分布是否匹配t1_valmod.tt[1]validTrueforjinrange(1,mod.nn1):k_digitt1_val%10t1_val//10mod.b[k_digit]1# 类似上述的分布检查ifk_digitn:ifmod.b[k_digit]mod.a[k_digit]:validFalsebreakelse:mod.b11ifmod.b1mod.nb[n]:validFalsebreak# 找到有效水仙花数ifvalidandjmod.nn:# jnn表示完整校验print(f{mod.tt[1]})# 输出结果mod.mm1# 增加计数# 主程序if__name____main__:start_timetime.time()init()# 初始化预计算数据# 搜索1到60位水仙花数fornninrange(1,61):print(f n {nn})mod.nnnn# 设置当前位数mod.tt[10]0# 初始化结果缓存mod.nb[10]nn# 初始化总数字数mod.b[0]*10# 重置统计mod.b10# 从数字9开始递归搜索cc(9)# 39位时输出中间计时ifnn39:elapsed_msint((time.time()-start_time)*1000)print(f\n time {elapsed_ms}ms\n)# 输出最终计时和总数total_elapsed_msint((time.time()-start_time)*1000)print(f\n time {total_elapsed_ms}ms)print(f total {mod.mm})附计算结果n1987654321n2n3407371370153n4947482081634n5930849272754748n6548834n79926315980081742108181741725n8885934772467805124678050n9912985153534494836472335975146511208n104679307774n119420459191482693916578493885506064470863567942678290603400283942253216404965132164049650n12n13n1428116440335967n15n1643382817693913714338281769391370n17358756990622500353564159420896413221897142587612075n18n194929273885928088826449812879116462486932895829844431870321517841543307505039n2063105425988599693916n21449177399146038697307128468643043731391252n22n233545259010403169193594328361281321319229463398279078650099770525678142787969489305407447140521887696841122916288858n24239313664430041569350093188451485447897896036875174088005938065293023722n2544220951180958996194579383706907995955475988644381370690799595547598864438015532421628937718506693781550475334214501539088894n26n27177265453171792792366489765174650464499531377631639254128851796696487777842012787121270696006801314328439376121204998563613372405438066n28n2923866716435523975980390369295190081741362542799950127347411900817413625427999501273474014607640612971980372614873089n30n31230909268261619030750969533891519278904571429606975806362366391145037275765491025924292050346n3217333509997782249308725103962772n33186709961001538790100634132976991186709961001538790100634132976990n341122763285329372541592822900204593n351267993778027227856630388559419692212639369517103790328947807201478392n36n371219167219625434121569735803609966019n3812815792078366059955099770545296129367n39115132219018763992565095597973971522401115132219018763992565095597973971522400time1630ms(gfortran1020ms)(Python30122ms)n40n41n42n43n44n45n46n47n48n49n50n51n52n53n54n55n56n57n58n59n60time2850ms(gfortran1750ms)(Python51376ms)total88