大数求阶乘Fortran实现初探

Jackdaw · 发表于 2018-2-4 20:21:44

本帖最后由 Jackdaw 于 2018-2-5 09:05 编辑

丁酉年癸丑月丙寅日，遇群友讨论大数阶乘，遂编制代码实现之，与朋友们分享（我不会说我是为了装X的

）。
此文分三部分内容，源码+Makefile+测试结果
测试电脑：Ubuntu 16.04.3 LTS 64-bit CPU: i5-4210M@2.6GHz

源码：

[Fortran] 纯文本查看 复制代码

001module factorial_mod
002  implicit none 
003  private
004  integer    ,parameter   :: id = 4
005  integer(id)             :: num_zeros 
006  integer(id),allocatable :: fact(:)
007  integer(id)             :: n = 12000
008  integer(id),parameter   :: n_digit  = 5
009  integer(id),parameter   :: critical = 10**n_digit
010  integer                 :: up
011   
012  public id
013  public factorial
014   
015contains
016   
017  subroutine init() ! 初始化
018    implicit none 
019    allocate( fact(n) )
020    fact = 0
021    fact( 1 ) = 1
022    num_zeros = 0
023    up = 1 
024  end subroutine init 
025 
026  subroutine finalize() ! 结束
027    implicit none
028    deallocate(fact)
029  end subroutine finalize 
030 
031  integer(id) function digit( a ) ! 几位数
032    implicit none
033    integer(id),value :: a 
034    digit = 1
035    do while ( a .ne. 0 )
036      a = a / critical
037      digit = digit + 1
038    end do
039  end function digit
040   
041  subroutine carry_forward() ! 进位
042    implicit none
043    integer(id)  :: i,j
044    integer(id)  :: tmp
045    do j = n , 1 , -1
046      if( fact( j ) .ne. 0 ) exit
047    end do 
048    up = j
049    do i = 1, up 
050      if( fact(i).ge.critical ) then 
051        tmp       = mod(fact(i),critical)
052        fact(i+1) = fact(i+1) + fact(i)/critical
053        fact(i)   = tmp 
054        if(i .eq. up) up = up+1
055      end if
056    end do  
057  end subroutine carry_forward
058    
059  subroutine carry_back() ! 退位
060    implicit none
061    integer(id) :: i   
062    do i = 1, up
063      if( fact(i) .ne. 0) exit
064    end do 
065    num_zeros = num_zeros + (i - 1) * n_digit
066    fact(1:up-i+1) = fact(i:up)
067    fact(up-i+2:)  = 0 
068    up = up -i +1
069  end subroutine carry_back
070 
071  subroutine display()  ! 打印结果
072    implicit none
073    integer(id) :: i   
074    do i = up, 1, -1
075      write(*,"(i5.5)",advance="no") fact(i)
076      if(mod(up-i+1,35) .eq. 0) write(*,*)
077    end do 
078    write(*,"(a,g0)") " E+",num_zeros 
079  end subroutine display
080   
081   
082  subroutine factorial( num ) ! 阶乘
083    implicit none
084    integer(id),intent(in),value :: num 
085    integer(id)                  :: i,j,k,npos
086    integer(id),allocatable      :: mod_i(:)
087    integer(id),allocatable      :: tmp_fact(:)
088    integer(id)                  :: tmp 
089 
090    call init()
091    if( .not. allocated(tmp_fact) ) allocate( tmp_fact(n) )
092     
093    do i = 2, num  
094      npos = digit(i)
095      if( .not. allocated(mod_i) )allocate( mod_i(npos) )
096      tmp = i
097      do j = 1, npos 
098        mod_i(j) = mod(tmp,critical)
099        tmp = tmp / critical
100      end do
101      if(npos .ge. critical ) then ! 非 n_digit 位数
102        tmp_fact = 0
103        do j = 1, npos 
104          if(mod_i(j) .eq. 0) cycle
105          k = up  
106          forall (k = 1:up ) tmp_fact(k+j-1) = tmp_fact(k+j-1) + fact(k) * mod_i(j)
107          do k = 1, up
108            if( tmp_fact( k ) .ge. critical) then 
109              tmp           = mod(tmp_fact(k),critical)
110              tmp_fact(k+1) = tmp_fact(k+1) + tmp_fact(k)/critical
111              tmp_fact(k)   = tmp
112            end if
113          end do  
114        end do
115        fact = tmp_fact 
116        call carry_forward()       
117        call carry_back()
118      else ! n_digit 位数 
119        forall( k = 1:up ) fact(k) = fact(k) * i 
120        call carry_forward()        
121        call carry_back()
122      end if 
123      if(npos.ne. digit(i+1) ) deallocate(mod_i)  
124    end do
125     
126    call carry_forward()
127    call carry_back()
128    call display()
129    call finalize()
130    if( allocated(tmp_fact) ) deallocate( tmp_fact )
131  end subroutine factorial 
132end module factorial_mod
133 
134 
135program main 
136  use factorial_mod
137  implicit none
138  integer(id) :: t(3)
139  integer(id) :: num   
140  write(*,*) "Please input a positive integer number(not greater than 13500):"
141  read(*,*) num
142  call system_clock(t(1),t(3))
143  call factorial(num)
144  call system_clock(t(2),t(3))
145  write(*,"(a,f7.3,a)") " Time :",(t(2)-t(1))*1.0/t(3) ," s"
146end program main

Makefile：

[Bash shell] 纯文本查看 复制代码

01#!/bin/bash
02FC = gfortran
03FF = -static -O2
04SC = factorial.f90
05XX = x.out 
06 
07$(XX) :
08        $(FC) $(FF) $(SC) -o $(XX) 
09          
10clean:          
11        rm -rf *.mod *.txt $(XX)