[数值问题] 大数计算

program Console1
    implicit none
    integer i
    real*8 a
    integer b,sum
    a=1
    do i=1,100
    a=a*i 
    end do
    sum=0
    do while (a>0)
    b=mod(a,10.)
    a=a/10.
    sum=sum+b
    end do
    print *, sum
    end program Console1

program www_fcode_cn
  use big_integer_module
  implicit none
  integer i
  type(big_integer) :: a
  integer :: sum , b
  a = 1
  do i=1,100
    a = a * i
    !if ( modulo(a,10)==0 ) a = a/10
    !// 上面一行添加后可减少位数
  end do
  sum=0
  call print_big (a)
  write(*,*)
  do while (a>0)
    b=modulo(a,10)
    a = a / 10
    sum=sum+b
  end do
  write(*,*) sum
end program www_fcode_cn

! 计算代码来自 [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