编辑代码

module abc
   type big_integer
      integer*1 a(63)
      integer*1 digit
   end type
   type(big_integer) s(9,60,0:60)
   type(big_integer) tt(10),t1,t2
   integer*1 b(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)
call init
do nn=1,39
   write(*,'(2x,a,g0)') 'n = ',nn
   call equ(tt(10),0_1)
   nb(10)=nn
   n=9
   call cc(n)
end do
call cpu_time(x)
nt=x*1000
write(*,'(/2x,a,g0,a)') 'time = ',nt,' ms'
write(*,'(2x,a,g0)') 'total = ',mm
end
recursive subroutine cc(n)
use abc
integer*1 n,i,j,k,ii
do i=nb(n+1),0,-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) 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)
         if(k.ge.n) then
            b(k)=b(k)+1
            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)
         if(k.gt.0) then
            b(k)=b(k)+1
            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(h.eq.'0') 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
if(aa%digit.eq.1.and.aa%a(1).eq.0) then
   cc=bb
   return
end if
if(bb%digit.eq.1.and.bb%a(1).eq.0) then
   cc=aa
   return
end if
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*60 hh
type(big_integer) aa
integer*1 ii,j
hh=' '
ii=aa%digit
do j=1,ii
   hh(ii-j+1:ii-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,60
      call mul(s1,s2,i)
      s1=s2
      do k=0,60
         call mul(s1,s(i,j,k),k)
      end do
   end do
end do
end