大数求阶乘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] 纯文本查看 复制代码

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：

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

#!/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的阶乘