大数计算
n! means n(n1)...321For example, 10! = 109...321 = 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!
我的程序:(计算结果有错误)。 求大神帮忙解答，为什么在除以十取余后有错误。有不有更好的算法。
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
kind=8的实数只有15位左右的精度，无法精确表示100的阶乘。 可以考虑使用大数模块。
经过计算。100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
其各项相加为 648
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 chuxf 发表于 2014-3-5 09:50
可以考虑使用大数模块。
经过计算。100! = 93326215443944152681699238856266700490715968264381621468592 ...
非常感谢啊 最有用的帖子！ 这个是一个实际应用呢，还是某个课程的课堂作业呢？
我觉得用大数模块是一种“暴力”手段，不至于成为课堂作业，否则对学生的要求未免太高了点。
是不是有更简便的计算方式？ fcode 发表于 2014-3-5 10:52
这个是一个实际应用呢，还是某个课程的课堂作业呢？
我觉得用大数模块是一种“暴力”手段，不至于成为课堂 ...
计算大数阶乘有很多算法（http://www.luschny.de/math/factorial/FastFactorialFunctions.htm），牟尼的“另一种求大数阶乘的算法”应该比较容易改成Fortran，实现楼主要求的累加应该也很容易，有时间可以试试。 所谓“另一种求大数阶乘的算法”其实很简单：用数组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码的计算或处理，则有很多现成的代码。 chuxf 发表于 2014-3-6 09:50
嗯，算法原理是没有什么难点。
每一个 a 其实可以存储4个或者8个数字的，这就有点类似BCD码了。
改自 http://blog.csdn.net/liangbch/article/details/7633574的代码。
! 计算代码来自 http://blog.csdn.net/liangbch/article/details/7633574
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="
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
页:
[1]
2