大数求阶乘Fortran实现初探
本帖最后由 Jackdaw 于 2018-2-5 09:05 编辑丁酉年癸丑月丙寅日，遇群友讨论大数阶乘，遂编制代码实现之，与朋友们分享（我不会说我是为了装X的:-lol）。
此文分三部分内容，源码+Makefile+测试结果
测试电脑：Ubuntu 16.04.3 LTS64-bitCPU: i5-4210M@2.6GHz
源码：
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
Makefile：
#!/bin/bash
FC = gfortran
FF = -static -O2
SC = factorial.f90
XX = x.out
$(XX) :
$(FC) $(FF) $(SC) -o $(XX)
clean:
rm -rf *.mod *.txt $(XX)
测试结果：
分别为10、100、 300、1000的阶乘
10000的阶乘（0.224秒）
ps：屏幕太小放不下了
牛逼牛逼
{:5_119:} lz为啥不考虑用FFT变换加速整数相乘呢？
http://numbers.computation.free.fr/Constants/Algorithms/representation.html
http://numbers.computation.free.fr/Constants/Algorithms/fft.html pasuka 发表于 2018-2-24 11:06
lz为啥不考虑用FFT变换加速整数相乘呢？
http://numbers.computation.free.fr/Constants/Algorithms/repres ...
直接算的，没考虑什么其他的算法，我试试这个:-lol 牛叉{:2_26:}{:2_26:} 你为何这么吊
页:
[1]