Fortran Coder

查看: 760|回复: 5

[数理统计] 大数求阶乘Fortran实现初探

[复制链接]

38

帖子

5

主题

0

精华

熟手

超凡脱俗

F 币
250 元
贡献
144 点
发表于 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)



测试结果:

r.png
分别为10、  100、 300、  1000的阶乘

r10000.png
10000的阶乘(0.224秒)
ps:屏幕太小放不下了













天下英雄出我辈,一入江湖岁月催。

鸿图霸业谈笑间,不胜人生一场醉。
回复

使用道具 举报

568

帖子

0

主题

0

精华

大师

F 币
396 元
贡献
263 点

规矩勋章元老勋章新人勋章水王勋章

QQ
发表于 2018-2-4 22:23:33 | 显示全部楼层
牛逼牛逼
回复

使用道具 举报

462

帖子

3

主题

0

精华

大宗师

F 币
3115 元
贡献
1848 点

水王勋章元老勋章热心勋章

发表于 2018-2-24 11:06:59 | 显示全部楼层

38

帖子

5

主题

0

精华

熟手

超凡脱俗

F 币
250 元
贡献
144 点
 楼主| 发表于 2018-4-6 10:02:41 | 显示全部楼层
pasuka 发表于 2018-2-24 11:06
lz为啥不考虑用FFT变换加速整数相乘呢?
http://numbers.computation.free.fr/Constants/Algorithms/repres ...

直接算的,没考虑什么其他的算法,我试试这个
天下英雄出我辈,一入江湖岁月催。

鸿图霸业谈笑间,不胜人生一场醉。

2

帖子

1

主题

0

精华

入门

F 币
73 元
贡献
30 点
发表于 2018-4-12 21:54:49 | 显示全部楼层
牛叉
回复

使用道具 举报

21

帖子

4

主题

0

精华

熟手

F 币
147 元
贡献
77 点

规矩勋章爱心勋章

发表于 2018-6-20 20:27:06 | 显示全部楼层
你为何这么吊
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2018-10-16 13:19

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表