FORTRAN 95准确度很大的除法

Question

我写了一个子程序来计算二项式系数（n选择k），它将用于n和中等k值的大值（因此接近给定n的最大值）。子程序代码是

! 
! Subroutine to calculate combinatoric term (n choose x) 
! 
subroutine COMBO(m,y,combin) 
implicit none 
integer m,y,b1,b2,maxx,minn 
real (kind=16) combin,temp1,temp2 
temp1 = real(1,16) 
temp2 = real(1,16) 
maxx = MAX(y,m-y) 
minn = MIN(y,m-y) 
do b1 = maxx+1,m 
temp1 = temp1*real(b1,16) 
enddo 
do b2 = 1,minn 
temp2 = temp2*real(b2,16) 
enddo 
combin = temp1/temp2 
end subroutine COMBO 

这适用于中等大的值。但是，如果我用N = 100和K = 55我得到以下

61448471214136179596720592959.998

小数部分显然是错误的，因为组合总是整数。主程序是

program chkint 
implicit none 
integer i,j,n,k 
real (kind=16) cmb 

print*,"what is n?" 
read*,n 
print*,"what is k?" 
read*,k 

call combo(n,k,cmb) 
print*,cmb 

end 

顺便说一句：kind = 16是我的机器上的四倍精度。

感谢

Answer 1

取决于你使用整数magniitude，有可能是系统上使用的扩展精度整数便携式手段。标准模块ISO_FORTRAN_ENV应该定义一些整数类型（INT8，INT16，INT32，INT64等）。从Fortran 2008开始，这些包含在数组INTEGER_KINDS中，但目前这似乎很实用。

如果您的房间还没有用完，您有两种选择：重新设定您的问题，使其在您的精度范围内工作，或寻找扩展或任意精度的数学库。令人高兴的是，许多这样的库都在LBNL中列出http://crd-legacy.lbl.gov/~dhbailey/mpdist/

更新：

selected_int_kind是创建扩展精度值的另一种动态的手段，但你仍然需要知道哪些“之类的价值观被编译器接受。例如：

program kind_int 
    use iso_fortran_env, only: output_unit, INT8, INT16, INT32, INT64!, & 
!  INTEGER_KINDS 

    implicit none 

    integer :: i 

1 format(A) 
2 format('Kind name = ', A5, T20, 'kind value = ',I2,    & 
     T40, 'maximum = ', I19) 
3 format('Available kind values: ', I2, *(', ', I2)) 
    continue 

    write(output_unit, 1) 'Maximum value of supported integer kinds:' 
    write(output_unit, 2) 'INT8 ', INT8, huge(1_INT8) 
    write(output_unit, 2) 'INT16', INT16, huge(1_INT16) 
    write(output_unit, 2) 'INT32', INT32, huge(1_INT32) 
    write(output_unit, 2) 'INT64', INT64, huge(1_INT64) 
    write(output_unit, 2) '16 ', 16, huge(selected_int_kind(16)) 
    write(output_unit, 2) '32 ', 32, huge(selected_int_kind(32)) 

! write(output_unit, 3) INTEGER_KINDS 

    stop 
end program kind_int 

给出LINUX64系统上gfortran以下：

Maximum value of supported integer kinds: 
Kind name = INT8 kind value = 1  maximum =          127 
Kind name = INT16 kind value = 2  maximum =         32767 
Kind name = INT32 kind value = 4  maximum =        2147483647 
Kind name = INT64 kind value = 8  maximum =      9223372036854775807 
Kind name = ik32 kind value = 16  maximum = 170141183460469231731687303715884105727 
Kind name = ikxx kind value = 16  maximum = 170141183460469231731687303715884105727 

Selected int kind value = 1   kind value = 1 
Selected int kind value = 2   kind value = 1 
Selected int kind value = 3   kind value = 2 
Selected int kind value = 4   kind value = 2 
Selected int kind value = 5   kind value = 4 
Selected int kind value = 6   kind value = 4 
Selected int kind value = 7   kind value = 4 
Selected int kind value = 8   kind value = 4 
Selected int kind value = 9   kind value = 4 
Selected int kind value = 10   kind value = 8 
Selected int kind value = 11   kind value = 8 
Selected int kind value = 12   kind value = 8 
Selected int kind value = 13   kind value = 8 
Selected int kind value = 14   kind value = 8 
Selected int kind value = 15   kind value = 8 
Selected int kind value = 16   kind value = 8 
Selected int kind value = 17   kind value = 8 
Selected int kind value = 18   kind value = 8 
Selected int kind value = 19   kind value = 16 
Selected int kind value = 20   kind value = 16 
Selected int kind value = 21   kind value = 16 
Selected int kind value = 22   kind value = 16 
Selected int kind value = 23   kind value = 16 
Selected int kind value = 24   kind value = 16 
Selected int kind value = 25   kind value = 16 
Selected int kind value = 26   kind value = 16 
Selected int kind value = 27   kind value = 16 
Selected int kind value = 28   kind value = 16 
Selected int kind value = 29   kind value = 16 
Selected int kind value = 30   kind value = 16 
Selected int kind value = 31   kind value = 16 
Selected int kind value = 32   kind value = 16 
Selected int kind value = 33   kind value = 16 
Selected int kind value = 34   kind value = 16 
Selected int kind value = 35   kind value = 16 
Selected int kind value = 36   kind value = 16 
Selected int kind value = 37   kind value = 16 
Selected int kind value = 38   kind value = 16 
Selected int kind value = 39   kind value = -1 

随着微小修改，该程序产生以下使用ifort在同一系统上：

Maximum value of supported integer kinds: 
Kind name = INT8 kind value = 1  maximum =          127 
Kind name = INT16 kind value = 2  maximum =         32767 
Kind name = INT32 kind value = 4  maximum =        2147483647 
Kind name = INT64 kind value = 8  maximum =      9223372036854775807 

Selected int kind value = 1   kind value = 1 
Selected int kind value = 2   kind value = 1 
Selected int kind value = 3   kind value = 2 
Selected int kind value = 4   kind value = 2 
Selected int kind value = 5   kind value = 4 
Selected int kind value = 6   kind value = 4 
Selected int kind value = 7   kind value = 4 
Selected int kind value = 8   kind value = 4 
Selected int kind value = 9   kind value = 4 
Selected int kind value = 10   kind value = 8 
Selected int kind value = 11   kind value = 8 
Selected int kind value = 12   kind value = 8 
Selected int kind value = 13   kind value = 8 
Selected int kind value = 14   kind value = 8 
Selected int kind value = 15   kind value = 8 
Selected int kind value = 16   kind value = 8 
Selected int kind value = 17   kind value = 8 
Selected int kind value = 18   kind value = 8 
Selected int kind value = 19   kind value = -1 

正如Vladimir F所建议的那样，通过使用selected_int_type可以获得一些额外的精度，但它会随着硬件和编译器的不同而变化，并且可能会或可能不会给您足够的精度在问题上。