大数计算

abandon

n! means n  (n  1)  ...  3  2  1
For example, 10! = 10  9  ...  3  2  1 = 3628800,
and the sum of the digits in the number 10! is 3 + 6 + 2 + 8 + 8 + 0 + 0 = 27.
Find the sum of the digits in the number 100!
我的程序:(计算结果有错误)。 求大神帮忙解答，为什么在除以十取余后有错误。有不有更好的算法。

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

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

jason388

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

楚香饭

可以考虑使用大数模块。

经过计算。100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
其各项相加为      648

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

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

大数模块，在本站主站“代码”栏目中有提供

http://www.fcode.cn/code_gen-46-1.html

abandon

chuxf 发表于 2014-3-5 09:50
可以考虑使用大数模块。

经过计算。100! = 93326215443944152681699238856266700490715968264381621468592 ...

非常感谢啊

lm_lxt

最有用的帖子！

fcode

这个是一个实际应用呢，还是某个课程的课堂作业呢？

我觉得用大数模块是一种“暴力”手段，不至于成为课堂作业，否则对学生的要求未免太高了点。

是不是有更简便的计算方式？

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 发表于 2014-3-6 06:44
所谓“另一种求大数阶乘的算法”其实很简单：用数组a存储阶乘结果的每一位数字（a(0)存储阶乘结果的位数，a ...

嗯，算法原理是没有什么难点。
每一个 a 其实可以存储4个或者8个数字的，这就有点类似BCD码了。

而BCD码的计算或处理，则有很多现成的代码。

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