FORTRAN 95 对非常大数的除法精度答案

【问题标题】：FORTRAN 95 accuracy for division of very large numbersFORTRAN 95 对非常大数的除法精度
【发布时间】：2014-12-06 05:59:26
【问题描述】：

我编写了一个子程序来计算二项式系数（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 在我的机器上是四精度。

谢谢

【问题讨论】：

你用的是什么编译器？我收到了61448471214136179596720592960.0000，这是一个预期的答案。小数点在那里，因为cmb 是真实的。由于浮点计算，您得到的是错误。
感谢您的回答。我正在使用 Absoft FORTRAN 95。如果我将 cmb 声明为整数，它根本不起作用。由于浮点计算导致的错误，我能做些什么吗？
我没有，用的是 Intel 的。它还可能取决于 x86 或 x64 操作系统/CPU。如果你有 x64 编译器和系统 - 为什么不使用 selected_int_kind(32)
如果您关心精确的数值精度，请不要使用不可移植的kind=16 来指定四倍精度，但可以使用许多其他可移植方式。它已经在这里重复了数百万次。即使real*16 比 IMO 更好，但最好是使用具有工作精度的命名常量，可以轻松更改。
在这种情况下，请尝试使用 Cheery 建议的最大整数类型。

标签： fortran

【解决方案1】：

根据您使用的整数的大小，可能有可移植的方法在您的系统上使用扩展精度整数。标准模块 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 获得一些额外的精度，但它会因硬件和编译器而异，并且可能会或可能不会根据问题为您提供足够的精度。

【讨论】：

我一般会同意你的观点，但在这种情况下需要int128，并且它不是由 ISO_FORTRAN_ENV 定义的，即使 128 位整数在某些编译器中实际上是可用的。这就是为什么我更喜欢 cmets 中已经建议的selected_int_kind，而不是这个问题的 ISO_FORTRAN_ENV 种类常量。
这种方法的一个问题是您需要知道哪些类型的值可以作为selected_int_kind 的参数；我将通过更新我的回复来澄清这一点
什么？你只需提供你需要的位数，如果你得到它。很简单，只要数出问题中的数字即可。
也许，也许不是。见上文。
当然！它不是动态的，它是完全静态的！关键是你可能从中得到可用的类型，但你肯定不会从iso_fortran_env命名常量中得到它（还有其他方法）。我从来没有说过你可以保证从selected_int_kind 获得所需的种类，这就是为什么我在之前的评论中写道你需要看看你是否能得到它以及为什么我在第一个评论中写道它可以在 some 编译器。