搜寻广义水仙花数（包含首位为0的情形）

刚刚更新了搜寻水仙花数的代码：
单线程代码 …… 2 秒
多线程代码 …… 0.4 秒

网友 wsktuuytyh 建议搜寻一下首位为 0 的情形。于是，把首位为 0 的一起纳入水仙花数，姑且把它们统称为广义水仙花数（广义阿姆斯特朗数）吧。

---

确定位数的上限：
9⁸⁵×85 ≈ 1.097×10⁸³，最大可达 84 位；
9⁸⁶×86 ≈ 9.985×10⁸³，达不到 85 位；
因此，超过 85 位不存在广义水仙花数。

调整代码：
搜寻位数，1～85；
递归退出两个条件，大于 nn 位，保持不变；小于 nn 位，改成小于 (nn-1) 位；
输出时，首位为 0 的，添加 0 一并输出。

单线程代码见文末。

---

计算结果：
广义水仙花数共 159 个（包括严格定义的水仙花数 88 个）。
最大的为 64 位（首位 0，+ 63 位正整数）：0129729706383218800940922971779520807379706064324417989786624799。

计算耗时：
Linux 在线系统 + 3.6GHz 主频 + gfortran 编译器，单线程，95 秒；
win7 + 3.6GHz 主频 + CVF 编译器，4 线程并行，36 秒；
win10 + intel core i9 29.1GHz 主频 (实测等效) + CVF 编译器，20 线程并行，12 秒。

注：并行拆分任务中，最小单位耗时最长的可能达 4～8 秒（n = 55～65），总体效能提升很有限。

  n = 1
  9
  8
  7
  6
  5
  4
  3
  2
  1
  n = 2
  01
  n = 3
  407
  371
  370
  153
  n = 4
  9474
  8208
  1634
  n = 5
  93084
  92727
  54748
  04151
  04150
  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
  0564240140138
  n = 14
  28116440335967
  n = 15
  n = 16
  4338281769391371
  4338281769391370
  n = 17
  35875699062250035
  35641594208964132
  21897142587612075
  n = 18
  n = 19
  4929273885928088826
  4498128791164624869
  3289582984443187032
  1517841543307505039
  n = 20
  63105425988599693916
  n = 21
  449177399146038697307
  128468643043731391252
  n = 22
  n = 23
  35452590104031691935943
  28361281321319229463398
  27907865009977052567814
  27879694893054074471405
  21887696841122916288858
  n = 24
  239313664430041569350093
  188451485447897896036875
  174088005938065293023722
  n = 25
  4422095118095899619457938
  3706907995955475988644381
  3706907995955475988644380
  1553242162893771850669378
  1550475334214501539088894
  0832662335985815242605071
  0832662335985815242605070
  0114735624485461118832514
  n = 26
  n = 27
  177265453171792792366489765
  174650464499531377631639254
  128851796696487777842012787
  121270696006801314328439376
  121204998563613372405438066
  077888878776432530886487094
  n = 28
  n = 29
  23866716435523975980390369295
  19008174136254279995012734741
  19008174136254279995012734740
  14607640612971980372614873089
  05022908050052864745436221003
  04716716265341543230394614213
  n = 30
  n = 31
  2309092682616190307509695338915
  1927890457142960697580636236639
  1145037275765491025924292050346
  0793545620525277858657607629822
  n = 32
  17333509997782249308725103962772
  n = 33
  186709961001538790100634132976991
  186709961001538790100634132976990
  032186410459473623435614002227248
  n = 34
  1122763285329372541592822900204593
  n = 35
  12679937780272278566303885594196922
  12639369517103790328947807201478392
  07673249664848285722449710136138169
  05250083909873201044638631458484846
  n = 36
  091097771122214850683543503173498149
  n = 37
  1219167219625434121569735803609966019
  0618670315011216949642642321868915308
  0418510620208153136884574959404115822
  n = 38
  12815792078366059955099770545296129367
  07403697806790834730831423191927508283
  07320233109580046612992702336326619665
  n = 39
  115132219018763992565095597973971522401
  115132219018763992565095597973971522400
  083281830613691836766959173718984508549
  083281823928125880164896079973522049472
  016427762135335641330720936105651700735
  n = 40
  n = 41
  n = 42
  036428594490313158783584452532870892261556
  024202117350726607855379216981761437616685
  n = 43
  0865428295996489702819913290126290251086057
  0649929577038606483022256247961798702282298
  0541494183189741758024440515698931971170565
  0540133062379646837091926683197423034975342
  0431693315336650106841122671394059061346259
  0216196771535553659064432337763021725280303
  0110488001357745108347376522372211155920487
  n = 44
  06810209536021751861114918348460992955509943
  05856845793349592707882266933841781895916377
  04870692736452280126562816349691961381966274
  03906378741950401404540845367463938860869507
  n = 45
  052629118994026288113690621458128419465866933
  052455000909452690385913874112931520316143509
  043988588628637966718625825202934800362393972
  017804804411551074034314974935268388104718848
  017674136814275371165896659566101835068144683
  017500125632575295611741402478151097427104361
  017500125632575295611741402478151097427104360
  n = 46
  0552652296628071929876997582254879810411842988
  0394500610649009608022685635958104085968134641
  0394500610649009608022685635958104085968134640
  n = 47
  04959921468176469608573534993418650725564922822
  03551561298137983224493440104432142398068825798
  02844606761993812649569688176852543331764572068
  02134869831132206002876377009473775323572880694
  n = 48
  n = 49
  0344121003015798953461109472522636869502957985402
  n = 50
  n = 51
  055683645994157961689595909141165132491193429749749
  027887539002613689375558346990932796323618371540363
  n = 52
  n = 53
  03009336689622424785094591221547283824126694123199598
  n = 54
  027051156903992754495077743355694541141634293209536961
  027051156903992754495077743355694541141634293209536960
  010185128017405712763334458116098808034687490694771371
  010185128017405712763334458116098808034687490694771370
  n = 55
  0152490314755876303649798312407542072870728583784707998
  0152349829250452484527824009552907045495378563176436265
  n = 56
  02193762240761908392137860899658607674401938496187046968
  01644854105813017601539897250795638996670193545173003361
  01644854105813017601539897250795638996670193545173003360
  n = 57
  014799195103184806214749814246033579630671509034462600766
  n = 58
  n = 59
  01199348308270467958191508402511086560689741786191274137702
  01199348017964116555212449342061478333395100875934086835818
  n = 60
  016177694250271739806081953963047779450935812520923227956595
  014386811900482199401154446588646400893796400907674489355497
  n = 61
  0145594628412421292850519939876219279596702430926364267115105
  0129495102611459020639899678988620521527176986805703581773808
  n = 62
  n = 63
  011793324696943658158328433258762969319709773412597747795160016
  n = 64
  0129729706383218800940922971779520807379706064324417989786624799
  n = 65
  n = 66
  n = 67
  n = 68
  n = 69
  n = 70
  n = 71
  n = 72
  n = 73
  n = 74
  n = 75
  n = 76
  n = 77
  n = 78
  n = 79
  n = 80
  n = 81
  n = 82
  n = 83
  n = 84
  n = 85

  time = 94571 ms
  (4线程并行，time = 36320 ms)
  (20线程并行，time = 11708 ms)
  total = 159

---

Fortran 代码：（恰好是 159 行）

module abc
   type big_integer
      integer*1 a(88),digit
   end type
   type(big_integer) s(9,85,0:85),tt(10),t1,t2
   integer*1 a(9),b(0:9),b1,nn,nb(10)
   integer mm
   data mm/0/
   character*85 h
end module
use abc
integer*1 n
call cpu_time(x)
call init
do nn=1,85
   write(*,'(2x,a,i0)') 'n = ',nn
   call equ(tt(10),0_1)
   nb(10)=nn
   n=9
   b=0
   b1=0
   call cc(n)
end do
call cpu_time(x)
nt=x*1000
write(*,'(/2x,a,i0,a)') 'time = ',nt,' ms'
write(*,'(2x,a,i0)') 'total = ',mm
end
recursive subroutine cc(n)
use abc
integer*1 n,i,j,k,ii
do i=nb(n+1)-b1+b(n),b(n),-1
   a(n)=i
   nb(n)=nb(n+1)-i
   call add(tt(n+1),s(n,nn,i),tt(n))
   if(tt(n)%digit.gt.nn) then
      cycle
   else if(tt(n)%digit.lt.nn) then
      call add(tt(n),s(n-1,nn,nb(n)),t2)
      if(t2%digit.lt.nn-1) exit
   end if
   if(n.gt.1) then
      k=s(n-1,nn,nb(n))%digit
      if(s(n-1,nn,nb(n))%a(k)+tt(n)%a(k).ge.9) k=k+1
      do ii=k,nn
         if(tt(n)%a(ii).ne.9) exit
      end do
      b=0
      b1=0
      do j=ii+1,nn
         k=tt(n)%a(j)
         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
      b=0
      b1=0
      do j=1,nn
         k=tt(n)%a(j)
         b(k)=b(k)+1
         if(k.gt.0) 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 chr(tt(n),h)
      if(tt(n)%a(1)+tt(n)%digit.eq.1) cycle
      write(*,'(2x,a)') trim(h)
      mm=mm+1
   end if
end do
end
subroutine equ(aa,n)
use abc
type(big_integer) aa
integer*1 n
aa%a=0
aa%digit=1
aa%a(1)=n
end
subroutine mul(aa,bb,n)
use 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
   call equ(bb,0_1)
   return
else if(n.eq.1) then
   bb=aa
   return
end if
bb%a=0
k1=0
ii=aa%digit
kn=n
do i=1,ii+3
   k=aa%a(i)*kn+k1
   bb%a(i)=mod(k,10)
   k1=k/10
end do
bb%digit=ii
if(bb%a(ii+1).ne.0) bb%digit=ii+1
if(bb%a(ii+2).ne.0) bb%digit=ii+2
end
subroutine add(aa,bb,cc)
use abc
type(big_integer) aa,bb,cc
integer*1 k,k1,i,ii
cc%a=0
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.ge.10) then
      cc%a(i)=k-10
      k1=1
   end if
end do
cc%digit=ii
if(cc%a(ii+1).gt.0) cc%digit=ii+1
end
subroutine chr(aa,hh)
use abc
character*(*) hh
type(big_integer) aa
integer*1 j
hh=' '
do j=1,nn
   hh(nn-j+1:nn-j+1)=char(48+aa%a(j))
end do
end
subroutine init
use abc
type(big_integer) s1,s2
integer*1 i,j,k
do i=1,9
   call equ(s1,1_1)
   do j=1,85
      call mul(s1,s2,i)
      s1=s2
      do k=0,85
         call mul(s1,s(i,j,k),k)
      end do
   end do
end do
end