大数计算

jason388

kind=8的实数只有15位左右的精度，无法精确表示100的阶乘。

jason388

fcode 发表于 2014-3-5 10:52
这个是一个实际应用呢，还是某个课程的课堂作业呢？

我觉得用大数模块是一种“暴力”手段，不至于成为课堂 ...

计算大数阶乘有很多算法（http://www.luschny.de/math/factorial/FastFactorialFunctions.htm），牟尼的“另一种求大数阶乘的算法”应该比较容易改成Fortran，实现楼主要求的累加应该也很容易，有时间可以试试。

jason388

所谓“另一种求大数阶乘的算法”其实很简单：用数组a存储阶乘结果的每一位数字（a(0)存储阶乘结果的位数，a(1),a(2),...,a(n)分别存储个位，十位，...，最高位），然后按n!=n*(n-1)!，用n乘数组中存储的(n-1)!阶乘的每一位，利用乘法原理重新生成a的每一个元素。

jason388

chuxf 发表于 2014-3-6 09:50
嗯，算法原理是没有什么难点。
每一个 a 其实可以存储4个或者8个数字的，这就有点类似BCD码了。

改自 http://blog.csdn.net/liangbch/article/details/7633574的代码。

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

! 计算代码来自 [url]http://blog.csdn.net/liangbch/article/details/7633574[/url]
program factorial
  implicit none   
  integer :: i,j,m,n,c,prod,rad,istat,len,buff_len
  integer,allocatable :: buff(:)
  real(8) :: pi=acos(-1.0_8)

  ! rad may equal other values but no output for n! 
  write(*,'(a)',advance='no')"radix[10,100,1000,10000]="
  read*,rad
  if(mod(rad,10) /= 0)then
     rad=10
  end if
  
  m=int(log10(real(rad,8)))
  
  ! n should less than 10^9 because range of integer is 9 by default.
  write(*,'(a)',advance='no')"n[<10^9]="
  read*,n

  if(n==1)then
     buff_len=1
  else
     buff_len=ceiling((n*log(real(n,8))-real(n,8)+log(2.0_8*n*pi)/2.0_8)/log(10.0_8))
  end if
  buff_len=buff_len/m+1

  allocate(buff(buff_len),stat=istat)
  if(istat /= 0)then
     stop "allocate array failed."
  end if
    
  len=1
  buff(1)=1

  do i=1,n
     c=0
     do j=1,len
        prod=buff(j)*i+c
        buff(j)=mod(prod,rad)
        c=prod/rad 
     end do

     do
        if (c==0) exit
        buff(len+1)=mod(c,rad)
        c=c/rad
        len=len+1
     end do
  end do

  if(m<5)then
     select case(m)
     case(1)
        print '(i0,"!=",i0,*(i1.1))',n,buff(len),buff(len-1:1:-1)
        print '("sum of digits = ", i0)',sum(buff(1:len))
     case(2)
        print '(i0,"!=",i0,*(i2.2))',n,buff(len),buff(len-1:1:-1)
     case(3)
        print '(i0,"!=",i0,*(i3.3))',n,buff(len),buff(len-1:1:-1)
     case(4)
        print '(i0,"!=",i0,*(i4.4))',n,buff(len),buff(len-1:1:-1)
     end select
  end if

  len=(len-1)*m+log10(real(buff(len)))+1
  print '("digits number of ",i0,"! = ",i0)', n,len

  deallocate(buff)
end program factorial