[Fortran] 纯文本查看 复制代码
module factorial_mod
implicit none
private
integer ,parameter :: id = 4
integer(id) :: num_zeros
integer(id),allocatable :: fact(:)
integer(id) :: n = 12000
integer(id),parameter :: n_digit = 5
integer(id),parameter :: critical = 10**n_digit
integer :: up
public id
public factorial
contains
subroutine init() ! 初始化
implicit none
allocate( fact(n) )
fact = 0
fact( 1 ) = 1
num_zeros = 0
up = 1
end subroutine init
subroutine finalize() ! 结束
implicit none
deallocate(fact)
end subroutine finalize
integer(id) function digit( a ) ! 几位数
implicit none
integer(id),value :: a
digit = 1
do while ( a .ne. 0 )
a = a / critical
digit = digit + 1
end do
end function digit
subroutine carry_forward() ! 进位
implicit none
integer(id) :: i,j
integer(id) :: tmp
do j = n , 1 , -1
if( fact( j ) .ne. 0 ) exit
end do
up = j
do i = 1, up
if( fact(i).ge.critical ) then
tmp = mod(fact(i),critical)
fact(i+1) = fact(i+1) + fact(i)/critical
fact(i) = tmp
if(i .eq. up) up = up+1
end if
end do
end subroutine carry_forward
subroutine carry_back() ! 退位
implicit none
integer(id) :: i
do i = 1, up
if( fact(i) .ne. 0) exit
end do
num_zeros = num_zeros + (i - 1) * n_digit
fact(1:up-i+1) = fact(i:up)
fact(up-i+2:) = 0
up = up -i +1
end subroutine carry_back
subroutine display() ! 打印结果
implicit none
integer(id) :: i
do i = up, 1, -1
write(*,"(i5.5)",advance="no") fact(i)
if(mod(up-i+1,35) .eq. 0) write(*,*)
end do
write(*,"(a,g0)") " E+",num_zeros
end subroutine display
subroutine factorial( num ) ! 阶乘
implicit none
integer(id),intent(in),value :: num
integer(id) :: i,j,k,npos
integer(id),allocatable :: mod_i(:)
integer(id),allocatable :: tmp_fact(:)
integer(id) :: tmp
call init()
if( .not. allocated(tmp_fact) ) allocate( tmp_fact(n) )
do i = 2, num
npos = digit(i)
if( .not. allocated(mod_i) )allocate( mod_i(npos) )
tmp = i
do j = 1, npos
mod_i(j) = mod(tmp,critical)
tmp = tmp / critical
end do
if(npos .ge. critical ) then ! 非 n_digit 位数
tmp_fact = 0
do j = 1, npos
if(mod_i(j) .eq. 0) cycle
k = up
forall (k = 1:up ) tmp_fact(k+j-1) = tmp_fact(k+j-1) + fact(k) * mod_i(j)
do k = 1, up
if( tmp_fact( k ) .ge. critical) then
tmp = mod(tmp_fact(k),critical)
tmp_fact(k+1) = tmp_fact(k+1) + tmp_fact(k)/critical
tmp_fact(k) = tmp
end if
end do
end do
fact = tmp_fact
call carry_forward()
call carry_back()
else ! n_digit 位数
forall( k = 1:up ) fact(k) = fact(k) * i
call carry_forward()
call carry_back()
end if
if(npos.ne. digit(i+1) ) deallocate(mod_i)
end do
call carry_forward()
call carry_back()
call display()
call finalize()
if( allocated(tmp_fact) ) deallocate( tmp_fact )
end subroutine factorial
end module factorial_mod
program main
use factorial_mod
implicit none
integer(id) :: t(3)
integer(id) :: num
write(*,*) "Please input a positive integer number(not greater than 13500):"
read(*,*) num
call system_clock(t(1),t(3))
call factorial(num)
call system_clock(t(2),t(3))
write(*,"(a,f7.3,a)") " Time :",(t(2)-t(1))*1.0/t(3) ," s"
end program main
）。
